Pressure in electronically excited warm dense metals

Non‐equilibrium two‐temperature warm dense metals consist of the ion subsystem that is subjected to structural transitions and involved in the mass transfer, and the electron subsystem that in various pulsed experiments absorbs energy and then evolves together with ions to equilibrium. Definition of pressure in such non‐equilibrium systems causes certain controversy. In this work we make an attempt to clarify this definition that is vital for proper description of the whole relaxation process. Using the density functional theory we analyze on examples of Al and Au electronic pressure components in warm dense metals. Appealing to the Fermi gas model we elucidate a way to find a number of free delocalized electrons in warm dense metals. (© 2015 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)


Introduction
Irradiation of solids with ultrashort laser pulses opened an exciting field of research (e.g. [1,2]). New emerging physics is connected with formation of warm dense matter (WDM) at the initial transient state of material evolution after energy deposition into electron subsystem. WDM is important not only for laser-matter interaction physics but in a wider context as well [3]. WDM in such ultafast phenomena is a non-equilibrium state that makes it very challenging for theory, modelling and simulation. Usually WDM can be described as a two-temperature (2T) system when electron and ion subsystems can be considered in quasi-equilibrium at T e > T i [4].
Physics of WDM is not a subject of pure fundamental interest. One of the major phenomena where WDM properties are crucially important in modelling and simulation is laser ablation [5]. Laser ablation is a multiscale phenomenon. Subpicosecond laser excitation transforms material under normal conditions into 2T-WDM in the isochoric way. Relaxation takes about tens of picoseconds before T i and T e become equal. This relaxation stage governs the details of ablation mechanism and is the focus of the on-going modelling effort. It can be described at the continuum level (e.g. [6][7][8][9]). However atomistic modelling gives the possibility to capture a richer spectrum of structural transitions and nucleation effects (e.g. [10][11][12][13][14][15][16][17][18]).
Both in continuum and in atomistic models it is assumed that the quasi-equilibrium 2T-WDM can be described using thermodynamic concepts. The free energy is represented as where E b is the binding energy, F fe is the free energy of free (ionized, delocalized) electrons (that vanishes at T e = 0), F i is the free energy of the ion subsystem. The corresponding representation of pressure is Usually the first term is called the cold or binding pressure and depends on density P b = P b (ρ), the second term is the kinetic electronic pressure P fe = P fe (n fe , T e ) and the third term is the thermal ionic pressure P i = P i (ρ, T i ).
It is the total pressure P that defines the mass transfer in continuum 2T-WDM models (e.g. [9]). In the atomistic 2T-WDM models sometimes the so-called blast force [6,19] is introduced in order to describe the influence of kinetic pressure of electrons on the ion subsystem [6,11]. The kinetic pressure of electrons is assumed to be a manifestation of free electrons. The difficulty connected with P fe calculations stems from the fact that for this purpose one needs the number density of free electrons n fe (i.e. the degree of ionization). The deployment of the finite-temperature Thomas-Fermi model is one way to overcome this problem [20]. Among the unusual properties of 2T-WDM is the change in lattice stability of some solids at high levels of electronic subsystem excitation [21]. The bond hardening effect for gold was initially predicted by finite temperature density functional theory (FT DFT) calculations in [21] and got the experimental support in [22]. This effect is not pronounced in simple metals like Al, but is quite high in d-metals like Au. It stems from the renormalization of the electronic structure at T e ∼ T F (T F is the Fermi temperature) and results in the dependence of the interatomic interaction in 2T-WDM material on electronic temperature T e (not only for metals [21,23,24]). Therefore T e changes mechanical, thermodynamic and other properties of solid, i.e. the equation of state (EoS). In [25] the dependence of the interatomic potential for tungsten on T e was considered. In [12][13][14] the electronictemperature-dependent (ETD) interatomic potential for gold in the framework of the embedded atom method (EAM) [26] was developed as well as a numerical scheme that incorporated this ETD-potential into atomistic modelling techniques [27]. This ETD EAM model was successfully used in the study of structural dynamics of laser-irradiated gold nanofilms [15].
The bond hardening manifests itself as a significant pressure build-up in the electronically exited metal. In [13] the comparison of the pressure build-up in ab initio calculations based on FT-DFT P DF T with the pressure buildup in the ETD EAM potential model (based on the same ab initio results) P EAM allowed to propose an ad hoc separation of the pressure increase due to localized and due to delocalized electrons in 2T-WDM gold. In the atomistic model for ablation of gold [12][13][14] the pressure build-up due to localized electrons was taken into account by the ETD EAM potential that effectively describes the dependence of the binding pressure on electron excitation P b = P b (ρ; T e ). In [13,14] the pressure build-up due to delocalized electrons in 2T-WDM gold was taken into account using the blast force fitted to match the difference (P DF T − P EAM ).
In this work we discuss two connected questions. Despite its fundamental significance, the calculation of pressure in the general context of atomistic modelling and simulation is a subject of on-going research (e.g. [28,29]). On examples of aluminum and gold we are making an attempt to analyze the electronic contribution to the total pressure in 2T-WDM metals and to clarify the representation (1). Another question is the separation of electrons into bound and free (or localized and delocalized) that is a general problem for non-ideal plasma physics (e.g. [30][31][32]) including EoS models (e.g. [33]), ion-ion correlations [34,35], continuum lowering [36] and plasma-plasma phase transition [37], for models of laser-matter interactions (see above), models of swift ion tracks (e.g. [38][39][40]) and other high energy density physics applications.

