Predicting the charged particle transport properties of warm dense matter/hot dense plasma mixtures is a challenge for analytical models. High accuracy ab initio methods are more computationally expensive, but can provide critical insight by explicitly simulating mixtures. In this work, we investigate the transport properties and optical response of warm dense carbon–hydrogen mixtures at varying concentrations under either conserved electronic pressure or mass density at a constant temperature. We compare options for mixing the calculated pure species properties to estimate the results of the mixtures. We find that a combination of the Drude model with the Matthiessen's rule works well for DC electron transport and low-frequency optical response. This breaks down at higher frequencies, where a volumetric mix of pure-species AC conductivities works better.

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.4 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.12 

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.11 In addition, recently-determined C–H equations of state for the giant icy planets,26 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 VAB 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 Pe 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)32 and the Stochastic and Hybrid Representation of Electronic Structure by Density functional theory (SHRED)33 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.34 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.

At each time step t for a periodically replicated cell of volume (V) containing Ne active electrons and Ni ions at fixed spatial positions R q ( t ), we first perform a FTDFT calculation within the Kohn–Sham (KS) construction to determine a set of electronic state functions [ Ψ i , k ( t ) | i = 1 , n b] for each k-point k as follows:
H K S Ψ i , k ( t ) = ϵ i , k Ψ i , k ( t ) ,
(1)
where ϵ i , k is the eigenenergy. A velocity-Verlet algorithm advances the ions, based on the force from the ions and electronic density, to obtain a new set of positions and velocities. Repeating these two steps propagates the system in time yielding a trajectory consisting of nt sets of positions and velocities [ R q ( t ) , V q ( t )] of the ions and a collection of state functions [ Ψ i , k ( t )] for the electrons. These trajectories produce a consistent set of static, dynamical, and optical properties. All molecular dynamics (MD) simulations employed only Γ ( k = 0) point sampling of the Brillouin Zone in a cubic cell of length L ( V = L 3 ).

1. Static and transport properties

The total pressure (P) of the system consists of the sum of the electronic pressure Pe, computed via the forces from the DFT calculation, and the ideal gas pressure of the ions at number density n = N i / V,
P = P e + n k B T .
(2)
The electronic pressure is an average over the pressures at different times along the MD trajectory once the system has equilibrated.
Diffusion of warm dense matter mixtures has been examined using effective potential models,35 classical MD and one-component plasma models,36,37 and quantum molecular dynamics and pseudo-ion in jellium models.38–42 For a detailed comparison of different modeling techniques, see the results of the first and second Charged-Particle Transport Coefficient Code Comparison Workshops.12 Following the standard prescription,38 the self-diffusion coefficient Ds is computed from the trajectory by either the mean square displacement (MSD) or by the velocity autocorrelation (VAC) function,
D s = 1 6 t | R i ( t ) R i ( 0 ) | 2 = 1 3 0 V i ( t ) · V i ( 0 )   d t .
(3)
The brackets denote statistical averaging over the trajectories. A similar formula yields the mutual diffusion DAB between the two species.38 Under warm dense matter conditions, the Darken approximation generally provides reliable results using only the self-diffusion coefficients:
D A B = x B D A + x A D B .
(4)
In other words, only interactions of a particle of a given species and itself at different times govern the mutual diffusion. From the e-folding time of the VAC function, we determine a correlation time τ. Time steps separated by τ are considered statistically uncorrelated, and statistical error is estimated from the Zwanzig formula.43 

2. Electrical and thermal properties

