Observations of non-linear plasmon damping in dense plasmas

The physics of dense plasmas is dictated by strong ion-ion correlations and partially degenerate electrons so that quantum statistical methods have to be applied in order to derive the equilibrium and non-equilibrium properties accurately. In particular, correlation and quantum effects become even more pronounced in the warm dense matter region (WDM), i.e., for matter at solid-state-like density and beyond with temperatures of a few to tens of eV.¹ WDM is relevant for astrophysical objects like giant planets^2–4 which have been observed in great number over the past two decades. To understand their striking diversity is a great challenge to current formation, evolution, interior, and dynamo models⁵ and perhaps inaccessible without a better knowledge of the properties of WDM. Furthermore, WDM is also important for another grand challenge of plasma physics, inertial confinement fusion,⁶ since cryogenic deuterium-tritium capsules have to be driven to ignition by powerful lasers like the National Ignition Facility (NIF) through the WDM state.

Besides equation-of-state data, the optical and transport properties are key for that purpose. Since the dielectric function or, equivalently, the dynamic structure factor (DSF) is closely related to the differential scattering cross section of x-rays focused onto WDM or a dense plasma, experiments using brilliant x-ray sources, generated either by powerful optical lasers⁷ or using free electron lasers (FELs),⁸ are perfectly suited to determine quantities like reflectivity and the absorption coefficient or the dynamic electrical conductivity. This has been shown impressively in recent highly resolved x-ray Thomson scattering (XRTS) experiments.^9,10

In particular, the electrical conductivity has been used for plasma diagnostics for decades. While the analytical Spitzer¹¹ and Ziman theory¹² are valid for low-density plasmas or liquid metals, generalized methods have been formulated for a wider range of densities and temperatures like, e.g., the Lee-More model^13,14 or correlation function approaches.^15,16 Nevertheless, the performance of these methods is limited due to the use of the relaxation time approximation or perturbation theory.¹⁵ 

Further, the Drude model is often used to describe the dynamic conductivity of simple metals. However, earlier absorption measurements¹⁷ suggested that even in aluminum with three quasi-free electrons per atom, the Drude model breaks down below the ionization energy. It has been shown recently¹⁰ that the understanding of XRTS spectra of dense plasmas requires a microscopic description beyond the conventional Drude model.

Ab initio methods like molecular dynamics simulations based on finite-temperature density functional theory (DFT-MD) offer a promising alternative. Although computationally expensive, this method avoids the use of plasma parameters like the ionization degree which is not well-defined for WDM states. Furthermore, this method does not rely on effective pair potentials or two-particle cross sections. Therefore, DFT-MD simulations are much better suited to describe the thermophysical properties of WDM.^18,19 As an alternative to the full Kohn-Sham schema, the density-functional neutral-pseudoatom (NPA) model has been proposed for dense plasmas and WDM, see Refs. 20–22.

DFT-MD is now widely used to compute the material properties needed for modeling astrophysical objects⁴ and for the design of inertial confinement fusion (ICF) targets²³ 

However, the accuracy of a DFT-MD calculation depends on the choice of the exchange-correlation (XC) functional. Additionally, thermal exchange and correlation effects have only recently been integrated into local density approximation (LDA)- and generalized gradient approximation (GGA)-type functionals, and the quality for the WDM regime is still an open question.²⁴ Comparison of DFT-MD based optical properties with newly obtained experimental results is of great interest in this context.

Following this introduction, the theory and details of the DFT-MD simulations are described in detail. In Sec. III, we discuss the dielectric function as derived from our simulations. In Sec. IV, we compare the frequency-dependent dielectric function to absorption and, via the fluctuation-dissipation theorem, to scattering measurements and, in particular, discuss the low-frequency limit.

We perform DFT-MD simulations using the Vienna ab initio simulation package (VASP),^25–27 which calculates the electronic structure within DFT. The forces onto the classical ions are derived in each time step from this electronic structure via the Hellmann-Feynman theorem within the Born-Oppenheimer approximation.

