Electronic nonequilibrium effect in ultrafast-laser-irradiated solids

This paper describes the effects of electronic nonequilibrium in a simulation of ultrafast laser irradiation of materials. The simulation scheme based on tight-binding molecular dynamics, in which the electronic populations are traced with a combined Monte Carlo and Boltzmann equation, enables the modeling of nonequilibrium, nonthermal, and nonadiabatic (electron-phonon coupling) effects simultaneously. The electron-electron thermalization is described within the relaxation-time approximation, which automatically restores various known limits such as instantaneous thermalization (the thermalization time τe−e→0 ) and Born-Oppenheimer (BO) approximation ( τe−e→∞ ). The results of the simulation suggest that the non-equilibrium state of the electronic system slows down electron-phonon coupling with respect to the electronic equilibrium case in all studied materials: metals, semiconductors, and insulators. In semiconductors and insulators, it also alters the damage threshold of ultrafast nonthermal phase transitions induced by modification of the interatomic potential due to electronic excitation. It is demonstrated that the models that exclude electron-electron thermalization (using the assumption of τe−e→∞, such as BO or Ehrenfest approximations) may produce qualitatively different results, and a reliable model should include all three effects: electronic nonequilibrium, nonadiabatic electron-ion coupling, and nonthermal evolution of interatomic potential.


I. Introduction
Ultrafast laser irradiation of materials plays an important role in both, fundamental and applied sciences [1][2][3][4][5].Understanding of basic phenomena in the physics of solids, nonequilibrium kinetics, and highly excited states of matter, benefits from experiments accessing the natural time window of the involved processes (femtosecond timescales) such as electron kinetics, electron-ion (electron-phonon) coupling, and atomic response [5][6][7].Laser irradiation experiments are vital for materials processing, nano and micro technology, and medical applications such as laser surgery [4,5,8,9].In turn, the interpretation of experimental results requires detailed theoretical and model descriptions of the processes involved.
When an intense ultrafast laser pulse arrives at a material, a sequence of processes takes place, ultimately leading to material modification and the formation of a new state [1,2,10].It starts with photoabsorption, which, for photon energies above the bandgap, is primarily by electrons of the matter [11].The photoelectrons are excited to higher energy states, which drives the electronic ensemble out of equilibrium [11,12].The electrons then scatter among themselves in the conduction band of material, with electrons from the valence band and deeper shell (for those shells whose ionization potential is below the electron energy), or elastically with the atoms of the target (also known as the electron-phonon scattering).Electron-electron scattering thermalizes the electronic ensemble towards equilibrium Fermi-Dirac distribution, whereas electron-atom scattering exchanges energy between the two systems, heating the lattice.Typical electronic thermalization in metals takes place at femtosecond timescales [1,13], while equilibration between the electronic and atomic temperatures may take up to tens of picoseconds [1,14].
At the same time, since the electrons form the interatomic potential of the material, the excitation of electrons modifies the potential energy surface [1,10,15].Atoms, in their former equilibrium positions, now experience new forces, which may drive them into a different material phasethe processes known as nonthermal phase transitions, the most famous example of which is nonthermal melting in covalently bonded solids [16,17].Furthermore, in the time window when the electronic temperature (or, more generally, parameters of the nonequilibrium electronic distribution) is different from the atomic one, this so-called two-temperature state departs from the equilibrium material phase diagram and may create states inaccessible under equilibrium conditions [18,19].
The variety of processes, taking place after ultrafast irradiation, may create synergy leading to nonlinear behavior [20].Previously, the nonthermal phase transitions were studied [21][22][23], in a few cases with nonadiabatic coupling included [10,24].Effects of the electronic nonequilibrium were also analyzed separately [1,13,25,26].The mutual effect of the electronic nonequilibrium and the electron-phonon coupling has also been demonstrated [27,28].However, a combined effect of the interplay of the electronic nonequilibrium, nonadiabatic coupling, and nonthermal effects, was not studied before.
In this work, the interplay of these effects is investigated using the developed combined model that simultaneously traces the kinetics of electrons, the evolution of the electronic structure, and the response of the atomic system.It is demonstrated that electron nonequilibrium slows down electron-phonon coupling, leading to prolonged heating of the atomic system.In semiconductors and insulators, electronic nonequilibrium enhances nonthermal effects, altering the damage threshold.