The basic electrical and thermal properties of the medium derive from the frequency-dependent Onsager coefficients30,31 given by
L n m ( ω ) = 2 π Ω i , j F i j | D i j | 2 [ ϵ i + ϵ j 2 h ] m + n 2 δ ( ϵ i ϵ j ω ) ,
(5)
where Ω is the atomic volume, and ϵi is the energy of the ith state. We have assumed an implicit summation over k-points. The summed-over quantities are the difference between the Fermi–Dirac distribution at energy levels ϵi and ϵj at temperature T,
F i j = [ f FD ( ϵ i ) f FD ( ϵ j ) ] / ω ,
(6)
and the velocity dipole matrix elements
| D i j | 2 = 1 3 α | Ψ i | α | Ψ j | 2 ,
(7)
with α representing the directions x, y, and z, and ψi is the wave function for the state with energy ϵi given by Eq. (1). For practicality, the δ function in Eq. (5) is approximated by a Gaussian of width Δ G. The enthalpy is defined as h = μ + T s with s, the entropy per particle, and μ, the chemical potential or Fermi energy ϵF. The zero-frequency values of the Onsager coefficients determine basic properties such as the DC conductivity σ dc, the thermal conductivity κ, the thermopower α, and the Lorentz factor L according to the following relations:
σ dc = L 11 ( 0 ) ,
(8)
κ = 1 T [ L 22 ( 0 ) L 12 2 ( 0 ) L 11 ( 0 ) ] ,
(9)
α = L 12 ( 0 ) T L 11 ( 0 ) ,
(10)
L = ( e 2 k B 2 ) κ T σ dc ,
(11)
where e is the electric charge, and k B is the Boltzmann constant. The Onsager coefficients satisfy certain symmetry rules: L n m ( ω ) = L n m ( ω ) and L n m ( ω ) = L m n ( ω ). When the Lorentz factor L is a constant, this relation yields the well-known Wiedemann–Franz law.
The frequency-dependent conductivity L 11 ( ω ) satisfies a simple selection rule of the form
S = 2 π V N e 0 L 11 ( ω ) d ω = 1 ,
(12)
which provides a check on the number of states (bands) employed in the calculation of the optical properties.

Given the behavior of the function Fij 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.

We can extract other optical properties from the frequency-dependent real σ 1 ( ω ) =  L 11 ( ω ) and imaginary σ 2 ( ω ) components of the electrical conductivity. The imaginary part derives directly from a Cauchy principal value ( P) of the integral over the real part,
σ 2 ( ω ) = 2 ω π   P 0 σ 1 ( ν ) ( ν 2 ω 2 )   d ν .
(13)
In terms of the complex conductivity, the components of the dielectric function ϵ ( ω ) = ϵ 1 ( ω ) + i ϵ 2 ( ω ) are written as
ϵ 1 ( ω ) = 1 4 π ω σ 2 ( ω ) ,
(14)
ϵ 2 ( ω ) = 4 π ω σ 1 ( ω ) .
(15)
Furthermore, the real n ( ω ) and imaginary k ( ω ) parts of the index of refraction,
n ( ω ) = 1 2 { | ϵ ( ω ) | + ϵ 1 ( ω ) } ,
(16)
k ( ω ) = 1 2 { | ϵ ( ω ) | ϵ 1 ( ω ) } ,
(17)
combine to give the reflectivity r ( ω ) and the absorption coefficient α ( ω ),
r ( ω ) = [ 1 n ( ω ) ] 2 + k ( ω ) 2 [ 1 + n ( ω ) ] 2 + k ( ω ) 2 ,
(18)
α ( ω ) = 4 π n ( ω ) σ 1 ( ω ) .
(19)
Finally, the Rosseland mean opacity (RMO) κR is given by
1 κ R = 0 B ( ν ) α ( ν )   d ν ,
(20)
where B ( ν ) is the derivative of the normalized Planck function with respect to temperature. Since the function B ( ν ) peaks around 4 k B T, we expect the computed opacities will be most sensitive to differences in the absorption coefficient around this energy.

1. Pressure matching (PM) or Amagat

The pressure matching (PM) scheme44,45 considers a composite sample of two constituents, A and B with particle numbers NA and NB, respectively, at a given total density nij with concentrations xA and xB and temperature T. Varying the densities (volumes Vi) of the pure species until the following conditions,
V A B = V A + V B ,
(21)
P A [ V A ] = P B [ V B ] ,
(22)
are satisfied establishes the PM prescription, where Pi is the pressure of species i at number density ni = Ni/Vi. The other composite properties (ΩAB), such as conductivities, are deduced from the relation
Ω A B = α A Ω A [ n A P ] + α B Ω B [ n B P ] ,
(23)
where α i V i / V A B. The determination of the pressure constraint Eq. (22) necessitates the independent construction of pressure–volume tables for the individual species (A, B) over a requisite range of densities for a particular T. In addition, the single-species properties require calculation at the matched densities, which may vary considerably.

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.

