Nonthermal phase transitions in semiconductors induced by a femtosecond extreme ultraviolet laser pulse

In this paper, we present a novel theoretical approach, which allows the study of nonequilibrium dynamics of both electrons and atoms/ions within free-electron laser excited semiconductors at femtosecond time scales. The approach consists of the Monte-Carlo method treating photoabsorption, high-energy-electron and core-hole kinetics and relaxation processes. Low-energy electrons localized within the valence and conduction bands of the target are treated with a temperature equation, including source terms, defined by the exchange of energy and particles with high-energy electrons and atoms. We follow the atomic motion with the molecular dynamics method on the changing potential energy surface. The changes of the potential energy surface and of the electron band structure are calculated at each time step with the help of the tight-binding method. Such a combination of methods enables investigation of nonequilibrium structural changes within materials under extreme ultraviolet (XUV) femtosecond irradiation. Our analysis performed for diamond irradiated with an XUV femtosecond laser pulse predicts for the first time in this wavelength regime the nonthermal phase transition from diamond to graphite. Similar to the case of visible light irradiation, this transition takes place within a few tens of femtoseconds and is caused by changes of the interatomic potential induced by ultrafast electronic excitations. It thus occurs well before the heating stimulated by electron–phonon coupling starts to play a role. This allows us to conclude that this transition is nonthermal and represents a general mechanism of the response of solids to ultrafast electron excitations.

Similar to the case of visible light irradiation, this transition takes place within a few tens of femtoseconds and is caused by changes of the interatomic potential induced by ultrafast electronic excitations. It thus occurs well before the heating stimulated by electron-phonon coupling starts to play a role. This allows us to conclude that this transition is nonthermal and represents a general mechanism of the response of solids to ultrafast electron excitations.
Laser-induced nonequilibrium phase transitions have been intensively studied over the last decade since qualitatively new effects manifested themselves in the experiments [18,19,28,34]. The solid-solid phase transition is of particular interest for basic research as well as for practical applications. A few theoretical studies have been performed on the kinetics of these ultrafast phase transitions [18-20, 22, 23]. The diamond-to-graphite phase transition, following intense visible light irradiation, is taking place within a hundred femtoseconds, before electron-phonon coupling starts to play a role [22]. The graphitization of diamond is an example of the socalled nonthermal phase transition (the most famous realization of which is nonthermal melting [18,19]). It was also observed for other semiconductors and dielectrics [18][19][20][21][22][23]35]. The nonthermal phase transition is caused by the transient changes of the interatomic potential, which stimulate atomic motion and the following transition of the whole electron-atom system into a new transient state [22,23,36]. After the relaxation of the electronic subsystem, and thus of the potential energy surface, the material may relax into a new structural equilibrium configuration, different from the initial one [36]. This process requires only electronic excitation and the corresponding change in the atoms' potential energy surface, without exchange of the kinetic energy between electron and phonon subsystems, and therefore it takes place on the time scale of an atomic vibration (a few tens of femtoseconds) [18].
However, until now it was not clear whether this nonthermal phase transition is only specific for irradiation with visible light. In this work, we perform a theoretical study of the transient response of diamond to femtosecond XUV laser pulse irradiation. For the first time, we show that the nonthermal phase transition of diamond into the graphite phase should occur also at the irradiation at these wavelengths on femtosecond time scales. Within our model we are able to distinguish between (i) the nonthermal phase transition, which is a result of interatomic potential modifications due to electronic excitations, and (ii) the heating of the atoms due to electron-phonon coupling. Our results clearly indicate that nonthermal phase transition is the mechanism of the ultrafast graphitization of diamond.
This paper gives a detailed description of our newly developed hybrid model. It starts with the description of the basics of the model and then discusses the assumptions and approximations made. Let us emphasize that our hybrid model is based on earlier approaches: for atomic motion and the tight-binding method, we used the same method as that in the work by Jeschke et al [22]; for tracing of high-energy electrons the Monte-Carlo (MC) method is used, which is based on the work of [13][14][15]. The combination of MC and the temperature equation relies on the work by Osmani et al [37]. The novelty of our approach consists in combining all these different approaches into one consistent model. This allows us to follow in a computationally efficient way the x-ray irradiation of the material. More technical details of the methods, as well as the specific parameters, are listed in the appendices. Finally, the results of the model, describing the FEL irradiation of diamond in the conditions used in experiments with the Free-Electron Laser in Hamburg (FLASH) [9,10], are shown. A discussion of the results concludes the paper.

