Equations of state, energy transport and two-temperature hydrodynamic simulations for femtosecond laser irradiated copper and gold

In this work, the model of thermodynamic and transport properties of copper and gold at electron-ion non-equilibrium state is presented. Accepted here ranges of electron temperature and pressure are enough to describe the experimentally achievable states. The changes in electron spectra due to electron heating and compression or expansion are taken into account using two-parabolic model. In our previous works, thermal conductivity and electron-ion coupling were considered as dependencies on electron and ion temperatures. Now the dependence on density for these coefficients is taken into account. To include exchange-correlation effects on electron-electron collisions we have found out how this effect can be included in electron screening. In addition, we have renewed our approach for heat conductivity calculation to include thermoelectric phenomena, which are significant at high electron temperatures. The effect of electron heating on sound velocities in aforementioned metals is investigated. The two-temperature hydrodynamics simulation of film expansion was provided with the use of the model presented here.


Introduction
The processes in metals irradiated by ultrashort femtosecond laser pulses with absorbed energy fluences higher than a few tens of mJ/cm 2 are investigated. For such fluences of laser energy the function of the Fermi distribution for electrons comes to equilibrium within a short period of time 10-100 fs. As the electronic subsystem has a temperature T e is rather greater than the temperature of the ion subsystem, well-known two-temperature hydrodynamic model (2THD) is used to estimate the propagation of heat associated with the diffusion of hot electrons from conduction bands and electron-ion energy exchange. Model coefficients are calculated for the electron thermal conductivity κ(ρ, T e , T i ) and electron-ion heat exchange α(ρ, T e ) with the two-temperature equation of state E(ρ, T e , T i ) and P (ρ, T e , T i ). Similar equations of state for aluminum and tungsten were obtained in [1,2], where various quantum ab initio methods were used. In the recent paper [3], quantum molecular dynamics calculation of the shock Hugoniot and release isentropes for aluminum was carried out in wide ranges of pressures (up to 1 TPa) and temperatures (up to 70 kK).
Simulation using 2THD code shows the transition of the electronic thermal wave from supersonic to subsonic regime of shock wave propagation. At the moment of transition heat wave emits a compression shock wave moving to the deep of metal. Nonlinear evolution of compression wave after its separation from subsonic thermal wave to the cleavage on the back side of the film layer can be seen in 2THD modeling and molecular dynamics simulations (MD).
To this date, none of the previous works has not been fully demonstrated the effect of electron temperature, and in particular, the density of metal on the 2T stage characteristics required for further calculations. In the paper [4] laser femtosecond irradiation of copper with moderate intensities up to 10 15 W/cm 2 was considered. Theoretical model presented in [4] was applied for isochorical electronic heating. In the present work, the choice of parameters for two-parabolic electron density of states was determined by the pseudopotential and FP-LAPW calculations, carried out using the computational codes VASP [5] and Elk [6], and, therefore, does not contain any arbitrary assumptions, which would require additional study.
Also it should be noted that in previous works [7,8] was not properly taken into account the dependence of used in 2THD modeling approach two-temperature electronic properties on the density of the metal. These properties should be considered as direct functions of the density and simultaneously as implicit functions of density due to the changes in the electronic density of states caused by deformation.
It should be pointed out that the results for the thermodynamic properties and transport coefficients used for the 2THD calculations have the form of user-friendly analytical dependences on electron and ion temperatures and density.

The model for calculations
The influence of high hydrostatic tension and compression was reviewed using the isotherms calculated by DFT where the number of free electron states was increasing with electron temperature.

Computational details
Calculation of electron temperature isotherms for internal energy E and pressure P at zero lattice temperature was provided for fcc copper and gold in the range of specific volume from 0.8 to 1.1 using DFT codes VASP [5] and Elk [6]. In the VASP calculations we use PAW approach with exchange-correlation functionals in the form of LDA and PBE. In the full-electron calculations only PBE exchange-correlation functional was used. The cut-off energy of plane waves was fixed at 500 eV, the Monkhorst-Pack grid of k-points had the size 21 × 21 × 21, the number of empty electronic levels was taken not less than 20. When the maximum of the considered electronic temperature T e = 55000 K was reached, the number of empty electronic levels increased to 30. These parameters were selected as a result of the convergence test so that the convergence error is estimated as less than 1%.

