Multistep transition of diamond to warm dense matter state revealed by femtosecond X-ray diffraction

Diamond bulk irradiated with a free-electron laser pulse of 6100 eV photon energy, 5 fs duration, at the ~19–25 eV/atom absorbed doses, is studied theoretically on its way to warm dense matter state. Simulations with our hybrid code XTANT show disordering on sub-100 fs timescale, with the diffraction peak (220) vanishing faster than the peak (111). The warm dense matter formation proceeds as a nonthermal damage of diamond with the band gap collapse triggering atomic disordering. Short-living graphite-like state is identified during a few femtoseconds between the disappearance of (220) peak and the disappearance of (111) peak. The results obtained are compared with the data from the recent experiment at SACLA, showing qualitative agreement. Challenges remaining for the accurate modeling of the transition of solids to warm dense matter state and proposals for supplementary measurements are discussed in detail.

Over the years, extensive theoretical studies of many-body systems lead to a development of two complementary approaches towards their description: (i) In case when the average kinetic energy of constituent particles is much larger than the average potential energy between them and than the Fermi energy of the system, the system is in the classical, ideal plasma state 1 . This state can be well described within semi-classical approximations 2,3 . (ii) In the opposite case, when the Fermi energy and the average potential energy are much larger than the average kinetic energy of the particles, the system is in the condensed matter state. For its description, the quantum methods including the description of interparticle interactions should be applied. In particular, the density-functional theory (DFT) methods have been extensively developed for the last few decades 4,5 .
At moderate temperatures and densities, when the two regimes (i) and (ii) meet, a new state, the so-called warm dense matter (WDM) emerges [6][7][8] . Neither of the above mentioned approaches can be rigorously applied in this specific regime, as the potential energy of interaction among the particles is of the same order as their kinetic energies. No rigorous theory or model has been proposed so far to describe properties and behavior of the exotic state of matter, in particular, the non-equilibrium transition from condensed matter into the plasma state. This is one of the most challenging problems bridging the solid-state and the plasma communities 6 .
Yet, the WDM state is common in the Universe: it exists in the inner core of large planets (such as Jupiter and Saturn), in white dwarf stars, and, supposedly, on the surface of neutron stars 6,7 . In the laboratories, it is produced as a transient state, following a deposition of high energy dose into a solid target 6,9 . Such energy deposition is often performed through the irradiation of the target sample with high-intensity ultrashort laser pulses, e.g., strongly focused X-ray pulses from the free-electron lasers (FEL) such as FLASH 10 , LCLS 11 , SACLA 12 , SwissFEL 13 , and European XFEL 14 . The modern free-electron lasers provide femtosecond intense pulses of X-ray photons, sufficiently bright to create a WDM state in a single shot [10][11][12] . These new experimental opportunities create a strong demand for a theoretical support and numerical modeling to describe experimental results, to suggest and guide new experiments, and ultimately to understand the fundamental physics of the warm dense matter.
In an ultrafast strongly driven matter, intrinsic processes occur at femtosecond timescales [15][16][17] . Excited electrons redistribute their energy among themselves and, later, also transfer it to the ions. The latter process may lead to material modifications, and to formation of transient nonequilibrium plasma 18 .
In the recent experiment at the free-electron laser facility SACLA, diamond was studied with an X-ray X-ray pump-probe scheme 19 . The experiment utilized the modern capability of the FELs: two color mode, enabling almost homogeneous deposition of an extremely high fluence into the sample with the pump pulse, and a measurement of the X-ray diffraction patterns with the probe pulse. The probe pulse arrived after a certain femtosecond delay from the pump pulse. The photon energy of the pump pulse was 6100 eV (2.03 Å), whereas the probe pulse had photon energy of 5900 eV (2.10 Å). Both pulse durations were 5 fs 19 .
The pump fluences achieved in the experiment were 2.3 × 10 4 J/cm 2 , 2.7 × 10 4 J/cm 2 , and 3.1 × 10 4 J/cm 2 . For 6.1 keV X-rays, such fluences result in average absorbed doses of 18.5 eV/atom, 21.7 eV/atom, and 24.9 eV/atom respectively, deposited in diamond within the photon attenuation length of 279 μm 20 , much larger than the sample thickness. After such energy deposition, the material turns into nonequilibrium warm dense matter state 16,21 .
Monitoring X-ray diffraction patterns of irradiated diamond, Inoue et al. observed that the integrated probe diffraction intensity of the (220) reflection decreased faster than that of the reflection (111). They both significantly decreased within 80 fs after the maximum of the pump pulse, which was the longest delay available between the pump and the probe pulses in that experiment 19 .
The interpretation of the results turned out to be challenging. The decreasing intensity of the diffraction peaks could be analyzed in ref. 19 only in the framework of the Debye-Waller factors, assuming thermal atomic oscillations. A possibility of nonthermal material destabilization due to electronic excitation was also discussed by the authors. Such nonthermal phase transition results from a change in the potential energy surface of atoms due to the electronic excitations. A well known example of such transition is nonthermal melting, see, e.g. 22,23 . In contrast, thermal phase transition is due to the kinetic energy exchange between hot electrons and atoms (e.g., electron-phonon coupling) through non-adiabatic coupling between the two subsystems.
However, the contribution of those mechanisms was addressed in ref. 19 only on the level of a general discussion, not supported by any quantitative theoretical analysis. It remained then unclear which physical mechanism did lead to the experimental observation, and why one of the peaks decreased faster than the other. These questions are of the main interest for the current study which uses a modeling tool that can address transient processes occurring on sub-ps timescale.
We apply the recently developed hybrid code XTANT (X-ray-induced Thermal And Nonthermal Transitions 24,25 ) to model X-ray irradiated diamond at the experimental conditions used in ref. 19 . The model consists of a few modules dedicated to simulate various processes induced by the incoming X-ray FEL radiation: Newton equations for nuclei, with the interaction potential evaluated within the TB module. (c) Electron occupation numbers within the TB-based band structure are assumed to follow Fermi-Dirac distribution with a transient temperature and chemical potential evolving in time. The electron temperature is changing due to interaction of band electrons with X-rays, high-energy electrons excited by X-ray photons, impact ionizations, or Auger-decays; or due to their nonadiabatic interaction with nuclei (through electron-phonon, or more generally, electron-ion scattering). (d) Non-equilibrium part of the electron distribution at high-energies is treated with a classical event-by-event Monte Carlo simulation. It stochastically models X-ray induced photoelectron emission from K-shell or from the valence band, Auger decays, and the scattering of high-energy electrons. (e) Electron-ion energy exchange, mentioned above, is calculated with a nonadiabatic approach, in which matrix elements are calculated as an overlap of TB wave-functions plugged into the Boltzmann collision integral.
More details are given in the Methods section. This combined approach allows to treat predominant processes within an X-ray FEL irradiated samples, including their non-equilibrium evolution stage, possible nonthermal effects and phase transitions 26 . With this approach, we will now address the formation of warm dense matter from an X-ray irradiated diamond under the conditions of the experiment 19 .

