Knowledge of the thermodynamic, transport, and optical properties of warm dense matter (WDM) is of outstanding importance for planetary science.1 However, a dense plasma (or WDM) under planetary interior conditions is strongly correlated and degenerate, and so simple models to calculate these properties fail. Therefore, material properties need to be determined preferably on the basis of ab initio simulations by using, for example, density functional theory (DFT) for the electrons coupled to classical molecular dynamics (MD) for the ions. This complex program has been executed recently, for example, along the pressure–temperature–density (P–T–ρ) profiles for the interior of Jupiter,2 Saturn,3 and brown dwarfs4 via large-scale DFT-MD simulations. Other extensive ab initio calculations have been performed for warm and hot dense matter, with the aims, for example, of generating an equation-of-state table of deuterium for inertial confinement fusion applications5 and of providing a so-called first-principles equation-of-state (FP-EOS) collection for a number of elements and compounds.6

Space missions such as Voyager and Galileo, the recent Cassini–Huygens mission to Saturn that ended in 2017, and the Juno spacecraft still orbiting Jupiter have provided a wealth of new insights into the properties of the gas giants in our solar system. In particular, gravity data of so far unprecedented accuracy have revolutionized our understanding of the interiors of Jupiter7 and Saturn,8 which is intimately connected with models for the formation and thermal evolution of these planets. For instance, there is now clear evidence for the existence of inhomogeneous regions in these interiors, which result from hydrogen–helium demixing on the one hand and dissolution of higher-Z material contained in the core region into the mantle of metallic hydrogen on the other.1 As well as a more detailed picture of the interior structure, advanced models of the dynamo action in the deeper layers that generates Jupiter’s magnetic field9 and of the zonal wind patterns10 have been established on the basis of the Juno data.11–13

Measuring material properties such as electrical and thermal conductivity and the viscosity of matter under such extreme conditions remains a great challenge for experimental high-pressure physics. Nevertheless, these quantities are an essential input into magnetohydrodynamic simulations describing the dynamo action in the deep interior of gas giants, i.e., the generation and evolution of their magnetic fields.14–17 These data are also important when calculating the influence of Ohmic dissipation by deep-seated zonal winds in the outer atmospheres of hot Jupiters.18–20

Therefore, numerous ab initio studies have been performed in the last decade to calculate material properties under such extreme conditions for the most abundant elements hydrogen and helium (see, e.g., Refs. 2–6 and 21–24). Comprehensive reviews of the properties of warm dense hydrogen and helium have been provided by McMahon et al.25 and Bonitz et al.26

The transport properties of fully ionized hydrogen plasmas have been calculated for a wide range of densities and temperatures, in particular in the context of studies of stellar structure and evolution (see, e.g., Refs. 27–30). To save computational costs, the development of elaborate conductivity models applicable to plasmas in a wide range of density and temperature is a very important task. These models should reproduce the known limiting cases: the Spitzer theory31 for low densities and high temperatures and the Ziman formula32 for high densities and low temperatures. Conductivity models include free parameters that are usually adjusted with reference to available theoretical and experimental data so that they can be applied for arbitrary densities and temperatures.

A well-known example is the quantum-statistical approach of Ichimaru and Tanaka33 based on a general solution of the Ziman formula for arbitrary degeneracy. Röpke and Redmer34 provided a virial expansion for the electrical conductivity that is valid for classical fully ionized plasmas. The underlying linear response approach35 was later generalized to the degenerate domain, and, in particular, the influence of electron–electron scattering was studied in detail.36 Other conductivity models were derived by, for example, Stygar et al.37 for nondegenerate Lorentz plasmas in magnetic fields and by Fortov et al.38 for strongly correlated plasmas on the basis of experimental data for hydrogen and noble gases. Theoretical approaches to the description of charged particle transport in dense plasmas and for applications in astrophysics and inertial confinement fusion research have been compared recently, spanning methods from semi-empirical models to ab initio calculations.39,40