Two-parabolic DoS parameters
The data shown in figure 1 indicates that for the noble metals the changes in the electron spectrum caused by both compression and electronic heating are equally important.
Indeed, at the experimental conditions for femtosecond laser irradiation of a metal surface the amplitude of the shock wave has the upper limit of about 100 GPa. In this case there is no need to consider the absolute values of change in volume of 0.1-0.2 greater than the magnitude of the initial volume. Then, as shown in figure 1a for gold, the effects of cold compression and electron heating on the discussed properties are comparable, although they are expressed in different ways. In the case of compression, we are dealing with the expansion of the d-band, and in the case of heating-with its contraction, followed by a shift towards lower energies. Trends discussed in more detail illustrated in figures 1b-1d. The further note is that the question about correctness of DFT calculation of maximum d-band energy is open until this time. It is known that calculated maximum d-band energy for copper and gold is overestimated by 0.4 and  In the former case the value ε 2 was fixed equal to experimental value (−2.5 eV). In the latter case green solid line corresponds to equilibrium density and green dashed line corresponds to volume compression V /V 0 = 0.9. (d) The numbers of valence s-band electrons for gold and copper. The cases of non-compressed and slightly compressed V /V 0 = 0.9 metals are considered. The s-electron number as a function of electron temperature for gold is shown with diamonds and s-electron numbers for copper are shown by lines with circles. 0.8 eV, respectively, in comparison with experiment [9,10]. From previous studies on the twotemperature electronic properties has been known that the characteristics of copper and gold show a qualitatively similar behavior. It was decided for copper, since in this case overestimation in DFT is not so great, to adjust the data for maximum d-band energy at different densities and electron temperature by shifting on −0.4 eV to the lower energies. In the case of gold under normal conditions the magnitude of the discrepancy of experimental and calculated data in 2 times higher, it was decided to fix the ε 2 energy equal to −2.5 eV regardless of density and electron temperature.

Two-temperature thermodynamics
The results of DFT calculations of electron heat capacities and electron pressures initiated by electron heating for copper and gold are shown on figure 2. Thermal capacity of gold is slightly higher than one for copper but, in the case of electron pressure, the situation is opposite. In the work [11] the role of delocalized electrons was discussed in gold and the total electron pressure obtained in [11] is in quantitative agreement with our data. Both electron heat capacities and electron pressures almost not depend on density in the range from V /V 0 = 0.9 to V /V 0 = 1.1 at electron temperatures up to 50000 K for the discussed metals.
In addition to calculations of the response of a two-temperature electron thermodynamics on hydrostatic compression/expansion, a study was conducted to investigate a response of electron thermodynamics to uniaxial compression along the [100] direction at different electron temperatures. This calculation is carried out using VASP, allowed us to obtain data on the behavior of the diagonal components of the stress tensor under uniaxial deformation in the aforementioned direction and when the electron subsystem is significantly heated. The purpose of this calculation was to obtain the dependence of the longitudinal sound velocity on density the result of phonon spectra computation provided by Recoules [12]. and electron temperature. This dependence is obtained simply by the definition of longitudinal sound velocity c 2 l through the square root of the ratio of the elastic modulus C 11 and density ρ: The elastic modulus C 11 is found from the expression for the diagonal stress-tensor component σ xx , where x is the direction of deformation and V -volume: According to Lindemann [13], the melting temperature of solids may be defined as the temperature at which the ion lattice vibrations reach some known in advance value for each material. According to Ubellohde [14] this amplitude oscillations x L of ions divided to the distance between the atoms of the lattice, is equal to 0.1 in the case of gold. Rewriting the expression for the Lindemann criterion in an convenient form to use obtained data on the dependence of the elastic modulus C 11 on density and electron temperature, we can write: In this formula, k B -Boltzmann constant, M I -mass of metal ion in lattice (6). Then we can obtain the dependence of the melting temperature of gold on electronic temperature. This evaluation is shown in figure. 3b. It is seen that the data obtained by the method described above in satisfactory agreement with the data calculated by V. Recoules [12], wherein the melting temperature was also evaluated by the Lindemann criterion, but the dependence of longitudinal sound velocity was obtained on DFT calculations of the phonon spectrum for gold.