Finite temperature Kohn-Sham DFT for 2T-WDM with fcc ionic lattice
Quantum mechanical description of electrons in condensed phase is a general problem of the highest importance. However the rigorous wave-functions based theory is very complicated even at T e = 0, the finite temperature case being much farther from applicability in practice. Only recently accurate quantum Monte Carlo approaches have been developed in this field [41]. That is why the Kohn-Sham (KS) DFT method in the FT formulation became a tool of choice (for some first examples of its applications for metals e.g. see [42,43]).
In the FT KS DFT framework the free energy of electron subsystem in the external potential of ions is given as where E k is the kinetic energy term, E H is the Hartree term, E xc is the exchange-correlation term, E ei is the electron-ion interaction term, E ii is the ion-ion interaction term and S is the entropy of the non-interacting electrons. KS electronic states are populated according to the Fermi- In this work we perform FT KS DFT calculations in the plane-wave basis using VASP [44][45][46][47] and ABINIT [48,49]. The electron-ion interaction for Al and Au is described by the projected-augmented wave (PAW) potentials [50][51][52]. The LDA approximation is used for E xc [53].
For fcc Al and fcc Au calculations are performed with the lattice constants a = 3.98Å and a = 4.06Å respectively. The k-mesh density and the plane-wave basis energy cutoff have been varied to converge pressure to 0.01 GPa. The 37x37x37 Monkhorst-Pack k-points mesh is used to sample the Brillouin zone. The plane-wave basis energy cutoff is set to 700 eV. The P e (T e = 0) values for different PAW models are Al VASP 3e -0.07 kbar and 9e 0.49 kbar, ABINIT 3e -6.6 kbar and Au VASP 11e -0.61 kbar and 17e -0.04 kbar, ABINIT 11e 24.7 kbar.
www.cpp-journal.org Fig. 1 Components of electronic pressure Pe as a function of Te from the FT KS DFT for 2T-WDM of aluminum and gold with fcc ionic lattice calculated with VASP (3e and 9e PAW potentials for Al and 11e and 17e PAW potentials for Au) and ABINIT (3e PAW for Al and 11e PAW for Au). Dashed lines show the 9e PAW data for Al and the 17e PAW data for Au. Color filling between the curves for different PAWs is used to guide the eye. Fig. 1 shows the dependencies on the electronic temperature of the electronic pressure P e components for fcc aluminum and gold (similar but not extensive data have been reported previously in [54]). The VASP results illustrate in absolute values how pressure components change when one switches from the PAW models with the minimum number of valence electrons (3e for Al and 11e for Au) to the PAW models that treat deeper electron levels explicitly. The Ewald, non-local and xc components go down as the number of electrons in the PAW model increases that is compensated by the rise of all other components. We see the strongest dependence on T e in the kinetic and non-local pressure components. The explanation is that only kinetic and non-local pressure components depend on KS wavefunctions and their partial occupancy values f ( i ). It is the kinetic and non-local pressure components in the case of Au that manifest the largest differences between VASP and ABINIT results (probably due to the differences in the PAW models). One can notice that the Hartree pressure component in Au depends on T e stronger that in Al. This fact implies that the electronic density redistribution inside the fcc unit cell is more pronounced in the case of Au than for Al. In order to analyze these density redistributions we obtain (using VASP) the radial dependencies for the total electron density and for the density that can be attributed to different angular components via the projected density of states (DOS) calculations (Fig. 2). This projected DOS analysis provides only qualitative results because i) the KS wavefunctions inside PAW spheres play an auxiliary role and ii) it is not possible to separate unambiguously the KS wavefunctions into s-, p-and d-components. However some conclusions can be made. In fcc Al the rise of electronic temperature does not essentially change the total electron density. But in fcc Au there is a significant increase of the electron density in the interionic region of the unit cell. The bond hardening effect in fcc Au [21] can be deciphered as a transition of electrons from d-like states into s-and p-like states as well as into the plane-wave like (i.e. free) states. The r s values illustrate that the electron gas in the interionic region of the unit cell is non-ideal in the 2T-WDM conditions considered.

