Theoretical Studies of Anisotropic Melting of Ice Induced by Ultrafast Nonthermal Heating

Water and ice are routinely studied with X-rays to reveal their diverse structures and anomalous properties. We employ a hybrid collisional-radiative/molecular-dynamics method to explore how femtosecond X-ray pulses interact with hexagonal ice. We find that ice makes a phase transition into a crystalline plasma where its initial structure is maintained up to tens of femtoseconds. The ultrafast melting process occurs anisotropically, where different geometric configurations of the structure melt on different time scales. The transient state and anisotropic melting of crystals can be captured by X-ray diffraction, which impacts any study of crystalline structures probed by femtosecond X-ray lasers.


Additional information
The structure of ice The length scales which describe the RDF peaks are shown in the table SI. 1.

Atomic model convergence for CR simulation
We model the interaction between photons and matter based on atomic physics, even though the bound electrons in the water molecule are more accurately described by molecular orbitals.This is motivated by the fact that core electrons are rather localized, especially when the atom is positively charged.Furthermore, the error in the ionization energy and cross-sections from neglecting molecular orbitals will be small since the molecule will be in a plasma, which causes many other perturbations to the orbitals of the molecule.
The screened hydrogenic model (SHM) in Cretin 2 is developed by averaging levels from a more complicated atomic model (which for instance can contain l quantum numbers) where the atom is only described by its n quantum numbers.The resulting model can still be reduced, by for instance reducing the total number of available n levels, limiting the number of core holes in higher n quantum number levels, constraining the total number of electronic excitations and lowering the largest principal quantum number used for excitations.
We studied the convergence of the free electron temperature as a function of the number of principal quantum numbers in the model.The convergence was conducted by simulating the SI. 2 interaction between the pulse and the sample, using the largest intensity in the manuscript (10 19 W/cm 2 ).The average electron temperature of the last 10 fs of the 50 fs pulse was used for analyzing the convergence.The result is shown in figure (SI.1).The figure shows that the electron temperature has a small change with increasing principal quantum number (PQN).Thus, the final model of oxygen was set to have a PQN of 10.We allowed for 8 electrons to be excited in the atom at the same moment.Furthermore, the model was allowed to contain two-core hole states.The convergence study allowed us to utilize a reduced model, which provided a speed-up in the computation of the production simulations which spanned multiple intensities.

Isotropic scattered intensity
The time-resolved scattered intensity from ice can be calculated by the formula 3 where the first term is the self-scattering term I sel f (q,t), which does not provide any information regarding the relative position of the atoms.The second part corresponding to the inter-molecular SI. 3 term I inter (q,t) encodes structural information, through the use of summation of the partial structure factors S α,β .The quantity S α,β is computed using the RDFs g α,β (r) in the hybrid plasma/molecular dynamics simulations Given that the two atomic species we have in our system is oxygen (O) and hydrogen (H), we get α, β = O, H. Therefore, we can write the structure factor due to the RDF of O-H being symmetric g OH (r) = g HO (r) as,

Coherent diffraction
We computed the coherent diffraction from single crystals using the following equation I(q,t) = I 0 (t) Ω(q) |F(q,t)| 2 (dσ /dΩ) T , (SI. 4) where I 0 (t) is the illumination fluence assumed to be homogeneous across the sample, Ω(q) is pixelwise detector solid angle, F(q,t) is the time-dependent structure factor, and is the Thompson scattering cross-section with ε and ε′ representing the polarization of the incident and scattered X-rays, respectively 4 .The structure factor is defined as Here c i (t) represent the population and f 0,i (q) the atomic scattering function of the electronic configuration i for example O-1s 1 2s 2 2p 4 .Similar methods to determine structure factors have been SI. 4 reported in the literature 5 .We assumed spherical symmetry on the atomic scattering function f 0,i (q) = 4 π dr |ψ i (r)| 2 sinc(q r) r 2 .(SI. 6) We computed the radial atomic wavefunctions ψ i (r) for each configuration with electronic structure code RSPt and their time-dependent populations with collisional radiative code cretin.The diffraction simulations ran in Julia 1.7.2 taking advantage of the non-uniform fast Fourier transforms package FINUFFT 3.1.0 6, which allowed us to sample the real space at atomic positions.We sliced the incident beams with 10 different energies within 1% of the beam's photon energy to represent the bandwidth and averaged the resulting diffraction patterns incoherently.We disregarded polarization, quantum efficiency effects, and discretization of the photons upon detection.We used a tolerance of 1 × 10 −10 to solve the fast Fourier transform.