Another way to construct a user-friendly conductivity model is by solving the Boltzmann equation in relaxation-time approximation (RTA) for arbitrary degeneracy, as Lee and More (LM)41 did in their well known study in 1984. By considering scattering cross sections in the Born approximation and applying simple cutoffs in the corresponding Coulomb logarithm, they provided expressions for the complete set of transport coefficients for plasmas in electric and magnetic fields for a given ionization state Z*, which can be used for many practical applications. For weakly correlated high-temperature plasmas, the Spitzer limit with the prefactors of the Lorentz plasma42,43 follows because electron–electron collisions are not accounted for in the RTA. The LM model was later generalized by Desjarlais44 to include the effects of partial ionization and of electron–electron scattering. However, on using the LM model for WDM, one finds huge deviations from DFT-MD data,23 which underlines the need to improve this and similar conductivity models for those conditions where strong collisions, electronic correlations, and the ionic structure factor have to be considered.

In this work, we describe an approach to improving the LM conductivity model for the case of a hydrogen plasma in the WDM regime. We use an extensive conductivity data set for the electrical and thermal conductivity based on DFT-MD simulations. Finally, we provide convenient interpolation formulas that can be used for applications in plasma and astrophysics as mentioned above.

We use an ab initio method based on the Born–Oppenheimer approximation. We describe our system by combining classic MD for the ions and DFT for finite temperatures45–47 for the electrons. For the DFT-MD simulations, we use the Vienna Ab initio Simulation Package (VASP)48–50 in version 5.4.4. Exchange and correlation (XC) effects of the electrons are approximated with the XC functional of Perdew, Burke, and Ernzerhof (PBE).51 In our simulation, we use a projector augmented wave (PAW) potential52–54 (labeled PAW H_h_GW) with a plane-wave cutoff of 1200 eV. For all simulations, we calculate the electronic wave function at the Baldereschi mean value point.55 The temperature is controlled by a Nosé–Hoover thermostat.56,57 Depending on the density and temperature of the system, we use up to 1400 atoms in our simulation box to get converged conductivity results and use a time step of 0.05–0.2 fs. The simulations run at least one ps up to several ps to provide a sufficient number of steps for the averages in equilibrium. The density and temperature vary from ρ = 0.1 to 8 g/cm³ and from T = 2000 to 500 000 K.

To calculate the electronic transport properties, we run static DFT calculations for 10–200 ionic configurations from each equilibrated simulation and average the results over the number of configurations. Here, we use the Coulomb potential (instead of the PAW potential) with a cutoff of 2000 eV, which has also been done in previous work to get well-converged conductivity data.23 We use several different Monkhorst–Pack k-point sets58 to check the convergence behavior with the number of k-points.

The following expression for the frequency-dependent Onsager coefficients L_n(ω) can be derived from linear response theory:21,22Ln(ω)=2π(−1)nVω∑kνμRe(⟨kν|v̂|kμ⟩⊗⟨kμ|v̂|kν⟩)×(fkν−fkμ)Ekμ+Ekν2−hen×δ(Ekμ−Ekν−ℏω),(1)where ω is the frequency, V is the volume of the simulation box, 〈kν|v̂|kμ〉 are the matrix elements of the velocity operator v̂, with |kμ⟩ the Bloch states, f_kμ are the Fermi occupation numbers, E_kμ are the energy eigenvalues, h_e is the enthalpy per electron, and ℏ is the reduced Planck constant. The matrix elements are calculated from the dipole matrix elements 〈kν|r̂|kμ〉, which are taken directly from the optical routines of VASP.22,59 To calculate transport properties from the Onsager coefficients, a Gaussian function is used to broaden the δ function to a small finite width. In our simulations, we select a broadening of 0.003–0.03 eV to get reasonable results.

The static electrical conductivity σ_e and thermal conductivity λ_e of the electrons are then given byσe=e2limω→0L0(ω),(2)λe=1Tlimω→0L2(ω)−L1(ω)L0−1(ω)L1(ω),(3)where e is the elementary charge. The coefficient L₀(ω) is also known as the ω-dependent Kubo–Greenwood formula.60,61 In practice, the static values are calculated by using a linear extrapolation to ω = 0 across an unphysical decrease at very small ω. This is typically done in the energy region from 0.1 to 1 eV.

