Various damage mechanisms in carbon and silicon materials under femtosecond x-ray irradiation

We review the results of our research on damage mechanisms in materials irradiated with femtosecond free-electron-laser (FEL) pulses. They were obtained using our hybrid approach, XTANT (X-ray-induced Thermal And Nonthermal Transitions). Various damage mechanisms are discussed with respect to the pulse fluence and material properties on examples of diamond, amorphous carbon, C60 crystal, and silicon. We indicate conditions: producing thermal melting of targets as a result of electron-ion energy exchange; nonthermal phase transitions due to modification of the interatomic potential; Coulomb explosion due to accumulated net charge in finite-size systems; spallation or ablation at higher fluences due to detachment of sample fragments; and warm dense matter formation. Transient optical coefficients are compared with experimental data whenever available, proving the validity of our modeling approach. Predicted diffraction patterns can be compared with the results of ongoing or future FEL experiments. Limitations of our model and possible future directions of development are outlined.

Often advanced simulation tools for treating FELexcited matter are based on the density functional theory molecular dynamics (DFT-MD) schemes, see, e.g., [23,[25][26][27]. Various improvements upon the standard DFT packages are being tested, such as incorporation of corehole excitations [28], or attempts to reconcile bound and free electrons into a unique formalism [29][30][31][32]. However, important non-adiabatic effects governing electron-ion (electron-phonon) energy exchange are difficult to incorporate into ab initio simulations. The state of the art in the solid state community is so far an approximate treatment of electron-phonon scattering, valid only near the room temperature of atoms [33,34]. In the ab initio femtochemistry, advanced techniques allow for a treatment of non-adiabatic electron-ion coupling [35]. Attempts have been made to introduce similar methods into the solid-state models [36]. But, at present, there is no standard methodology enabling to incorporate electron-ion interaction and the resulting energy exchange within ab initio approaches for highly excited many-body systems.
Molecular dynamics (MD) methods with classical potentials or force fields are among the most commonly used numerical tools for large-scale simulations in the solid-state modeling, see, e.g., [37][38][39]. Utilization of classical interatomic potentials makes MD simulations very efficient and capable of treating large systems [39,40]. However, such methods require predefined force fields, whose applicability is limited to solids under low electronic excitation. A few attempts have been made to develop potentials depending on electronic temperature and/or number of excited electrons, e.g., [41,42].
Classical Monte Carlo (MC) methods are often applied to simulate FEL-induced transport of electrons, photons or other particles in matter [43]. MC methods are suitable for simulating behavior of high-energy classical particles, whereas the propagation of low-energy particles involves strong quantum-mechanical effects which can hardly be treated. This limits the applicability of asymptotictrajectory MC schemes.
Apart from the integral methods, differential methods are also used to treat systems with electronic excitations. Models operating with ensembles instead of individual particles can be derived from the Liouville equation. They are typically based on the single-particle distribution equations, such as kinetic equations and collision integrals [44,45]. Kinetic equations are capable of treating nonequilibrium evolution of solids. However, they are usually relying on a free-electron approximation or on a predefined band structure of a material [46]. In addition, the electronion coupling, defining energy exchange between the electronic and the atomic system, is one of the least known key parameters in these approaches.
The next level of approximation is hydro- [47][48][49][50] or thermodynamic methods [51,52]. The most commonly used is the two-temperature model (TTM), which treats electrons and phonons as separate interacting subsystems in local equilibrium. Thermodynamic approaches are widely used due to their simplicity, but are limited to equilibrium conditions (see, e.g., review [53]). TTM is inapplicable under non-equilibrium within electronic or atomic systems; in such cases, kinetic methods must be employed.
Hybrid approaches are an actively developing field in numerical simulations. A hybrid approach combines two or more simulation techniques into a unified model. In such combinations, various approaches strengthen and complement each other, compensating mutually for their shortcomings, in particular, with respect to the overall computational efficiency. A proper combination of different models can alleviate limitations of each individual approach, thereby significantly extending the applicability of the combined model keeping the implementation simple enough. For example, a combination of the TTM with molecular dynamics, TTM-MD [54], was a major step that allowed to simulate experimental data on laser irradiation of solids with a high accuracy [55]. Also, MC-MD combination allowed to treat accurately both atomic and electronic kinetics in case of plasmas [56].
In this paper, we review our hybrid approach designed to treat solids under FEL irradiation. We demonstrate that such an approach is capable of capturing the essential stages of the material evolution under femtosecond X-ray irradiation, starting from photo-induced non-equilibrium electron kinetics on femtosecond timescales, and progressing later towards atomic dynamics of thermal and nonthermal damage formation up to ps timescales. Thermal phase transition is a result of a kinetic energy exchange between hot electrons and atoms (e.g., electron-phonon coupling), occurring due to non-adiabatic coupling between the two systems. Non-thermal phase transition is a consequence of the changes of the potential energy surface for atoms due to the electronic excitations (e.g., non-thermal melting [57,58]).
In what follows, we present a detailed analysis of different damage channels in various carbon and silicon allotropes. Those materials are widely used in FEL-related optical elements, such as, e.g., X-ray mirrors. Their radiation tolerance plays a crucial role for this application.
Our results indicate that different allotropes of carbon and silicon follow different damage mechanisms under FEL irradiation.
With X-ray FELs, one can experimentally observe the evolution of diffraction patterns from irradiated samples with femtosecond resolution (e.g., Ref. [14,59]). Such a scheme is particularly interesting, because it can resolve the changing molecular structure of the sample, resulting from its progressing damage. Thus, throughout the paper, we present diffraction patterns calculated for each material and damage channels discussed. At the end, we discuss the applicability of the presented hybrid approach and possible pathways of its further development.

Hybrid approach
Our recently developed hybrid code X-ray-induced thermal and non-thermal transitions (XTANT) [60] is a combination of a few different schemes interconnected and executed in parallel. These approaches include: (a) MC module tracing X-ray photon absorption, high-energy electron and core hole kinetics; (b) a module describing the evolution of the low-energy electrons, using rate equations and thermodynamic modeling (in a similar manner to the TTM); (c) the Boltzmann collision integral module following the non-adiabatic electron-ion energy exchange; (d) transferable tight binding (TB) model for calculations of the transient electronic band structure and atomic potential energy surface; and (e) MD simulation tool to follow atomic motion. The scheme showing the interconnection of the most important modules is presented in Figure 1. This diagram also indicates data flows between different modules at each time-step of the simulation. They will be described below in more detail.
The chosen combination of approaches relies on the fact that the time-dependent electron distribution function, f e ðE i ; tÞ, affects atomic motion. That is because transient electron distribution enters equations for the atomic potential energy surface (cf. Eq. (10)). To efficiently trace the evolution of the electron distribution function in time, we notice that the typical transient electron distribution after X-ray irradiation has the shape of the so-called 'bump on hot tail' [61]. It combines (i) a (nearly) thermalized fraction of low-energy electrons within the valence band and the bottom of the conduction band, and (ii) a high-energy non-equilibrium tail containing a few highly energetic electrons left after photoionization and Auger decays. It has been proven theoretically that such a transient electron distribution is typical in various materials after an FEL pulse [61][62][63][64]. Recently, predictions on the electron distribution function in FELirradiated aluminum from reference [63] were confirmed experimentally in reference [65].
This particular shape of the electronic distribution function allows to simplify the model by combining two efficient approaches for low-energy and high-energy fractions of the electronic distribution function. It is used in the presented XTANT code as described in detail in Section 2.2.
Our hybrid model was specifically designed to follow the processes occurring in a solid target under irradiation with a femtosecond free-electron laser pulse. It is applicable in a broad range of photon energies: from extreme ultraviolet (XUV, photon energy above a few tens of eV, depending on material) to hard X-rays (∼40 keV photon energy) [13,60,[66][67][68].

Creation and kinetics of high-energy electrons and core holes modeled with MC method
For the modeling of photon absorption, the kinetics of high-energy electrons, and Auger decays of core holes, we apply an event-by-event individual-particle MC scheme [43,69,70]. For XUV and X-rays at fluences currently accessible at FELs, single photon absorption is the dominant interaction channel [71]; other interaction channels such as elastic scattering, multiphoton absorption and inverse Bremsstrahlung (in XUV regime) are at least by two orders of magnitude less probable. An absorbed photon initiates the release of a high-energy photoelectron, and leaves a core hole in an atom.
In the presented approach, the choice of whether the photoabsorption occurs by an excitation of an electron from the valence band or a deep shell is made, using the subshell photoabsorption cross sections taken from the EPDL97 database [72]. Atomic cross sections for photoabsorption are applied, which are a good approximation for core shells. For the valence band, the photoabsorption cross section is obtained from the experimental optical coefficients [73].
The released photoelectron has initial kinetic energy equal to the difference between the photon energy and its ionization potential, I p . Ionization potentials for core shells are taken from the atomic EADL database [74], whereas for the valence or conduction band photoabsorption, the TB calculated transient energy levels are used (Eq. (9) in Sect. 2.4).
High-energy electrons, i.e., those populating states at energies above a certain energy threshold, are treated as classical individual particles within the MC routine. In the presented calculations the threshold is chosen to be E cut = 10 eV, counted from the bottom of the conduction band. Electrons with lower energies are attributed to the ʻlow-energy' domain (see Sect. 2.3) [75]. Alternatively, the threshold can be chosen as equal to the uppermost energy level of the conduction band produced by the TB calculations. Influence of the cut-off value on the results was analyzed in reference [66], showing almost no effect on the results, if varied by a few eV around the 10 eV value.
Each high energy electron can scatter inelastically on the atomic core-shells, if its kinetic energy is higher than the respective binding energy, E e > I p ; otherwise, only the scattering on the valence-band electrons is possible. The electron mean free paths are estimated via the scattering cross sections, s i (E e ), within the first Born approximation in terms of the complex dielectric function (CDF), є(v e , q) [70,76,77]: Here the cross section also depends on the energy hv e gained by the secondary electron in a collision, and is integrated over the transferred momentum q; e denotes the electron charge; m e is the free electron mass; n a is the atomic density; h is the Planck constant; and v is the incident electron velocity corresponding to the kinetic energy E e . The CDF is parameterized for each shell of each element; all the parameters and accuracy checks of the cross sections can be found in reference [77]. The calculated electron inelastic mean free paths showed a very good agreement within a few percent with the NIST database [78] and available experimental data for electron energies above ∼50 eV.
When a high-energy electron collides with the valence or deep shell electrons, the probabilities for these collisions are estimated using the scattering cross sections from equation (1). The initial energy of the secondary electron emitted during an inelastic collision is calculated from the energy conservation, i.e., it equals to the difference between the energy lost by the incident electron and the binding energy of the level from which this electron is being ionized. If the energy of any electron falls below the E cut , this electron is removed from the MC domain and added into the low-energy domain. Detailed description of this inter-domain energy and particle exchange is presented in Section 2.3.
The elastic scattering of electrons on ions is modeled using Mott's cross section with modified Moliere screening parameter [43]. However, this interaction channel produces only a minor energy loss for electrons (and correspondingly negligible energy increase for atoms) [60,66,79]. Most of the energy transfer to the ions is produced by low-energy electrons via non-adiabatic coupling, as will be described in detail (Sect. 2.3).
After a photoionization or impact ionization of a core shell, a hole is left behind. Predominantly, for light elements, this hole will decay via Auger processes [80]. We use the Poissonian probability distribution to model Auger decays. The characteristic hole decay times are taken from the EADL database [74]. When a core hole relaxes during an Auger process, one electron is promoted from an upper shell or the valence band into high energy states of the conduction band, leaving another hole. In case of the valence band, the energy level from which the Auger electron is emitted is chosen randomly among all currently populated levels E i . The Auger electron receives the excess energy and is then treated in the same way as other secondary electrons.
The MC tracing of photons, electrons, and core holes is split into time steps dt equal to the MD time steps. The number of iterations of the MC subroutine is proportional to the photon energy. This is necessary in order to obtain sufficient statistics to follow accurately small numbers of high-energy photo-electrons. E.g., for intermediate photon energies of a few keV, each time step of MC is iterated for more than 30 000 times for reliable statistics. Trajectories of all electrons are propagated simultaneously during each iteration. The calculated electron distributions are then statistically averaged.
As we assume homogeneous excitation, we neglect a contribution of electron transport and heat diffusion to the overall electron kinetics. This can be justified as an X-ray irradiation homogeneously heats up the sample down to a few micrometers depth within the laser spot of typically a few microns size. The approximation of homogeneity allows us to use periodic boundary conditions, but excludes from our considerations the thin near-surface layer, from which the high-energy electrons could escape outside the material. To account for the effects of electron emission whenever needed, we artificially remove an electron from the simulation box after a certain number of collisions, introducing charge non-neutrality. This charge non-neutrality then contributes to additional Coulomb forces acting on atoms (Sect. 2.4).

Low-energy electrons and non-adiabatic coupling modeled with rate equations including Boltzmann collision integral
At each time step, we track how the electrons are distributed between the high-and low-energy domains.
The total number of low-energy electrons, N low e , and their energy, E low e , are calculated knowing how many electrons were excited to the high-energy domain by the incoming photons from the laser pulse, N ph , by the secondary electron collisions, N imp , by Auger decays of core-shell holes, N A , and by electron emission in case of thin films, N em . Electrons which are excited to the energy levels above the cut-off energy, are transferred to the high-energy domain and treated with the MC algorithm (Sect. 2.2). Vice versa, when an electron from the highenergy fraction loses its energy below the cut-off energy, it joins the low-energy domain. Thus, the total number of low energy electrons is calculated as: where N high e ðtÞ denotes the fraction of high-energy electrons that fell into or were excited off the low-energy domain, as counted within the MC module. The total number of electrons in high-and low-energy domains, including those originating from core-shell ionizations and electron emissions, is conserved.
The total energy of low-energy electrons, E low e , is calculated with an equation similar to (2), in the following way: where E low e is the total energy of low-energy electrons; E high e ðtÞ is the energy brought in or out by high-energy electrons that fell into or jumped off the low-energy domain (e.g. by photoabsorption); E imp is the energy delivered during the impact ionization events by the highenergy electrons; E A is the energy delivered by Auger decays of core-shell holes that involve valence or conduction band electrons; and dE elÀion is the energy transferred to (or from) ions.
The energy flux between the electrons and ions at each time step in equation (3) can be calculated via a collision integral I eÀat i;j as: where the summation runs over all the electronic orbitals for transitions between each pair of levels [68], and I eÀat i;j used in this work is the Boltzmann collision integral: where w i,j is the rate for an electron transition between the energy levels i and j; f e (E i ) is a transient electron distribution function (normalized to 2 due to spin degeneracy); and G at (E) is the integrated Maxwellian function for atoms [67]. One could, in principle, use the transient atomic distribution function obtained from the MD simulations, but due to a relatively small number of atoms in the supercell, the fluctuations do not allow us to obtain a smooth function. This, in turn, would introduce a numerical dependence of the energy transfer rate on the number of atoms. For this reason we approximate the atomic distribution in equation (5) as an equilibrium one, using a Maxwellian function.
The following expression suitable for finite-difference implementation for the transition probability is used [68]: where |i(t)⟩ and |j(t)⟩ are eigenstates of the system Hamiltonian.
Note that this expression for the transition rate only reduces to the Fermi's Golden Rule in case of periodic atomic vibrations (phonons), and at long timescales in comparison to a duration of each individual act of scattering [68]. None of these assumptions holds true at subpicosecond timescales in case of irradiation with femtosecond laser pulses; thus, more general expression, equation (6), has to be used in this case.
Knowing the total energy and number of the lowenergy electrons at each time-step, we can estimate their temperature and chemical potential. They are calculated from the 0th and the 1st moments of the Fermi distribution function by solving the inverse problem [81]: where the summations include all energy levels, E i , corresponding to the current band structure of the material (that evolves in time), and E min is the lowest energy level of the valence band. These energy levels are calculated with the TB method (Sect. 2.4, Eq. (9)). The factor 2 in the Fermi-distribution function f e (E i ) accounts for the electron spin; m is the transient chemical potential of the electrons and T e is their temperature (in energy units). We solve the system of equation (7) for the known values of N low e and E low e at each time step by the bisection method. In this way, we find the transient values of m and T e [60,66].
Note that, in principle, the model does not require the low-energy electron distribution to necessarily obey the equilibrium Fermi function. It can work for any transient non-equilibrium distribution function, such as provided by the Boltzmann equation [45,82]. For example, we tested the possibility of independent Fermi distributions of holes in the valence and of electrons in the conduction band in GaAs [83]. For this case, the system of equation (7) can be written for each band independently, and each band can independently exchange the energy with the high-energy electrons in the MC domain. The work on this scheme is in progress, and will not be discussed in the current paper. Correspondingly, below we assume instant low-energy electron thermalization within the entire lowenergy domain (Eq. (7)), which significantly simplifies the calculations.

Transferable tight binding molecular dynamics
To trace material modifications on the level of both electronic and atomic processes, we employ tight-binding molecular dynamics [84][85][86]. This method relies upon transferable TB Hamiltonian to evaluate electronic energy levels (band structure) and the interatomic potential energy surface. Transferable TB means that for given atomic species the parameterizations of the hopping integrals and the repulsive potential are constructed to reproduce several material phases.
The transferable TB Hamiltonian is written as follows [84][85][86]: is the repulsive part describing the effective repulsion of atomic cores and H TB is the attractive part calculated with the TB Hamiltonian. The TB part is constructed of the on-site energies, ih , and the pairwise overlap integrals, t hn ij , within the sp 3 basis set. Those are parameterized functions, which for the case of carbon-based materials can be found in reference [84], and for silicon in reference [85].
In the case of an orthogonal Hamiltonain parameterization, the electron energy levels are obtained by a direct diagonalization: The potential energy surface F({ r ji } , t), needed as an input to the equations of motion for atoms, can be derived from equation (8) within the Born-Oppenheimer (BO) approximation from the Hellmann-Feynman theorem as follows: Here, f e (E i , t) is the transient electron distribution function indicating fractional electron population numbers on the transient energy levels E i given by equation (9). Additional terms beyond the BO approximation result from the energy transferred from the electron in nonadiabatic transitions between the energy levels mediated by the atomic displacements [68]. This transferred energy N. Medvedev et al.: 4open 2018, 1, 3 is then distributed among all the atoms in the simulation box by the appropriate velocity scaling. The calculation of such energy transfer was described in Section 2.3.
The typical laser spot radius for an FEL laser is of the order of a few micrometers, and a photon penetration depth may be also on a few micron scale. This volume corresponds to at least a few billion atoms. Thus, we can choose only a small simulation box (supercell) inside the laser spot with a size much smaller than the laser spot, and apply periodic boundary conditions. The periodic boundary conditions in XTANT can be introduced in two ways. First, we can keep the super-cell vectors constant during the simulation, modeling the NVE ensemble (i.e. constant volume simulation). Alternatively, within the Parrinello-Rahman method [87], we can account for the changing geometry of the supercell. It is traced via additional variables entering the Lagrangian of motion [86,87], for the NPH ensemble (i.e. assuming constant pressure simulation). We use the velocity Verlet algorithm for propagating atomic coordinates and velocities in time [60]. The applied time-step is usually 0.01 fs, ensuring a stable numerical scheme [68].
For the van der Waals forces acting on carbon atoms in case of C 60 crystal or graphite, we employ additional semiempirical Lennard-Jones (6-12) potential, softly cut at short and large distances [88]. Soft cut-off at short distances ensures that it does not overlap with the shortrange forces treated within the TB approach.
In case of charge non-neutrality, which may occur in thin films after electron emission (see Sect. 2.2), unbalanced positive charge is then accounted for as an additional fractional charge equally distributed among all the atoms in the simulation box [88]. It produces additional long-range Coulomb potential, which is not accounted for in the TB.

Data analysis based on optical properties, autocorrelations and diffraction patterns
Throughout the paper, we define the absorption dose as the absorbed energy per atom at the depth, d, equal to the attenuation length of the considered photon, l (i.e. d = l): where D abs is the absorbed dose (eV/atom); F is the incoming fluence into the material; R is the reflectivity of the sample at the wavelength of the incoming pulse; n at is the atomic density. For comparison with experiments, one may use equation (11) to evaluate the incoming fluence from a given dose. Such a connection assumes linear photoabsorption (without multiphoton effects), which is generally a good approximation for photon energies above a few tens of eV, typical for current FELs.
Knowing atomic positions and the supercell vectors from the MD module, one may obtain powder diffraction patterns with available software. We calculate the patterns with help of the Mercury software [89], renormalizing them to the area under the peaks instead of the highest peak height.
We obtain the vibration spectrum from the time evolution of the atomic velocities v k ! ðtÞ in excited solids, with consideration of the large amplitude motion [90][91][92].
Iðt; vÞ ¼ ℱ t ½Iðt; tÞðvÞ where ℱ t denotes the Fourier transform; N at is the number of atoms, and the parameter a is chosen to suppress the autocorrelation function within 200 fs to effectively select the vibrational modes that are present within a given time interval t to ∼t + 200 fs. In this way, the temporal evolution of vibrational modes can be revealed by their frequencies and amplitudes. In particular, it allows to trace disappearance of a harmonicity in the atomic motion, indicating loss of the structure in the crystal, i.e., phase transition to a disordered phase. Transient optical properties of materials, such as reflectivity or transmissivity, can also be measured in pump-probe experiments with a time resolution down to 10 fs [9,10,93]. Thus, along with diffraction patterns the latter can be used for comparison between the experimental and simulation results. Within the linear response theory, optical properties are defined by the CDF. The random phase approximation (RPA) provides the following expression for the dielectric function (the Lindhard formula) [94,95]: Here, E hn = E n À E h is the transition energy between two eigenstates |h⟩ and |n⟩; f e (E h ) and f e (E n ) are the transient occupation numbers of the corresponding states (electron distribution function) as defined above; F hn are the diagonal elements of the oscillator strength matrix [95,96]; V is the volume of the supercell; and 0 is the vacuum permittivity in SI units. A particular choice of the (small) parameter g does not affect the results beyond the broadening of peaks in the CDF [97]. The dependence of the results on g was investigated in detail in reference [98].
The real, e 0 , and imaginary, e 00 , parts of the CDF define the components n and k of a complex index of refractioñ nðvÞ ¼ n þ ik by relations: Using Fresnel's and Snell's laws, the reflectivity coefficient can be expressed as follows [99]: where u is the angle of incidence of the probe pulse. The transmission coefficient of the material also depends on the material thickness d and the wavelength of the incident probe pulse l. In the case of a bulk or a thick layer of a material and ultrashort probe pulse, we assume the first ray propagation with no interference effects included from multiple reflections on the material boundaries [96,99]: The absorption coefficient can then be obtained from the normalization condition: 3 Results

Low fluence 3.1.1 Non-thermal graphitization of diamond
In a series of papers, we modeled diamond under femtosecond FEL irradiation in a wide range of photon energies, using XTANT code described above [10,13,66,75,100]. We showed that at the absorbed doses above the damage threshold of ∼0.7-0.75 eV/atom, diamond undergoes a non-thermal solid-to-solid phase transition into graphite phase on an ultrashort timescale of ∼150 fs. An average dose of above 0.7 eV per atom in diamond leads to the excitation of over 1.5% of electrons from the bonding states of the valence to the antibonding states of the conduction band. This is sufficient to trigger the graphitization.
Non-thermal graphitization of diamond proceeds in the following steps [100]: -Initial electronic excitation occurs during the FEL pulse.
In case of X-ray pulses, photoelectrons relax to the bottom of the conduction band within a few up to a fewtens of femtoseconds (depending on the photon energy, see, e.g., [77]) via collisional processes and Auger recombinations of K-shell holes. During this step, lowenergy electrons in the valence and the bottom of the conduction band receive energy from the high-energy electrons, and are starting to exchange it with the lattice. -Electronic excitation triggers a band gap collapse (see, e.g., [100]). It occurs within ∼50 fs (for soft X-rays) after the pulse maximum of the FEL pulse at the time instant when the density of conduction band electrons overcomes the threshold value of ∼1.5%, as mentioned above. This is accompanied by the interatomic sp 3 bonds breaking [100]. For higher photon energies, the electron cascades last longer, thereby delaying all the ensuing processes [75]. -These processes are followed by the atomic relocation (occurring at ∼150-200 fs), which significantly changes the material properties: from insulating diamond to semi-metallic graphite. The electronic density in the conduction band further increases, leading to the final irreversible atomic relocation [100]. In this rearrangement, atoms settle at the new positions corresponding to overdense graphite.
All these stages can be seen in Figure 2 showing atomic snapshots in diamond under irradiation with an FEL pulse with the absorbed dose of 0.8 eV/atom, photon energy of 47.4 eV, and FWHM pulse duration of 52.5 fs. One can clearly observe the formation of the graphite planes. Note that the overdense graphite planes are changing orientations at the edges for a sufficiently large simulation box: in this figure 512 atoms in the box was used, whereas for 216 atoms box formed planes were always perfectly aligned in the simulations. This indicates a presence of multiple nucleation centers of the new phase, as in the case of the homogeneous phase transition [101,102].
These changes of the atomic structure are reflected in the evolution of the diffraction patterns. In Figure 2, the diamond reflections (1 1 1) are present at around 42°, (2 2 0) at around 74°, and (3 1 1) at around 91°at the beginning of the simulation. Then, the peaks (2 2 0) and (3 1 1) become reduced significantly during the graphitization. Later, they almost disappear. In contrast, the peak (1 1 1) is only slightly reduced. The reason for that is that this particular reflection of diamond coincides with a reflection of the overdense graphite, due to their equal density. At later timescales, when expansion of the overdense graphite takes place, the peaks shift towards smaller angles (not shown). The overlap of the diamond and graphite peaks may complicate the analysis of the experimentally recorded diffraction patterns. Luckily, as Figure 2 shows, an additional peak (0 0 2) at around 50°e merges. It corresponds to graphite structure, clearly marking the phase transition. Additional peaks at 40°are observed here, corresponding to the defected structure of the bending graphite planes À those peaks are absent in case of perfectly oriented graphite planes.
Apart from the structure of the non-thermally created graphite phase, we can also analyze its dynamical properties. The vibration spectrum, equation (12), calculated for irradiated diamond is shown in Figure 3a. During the non-thermal graphitization, coherent acoustic phonon excitation takes place (see spectra at times after t ∼ 200 À 300 fs). Such coherent excitation of phonons is in agreement with the model of displacive excitation of coherent phonons (A 1 symmetric phonons) [103], and was known previously for non-thermal melting in different materials [104]. The excited modes correspond to A 2u mode at ∼900 cm -1the characteristic frequencies of D 4 6h symmetric graphite [105,106].
At this time period, the optical phonon modes are inhibited in comparison to the density of states (DOS) characteristic to the equilibrium graphite [107]. As the material tends to a new equilibrium, the energy is transferred to the optical phonons (t = 700 fs), leading to the relative intensity ratio approaching that of the equilibrium phonon spectra. This step of the diamond to graphite transition can be characterized by the onset of optical phonons at ∼1600 cm -1 for the lattice modes of E 2g2 and E 1u symmetry [105,106]. This spectrum already closely resembles known graphite phonon DOS, although at high temperature (these phonons have much larger amplitude than the equilibrium ground state ones) [107]. This is reflected by the autocorrelation function, Figure 3b, where a sudden enhancement of vibrational amplitude ∼200 fs after the pump pulse can be seen. The predicted temperature increase is consistent with the increase of the atomic kinetic temperature estimated from the MD simulation (up to ∼1600 K) [60].
This ultrafast solid-to-solid phase transition À graphitization of diamond under an FEL pulse irradiation À is unambiguously reflected in the evolution of the transient optical properties [10]. Comparison between the experimentally measured transmittance of the optical probe pulse (wavelength of 630 nm) with theoretical predictions is shown in Figure 4. The transmission curve is normalized to the initial transmission of non-irradiated diamond. It exhibits the characteristic multistep process of graphitization described above: (i) initial electronic excitation, (ii) band gap collapse, and (iii) atomic relocation [10]. The remarkable agreement between the theoretical calculations and the recent experiment confirms the transient timescales of the ultrafast graphitization, completed within ∼150-200 fs À to our knowledge, the fastest solid-solid transition observed up to now.

Thermal melting of silicon
As we discussed in the previous papers [67,68,96], silicon can undergo various phase transitions depending on the absorbed dose: a thermal one, occurring for the doses about ∼0.65 eV/atom, and a non-thermal one, for doses above ∼0.9 eV/atom. Thermal melting turns crystalline silicon into low-density liquid by heating the lattice via electron-ion coupling mechanism (non-adiabatic energy exchange discussed in Sect. 2.3). Non-thermal melting quickly leads to high-density liquid phase via an interplay of the lattice heating and non-thermal modification of the interatomic potential, similar to the case of diamond graphitization. The latter case will be studied below in Section 3.2.2.
Irradiation of silicon with X-ray radiation providing an absorbed dose 0.65-0.9 eV/atom induces the following processes. Firstly, high-energy electron cascades deliver photoabsorbed energy to low-energy domain via impact ionizations [67]. Later, low-energy electrons couple to the ions, providing them with energy via non-adiabatic coupling, equations (5)-(4). During this phase, electronic structure of the material is changing. A band gap collapse follows, indicating a transition into a semimetallic phase. Once the lattice is heated enough, the atomic structure transforms into the new phase. Figure 5 shows the evolution of the atomic structure of silicon after irradiation with an FEL pulse of 30 eV photon energy, 60 fs FWHM duration, and 0.72 eV/atom absorbed dose. After approximately ∼300-500 fs, one can see an onset of the melted phase. Diffraction patterns demonstrate that the short-range order is still preserved in the new state, persisting after equilibration of electronic and atomic temperatures on the timescale of ∼1 ps [67,68].
As the stable liquid phase of silicon is a high-density liquid, we expect that LDL phase will densify at longer timescales, unless it is quickly cooled down to be frozen in this low-density state. Presumably, transition to the HDL would take at least a few tens of ps at the considered near-threshold absorbed doses [108]. However, such timescales are too long to access with our present approach, and require dedicated investigations. These aspects are beyond the scope of the current discussion.
To analyze the dynamical properties, the temporal evolution of silicon spectra is shown in Figure 6 for the case of the absorbed dose of 0.72 eV/atom. At the beginning t = 0 fs, the initial phonon spectra are close to the equilibrium spectra of silicon. The peak around 600 cm -1 corresponds to the optical phonons in the vicinity of the G point [110,111].
At later times after ∼300-500 fs, silicon crystal enters a regime of much higher lattice vibration amplitude, which is a signature of melting. The phonon peak at around 600 cm -1 is then shifting to lower frequencies at times around 300-500 fs, indicating phonon softening [112]. It occurs due to the ongoing transition to the low-density liquid state [67]. This observation is consistent with the diffraction patterns shown above in Figure 5. The presence of the local order in the melted phase is also consistent with the experimentally observed low-density liquid phase created in solid silicon by FEL-irradiation [113]. The increase of the atomic temperature can be seen in the gradual increase of the amplitude of the autocorrelation function in Figure 5b.
We can compare the calculated optical properties of irradiated silicon with experimentally available data. Unfortunately, such data exists only for femtosecond optical pulses [109]. However, as we discussed in references [96,114], the electrons after VUV irradiation as well as after an optical one À at the same near-threshold absorbed doses À relax quickly to an equilibrium Fermi distribution. The two cases become nearly identical within a few fs after the exposure, as long as we consider bulk material with periodic boundaries, without any essential contribution of the particles and energy transport. Such a comparison is shown in Figure 7.
This figure also clearly shows that the observed reflectivity overshooting (i.e., the fact that final reflectivity of the irradiated sample is higher than the initial one) is a result of the thermal lattice heating. If we exclude non-adiabatic electron-ion coupling, overshooting does not show up [68]. These results confirm the idea presented earlier in reference [115]: the overshooting effect observed in experimental data [7,93] is a consequence of the ion heating and the resulting band-gap shrinkage. This observation proves that it is essential to step out beyond the BO approximation when modeling the evolution of irradiated solids. Otherwise, important non-adiabatic mechanisms may be missing, and a proper description of the solid evolution can never be achieved.

Coulomb explosion of C 60 crystal layer
In the recent paper, we studied behavior of thin layers of C 60 crystal under irradiation with FELs X-ray pulses [88]. Our predictions indicate that irradiated C 60 crystal disintegrates into single intact fullerenes. The observed fullerene behavior is caused by a Coulomb explosion induced by the charging of fullerene cages. It was confirmed by a comparison of these calculations with the calculations performed assuming all electrons to be confined within the system, which preserves charge neutrality. In the latter case, no breaking of C 60 crystal was predicted [88].
The unbalanced charge is produced due to the photoabsorption by extreme ultraviolet/soft-X-ray (XUV/SXR) laser radiation and impact ionization by photo-electrons and secondary electrons. When the energy of an excited electron is above the work function of C 60 (which is 7.6 eV), the electron can be emitted leaving a positive charge behind. The repulsive forces between neighboring fullerene    [109] reflectivity of silicon at 625 nm probe wavelength under 70.5°i ncidence for the fluence corresponding to 0.72 eV/atom absorbed dose. The data were convolved with 60 fs gaussian probe pulse. Two cases are compared: with electron-ion coupling, and without it (BO approximation). The figure is reproduced from reference [68]. cation radicals then decompose the molecular crystal structure, releasing fullerenes into the vacuum. Figure 8 shows calculated snapshots of the sample at different times following the FLASH irradiation of C 60 crystal (for the photon energy of 92 eV, 30 fs FWHM, absorbed dose of 0.28 eV/atom), reference [88]. This figure shows the scc structure. The fcc structure simulation looks nearly identical. The calculated damage threshold of fcc structure is lower only by ∼10% as expected from the considerations of their cohesive energies [116]. The decomposition of the C 60 layer can clearly be observed. We notice that here, the C 60 cages start to separate from each other on the timescale of ∼2 ps, although the unbalanced charge was created within the first ∼100 fs during the electron cascading. This timescales mismatch is due to a great inertia of the fullerene molecules À such massive objects repeal each other slowly in comparison with the atomic repulsion in case of diamond and silicon phase transitions.
In Figure 8, the diffraction reflections related to the intermolecular distances shift to lower angles, which is a signature of the material expansion. Although the intensity of these peaks is very low, in comparison to the intramolecular peaks at small angles, they clearly indicate a separation of C 60 cages. The fact that the intramolecular peaks stay practically unchanged confirms that intact fullerens are emitted from the irradiated layer, as seen in the atomic snapshots.
We do not perform here an autocorrelation analysis, since the van der Waals potential used [88] was not specifically designed to reproduce the vibrational frequencies but only to yield the correct cohesive energy and structure.
The damage threshold for the molecular Coulomb explosion in thin layer depends on the produced unbalanced charge due to electron emission. This, in turn, depends on the FEL photon attenuation length, and the layer thickness. Thus, we cannot present here a universal damage threshold dose. For the particular photon energy studied, hv = 92 eV, the threshold charge was estimated to be 0.0018 electrons/atom for fcc structure of C 60 crystal (or 0.002 for scc C 60 molecular arrangement). In this case, the absorbed dose of 0.18 eV/atom for fcc structure (or 0.21 eV/atom for scc) produces the corresponding unbalanced charge leading to breaking of bonds between C 60 cages. For comparison, experimental measurement produced the damage threshold of ∼0.15 eV/atom for the identical FEL pulse conditions [88].

Medium fluence 3.2.1 Graphitization of amorphous carbon
Simulated samples of amorphous carbon (a-C) were prepared by quenching of the melted phase starting from underdense diamond (with the density equal to that of the desired a-C, r = 2.6 g/cm 3 ), heated up to the temperatures of a few thousand Kelvins. Quenching was performed by artificially setting atomic velocities to zero every few femtoseconds (similar to the standard zero-temperature MD schemes for quenching [117]) during ∼5 ps until the total energy change in the system became negligible. It indicates that a stable configuration was reached. In a series of simulations multiple initial a-C states were created and then checked for their quality. The most homogeneous and stable one was chosen for the simulation of X-ray irradiation.
Calculations with XTANT for a-C showed that for the considered parameters of FLASH irradiation ( hv = 92 eV, 30 fs FWHM), the spallation threshold is ∼0.85 À 0.9 eV/ atom (cf. the experimental dose of 0.88 eV/atom [88]). Here, the spallation threshold is identified by modeling within Parrinello-Rahman method for NPH ensemble. This allows to trace material expansion and its eventual fragmentation. For an above-spallation dose, the irradiated sample breaks apart into molecular fragments, and the volume of the modeled supercell expands indefinitely (see Sect. 3.3.1). For the below-spallation case, no fragmentation was observed.
In the below-spallation absorbed dose case (<0.85 eV/atom), Figure 9 shows expansion of the irradiated material, which saturates after ∼2.5 ps, with formation of graphite-like structures inside. As atomic snapshots in Figure 9 demonstrate, this process is similar to graphitization, although the formed graphite-like planes are bent and highly defected. This belowspallation expansion reproduces the experimental finding [88].
Analysis of the diffraction patterns for the corresponding atomic snapshots shows emergence of sharp peaks, indicating ordering of the material. This supports the above-mentioned scenario of graphitization of a-C under FEL irradiation.
The good agreement of the calculated spallation threshold with the experimental one confirms the reliability of the model [88]. To the best of our knowledge, there is no time-resolved experimental data available yet on the process of graphitization of amorphous carbon. However, we expect that the timescales of material expansion reported here are probably underestimated, as is typical for Parrinello-Rahman MD simulations.

Thermal and non-thermal melting of silicon
We studied non-thermal melting of silicon under a fs laser pulse irradiation in detail in references [67,68,96]. We performed the silmulations of irradiated silicon for the following FEL parameters: 1.08 eV/atom, photon energy of 30 eV, and FWHM pulse duration 60 fs. This absorbed dose is above the threshold for the non-thermal melting of ∼0.9 eV/atom, calculated in [67]. Note that this dose is significantly lower than the non-thermal melting threshold predicted earlier, e.g., [118], estimated to be ∼2.1 eV/atom. This is due to the fact that in our approach we included not only the changes of the potential energy surface due to excitation of electrons (BO approximation), but also the electron-ion (electronphonon) coupling via non-adiabatic effects [67]. This heating of the lattice by electrons significantly lowers the damage threshold [67].
The damage threshold can also be expressed in terms of the number of excited electrons, which is lowered from ∼9% within BO approximation, to ∼4.5-5% if electronion coupling is included. For comparison, experimental estimation of the electron density threshold is ∼6% [119]. Thus, one can conclude that the interplay of thermal heating with the non-thermal evolution of the potential energy surface plays an important role in the damage formation in silicon [67,68].
Calculations with Parrinello-Rahman MD allowed us to demonstrate that this phase transition proceeds via low-density liquid phase at a picosecond timescales into the ultimate high-density liquid phase, in a good agreement with experiments [67,113]. Atomic snapshots of the material evolution are shown in Figure 10. In the same figure one can also see that the material disorders on the scale of ∼500 fs, with diffraction peaks almost completely disappearing. By the time of 1 ps, only the diffuse scattering background is visible in the powder diffraction.
The evolution of vibrational spectra in silicon after absorbed dose of 1.08 eV/atom is demonstrated in Figure 11. Again, as in the low-dose case discussed above (Sect. 3.1.2), at t = 0 fs, the initial phonon spectra are close to the equilibrium spectra of silicon. After the FEL irradiation, the optical phonon peak (∼600 cm -1 ) completely disappears after ∼500 fs, in agreement with the conclusions drawn from the diffraction patterns. Since the optical phonons correspond to the relative motion of silicon atoms inside the primitive  cell, the disappearance of the relevant peak reflects disordering of the original structure. The non-thermal transition to the high-density liquid state is thus completed [67]. Again, increase of the amplitude in the autocorrelation function in Figure 11b indicate heating of the lattice, however, one can also see that at times after ∼500 fs disorder in the system changes the dynamics qualitatively.
The timescales of damage can be compared with experiments by tracing evolution of the optical properties. In the case of the absorbed dose above the nonthermal melting threshold, the band-gap collapse is induced via non-thermal melting [67,68]; significant heating of the lattice is not necessary for that. For such doses, the overshooting effect allows to extract timescales of the predominant non-thermal melting. We again use the experimental data on the optical pulse irradiation [109], as no FEL-pump time-resolved data exists as of yet to the best of our knowledge. Such a comparison is shown in Figure 12. One can see a very good agreement between the calculated and experimental reflectivities, indicating that the predicted timescales of damage are correct.

High fluence 3.3.1 Spallation of amorphous carbon
For studying disintegration of amorphous carbon, we start with the same initial sample as discussed in Section 3.2.1. Then, we model the evolution of the sample after its irradiation with pulses of different fluences (absorbed doses) with NPH ensemble MD. After running a simulation for 10 ps, we can determine whether the super-cell volume expansion saturates or proceeds continuously and the material breaks apart. From the set of simulations, we can identify the damage threshold dose.
The results show that after irradiation with a dose above ∼0.85-0.9 eV/atom, a-C spallates. This dose corresponds to the peak density of the CB electrons of ∼5.5%, or the electronic temperature of ∼13 000 K.
The spallation proceeds via graphitization on a timescale of a few picoseconds, similar to the abovementioned case from Section 3.2.1. However, in the abovethreshold case, one can clearly see spallated parts of the disintegrating sample at the time instant of 5 ps, see Figure 13.
Emergence of the diffraction peaks from the diffuse scattering pattern clearly shows some ordering in the structure before material disintegration. Then, material breaks apart into a few fragments, but not into atomic species, which is indicative of a spallation regime of the material removal rather than ablation [37]. Note that in contrast to a thin layer of C 60 here we modeled a bulk material with no open surfaces, and thus with no unbalanced charge in the system. The spallation, thus, can be considered as a thermal effect due to atomic heating via electron-ion coupling.

Ablation of silicon
An ablation threshold of silicon was determined within the Parrinello-Rahman MD. From a set of calculations for different deposited doses, the damage threshold was estimated to be ∼2.6 eV/atom. Below this threshold, Si samples demonstrated only (non-thermal) melting, without disintegration. For doses above that threshold, material disintegrated into fragments. This damage threshold dose corresponds to the maximal excited electron density of ∼12À13%, and the electronic temperature of ∼20 kK.
An example of such a simulation is shown in Figure 14 for 3 eV/atom absorbed dose, delivered with an FEL pulse of 92 eV photon energy and 10 fs FWHM. The material disintegration appears to be in the ablation regime, as even individual atoms and small molecular fragments are observed.
Silicon ablation occurs via transient non-thermal melting on sub-ps timescales. This is consistent with the reported above non-thermal melting regime described above (Sect. 3.2.2). The initial changes in the electron structure take place on sub-100 fs timescales, the band gap collapses within ∼40 fs, producing semi-metallic silicon. Already within ∼200-300 fs, Figure 14 shows that powder diffraction reflections vanish into the rising diffuse scattering background. Thus, a complete loss of structural order takes place at this ultrashort timescale. After that, however, new diffraction peaks emerge. They confirm formation of small molecular fragments, emitted from the disintegrating bulk sample. This process is a non-thermal ablation with the contribution from thermal electron-ion heating. We note again that due to inclusion of non-adiabatic effects beyond the BO approximation, our predicted ablation threshold is significantly lower than previously predicted purely non-thermal thresholds within the Fig. 12. Comparison of the calculated and experimental [109] reflectivity of silicon at 625 nm probe wavelength under 70.5°i ncidence for the fluence corresponding to 1.08 eV/atom absorbed dose. The calculated data were convolved with 60 fs gaussian probe pulse. The figure is reproduced from reference [68].
However, again, due to application of the Parrinello-Rahman MD scheme, we expect that the timescales of ablation are underestimated. The experimental studies would probably reveal longer times of silicon ablation.

Ultrahigh fluence: warm dense matter formation
In order to test the limits and capabilities of our developed hybrid approach, we performed a set of simulations at very high deposited doses, leading to WDM formation, reference [121]. The simulations of diamond irradiated with an extremely intense X-ray FEL pulses (at average absorbed doses 18.5-24.9 eV/atom; 6.1 keV photon energy; 5 fs FWHM) performed with XTANT code demonstrate that the atomic structure quickly disorders, on a timescale of a few tens of fs, see atomic snapshots in Figure 15. Diamond transiently undergoes through a stage similar to the graphitization discussed above (Sect. 3.1.1), lasting only for a few fs. From there on, a complete disordering of the sample proceeds.
These effects are visible in the simulated diffraction patterns: the nearest neighbour peak, corresponding to the reflection (2 2 0), disappears faster than the (1 1 1) peak. That is due to the fact that the peak (1 1 1) is present in diamond as well as in overdense graphite formed before material expansion, see Figure 2. Thus, a  This effect is confirmed by the powder diffraction patterns in Figure 15 (e.g. see the snapshot at 20 fs). After that phase, a quick atomic rearrangement follows, leading to the sample disordering at times >20 fs. This is due to the fact that the absorbed dose lies so much above the graphite damage threshold that it triggers atomic disordering in the graphite-like state as soon as it is formed.
This process is clearly of non-thermal nature at its early stage, i.e., until 15-20 fs. The ion temperature increases due to the non-thermal changes of the potential energy surface, and not due to the electron-ion coupling, as confirmed by a comparison with a dedicated simulation within the BO approximation. Such a simulation produced nearly identical result, showing a negligible contribution of the electron-ion coupling at such extremely short timescales [121].
The calculated diffraction peaks in Figure 15 show qualitative agreement with the recent experiment [59]: the intensity of peak (2 2 0) decreases faster than the peak (1 1 1) both in our simulation and in experiment. However, quantitatively, the simulated diffraction peaks vanish significantly faster (by the time of ∼25 fs) than the experimental ones (>80 fs). This observation indicates the limit of validity of our approach reached in this particular case.
The shifts of diffraction peaks (1 1 1) and (2 2 0) measured in experiments were smaller than 0.15° [59]. The lack of any significant shifts of maxima positions of both Bragg peaks indicates that the material expansion due to ablation was insignificant at the experimental timescale of 80 fs. It justifies the usage of MD simulation scheme at a constant volume (V = const, NVE ensemble), which we applied. We note that application of the NPH Parrinello-Rahman scheme in this case resulted in an explosive expansion of the material at the modeled timescales of 30 fs, yielding a clear disagreement with experimental data. This, again, is in line with our argumentation above on the underestimation of the timescales within this MD scheme.

Damage thresholds vs. photon energy
As we have discussed in the course of this paper, in many cases for materials irradiated with fs X-ray pulses, there are universal damage threshold doses. Knowing the photon attenuation length in the material, one can evaluate the corresponding threshold fluence with the inverse of equation (11). It relies upon the assumption of the negligible electron transport, so that the delivered energy is distributed in the material according to the photon penetration profile. This assumption is generally satisfyed in case of X-rays irradiation under normal incidence, since the electron ranges [123] are then significantly shorter than the X-ray attenuation length [124]. It also assumes the thickness of the sample to be significantly larger than the electron range, so that the electron emission from the surfaces could be neglected.
The calculated threshold fluences for various damage channels in allotropes of carbon and in crystalline silicon are presented in Figure 16. They could be used as guidance for preparation of experiments at FEL facilities, which utilize photon energies from few tens of eV to tens of keV. Note, however, that they are yet to be validated by experiments at hard X-rays.

Limitations of the model
Our results from high-fluence regime indicate a few limitations of the XTANT model. We can estimate the limiting fluence or an absorbed dose. Our results discussed in Section 3.4 showed only qualitative agreement with experiment, implying the absorbed dose of ∼18.9 eV/atom is already too high for quantitative model application. On the other hand, comparison with the optical properties of diamond after irradiation with ∼5.4 eV/atom produced a reasonably good agreement with experiment in reference [96]. Thus, the limitation of the model is somewhere in the range of 5.4-18.9 eV/atom. A few reasons for that can be identified.
Firstly, we can expect that the duration of cascading of created high-energy photo-electrons, in reality, may be different from that estimated with the model. This is due to the fact that the XTANT code currently uses cross section for electron impact ionization, calculated for neutral medium with the CDF [60,75]. In case of a strong ionization of the sample by a high radiation dose deposited, the sample neutrality quickly breaks down. To the best of our knowledge, there are no rigorously derived impact ionization cross sections in highly excited solids yet available. However, the recently proposed semiempirical expression based on the experimental data suggests that the cross section of electron impact ionization in dense plasma increases, which may speed up the electron cascading [125].
In reference [121], we showed that (artificially) reducing impact ionization cross sections would be needed to slows down the diamond damage, but even then, it was not sufficient to achieve a quantitative agreement with experiment on WDM formation.
Secondly, created K-shell holes may strongly perturb the electronic band structure, which is not taken into account in the present model. However, in reference [121] the number of K-shell holes was shown to be small enough not to influence significantly the sample evolution.
The most important effect which is not accounted in XTANT seems to be related to the fact that the model relies on the transferable TB parameterization, whose parameters were fitted to the equilibrium configurations of different carbon phases [84]. This approximation misses the effect of the shifts of the electronic energy levels due to the presence of excited electrons. This problem is also known in the plasma community in the context of the ionization potential depression (IPD) [126][127][128]. With the increasing temperature within the heated solid, higher charges appear within the sample (cf. Fig. 5 in Ref. [127]). The energy levels within the band correspondingly move down. Electrons occupying the valence levels below the Fermi level form bonding states, whereas electrons populating the levels in the conduction band above the Fermi level contribute to antibonding states. Thus, lowering of the conduction band levels beyond the Fermi level in the strongly heated diamond may temporally change the bonding from repulsive to attractive. This effect may stabilize diamond on the way to the WDM state and prolong the timescales of WDM formation. Since transient amount of excited electrons in diamond is close to a metal density, perhaps similar effects to the bond hardening observed in metals can be transiently expected in diamond during WDM formation [129,130]. This conjecture may be tested in the future by a dedicated ab initio modeling, such as, e.g., finite-temperature DFT.
Another kind of a limitation of the model is related to the periodic boundary conditions applied. In case of the normal incidence hard X-rays irradiation, when heating is homogeneous within the material down to the depth of a few microns or more, periodic boundary conditions can be applied with a good accuracy as it was done throughout the paper. For lower photon energies, or grazing incidence angles, the photon attenuation length is small and the periodic boundary conditions are inapplicable. Near the plasmon resonance, the attenuation length may be as short as a few nanometers, which induces strong gradients in the absorbed energy, and triggers strong particle and heat transport. These effects hinder the application of the and (right panel) for crystalline siliconmelting into low-density liquid, high-density liquid, and ablation, calculated with XTANT. Experimental data for graphitization of diamond are shown for comparison; points for photon energies below 300 eV are reproduced from reference [13], while the point at 830 eV is from reference [122]. model to the photon energies below about 50 eV. For very low energies, in the ultraviolet or optical regime, additional nonlinear effects such as photoabsorption and inverse Bremsstrahlung must be taken into account [71].
And last but not least of the discussed issues, the electron-ion coupling scheme which we currently apply in XTANT seems to have limitations. In our attempt to apply it to GaAs, the calculated non-adiabatic electronion coupling rate seemed to underestimate the experimental values extracted from the optical probe data [83]. The reason for this discrepancy is likely that GaAs has a very narrow conduction band minimum. For its proper sampling, one would need a large number of k-points [131]. Thus, our electron-ion energy exchange, so far written for the Gamma-point (Eq. (6)), must be extended to include multiple k-points. This will be a topic of a separate study.

Future development
The reliability of the XTANT code, confirmed by its quantitative agreement with various experimental results, proofs the correctness of the approach and demonstrates powerful capabilities of hybrid models. Consequently, one can think of further improvements of the model in order to extend its applicability.
Firstly, in order to enable the usage of the model for other materials, a non-orthogonal TB parameterizations can be implemented. As was mentioned, e.g., in [132], a non-orthogonal TB allows for easier construction of transferable parameterizations. It is already available for many elemental solids [132].
Secondly, in order to go beyond the simple periodic boundary conditions, one can incorporate particle and energy transport effects into the model. One way is to include additional source and sink terms into the electronic and atomic equations. For example, a proper tracing of the electronic transport can be done by extending the MC module and the rate-equation/thermodynamics module to account for spatial dimensions, see, e.g., [133]. In this way, a sample can be discretized into a set of simulation boxes. Each of them would have its specific absorbed dose. Exchange of particles and energies between each other could also be accounted for.
As shown, e.g., in reference [45], at low fluences the thermalization of low-energy electrons may take long times. In order to take this effect into account, one could replace the rate-equations and the TTM with a proper Boltzmann electron-electron collision integrals [45], or a full Boltzmann transport equation [44]. As non-equilibrium electron distribution couples differently to the ions, this may affect the atomic dynamics in the case of long electron thermalization times. Non-equilibrium electron distribution could also affect optical coefficients (cf. Eq. (13)), which may improve their agreement with experiments (e.g. in Fig. 7).
Thirdly, the calculations presented in Sections 3.3.1-3.4 show that the high-fluence regime lies on the border of applicability of the TB method. The latter relies on the ground-state parametrization and therefore should be applied with care under such conditions. In order to reliably describe irradiation at higher fluences, one would need to replace the TB module with a more suitable approach. Since at higher material excitations we expect to reach states far from equilibrium, we require a robust method not relying on the ground-state approximation, which casts doubts on the applicability of the DFT schemes.
In reference [134], an ab initio scheme X MOLECULE for the calculation of the electronic structure in molecular ensembles within the Hartree-Fock (HF) approximation was developed. It includes multiple-hole configurations of molecules or solids formed during XFEL pulses and is able to provide necessary information about the electronic states, also including the influence of core holes on the interatomic forces under strongly non-equilibrium conditions. Incorporation of such a method into XTANT would allow to simulate previously unreachable conditions, such as X-ray-generated warm or even hot dense matter. A disadvantage of such an approach is its high computational costs. In order to decrease them, one could try to use XMOLECULE to obtain classical interatomic force-fields with the parameters adjusted on-the-fly. One of the fitting methods suitable for such an implementation is the force-matching method [135]. Another option would be to create a HF-based TB parameterizations on-the-fly, similarly to the density-functional-based tight binding (DFTB) methods [136,137].

Summary
We proposed a hybrid model specifically designed to treat femtosecond free-electron-laser irradiation of solids. The corresponding code XTANT has been used here for carbon-and silicon-based materials. It provides a good quantitative description of various damage channels: nonthermal graphitization of diamond and amorphous carbon (a-C), thermal spallation of a-C, molecular Coulomb explosion of C 60 crystal, thermal melting of silicon into low-density liquid phase, non-thermal melting of silicon into high-density liquid phase, ablation of silicon. It also gives qualitative insights into WDM formation in diamond.
Wherever available, a comparison with time-resolved experimental data was provided. This was done by extracting experimental observables from XTANT, such as transient transmission and reflectivity coefficients for an optical probe pulse. The reliability of the XTANT code, confirmed by its quantitative agreement with various experimental results, proves the correctness of the approach and demonstrates powerful capabilities of hybrid models in general.
For each studied case, we presented transient atomic snapshots and powder diffraction patterns. The timeresolved diffraction patterns can be used for a comparison with future X-ray-pump X-ray-probe experiments at FELs. Damage thresholds calculated for a wide range of photon energies were also shown.
Finally, we presented a discussion on the limitations of the developed model and outlined directions for its future improvements.