The central approximation in DFT is the choice of the XC functional for which various approximations exist. They are suited for different types of applications²⁸ and can be ordered in a hierarchy ladder after Perdew et al.^28,29 The first ladder rung is the local density approximation, followed by the generalized gradient approximations [e.g., PW91,³⁰ AM05,³¹ Perdew, Burke, Ernzerhof approximation (PBE),³² and other functionals], to meta-GGA (SCAN³³ and others), hybrid functionals [Heyd, Scuseria, Ernzerhof (HSE)^34,35 and others], and fully nonlocal functionals.²⁸ In several other papers, electrical conductivities for warm dense matter were calculated with the PBE^36–42 functional, without studying the influence of the XC functional in further detail. There are only few examples where the HSE functional was used to improve the theoretical predictions for the electrical conductivity compared with calculations based on PBE.^3,10,43–45

In this work, we perform comparative electronic structure calculations with three different XC functionals: first, the “Perdew, Burke, Ernzerhof approximation” (PBE),³² second, the “Strongly Constrained and Appropriately Normed Semilocal Density Functional” (SCAN),³³ and, third, the hybrid functional from “Heyd, Scuseria, Ernzerhof” (HSE).^34,35 HSE contains parts of the non-local Fock exchange energy and is explicitly dependent on occupied states. The latter calculations are approximately two orders (one order) of magnitude more expensive than PBE (SCAN) and therefore only feasible for static DFT calculations yet.

We apply a GW-labeled 11-electron projector-augmented wave (PAW) pseudopotential^46,47 for aluminum with a plane-wave cutoff energy of 500 eV. Our DFT-MD simulations were performed at several temperatures using a Nosé thermostat.⁴⁸ The timesteps in the MD simulations had a duration of 0.8 fs to 1.5 fs, and the simulations were run in total for 10 ps to 15 ps after equilibration. For the evaluation of the electronic transport coefficients, we selected 10 to 20 uncorrelated ionic configurations generated with the DFT-MD simulation runs. Importantly, all DFT-MD simulations were run with the PBE functional. The subsequent calculations of the electronic transport coefficients with the HSE, SCAN, and PBE functionals were then made with these ion distributions. We checked convergence of our results against the number of particles and the k-point sets used for the evaluation of the Brillouin zone integrations. At temperatures above 0.3 eV our simulations yielded well-converged results for Baldereschi mean-value point and 64 atoms, which allowed us to reproduce results from earlier work.⁴¹ Particle numbers were checked from 16 to 108 and a number of 4500 bands were needed for 64 particles to achieve convergence for high frequencies of the dielectric function. For smaller temperatures in the liquid phase, a denser k-points grid was necessary, i.e., 3 × 3 × 3 or 4 × 4 × 4 Monkhorst-Pack grids.⁴⁹ 

For the calculation of transport properties, we employ the following Onsager coefficients from linear response theory:^40,50–52

The transition matrix elements of the velocity operator $v̂$ with the Bloch states $|kμ〉$ with the corresponding Fermi occupation number $fkμ$ and eigenenergy $Ekμ$ are calculated from the optical routines of VASP.^44,53 Additionally, ω, Ω, e, and h_e describe the frequency, simulation cell volume, elementary charge, and enthalpy per electron, respectively. The real parts of the dynamic electrical and thermal conductivities, σ(ω) and λ(ω), respectively, are then related to the Onsager coefficients as follows:

The imaginary part of the electrical conductivity $Im[σ(ω)]$ has to be determined from a Kramers-Kronig relation.⁵⁴ Consequently, we have direct access to the dielectric function ϵ(ω) using

In the following paragraph, the dielectric function will be connected to experimentally accessible quantities like the absorption coefficient or the dynamic structure factor for XRTS measurements.

We calculate the frequency-dependent absorption coefficient α(ω) via

For the description of XRTS measurements, the key quantity is the dynamic structure factor,⁵⁵ which can be calculated in reciprocal space k from an improved Chihara⁵⁶ approach

The first term is the ion feature, which describes scattering on the tightly bound and screening electrons with the total form factor N(k). These electrons follow the ionic motion described by the dynamic ion structure factor S_ii(k, ω). The latter reveals access to sound speeds of acoustic modes and other properties of the ionic system.^57,58 In our modified Chihara approach, the electronic transition feature S_et(k, ω) contains the physics of free-free (i.e., plasmons), bound-free, and bound-bound transitions via the fluctuation-dissipation theorem using the complex dielectric function