Another way to calculate transport properties of partially ionized plasmas is to develop an appropriate multicomponent plasma model.22 Solving the Boltzmann equation in the RTA, Lee and More41 derived the following expression for the static Onsager coefficients:Kn=−2me3π2ℏ3∫0∞dEEn+1∑snsQes(E)∂fe0∂E,(4)which are directly connected to the electrical and thermal conductivities byσe=e2K0,λe=1TK2−K12K0.(5)The Onsager coefficients K_n depend on the sum of transport cross sections Q_es(E) and the derivative of the Fermi function. The transport cross sections describe how strongly electrons (e) are scattered by a specific scattering partner (s) for a given energy. For each scattering partner, one can write the specific Onsager coefficient in the form22Kn,es=−2me3π2ℏ3ns∫0∞dEEn+1Qes(E)∂fe0∂E(6)or, by using the relaxation time τ,23Kn,es=−23/2me1/23π2ℏ3∫0∞dEEn+3/2τes(E)∂fe0∂E,(7)where the specific relaxation time τ_es = 1/(vQ_esn_s) for the scattering partner s. In the following, we use Eq. (7) for electron–ion scattering and Eq. (6) for electron–neutral scattering to calculate the specific electrical conductivity σ_es and thermal conductivity λ_es. The total conductivity values are then calculated according to the Matthiessen ruleσe−1=∑sσes−1,λe−1=∑sλes−1.(8)This Matthiessen-rule-like ansatz for the thermal conductivity is motivated by theories from liquid mixtures62 and a previous approach to add electron–electron scattering to DFT results.63

In most conductivity models, electron–ion scattering is treated with a screened Coulomb potential using, for example, the Debye–Hückel screening length to derive an expression for the relaxation time in the Born approximation.22,23,64 This leads to a dependence on the so-called Coulomb logarithm, which is approximated via certain points that describe, for example, the effective minimal and maximal scattering distance from the ions.65

Here, we adopt a simple expression for the relaxation time,τeI(E)=E3/2P(T,ρI)nIZI2,(9)which has the same dependence on energy, particle density and charge as the relaxation time in the Born approximation. The fitting parameter P(T, ρ_I) is used to provide an accurate match of the LM model, which describes weakly coupled systems at high temperatures, to our DFT results for stronger coupled systems at lower temperatures so that both cases are included. This parameter depends only on the temperature T and the ion density ρ_I and not on the energy, and so we can take it out of the integral in Eq. (7) and get a simple expression for this Onsager coefficient:Kn,eI=−23/2me1/2P(T,ρI)3π2ℏ3nIZI2∫0∞dEEn+3∂fe0∂E.(10)

For the scattering of electrons by atoms and nonpolar molecules, we adopt the expression for the Onsager coefficients given by French and Redmer:22Kn,eN=211/2π1/2(n+3)!ε02(kBT)n+3/2ne3e4me1/2nN∑iZN,i2⁡lnAN,i(xn,N,i),(11)with the atomic logarithm functionlnAN(E)=ln(xN+1)2−3aN2+6aNbN−2bN26aN2+4aNbN+aN22aN2(xN+1)−bN2+aNbNaN2(xN+1)2+2bN23aN2(xN+1)3.(12)The parameters Z_N, a_N, and b_N are used to fit the transport cross section of electrons on neutral species N; they are also taken from Ref. 22. The logarithm function is evaluated atxn,N,i=8(n+2)(n+3)bN,i2mekBT/ℏ2.(13)With this expression, one can calculate all possible processes of scattering of electrons by neutral atoms and nonpolar molecules. Here, we assume scattering by hydrogen atoms only. The influence of scattering by hydrogen molecules becomes important for low-temperature plasmas, which are not studied here. In the temperature and density region considered in this work, hydrogen atoms are either ionized or single atoms, which are mostly free or bonded for very short times with a second hydrogen atom (transient molecules). This feature permits us to neglect dissociation processes. Thus, the number of neutral scattering partners is defined by the ionization degree α.