II. Model
To model the response of materials to laser irradiation, XTANT-3 hybrid code is modified to account for electronic nonequilibrium among the low-energy electrons populating the valence and the conduction band [10].The code combines a few different models to trace the essential effects in both, the electronic and the atomic systems of the target.(i) The electrons with energy above a certain chosen energy cutoff are modeled with the Monte Carlo (MC) simulation.(ii) The fractional populations of low-energy electrons on the valence and conduction band energy levels are traced with the Boltzmann equation (BE).(iii) The interatomic forces are calculated from the transferable tight-binding (TB) formalism, which also traces the evolution of the electronic energy levels (molecular orbitals).(iv) The motion of atoms is traced with help of the molecular dynamics (MD).Below, we will discuss the relevant details of each model, and their interconnection on-the-fly, enabling to model laser irradiation of matter.
(i) In the current implementation, photoabsorption is modeled in the linear regime only.No nonlinear or simultaneous multiphoton effects are included; although sequential multiphoton absorption is possible in this model.The photoabsorption is modeled with the Monte Carlo simulation, which uses photoabsorption cross sections from the EPICS2017 database for deep atomic shells (possible for X-ray irradiation) [29], whereas for the valence or conduction band photoabsorption may be included via complex dielectric function [30].
Electrons above a chosen cutoff energy (typically 10 eV counted from the bottom of the conduction band) are also treated with the MC module as free electrons [31].Such a separation of electrons into the lowenergy fraction and the high-energy fraction is necessary within the model for several reasons.As TB energy levels are calculated as the linear combination of atomic orbitals, they are typically limited to energies of about ten eV, and higher energy electrons cannot be reliably traced.On the other hand, modeling high-energy electrons as free particles is a standard approximation in transport MC methods and is significantly easier and faster than more advanced methods such as time-dependent ab initio simulations [32].
The event-by-event MC scheme relies on the asymptotic trajectory method.The impact ionization and (quasi-)elastic scattering processes are included.The impact ionization is modeled with the binaryencounter Bethe (BEB) cross-section [33].The elastic scattering is modeled with the screened Rutherford (Mott) cross-section with the modified Molier screening parameter [34].
Upon impact ionization, a new electron may be excited above the cutoff energy.In this case, such an electron is traced in MC in the same way as the primary photoelectron.When an electron loses energy below the cutoff, it disappears from the MC module and joins the low-energy electrons populating the valence and the conduction band traced with the BE, as will be described below.
In an elastic scattering event, no ionization is produced but the energy is transferred to the atomic system as an increase in the kinetic energy of atoms.The energy transfer is traced from all electrons scattering within the given MD timestep and averaged over the MC iterations.
Core shells, created by X-rays or impact ionization by fast electrons, decay via the Auger channel with the characteristic times from the EPICS2017 database [29].In such a process, two new holes are created, which, if in core shells, will undergo their own Auger decays.High-energy Auger electrons, emitted in decays, are traced in the MC module.
(ii) To describe the dynamics of (low-energy) electrons and atoms, one may start with the Ehrenfest approximation [35].It is based on the mean-field approach, allowing a reduction of the problem to quantum electrons and classical atoms.Electrons are then described with the equation for the evolution of the density matrix, whereas expectation values are tracked for atoms, which reduces to the classical (Newton's) equations with the potential defined by the transient state of the electronic system.It is further assumed that the diagonal elements of the density matrixthe electron populations or distribution function,   (  , ),may be traced with a semiclassical Boltzmann equation.This assumption also requires statistical averaging over the electronic ensemble, which makes the electron populations fractional.The off-diagonal elements are then only represented as the scattering integrals, defining the evolution of the electronic distribution: The distribution function describes fractional electronic populations on discrete energy levels   = ⟨|  |⟩, which are the eigenfunctions of the TB Hamiltonian (described below) at the current MD timestep.Here  − is the electron-electron collision integral responsible for electronic thermalization;  − is the electron-atom collision integral defining the nonadiabatic energy exchange between electrons and atoms induced by atomic motion (e.g., electron-phonon coupling);   is the source term describing the change of the distribution function induced by the processes accounted for in the MC module (photoabsorption, Auger-decays involving valence/conduction band, high-energy electrons scattering and influx).The standard Ehrenfest dynamics, being a mean-field approximation, does not by itself account for electronic thermalization.Thus, in Eq.(1), the term  − is added ad hoc, knowing that the Boltzmann equation must lead to thermalization.
For the electron-atom scattering,  − , the model developed in Ref. [14] is employed.It is a semiclassical Boltzmann collision integral, in which the probability of electronic transition between energy levels induced by atomic displacement at each MD timestep is calculated from the overlap of the electronic wavefunctions, obtained within the TB (described below).The classical Maxwellian distribution is used for atoms with the kinetic temperature calculated directly from the MD velocities, whereas for electrons, Pauli blocking factors are included [14].
For the electron-electron scattering, the relaxation time approximation is used in the current implementation [36]: Here  − is the characteristic electron-electron relaxation time;   ( Eqs. (3) define the equivalent electronic temperature (also called the kinetic temperature,   [37]) and the equivalent chemical potential (µ).Within this ansatz, the total number of low-energy electrons and the total energy (in electrons and atoms) are conserved within an MD timestep (changes in the number of electrons may only occur via the   term).
Eqs. (1,2) naturally unify various widely used approaches to quantum-classical dynamics.It recovers the limiting cases of the Born-Oppenheimer (BO) molecular dynamics (in the limit of infinite electronic thermalization time,  − → ∞ , and no nonadiabatic electron-atom coupling,  − = 0 ) in which the electronic populations are fixed; the Ehrenfest dynamics which includes average electron-atom energy exchange but no electron thermalization (  − → ∞ ,  − ≠ 0 ); instantaneous thermalization in the adiabatic microcanonical ensemble ( − = 0,  − = 0, used e.g. in Refs.[31,38]); and nonadiabatic dynamics with instantaneous electron thermalization ( − = 0,  − ≠ 0, used in the previous versions of XTANT [10,14]; the same assumptions are used in the two-temperature based models, TTM-MD [39][40][41]).Let us also note that the BO approximation only assumes the decoupling of the electronic and atomic wave functions, but not necessarily the ground stateexcited states adiabatic states may also be calculated in many cases, where potential energy surfaces are far apart, thus suppressing electronic transitions between them (far from such situations as conical intersections or avoided crossings) [42].The ground-state BO molecular dynamics is a separate additional approximation, which can be reproduced within the relaxationtime formalism as the zero-temperature instantaneous thermalization with no coupling ( − → ∞,  − = 0, and   = 0).
The influx or outflux of electrons from the valence and conduction band (low-energy electrons) via such processes as photoabsorption, scattering of high-energy electrons, and Auger decays, are traced in the MC.In each scattering act involving the valence and conduction electrons, the energy level involved is sampled and recorded into the   .The levels are sampled according to probability following the Pauli blocking term in the Boltzmann collision integral:   ~(  )(2 − (  + ∆)), where ∆ stands for the energy transferred in the collision under consideration (factor 2 is due to spin degeneracy).All the changes in populations in each energy level within the given timestep are then averaged over the MC iterations.
A single discrete level from which an electron is emitted (in photoabsorption, impact ionization, or Auger decay of a core hole) can be chosen.However, an incoming electron generally comes with energy somewhere in between the discrete energy levels (since the energy grid cannot be chosen arbitrarily but is unambiguously defined as the eigenstates of the Hamiltonian).In this case, the influx of the particles and energy is distributed between the two closest levels under the condition of conservation of the number of particles and energy: increasing the total number of low-energy electrons by one out of   iterations and bringing energy ∆ in the scattering event: Where for one incoming electron out of   , the change in the number is ∆ = 1/  , and the corresponding energy change is ∆ =   /  (  is the energy brought by this electron).This way, the total number of electrons and energy, in the low-and high-energy fractions (MC and BE modules), are conserved.In the case when one of the levels (i or j=i+1) is fully occupied, another set of levels can be chosen -Eqs.(4) hold for arbitrary levels i and j.The total change of the population on each level is then summed over all scattering acts within the given timestep:   = ∑ ∆  (  , ).
(iii) The tight-binding method is used to calculate energy levels and interatomic forces (potential energy surface) with the Slater-Koster approximation [43].To make the method transferable to multiple atomic structures, the radial part in the hopping integrals is a function depending on the interatomic distance [44].For elemental metals, we use NRL TB parameterization [45,46], for Si the parameterization from Ref. [47] is used, and for compounds, we employ DFTB matsci-0-3 parameterization [48].Nonorthogonal TB Hamiltonian is constructed on each MD timestep of the simulation and diagonalized by solving the secular equation [44].TB wavefunctions (eigenfunctions of the transient Hamiltonian) can be used to calculate the probabilities of the nonadiabatic electron coupling to atomic motion [14].Unfortunately, the tight binding method does not allow evaluating electron-electron scattering probabilities at the same level of theory, hence the characteristic relaxation time  − is an external parameter.Note, however, that further approximations can, in principle, be made for calculations of the electron-electron scattering probabilities or cross sections (such as linear response theory), which may be the subject of a future work.
(iv) Molecular dynamics module in XTANT-3 uses Martyna-Tuckerman's 4 th order algorithm to propagate classical trajectories of atoms typically with the timestep smaller than 1 fs [49].Over 200 atoms are used in each simulation box with periodic boundary conditions to mimic bulk materials.No energy sinks associated with particles and energy transport out of the simulation box are included, assuming homogeneous material excitation in all presented simulations.
The forces acting on the atoms are calculated as the gradients of the potential energy defined within the TB formalism: where the potential V depends on pairwise distances between all the atoms in the simulation box {  ()};   is the effective ion-ion TB repulsion term (containing all contributions apart from the electronic band energies, if required in the given TB parameterization).Additionally, the nonadiabatic energy transfer from the Boltzmann collision integral and the energy transfer from high-energy electrons elastic scattering are fed to atoms via velocity scaling at each timestep, which ensures the energy conservation in the entire system: all electrons and atoms (microcanonical ensemble).
Using Eq.( 1) allows directly tracing an effect of the nonequilibrium electronic distribution function on the interatomic forces in Eq.( 5) (as well as on the electron-atom coupling,  − [14]), and thus the dynamics and stability of the material.As mentioned above, various approximations for the relaxation time also enable comparisons of different standard methods (e.g., BO, instantaneous thermalization approximation).
Note that the formalism, at its core, relies on the Ehrenfest dynamics, and not on the finite-temperature extension to ab initio simulations [50].The difference is, the underlying physical picture of the finitetemperature models is that the constant electronic temperature is enforced by the interaction with the bath.In our case, however, even in the limit of the instantaneous electron thermalization (  − → 0), no interaction with the bath is assumed, and the electronic thermalization is an intrinsic process (non-meanfield effect added).The electronic temperature and chemical potential are not variables in this formulation but parameters.Let us also note that Eq.( 2) could be replaced with the electron-electron Boltzmann collision integral, which does not employ any electronic temperature or chemical potential, even as parameters of the equivalent equilibrium distribution [13].

Metals
We start with the study of the evolution of the electron distribution function in irradiated metals in different approximations for the electron-electron characteristic relaxation time.For the study, we chose the photon energy of 2 eV (wavelength of 620 nm) for metal irradiation; different photon energies will be compared below for the case of semiconductors.In all simulations, the temporal profile of the laser pulse is assumed to be Gaussian with 10 fs FWHM centered at 0 fs. Figure 1 shows the electronic distribution function in aluminum under irradiation with the optical pulse simulated with various approximations:  − = 0 fs (instantaneous electron equilibration),  − = 15 fs,  − = 150 fs, and  − → ∞ (no electron-electron relaxation; Ehrenfest-like approximation).In the instantaneous thermalization limit ( − = 0 fs, Figure 1a), the distribution function coincides with the equilibrium Fermi distribution at all times.In simulations with finite thermalization times, the nonequilibrium electronic distribution function transiently forms a staircase-like shape, in agreement with the previous studies [13].When fast electronic thermalization is assumed ( − = 15 fs, Figure 1b), the distribution function quickly turns into equilibrium Fermi function.With longer characteristic electronelectron thermalization times ( − = 150 fs, Figure 1c), the deviations from the Fermi functions are stronger and the equilibration takes proportionally longer.

Figure 2. Evolution of the electronic entropy in aluminum irradiated with a pulse of 2 eV/atom, 2 eV photon energy, 10 fs FWHM. Simulations performed with different characteristic electronic relaxation times: (a) τ=0 fs (instantaneous relaxation); (b) τ=15 fs; (c) τ=150 fs; (d) τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation). Solid lines are the entropies corresponding to the transient distribution functions, dashed lines are the entropy associated with the equivalent Fermi distributions (maximal entropy).
The most surprising is that the Ehrenfest-like dynamics (with no electron-electron scattering,  − → ∞, Figure 1d) also leads to the eventual thermalization of electrons via electron-phonon scattering, when the kinetic temperatures of electrons (the equivalent temperature of the equilibrium distribution function according to Eq.( 3)) and ions equilibrate.This effect is quantified in Figure 2, where the electronic entropy (  = −  2 ∑[(  /2) ln(  /2) + (1 −   /2) ln(1 −   /2)], with   being the Boltzmann constant, and factors of 2 due to the normalization of the distribution function according to the spin degeneracy) is shown.For short electron-electron thermalization times, the entropy quickly reaches the maximal value (Figure 2a  and b), which corresponds to the equilibrium distribution function.In all cases, the maximal entropy is decreasing with time due to the cooling down of the electronic ensemble via electron-phonon coupling.In the Ehrenfest-like approximation, the entropy coincides with the maximal one at the time of ~1.5 ps, indicating electronic thermalization.
It is important to keep in mind that this is only an effect of one channel of electron interaction, and, despite the apparent thermalization, the actual channel responsible for the equilibration of electrons is missing in this approximation.Including electron-electron interaction leads to much faster electronic thermalization (compare, e.g., Figure 2b vs. Figure 2d).
The nonequilibrium electron distribution affects the electron-phonon coupling, and thereby the atomic temperature evolution in irradiated aluminum, see Figure 3.The thermalized electrons ( − = 0 fs) couple to ions/phonons most efficiently, resulting in the fastest increase of the atomic temperature, as was also noted in previous works [13].Strongly nonequilibrium electronic ensemble ( − → ∞) produces the smallest electron-phonon coupling and, correspondingly, the slowest atomic heating.The coupling parameter depends on both, the electronic and atomic temperature (among other parameters such as material structure and density), and thus evolves in time [14].The kinetic (equivalent) temperatures in the case of nonequilibrium distributions are used to define the coupling parameter in Figure 3b; also note that it is only defined for   ≠   , thus the timescale in the figure is limited to before equilibration of the electronic and atomic temperatures.The difference in atomic temperatures, however, is not large among all the simulations, and by the time of ~1500 fs, the atomic temperature reaches the plateau (equilibrated with the electronic kinetic temperature) in all cases in aluminum after deposition of 2 eV/atom absorbed dose.The difference at other deposited doses is comparably small, see Appendix.Thus, in the case of optical irradiation, the effect of electronic nonequilibrium on the electron-phonon coupling and atomic heating in aluminum appears to be relatively small.The same can be concluded for other studied metals: gold and tungsten (see Appendix).
Since in bulk metals electronic excitation leads to a phonon hardening effect [51], no nonthermal damage is observed in the current simulations.The situation may be different in nano-sized metallic samples where nonthermal phase transitions can take place in response to an increase of the electronic pressure and ensuing material expansion [51], but finite-sized samples are beyond the scope of the present work and maybe a subject of future dedicated research.
Having studied the effect of the nonequilibrium in the wide range of possible relaxation times, let us now consider the experiment measuring the photon spectra emitted from L-shell decays in aluminum following irradiation with a 35-femtosecond free-electron laser (FEL) pulse of 92 eV photon energy [52].In the soft X-ray regime, the photons are predominantly absorbed by the L-shell (2p atomic shell) of aluminum, with only a minor contribution of the conduction band photoabsorption.The L-shell holes formed during photoabsorption subsequently undergo decay, which is primarily Auger decay, with radiative decay contributing to a lesser extent.The radiative decays produce photons with energies defined by the ionization potential of the L-shell (~72 eV), which are detected experimentally.The detected spectrum is naturally averaged over the lifetime of the L-shell holes, which is ~40 fs in solid aluminum [52].

(b) Electronic entropy, simulated with various characteristic relaxation times of low-energy electrons. (c) The density of holes in various atomic shells. (d) Energies in various systems: total energy in the simulation box, energy of atoms and low-energy electrons, and potential energy of atoms (without the kinetic energy).
We conducted simulations of aluminum irradiated with an FEL pulse with the parameters corresponding to those of the experiment and 500 atoms in the simulation box. Figure 4a illustrates the calculated electronic distribution function.The characteristic shape of the bump-on-hot-tail distribution, typical for FEL-irradiated materials, is emerging during and after the FEL pulse [25,53].The high-energy tail originates from the conduction-band photoabsorption during the pulse and is subsequently sustained by the Auger-decays of L-shell holes (see a detailed discussion in [25]).This tail lives as long as the L-shell holes are still present in the material, which may last for a few hundred femtoseconds (exponentially decaying with the characteristic time of ~ 40 fs, see Figure 4c) [54].However, the low-energy fraction of electrons around the Fermi level is much closer to partial equilibrium, as was also discussed in [25].These electrons participate in the radiative decays of L-shell holes, forming the experimentally observed spectrum.
The deviation of the electronic ensemble from equilibrium is not large, as seen by the electronic entropy (Figure 4b), even for the relaxation time of the low-energy electrons τ=100 fs, however, the minority of high-energy electrons transiently stores some energy; a larger amount of energy is transiently stored in Lshell holes.Figure 4d shows that the energy of L-shell holes and high-energy electrons (the difference between the total energy plotted and the energy of atoms and slow electrons) transiently may be dominant.This figure also demonstrates that the potential energy of atoms starts to decrease due to heating, and the corresponding energy transfers to the kinetic energy of atoms, heating the system.It also can be seen that the total energy in the simulation box is conserved (apparat from the deposition by the FEL pulse, centered around zero), validating the developed simulation scheme.
The calculated distribution function averaged over time with the transient density of L-shell holes (Figure 4), is then compared with the distribution function extracted from the experimental spectra (by division of the spectrum by the density of states).The results of this comparison are depicted in Figure 5, for the few relaxation times employed in XTANT-3 simulation.It is apparent from the figure that the closest match between the two distribution functions occurs for the characteristic time of ≤10 fs.

Figure 5. Low-energy parts of the electron distribution in aluminum irradiated with 35 fs FWHM 92-eV photon energy laser pulse, simulated
with XTANT-3 with various characteristic relaxation times, for the deposited doses of 1.8 eV/atom (a) and 5.8 eV/atom (b), compared with the experimental distributions extracted from the corresponding spectra from Ref. [52].
Note that this is the characteristic relaxation time of only the low-energy fraction of electrons, whereas the total thermalization of the electronic system requires considerably more time (cf. Figure 4).The lifetime of the long nonequilibrium tail is defined by the slowest process, which in the case of excitation of the Lshell of aluminum is the Auger-decay time.This process continuously supplies the distribution function with new out-of-equilibrium electrons at the energies of ~72 eV for several hundred femtoseconds, as was previously investigated in [25] and experimentally confirmed in [54].
This comparison allows us to conclude that the relaxation time approximation, combined with MC for high-energy electrons, is a reliable model, capable of describing experimental data with reasonable accuracy.

Semiconductors
We use silicon as the main target to study electronic nonequilibrium effects in semiconductors.Here, again, the laser pulse of 10 fs FWHM with the photon energy of 2 eV is used.The deposited dose is set to 0.9 eV/atom, which is close to the nonthermal melting threshold in silicon [55].Similar to the case of aluminum above, the distribution function of electrons in silicon forms a nonequilibrium shape during the laser pulse and secondary electron cascades, see Figure 6.The bandgap, seen by the positions of the points in the plotted distribution function, closes after a few hundred femtoseconds (for details, see e.g.Ref. [55]).The thermalization takes longer for larger electron-electron relaxation times, as seen by the behavior of the electron entropy in Figure 7. Due to the presence of the bandgap at the beginning, before silicon nonthermally melts and turns metallic with the collapse of the gap, some energy is stored in electrons as the potential energy, which slows down electron thermalization.By the time of 1 ps, the thermalization of the electronic ensemble is almost completed even in the simulation with no electron-electron relaxation included ( − → ∞)the same effect as described above for aluminum irradiation.Nonequilibrium in the electronic system affects the electron-phonon coupling in silicon stronger than in metals, see Figure 8: the electronic kinetic temperatures for large electron-electron relaxation times show a noticeable difference from the short-relaxation-time simulations.However, the atomic temperatures differ only mildly.
Nonthermal damage mechanisms are also affected by electronic non-equilibrium.Figure 9 shows that the nonthermal melting threshold first lowers with the increase of the electron-electron relaxation time, up to the times of  − ~150 fs, where the damage threshold is ~0.85 eV/atom vs. ~0.9eV/atom within the instantaneous thermalization (  − = 0 fs) reported earlier [55].For the thermalization times over  − ~250 fs, the nonthermal damage threshold increases again, reaching ~0.9-1.0 eV/atom in the limit of  − → ∞.
In germanium, the situation is qualitatively similar, albeit the change in the damage threshold dose becomes more noticeable at larger thermalization times, see Figure 9.This could be expected, considering that Ge atoms are heavier than Si, and it takes longer for them to react to electronic excitation and melt the lattice nonthermally.So, longer-lasting electronic nonequilibrium is required to influence this process.The nonlinear dependence of the damage threshold on the electron-electron thermalization time is mainly defined by the transient density of the electrons excited to the conduction band, see an example in Figure 10.Electrons, promoted to the conduction band by photoabsorption and relaxing via electron-electron scattering, excite secondary electrons into the conduction band, thereby increasing their density.When this channel is absent, fewer conduction band electrons are excited.With the increase of the relaxation time up to  − ~250 fs, the density of the conduction band electrons holds over-the-threshold values (above ~4 % of the valence band electrons [55]) for a longer time, allowing the atoms to react and form nonthermal damage.However, the peak density simultaneously reduces, and for some values of the relaxation times ( − ~250 fs) becomes insufficient to induce nonthermal melting (Figure 10).Thus, to reach the threshold electron density in simulations with larger relaxation times, a higher dose is required.It is interesting to note that the earlier results obtained with the instantaneous thermalization approximation ( − = 0) predicted that the damage threshold in terms of the absorbed dose is independent of the photon energy (under the assumption of homogeneous photon absorption and no energy sinks from the sample) [55].But nonequilibrium simulation produces a different damage threshold that depends on the particular shape of the transient electron distribution function.Therefore, the threshold absorbed dose depends on the photon energy.We demonstrate it in the example of the largest effect of nonequilibrium, the limit of  − → ∞, see Figure 11.[56] and references therein for XUV and x-rays, and from [57] for the optical pulse.

All simulations are performed with the characteristic electronic relaxation time τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation). The dashed horizontal line is the threshold dose in simulations with instantaneous electron thermalization (τ=0 fs). The crosses are experimental data from
For photon energies in the optical and ultraviolet range, there is a nonlinear dependence of the ultrafast nonthermal damage threshold.So, each particular photon energy requires a dedicated study.In contrast, for photon energies above the width of the valence band (~ 15 eV), the secondary cascades of high-energy electrons (photo-, Auger-, and impact-ionized) only differ in duration (see Ref. [58]), but the distribution of the secondary excited electrons does not differ significantly.Thus, with the increase of the photon energy above some ~15 eV, the damage threshold dose stays nearly constant up to the highest photon energy of 10 keV studied here.This result suggests that an indirect validation of the electron-electron relaxation time is potentially possible by experiments measuring the damage threshold.If all the other uncertainties of the simulation and experimental procedure were eliminated, the damage threshold dose would enable us to extract average electron-electron thermalization time by comparison with simulations.Unfortunately, the existing experimental data do not allow for discrimination between different characteristic thermalization times in the present calculations (see Figure 11).In Figure 11, the comparisons are shown for XUV and X-ray experiments with only one point for the optical pulse, which was converted from the excited electron densities from Ref. [57], since in most of the optical-laser experiments, the threshold is reported only in terms of the incident fluence.The nonlinear photoabsorption and other ensuing effects in the optical case make it difficult to convert the thresholds from incoming fluence to the absorbed dose [59].