as successfully applied in Ref. 10. The advantages of this approach are, first, that the total electron structure factor S_ee(k, ω) can be consistently calculated from DFT-MD simulations. In particular, we avoid using an artificial composition of different theoretical approaches by consistently calculating all transition matrix elements in the long wavelength limit, which can be interpreted as free-free, bound-free, and bound-bound transitions; i.e., best-fit procedures for plasmon and bound-free features are unnecessary.^59,60 The plasmon dispersion for finite k vectors can be calculated by adjusting the Mermin dielectric function to the DFT dielectric function,¹⁰ see Eq. (3). Second, no distinction between bound and ionized electrons for the total form factor N(k) is needed because the wave functions of all electrons are calculated with DFT.^58,61 Third, the dynamics of the ionic system is included, which can be measured by using a narrow bandwidth instrumental function (i.e., small probe laser width and small detector broadening) in XRTS.

Figure 1 shows the dynamic electrical conductivity for the HSE and PBE XC functionals. We observe a Drude-like conductivity for energies $ℏω<21 eV$⁠. Therefore, we apply a Drude fit

to our results for the real part of the electrical conductivity. This determines the zero-frequency value for σ(ω).^37,41,64 In the same way, we determine the zero-frequency value of the dynamic thermal conductivity λ(ω),⁴¹ which will be discussed in Sec. IV C. Between 21 eV and the 2p ionization energy, we observe the expected non-Drude behavior⁶⁵ of aluminum. The local minimum around $ℏω=20–25 eV$ is known as the Cooper minimum^66,67 for transition matrix elements. At T = 0.5 eV, the minimum is located at 22 eV (24 eV, 25 eV) for the PBE (SCAN, HSE) XC functional, respectively. Above the ionization energy of 2p electrons ΔE_2p–free, defined as the energy difference of the mean 2p state eigenenergy to the Fermi energy, we observe a sharp increase in the conductivity. We find different ionization energies for the HSE, SCAN, and PBE functionals as the 2p eigenenergies are shifted downward, respectively. By distinguishing between the different states in the transition matrix [see Eq. (1)] contributing to the dynamic electrical conductivity, we also show the portion of 2s ionization which occurs for $ℏω>ΔE2s−free$ and makes up roughly 10% of the bound-free transitions.

For cold solid aluminum, the ionization energies are given in the literature¹⁷ by ΔE_2p–free = 72.5 eV, ΔE_2s–free = 117.8 eV, and (not studied in this work) ΔE_1s–free = 1559.6 eV. The transition energies discussed in this work are shown in Table I. In Ref. 10, the HSE ionization energy was found to give good agreement for warm dense aluminum, in contrast to those from PBE. The origin of this difference is the so-called bandgap problem, related to the lack of a derivative discontinuity in PBE.^30,68 Nevertheless, we find for the conditions of this work that the ionization energy for the 2p (2s) state from HSE is smaller by 3.0 eV (7.9 eV) than given in Ref. 17 for the cold solid.

In this work, we also compare with the dynamic electrical conductivity from an Average Atom (AA) model calculated with PBE.^69,70 AA calculations are computationally less expensive than DFT simulations, which is an advantage when large data sets of spectra are needed for astrophysical or ICF simulations. However, critical tests against the more accurate DFT method are needed for validation. The AA model also shows the Cooper minimum in warm dense aluminum. We find a good overall agreement with the DFT simulations using the PBE XC functional. However, deviations occur at small frequencies, but these tend to vanish at higher temperatures (above 6 eV). Using an AA model which is capable of predicting accurate band energies for the calculation of the dielectric function in a less correlated, the warmer regime will therefore be very useful as DFT-MD calculations become more demanding with increasing temperature.