Fermi gas pressure as a tool for revealing free electrons
The concept of the homogeneous electron gas is a starting point of the DFT-like theories of inhomogeneous systems (e.g. [55]). Conduction band electrons in metals are usually considered to be close to the homogeneous electron gas with negligible electron-electron interaction, the so-called Fermi gas (FG). Here we expect that the comparison of the FG model with the FT KS DFT results could allow us to understand the 2T-WDM properties better. The FG model is a simple quantum theory that have one input parameter n fe . If the free electron density is known one can calculate all thermodynamic properties of the FG model. We get the chemical potential μ by numerical solution of n fe V =    3 shows the comparison of FT KS DFT results with the FG model. The kinetic pressure component for Al is in a good qualitative and quantitative agreement with P F G fe with 3 electrons per atomic volume that manifests week ion scattering of conduction band electrons in Al. We notice that P F G fe for 11 electrons per atomic volume is close to the VASP data on P e in fcc Au by the absolute values despite a somewhat different trend (the strong difference of the similar ABINIT data requires additional consideration of the corresponding PAW potential). An interesting fact is that the sums of kinetic and non-local pressure components for VASP and ABINIT look very www.cpp-journal.org similar. We can clearly see that the non-local part of the PAW potentials plays an important role in the pressure build-up in 2T-WDM metals (it would be instructive to check this effect in orbital-free DFT models [56]).
However the kinetic component of electronic pressure calculated in FT KS DFT includes contributions not only from the delocalized electron states but from all the electrons considered. The correspondence of the absolute values of P F G fe for 3e with PAW of Al=[Ne]2s 2 3p 1 and for 11e with PAW of Au=[Xe4f 14 ]5d 10 6s 1 provides an illustration of this statement. It is possible to notice that the local kinetic energy and the closely related local electronic stress tensor are commonly used to elucidate chemical bonding patterns, especially for covalent bonds [57]. That is why we still need to find a way for the separation of the pressure build-up due to localized and due to delocalized electronic states.
In the case of aluminum Fig. 3 gives an a proiri evident answer: since the P e build-up is equal to the increase of P F G fe with 3 electrons per atomic volume then all 3 electrons occupy the free-like states for T e = 0 − 6 eV. However in the case of gold Fig. 3 gives no evident insight since the trends of P F G fe (T e ) from 1e to 11e per atomic volume essentially differ from the P e (T e ) data.
The results for aluminum show that there is no need to consider any change in the binding pressure P b with the increase of T e . This fact is in agreement with the very week dependence of the Hartree and local components of P e (Fig. 1) and with the vanishing redistribution of electron density (Fig. 2) for fcc Al at T e = 0 − 6 eV. However for fcc Au it is not the case. The ETD EAM model [12][13][14] gives us the data on the P b (T e ) dependence. The ETD EAM potential was obtained using the force matching procedure [58], i.e. its empirical form was optimized to provide the best possible match with the forces acting on ions in the FT KS DFT calculations for a large set of different ionic configurations (for other examples of force matching method applications see [59,60]). The significant difference (P DF T − P EAM ) [13] suggests that the FG component of electron subsystem does not contribute to the ionic forces and the pressure of the ETD EAM potential corresponds only to the binding component of the electron pressure determined by the localized electrons. The best fit of the pressure build-up by the sum P EAM b + P F G fe gives n fe = 3 ± 0.3 electrons per atomic volume for the range T e = 1 − 6 eV (Fig. 4).
The ETD EAM model [12][13][14] is based on three EAM potentials created for T e = 0.8, 3 and 6 eV. Data for other T e values are interpolated. That is why the accuracy of the P b (T e ) dependence at low T e is currently not sufficient to determine the best n fe at T e = 0 − 1 eV by the (P e − P b ) fitting with P F G fe . It seems that n fe = 1 electron could be a more accurate choice than 3 electrons. That could make sense since in this case the free electron number n fe per one gold atom would change form the conventional value of 1 electron at T e = 0 to 3 electrons at T e > 1 eV. Recent experimental data [61] suggest that 5d electrons have prominent contributions to the conduction electrons. These findings agree perfectly with our results if we assume that the electron temperature in [61] reaches T e ∼ 1 eV. In different experiments different electron temperatures can be reached in 2T-WDM state resulting in different n fe values.
In the recent work of Bévillon and coauthors [62] another approach for the calculation of n fe has been used based on the fitting of the FT KS DFT electronic DOS by the square-root FG DOS. This approach gives n fe = 2.4 − 3.5 electrons per ion for fcc Au for T e = 0 − 6 eV. Their approach does not require the knowledge of P b (T e ), i.e. the ETD EAM model for the warm dense metal. However the analysis of electron pressure presented in this work allows us to elucidate the interplay between localized and delocalized contributions to pressure that is important for the applications in hydrodynamic and atomistic simulations. The dependence of the binding pressure on electronic temperature T e manifests a corresponding dependence of the EoS that means the representation (1) provides simplification only in the special cases like Al when one can ignore the dependence of P b and P i on T e . In the general case all three terms in (1) depend on T e and therefore special tools are required for modelling and simulation of 2T-WDM (e.g. the ETD models).