2. Ionization matching

Rather than a pressure balance, another set of possible mixing rules focus on the effective ionization within the mixture, such as the ones proposed in Appendix A of a recent paper by Starrett et al.46 In this case, the mixing rule for a given composite property Ωij becomes
Ω A B = α A Ω A [ ρ ] + α B Ω B [ ρ ] .
(24)
with α i = C i / C , C i = x i ( Z ¯ i p ) 2, and C = i C i, where an effective species charge Z i p within the plasma mixture determines the single-species mixing coefficient, αi. Here, ρ is the same mass density for both the pure and mixed systems, i.e., the isodensity case. Many possibilities exist for the choice of Z ¯ i p, some through average atom formulations.46 
Additional flexibility in the mixing rules comes from adjustment of the composite property, which is mixed. For example, solving
Ω A B 1 = α A Ω A 1 [ n A ] + α B Ω B 1 [ n B ] ,
(25)
rather than Eq. (23) for ΩAB directly. To illustrate, consider the conductivity, σ ( ω ) and its inverse, the resistivity, R ( ω ). Ideally, the inverse of the mixed resistivity would be equal to the mixed conductivity, but, as we will see, this is not the case.
Assuming some effective number of contributing electrons per atom Z ¯, we can define the “valence” electron density of the mixture: n ¯ e = ( x A Z A ¯ + x B Z A ¯ ) / V A B. We can define an effective frequency-dependent scattering rate as γ ( ω ) = ( σ ( ω ) / n ¯ e ) 1. Matthiessen's rule states that the total scattering rate is the sum of all scattering rates and is the basis of the ionization matching procedure.46 We will thus also consider the direct mixing of the effective scattering rates and the resulting conductivity as follows:
γ A B = α A γ A [ n A ] + α B γ B [ n B ] ,
(26)
σ A B = n ¯ e [ x , V A B ] × γ A B 1 .
(27)

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/cm3 and for the isobaric case Pe = 5580 GPa, the value for x H =  x C = 0.5. As an example of the isobaric, a density of 8 g/cm3 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.47 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.48 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.

FIG. 1.

Multi-component QMD-calculated electronic pressures as a function x C for isodensity mixtures of density ρ HC = 10 g/cm3.

FIG. 1.

Multi-component QMD-calculated electronic pressures as a function x C for isodensity mixtures of density ρ HC = 10 g/cm3.

Close modal

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.

FIG. 2.

Multi-component QMD-calculated density as a function x C for isobaric mixtures with pressure equal to Pe = 5580 GPa, which is the value for x H =  x C = 0.5.

FIG. 2.

Multi-component QMD-calculated density as a function x C for isobaric mixtures with pressure equal to Pe = 5580 GPa, which is the value for x H =  x C = 0.5.

Close modal

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.40 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.49 For the isodensity case, the volume expansion required to maintain mass density dominates; the hydrogen diffusion increases as total collisions are diminished.

FIG. 3.

Diffusion as a function x C for (a) isodensity mixtures and (b) isobaric mixtures. The blue and red markers are, respectively, the calculated results for D H and D C. The purple square markers are the D CH results. The black square markers connected by the dashed line are the results given by the Darken relation.

FIG. 3.

Diffusion as a function x C for (a) isodensity mixtures and (b) isobaric mixtures. The blue and red markers are, respectively, the calculated results for D H and D C. The purple square markers are the D CH results. The black square markers connected by the dashed line are the results given by the Darken relation.

Close modal

In both cases, the Darken relation given by Eq. (4) gives strong agreement with the measured mutual diffusion values. In the Darken approximation, DH and DC 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.

FIG. 4.