Transport properties
With regard to the calculations of the electronic kinetic coefficients based on the solution of the kinetic Boltzmann equation in the relaxation time approximation (electronic conductivity) and the Kaganov, Lifshitz and Tanatarov theory (electron-ion heat exchange) several changes were made in the calculation scheme presented in [15]: (i) Fully integrated effect of density and electron temperature on the electronic density of states in two-parabolic approximation. This means that since the definition of energy and the occupation numbers of electron valence bands at all the following stages of the calculation of the kinetic coefficients changes in the spectrum are taken into account. (ii) The effect of screening in electron-ion and electron-electron interactions on the electronic thermal conductivity and the electron-phonon coupling is investigated for noble metals. Based on this data, it was found that among the four considered types of screening (Thomas-Fermi, Lindhard, Hubbard, Singwi-Sjoelander) the latter is the most correct choice of screening [16]. In this type of screening the effects of exchange and correlation in collisions of electrons with electrons and ions are taken into account. In the case of electron-electron collisions Hubbard and Singwi-Sjoelander dielectric functions take the form ε ee = 1 + (ε L − 1)(1 − G q ), in the case of electron-ion collisions in the same cases we have Here, ε L is the dielectric function in random phase approximation (Lindhard). The functions G q is special factors which gives the form of correction to the initial RPA dielectric function. (iii) The effect of thermoelectric phenomena on the contribution in electronic thermoresistivity due to electron-electron collisions has been evaluated in this work. It was already noted [7] that used up to now our method for calculation of electronic thermoresistivity is valid neglecting the contribution of the electron-electron collisions in resistivity that is incorrect for d-metals at high electron temperatures. In this paper, evaluation of thermoelectric effects was based on the definition of electronic coefficients through the matrix of the Onsager coefficients in the presence of such thermodynamic forces as the electric field and the electron temperature gradient. It was found that all Onsager coefficients can be obtained by using previous expression [7] of the electron-electron effect in the electron heat conductivity. Here, v 2 , τ −1 and µ are average square of s-electron velocity, effective rate of electron collisions and chemical potential, correspondingly. Discussed here electron heat conductivity can be written in the form κ In accordance with these changes the partial contributions of s-s and s-d electronic collisions were calculated. It is convenient to represent these contributions in the form of reciprocal value which is thermal resistance.
Dependence of the contribution of electron-electron collisions in the thermal resistance for gold or copper is qualitatively behaves the same when the densities and parameters of twoparabolic model are in the accepted ranges. For electron temperatures below about 2 eV the discussed contribution in the electronic thermorestivity shows sharp and quasilinear growth that is followed by the slow decline in the remainder range of considered electron temperatures. It is worth to note that in the case of gold, the influence of the d-band maximum energy correction is more significant than the influence of density. In the case of copper, the situation changes noticeably and the effect of density change and the effect of shift of maximum d-band energy have the same importance. Also we see that the effect of density changes is relatively small for any electronic temperatures and any considered metal, so the electron-electron contribution to thermoresistivity for these metals can be considered as independent on density. Also some remarks about the changes in the previously used model for evaluation of electronion collision effect on thermoresistivity should be made. In the previously used model [15] our calculation of this effect was based on experimental data for electroresistivity of considered metals at atmosphere pressure.
In the sense of describing the dependence of this contribution on the ion temperature scheme has not changed, but the available experimental data were extrapolated to treat the considered range of densities for each metal. The initial dependencies of the contribution of electron-ion collisions in thermoresistivity and the effective frequency of electron-ion collisions on the ion temperature were obtained using Drude relations: Here, ν ei -the effective frequency of electron-ion collisions, C s -heat capacity of s-electron gas, v 2 s -average square of s-electron velocity, e 2 -square of electron charge (6), n s and m s -selectron density and mass, respectively. To obtain the electrical resistivity S ie as function of equilibrium temperature the known experimental data describing the state of metal up to several 1000 K [17] were fitted and extrapolated using the Mott's limit In this expression the following notation is used: L is the number of atoms in a cube with edges which are equal to the lattice parameter, a is a constant, selected using the data of quantum molecular dynamics calculations based on Kubo-Greenwood theory. This constant is about 0.2 for gold and copper. The Planck constant and electron temperature as well density dependent number of free s-electrons per atom z s are also used (6).
On a step of extrapolation to treat the considered ranges of densities the formula (7) is preserved, but the low-temperature and high-temperature terms are differently determined as functions of density. For the low-temperature limit, using the expressiion of the electrical conductivity by the Drude formula where l-mean free path of the s-electron, V F -Fermi velocity of s-electron, one can obtain the following estimation for the dependence of the conductivity on density: In this derivation the relations were used between the path length l = (nσ) −1 and the scattering cross section σ, and between the scattering cross section and the mean square amplitude of phonon thermal vibrations Σ ∼ u 2 T . Also the definitions of the Debye temperature T D = ηω D /k B = η/k B ζ/M i and the Grüneisen parameter G = d ln T D /d ln V are used here.
According to this estimate, the electrical conductivity as function of density was obtained by using the extrapolation formula: Here, in the second bracket the numbers of s-electrons at the fixed electron temperature 1000 K are taken. Used in this paper dependence of Grüneisen parameter on volume for gold and copper was obtained [18].
In the case of high temperatures we consider Mott's limit (7) as a function of volume. The dependencies of the concentration of atoms and the number of s-electrons per atom on the density are taken into account. In this case, used numbers of s-electron per atom were taken at the maximum electronic temperature which is 55000 K. The number of s-electrons per atom at given electronic temperature can be considered as upper limit of this number in our calculation.
As can be seen from figure 5, for each of the considered configurations, the contribution of electron-ion collisions in thermoresistance monotonically decreases with electron temperature, going on close to high-temperature asymptotic, which is already near 20 000 K. It is also worth noting that similar to the electron-electron contribution shown in figure 4, in this case the dependence on the density is also sufficiently weak. Therefore, both the contributions of two types of collisions and the total resistance are weakly changed with compression and strain achieved in experiments, but the influence of electronic temperature is significant. This statement is illustrated in figure 6, according to which total electron thermal conductivity of gold, which was obtained according to the Mathiessen rule, the influence of two-parabolic approximation parameters is more appreciable than the effect of compression.   [19,20].
Similar to the changes of the computational scheme for electron heat conductivity, in the calculation of the electron-ion coupling the evolution of the electronic structure and the screening of the electron-ion interaction are also considered. As a result the electron-ion coupling as function of electron temperature at different densities is obtained and shown in figure 7. It can be noticed that in contrast to the electron thermal conductivity, the change of density significantly affects the value and even behavior electron-ion coupling. Indeed as shown in figure  7 at low densities the electron-ion coupling of copper has a steeper growth than in the case of compression. It is also important that the values of electron-ion coupling for gold and copper at low electron temperature are in good agreement with the available experimental data [21,22].

