Tracing X-ray-induced formation of warm dense gold with Boltzmann kinetic equations

In this paper, we report on the Boltzmann kinetic equation approach adapted for simulations of warm dense matter created by irradiation of bulk gold with intense ultrashort X-ray pulses. X-rays can excite inner-shell electrons, which triggers creation of deep-lying core holes. Their relaxation, especially in heavier elements such as gold (atomic number Z=79\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z= 79$$\end{document}) takes complicated pathways, involving collisional processes, and leading through a large number of active configurations. This number can be so high that solving a set of evolution equations for each configuration becomes computationally inefficient, and another modeling approach should be used instead. Here, we use the earlier introduced ’predominant excitation and relaxation path’ approach. It still uses true atomic configurations but limits their number by restricting material relaxation to a selected set of predominant pathways for material excitation and relaxation. With that, we obtain time-resolved predictions for excitation and relaxation in X-ray irradiated bulk of gold, including the respective change of gold optical properties. We compare the predictions with the available data from high-energy-density experiments. Their good agreement indicates ability of the Boltzmann kinetic equation approach to describe warm dense matter created from high-Z materials after their irradiation with X rays, which can be validated in future experiments.


Introduction
Highly excited matter created after intense X-ray radiation is an object of intense experimental studies with high power laser sources, in particular, with free electron lasers (FELs), see e.g., [7,9,21]. The experiments trace non-equilibrium mechanisms of plasma formation, also through a transient state of warm dense matter [12,22,24]. Ultrashort duration of FEL pulses also makes it possible to probe the unexplored regime of electronic thermalization which takes up to several tens of femtoseconds since the beginning of X-ray exposure. Further dedicated experiments are underway. However, the analysis of such experimental results requires development of theoretical tools able to describe plasma formation under strong non-equilibrium conditions. For example, modeling of the ionization dynamics can con-veniently be done in many applications with a continuum approach [5,17,25]. With such an approach, evolution equations in phase-space are formulated for distributions of electrons, atoms and ions. They are solved on a phase-space grid. As computational costs depend only on the grid size, one then avoids a direct scaling of the computational costs with particle number, O(N 2 ), typical for particle approaches. This feature of continuum models allows one to also treat large samples efficiently.
For the simulation of the early non-equilibrium stages of the sample evolution, full kinetic equations should be applied. They follow sample evolution during this stage, delivering the information on transient electron and ion distribution. Such equations should treat every active atomic configuration appearing during sample excitation and relaxation.
In this paper, we report on the application of the Boltzmann kinetic equation approach for simulations of bulk gold irradiated with X-ray pulses, performed with Boltzmann equation solver combined with the 'predominant excitation and relaxation path' approach (PERP) introduced in [32]. The Boltzmann equation solver was originally designed to describe the evolution of 'plasma-like' samples composed of atoms, ions and free electrons [25]. Therefore, the bulk Au is assumed to be initially composed of (unbound) atoms. After X-ray excitation, ions and free electrons appear. Whereas for samples irradiated with VUV radiation the number of active atomic configurations is small, as the incoming photons can excite only outer-shell electrons, for samples excited with X-rays this number rapidly increases, as the Xrays can excite inner-shell electrons. This triggers the creation of deep-lying core holes, which complex relaxation, especially in heavier elements, involves a large number of active configurations. To illustrate, in carbon (Z = 6) the total number of possible atomic configurations, with 6 electrons distributed on levels (1s, 2s, 2p), is 27. It is then still possible to solve the set of evolution equations for all the configurations. The same can be done for other light elements. However, already for noble gas atom, neon (Z = 10), the total number of atomic configurations is 63. For argon (Z = 18), it amounts to 1323. This steep increase of the configuration number gives a clear limitation for kinetic simulations of high-Z materials, if including all active configurations.
To overcome this difficulty, the existing approaches such as, e.g., [4] use a superconfiguration approach, i.e., they do not follow the evolution of individual configurations but, instead, use a simplified set of 'averaged' configurations [11,14]. In [32], we introduced an alternative approach which still uses true atomic configurations but limits their number by restricting the sample relaxation to the predominant excitation and relaxation pathways. Applying this approach in what follows, we will obtain time-resolved predictions for excitation and relaxation of X-ray-irradiated bulk of gold, including the change of gold optical properties in response to an X-ray pulse.
The paper is organized as follows. First, we will recall kinetic equation formalism for X-ray-irradiated samples. Especially, simulations of irradiated bulk material will be discussed. With this scheme, we will follow Xray irradiation of bulk gold, leading to the formation of warm dense matter. Bulk gold has been selected due to many available theoretical and experimental data on its behavior in the warm dense matter regime (see, e.g., [16]). The results obtained for electronic and optical properties of the irradiated sample will then be discussed, and compared to the available data from highenergy-density experiments. Finally, we will present our conclusions and outlook.