Insulators
To study the effects of electron-electron equilibration in an insulator, we use Al2O3.Since, in the current implementation, XTANT-3 only includes linear photoabsorption, we cannot model photon energies smaller than the bandgap of an insulator, which is Egap~9 eV in Al2O3.Thus, we chose an XUV pulse with a photon energy of 30 eV (wavelength of 41 nm).In this case, the high-energy electrons created due to photoabsorption, relax to the energies below the cutoff of 10 eV within ~1 fs, joining the low-energy fraction of electrons (utilizing Eqs. ( 4)) already during the laser pulse [58].
The evolution of the electronic distribution in irradiated Al2O3 is shown in Figure 12.Both, the transient nonequilibrium distribution and the equivalent equilibrium distribution function plotted have initially a large discontinuity at the bandgap.As the material starts to respond to the deposited energy, some energy levels are shifting into the gap, see Figure 13, and the electronic populations are changing [60].In the case of fast thermalization ( − = 0 fs and  − = 15 fs), the kinks and spikes in the distribution function relax fast towards the Femi distribution, as expected.Slower thermalization time of  − = 150 fs demonstrates persisting nonequilibrium at times of 1 ps.In the absence of electron-electron scattering ( − → ∞ ), the electron interaction with phonons leads to only partial thermalization in the electronic ensemble, in contrast to the case of metals and semiconductors discussed above.In an insulator, phonons cannot provide or accept sufficient energy for electrons to cross the bandgap, thus only separate intraband electron thermalizations are observed.Each bandconduction, valence, and subvalenceestablishes its own separate distribution, close to the Fermi function (see Figure 12d).There is only a slow small leak of electrons from the conduction band to the valence band due to the lowering of the bands and occasional crossing of the valence and the conduction levels in the highly excited material, see Figure 13.Crossing or avoided crossing of energy levels allows for nonadiabatic electron transitions between them [61,62].These processes are reflected in the electronic entropy being noticeably below the maximal value corresponding to the complete thermalization of all electrons, see Figure 14.The full electron thermalization within the Ehrenfest-like approximation does not take place in insulators at considered doses.In the case without electron-electron relaxation, a large amount of energy is stored as the potential energy of electrons: the energy that could be released if electrons had a mechanism to relax across the band gaps.The effect of it is twofold.On the one hand, the potential energy stored cannot be transferred to atoms/phonons via the nonadiabatic electron-ion coupling.This leads to slower heating of atoms and a lower final atomic temperature, see Figure 15.In the Ehrenfest-like approximation, the kinetic electronic and atomic temperatures reached saturation by the time of ~1 ps, and do not equilibrate.
On the other hand, without direct electron-electron relaxation, valence holes and conduction band electrons last longer than in the simulation where electron-electron relaxation is included.That means, the electronic population on the bonding levels in the valence band is reduced, whereas the population of electrons on the antibonding states in the conduction band is increased, as can be seen by the persisting nonequilibrium shape of the distribution, recall Figure 12.The interatomic potential is directly affected by the electronic distribution on the levels (see Eq.( 5)) [18,60].This nonequilibrium electron distribution leads to ultrafast nonthermal phase transition in Al2O3 at doses lower than the equilibrium one: ~2.4 eV/atom in the case of the instantaneous thermalization ( − = 0 fs) vs. ~2.2eV/atom without electron-electron thermalization ( − → ∞).The state after the phase transition seems to be the same independently of the electron-electron thermalization timetransiently superionic state with the liquid oxygen subsystem and the solid aluminum lattice, for details see Refs.[18,60].
Thus, it can be concluded that the electronic nonequilibrium in insulators plays an important role in influencing both, nonadiabatic (electron-phonon) coupling and nonthermal phase transitions.The same effects are expected to take place in other covalent insulators where the nonthermal bandgap collapse is induced by ultrafast irradiation.The ionic crystals may behave differently since they typically do not experience bandgap collapse under irradiation [63].This topic is beyond the scope of the present work and may be studied in the future.