DC electrical conductivity, σ ( ω = 0 ), for (a) isodensity and (b) isobaric mixtures for different concentrations of carbon χC. From highest to lowest, the lines are (blue) volumetric mix of conductivity, (red) volumetric mix of resistance, (green) Starrett ionization for three different ionization option (see the legend) (purple) volumetric mix of effective scattering rate assuming three different ionization options (see the legend). Black dots are the fully atomistic calculation.

FIG. 4.

DC electrical conductivity, σ ( ω = 0 ), for (a) isodensity and (b) isobaric mixtures for different concentrations of carbon χC. From highest to lowest, the lines are (blue) volumetric mix of conductivity, (red) volumetric mix of resistance, (green) Starrett ionization for three different ionization option (see the legend) (purple) volumetric mix of effective scattering rate assuming three different ionization options (see the legend). Black dots are the fully atomistic calculation.

Close modal
FIG. 5.

DC electrical component of thermal conductivity κ ( ω = 0 ) for (a) isodensity and (b) isobaric mixtures for different concentrations of carbon χC. From highest to lowest, the lines are (blue) volumetric mix of conductivity, (red) volumetric mix of resistance, (green) Starrett ionization for Z C = 4 / Z H = 1 (purple) volumetric mix of effective scattering rate assuming Z C = 4 / Z H = 1. Black dots are the fully atomistic calculation.

FIG. 5.

DC electrical component of thermal conductivity κ ( ω = 0 ) for (a) isodensity and (b) isobaric mixtures for different concentrations of carbon χC. From highest to lowest, the lines are (blue) volumetric mix of conductivity, (red) volumetric mix of resistance, (green) Starrett ionization for Z C = 4 / Z H = 1 (purple) volumetric mix of effective scattering rate assuming Z C = 4 / Z H = 1. Black dots are the fully atomistic calculation.

Close modal
FIG. 6.

Absorbance as a function of photon energy for (a) isodensity mixtures and (b) isobaric mixtures. The different color curves correspond to various mixture ratios shown in the legend. The curves from top to bottom correspond to 0 % C (red), 10 % C (blue), 25 % C (black), 50 % C (green), 75 % C (orange), and 100 % C (purple). The corresponding top to bottom ordering in (b) is taken at photon energy 0.5 on the x-axis.

FIG. 6.

Absorbance as a function of photon energy for (a) isodensity mixtures and (b) isobaric mixtures. The different color curves correspond to various mixture ratios shown in the legend. The curves from top to bottom correspond to 0 % C (red), 10 % C (blue), 25 % C (black), 50 % C (green), 75 % C (orange), and 100 % C (purple). The corresponding top to bottom ordering in (b) is taken at photon energy 0.5 on the x-axis.

Close modal
We compare different options for mixing rules, including the volumetric conductivity mix,
σ C H σ = V H V C H × σ H ( χ C = 0 ) + V C V C H × σ C ( χ C = 1 ) ,
(28)
resistivity mix,
σ C H R = ( V H V C H × σ H 1 ( χ C = 0 ) + V C V C H × σ C 1 ( χ C = 1 ) ) 1 ,
(29)
the mixing rule proposed by Starrett,
σ C H Sta = C H C × σ H ( χ C = 0 ) + C C C × σ C ( χ C = 1 ) , C I = x I ( Z I ¯ ) 2 ; C = I C I
(30)
and, finally, the mixing of the effective scattering rates
σ A B = n ¯ e [ x , V A B ] × γ A B 1 , n ¯ e = ( x A Z ¯ A + x B Z ¯ B ) / V A B , γ A B = V H V C H × γ H ( χ C = 0 ) + V C V C H × γ C ( χ C = 1 ) .
(31)
For the isodensity case, we see that the volumetric mix overestimates the conductivities of the mixtures, while the resistive mix underestimates it. We compare the Starret approach46 for a variety of Z ¯ I options. The Z C = 2.6 / Z H = 0.82 case corresponds to the effective charge from average atom Kohn–Sham calculations from the Tartarus code.50 We also include the valence charges Z C = 4 / Z H = 1 to calculate the Fermi momentum. When using a Thomas Fermi model to calculate Z C / Z H, we note that these results agree better with Drude fits of the AC conductivities,19,51 and the fully ionized charges Z C = 6 / Z H = 1. We see that the fully ionized charges do provide good agreement, but the sensitivity to the charges is large, and full ionization is unreasonable in this temperature density regime, so we expect the agreement is simply fortuitous. In Starrett et al.46 this disagreement between atomistic mixtures and the mixing rule was considered as a consequence of the low temperature. However, mixing the scattering rates, γ mix, gives much less sensitivity to the charges used to define the electron density, with both the average atom and valence charges giving good agreement.

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.52 

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”53 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.