Results
The simulation of diamond irradiated with an X-ray FEL pulse (Gaussian pulse of 5 fs FWHM duration; 6.1 keV photon energy; average absorbed dose 24.9 eV/atom) performed with XTANT code demonstrates that the atomic structure quickly disorders, on a timescale of a few tens of femtoseconds, see atomic snapshots in Fig. 1. This damage proceeds as a sequence of steps, similar to the graphitization of diamond occurring at lower fluences (reported earlier in, e.g., ref. 27 and confirmed in ref. 28 ), however, on much shorter time-scales. During the simulation, firstly, valence electrons are excited into the conduction band during and for some time after the pulse. The number of the excited electrons is high, overcoming the damage threshold value (~1.5%) 24 already early during the exposure to the pump pulse 27 . This leads to the band gap collapse at the time ~15 fs. This evolution stage is documented in more detail in the Methods section.
Diamond transiently undergoes the graphitization-like stage, lasting only for a few femtoseconds. It is indicated by a transient presence of the graphite nearest-neighbor peak in the pair correlation function in Fig. 1, and a small peak in the diffraction patterns ( Fig. 2) at 20 fs and 23 fs. The intensities of diffraction peaks from our simulated atomic snapshots were extracted using an open software reciprOgraph 29 . Zoom into the graphite-like peak and an extended discussion on its occurrence are presented in the Supplementary Information. Further material disordering continues from this graphite-like sp 2 -carbon phase. A quick atomic rearrangement follows which leads to the sample disordering at times >20 fs, see to the material is sufficiently high not only to convert diamond into a graphite-like phase, but also to damage the forming phase on femtosecond timescale. In particular, the previously reported ablation dose in graphite (~3.3 eV/atom 30 ) is overreached already during the beginning of the FEL pulse, triggering a rapid damage of the graphite-like state while it is still forming.
This process is clearly of nonthermal nature at its early stage, i.e., until 15-20 fs. The temperatures of electrons and ions are shown in Fig. 3. It is important to emphasize that the ion temperature increases here due to nonthermal changes of the potential energy surface, and not due to electron-ion (electron-phonon) energy exchange.
To justify this statement, we performed a simulation within the Born-Oppenheimer approximation that naturally excludes nonadiabatic electron-ion coupling. Figure 4 compares the full calculations with the BO approximation. The BO-results are nearly identical to the case with electron-ion coupling included. This confirms a  negligible contribution of electron-ion coupling at such extremely short timescales, and proves that the ion temperature increase and the material damage are of nonthermal origin.
Electronic ensemble reaches temperatures of 7 to 9 eV at the time of ~20 fs after the maximum of the pump pulse, due to photoabsorption and further secondary high-energy electron scattering. Such an extreme level of electronic excitation severely modifies potential energy surface that keeps ion lattice stable 24 . Ions remain cold until ~15 fs. Their temperature starts to raise quickly when the modified interatomic potential pushes them out from their former equilibrium positions, at the same time increasing their temperature up to ~3 eV. The timescale of the ion temperature increase is in agreement with the estimates of ion displacements obtained with the Debye-Waller fit in ref. 19 (Fig. 4 therein). The predicted displacements start to increase fast at the times above 20 fs, in agreement with our predictions on ion temperature from Fig. 3. This way, a warm dense matter state is formed on a ten-femtosecond timescale.
To support the proposed mechanism of diamond's structural disordering, we will now compare our simulation results to the experimental data. In ref. 19 , the authors presented detailed data on the decrease of (111) and (220) diffraction peak intensities in diamond after its irradiation with an X-ray FEL pulse.
The calculated results in Fig. 5 show qualitative agreement with the experiment. The intensity of peak (220) decreases faster than the peak (111) both in our simulation and in the experiment. However, quantitatively, the simulated diffraction peaks vanish significantly faster (by the time of ~25 fs) than the experimental ones (>80 fs). The reasons for this discrepancy will be discussed later.
To interpret those findings, we first discuss the features of the diffraction peaks of equilibrium diamond, overdense and equilibrium graphite. Figure 6 shows the diffraction peaks at the probe wavelength of 2.1 Å  (5900 eV). The corresponding atomic structures are shown in the insets on the right. One can note that the intensity of the diamond reflection (111) is only slightly reduced in the overdense graphite (with the density equal to that of diamond, as transiently formed during the graphitization-like process before material expansion), whereas the intensity of the (220) reflection is reduced by half. The reason is that the changes of the (220)-peak indicate a change of the nearest neighbor distance in diamond which is a primary indication of a progressing structural transition. Thus, the strong reduction of the (220)-peak, with the (111)-peak still present, indicates a progressing graphitization-like process.
Indeed, in the experiment the intensity of the reflection (220) decreases faster than of the reflection (111). This indicates that the damage proceeds via a transient graphitization-like phase, similarly to that described in ref. 26 . After that, the graphite-like structure disorders, which leads to a disappearance of the (111) peak. The simulation results indicate that the femtosecond time delay between the disappearance of the peak (220) and the peak (111) corresponds to the duration of the transient graphite-like sp 2 -carbon state. This prediction can be tested in the future provided that the temporal resolution of pump-probe experiments improves enough to measure such short-living states.
Additionally, ref. 19 presents angular positions of the diffraction peaks. Only slight shifts (<0.15 o ) of the maxima positions for both (111) and (220) peaks towards smaller angles were observed at the timescale of 80 fs. Although the absorbed dose exceeded: (i) the threshold for the graphitization of diamond, ~0.7 eV/atom, corresponding to 1.5% of excited valence electrons, (ii) the threshold for damage of graphite, ~2 eV/atom, corresponding to ~9% of excited electrons 30 , and (iii) the ablation threshold of carbon, ~3.3 eV/atom, corresponding to ~12% of excited valence electrons 30 , estimated with XTANT, 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 timescale of 80 fs. It justifies the usage of MD simulation scheme at a constant volume (V = const, NVE ensemble) which we applied here.