The model
XUV photons, when penetrating materials, excite electrons into high-energy states. At the range of intensities relevant for solid-solid laser-induced phase transitions, single-photon absorption is the dominant channel of photoabsorption [15]. A photoexcited electron then loses its energy via collisions with atoms. If the electron energy is sufficiently high, it can excite secondary electrons via inelastic scattering (the secondary impact ionizations). It can also scatter elastically on atoms of the lattice. Both processes are taking place at the femtosecond time scale [38]. Losing its energy, the electron thermalizes and finally falls into the Fermi sea of the valence or conduction band electrons. For the specific case of XUV excited electrons, it has been shown that partial thermalization of the low-energy electrons occurs already within a few femtoseconds (during the laser pulse itself) [12][13][14][15]. The deviation from the thermal equilibrium distribution for this fraction of low-energy electrons is fairly small, even during the transient nonequilibrium state. However, the second fraction of electrons-the high-energy electrons produced by the XUV photons directly, or indirectly via a subsequent relaxation of a deep-shell hole-remain nonthermalized during a few tens or even hundreds of femtoseconds. The number of these electrons is small, compared to the total number of valence electrons, and their energy is high. The characteristic two-component shape of the electron distribution function after intense XUV irradiation [12][13][14][15] is then referred to as the 'bump on hot tail' distribution [17]. That particular shape encourages us to treat those two fractions of electrons with different dedicated (numerically efficient) methods.
Electrons of the valence band act as an attractive potential for ionic cores in semiconductors and dielectrics. Their distribution directly affects the potential energy surface of atoms and, therefore, it influences the structure of the material. This is of importance for the solid-solid phase transitions, because the electrons excited are then also affecting the interatomic potential. These fast radiation-induced changes of the potential energy surface of atoms define the kinetic pathways and can stimulate phase transitions. Consequently, the atomic disorder affects the electronic band structure and the kinetics of the valence and conduction band electrons. Both these effects must be addressed if we want to model the ultrafast material response to XUV irradiation. To treat the electronic and atomic dynamics efficiently, we have developed a dedicated hybrid model whose details are described below.

Overview of the hybrid model
As described before, the radiation-induced electron transient distribution during and after an XUV laser pulse has the shape of the so-called 'bump on hot tail' [17] distribution consisting of: (i) a thermalized bath of low-energy electrons within the valence and conduction bands and (ii) a high-energy nonequilibrium tail [12,14,15]. In the case of a correlated many-body system (such as the low-energy-electron fraction), it is efficient to treat its kinetics with statistical approaches, such as, e.g. the Boltzmann kinetic equation [12,39,40]. If the electronic system is close to equilibrium, then the simpler statistical concept of a temperature can be applied.
For the high-energy nonequilibrium electrons, which are very few at a low laser fluence, methods for tracing individual particles are more efficient. The MC method, simulating the trajectories of individual particles event by event, does not require a numerical grid in momentum or energy space, in contrast to techniques based on Boltzmann equations. The extensive grid necessary for describing energetic particles would make the Boltzmann method inefficient. Therefore, in our model we apply a combination of both the methods: the temperature approach for low-energy electrons and the MC modeling of the high-energy nonthermalized electrons (similar to [37]). A proper tracing of the high-energy electron cascades is important to interpret the experimental data [11].
For tracing the lattice dynamics, we employ a molecular dynamics (MD) technique [22]. The atomic motion is calculated on a potential energy surface which is updated at each time step according to the transient electron distribution. After moving the atoms, we recalculate the electronic band structure using the tight-binding method [24], which accounts for the transient spatial positions of all atoms and the distribution function of the electrons. The numerical scheme is presented in figure 1, and further details are given in the following subsections.