Molecular dynamics bonded interactions
For running classical molecular dynamics simulations, a force field, which describes the interactions between the atoms in the system, is needed.In our simulations we assume that the sample is a plasma, and therefore we must alter the interaction potentials used.Since the assumption of Cretin is that the sample is plasma, the method is only valid if the native sample transitions into the plasma phase in a time-scale which is short compared to the time-scale of atomic motion.We note that during the first 1-2 fs of the interaction between the photons and the matter the sample is not described by a plasma, but is so after, as can be observed in figure (SI.11).
Since the atoms occupy essentially their native positions in the time it takes to reach the plasma phase, we are sure that we can base our physical description of the system on plasma physics.
Even though a plasma usually describes a system where bonds do not exist, and the probe will turn the sample into a plasma, the charge distribution will depend on the parameters used, and thus it is likely that there will be a fraction of atoms which are neutral, as is proven for both hydrogen and oxygen in figures (SI. 9) and (SI.10).In this case, one still wants the bonds to be intact.
Therefore, we decided to use Morse potentials 7 for atoms which are bonded in the native system.
The potential errors introduced by the parameters used are essentially removed by starting the simulations at 0 K.This results in that the propagation of the molecular dynamics is dictated solely by the collisional-radiative calculations.The use of the Morse potential is further motivated by the fact that one needs to make sure that the atoms' electronic clouds repel as the atoms approach each other.The Morse potential is defined as where the depth of the well is D e , the well steepness β and the equilibrium distance r e all depend on the two atoms involved.We validated the parameters by checking the stability of a simulation starting at 0 K using only the Morse and Lennard-Jones potential.We observed that the temperature did not increase and therefore concluded that the force field would be adequate for running the plasma simulations.The converged parameters were then used in the production simulations for running the hybrid plasma/MD simulations of ice.

The screened non-bonded interactions
In this section, we present the potentials that are used to compute the non-bonded interactions between the atoms.The screened Coulomb potential was defined as 8 where f = 138.935458kJ mol −1 nm e −2 , and the scalar force was given by F = −∇V (r) (SI. 9) SI. 6 The screened Lennard Jones potential is defined as and the force . (SI. 11) The parameters c 6 and c 12 depend on the two atom types involved.In our simulations, they dynamically change due to the reduction of bound electrons in the atoms.As one removes electrons from the atoms, the ionic radius will reduce and the atoms can reach closer to each other.Therefore, the parameters are changed such that this effect is present.Furthermore, since the Lennard-Jones interaction is an interaction originating from the electromagnetic force, we apply the same screening factor like in the Coulomb potential case, in order to get the effects of the electronic environment.
Given the important length scales as described by the peaks in the RDFs, we compute the screened force at the first two peaks (r x = r 1 , r 2 ) as a function of time, normalized by the regular Coulomb interaction, 2) shows how the forces at these length scales are greatly reduced due to the screening of the free electrons, as computed by equation (SI.12).Leaving out the effect of the electron cloud would result in an increased Coulomb force and thus a much more rapid expansion of the system.

Additional results
The