FIG. 7.

Reflectivity as a function of photon energy for (a) isodensity mixtures and (b) isobaric mixtures. The different color curves correspond to various mixture ratios shown in the legend. The curves from top to bottom correspond to 0 % C (red), 10 % C (blue), 25 % C (black), 50 % C (green), 75 % C (orange), and 100 % C (purple).

FIG. 7.

Reflectivity as a function of photon energy for (a) isodensity mixtures and (b) isobaric mixtures. The different color curves correspond to various mixture ratios shown in the legend. The curves from top to bottom correspond to 0 % C (red), 10 % C (blue), 25 % C (black), 50 % C (green), 75 % C (orange), and 100 % C (purple).

Close modal
Electronic transitions at different frequencies are of different nature (e.g., bound–bound, bound–free, and free–free). Thus, we investigate how the efficacy of the mixing rules changes in optical properties between low and high frequencies. In Fig. 8 (Fig. 9), we show the absorbance (reflectance) of the 10% C mixture for the isodensity (top) and isobaric (bottom) cases. Other mixtures are shown in the supplementary material. Generally, we see that the γ mix works well at low frequencies, in both the isodensity and isobaric case, while the σ mix works better for higher frequencies. For the isobaric case, the direct mix of the optical property and the optical property calculated via mixed conductivities gives similar results. The γ mix is based on Matthiessen's rule for adding scattering rates. In this picture, the electron diffuses through the system interacting with different scattering centers. This is appropriate at low frequencies when an electron can diffuse through the system. The insets in the absorbance plots, Fig. 8, expand the low-frequency regime, showing the superior agreement of the γ mix at low frequencies. In a classical picture, at high frequencies, the electron is oscillating rapidly, with limited ability to traverse between scattering centers. Thus, the direct volumetric mix works well for high frequencies. For the reflectance, we can see the agreement more easily. We also plot a transitional mix, where a linear combination of both γ and σ mix is used, weighted by a Fermi–Dirac factor to interpolate between the two,
Ω C H T = F D ( ω ) × Ω γ ( ω ) + ( 1 F D ( ω ) ) Ω σ ( ω ) , F D ( ω ) = ( 1 + e ( ω 1.0 ) / 0.5 ) 1 .
(32)
We have chosen the values of the crossover point and smearing to be 1.0 and 0.5 atomic units, respectively. We can estimate this crossover by considering the following simple model. We approximate that the frequency of a transition given by the change in the kinetic energy of the electron and that the transitions are centered around the thermal electron kinetic energy, ω v e × δ v × m e. Here, v e = 3 k B T + ( k F / m e ) 2 is the thermal electron velocity, and me is the electron mass. We then assume that the crossover occurs when δ v / ( 2 × r W S × m e ), where rWS is the Wigner–Seitz radius of the averaged ion. This leads to a range of crossover frequencies from 1.3   ( 1.7 ) to 0.9   ( 0.8 ) atomic units for the isobaric (isodensity) case when using Z C = 4 / Z H = 1 to calculate the Fermi momentum. When using an average-atom Thomas Fermi model to calculate Z C / Z H, we see a similar range of 0.8   ( 1.3 ) to 0.7   ( 0.6 ) atomic units for the isobaric (isodensity) case. Given the simplicity of this model, we simply use the fixed crossover frequency of 1 a.u. in the plots, and only use this analysis to build a preliminary understanding. Empirically, we see that the transition is broad compared to these differences, and thus the interaction of the electrons and ions is important.
FIG. 8.