Description of Boltzmann equation solver
For the description of X-ray excited bulk Au, we have applied our solver of Boltzmann kinetic equations developed in [6,[25][26][27][28][29][30][31][32]. The classical Boltzmann equations originate from the reduced N-particle Liouville equations and include only single-particle phasespace densities of ions and free electrons in the sample. The equations then model all systems as samples built of atoms/ions and of free electrons represented by their classical phase-space densities. In the atomistic approximation, the resulting kinetic equations for the electron distribution in phase-space, ρ (e) , as well as the kinetic equations for various Au ion configurations, ρ (i,j) are formulated as follows: for electrons, and: for ions, with the index, i = 0, . . . , N J , describing ion charge (where N J is the highest charge state present in the system), and the index, j = 0, . . . , N C (i), being the active configuration number (where N C (i) is the maximal number of ion configurations considered for a given ith ion charge). Here, m is the electron mass, and M is the ion mass.
In the most general case, the electromagnetic force acting on a single electron, , with e being the magnitude of electron charge, has two components, the electric one, E (r, t) and the magnetic one, B(r, t). They describe the interaction of charges within the sample with the external laser field as well as mutual electromagnetic interactions between charges. In our specific case of X-ray irradiation, the magnetic component of the electromagnetic interaction, B(r, t), as well as the driving effect of the laser field, E (r, t) on charged particles can be neglected (for details, see [25]). The force then reduces to the component describing mutual electrostatic interactions between electrons and ions in the sample. This component is a non-local function of electron and ion densities and was described in detail in [25] (Eq. (5) therein).
The kinetic equations follow non-equilibrium evolution of the emerging free-electron and ion densities due to photoexcitations, Auger decays of core holes and the subsequent electronic collisional processes such as elastic electron-ion collision, electron impact ionization and three-body recombination. The respective collision terms, Ω (e) and Ω (i,j) , then describe the change of the electron and ion densities, respectively, due to: (i) the creation of the secondary electrons and highly charged ions via photo-and collisional ionizations of atoms and ions, and Auger decays of core holes, (ii) elastic and inelastic collisions of electrons and ions, (iii) recombination processes, and (iv) short range electronelectron scattering. The cross sections and rates of those processes are derived in the atomistic approximation including the interaction of isolated atoms with impact particles, and implemented within the two and three body Boltzmann collision integrals (for details see, e.g., [3,20]). The cross sections and rates of atomic processes induced by X-ray photons are obtained with the XATOM code [10,18] based on Hartree-Fock-Slater scheme. The impact ionization cross sections (and the respective recombination rates) are calculated with Lotz formulas [13]. The short range electron-electron scattering are modeled with the Fokker-Planck collision integral [20]. The Pauli blocking is not included, as the electron system is assumed to be classical in this model.
Details on the collision terms can be found in Refs. [25]. If the collision terms are neglected, Boltzmann equations, Eqs. (1-2), reduce to the Vlasov equation [20] describing the evolution of a collisionless plasma. Initial input of Eqs. (1-2) is given by atomic density function, which represents the neutral sample.
Both the collisionless part of the equations and the collision terms conserve per construction the number of particles and the total energy in the system. This has been proven also for the numerical algorithms applied [25]. The numerical conservation of particle number and energy have been checked by dedicated tests in all numerical implementations of the Boltzmann Eqs.
The results obtained with the first ('core') version of the Boltzmann equation solver were published in [25]. The sample studied there was a noble-gas cluster irradiated with VUV/XUV radiation. The system of equations contained only ground-state ion configurations as the incoming photons excited electrons from valence levels only. The inverse bremsstrahlung process was also included as it contributed to the electron heating in this photon energy regime. On the femtosecond timescales considered, electron recombination was neglected. Further tests and applications of the code to the similar systems (weakly bonded noble gas clusters) then followed in Refs. [26][27][28]. The observables studied there were usually ion time-of-flight spectra.
An extension of the code was performed in [29], where the three-body recombination (in solids known as the Auger recombination) was added to the code. In [30], the Boltzmann equation solver was successfully applied for the analysis of electron spectra. It followed the evolution of the electron distribution through its nonequilibrium stage toward local thermodynamic equilibrium, without making any 'instantaneous thermalization' assumption, quite frequently used at that time (see, e.g. [8]). This ability encouraged the first application of the code in the context of plasma and WDM studies: in [6] the code was used to investigate femtosecond thermalization of electrons created after X-ray irradiation of hydrogen. In [31], for the first time heterogeneous samples, atomic clusters, containing atoms of two noble gases, were successfully studied.
The Boltzmann equation solver, originally prepared to follow non-equilibrium evolution of finite systems (e.g., atomic clusters), was adapted in [32] to study creation of plasma from bulk materials irradiated by Xrays. Irradiation with high fluence X-ray pulses enables to achieve transiently an almost homogeneous ionization dynamics within a large volume inside the irra-diated solid, provided wide beam focusing and thickness of the target layer comparable to the absorption depth of the X-rays used [23]. When compared with the simulations of finite samples, for which particle distributions in full phase space should be followed, the spatially-uniform ionization dynamics within the bulk material allows for a significant simplification of kinetic equations. One can then assume that both electron and ion distributions are approximately uniform within the bulk, which implies that the evolution of these distributions is position-independent: Consequently, the quasineutrality condition is naturally imposed on the sample, as at each space point the number of created electrons and ions is identical. The net electrostatic force responsible for mutual interactions of electrons and ions is then equal to zero. The equations still account for the fast electron thermalization, as the short-range electron-electron interactions remain unaffected.
Further, any directionality imposed by the initial photoionization process is quickly lost, due to a large number of collisional ionization events with isotropic emission of secondary electrons, occurring soon after Xray exposure begins. In this way, one can with a good accuracy restrict to the isotropic component of electron and ion distributions [25]. The evolution of the distributions only in momentum space is then followed. These simplifications, as well as the net electrostatic force equal to zero, extensively speed up calculations, as no adaptive stability condition for the evolution in real and momentum space, restricting computational timesteps, is then necessary. The simplified equations can then be conveniently applied to analyze ionization dynamics within a bulk material after its X-ray irradiation (for further details see [32]).
In parallel, the extension of the code applicability to hard X-ray regime was performed. To circumvent the 'bottleneck' of very high number of active configurations involved in the excitation and relaxation of an Xray-irradiated sample, in [32], an alternative 'predominant excitation and relaxation path' (PERP) approach was proposed. It still uses true atomic configurations but limits their number by restricting the sample relaxation to the predominant relaxation paths, determined by largest cross sections and transitions rates. The current scheme includes the most probable photoionization and the most probable Auger decay from each configuration within the path. Consistently with the treatment of only the predominant photoinduced processes, we included into the code the predominant collisional ionization processes (from the outermost shell of all considered atoms and ions) and the corresponding three-body recombination rates. Let us emphasize that including collisional processes extends the number of active configurations by those that can be obtained from each 'photoinduced' atomic configuration by a sequence of collisional ionization or recombination from outermost atomic shells. This forms a ladder of additional active configurations with charges from 0 to the highest one allowed in the system for each active configuration obtained through a photo-or Auger ionization. The reliability of the scheme was tested by performing respective calculations for a bulk material consisting of light atoms (carbon) and comparing their results with a full calculation including all relaxation paths. Later, these results were also benchmarked with an independent Molecular Dynamics calculations and found to be in a good agreement [1]. The selection of the predominant excitation and relaxation paths is now performed automatically by a dedicated code. This allows to apply the approach also to heavy elements, and at any X-ray photon energies.

Application of Boltzmann equations to follow evolution of X-ray-irradiated Au layer into warm-dense-matter state
The Boltzmann equation solver was originally designed to follow ultrafast electron and ion dynamics on a few 100 fs timescale. As ion masses were much larger than electron masses, the recoil momenta of the ions during electron-ion collisions could be were neglected on those ultrashort timescales. Therefore, the ions remained 'cold', i.e., they retained their initial kinetic energy throughout the entire (a few 100 fs long) simulation.
However, on longer simulation timescales the energy exchange between ions and electrons and the resulting electron-ion thermalization cannot be neglected. For the purpose of the current project, the original Boltzmann equation solver has been extended so as to also include the electron-ion coupling, modeled here with the Spitzer rate [20]. With this extension, the code can now be used on picosecond timescales. However, as the electron-ion coupling in bulk Au is weak, it does not visibly affect the sample evolution on a picosecond timescale. It should be also emphasized that the plasma code per construction does not include band structure. It describes a sample as an ensemble of unbound atoms, and models their interaction with X-rays (including subsequent Auger decays), and with the Auger and secondary electrons, using atomistic cross sections and rates.
However, due to the presence of numerous electrons and ions in the dense plasma, the ions cannot be treated any longer as isolated particles. In order to account for this effect, we introduce the respective continuum lowering [15], which can be estimated with the code XATOM for plasmas within a wide regime of electron density and temperature [10,19]. We then calculate the excitation energy for each level in respect to the lowered continuum. This yields 6s electrons to be delocalized (free) in ions already under the ambient conditions (see [16], Fig. 3a therein) and after X-ray excitation. Consequently, the 6s level is then the limiting energy level splitting the domains of localized (bound) and delocalized (free) electrons. The initial state of the irradiated sample, the neutral Au, is then an 'ensemble' composed of atomic Au +1 ions with delocalized 6s electrons as indicated by experiments [16].
The code can treat atomic Au ions emerging after X-ray excitation up to a charge of +9. The number of photoinduced transition between 144 various ion configurations of Au considered in the PERP approach [32] is ∼140. As mentioned earlier, their cross sections and rates were obtained with the XATOM code [10,18].

Results for electronic and optical properties of X-ray irradiated gold
We applied our code to a test case, using a set of pulse parameters available at the FLASH experimental facility for high-energy-density experiments. We assumed that X-ray pulses with photon energy of 245 eV and FWHM pulse duration of 60 fs irradiated a 30-nm-thick layer of Au. The thickness of the Au layer is comparable with the attenuation length of 245 eV photons in bulk gold (∼37 nm) which prevents forming a strongly nonuniform distribution of X-ray energy absorbed within the layer. Pulse intensity (of a Gaussian temporal profile) was adjusted so as to yield an average dose of ∼3 eV/atom which corresponded to a deposited dose of 1.6 MJ/kg.
During the irradiation of the Au layer with 245 eV X-ray photons, 4f-shell electrons are predominantly excited, leaving behind 4f holes with lifetime of ∼6.6 fs. The 4f core hole is later filled with an electron, most probably from 5d shell, and another 5d Auger electron is emitted. Further photoionization of the 4f core hole configuration is also possible to a double core hole state. Photoinduced transitions occur further but at the same time the released electrons start to ionize Au ions through impact ionization or to recombine with the ions. The complex ionization dynamics leads to a fast increase of the average ion charge and fast energy exchange within the electronic system, establishing a local thermodynamic equilibration on tens femtosecond timescale, as it will be shown below.
Under the current X-ray pulse conditions, the model predicts a presence of atomic ions of maximal charge +4. Au +1 ions build the ground-state configuration (with delocalized 6s electrons). For Au +2 ions, we have 2 active configurations, for Au +3 and Au +4 ions we have 3 active configurations each. The total number of active configurations with charges up to +4 is 9, and the number of the respective photoinduced transitions involved is 8. Within the PERP approximation used in the Boltzmann equation solver, established for the current experimental conditions, subdominant rates for other transitions were consistently neglected. The list of the selected active configurations (of charge up to +4) is given in Appendix A. In Fig. 1, we show the corresponding average ionization degree of the bulk gold (per atom, Fig. 1a) and the relative charge content, i.e., the number of specific Au ions (summing contributions from all active ionic configurations) normalized to the initial total number of Au atoms (Fig. 1b) as a function of time. Please note that in this picture already the 'neutral' bulk contains only ions, Au +1, which are not included in the calculation of < Z >. We show our predictions on timescales up to 200 fs after the maximum intensity of FEL pulse (i.e., time zero), when ionization process is rapid. At later times, there is not much change in the plotted quantities.
The average charge < Z > per atom increases shortly after the exposure start (Fig. 1a). It saturates at around 100 fs to increase later slightly due to electron production in long-timescale Auger decays. At longer times, also three body recombination starts to play a role, making the relaxation dynamics even more complex. The largest relative ion content is observed for Au +1 with one (delocalized) 6s electron (neutral bulk Au). The ion content for Au+2 with delocalized 6s electron is the second largest. Au+2 charge states can be achieved either after a core hole excitation 4f or 5p, or a photoor impact ionization from level 5d. The ion content Au +3 constitutes only a small fraction of the overall number of ions and originates either from Auger decay of core hole (here 4f -5d 5d), or from a photoionization or the impact ionization of Au +2 ion. Au +4 ions are created in a similar way (see Appendix B).
In Fig. 2, we show the prediction on the kinetic electron temperature and transient electron-ion collision time, τ , calculated selfconsistently from the actual electron-ion collision rates in the sample, involving transient electron and ion distributions. The kinetic electron temperature is obtained as 2/3 of the total kinetic energy of all free electrons above the 6s level (calculated in respect to the 6s level) divided by the number of free electrons above the 6s level. The factor of 2/3 originates from a standard definition of kinetic temperature in 3D hot electron gas: E kinetic = 3/2 kT kinetic . When electron gas becomes thermalized, kinetic temperature equals to the Maxwell-Boltzmann temperature. These predictions do not include the contribution from the delocalized 6s electrons, which cannot be extracted within the framework of our model.
Kinetic electron temperature corresponds initially to the temperature of the emitted photoelectrons (∼93 eV) after the excitation of a 4f electron (Fig. 2a). Due to fast energy exchange among numerous secondary electrons, the electron distribution quickly thermalizes, and the system enters a local thermodynamic equilibrium regime at ∼50 fs, when kinetic temperature almost stops to change. This thermalization is also reflected in the predicted electron-ion collision time defined as an inverse of the total electron-ion collision rate normalized per one electron (derived in the same way as described in Section 2.3 of [4]). The calculation of the total electron-ion collision rate takes into account the total collisional cross section, including the scattering of hot electrons on ions in any possible atomic configuration included in the model, i.e., all free electronbound electron collisions. The collision time strongly increases, following the decrease of the kinetic temperature. It peaks around −35 fs, i.e., at the same time, when the kinetic electron temperature reaches its minimum. Later, after the thermalization of free electrons, it saturates at the value of ∼ 1.2 fs, in a very good agreement with the equilibrium value measured in Ref. [16] after relaxation of optically excited Au. In Fig. 3, we show the evolution of transient freeelectron-energy distribution, n e (E) (normalized per electron) as a function of energy. Snapshots for the times −90 fs, −50 fs, 0 fs (FEL pulse maximum), 50 fs and 90 fs are shown. The normalized distribution is multiplied by the actual value of the transient ionization degree, < Z >, i.e., the number of free electrons above the 6s level divided by the total number of atoms (Fig. 1a), in order to show the increase of the free electron density per atom in the sample. The transient electron energy distributions are compared with the corresponding Maxwell-Boltzmann (M-B) distributions calculated, using the transient kinetic electron temperature, and multiplied by the value of the transient ionization degree. Initially, the transient electron energy distributions show a significant deviation from the M-B distributions both in the low and in the high energy range, where the contribution of photo-and Auger electron peaks to the spectra at ∼100 eV is visible. After time zero, the latter contribution continuously decreases and the low energy part of the transient spectra approaches the M-B distribution. This shows that thermalization of the electrons progresses toward a local thermodynamics equilibrium. The M-B curves for times 50 fs and 90 fs are almost overlapping, coinciding with the respective transient free-electron-energy distributions. This confirms that the thermalization has been completed by those times.
The free-electron density and electronic temperature can be used to obtain predictions on the transient optical properties of WDM Au. For the conversion, we use the Drude model described in [16], including contributions of interband transitions. The scheme is the following: using known experimental data on gold from [16], we can identify the initial value of the complex dielectric function, ε(ω, t 0 ). Taking the Eq. (1a) and (1b) from [16], we arrive at the expression: As we know the both components of the complex dielectric function ε(ω, t 0 ) (from the experimental data), as well as the Drude terms, ε f 1 (ω, t 0 ), ε f 2 (ω, t 0 ), calculated with the equilibrium value of the collision time and the initial plasma frequency (obtained with 6s electron density), we can calculate ε ib 1 (ω, t 0 ) and ε ib 2 (ω, t 0 ) from Eq. (4). Please note that in this way the interband contributions remain temperature-independent which may be a potential source of error. However, the analysis provided in [16] indicates that the temperature dependence is rather weak under the conditions of the actual simulation.
The Drude model applied here a parameter-free model, as it uses the electron-ion collision time calculated on-the-fly from the Boltzmann equations. For the initial condition, let us recall that within the atomistic picture, the 'cold' Au corresponds to an ensemble of Au+1 ions with the delocalized 6s electrons, i.e., the valence electron density of 'cold' bulk Au (with 1 valence electron per atom) is then equal to the number density of Au.
Below we plot the DC conductivity, σ 0 as a function of time, σ 0 = ω 2 p τ/4π, where ω p is the electron plasma frequency, and τ is the electron-ion collision time (Fig.  4). As our atomistic calculation cannot include the contribution of the delocalized 6s electrons per se, at early times (up to −80 fs) we use an ambient value of τ for gold from [2], and switch to τ from Fig. 2b,   collision time, Fig. 2b), until free electrons thermalize (∼50 fs) (cf. Fig. 2) and then saturates at an equilibrium value of 1.94 × 10 16 1/s which is in a very good agreement with the equilibrium value of ∼ 2 × 10 16 1/s measured in Ref. [16].
The predictions on transient optical properties, reflectivity and transmissivity are shown in Fig. 5. They were obtained with the optical probe of wavelength: (a) 500 nm (2.48 eV), and (b) 900 nm (1.38 eV). The 500 nm pulse probes bulk Au directly above the 'cold' interband transition threshold (∼ 2 eV), where interband contributions dominate. The 900-nm probe pulse restricts to the regime dominated by intraband contributions. Both regimes can be tested with our Drude model which includes both the interband and intraband contributions. In Fig. 5, we show predictions obtained with the Drude model, using the transient electronion collision time, τ, from Fig. 2b (green curves). As in case of DC conductivity, at early times (up to −80 fs) we use an equilibrium value of τ for gold from [2], and switch to τ from Fig. 2b, as soon as the transient electron-ion collision time approaches this equilibrium value.
Generally, both the initial increase or decrease of transmissivity and reflectivity follow the increase of the FEL pulse intensity. The corresponding extrema of optical properties are slightly delayed in respect to the time zero. This is due to the still developing electron cascades at this time. Only after the cascading stops, optical properties stop to change. The predictions show also a strong effect of the finite collision time which is especially pronounced for reflectivity with intraband contributions only (for 900-nm probe pulse). Neglecting collisions induces a reverse behavior of reflectivity, observed also for other probe pulse wavelengths (not shown). Further validations of these predictions require a dedicated time resolved experiment. The respective efforts are in progress.

Conclusions and outlook
The predictions on the evolution of X-ray-irradiated gold obtained with atomistic Boltzmann kinetic equations have been reported. The equations followed complex relaxation path of bulk Au after the impact of 245 eV photons, using the 'predominant excitation and relaxation path' approximation introduced in [32]. The equilibrium values of free-electron density, electron collision time, as well as the DC conductivity obtained after electronic thermalization show good agreement with the measurements from [16] obtained for optically excited gold at the similar absorbed dose. This demonstrates the ability of the Boltzmann equations approach to trace formation of warm dense matter created from materials composed of heavy elements, irradiated with X rays. In addition, the self-consistent parameter-free Drude model, including interband transitions and transient collision time calculated on-the-fly from the actual state of the Au sample, has been used to obtain predictions on transient optical reflectivity and transmissivity of warm dense gold. These predictions can be validated by dedicated experiments in future.