Discussion
Our results achieve only qualitative agreement with experiment. In particular, the simulated timescales of WDM formation are much shorter than those observed in experiment 19 . This is due to a few limitations of the XTANT model, which manifest here due to the extreme irradiation conditions of the simulated experiment. Below we discuss in detail the reasons for the discrepancy observed.
Firstly, it is to be expected that the created high-energy photo-electrons are, in reality, cascading for longer times than estimated with the model. The XTANT code currently uses cross section for electron impact ionization, calculated for neutral diamond from the complex dielectric function 17,24 . The calculations with different impact ionization cross sections for neutral sample produce very similar cascading times 31 . However, here, due to a strong ionization of the sample by a high radiation dose deposited, the sample neutrality quickly breaks down. There are no rigorously derived impact ionization cross sections in highly excited solids yet available. Nevertheless, it can be expected that the cross sections calculated within highly excited solids should be smaller than those obtained for neutral samples, as the ionized nuclei attract the valence electrons stronger, and their binding energies would increase as already shown for plasmas, e.g., in ref. 32 . Including the proper cross sections would then lower the impact ionization rate and slow down the progress of impact ionization, when compared to the current prediction with our code. This would, in turn, slow down diamond damage, pushing our predictions towards the experimentally observed longer disordering timescales.
We performed a test case study, in which we artificially reduced the electron inelastic scattering cross sections by a factor of two. This change slightly prolonged the electron cascading times, and delayed the onset of the damage by a few femtoseconds (see Fig. 4, the curve marked as σ/2). However, it did not significantly reduce the temporal discrepancy between our predictions and the experiment.
At this point, we can also evaluate the impact of the incoming X-ray energy on the graphitization-like progress. Figure 7 shows the decrease of intensity peaks as a function of the X-ray photon energy. Reducing the X-ray photon energy reduces the formation time of secondary electron cascades, as the corresponding energy of the photoelectron triggering the cascade is significantly lower. The transient graphitization-like process would then still occur, however, it would start much earlier than in the reference case (X-ray energy ~6.1 keV). In the opposite case of much higher photon energies, the graphitization-like onset would be delayed. These observations could be used when planning any future measurements to unveil the details of this transient process.
Secondly, the majority of the emitted high-energy photoelectrons leave positively charged ions with K-shell holes. The presence of K-shell holes may affect the sample dynamics. Two effects can contribute: modification of the electronic band structure and the additional potential from an effective charge non-neutrality. To estimate their potential impact, we calculated the densities of K-shell holes, Fig. 8. The calculation shows that the number of K-shell holes per atom is only on the order of 0.3-0.4% at most for the considered experimental fluences, in agreement with the estimates from ref. 19 . Thus, one can assume that K-shell holes do not influence the atomic dynamics significantly -in this particular case.
The number of high-energy electrons (with energies above 10 eV counted from the bottom of the conduction band) is higher than of the core holes, Fig. 8, ~18-25% per atom. Thus, we analyzed their influence on the atomic motion, performing an additional simulation with the incoming soft X-ray photons of energies ~50 eV, at the pulse parameters selected to yield exactly the same average absorbed dose per atom as the doses from the SACLA experiment. Thereby, we enforced the photoexcited and impact-ionized electrons to stay within the low energy regime only, and relax quickly to the bottom of the conduction band. In such a simulation, charge neutrality restores fast. The results of this simulation show that the electrons at the bottom of the conduction band are speeding up the atomic dynamics, rather than slowing it down -the diamond disordering in this case occurs already at ~15 fs (Fig. 7). Thus, we can confirm that the numerous presence of high-energy electrons seems not to be the reason for the faster atomic disordering predicted by our model, when compared to the experimental results.
Thirdly, the XTANT model relies on the transferable tight binding parameterization, whose parameters were fitted to the equilibrium configurations of different carbon phases 33 . 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) 32,[34][35][36] . With the increasing temperature within the heated solid, higher charges appear within the sample (cf. Fig. 5 in ref. 32 ). The energy levels within the band correspondingly move down. This, again, supports the expected lowering of the impact ionization rate as the ionized nuclei attract the valence electrons stronger, and the corresponding binding energies increase. Electrons occupying the valence levels below the Fermi level form attractive bonds, whereas electrons populating the levels in the conduction band above the Fermi level contribute to repulsive bonds. 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 warm dense matter state and prolong the timescales of WDM formation. Since the density of transiently excited electrons in diamond is close to the conduction electron density in metals (see Fig. 9 below), the effects similar to the bond hardening observed in metals may be expected in diamond underway to warm dense matter state 37 .
In order to estimate the regime of applicability of the tight binding method, let us recall the observation made about the weak influence of high-energy electrons on the dynamics within the irradiated sample. As long as the fast photoexcited electrons stay within the high energy regime, i.e., are weakly coupled, they perturb weakly the electronic structure of the solids. The tight binding method for equilibrium can still be then applied with a good accuracy, if the number of high-energy electrons is not too high to significantly affect the charge neutrality. This time span corresponds approximately to the cascading time of high-energy electrons, here ~15-20 fs (Fig. 9). As soon as the electrons start to fill in the low-energy domain by secondary impact ionizations, the equilibrium tight-binding approximation breaks down due to the presence of so many excited electrons in the low-energy domain.
Indepedently of these considerations, it would be good to check the influence of the specific TB parametrization 33 itself on our predictions, replacing it with another one. However, a construction of another fully transferable TB parametrization or an alternative ab-initio model is a nontrivial task, which is far beyond the scope of  Finally, in these simulations, the effect of the probe pulse on the sample was neglected. Although the fluence of the probe pulse is twice as high as the fluence of the pump pulse, the probe pulse takes only 5 fs FWHM. This timescale of the diffraction measurement is too short to record any atomic relocations induced by the probe pulse and reflected by Bragg peaks 38 . We can, therefore, neglect its influence on the simulated results.
In conclusion, with our simulations, we have followed the multistep pathway of a diamond crystal exposed to intense hard X-ray FEL pulses (6.1 keV photon energy, 5 fs FWHM duration) to the warm dense matter state. Our simulations indicate a few evolution stages of the irradiated sample. Diamond transition starts with the strong electronic excitation, leading to nonthermal band gap collapse and the subsequent transient formation of a graphite-like sp 2 -carbon phase, followed by the sample disordering -within ~20 fs in total. The results are compared with the recent experimental data on diffraction peaks of diamond recorded at SACLA with the X-ray X-ray pump-probe. In particular, our simulations confirm that the intensity of (220) reflection decreases faster than the intensity of (111) reflection, in a qualitative agreement with experiment. This effect is caused by the transient graphitization-like process of a femtosecond duration, occurring underway to warm dense matter state. Our finding then allows to explain the puzzling experimental observation on the disappearance of the diffraction peaks at different timescales.
Quantitative discrepancy between the theoretically estimated timescales of diamond disordering and the experimental ones can be attributed to the shortcomings of the tight-binding parameterization applied here under the conditions of highly excited electronic system. We estimated that the validity of this approximation under the considered experimental conditions (at average absorbed doses of 19-25 eV/atom) maintains on the timescale of up to 20 fs at most. In order to extend the applicability of our code to longer times, the TB parametrization should be replaced with another ab-initio module which can self-consistently account for changes of the electronic band structure and the impact ionization rate within such a highly excited system, out of equilibrium. Such approaches are not yet available, however, promising developments in this direction from TDDFT calculations 39 and Hartree-Fock schemes 40 are underway.