Combined Monte-Carlo and temperature approaches for electrons
As schematically shown in figure 1, we divide all electrons into two fractions: (i) low-energy electrons (with energies below a cut-off energy ∼E cut , which is initially equal to the band gap energy: around 5 eV for diamond) are treated with a temperature model, while (ii) highenergy electrons capable of performing secondary ionization of the valence band electrons (i.e. with energy above the band gap) are followed individually with the MC scheme. At each time step, we calculate the total number of low-energy electrons, N low e , knowing how many electrons got excited by the incoming photons from the laser pulse, N ph , and by the secondary collisions of the high-energy electrons, N imp . Electrons that reach the energy above the cut-off energy are transferred to the high-energy domain; conversely when an electron from the highenergy domain 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 N low (t) denotes the fraction of high-energy electrons that fell into the low-energy domain. The total energy of low-energy electrons, E low e , is calculated from the conservation of total energy among all the electrons (in both domains), and the potential and kinetic energy of atoms. More details will follow in section 2.4.
Knowing the total energy and the number of the low-energy electrons, we can find their temperature and chemical potential. They are calculated from the zeroth and the first moments of the Fermi distribution function by solving the inverse problem [37] N low where the summations include all energy levels, E i , corresponding to the actual band structure of the material that changes in time, and E min is the lowest energy level of the valence band. These energy levels are calculated with the tight-binding method described in section 2.4. The factor 2 in the Fermi-distribution function f e (E i ) accounts for the electron spin; µ is the transient chemical potential of electrons and T e is their temperature. We solve equations (1) for the known values of N low e and E low e at each time step by the bisection method, finding the corresponding values of µ and T e . Note that at the partial occupation of the conduction band (e.g. at nonzero electronic temperatures) the cut-off energy, E cut , must be extended above the band gap value so as to account for the exponential tail of the electron distribution function. This gives an additional constraint for E cut , depending on µ and T e . The choice of E cut made above satisfies this constraint.
For high-energy electrons, we perform an MC tracing of individual particles with eventby-event simulations, similar to those in [14][15][16]. First, we calculate how many photons are absorbed at each time step. We assume that each photon excites only one electron. As discussed earlier, this is a good approximation for photons at XUV energies [15]. The shell from which an electron is excited is chosen according to the relative photon attenuation length for shell photoabsorption [41]. In the case of valence-band photoabsorption, we choose a particular energy level E i randomly.
After the electron is excited from the valence band, it 'disappears' from the low-energyelectron domain, and joins the high-energy fraction. All photons absorbed at this time step undergo the same procedure, creating highly excited electrons. If the energy of the excited electron falls below E cut , this electron does not appear in the high-energy-electron domain any more, but only contributes to the source term in the low-energy-electron domain,Ñ high e , and its total energy, E low e . The excited high-energy electron starts to propagate and to perform secondary scattering events. In order to model this, we calculate the electron mean free path, depending on the current electron energy. We use the method of the generalized complex dielectric function [42][43][44][45]. With this method, the cross-section of inelastic electron scattering, σ i (E e ), is calculated from the complex dielectric function, (ω, q), as with Here the cross-section also depends on the transferred energyhω, and is integrated over the transferred momentum q. The charge e denotes the electron charge and n e = N low e / is the current electron density of the low-energy domain. The volume is the volume of our simulation box, and it is changing in time as will be described in section 2.3. The constant,h, is the Planck constant, and v is the incident electron velocity corresponding to the energy E e . Using the experimental data, the inverse complex dielectric function can be parameterized as, e.g., proposed in [42,45]: The coefficients used, A j , i and E 0 j ( j = 1, . . . , 6), are listed in table A.1 in appendix A. The corresponding free flight distance, l e , of an electron can then be sampled according to the Poisson distribution with the mean free path l e , which is linked to the calculated total cross-section, σ i (E e ), as where the random number, γ , is uniformly distributed in the range of (0,1]. During collisions, the transferred energy, E, is calculated from the differential and the total cross-sections, using another random number γ [13] 5 : An inverse problem must be solved to find the transferred energy E out of the given values of γ , dσ i and σ i (E e ) in equation (5). This can be done again with the bisection method.
After the collision we update the electron energy, E e → E e − E, and calculate the new mean free path for the electron and the time of the next collision. If this time is smaller than our current time step, we simulate the next collision for the electron following the same procedure. In the case of the next collision occurring during the next time step, we save the electron parameters, and follow the next electron [14].
The secondary electron to which the energy, E, is transferred, is then chosen randomly among all the electrons of the valence band, with the ionization potentials E i not higher than the transferred energy [15]. After that, from the energy of the secondary electron, E − E i , we determine whether it belongs to the low-energy domain or to the high-energy domain. In the latter case, we follow the kinetics of this electron in the same manner as in the case of the primary electron. In the same way, we follow all the produced primary and secondary electrons.