We study the influence of temperature on the dynamic electrical conductivity in Fig. 2. With increasing temperature from T = 0.5 eV to T = 6 eV, the Cooper minimum becomes less pronounced and disappears at T = 12 eV. We also observe a pre-edge occurring near $ℏω=60 eV$⁠, as observed in XANES spectra¹⁸ or in Ref. 71, which arises from unoccupied s-type states in the conduction band below the Fermi energy (see Sec. III C). Our calculation at T = 12 eV ≈ 0.2 ΔE_2p–free describes a regime, where the outermost 2p electrons are partially ionized (approximately 1%) due to the high temperature entering the Fermi-Dirac distribution function. This allows for 2s–2p bound-bound transitions, which produce a signal in the dynamic conductivity at 40.4 eV, 40.1 eV, and 38.9 eV in HSE, SCAN, and PBE, respectively. In similar spirit, Desjarlais et al.³⁶ found 3s–3p bound-bound transitions in hot low-density aluminum vapor, where continuum states localize at atoms.⁷² 

The ionization energies in the warm dense matter regime need further validation from experiments since models for ionization potential depression (IPD) yield strongly different results.⁷³ In even more compressed, hotter and more highly ionized aluminum, the 2p state will be even less populated, which means the 2s–free transitions will contribute to the total bound-free spectrum in a larger portion. As a consequence, the 2s excitation energy will possibly become measurable by XRTS in analogy to Ref. 10. In general, these intra-atomic transitions contribute to the understanding of IPD and the population of electronic states in the warm dense matter regime.^10,73,74

We project the plane wave functions of the electrons in aluminum onto s–type, p–type, and d–type wave functions to understand the characteristics of electronic transitions contributing to the electrical conductivity and discuss our results using the analogy to well-known selection rules, see Fig. 3.

The dc value in the Onsager coefficient L₁₁(ω) is determined by transitions between partially occupied states with ΔE → 0 around the Fermi edge, see Eq. (1). In the case of Fig. 3, the relevant states are located at energies of E_F ± 0.5 eV. The probability of transitions is influenced by the shape of the orbitals and can be interpreted with selection rules for specific transition matrix elements.⁶⁷ For instance, a transition from a s–type (blue) state to a p–type state (green) is more likely than a transition from the s–type state to another s–type state as the transition matrix element goes to zero. As we showed above, the HSE XC functional does not shift all electronic states by the same energy, i.e., the 2s states more than the 2p states. Consequently, the shifts in energy eigenvalues due to the XC functional (compare DOS given in Fig. 1 of Ref. 10) depend on the type of the electronic state, and this also changes the position of the Cooper minimum, the DOS, and the dc conductivity. In particular, HSE shifts the s–type and p–type DOS in the conduction band relatively to each other, leaving less s–type states at the Fermi energy and decreasing the dc value, consistent with the prediction of Karasiev et al.²⁴ 

In this section, we converted the dynamic conductivity to the absorption coefficient by using Eqs. (3)–(5). Also the Mermin approach^10,75 was used to derive the structure factor with Eq. (6) for XRTS. We compare to experimental results and discuss the zero-frequency limit of the thermal and electrical conductivity.

In Fig. 4, we show the absorption coefficient, see Eq. (4), from DFT using the HSE functional. The dielectric function used at T = 0.3 eV has been benchmarked against XRTS spectra measured at the Linac Coherent Light Source (LCLS), where good agreement was found for the free-free and bound-free transitions.¹⁰ Here we compare with different measurements^65,76,77 in cold aluminum and find reasonable agreement. We mention that the measurements and our calculations differ from the DFT calculation of Vinko et al.,³⁸ who found the Cooper minimum at ∼30 eV using a three electron potential with PBE. The agreement of the DFT-MD dielectric function, including the Cooper minimum, with XRTS measurements¹⁰ and absorption measurements shows the importance of a correct quantum description of transition matrix elements in the dielectric function. In particular, constant collision frequencies [in the Random Phase Approximation (RPA), or in the Born-Mermin Approximation (BMA)], which do not capture the microscopic physics of the free-free transitions, are in our cases not suited to describe absorption³⁸ and XRTS measurements¹⁰ of warm dense aluminum correctly.