Methods
Simulation tool. Time evolution of X-ray irradiated diamond was modeled with the hybrid code XTANT 24 .
The model employs tight-binding (TB) molecular dynamics (MD) to follow trajectories of atoms within a supercell. TB module also calculates transient electronic structure within the valence band and bottom of the conduction band via a direct diagonalization of the transient Hamiltonian. The so-called hopping integrals, entering TB Hamiltonian, are functions of interatomic distances and positions, adjusted to reproduce atomic structures of various carbonaceous materials 33 . Thus, the method is capable of reproducing different phase transitions.
This enables us to calculate the potential energy surface, Φ({r ij (t)}, t), depending on all the atomic positions pairwise: Here, the first term is the attractive part formed by the electrons represented by the electron distribution function f e (E i ), populating the transient energy levels E i ; E rep is the core-core repulsive potential of ions. Derivative of the potential energy then enters equations of motion. They are solved for all atoms in the simulation box with periodic boundary conditions by applying Verlet algorithm in its velocity form.
To trace the electron distribution function, f e (E i ), which enters Eq. (1), XTANT splits electrons into two fractions, according to their energy: (i) High-energy electrons (with energies above a threshold of ~10 eV 24 ) are modeled as individual particles with the event-by-event Monte Carlo (MC) simulation. Their time propagation is sampled as a sequence of collisions, defined by the cross sections of interactions, which are derived from the complex dielectric function 17,24 . High-energy electrons perform secondary ionizations, exciting new electrons from the valence to the conduction band of diamond. Photoabsorption and decays of core K-shell holes are also traced within the MC scheme. When an electron loses its energy below the low-energy threshold, it joins the electronic low-energy fraction (see below), thereby heating up the low-energy electronic distribution. All the details of the MC scheme including cross sections of scattering were described in ref. 24 . (ii) Low-energy electrons (with energies below the threshold) are assumed to obey Fermi-Dirac distribution at all times. Their evolution is traced with a simplified Boltzmann equation 31 . Boltzmann collision integral for the energy exchange with ions (electron-phonon coupling) reads (as in ref. 25 ): ( ) (2 ( )) ( ) (2 ( )) ( ), for , ( ) (2 ( )) ( ) ( ) (2 ( )), for , where g at (E) is the integrated Maxwellian distribution for atoms with a transient ion temperature, and w i,f is the electron-ion scattering rate during the actual time-step, calculated as follows 25 : where 〈i(t)| and | j(t + δt)〉 are the i-th and j-th eigenfunctions of the Hamiltonian at the time instants t and t + δt, respectively; δt is the simulation time-step 25 .

Evolution of electronic subsystem in X-ray irradiated diamond. Simulations of diamond damage
show that this process proceeds in multiple steps, as discussed in the Introduction. Firstly, valence electrons are excited into the conduction band during and for some time after the pulse, see Fig. 9. The number of the excited electrons is high. It overcomes the damage threshold value (~1.5%) 24 already early during the exposure to the pump pulse 27 . This leads to the band gap collapse at the time ~15 fs since the maximum of the pump pulse, as shown in Fig. 10. It is known from our previous studies 27 that the band gap collapse indicates an onset of the irreversible phase transition, a graphitization-like process. In case of all three absorption doses considered in the experiment, the electronic processes have a similar temporal characteristics. As expected, the fraction of electrons excited increases with the increasing absorption dose.