Ionic forces in 2T-WDM
The bond hardening in fcc Au manifests itself by the increase of the phonon frequencies and the increase of P b [13,21]. These facts illustrate certain aspects of the changes with T e that experience the potential energy surface governing ion dynamics of warm dense metal. However, because of the lattice symmetry we have no information about the effects of T e on the ionic forces (they always vanish in a stable crystal lattice). Fig. 5 shows that the forces in 2T-WDM disordered gold increase with T e essentially stronger than in 2T-WDM disordered aluminum (data presented on Fig. 5 correspond to the disordered structures of Al and Au equilibrated via ab initio molecular dynamics at T i = T e = 2000 K). This is another manifestation of the bond hardening in Au.
These forces are the Hellmann-Feynman forces acting on ions (with the contributions from the non-local part of the PAW potentials) that are known to be consistent with the FT KS DFT [63]. However since the ETD EAM model for 2T-WDM Au captures only the P b + P i part of the total pressure, there still stands a question about www.cpp-journal.org the necessity of including the effect of P fe into the ionic forces (i.e. via a blast force [6,11,19]). This question requires further analysis.

Conclusions
The sum of kinetic and non-local pressure components governs the build-up of the electron pressure in FT KS DFT models of warm dense metals. The non-local component of the electron-ion interaction plays an important role in the changes of warm dense metal properties with the increase of electronic temperature.
In the general case the increase of electron temperature in warm dense metal not only causes the free-electron pressure build-up, but also changes the effective interionic forces and hence the binding pressure and the EoS. Therefore in the representation (1) of the total pressure for 2T-WDM as P = P b + P fe + P i all three components depend on T e that should be taken into account in hydrodynamic and atomistic models.
Using effective ion-ion potentials obtained by force-matching for different T e it is possible to separate the binding component of electron pressure and the free electron like component. This is a method for calculation of the number of free electrons, i.e. the ionization degree of WDM.
Additional studies are required for answering the question whether the ion dynamics and mass transfer should be associated with P b + P i pressure and how P fe should be taken into account (i.e. as a blast force).