We assume a system of hydrogen atoms H, protons H⁺, and electrons e, and so the ionization degree is defined asα=NefreeNetotal=NH+NH++NH,(14)where NH++NH=NH,total is the total number of protons or hydrogen atoms. The density of the system is also directly connected to the total number of hydrogen atoms:ρ=ρH+ρH++ρe≈ρH+ρH+=ρH,total.(15)We use these two quantities to define the densities of ions and of neutral atoms, which are needed in Eqs. (10) and (11):ρI=ρH+=αρH,total,(16)ρH=(1−α)ρH,total.(17)The number density is then given by n_i = ρ_i/m_H, with i = {H, H⁺} and m_H the mass of the hydrogen atom. For hydrogen, the number density of free electrons n_e is equal to the number density of ions (protons). In this work, α is calculated by matching our conductivity model to the DFT data and not by applying a chemical model such as the Saha equation to determine the plasma composition separately (see Sec. III B).

First we investigate the simplest case of a fully ionized hydrogen plasma, where all atoms are ionized and only electron–ion scattering occurs. We use the electrical conductivity results from the DFT simulations using the Kubo–Greenwood formula and the results from the LM model with a screened Coulomb logarithm that includes degeneracy23 for high temperatures to fit Eq. (10). For that purpose, we employ the following expression for the fitting parameter P(T, ρ_I):P(T,ρI)=a0Ta1ρa2+a3ρa4+a5+a6ρa7⁡exp−log⁡T−a8a92,(18)where we insert the values for T in kelvin and ρ in g/cm³. The coefficients a₀, …, a₉ are given in Table I. The first three terms of the fitting parameter P(T, ρ_I) describe the high-temperature behavior of the LM model, while the last term describes the differences between the high-temperature LM model and the DFT results for lower temperatures where stronger interactions have to be taken into account. The resulting electrical and thermal conductivities are shown in Figs. 1 and 2.

For temperatures above 10⁶ K, we consider a weakly coupled plasma and use the results from the LM model. If we ignore electron–electron scattering, the electrical conductivity exhibits Lorentz plasma behavior for a weakly coupled plasma at high temperature. For temperatures below 10⁶ K, we use the DFT conductivity to capture stronger correlations in our fit. The curves in Fig. 1 representing the results using the fitting parameter P(T, ρ_I) nearly overlap the circles representing the DFT and RTA results. This indicates that it is possible to modify the LM model for a fully ionized one-component plasma by using an improved relaxation time that includes the effects of strong correlations, which are essential in WDM. Nevertheless, it is questionable whether the fitting parameter (18) also holds for the very high densities typical of stellar cores or the interiors of compact objects such as white dwarfs. The limits of our model are discussed in Secs. III B and III C. In general, the model describes metallic behavior for lower temperatures, where the stronger vibrations of the ions with temperature lower the conductivity until the system runs into the Spitzer or, here, the Lorentz limit, where the conductivity increases with temperature again.

We can also obtain correct results for the thermal conductivity. In Fig. 2, we show only the region where we compare our fitting parameter with the DFT data. The agreement is again very good. Note that our model slightly overestimates the results for both conductivities for lower densities, but all deviations are within 10%, which is a typical estimate for the uncertainty of DFT conductivities in the WDM domain.

Ionization or dissociation effects have not been included in our model so far, although these are essential for lower temperatures and lower densities. The EOS and the composition of partially ionized hydrogen plasmas have long been treated by chemical models. For instance, Ebeling and Richert constructed efficient Padé formulas valid for the thermodynamic functions of fully ionized hydrogen,66 which were used to predict the plasma composition via chemical models (see, e.g., Ref. 67). Fluid variational theory (FVT) was applied successfully to infer EOS data for shock-compressed liquids68 and, in particular, for the description of the dissociation equilibrium in dense hydrogen.69,70 The free-energy model of Saumon, Chabrier, and van Horn71 treats the dissociation and ionization equilibrium in hydrogen self-consistently and has found widespread applications in plasma physics and astrophysics. Däppen et al.72–75 derived an EOS for hydrogen plasmas to model stellar envelopes. Chemical models have also been applied to evaluate explosively driven shock-wave experiments.76 The corresponding transport coefficients of partially ionized hydrogen plasmas have been calculated with the aim of locating the nonmetal-to-metal transition in density–temperature space (see, e.g., Refs. 77–79).

