Electrical conductivity of hydrogen plasmas: Low-density benchmarks and virial expansion including e–e collisions Open Access

Equilibrium correlation functions can be evaluated using the methods of quantum statistics, such as the Green's function method.^4,5 This is a perturbative approach in which an expansion with respect to the Coulomb interaction is performed.^6–8 By applying partial summations describing, e.g., quasiparticles, screening, and the formation of bound states, useful results are obtained for a wide range of plasma parameters. The application to electrical conductivity was also shown in Ref. 9. For the relation to kinetic theory, see Ref. 3. However, exact results are only obtained in special limiting cases. In the following, we discuss the virial expansion of the inverse conductivity of the hydrogen plasma, which contains these exact results. They can be employed as benchmarks for any approach.

A possibility to avoid the perturbation-theoretical treatment of the electron–ion interaction is provided by the density functional theory (DFT), which is already successfully used in condensed matter physics. The application to the calculation of electrical conductivity was made within the framework of the DFT-MD (molecular dynamic) simulations.¹⁰ The ion dynamics are treated classically using molecular dynamics (MD) simulations, while the electrons are treated as quantum particles, with the wave equations being solved using the DFT formalism. The evaluation of the correlation function (4) is performed numerically using the Kubo–Greenwood formula. However, the treatment of the electron–electron (e–e) interaction is approximated by a suitable choice of the exchange-correlation functional E_xc.¹¹ This problem is also a subject of the present work.

The rigorous treatment of the e–e interaction is possible with the help of path-integral Monte Carlo (PIMC) simulations (see Refs. 12–14 and references therein). Note that the equilibrium correlation function in Eq. (4) is determined by the corresponding spectral function, which can be calculated using the imaginary time (Matsubara) technique.^4,5 However, the calculations are very complex, and the current shortcomings of this approach include the relatively small number of particles and the sign problem for fermions.

In this work, we focus on two aspects: the virial expansion of the resistivity of hydrogen plasmas in the low-density region and how DFT-MD simulations can be improved to describe the contribution of electron–electron collisions to plasma conductivity. In Sec. II, we discuss the virial expansion of the resistivity. The DFT-MD simulations of plasma conductivity are presented in Sec. III, and ways to account for e–e collisions are discussed in Sec. IV. Finally, we draw conclusions in Sec. V.

For the sake of simplicity, we consider hydrogen plasmas. The results can be transferred to other plasmas with ionic charge Ze, since the low-density limit considered here is dominated by the Coulomb force and short-range interactions contribute to higher orders of the virial expansion of the resistivity.

This method was applied to several approaches for determining the electrical conductivity of hydrogen plasmas, in particular to the data presented in Ref. 20. As shown in Ref. 19, several semi-empirical expressions for plasma conductivity fail to fulfill this benchmark.

The choice of the reference density influences the higher virial coefficients. In principle, $n ̂ ref ( T )$ can be chosen such that the next virial coefficient disappears, $ρ ̂ 2 ( T ) = 0$⁠, as shown at the end of this section [see Eq. (28)]. The logarithmic term in the virial expansion gives us the freedom to choose $n ref ( T )$ in such a way that the contribution of the higher order virial terms is minimized.

As example, the values for the conductivity for specific values of T, n are shown in Table I. The values for low densities and high T are considered where the parameter $x , x ̃ = 1 / Λ$ or $x ̂ = 1 / Λ ̂$ is small so that the virial expansion with error bars below few percent is valid. The largest discrepancy occurs for the relatively high density 10 g/cm³ and the low temperature 100 eV. The method of virial plots can also be used to extract the virial coefficients $ρ 1 ( T ) , ρ 2 ( T ) , ρ ̃ 2 ( T )$⁠, or $ρ ̂ 2 ( T )$ from experimental data or numerical simulations. An example is discussed in Sec. III, see also Ref. 19.

The method of virial plots has also been applied to thermodynamic variables describing the equation of state of plasmas (see Ref. 23). There, higher order virial coefficients and generalized virial expansions are discussed. Higher order virial coefficients are of interest to describe degeneracy as known from the Lee–More model (see, e.g., Ref. 22).

Due to the finite size of the simulation box, the δ-function in Eq. (29) must be approximated by a Gaussian with finite width, which also prevents the direct calculation of the dc conductivity at ω = 0. Therefore, the dynamic conductivity is extrapolated to the limit $ω → 0$ using a Drude fit. The extrapolation procedure to ω = 0 can be improved by using a frequency-dependent collision frequency.⁸ 