IV. Conclusions
Materials response to ultrafast laser irradiation was studied theoretically using a comprehensive hybrid simulation tool XTANT-3.The evolution of the nonequilibrium electronic distribution function on the transient energy levels was included, which enabled accounting for the interplay of nonequilibrium, nonadiabatic and nonthermal effects.The solution of the Boltzmann collision integrals developed is particle and energy conserving by construction.The model was applied to various classes of solids: metals, semiconductors, and insulators.
The results suggest that models based on Ehrenfest-like dynamics may be applied to describe the atomic response to irradiation of metals and semiconductors at doses above nonthermal melting, but they fail to describe insulators since electronic thermalization across the bands is excluded.Such models exhibit equilibration within the electronic ensemble mediated by the electron-phonon coupling, but on much longer timescales than the actual thermalization via electron-electron interaction.The non-equilibrium distribution of electrons influences both, nonadiabatic electron-phonon coupling and nonthermal modification of the interatomic potential.
The nonequilibrium electron distribution leads to smaller electron-phonon coupling and, correspondingly, slower atomic heating in all studied cases (all materials, irradiation doses, and photon energies).In metals, the difference in the atomic response simulated with different models is relatively small: only slightly slower atomic heating is observed for larger electron-electron relaxation times.Comparison of the simulated electron distribution function with the experimental emission spectra enabled estimating the thermalization time in aluminum as ≤10 fs.However, longer-lasting non-equilibrium may be expected in some cases [64], or, perhaps, under long-pulse irradiation continuously driving the electronic system out of equilibrium.
In semiconductors, the effect of electronic nonequilibrium on the electron-phonon coupling is similar, but it also plays a role in nonthermal melting.Depending on the electron-electron relaxation time, the nonthermal damage threshold may be lowered or increased, since electron-electron scattering affects the excitation of electrons across the bandgap, and the conduction-band electrons and valence holes are the main driving force behind the nonthermal melting.
In insulators, electronic nonequilibrium leads to a lower nonthermal damage threshold, since unrelaxed electrons excited to the conduction band affect the interatomic potential stronger than that in the case of assumed instantaneous electronic thermalization.Without electron-electron thermalization, there is no mechanism for electrons to hop between the conduction and the valence band, which strongly affects both: slows down electron-phonon coupling precluding electron-ion thermalization, and changes the threshold of nonthermal phase transitions.Thus, an appropriate model of the insulator's response to ultrafast irradiation should include all three effects: electronic nonequilibrium, nonadiabatic electron-ion coupling, and nonthermal evolution of interatomic potential.resources were supplied by the project "e-Infrastruktura CZ" (e-INFRA LM2018140) provided within the program Projects of Large Research, Development and Innovations Infrastructures.