The experimental campaign at LCLS^9,10 measured collective XRTS from 50 μm thick aluminum foils using 7.98 keV photons. Focal spot sizes of 1 μm and 10 μm of the 0.1 mJ x-ray pulse heated the system isochorically to T_e = 0.3 eV and T_e = 6 eV. The system was probed under different scattering angles.¹⁰ In Fig. 5, we show a XRTS spectrum at 24° for a 10 μm spot size. We calculate the plasmon signal at T_e = 6 eV as determined in Ref. 9 from the dynamic electrical conductivity shown in Fig. 2 and find good agreement with the experiment. In this regime, the non-Drude behavior of the conductivity (also present in AA calculations⁷⁰) is still causing strong frequency-dependent damping not included in the BMA. This, again, verifies our calculation of the dielectric function. Furthermore, from analyzing parts of the absorption data from Henke et al.,¹⁷ we propose that plasmons in zinc (Z = 30), gallium (Z = 31), molybdenum (Z = 42), indium (Z = 49), tantalum (Z = 73), gold (Z = 79), and other materials might also be non-linearly damped under warm dense matter conditions, where analytical models for the collision frequency break down.

The ion temperature for this XRTS spectrum is inferred¹⁰ from a DFT-MD fit procedure.⁷⁸ At solid density, we run several simulations at different ion temperatures and compare the ionic structure factor with the experimental value. We found best agreement by using 930 K for all data at different wave numbers.¹⁰ This is below the melting temperature, and the ionic structure factors for these cases show Laue diffraction peaks with temperature dependent diffusive background.

In Fig. 6, we show the plasmon width for different wavenumbers from the experiment⁹ at 1 μm focal spot size of the LCLS (T = 6 eV) and compare with theoretical calculations. Good agreement is found with our DFT-MD collision frequency in the Mermin approach. The BMA and RPA do not reproduce the measured plasmon dispersion.

In 2015, Sperling et al.⁹ showed, for the first time, that the dc conductivity value can be extracted from the frequency dependent scattering signal. They used a generalized Drude expression to extrapolate the dielectric function below $Δℏω≈5 eV$ to $Δℏω=0$⁠. Due to the presence of the ion feature around $Δℏω≈0 eV$ in XRTS, an extrapolation to the dc value is generally necessary and dependent on the extrapolation model used. With a narrower instrument function, the quality of the extrapolation will improve, resulting in a decrease in the error bars. In this work and in Ref. 10, we have shown that XRTS spectra of warm dense aluminum are described well by DFT-MD simulations all across the accessible frequency range for different temperatures and different scattering angles. Best agreement was found by applying the HSE XC functional. Therefore, $σdcHSE$ from our DFT-MD simulations (shown in the next paragraph) represents an improved extraction of dc conductivities from XRTS compared to the evaluation of the generalized Drude-expression as done in Ref. 9.

As pointed out above, zero-frequency values of the dynamic electrical and thermal conductivity in the WDM regime are of great interest for astrophysical applications or fusion science. In addition, also the behavior of the Lorenz number^79,80

is not fully understood in the WDM regime yet. For highly degenerate electrons, the Wiedemann-Franz law^79,80 holds and states L = π²/3, in contrast to the Spitzer limits for fully ionized plasmas at high temperatures.^11,16 Following the summary of Konôpková et al.,⁸¹ the Lorenz number varies strongly around the value given by the Wiedemann-Franz law: large deviations are reported in the literature for different temperature and density conditions in iron (±30%), platinum (±30%), molybdenum (–30%), copper (–30%), or in tungsten (±31%).

In the upper panel of Fig. 7, we show the thermal and electrical conductivity from the HSE and PBE functionals as a function of temperature. In the lower panel, we show the Lorenz number from Eq. (8). Our PBE results in agreement with the Bloch-Grüneisen fit of Ref. 41, which was obtained from DFT-MD results using the PBE functional at temperatures from 0.1 eV to 0.8 eV.⁴¹ For temperatures from 0.3 eV to 1.0 eV, the HSE functional predicts Lorenz numbers higher by up to 14% than those from PBE.

In addition, both our DC conductivity values and those from Ref. 41 below T = 0.2 eV and above T = 1 eV disagree with corresponding results calculated with an orbital-free DFT implementation combined with electronic transport calculations using MD simulations,⁸⁴ which might be caused by the applied pseudopotential.