One of the main shortcomings of the DFT-MD approach is that the many-particle interaction is replaced by a mean-field potential. When using product wave functions for the many-electron system, correlations are excluded. The exchange-correlation energy density functional reflects the Coulomb interaction in a certain approximation, e.g., as it exists in the homogeneous electron gas. However, it becomes problematic in the low-density limit, where correlations, especially e–e collisions, are important.

DFT-MD simulations have been successfully used to calculate the transport properties in the warm–dense matter region for various materials. The transition to condensed matter, the liquid or solid state, is adequately described, in particular for the degenerate electron system. The question whether the contribution of e–e collisions to the conductivity is correctly described also in the low-density range has long been the subject of controversial debate. The e–e Coulomb interaction is included in the exchange-correlation energy of the DFT functional so that it was argued that the ab initio approach also describes the Coulomb interaction of the electrons.^10,30 However, it was pointed out⁸ that DFT-MD is not able to reproduce the Spitzer limiting value of plasma conductivity, what leads to significant deviations in the low-density range. The correct behavior of the plasma conductivity at low density was determined using generalized linear response theory (gLRT).

We would like to mention that in the case of thermal conductivity, it has been shown that the contribution of e–e collisions is not taken into account in DFT-MD simulations^10,22,30 and results in an additional term. Other approaches such as generalized linear response theory can be used to solve this problem.

The DFT-MD simulations have been used very successfully to calculate the transport properties of warm and dense matter (WDM) in the region of degeneracy. However, electron–electron collisions are not treated adequately. For degenerate systems, e–e collisions are suppressed by Pauli blocking so that their contribution can be neglected and the results of DFT-MD are applicable. The systematic inclusion of e–e collisions is indispensable for the description of plasmas at $Θ > 1$⁠.

We refer to a Lorentz plasma as a model system of electrons and ions, in which the interaction of electrons with ions is considered [see Eq. (32)]. The electron–electron interaction is only considered in the mean-field approximation so that the electron–ion interaction is screened, but the effect of electron–electron collisions on the transport coefficients is neglected. For this Lorentz model plasma, the solution of the linearized Boltzmann equation is found using the relaxation time approach. The relaxation time ansatz is possible for elastic collisions, which is the case in the limit $m i / m e ≫ 1$⁠, but not for e–e collisions.

For the Lorentz model plasma, the electrical conductivity $σ Lorentz ( T , n )$ can be calculated beyond perturbation theory. The effective interaction of electrons with the ion subsystem is treated, for example, using the density functional theory. Similar approaches are given by the average atom model or the Lee–More model. These approaches are applicable in the high-density, degenerate case (⁠ $Θ ≤ 1$⁠), in which electron–electron collisions are suppressed due to Pauli blocking. They fail to describe real plasmas such as hydrogen plasmas in the low-density range, as e–e collisions are not taken into account.

Exact expressions for the conductivity $σ ( T , n )$ of plasmas are obtained from the generalized linear response theory (gLRT) in the form of correlation functions. The kinetic theory (KT) follows as a special case.³ However, the analytical evaluation of the correlation functions is only possible in some approximations using Green function diagram techniques. For example, the dynamic screening of the electrons is treated in the random phase approximation (RPA), the ionic structure factor is also approximated, and the Born approximation can be improved by a T-matrix approach.

Within the framework of a refined calculation, we must investigate the effects of dynamical screening, ion–ion structure factor, and strong collisions (see Refs. 7, 15, 33, and 34). These effects are of minor importance for the correction due to e–e collisions, both in the high-density limit where $lim Θ → 0 R e e ( Θ ) = 1$ and in the low-density limit, where they only occur in higher orders of the virial expansion. Although these corrections are small, we analyze them in Sec. IV B.

As shown in Fig. 1, the DFT-MD simulations for the effective resistivity $ρ 1 eff , DFT − MD ( T , n )$ are on a straight line pointing to the Lorentz limit (see Ref. 19). The corrected values $ρ 1 eff , ee − corr ( T , n ) = ρ 1 eff , DFT − MD / R e e$ are shown in Table II. The correction factor R_ee, Eq. (35), leads to a considerable increase (up to a factor $≈ 2$ in the low-density range) of the resistivity and solves the drawback of the correct ρ₁ benchmark by construction. However, the slope determining the second virial coefficient seems to remain too large compared to the high-temperature limit, Eq. (11), obtained from QLB or gLRT approaches.^15,17 While the slope of the DFT-MD simulations is 0.9965, the slope of the e–e corrected virial line is 0.732. The slope of the virial line in Fig. 1 determines the second virial coefficient $ρ 2 ( T )$⁠. The e–e corrected value is closer to the high-temperature benchmark 0.4917, but remains different.