An important application of models of partially ionized hydrogen plasmas is to the outer layer of gas giant planets, which are cold and thin enough to consist of molecular hydrogen plus helium and a small amount of heavier elements. Increasing pressure and temperature toward the center of the planet drives dissociation and ionization—hydrogen transforms to a metal. This nonmetal-to-metal transition is of first order below a critical temperature of about 1500 K and a pressure of a few Mbar: a liquid–liquid phase transition (LL-PT). Above the critical temperature, the transition is abrupt but continuous. The location of the corresponding coexistence line has been discussed for many years; for a recent review, see Ref. 26. Calculations based on DFT locate the continuous nonmetal-to-metal transition at about 90% of Jupiter’s radius2 and 70% of Saturn’s radius,3 respectively. The conductivities and other material data are important to better understand the dynamo action in the deep interior of gas giants via MHD simulations (see, e.g., Refs. 16 and 17). Another example is Ohmic dissipation, which might occur in the atmospheres of hot Jupiters. This process is a candidate to explain the radius anomaly of many hot Jupiters.18,19,80

Therefore, we generalize the conductivity model derived in Sec. II B to a partially ionized hydrogen plasma. As already mentioned in Sec. II, we assume that the hydrogen plasma consists of atoms, protons, and electrons. This means that we treat all neutrals as atoms. Hydrogen molecules are consequently described by two individual hydrogen atoms. The ionization degree α is calculated from the formulaα=1−1b0Tb1ρ(b2Tb3+b4)+1(19)by fitting its value to the electrical conductivity from DFT calculations using the model for a fully ionized hydrogen plasma. Equation (11) for the contribution of electron–neutral interactions is then used to determine the total conductivity via Eq. (8). The parameters b_i in Eq. (19) are given in Table II. We insert again the values for T in kelvin and for ρ in g/cm³.

In both terms, the ionization degree α defines the number of electrons and the number of scattering partners. The calculated ionization degree that reproduces the DFT results within our model is shown in Fig. 3. The curves show the ionization degree α fitted using Eq. (19). Note that we have used some additional conductivity data for 1 g/cm³ at 1000 and 2000 K from Holst et al.21 to achieve an improved behavior for the ionization degree in that domain.

Applying the widely used but approximate Saha equations leads to too low ionization degrees for such dense plasmas, as can be seen from the dashed curves in Fig. 3. A rigorous method to determine the ionization degree α is based on the Thomas–Reiche–Kuhn sum rule for the dynamic electrical conductivity, which is calculated via the Kubo–Greenwood equation using DFT. This method provides reliable ionization degrees for high-density carbon81 and beryllium plasma82 and was also tested here. However, the sum-rule method does not hold for conditions where the bandgap starts to close or vanishes. For hydrogen, we found an unphysical decrease of the ionization degree with temperature for high densities just in that region. This problem can probably be resolved by using improved meta-generalized gradient approximation (meta-GGA) XC functionals like SCAN (Strongly Constrained and Appropriately Normed)83 or hybrid functionals like HSE (Heyd–Scuseria–Ernzerhof),84 which yield more realistic bandgaps compared with the PBE GGA XC functional used in this study.

Therefore, the strong deviations from the target values for the ionization degree in that specific region motivated us to use the fitting formula (19). This provides reliable results for the ionization degree for most densities. Note that the fit starts to deviate strongly in the low-density region where the Saha equation approach holds and only temperature-driven ionization is relevant. Therefore, we use the results of the Saha equation approach whenever they are higher than the ionization degree predicted by the fitting formula (19).

Using the values for the ionization degree determined from Eq. (19), we calculate the electrical and thermal conductivities, as shown by the curves in Figs. 4 and 5, where the circles show the results of the DFT simulations. We very accurately resolve the region where the hydrogen plasma is metallic and the ionization degree is high within our model. We also capture the transition to the nonmetallic region in the low-temperature limit. One can see that the deviations become more pronounced at low temperatures and low densities. Note that the densities are still too high for the Saha equation approach, which results in some order-of-magnitude differences. Nevertheless, the plasma has a very low conductivity in both cases, i.e., the number of free electrons and also the corresponding ionization degree α are very low.