Appendix
Figure 16 shows the evolution of the electronic and atomic kinetic temperatures in aluminum irradiated with 10 fs FWHM laser pulse, 2 eV photon energy, and various absorbed doses.Two simulations are compared: with the electron-electron relaxation time of  − = 15 fs vs.  − → ∞.For all deposited doses, the difference between the two simulation schemes is comparably small.Due to the dependence of the electron-ion (electron-phonon) coupling on both, electronic and atomic temperatures, nonlinear behavior of the temperatures is observed: for higher deposited doses (peak electronic kinetic temperatures), the temperatures equilibrate faster than for lower doses, see discussions in [24,65].Figure 17 presents the kinetic electronic and atomic temperatures in gold and tungsten irradiated with 10 fs FWHM laser pulse, 2 eV photon energy, and 2 eV/atom absorbed dose, using various electron-electron relaxation times.The difference among all the cases is relatively small, indicating that the electronic nonequilibrium does not influence strongly the electron-ion coupling in these metals, supporting the conclusions in the main text.It appears to be the smallest in the case of the slowest electron-phonon coupling (typically, heavier elements [14]).

Figure 3 .
Figure 3. (a) Electronic and atomic temperatures, and (b) electron-phonon coupling parameter in aluminum irradiated with a pulse of 2 eV/atom, 2 eV photon energy, 10 fs FWHM.Simulations performed with different characteristic electronic relaxation times: τ=15 fs; τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation).Solid lines are the electronic kinetic (equivalent) temperatures, and dashed lines are the atomic kinetic temperatures.