The e–e correction factor, Eq. (35), was obtained from a two-moment approximation. The inclusion of higher moments of the distribution function would change the e–e correction factor and the corresponding values of $ρ 1 , ρ 2 ( T )$⁠. These corrections are expected to be small, see the good convergence with increasing number of moments reported in Ref. 15.

The virial coefficient $ρ 2 ( T )$ is determined by the screening of the Coulomb potential. The Born approximation can be applied in the high temperature range $T Ha ≫ 1$ so that strong collisions leading to a T-matrix approach are not relevant. A reason for the ρ₂ discrepancy can also be found in the treatment of screening, see the analysis of the slope parameter in Ref. 31. The electrons should be treated as dynamic screening, in lowest approximation by the RPA expression for the polarization function, to obtain the QLB benchmark Eq. (11). It is not clear to what extent the dynamic screening by electrons is already implemented in the DFT-MD simulations by the choice of the exchange-correlation part E_xc of the energy-density functional.

If $σ Lorentz ( T , n )$ which is required for the calculation of R_ee, Eq. (35), treats the dynamical screening of the electron–ion interaction differently than the DFT-MD approach $σ DFT − MD ( T , n )$⁠, the use of the correction factor is not consistent since the slope parameter $ρ 2 ( T )$ is not corrected. With respect to the contribution of the ions to the screening, the ion structure factor describes the static screening, which is well approximated in the DFT-MD simulations.

In conclusion, the correction factor (35) given in Ref. 8 reproduces the general behavior for $Θ ≤ 1$ on $R e e ≈ 1$ and is a reasonable approximation to the correct Spitzer result in the low-density limit. Only a two-moment approximation of the electron distribution function is considered to calculate R_ee, Eq. (35). It can be improved by considering higher moments of the electron distribution function so that the exact value of the Spitzer benchmark is reproduced. With respect to the second virial coefficient, the specific choice of the e–i interaction, which defines the Lorentz model plasma as a reference system, is discussed in Sec. IV C.

While electron–ion collisions are well understood to calculate the conductivity of plasmas, the inclusion of electron–electron collisions is notoriously difficult. One reason for this is that a relaxation time approach is only valid for elastic collisions where the single-particle energy is conserved, what holds to a good approximation for e–i collisions, but not for e–e collisions. Several approaches have been worked out to implement e–e collisions. At low densities, this problem is solved within the framework of kinetic theory. There, a linearized Boltzmann equation can be derived for a weak electric field.^1,5 The perturbation of the single electron distribution function results from the Kohler variational principle. In dense plasmas, gLRT is used to express the conductivity in terms of equilibrium correlation functions.^3,6

To calculate the conductivity of degenerate electrons in condensed matter, the Ziman formula is known, which considers the force–force autocorrelation function. Calculations of the collision frequency based on the solution of the electron–ion scattering problem are improved by a renormalization factor $r ( ω )$ to obtain the correct dynamic collision frequency of the plasma.^8,34,35 While the Ziman formula is frequently used to describe degenerate Coulomb plasmas such as electrons in liquid metals or solids, the renormalization factor allows to treat also non-degenerate plasmas, especially in the low-density range. The method of including e–e correlations by means of a renormalization factor $r ( ω )$ is equivalent to using the correction factor R_ee, Eq. (34), and requires analogous approximations.

An improved approach to treat dense plasmas is to start from the Lorentz model plasma, in which the e–i interaction is solved, and the electron distribution function is rigorously obtained from the relaxation time ansatz for any degree of degeneracy. Starting from the linearized Boltzmann equation, thermoelectric transport coefficients for Lorentz model plasmas are calculated, see the Lee–More model.³⁶ For a dense system, the averaged atom model was worked out, which describes the average influence of the surrounding plasma on the atom. The interaction of the electrons with the ion system is better described by the DFT approach, and the transition to condensed matter is well established. The electron–electron interaction is only considered in mean-field approximation to introduce an effective potential for the single-electron quasiparticle states.

The success of relaxation-time models for transport coefficients to treat the electron–ion interaction has led to develop a generalization of these models which includes e–e collisions. For example, as in Refs. 30 and 37, the relaxation-time model by Starrett³⁸ includes Pauli blocking and accounts for correlations using a mean-force scattering cross section for electron–ion collisions, but electron–electron collisions are accounted for only through the correction formula for R_ee proposed in Ref. 8.

