Optical and transport properties of plasma mixtures from ab initio molecular dynamics

Understanding the material properties of matter in extreme conditions is a critical task for predicting the behavior of complex high energy density physics experiments, e.g., inertial confinement fusion (ICF),^1–3 as well as astrophysical systems, e.g., the dependence of magnetic fields on planetary composition.⁴ Numerous models exist to predict these material properties, from kinetic plasma models to average atom density functional theory.^5–9 However, in the warm dense matter and hot dense plasma regimes, it remains difficult to develop accurate models.^10,11 Transport properties of electrons and ions in multi-species mixtures are particularly difficult to accurately predict.¹² 

Ab initio calculations, through Kohn–Sham density functional theory, (DFT) of transport from quantum molecular dynamics and the Kubo–Greenwood approach have become a gold standard for matter in extreme conditions.^13–21 However, the computational cost is large. Generating a table with a wide range of densities and temperatures at a variety of concentrations may be prohibitively expensive. Spherically symmetric average-atom models for DFT are much more efficient but can be less accurate and are not directly applicable to mixtures.^22–24 

Mixing rules are used to estimate the properties of a mixture based on knowledge of the pure component systems. They provide a route to calculate transport properties of mixtures from either rapidly generated average atom data or existing or more readily calculated single-species atomistic data. However for optical properties there has been a rather scant testing of the mixing rules. As a testbed, we investigate warm dense carbon (C) hydrogen (H) mixtures at the atomistic level using ab initio (many-atom) DFT and compare it to mixing rule estimations done at the same level of theory. Warm dense CH mixtures are of critical importance to ICF due to the use of high-density carbon or styrene as an ablator material and deuterium/tritium (DT) as fuel.^1,2,25 Thermal conductivities of warm dense C–H and Beryllium have recently been measured at the OMEGA laser facility.¹¹ In addition, recently-determined C–H equations of state for the giant icy planets,²⁶ Uranus and Neptune, may help explain the puzzling differences in their luminosities giving rise to exothermic and endothermic between the similar planetary structures. Here, we consider isobaric, representative of a pressure–temperature equilibrated interface, and isodensity mixtures.

The rest of the article is organized as follows: Sec. II contains the details of theoretical formalism that is used to predict various transport and optical properties of WDM mixtures with a focus on CH mixtures. In Sec. II B, we compare the results of several mixing rules to thermal and electrical conductivity of CH mixtures using pure species values. Section III contains the details of the computational methods, including the quantum molecular dynamics simulations, that are used to calculate properties of WDM mixtures. The results of the applied workflow are shown in Sec. IV. Conclusions and future outlook are given in Sec. V.

We consider a binary mixture of atoms of type A and B at a constant temperature T with a fixed total number of atoms $N A B = N A + N B$ at concentrations $x A = N A / N A B$ and $x B = N B / N A B$ with a volume V_AB and total number density $n A B = N A B / V A B$⁠. Our study examines the trends in static, dynamical, optical, and thermal properties over a range of concentrations from a pure A to a pure B system within two environmental conventions: (1) isodensity with the volume varied to produce the same total mass density for each choice of concentrations $[ x A , x B ]$⁠, including the pure cases $[ x A = 1 , x B = 0 ]$ and $[ x A = 0 , x B = 1 ]$ and (2) isobaric with the volume varied to produce the same electronic pressure P_e for each choice of concentrations $[ x A , x B ]$⁠.

Since the basic formulation and implementation of the molecular dynamics and optical properties simulations appear in a set of earlier papers,^13,27–31 we shall present only a brief overview of the procedures. We have performed quantum molecular dynamics (QMD) simulations employing the Vienna ab initio Simulation Package (VASP)³² and the Stochastic and Hybrid Representation of Electronic Structure by Density functional theory (SHRED)³³ codes within the isokinetic ensemble (constant NVT). The electrons are treated quantum mechanically through plane wave, finite-temperature-density-functional theory (FTDFT) calculations within the generalized gradient approximation (GGA) for the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional having the ion–electron interaction represented by projector augmented wave (PAW) pseudopotentials.³⁴ The nuclei evolve classically according to a combined force provided by the ions and electronic density. The system was assumed to be in local thermodynamical equilibrium (LTE) with equal electron $( T e )$ and ion $( T i )$ temperatures $( T e = T i )$⁠, in which the former was fixed within the FTDFT and the latter kept constant through simple velocity or force rescaling.