Charge dynamics in CR simulations
The main data collected from CR simulations to run classical molecular dynamics (MD) in GROMACS 9 is the fractional charge state of each atomic species (as seen in figure (SI.9), the free electron temperature and density.Blow we show the charge dynamics for oxygen and hydrogen atoms at different conditions.

Debye screening to model Coulomb interactions
For the plasma conditions studied here, the Debye screening model is not valid for all length scales.
Debye screening is appropriate for length scales much smaller than the dimensions of the system, which is always true for our simulations.However, we may not have enough charged particles in the Debye sphere to make the concept statistically significant as the Debye length becomes shorter.
Therefore, the concept of Debye screening breaks down close to the ion.An ion sphere model instead is most appropriate for a strongly coupled condensed matter system closer to the ion.For long range screening where we approach long distance away from the ion, the Debye length will be statistically valid 10 .We believe the way we incorporate the effects of the electrons on the ion-ion Coulomb interaction is a good approximation for this system.The complete mathematical form of the screening is something else than we use 10 , but our simplified formula will can capture well the essence of the screening phenomena.Comparing these two screening models in water close to the ion, we do not see enough difference to motivate using a hybrid screening model 10 .Finally, when SI. 16 simulating the same system with and without screening, the only notable difference in the dynamics is the time-scale for the melting, and not the process of melting.

Stewart-Pyatt continuum lowering
In our collisional-radiative calculations used to compute the photon-matter interactions, we utilize the Stewart-Pyatt model for the continuum lowering.There is some debate in the literature about the appropriateness of the Steward-Pyatt model to capture continuum lowering effects in favor of the Ecker-Kroll model 11,12 .Other papers suggest that Steward-Pyatt is most appropriate or that lowering is well captured by these parametric models 13,14 .Recent developments in the field involve the use of density functional theory or other methods 15,16 .To our knowledge there is no simple model that can be implemented in collisional-radiative simulations.The relevant physics to capture the transition from a degenerate plasma (condensed matter like) to warm dense matter implemented in Cretin 2 include: inclusion of the effects of electron degeneracy on the rates, opacities, emissitives and collisional ionization.This code has previously been employed to study x-ray induced dynamics in matter with comparison to experiments 17 and has along with appropriate atomic data been benchmarked to other similar codes 18 .Similar to our pipeline, other codes 19,20 use MD simulation supplemented by ionization dynamics from rate equations to model x-ray induced atomic damage in crystalline systems.While the validity of the SP model is been debated, it still gives adequate results which can be used to make conclusions in this study.Our observable is not the ionization state (or free electron dynamics), but the diffraction pattern, which depends on electronic occupation and atomic positions.We therefore do not expect much differences compared to a different continuum lowering model.

SI. 17
Temperature and root-mean-square displacement comparison between CR

and MD simulations
Information regarding structure is not available in Cretin.The dynamics of the structure are instead studied by taking data from Cretin to run MD in GROMACS.A advantage of using MD is that we can study dynamics for a system which is not isotropic.Thus, one should expect that a quantity such as temperature should be at least on the same order between the two methods.In figure (SI.14)   we compare the temperature from the MD and CR calculations.We note that the temperatures are indeed of similar order of magnitude with time for the highest intensity.Even though temperatures are only of the same order of magnitude, the differences are not expected to change the interpretation of the MD results.Furthermore, we note that the temperature evolution of different trajectories in the MD simulations, where the charge distribution is assigned differently in space, are reproducible.
We started simulating 10 different trajectories with different random seeds.Each simulation distributed the charge states as seen in figure (SI.9) differently.Out of the 10 trajectories, some of the simulations had to be removed since the temperatures increased in a single time-step in a non-physical way.We are not completely certain of its origin, but we theorize that since electron recombination causes a change in the Lennard-Jones parameters, it may be that two atoms that are at separation d t 1 at time t 1 , will be stable (i.e.no large LJ-forces) but at time t 2 , when one or both atoms have had recombination of an electron, will make this separation for the new LJ-parameters more repulsive and thus result in large forces.This would, in turn, lead to large spikes in the temperature.

Bond breaking with ab-initio methods
In our previous work, where the dynamics of charged amino-acids and peptides were studied, we concluded that some bonds could be robust to ionization 21,22 .For the water molecule, we justify the use of Morse potentials (bonds) by running density functional theory (DFT) simulations scanning different net charges (+0, +1, +2, +3, +4) and detecting the threshold charge for breaking the bonds.
We performed all-electron unrestricted Kohn-Sham calculations in ORCA 23 , with B3LYP as the exchange correlation functional and TZVPP as the basis set for both energy minimization and the molecular dynamics calculations.A single trajectory starting at 300 K was computed for 75 fs, where each time-step was 0.5 fs.The results in figure (SI.15) show that both O-H bonds are intact for up to +2 charge, while for +3, one of the bonds breaks.For higher charges than +3, both bonds break.Thus, we show using ab-initio molecular dynamics, that the bonds in water can be intact, even though the molecule is charged.

Fig. SI. 1 :
Fig. SI.1: Electronic temperature as a function of principal quantum number n used in the atomic model of oxygen.

Fig. SI. 2 :
Fig. SI.2: Screened Coulomb interaction normalized by the regular Coulomb force at different length scales.The data for the Debye length is from a simulation with an intensity of 1e18 Wcm −2 .The filled line is for O-O, dashed is for H-H and dotted for O-H.

Fig. SI. 8 :
Fig. SI.8: Average charge of oxygen (filled lines) and hydrogen (dashed lines) for a 50 fs flat shaped pulse XFEL pulse with 8 keV photon energy, shown for different XFEL intensities.The range of intensity is chosen to achieve an average ionization of oxygen that will allow for coherent scattering during the X-ray pulse.

Fig. SI. 9 :
Fig. SI.9: Example of the fractional charge distribution for an intensity of 10 19 Wcm −2 and photon energy 8 keV.

Figure
Figure (SI.SI.13) shows in particular the PADF contact r 2 from different angle (top view, hexagonal ice).

Fig. SI. 14 :
Fig. SI.14: The ion temperature for an intensity of 1e19 W/cm 2 and photon energy 8 keV.

Fig. SI. 15 :
Fig. SI.15: Distance between O-H as a function of time from ab-initio molecular dynamics of charged H 2 O, as calculated by ORCA.The dotted/non-dotted line each show one of the bonds to the hydrogens.

Table SI . 1 :
Table of the first three peaks for the radial distribution function of the ice structure, as computed using GenIce 1 .The values are in nanometers.