Figure 4 .
Figure 4. Parameters of aluminum irradiated with 92 eV photon energy, 35 fs FWHM pulse, 6 eV/atom absorbed dose.(a) Electron distribution function.(b)Electronic entropy, simulated with various characteristic relaxation times of low-energy electrons.(c) The density of holes in various atomic shells.(d) Energies in various systems: total energy in the simulation box, energy of atoms and low-energy electrons, and potential energy of atoms (without the kinetic energy).

Figure 8 .
Figure 8. Electronic and atomic kinetic temperatures in silicon irradiated with a pulse of 0.9 eV/atom, 2 eV photon energy, 10 fs FWHM.Simulations performed with different characteristic electronic relaxation times: 0 fs, 15 fs; 150 fs, 1 ps, and τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation).Solid lines are the electronic kinetic (equivalent) temperatures, and dashed lines are the atomic kinetic temperatures.

Figure 9 .
Figure 9. Damage threshold dose in silicon and germanium irradiated with a pulse of 10 fs FWHM as a function of the characteristic electronic relaxation times.

Figure 11 .
Figure 11.Nonthermal melting threshold deposited dose in silicon irradiated with a pulse of 10 fs FWHM as a function of the photon energy.All simulations are performed with the characteristic electronic relaxation time τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation).The dashed horizontal line is the threshold dose in simulations with instantaneous electron thermalization (τ=0 fs).The crosses are experimental data from[56] and references therein for XUV and x-rays, and from[57] for the optical pulse.

Figure 15 .
Figure 15.Electronic and atomic kinetic temperatures in Al2O3 irradiated with a pulse of 2 eV/atom, 30 eV photon energy, 10 fs FWHM.Simulations performed with different characteristic electronic relaxation times: 0 fs, 15 fs; 150 fs, and τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation).Solid lines are the electronic kinetic (equivalent) temperatures, and dashed lines are the atomic kinetic temperatures.

Figure 16 .
Figure 16.Electronic and atomic temperatures in aluminum irradiated with a pulse of various deposited doses, 2 eV photon energy, 10 fs FWHM.Simulations performed with different characteristic electronic relaxation times: (a) τ=15 fs; (b) τ=∞ (no electron-electron relaxation; Ehrenfest-like approximation).Solid lines are the electronic kinetic (equivalent) temperatures, and dashed lines are the atomic kinetic temperatures.

Figure 17 .
Figure 17.Electronic and atomic temperatures in (a) gold and (b) tungsten irradiated with a pulse of various deposited doses, 2 eV photon energy, 10 fs FWHM, and different characteristic electronic relaxation times.Solid lines are the electronic kinetic (equivalent) temperatures, and dashed lines are the atomic kinetic temperatures.