Molecular dynamics method for atoms
For modeling the atomic motion we use the classical MD method [46,47], where the forces acting between the atoms are calculated with the help of semi-empirical simulations, namely, the transferable tight-binding method [22,48]. This scheme relies on the Born-Oppenheimer approximation, assuming that the fast moving electrons follow adiabatically the motion of ions. In this way, the electronic subsystem only creates an additional potential for the atomic motion, which is calculated within the tight-binding approximation described in the next section. The adiabatic approximation allows us to solve the classical equation of motion for all atoms of our sample, which are described by spatial coordinates and velocities; however, the potential energy surface is calculated quantum mechanically.
At present, it is hardly possible to apply the tight-binding MD method when the number of ions exceeds a few hundred to a thousand. The typical laser spot radius for the FLASH laser is approximately 10 µm, and the XUV-photon penetration depth may be above 100 nm. 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 its size much smaller than the laser spot, and apply periodic boundary conditions. In this work, we considered three different numbers of atoms contained in the simulation cell: 64, 144 and 216, in order to investigate the effect of cell size on the simulation results. The changing number of atoms had no significant influence on the results, only reducing the statistical fluctuation of the predicted quantities. The predicted physical effects, as well as the related quantitative time scales and the average values of simulation observables remained unchanged with increasing the number of atoms. Note that it is computationally feasible to further increase the number of atoms within the unit cell but in the present case the relatively small number was sufficient to obtain convergent results.
For MD modeling, the periodic boundary conditions are introduced with the Parrinello-Rahman method [49], accounting for the changing geometry of the supercell. That means that the size of the simulation box is assumed to be an additional variable entering the Lagrangian of motion [22,[49][50][51]. This allows us to include the volume fluctuations of the material at a given external pressure, its expansion due to heating, and to properly model a phase transition, which changes the material density. Thus, we write the Lagrangian in the following form [49]: where the first two terms are responsible for the atomic motion; M i is the atomic mass, r i is the spatial coordinate of the atom i; s i are the relative coordinates of the atoms within the simulation box; vectors z run through all neighboring supercells around our simulation box and are responsible for the periodic boundaries; the matrix h of (3 × 3) size is constructed out of the three vectors h = {a, b, c} that span the MD simulation box [36,49,50]; the superscript T denotes the transposition of the matrix; and is the potential energy surface, which depends on the coordinates of all the atoms within the simulation box. The last two terms in equation (6) are responsible for the kinetics of the supercell: the stretching and the shape variation of the simulation box, changing length and orientations of the three vectors h = {a, b, c}. The constant W PR is the effective mass of the supercell [49,50]. With the demand that the volume oscillations in carbon-based materials occur on physically reasonable time scales, it can be chosen as W PR = M i /25.5 [22,36]. The notation Tr(X ) is used for the trace of the matrix X . The external pressure acting on the supercell is denoted as P ext . In our case it is equal to the normal atmospheric pressure. The volume = det(h) is the transient volume of the simulation box, defined as the determinant of the matrix h. As we do not use atomic pair potentials, but a full many-atom potential energy surface, the equations of motion for the relative atomic coordinates s i and for the supercell vector components of h cannot be simplified as in the original work [49], but should be written explicitly as [36] where g = h T h and σ = ∂ /∂h . Details of derivations are presented in appendix C. The potential energy surface, , is calculated within the tight-binding method, described in the next section.
We use the Verlet algorithm for propagating the atomic coordinates and velocities in time, as well as for the supercell coordinates and velocities [46,47]. This ensures a stable numerical scheme.

The transferable tight-binding method applied to carbon
Within the tight-binding method, the atomic Hamiltonian can be written as follows [22,36,48]: It consists of two parts: the attractive part, H TB , which is calculated with the tight-binding method, and the repulsive part describing the repulsion of atomic cores, E rep ({r i j }). This Hamiltonian is used to calculate the electronic energy levels at each time step, E i (band structure), entering equations (1), and to determine the potential energy surface ({r i j }, t), needed as the input to the equations of motion, equations (7). The coefficients, iη , include the on-site energy of atoms [48]; the hopping integrals, t ην i j , are calculated using the Slater and Koster approach [48,52]. In this method, the hopping integrals are considered as fitting parameters to reproduce the band structure of the material [52]. They interpolate between the accurate energies at selected k points. These energies are known from ab initio calculations. The method takes into account only the s, p x , p y and p z valence orbitals of carbon, but this is enough to reproduce the band structure with the required accuracy [22,48].
In order to study structural transformations within carbon materials, such as the diamondto-graphite phase transition, the hopping integrals have to be 'transferable' between the different structures. This means that these parameters have to be functions of the atomic coordinates of all atoms in our simulation box. Therefore, the following scheme is used to calculate the transferable hopping integrals-first, following the Koster-Slater procedure, we separate the angular and the distance parametrizations as follows [52]: where the distance-dependent functions V ξ (with ξ = {ssσ, spσ, ppσ, ppπ }) are now introduced [48]. The other hopping terms t can be written analogously [36]. Their dependence on the interatomic distance for each pair of atoms within the simulation box and from the neighboring supercells is described by some fitting functions. Also the repulsive potential energy is fitted. Following [48], these quantities are parameterized as follows: All the corresponding coefficients and further details of equation (10) are given in appendix B. The electronic energy level spectrum E i is then obtained by diagonalization of the Hamiltonian matrix H TB , constructed from equations (8)- (10). The potential energy surface can then be calculated as Its corresponding derivatives for equation (7) and the coefficients for these derivatives can be found in [36] and are also presented in appendix C. Let us emphasize that the transient electron distribution function, f e (E i , t), calculated in equations (1), enters in equation (11). Thus, the time-dependent electronic distribution affects the atomic motion. In this manner, the problem becomes self-consistent, closing the circle of the interconnections between the electrons and atoms. It must also be noted that the distribution function of electrons is split into two parts, as discussed in section 2.2. Within our approach, removing electrons into high-energy states reduces the number of low-energy electrons within the bands. This decreases the attractive part of the inter-atomic potential, creating an effective transient charge non-neutrality. The changed inter-atomic potential causes atom displacement which, in the case of high irradiation, may even lead to a Coulomb explosion of the atomic system. Other modifications of the band structure due to the presence of high-energy electrons are not taken into account. This is an approximate way to account for these effects, but, going beyond the scheme of slow energy equilibration between the high-energy electrons and the bands. We assume that this approximation works fine in the case of low radiation fluences, when the number of excited high-energy electrons is small in comparison with the number of valence electrons. However, it is not valid in the case of higher fluences. This is a limitation of the current model. Similarly, in the case of core-hole excitations, we also assume that the created core holes do not significantly affect the structure of the valence and conduction bands and the inter-atomic potential if the density of core holes is much lower than the total electron density. In the case of higher fluences, i.e. higher density of core holes, this assumption breaks down [8]. At a large number of highly excited electrons or holes, the tight-binding method with predetermined tight-binding parameters cannot be applied to describe the band structure, as the tight-binding parameters are adjusted to the ground state, and they are no longer able to reproduce the electronic structure when it is extremely excited. In both cases the correct description of perturbed band structure would then require ab initio methods (e.g. density functional theory and time-dependent density functional theory).

Results and discussions
We have modeled the irradiation of diamond with a femtosecond laser of τ = 10 fs duration at the full-width at half-maximum of a Gaussian temporal profile, and a photon energy of 92 eV. The chosen parameters correspond to parameters of the free-electron laser FLASH used in experiments with solids [1,6,7,9,10]. First, we place all the atoms at positions corresponding to the diamond structure, and assign to them initial velocities according to a thermal distribution (assuming room temperature). We then let the atoms relax for a number of time steps in order to reach an equilibrium distribution. We also place electrons according to the Fermi distribution at 300 K. For the case of a wide-band-gap semiconductor, this implies full occupation of the valence band and an empty conduction band. The time step used in our simulation is 0.1 fs. This guarantees energy conservation with an accuracy of ∼0.1% with the Verlet algorithm used for the simulation of atomic motion [36]. Figures 2(a) and (b) show the calculated total and potential energies in the diamond target irradiated with the laser pulse. The maximum of the Gaussian pulse is at 30 fs. After the irradiation, the total energy of the system is conserved. In order to check it, we performed longer calculations for up to 50 ps time for unexcited and weakly excited systems, and observed no changes in the energy conservation properties: the energy was conserved with the same accuracy of 0.1%. As expected, the figures show a direct increase of the total energy during the laser pulse (from ∼20 to ∼40 fs). The potential energy of the atoms is fluctuating, being interchanged with the kinetic energy of the atoms. The total energy of the atoms during and shortly after the pulse is slightly lower than the total energy in the system; this reflects the fact that a part of the energy is transiently stored in the high-energy electrons and then is relaxing back via secondary electron scatterings. These high-energy electrons transiently have energies between the bottom of the conduction band and 92 eV, but are relaxing fast to the low-energy states (see below, figure 3), returning their energy to the potential energy of atoms. As was explained in the previous section, the atomic potential energy contains the repulsive contribution of ionic core-core interaction, and the attractive contribution of low-energy electrons; thus, the low-energy-electron energy is already included in the potential energy of atoms. Therefore, the sum of high-energy electrons' energy and the total atomic energy is the conserved quantity in our simulation.
As can be observed in figure 3, the high-energy electrons are relaxing fast to the low-energy state, inducing secondary electron cascades until approximately 100-150 fs. The majority of the high-energy electrons fall into the low-energy domain even within some 50-70 fs. So, transiently, the contribution of the electrons to the interatomic attractive potential is reduced, but is then restored, still before the phase transition starts.
One can see a clear difference between the 'below damage threshold' and 'above damage threshold' cases in figures 2(a) and (b), respectively: the potential energy fluctuates much more strongly in the latter case. This occurs because in this case the atoms are already undergoing a phase transition from diamond to graphite. The predicted phase transition can be followed in figure 4, where snapshots of atomic positions are recorded at different times for the case of higher energy deposition (0.95 eV per atom). In an unexcited diamond, atoms are located at their equilibrium positions at t = 0 fs. Later, they absorb the energy from the laser pulse (see the snapshots at 20 and 40 fs) and start to move extensively. At 60 fs, the structural rearrangement already starts, and one can see the breaking diamond bonds. It is completed by the time of 80-100 fs (similar time scales were recently reported in [11]). The atomic structure is then graphite-like with the characteristic parallel planes. The plotted pair correlation function in figure 5 also indicates the occurrence of a structural transition from the diamond to the graphite phase. At these ultrashort times, the newly formed graphite is not in its natural state, as its density is still diamond-like, i.e. higher than the normal graphite density. On later time scales, the supercell is relaxing with its increasing size, and the natural graphite density is then reached. Figure 5 shows four snapshots of the calculated pair correlation function of atoms: before irradiation (0 fs); during irradiation (50 fs); shortly after the pulse, when the phase transition to the graphite state is in progress (75 fs); and when the phase transition is already completed (100 fs). At the first snapshot, the peak at 1.53 Å corresponds to the nearest-neighbor distance of atoms in solid diamond at room temperature. On later time scales, this peak is shifting to the position of 1.41 Å, which is the nearest-neighbor distance of atoms in graphite. Such a shift is a clear indicator of the phase transition. This transition is stable and irreversible, as the graphite phase is more stable than the diamond one [54]. In our case we have checked the stability of the new configuration by running the code for longer time scales (several ps). Figures 6(a) and (b) show the transient behavior of the predicted instantaneous temperatures of low-energy electrons and atoms. As expected, in both cases, below the damage threshold and above the damage threshold, the temperatures fluctuate, following the oscillations of kinetic energy. For the below-damage fluence, the lattice temperature almost does not change on a time scale of half a ps, similar to the results reported in [55]. The high-energy electron temperature increases during the pulse and then oscillates around some average value until the diffusive relaxation takes place on a picosecond time scale [22]. Such a relaxation is, however, not included in our model. In contrast, for the above-damage fluence, we can see a sudden jump in the lattice temperature after the pulse is finished. The average lattice temperature then increases further. This is correlated with the average decrease of the potential energy of the atoms ( figure 2(b)). Interestingly, although the final temperature of the atoms, reached at ∼80 fs, corresponds to the graphitization temperature of diamond [54], its increase coincides with the phase transition itself, as can be seen from the comparison of the time scales in figures 2(b) and 5. This implies that the temperature increase is not the cause of the phase transition, but its consequence. As atoms move to their new equilibrium positions on the changed potential energy surface (the nonthermal phase transition process), they gain kinetic energy and, thus, temperature. In this way it is possible for the atoms to undergo the structural phase transition much faster than at the normally expected pico-to nano-second time scales [36].
Our simulations were performed for different values of absorbed energies, and the results were analyzed in order to determine the threshold fluence for the phase transition. We found that the energy, which is enough to induce graphitization, is ∼0.69 eV per atom for 92 eV photons. Note that the nonthermal phase transition of diamond to graphite predicted with our model is induced by the electronic excitation, when the density of electrons in the conduction band (both in low-and high-energy domains) overcomes a threshold of ∼1.5% of the valence electron density. A preliminary comparison with experimental data obtained at FLASH [9,10] shows the agreement of the calculated damage threshold for diamond with the data. A detailed analysis is in progress.
It is also interesting to note that there can be a delay between the laser pulse and the phase transition for near-damage-threshold fluences. For absorbed fluences from 0.69 to approximately 0.7 eV per atom, the calculations indicate the start of the phase transition after an arbitrary delay from a few to a few hundred femtoseconds (the lower the fluence, the longer the delay). This indicates the stochastic origin of the process, showing that the atomic system might require some time to realize its kinetic relaxation pathway. This topic of near-threshold behavior requires additional investigations, which are, however, beyond the scope of the present paper.
The predicted nonthermal ultrafast phase transition under XUV irradiation of diamond is very similar to that observed for visible light irradiation [22]. The estimated damage threshold agrees well with the graphitization threshold for the visible light irradiation [22,36]. Although the transient electron behavior is different, the atomic kinetics on the changing potential energy surface manifests itself in a similar manner: the irradiated diamond undergoes a phase transition to the graphite phase on a femtosecond time scale. This allows us to conclude that the nonthermal phase transition can be a general mechanism of a semiconductor or a dielectric response to ultrafast electron excitation. As the nonthermal phase transition can be monitored directly with femtosecond resolution by time-resolved measurements of atomic kinetics [11,18,19], a dedicated experiment with FELs should verify our predictions.
Note also that transferable tight-binding parameters for different materials (semiconductors and dielectrics) can be found in the literature (see, e.g., [36,[56][57][58][59][60]), as well as electron crosssections of scattering (e.g. [43,45,61,62]). Thus, the extension of our method to different materials is straightforward. In this way a wide range of applications can be addressed with the developed model.

Conclusions
In this work, we have described our newly developed hybrid model of semiconductor response to XUV femtosecond laser pulses. The model consists of interlinked modules, dedicated to simulating electron and lattice dynamics. The transient electronic state is treated with the two methods: highly excited electrons are traced with the MC scheme of event-by-event individual particle simulation; the low-energy electrons, populating the valence and the conduction bands, are described by an electronic temperature, with electrons populating a transient electronic band structure. Both domains are, of course, interconnected, allowing the electrons to interchange between low-and high-energy sectors via excitation and relaxation mechanisms.
The lattice dynamics is evaluated with the Parrinello-Rahman MD method within a supercell, where the potential energy surface is updated at each time step of the simulation with a transferable tight-binding method. This method naturally accounts for the transient electronic band structure, which follows the current electron distribution and the changes of the atomic positions of all the atoms within the simulation box.
We applied our hybrid model to solid diamond, irradiated with a laser pulse of 10 fs duration and a photon energy of 92 eV. This corresponds to the typical FLASH laser parameters used in solid state and plasma experiments. The model has proved to be robust and numerically stable. The results of the modeling demonstrated that the ultrafast nonthermal diamond graphitization is completed within a few tens of femtoseconds after the exposure. This transition is followed by the increase of the lattice temperature up to the 'graphitization' temperature. The time scale of this transition is comparable with the electronic relaxation time. That allows us to conclude that nonthermal phase transitions are a general mechanism of the response of solids to ultrafast electronic excitations. The calculated graphitization threshold of diamond for a femtosecond laser pulse with 92 eV photon energy is found to be ∼0.69 eV per atom, which corresponds to the excitation of ∼1.5% of valence electrons into the conduction band.  atom i. The off-diagonal blocks, connecting atoms i and j, contain the hopping integrals [36]:   [48].