Note that the ionic contribution to the currents has to be taken into account for low temperatures and low densities as are typical of, for example, the outer layer of a gas giant planet. Hydrogen is mostly atomic or molecular in this layer, and the thermal conductivity is dominated by the ionic part λ_i and not by the electronic contribution λ_e (see Refs. 2 and 3).

In this subsection, we benchmark our conductivity model with other computational methods and available experimental data for dense hydrogen over a wide range of temperatures and densities.

First, we compare our results for the electrical conductivity with those obtained by Lambert et al.85 They used an orbital-free treatment for the electrons in MD simulations to generate atomic configurations for high densities and temperatures. Then, they calculated the transport coefficients within the Kubo–Greenwood formalism [see Sec. II A and Eqs. (1)–(3)] for five snapshots by using DFT in the local density approximation (LDA). To include electron-electron interaction, they used a correction factor given by Bespalov and Polishchuk86 as described in Ref. 87. The comparison is presented in Fig. 6.

We find very good agreement with the data of Lambert et al.85 for a density of 10 g/cm³ and the LM model for high temperatures. For higher densities and lower temperatures, we get considerably higher conductivity values. This deviation decreases with temperature, where the Spitzer limit is fulfilled for both methods. Note that we have included the same electron–electron correction factor as chosen by Lambert et al. in our model and the LM model, and so this is not the source of the deviations obtained. The deviations for very high densities, especially in the low-temperature region, may be traced back to a different treatment of the broadening or smearing parameter when calculating the Onsager coefficients. If the broadening parameter is chosen too high, the dynamic conductivity has a nonzero Drude-like behavior for ω → 0 (see Fig. 6 in Ref. 85), while a linear extrapolation should be performed toward the static limit in the region of the unphysical decrease (this paper). This would explain why the results of Lambert et al. are systematically lower than those obtained using a linear extrapolation method. Nevertheless, a deviation by a factor of 2 for lower temperatures can probably be traced back not only to the different extrapolation methods. Therefore, we have performed additional simulations for higher densities (20, 50, and 70 g/cm³) at T = 200 000 K and for 70 g/cm³ at T = 1 000 000 K. We find that our model overestimates the DFT conductivity values for higher densities. The differences amount to 50% at ρ = 70 g/cm³. The deviation decreases at higher temperatures. This trend indicates a stronger deviation for even higher densities at lower temperatures, where our model has to be used with caution. For high temperatures, we find better agreement with the LM model that describes the right Spitzer limit by including electron–electron collisions for all three densities. Thus, our model can be used for high densities in the high-temperature limit.

A rigorous formalism to include electron–electron collisions far beyond the simple model of Bespalov and Polishchuk86 was developed by Reinholz et al.36 within generalized linear response theory. The resulting correction factor R_ee that accounts for electron–electron scattering is shown in Fig. 7 for our density–temperature grid. We find that the electron–electron correction factor becomes important for very high temperatures at the densities considered here.

Electron–electron collisions are also important for the conductive opacity κ_c of astrophysical objects, which is directly connected to the thermal conductivity viaκc=16σBT33ρλ,(20)where σ_B is the Stefan–Boltzmann constant. For many applications in astrophysics, the conductive opacity is calculated with an effective collision frequency that depends on the Coulomb logarithm. Cassisi et al.88 and Blouin et al.89 provide fitting formulas for the opacity of fully ionized hydrogen plasmas over a wide range of densities and temperatures, including electron–electron scattering effects. Our model does not include electron–electron scattering for the thermal conductivity. Cassisi et al.88 provide opacity results for both cases, i.e., electron–ion and electron–ion + electron–electron scattering, while Blouin et al.89 give a correction factor for the Cassisi model at T > 0.1T_F, where T_F is the Fermi temperature. Combining the results from Figs. 6 and 7, a comparison is only applicable for temperatures T < 10⁵ K for lower densities and for very high densities around 100 g/cm³ where the electron–electron corrections are very small also for high temperatures. Much higher densities as are relevant in white dwarfs cannot be investigated within the current model, because in such a dense plasma relativistic effects have to be taken into account. Our results in the region where we trust our model provide good agreement with the results of Cassisi et al.88 for electron–ion scattering and also agree with the results of Blouin et al.,89 where we expect no or weaker electron–electron correction effects. Note that the decomposition for λ according to electron–ion and electron–electron scattering processes as introduced in Ref. 63 and used by Blouin et al.89 is not valid, as shown by French et al.23