Absorbance (α) as a function of photon energy for (top) isodensity mixtures and (bottom) isobaric 10% C mixture. Black line is the atomistic calculation result, green is the volumetrically mixed reflectivity, purple (blue) is the reflectivity calculated from volumetrically mixed effective scattering rate (conductivity). Yellow dotted line is the transitional mix from γ to σ mix at 1 a.u. The inset highlights the low photon energy range from 0 to 0.3 a.u.

FIG. 8.

Absorbance (α) as a function of photon energy for (top) isodensity mixtures and (bottom) isobaric 10% C mixture. Black line is the atomistic calculation result, green is the volumetrically mixed reflectivity, purple (blue) is the reflectivity calculated from volumetrically mixed effective scattering rate (conductivity). Yellow dotted line is the transitional mix from γ to σ mix at 1 a.u. The inset highlights the low photon energy range from 0 to 0.3 a.u.

Close modal
FIG. 9.

Reflectivity (r) as a function of photon energy for (top) isodensity mixtures and (bottom) isobaric 10% C mixture. Black line is the atomistic calculation result, green is the volumetrically mixed reflectivity, purple (blue) is the reflectivity calculated from volumetrically mixed effective scattering rate (conductivity). The yellow dotted line is the transitional mix from γ to σ mix at 1 a.u.

FIG. 9.

Reflectivity (r) as a function of photon energy for (top) isodensity mixtures and (bottom) isobaric 10% C mixture. Black line is the atomistic calculation result, green is the volumetrically mixed reflectivity, purple (blue) is the reflectivity calculated from volumetrically mixed effective scattering rate (conductivity). The yellow dotted line is the transitional mix from γ to σ mix at 1 a.u.

Close modal

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.