Given the behavior of the function F_ij in Eq. (5), the principal contributions to the Onsager coefficients arise from transitions between occupied and unoccupied eigenstates, requiring the determination of a much larger number of states (bands) than for the MD simulation. Fortunately, only between five and ten snapshots along the trajectory are required to converge the optical properties to within a few percent. The separation between sequential snapshots though should exceed the longest correlation time τ determined from the VAC.

This general procedure can be applied considering a constraint of any material property, e.g., the electronic pressure or the electronic chemical potential. If the parameter is similarly sensitive to the plasma environment, then matching will be less sensitive to the choice of the matched quantity. If other factors play a considerable role, such as the nuclear mass as is the case when mass density is constrained (Dalton's law), then the match can be significantly different.

The mixing of optical properties can become even more complex. As all optical properties [Eqs. (13)–(19)] can be computed from the real part of the frequency-dependent conductivity $σ 1 ( ω )$⁠, we may consider whether the derivative property, e.g., the absorption, should be mixed directly or whether the mixed conductivity should be transformed to the derivative property. Since the relationship between the optical properties is not linear, the results will differ.

We focus on a carbon–hydrogen system at T = 10 eV and for the isodensity case $ρ HC$ = 10 g/cm³ and for the isobaric case P_e = 5580 GPa, the value for $x H$ =  $x C$ = 0.5. As an example of the isobaric, a density of 8 g/cm³ recovers this pressure for $x H$ = 0.75 and $x C$ = 0.25 with $N HC$ =128. We employ ten concentration combinations: $x H$ = 1.0, 0.90, 0.80, 0.75, 0.63, 0.50, 0.37, 0.25, 0.10, 0.0, with $x C = 1 − x H$⁠. Since some of these combinations involve small numbers of atoms for a given species, we examine the convergence of various properties as a function of the total number of atoms $N HC$ = 64, 128, 192, 256, and 384. We find for the electronic pressure, diffusion coefficients, the DC electrical conductivity, and the thermal conductivity that $N HC$ = 128 gives values within better than 5% when compared with the $N HC$ = 192 and 256 simulations for most of the concentration combinations. For example, for the $x H$ = 0.75, $x C$ = 0.25 case, the thermal conductivity κ takes the following values: 6600, 7100, 7150, and 7100 (W/m/K) for $N HC$ = 64 (48/16), 128 (96/32), 192 (144/48), and 256 (192/64). An exception is the electron transport properties of pure hydrogen, and very low concentration (1% and 5%) carbon mixtures. We thus perform large 1000-atom simulations in these cases using the SHRED code, which has an efficient orbital and grid parallel implementation of the Kubo–Greenwood approach. This enables us to obtain a highly converged calculation with respect to the system sizes for a particularly sensitive set of CH systems.

For VASP calculations, we apply one- and four-electron “hard” PAW with a cutoff of 700 eV for hydrogen and carbon, respectively. The QMD trajectories consist of 2000–5000 time steps of length 0.1 fs. We generally employ four $k −$ points in the sampling although we have tested with 14, which makes a change of less than 5% in the various properties.⁴⁷ When simulating the dynamic properties, we include states with occupations >10^–5. The calculation of electronic transport properties, via the Kubo–Greenwood approach requires two to three times as many Kohn–Sham states as required to converge the electronic density. For SHRED calculations, the one- and four-electron hydrogen and carbon PAW potentials are utilized.⁴⁸ A long QMD of $∼ 8000$ steps is performed with a time step of 0.02 fs. The optical and transport properties are obtained by calculating the respective values at equally spaced uncorrelated static configurations from the trajectory, and then averaging over multiple configurations. We found that in most cases, averaging over 10 configurations was enough to estimate a converged value, with the exception of pure H system where we doubled the number of configurations in the average to 20.

The electronic pressure for different mixture ratios in the isodensity case are shown in Fig. 1. As the concentration of carbon is increased, a sharp decrease in the pressure is observed. This is largely due to the dramatic change in the volume required to maintain the constant mass density as hydrogen nuclei are replaced with carbon.

For the pressure matched system, the densities as a function of mixture ratio are shown in Fig. 2. The densities are taken from fitting three multi-component QMD calculations at total densities, which are ±10% of a guessed density (taken from a Thomas Fermi model), and then interpolating/extrapolating the results to get the matching density. For these densities, the electronic pressures are all within 2% of the target pressure (5580 GPa), with most cases being within 1/2%. The mass density changes by a factor of $∼ 3.2$ from pure hydrogen to pure C, in contrast to their atomic mass ratio of $∼ 12$⁠. The carbon 1s electrons do not contribute significantly to the pressure in this temperature density regime. Assuming these 1s core carbon electrons are frozen, the change in mass density is readily observed to be roughly equivalent to the change, which would be required to preserve a constant valence electron number density when changing from pure H to pure C, $( ∼ 2.98 = 12 1.007 · 1 4 )$⁠. We note that this agreement holds under these conditions and may not generally be true for other densities and temperatures.

The mass diffusion results are shown in Fig. 3. As expected in both the isobaric and isodensity cases, the carbon self-diffusion is much less sensitive to the change in concentration than the hydrogen self-diffusion.⁴⁰ This is due to the larger mass, which leads to a Brownian-type temperature dominated diffusion in the low concentration regime. The hydrogen self-diffusion is more sensitive to concentration change. Under isobaric conditions in asymmetric mixtures, the hydrogen transport crosses over to a Lorentz gas diffusion, where the hydrogen transport is nearly ballistic in between collisions with the higher charge species.⁴⁹ For the isodensity case, the volume expansion required to maintain mass density dominates; the hydrogen diffusion increases as total collisions are diminished.

In both cases, the Darken relation given by Eq. (4) gives strong agreement with the measured mutual diffusion values. In the Darken approximation, D_H and D_C are the self-diffusion calculated in the mixed system. Thus, the Darken relation should not be considered a “mixing rule.” Rather, it derives explicitly from neglecting inter-species correlations in the full Maxwell–Stefan mutual diffusion, which non-trivially extends to higher numbers of species.^38,39

The direct conduction conductivity $σ 1 ( ω = 0 )$ and thermal conductivity $κ ( ω = 0 )$ are shown from the mixed “atomistic” simulations in Figs. 4 and 5, respectively. The top (bottom) panels are the isodensity (isobaric) results. In both cases, the increase in carbon concentration yields a dramatic reduction in the conductivity. For the isodensity case, the increased volume required to maintain mass density leads to a decreased valence electron density. Thus, the DC conductivity follows similar behavior to the electronic pressure, Fig. 1. For the isobaric case, the valence electron density is nearly conserved (⁠ $2.18 − 2.35 × 10 24 e − / cc$⁠), but there is a drop in conductivity of $∼ 3.5$ from pure hydrogen to pure carbon. This indicates an increased electron scattering due to the carbon ions, which have higher effective charge. Thermal conductivity follows the same behavior; in fact, we observe that the Wiedemann–Franz law works well for these system, with Lorentz numbers only ranging from 2.35 to 2.5  $V 2 / K 2 × 10 − 8$⁠.

For the isobaric case, we see significantly different behavior. All mixing rules, except the direct conductivity mix, provide excellent agreement. This is because, in the isobaric case, the conducting electron density is nearly conserved, a consequence of electronic pressure match, thus the $R −$ and $γ −$ mix are nearly identical. The approximations involved in Starrett's mix also reduce to the γ mix when applied to a constant conducting electron density. The good performance of the σ mix in the isobaric case was also seen previously in gold–aluminum mixtures.⁵² 

Given that the Wiedemann–Franz law works well for these systems, L in Eq. (11) is nearly constant, and the mixtures are all isothermal, the mixing rules all behave similarly when applied directly to the thermal conductivity (κ), replacing σ with κ in all mixing rules. This is shown in Fig. 5.

Optical properties (dielectric function, index of refraction, reflectivity, and absorption) are related to the real frequency-dependent (AC) conductivity through Eqs. (13)–(19). In Figs. 6 and 7, we respectively plot the absorbance and reflectivity of the isodensity (top) and isobaric (bottom) mixes. The isodensity cases are again dominated by the drop in the valence electron density, with increasing concentration of carbon. This leads to a drop in reflectivity and absorbance across frequencies. The plasma frequency being largely dependent on the valence electron density, we see that the “plasma edge”⁵³ of the reflectivity (a measure of the plasma frequency) drops as the concentration of C increases. We also see that the reflectivity begins to drop significantly below the plasma edge as the average scattering rate increases. The transition in the isobaric mix case shows a characteristic increasing of the scattering rate (dampening) in a Drude–Lorentz model as the plasma transitions from hydrogen to carbon under a constant electron density. The absorbance peak shifts to higher frequencies and broadens, while reflectivity drops across frequencies without changing the plasma edge.

We have calculated and analyzed the ion and electron transport properties of warm dense CH mixtures across concentrations from pure hydrogen to pure carbon. We considered two cases, namely, isobaric and isodensity, and tested different mixing rules. We applied these mixing rules to mix atomistic Kohn–Sham DFT results for the pure species and compared it to similar atomistic calculations of the mixtures. Thus, we tested the accuracy of the mixing rules themselves, applied to the best reasonably achievable level of theory in these conditions, without convolution of error from the underlying theory and the mixing rule.

Under the warm dense CH plasma conditions we consider here, a volumetric mixing of the “effective” scattering rate provides good agreement in the isodensity case and excellent agreement in the isobaric case, for electrical conductivity, thermal conductivity, and low-frequency optical properties. This agreement is based largely on the accuracy of the Drude model for these systems and the Matthiesen rule for adding scattering rates. The isobaric mixing is much less sensitive to the mixing model and to choices of ionization required to determine a valence electron density for isolating the effective scattering rate. This is because the isobaric matching procedure naturally produces nearly uniform valence electron density in order to achieve electronic pressure match. Furthermore, we demonstrate that the mixing rule, which best applies to low-frequency excitation, is not the same that applies to higher frequencies. We postulate that the length-scale associated with higher frequency excitations is shorter than the interatomic distance, and thus the volumetric mix of conductivity outperforms the volumetric mix of scattering rates due to the breakdown of the diffusive electron transport picture (Matthiesen's rule).

While we have provided analysis here based on two cases, isobaric (∼5580 GPa) and isodensity (10 g/cc), we believe the general findings of this article will apply to other cases where Drude models reasonably fit to the low-frequency behavior, and atomic atom-free transitions dominate higher frequencies, i.e., partially ionized systems. In the case of colder condensed matter, where bonding plays a significant role, inter-molecular interactions will be more prominent, and mixing rules for optical properties will likely fail.

We have demonstrated that different mixing rules can achieve a wide range of results, and thus one should be careful to consider the accuracy of both the mixed quantities and the rules before applying any procedure. Fully atomistic calculations are thus an invaluable tool to explicitly evaluate the accuracy of these mixing rules.

See the supplementary material for additional plots of reflectivity and absorbance for different concentrations of C ( 1 % ,   10 % ,   25 % ,   50 % ,   75 %, and 90%). These are similar to Figs. 8 and 9 of the main article.

This work was supported by the U.S. Department of Energy through the Los Alamos National Laboratory (LANL). Research presented in this article was supported by Science Campaign 4 and the Laboratory Directed Research and Development of LANL under Project Nos. 20210233ER and 20230322ER. We gratefully acknowledge the support of the Center for Nonlinear Studies (CNLS). This research used computing resources provided by the LANL Institutional Computing and Advanced Scientific Computing programs. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

The authors have no conflicts to disclose.

Alexander J. White: Conceptualization (equal); Data curation (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Galen T. Craven: Data curation (lead); Investigation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Vidushi Sharma: Data curation (supporting); Investigation (equal); Writing – original draft (supporting); Writing – review & editing (equal). Lee A. Collins: Conceptualization (equal); Data curation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.