A possible solution for this particular problem of including e–e collisions in the determination of plasma conductivity can be the use of path-integral Monte Carlo (PIMC) simulations to calculate $σ ( T , n )$⁠. In PIMC simulations, e–e correlations are treated in an exact way. This enables the immediate calculation of $R e e ( T , Θ )$ according to (34) and the construction of consistent interpolation formulas. In particular, the determination of virial coefficients would be of interest. Such PIMC simulations are not yet possible. Highly accurate results are only available for the uniform electron gas.²³ In addition, there is the unsolved sign problem, and PIMC simulations are currently performed for too small numbers of particles (a few tens). It is expected that significant results for the conductivity of hydrogen plasmas from PIMC simulations will be available in the near future. This will overcome the drawback of DFT-MD simulations where e–e correlations are not treated precisely.

We propose to use the virial expansion (8) for the determination of the conductivity of the hydrogen plasma in order to analyze the consistency of different theoretical approaches using the virial coefficients as a benchmark. Having $ρ 1 , ρ 2 ( T ) ,$ etc. available, we can find improved versions of the correction factor $R e e ( T , n )$⁠. Since these correction factors relate the plasma conductivity $σ ( T , n )$ to the conductivity of a Lorentz model plasma $σ Lorentz ( T , n )$⁠, see Eq. (34), these correction factors depend on the definition of the Lorentz model Hamiltonian. They depend, in particular, on the definition of the effective electron–ion potential, which includes dynamic or static screening by electrons, or on the particular choice of the exchange-correlation energy E_xc in DFT-MD calculations. If the determination of the conductivity and virial expansion for the Lorentz model plasma is known, a corresponding correction factor can be defined.

We propose an exact virial expansion (8) for the plasma conductivity to analyze the consistency of theoretical approaches. The virial coefficients serve as a benchmark for different approaches and can be better visualized using virial plots. For example, various analytical calculations of the dc conductivity $σ ( T , n )$⁠, which were presented in Ref. 20, do not fulfill this exact requirement and do not provide accurate results in the range of small densities. The results of DFT-MD simulations, which are currently considered to be the most reliable, are checked by virial plots, and shortcomings are found in the region of small densities. By benchmarking with the virial expansion (8) for $x → 0$⁠, future PIMC simulations can also be tested. These ab initio simulations become a computational challenge in the low-density region, but the virial expansion enables extrapolation into this range. For the construction of interpolation formulas, see Ref. 18, knowledge of the virial coefficients is an important ingredient.

A particular problem concerns the treatment of the e–e interaction in DFT-MD simulations. To a certain extent, the e–e interaction is included in the exchange-correlation energy E_xc in order to obtain an effective mean field for optimal single-particle states. However, with the help of virial expansion, it was shown that e–e collisions are not included in the DFT-MD calculations. By introducing a correction factor R_ee, which is obtained from a model plasma with statically screened Coulomb interaction, the first virial coefficient ρ₁ can be reproduced. In this work, the problems with the second virial coefficient $ρ 2 ( T )$⁠, which is determined by the treatment of dynamic screening, are emphasized. Since the dynamic screening is not considered in the DFT-MD simulations, the second virial coefficient $ρ 2 QLB$⁠, Eq. (11), is also not obtained by the correction factor R_ee. A numerical calculation can be considered in the context of PIMC simulations to determine the second virial coefficient $ρ 2 ( T )$⁠.

A correction factor $R e e DFT − MD ( T , n )$ was introduced, which takes into account the specific choice of exchange-correlation energy in the DFT approach. The high-temperature limit of the virial expansion was given in Eq. (38), and the full T dependence would be of interest. The DFT-MD simulations are very accurate at higher densities since e–e correlations are strongly suppressed due to Pauli blocking in the degenerate range (⁠ $Θ ≤ 1$⁠). With the knowledge of $R e e DFT − MD ( T , n )$⁠, results for the conductivity of hydrogen plasmas can also be obtained from DFT-MD calculations for the non-degenerate range. In addition to the analytical evaluation of correlation functions for higher moments of the distribution function in the framework of gLRT, the numerical calculation in the context of PIMC simulations can be considered to determine the second virial coefficient $ρ 2 ( T )$⁠.

The approach described here can also be applied to other transport properties such as thermal conductivity, thermopower, viscosity, and diffusion coefficients.²² It is also interesting to extend the virial expansion to substances other than hydrogen, where different ions can be formed and the electron–ion interaction is no longer purely Coulombic.

The author would like to thank M. Bethkenhagen, M. French, R. Redmer, H. Reinholz, and M. Schörner for valuable and fruitful discussions.

The author has no conflicts to disclose.

G. Roepke: Conceptualization (lead); Investigation (lead); Writing – original draft (lead).

The data that support the findings of this study are available within the article.