1.
S. X.
Hu
,
L. A.
Collins
,
V. N.
Goncharov
,
J. D.
Kress
,
R. L.
McCrory
, and
S.
Skupsky
,
Phys. Plasmas
23
,
042704
(
2016
).
2.
S. X.
Hu
,
L. A.
Collins
,
T. R.
Boehly
,
Y. H.
Ding
,
P. B.
Radha
,
V. N.
Goncharov
,
V. V.
Karasiev
,
G. W.
Collins
,
S. P.
Regan
, and
E. M.
Campbell
,
Phys. Plasmas
25
,
056306
(
2018
).
3.
P.
McKenna
,
D. A.
MacLellan
,
N. M. H.
Butler
,
R. J.
Dance
,
R. J.
Gray
,
A. P. L.
Robinson
,
D.
Neely
, and
M. P.
Desjarlais
,
Plasma Phys. Controlled Fusion
57
,
064001
(
2015
).
4.
V. B.
Prakapenka
,
N.
Holtgrewe
,
S. S.
Lobanov
, and
A. F.
Goncharov
,
Nat. Phys.
17
,
1233
(
2021
).
5.
M. W. C.
Dharma-wardana
,
Phys. Rev. E
73
,
036401
(
2006
).
6.
S. D.
Baalrud
,
Phys. Plasmas
19
,
030701
(
2012
).
7.
S. D.
Baalrud
and
J.
Daligault
,
Phys. Rev. Lett.
110
,
235001
(
2013
).
8.
C. E.
Starrett
,
High Energy Density Phys.
25
,
8
(
2017
).
9.
P. A.
Sterne
,
S. B.
Hansen
,
B. G.
Wilson
, and
W. A.
Isaacs
,
High Energy Density Phys.
3
,
278
(
2007
).
10.
S.
Malko
,
W.
Cayzac
,
V.
Ospina-Bohórquez
,
K.
Bhutwala
,
M.
Bailly-Grandvaux
,
C.
McGuffey
,
R.
Fedosejevs
,
X.
Vaisseau
,
A.
Tauschwitz
,
J. I.
Apiñaniz
et al,
Nat. Commun.
13
,
2893
(
2022
).
11.
S.
Jiang
,
O. L.
Landen
,
H. D.
Whitley
,
S.
Hamel
,
R.
London
,
D. S.
Clark
,
P.
Sterne
,
S. B.
Hansen
,
S. X.
Hu
,
G. W.
Collins
et al,
Commun. Phys.
6
,
98
(
2023
).
12.
P. E.
Grabowski
,
S. B.
Hansen
,
M. S.
Murillo
,
L. G.
Stanton
,
F. R.
Graziani
,
A. B.
Zylstra
,
S. D.
Baalrud
,
P.
Arnault
,
A. D.
Baczewski
,
L. X.
Benedict
et al,
High Energy Density Phys.
37
,
100905
(
2020
).
13.
S.
Mazevet
,
M. P.
Desjarlais
,
L. A.
Collins
,
J. D.
Kress
, and
N. H.
Magee
,
Phys. Rev. E
71
,
016409
(
2005
).
14.
F.
Lambert
,
V.
Recoules
,
A.
Decoster
,
J.
Clérouin
, and
M.
Desjarlais
,
Phys. Plasmas
18
,
056306
(
2011
).
15.
G.
Faussurier
,
C.
Blancard
, and
P.
Cossé
,
Phys. Rev. E
91
,
053102
(
2015
).
16.
Z.
Chen
,
B.
Holst
,
S. E.
Kirkwood
,
V.
Sametoglu
,
M.
Reid
,
Y. Y.
Tsui
,
V.
Recoules
, and
A.
Ng
,
Phys. Rev. Lett.
110
,
135001
(
2013
).
17.
K.-U.
Plagemann
,
P.
Sperling
,
R.
Thiele
,
M. P.
Desjarlais
,
C.
Fortmann
,
T.
Döppner
,
H. J.
Lee
,
S. H.
Glenzer
, and
R.
Redmer
,
New J. Phys.
14
,
055020
(
2012
).
18.
Y.
Cytter
,
E.
Rabani
,
D.
Neuhauser
,
M.
Preising
,
R.
Redmer
, and
R.
Baer
,
Phys. Rev. B
100
,
195101
(
2019
).
19.
M.
Bethkenhagen
,
B. B. L.
Witte
,
M.
Schörner
,
G.
Röpke
,
T.
Döppner
,
D.
Kraus
,
S. H.
Glenzer
,
P. A.
Sterne
, and
R.
Redmer
,
Phys. Rev. Res.
2
,
023260
(
2020
).
20.
D. E.
Hanson
,
L. A.
Collins
,
J. D.
Kress
, and
M. P.
Desjarlais
,
Phys. Plasmas
18
,
082704
(
2011
).
21.
D. A.
Horner
,
F.
Lambert
,
J. D.
Kress
, and
L. A.
Collins
,
Phys. Rev. B
80
,
024305
(
2009
).
22.
C. E.
Starrett
,
High Energy Density Phys.
19
,
58
(
2016
).
23.
N.
Wetta
and
J.-C.
Pain
,
Phys. Rev. E
108
,
015205
(
2023
).
24.
T. J.
Callow
,
E.
Kraisler
, and
A.
Cangi
,
Phys. Rev. Res.
5
,
013049
(
2023
).
25.
F.
Lambert
and
V.
Recoules
,
Phys. Rev. E
86
,
026405
(
2012
).
26.
B.
Cheng
,
S.
Hamel
, and
M.
Bethkenhagen
,
Nat. Commun.
14
,
1104
(
2023
).
27.
I.
Kwon
,
L.
Collins
,
J.
Kress
, and
N.
Troullier
,
Phys. Rev. E
54
,
2844
(
1996
).
28.
L. A.
Collins
,
S. R.
Bickham
,
J. D.
Kress
,
S.
Mazevet
,
T. J.
Lenosky
,
N. J.
Troullier
, and
W.
Windl
,
Phys. Rev. B
63
,
184110
(
2001
).
29.
M. P.
Desjarlais
,
J. D.
Kress
, and
L. A.
Collins
,
Phys. Rev. E
66
,
025401(R)
(
2002
).
30.
V.
Recoules
,
F.
Lambert
,
A.
Decoster
,
B.
Canaud
, and
J.
Clérouin
,
Phys. Rev. Lett.
102
,
075002
(
2009
).
31.
B.
Holst
,
M.
French
, and
R.
Redmer
,
Phys. Rev. B
83
,
235120
(
2011
).
32.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
33.
See https://github.com/lanl/SHRED/ for “
SHRED: Stochastic and hybrid representation electronic structure by density functional theory, a plane-wave DFT code employing Kohn-Sham, orbital-free, stochastic, and mixed stochastic-deterministic DFT methods
.”
34.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
35.
N. R.
Shaffer
,
S. D.
Baalrud
, and
J.
Daligault
,
Phys. Rev. E
95
,
013206
(
2017
).
36.
37.
T.
Haxhimali
,
R. E.
Rudd
,
W. H.
Cabot
, and
F. R.
Graziani
,
Phys. Rev. E
90
,
023104
(
2014
).
38.
P.
Arnault
,
High Energy Density Phys.
9
,
711
(
2013
).
39.
A. J.
White
,
C.
Ticknor
,
E. R.
Meyer
,
J. D.
Kress
, and
L. A.
Collins
,
Phys. Rev. E
100
,
033213
(
2019
).
40.
A. J.
White
,
L. A.
Collins
,
J. D.
Kress
,
C.
Ticknor
,
J.
Clérouin
,
P.
Arnault
, and
N.
Desbiens
,
Phys. Rev. E
95
,
063202
(
2017
).
41.
C.
Ticknor
,
E. R.
Meyer
,
A. J.
White
,
J. D.
Kress
, and
L. A.
Collins
,
Phys. Plasmas
29
,
112703
(
2022
).
42.
J.
Clérouin
,
P.
Arnault
,
B.-J.
Gréa
,
S.
Guisset
,
M.
Vandenboomgaerde
,
A. J.
White
,
L. A.
Collins
,
J. D.
Kress
, and
C.
Ticknor
,
Phys. Rev. E
101
,
033207
(
2020
).
43.
R.
Zwanzig
and
N. K.
Ailawadi
,
Phys. Rev.
182
,
280
(
1969
).
44.
D. A.
Horner
,
J. D.
Kress
, and
L. A.
Collins
,
Phys. Rev. B
81
,
214301
(
2010
).
45.
R. J.
Magyar
,
S.
Root
, and
T. R.
Mattsson
,
J. Phys.: Conf. Ser.
500
,
162004
(
2014
).
46.
C. E.
Starrett
,
N. R.
Shaffer
,
D.
Saumon
,
R.
Perriot
,
T.
Nelson
,
L. A.
Collins
, and
C.
Ticknor
,
High Energy Density Phys.
36
,
100752
(
2020
).
47.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
48.
F.
Jollet
,
M.
Torrent
, and
N.
Holzwarth
,
Comput. Phys. Commun.
185
,
1246
(
2014
).
49.
J.
Clérouin
,
P.
Arnault
,
N.
Desbiens
,
A.
White
,
C.
Ticknor
,
J. D.
Kress
, and
L. A.
Collins
,
Contrib. Plasma Phys.
57
,
512
(
2017
).
50.
C. E.
Starrett
,
N. M.
Gill
,
T.
Sjostrom
, and
C. W.
Greeff
,
Comput. Phys. Commun.
235
,
50
(
2019
).
51.
A. J.
White
,
L. A.
Collins
,
K.
Nichols
, and
S. X.
Hu
,
J. Phys.: Condens. Matter
34
,
174001
(
2022
).
52.
J.
Clérouin
,
V.
Recoules
,
S.
Mazevet
,
P.
Noiret
, and
P.
Renaudin
,
Phys. Rev. B
76
,
064204
(
2007
).
53.
H.
Raether
,
Solid State Excitations by Electrons
(
Springer Berlin Heidelberg
,
Berlin, Heidelberg
,
1965
), pp.
84
157
.

Supplementary Material