Statement of the problem
The final stage of our study was the two-temperature hydrodynamic modeling of the electronion relaxation for thin 70 nm copper foil irradiated by femtosecond laser pulse which has the parameters, similar the experimental data [10]. Indeed, the used initial parameters for modeling is in accordance with performed experiment: the duration of pump pulse is 150 fs, wavelength is chosen equal to 1060 nm and the absorbed laser fluence is 160 mJ/cm 2 . The latter value is taken as the lower limit of the energy deposition in the experiment [10], because there were serious doubts in so small electron temperature at the peak of the electron heating. The calculation was based on the model of two-temperature hydrodynamics introduced by Anisimov [23] in its modern formulation for the simulation of femtosecond laser heating of metals [7]. The similar model established on two-temperature equations was used in the paper [4] where the loss of ionic subsystem energy due to hydrodynamics motion initiated by pressure action has been neglected. Considered in addition to the model considered here, the model presented in [15], the model of electron thermodynamics and the electron-ion coupling obtained by Zhigilei [24] in conjunction with the well-known phenomenological expression for the electron thermal conductivity "5/4" [25]. In this new model, two assumptions were made about the depth of the skin layer at the fixed heating laser wavelength were considered. In first case this depth was 14 nm in accordance with reference data [9]. In the latter case the suggestion where the skin depth is equal to foil thickness is considered [22]. Such suggestion is based on the allowable possibility of superfast fjil heating by excited electrons in the early moment of thermalization when the lattice has a temperature about 300 K may occur. Such superfast heating can take place only if the path length is about several tens of nanometers.
In the simulation it was found that for the calculation with the normal thickness of the skin layer, the shape of the electron temperature profile within the depth of foil is retained during a few picoseconds and creates the similar ion temperature profile which is shown in figure 8. Such similarity is also caused by weak dependence of electron-ion coupling on electron temperature.
In addition, it was found out that even with the equality of the skin layer to the thickness of the foil the temperature obtained on the surface in our modeling and the temperature obtained by the authors [6] are not in agreement. This fact is illustrated in figure 9, where for two initial moments of foil evolution the electron temperatures of 2THD simulation provided in [10] are shown by the horizontal lines intersects at 10 nm depth from the front surface the profiles of the    [10], which shows in the original paper the dependence from time without any reference to spatial localization, is based on the fact that the authors of the work [10] consider the simplified 2THD model where the term for energy transfer by heat conductivity is not introduced.