Our calculations of the thermal conductivity agree, again, with the Bloch-Grüneisen fit of Vlček et al. at small temperatures. We find the HSE functional to yield smaller thermal conductivities (by up to 15%), and we observe a monotonous increase in the thermal conductivity with temperature. Our calculations agree qualitatively with experimental data given by Mills et al.⁸³ for the expanded metal. Recently, McKelvey et al.⁸² published thermal conductivities from proton-heated warm dense aluminum. The data are accurate within one order of magnitude in the error bars at densities of 1.7 g/cm³ to 2.7 g/cm³ as shown also in Fig. 7. For the experimental data,⁸² also an EOS had to be constrained from hydrodynamic simulations. Our calculations agree with rather uncertain published data.

We compare our dc conductivity values in Fig. 8 to experiments and theories from the literature. We observe the following trend that the XC functionals influence the dc conductivity in the following way: $σdcPBE>σdcSCAN>σdcHSE$⁠. Similar to the result found in warm dense helium,⁴⁵ these differences decrease with increasing temperature. In our case of solid-density aluminum, they vanish above a temperature of T ≈ 2 eV, where the differences in the DOS of the continuum become rather negligible. It is nevertheless possible that differences in the conductivity reappear at higher temperatures where higher ionization states of aluminum are present. Note that the conductivity values reported here at T = 0.3 eV for the HSE functional are slightly higher than those published earlier in Ref. 10, where a part of the non-local contributions in the transition matrix elements was not taken into account.

In the temperature range below T = 0.5 eV, we compare with the precise liquid metal data of Gathers.^86,89 Two data sets are given by Gathers: First the direct measurement, density corrected, which means it follows an isobaric expansion (Fig. 9), and second, density uncorrected data, which is a simple extrapolation to solid density in Fig. 8.⁸⁶ For temperatures above T = 2 eV, we find good agreement with the measurement of Milchberg et al.⁸⁵ (blue squares). Furthermore, our DFT-MD results lie in between two analytical average atom based calculations.^87,88

In Fig. 9, we compare our calculations for the densities along the isobaric expansion with PBE to the calculations of Vlček et al.⁴¹ and find agreement. However, these values overestimate the direct conductivity measurements of Gathers considerably (red filled squares). The direct measurements (isobaric) of Gathers are labeled as density corrected, so are isobaric measurements of Desai et al.,⁹⁰ whereas the isochoric values are labeled as density uncorrected. When we apply the SCAN or HSE functional for those conditions, we find reduced conductivities and better agreement with the direct isobaric measurements which again demonstrates the high predictive power of the DFT-MD method but also the influence of the XC functional.

In conclusion, we have calculated the dynamic electrical and the thermal conductivity of warm dense aluminum at solid density from DFT-MD simulations for various temperatures. We find a strong influence of the exchange-correlation functional on the electronic transport coefficients by comparing PBE and HSE with the recently developed SCAN functional. By benchmarking the calculated dielectric functions against two different types of experiments, absorption and scattering measurements, we gain confidence that HSE is best suited to describe warm dense aluminum plasmas. Additionally, we calculate the dc conductivity and the thermal conductivity values. We find again that HSE yields the best agreement with experimental measurements. These noticeable differences due to the XC functional can be attributed to shifts of electronic states inside the continuum. Details on characteristic deviations from the Drude-like conductivity are discussed as well. We show the direct relation of absorption and XRTS spectra and thus propose non-linear plasmon damping for other materials. Finally, our DFT calculations confirm that the Lorenz number should not be assumed constant in the warm dense matter regime as suggested by the Wiedemann-Franz law, which is valid only in highly degenerate plasmas.

We thank A. Baczewski, S. B. Hansen, M. P. Desjarlais, D. Saumon, and K. Burke for helpful discussions. We thank C. Starrett for providing AA calculations to us and for helpful discussions. We acknowledge frequent and helpful discussions with C. Dharma-wardana concerning the interpretation of the dc conductivity values given by Gathers.^86,89 The ab initio calculations were performed at the North-German Supercomputing Alliance (HLRN) facilities. This work was supported by the DFG within the SFB 652 and FOR 2440, and by the DOE Office of Science, Fusion Energy Science under FWP 100182.