On the contrary, corrections due to electron–electron scattering can be neglected for low temperatures of only few thousand kelvin, owing to Pauli blocking (degenerate region). In Fig. 8, we compare the predictions of our conductivity model with (i) the gas-gun experiments of Weir et al.90 and Nellis et al.91 inspecting the low-temperature region of a few thousand kelvin and (ii) the virial expansion of Röpke et al.,64,92 which represents the correct high-temperature limit. We find very good agreement with the high-temperature limit given by the virial expansion and also match the experimental results in the lower-temperature region quite well. We can also reproduce the density between 0.6 and 0.7 g/cm³ as determined in the experiment.

We conclude that our conductivity model predicts the behavior of dense hydrogen over a wide range of temperatures. Other experiments by, for example, Ternovoi et al.93 and Fortov et al.38,76 show the same jump in conductivity with density, but at slightly higher temperatures. This deviation can be traced back to the determination of the temperature in these experiments.

We also compare our model with results from an average atom (AA) model94 in the low-density regime, which are shown in Fig. 42 of Ref. 26. The DFT results and the AA data match nearly perfectly for densities above 0.1 g/cm³, which was also shown in previous work.39 For lower densities, it is very difficult and numerically expensive to converge the DFT-MD simulations. Here, we can use our model to extrapolate to very low densities. In Fig. 9, we show the low-density behavior of our model compared with the AA model for temperatures of 1.35, 2.69, and 5.39 eV. We find the same general trends as predicted by the AA model, but the deviations are greater than for higher densities. For all three temperatures, one can see a kink in the data set of our model. At this point, we change the value for α from our fit via Eq. (19) to the value given by the Saha equation. This leads to higher ionization with decreasing density. The same behavior is observed in the AA model at nearly the same densities. Combining the results from Figs. 8 and 9, we can predict the general behavior of the ionization degree and the conductivity in the low- and high-density regions.

We have calculated an extensive DFT-MD data set for the electrical and thermal conductivity of hydrogen in a wide range of density ρ and temperature T. Based on these ab initio results, we have proposed a modified LM model that provides improved predictions for the electrical and thermal conductivities for a dense hydrogen plasma in the fully and partially ionized region by just changing the relaxation time. In addition, we have introduced an expression for the ionization degree α to generalize our model from the fully to the partially ionized plasma region.

Our new conductivity model shows in general very good agreement with other computational results and available experimental data over a wide range of densities and temperatures. The underlying ab initio data and the proposed fitting formula can be further benchmarked by future experiments to provide a better understanding of the behavior of warm dense hydrogen. This is essential for improved predictions for the interior structure, thermal evolution, and dynamo action in, for example, gas giant planets. Our model can also be used for conditions in stellar interiors exceeding densities of ρ = 100 g/cm³ and for temperatures up to T = 1000 eV.

Acknowledgment. We thank Martin French for motivating this study and for continual interest in this project. Furthermore, we thank Armin Bergermann, Martin Preising, Argha Roy, and Gerd Röpke for many helpful discussions. We thank Jan Vorberger and Stephanie B. Hansen for providing the AA results for the electrical conductivity shown in Fig. 9. This work was supported by the Priority Program SPP 1992 of the German Science Foundation (DFG) The Diversity of Exoplanets under project number 362460292. The ab initio calculations were performed at the North-German Supercomputing Alliance (HLRN) facilities and at the IT and Media Center of the University of Rostock.