Results
Dynamics of two-temperature relaxation using mentioned above three models for thermodynamic and kinetic characteristics at the two-temperature state is shown in figure 10. The left figure shows a comparison of the present model (solid line) with the model used by the authors previously (dash-dotted line) and with the model based on known from the literature data (dashed line). The electron and ion temperatures are given for the front surface of the foil during the first 15 ps. In all presented here calculations, the skin layer is taken equal to the thickness of the foil. Note that in many ways the serious difference between our old and new models results is due to the accounting for the first time dependence of electron two-parabolic density of states on density of metal. In addition, it is worth noting that the use of the model with the data given from the papers [24,25] is the reason of observed qualitatively different from the other cases the behavior when ion temperature exceeds electron temperature on reaching the 7 ps that would cause a reverse ion-electron heat. On the right side of figure 10 it is shown that in spite of the strongly increased the initial electron temperature on the surface, thermalization of two subsystems is achieved much faster if the skin depth is equal to reference data [9] than in the case of overestimated skin depth. Noticeable oscillations with time of electron temperature on the right side of this figure in the early moments of thermalization is based within the changes of taken into account changes of electron thermodynamics and kinetics with density.
Based on the results presented in figures 9 and 10, the authors conclude that shown in [10] dependencies of electron temperatures on the time can not be considered as the temperatures at any local point of the foil. The first reason is fact that in [10] investigated absorption of x-ray pulses with a wavelength of few nanometers so one can obtain a value of 125 nm for the extinction length. With this value of the extinction length x-ray laser pulse is passed through the foil only slightly losing its intensity. The collected data is accumulated information by pulse The electron and ion temperatures as functions of time calculated with the use of presented here model for two cases. In the first case skin depth equals to the thickness of submicron nanofilm (solid black lines) and the skin depth equals to 14 nm which is found using Palik book [9] was used in the second case (dashed green lines).
which is passing through the entire thickness of the foil, and thus provides an general picture of the temperature in the foil. Secondly, during the time of simulation the authors of [10] did not take into calculation mechanism of electron heat conductivity, and that is generally incorrect because, as seen from figure 10, for the more homogeneous distribution of electron temperature the greater thermalization time is required. However, according to the available experimental data [26,27], the duration of two-temperature thermalization in gold is about 5-10 ps. In copper such thermalization should be faster due to the greater electron-ion coupling.
Based on the assumption that the extinction length for probe x-ray laser pulse is not depend on density and obtained experimentally [10] electronic temperature T exp is integral value for the foil at fixed time the simplest method using the calculated profiles T calc to represent such an integral value can be described by the relation: Here, d-the thickness of the plate, l ext -the extinction length. In this approach, we have ignored the fact that the density of the foil with time and depth may vary. A more accurate interpretation of the experimental data may be performed using a known [28] mass attenuation coefficient µ, when considered the fact that the optical path is dependent on the local density 15 T e , kK Figure 11. Integral electronic temperature of copper foil irradiated by femtosecond laser and heated as a function of time. Experimental data are shown by points with error bars and obtained in the paper [10]. The result of the 2THD data evaluation using expression (8) is presented by blue circles, and the similar data obtained by the formula (8) are shown by red diamonds. The solid line correspond to the 2THD modeling performed by authors of paper [10].
As a result, the convolutions of calculated data from formulas (8) and (9) give the possibilities to carry out a valid comparison with the experimental results and 2THD modeling are made in [10].
As can be seen from figure 11, the dependence of the average electron temperature in the proposed interpretation on the basis of the relations (7) and (8) is monotone. The difference between these dependencies is only noticeable at the beginning of the two-temperature stage, when the electron heating, and as a result, the electron pressure leads to significant changes in density. It is worth noting that the average temperature of the foil with the changes in density in better agreement with the experimental values. On the other hand, data averaging according to the relation (7) in the first 6 ps are consistent with the results of calculation in [10], and then reduced with respect to the cited paper. According to the authors, the most likely explanation is related to the fact that the two-temperature model of the dynamics used in [10] does not consider the process of energy transfer by the action of pressure in foil. In other words the used in the paper [10] suggestion that electron-ion relaxation is isochoric may be valid according to our explanation only during some picoseconds after irradiation.

Conclusions
In conclusion of this work, it should be noted that the extended model of two-temperature electron characteristics is obtained. Two-temperature equations of state for gold and copper are evaluated with the calculations of electron heat conductivity and electron-ion coupling. For the first time dependence of these coefficients on the density is fully taken into account, including changes in the electronic density of states with density. The developed model was tested by the simulation of the experiment described in the work [10], and showed a satisfactory agreement with the results of the experiment in the cited work.