A new equation of state for dense hydrogen-helium mixtures

TL;DR

This paper introduces a comprehensive new equation of state for dense hydrogen-helium mixtures, integrating multiple models and ab initio calculations to cover a wide range of astrophysical conditions and providing entropy data crucial for stellar evolution.

Contribution

It presents a novel EOS for H/He mixtures that combines existing models with ab initio simulations and includes entropy calculations across extensive density and temperature ranges.

Findings

The new EOS aligns well with experimental shock data.

It matches first-principles simulations for pure elements and mixtures.

The EOS is applicable from planetary to stellar densities.

Abstract

We present a new equation of state (EOS) for dense hydrogen/helium mixtures which covers a range of densities from $10^{-8}$ to $10^6$ g.cm$^{-3}$, pressures from $10^{-9}$ to $10^{13}$ GPa and temperatures from $10^{2}$ to $10^{8}$ K. The calculations combine the EOS of Saumon, Chabrier & vanHorn (1995) in the low density, low temperature molecular/atomic domain, the EOS of Chabrier & Potekhin (1998) in the high-density, high-temperature fully ionized domain, the limits of which differ for H and He, and ab initio quantum molecular dynamics (QMD) calculations in the intermediate density and temperature regime, characteristic of pressure dissociation and ionization. The EOS for the H/He mixture is based on the so-called additive volume law and thus does not take into account the interactions between the two species. A major improvement of the present calculations over existing ones is that we calculate the entropy over the entire density-temperature domain, a necessary quantity for stellar or planetary evolution calculations. The EOS results are compared with existing experimental data, namely Hugoniot shock experiments for pure H and He, and with first principle numerical simulations for both the single elements and the mixture. This new EOS covers a wide range of physical and astrophysical conditions, from jovian planets to solar-type stars, and recovers the existing relativistic EOS at very high densities, in the domains of white dwarfs and neutron stars.

Figures

Tables

Table 1. Table 1 Example of the EOS table. Units are T 𝑇 T in K, P 𝑃 P in GPa, U ~ ~ 𝑈 {\tilde{U}} in MJ kg -1 , S ~ ~ 𝑆 {\tilde{S}} in MJ kg -1 K -1 . Each ( log ⁡ T , log ⁡ P 𝑇 𝑃 \log T,\log P ) or ( log ⁡ T , log ⁡ ρ 𝑇 𝜌 \log T,\log\rho ) table includes N T = 121 subscript 𝑁 𝑇 121 N_{T}=121 isotherms, each with N P = 441 subscript 𝑁 𝑃 441 N_{P}=441 pressure values or N ρ = 281 subscript 𝑁 𝜌 281 N_{\rho}=281 density values, with step values 0.05 for Δ ​ log ⁡ T Δ 𝑇 \Delta\log T , Δ ​ log ⁡ ρ Δ 𝜌 \Delta\log\rho and Δ ​ log ⁡ P Δ 𝑃 \Delta\log P . Note: the number of digits after the point has been truncated in this exemple to fit the journal format. In the online table, all quantities are given with 6 digits after the point.

$\log T$	$\log P$	$\log \rho$	$\log \tilde{U}$	$\log \tilde{S}$	${\left(\frac{\partial \log \rho }{\partial \log T}\right)}_P$	${\left(\frac{\partial \log \rho }{\partial \log P}\right)}_T$	${\left(\frac{\partial \log S}{\partial \log T}\right)}_P$	${\left(\frac{\partial \log S}{\partial \log P}\right)}_T$	${\nabla }_{ad}$
#iT=1 log T= 2.000
0.200E+01	-0.900E+01	-0.909E+01	-0.628E-01	-0.100E+01	-0.114E+01	0.983E+01	0.242E+02	-0.355E+00	0.300E+00
0.200E+01	-0.895E+01	-0.870E+01	-0.628E-01	-0.102E+01	-0.591E+01	0.572E+01	0.190E+02	-0.297E+00	0.300E+00
…	…	…	…	…	…	…	…	…	…
…	…	…	…	…	…	…	…	…	…
0.200E+01	0.129E+02	0.542E+01	0.781E+01	0.547E+01	0.454E+00	0.500E+01	-0.152E+01	0.506E+00	0.330E+00
0.200E+01	0.130E+02	0.544E+01	0.784E+01	0.550E+01	0.455E+00	0.498E+01	-0.152E+01	0.506E+00	0.330E+00
0.205E+01 -	-0.900E+01	-0.978E+01	-0.174E-01	-0.402E+00	-0.160E+02	0.161E+02	-0.260E+01	-0.315E+01	0.300E+00
0.205E+01	-0.895E+01	-0.911E+01	-0.174E-01	-0.548E+00	-0.105E+02	0.104E+02	-0.168E+01	-0.267E+01	0.300E+00
…	…	…	…	…	…	…	…	…	…
…	…	…	…	…	…	…	…	…	…
#iT=2 log T= 2.050
0.205E+01	-0.900E+01	…	…
…	…	…	…	…	…	…	…	…	…
…	…	…	…	…	…	…	…	…	…
0.205E+01	0.130E+02	0.547E+01	0.781E+01	0.542E+01	0.451E+00	0.498E+00	-0.152E+01	0.507E+00	0.331E+00
#iT=3 log T= 2.100
0.210E+01	…	…
…	…	…	…	…	…	…	…	…	…
…	…	…	…	…	…	…	…	…	…

Equations

$$T_m=T_0(1+P/a{)}^{b}exp(-cP),$$

$$\log T_m\approx 3.1+\frac{1}{3}\log \rho +2\log Z-\frac{1}{3}\log A,$$

$$\log T\approx 3.3+\frac{1}{2}\log \rho +\log Z-\frac{5}{3}\log A,$$

$$(\tilde{E}-{\tilde{E}}_0)+\frac{1}{2}(P+P_0)\times ({\rho }^{-1}-{\rho }_0^{-1})=0$$

$$\widetilde{C_P}$$

$$\widetilde{C_V}$$

$$\frac{{\kappa }_T}{{\kappa }_T^0}=\frac{\rho k_{B}T}{\mu P}{\chi }_{\rho }^{-1},$$

$$T_m=61.0P^{0.639}\mathrm{K},$$

$$W(T,P)=\sum _i{X}_i{W}_i(T,P),$$

$$\frac{1}{\rho (T,P)}=\frac{1-Y}{{\rho }_H(T,P)}+\frac{Y}{{\rho }_{He}(T,P)},$$

$$\tilde{S}(T,P)=\sum _i{X}_i{\tilde{S}}_i(T,P)+{\tilde{S}}_{mix}^{id}(T,P).$$

$$\frac{{\tilde{S}}_{mix}^{id}}{k_B}$$

$${\nabla }_{ad}$$

$${\nabla }_{ad}(X,P,T)=\frac{{\sum }_i{X}_i{S}_i(\frac{\partial \log S_i}{\partial \log P}{)}_T}{{\sum }_i{X}_i{S}_i(\frac{\partial \log S_i}{\partial \log T}{)}_P},$$

$$\alpha$$

$$C_S$$

$$\frac{{\kappa }_T}{{\kappa }_T^0}$$

$${\kappa }_S$$

$$\gamma$$

$$2.0$$

$$-9.0$$

$$-8.0$$

Full text

1. Introduction

Understanding the thermodynamic properties of hydrogen (H) and helium (He) at high density, characterized by their equation of state (EOS), is at the heart of numerous physical and astrophysical problems. From the point of view of fundamental physics, understanding the metalization of hydrogen has remained a major challenge since the pioneering work of Wigner & Huntington (1935), more than 80 years ago. The quest for its observational evidence remains so far unachieved but is in reach with both static and dynamic high pressure experiments, thanks, for these latter, to the achievement of modern techniques such as Z pinch magnetically driven shock experiments (Knudson et al. 2004), spherically converging shock wave experiments (Belov et al. 2002, Boriskov et al. 2003) and intense laser driven planar shock wave experiments (Collins et al. 1998, Hicks et al. 2009, Sano et al. 2011, Loubeyre et al. 2012, Brygoo et al. 2015). These experiments have revealed the principal Hugoniot of dense deuterium up to 200 GPa. Knowledge of the hydrogen and helium EOS is also central for inertial confinement fusion (ICF) and of course for the characterization of the interior or outer mechanical and thermal structures of dense astrophysical bodies. These latter include low-mass stars (generically stars smaller than the Sun, for which the perfect gas EOS or the Debye-Hückel expansion is no longer valid), brown dwarfs (objects not massive enough to sustain or even ignite hydrogen fusion in their core, whose mass distribution extends from about 0.07 ${\rm M}_{\odot}$ down to a few Jupiter masses), giant (solar and extrasolar) planets, but also the envelope of white dwarfs and the outer envelope and atmosphere of neutron stars.

In the meantime, ab initio numerical calculations of the properties of dense H and He, based either on quantum molecular dynamics (QMD), which combines molecular dynamics (MD) for the heavy classical particles and density functional theory (DFT) to treat the quantum electrons, or Path Integral Monte Carlo (PIMC) or Molecular Dynamics (PIMD) can now be performed in the density-temperature domain of interest (e.g. Militzer & Ceperley 2000, Holtz et al. 2008, Militzer 2009, 2013, Militzer et al. 2001, Morales et al. 2010ab, Lorenzen et al. 2009, 2011, Becker et al. 2014, Mazzolla et al. 2018, Schöttler & Redmer 2018), thanks to the enormous improvement in computer capacities. The widely used semi-analytical H/He model of Saumon & Chabrier (Chabrier 1990, Saumon & Chabrier 1991, 1992, Saumon et al. 1995 (SCvH)) can thus now be replaced by these calculations in the crucial domain of pressure ionization. Such an approach, combining ab initio calculations with the SCvH EOS in the low (mainly atomic/molecular) and high (fully ionized) domains, has been used by various authors (Caillabet et al. 2011, Militzer & Hubbard 2013, Becker et al. 2014, Miguel et al. 2016). These calculations, however, remain so far limited in two aspects. Either they cover only a limited density-temperature range, precluding the use of an EOS over a significant physical or astrophysical domain, or they do not provide the entropy. Indeed, while the pressure and internal energy are directly accessible to QMD or PIMC/PIMD calculations, the entropy is a much more cumbersome task, requiring a so-called thermodynamic integration over a large number of temperature and density points. The knowledge of the entropy, however, is central in stellar evolution calculations (as the cooling history of a star directly derives from the first principle of thermodynamics, $Q=dS/dt$) and even to determine the thermal structure of dense astrophysical bodies since their interior is quasi isentropic, due to the onset of convection to carry out their internal heat flux111It must be kept in mind that only for adiabatic reversible process, such e.g. as convection, is an adiabat ($dQ$=0) equivalent to an isentrope ($dS=0$). . The thermal profile and the contraction rate, thus the evolution of low-mass stars, brown dwarfs and giant planets is indeed entirely determined by their entropy profile (Chabrier & Baraffe 2000).

In the present paper, we follow the same method as mentioned above, by combining QMD calculations for pure hydrogen and helium with the SCvH (1995) and the Chabrier-Potekhin (1998) EOSs. As just mentioned, a striking advantage of the present calculations is that they provide the entropy over a wide temperature-pressure-density range, namely $10^{-8}$ to $10^{6}\,\,{\rm g}\,{\rm cm}^{-3}$, $10^{-9}$ to $10^{13}$ GPa and $10^{2}$ to $10^{8}$ K, covering essentially the domain of all dense astrophysical bodies. The paper is organized as follows. Sections 2 and 3 describe the H and He EOS, respectively, and make comparison with available experimental data or ab initio calculations. The calculations for the H/He mixture are described in §4, examples of the tables are presented in §5 while §6 is devoted to the conclusion.

2. The Hydrogen equation of state

2.1. Construction of the EOS model

Following the same procedure as Becker et al. (2014), our hydrogen EOS combines different calculations. For $T\geq 1.1\times 10^{5}$ K and/or $\rho>10\,{\rm g}\,{\rm cm}^{-3}$, the hydrogen becomes fully ionized and we use the Chabrier & Potekhin (1998, CP98) EOS model, based on the linear response theory to treat ion-electrons interactions. The CP98 EOS extends to very high temperatures or densities, when electrons become relativistic and recovers the Potekhin & Chabrier (2000) model EOS which handles the solid phase. The relativistic domain concerns essentially neutron stars or white dwarfs. For $T<1.1\times 10^{5}$ K, the EOS is divided into 3 density regimes, where we use 3 different EOS calculations:

• •

$\rho\leq 0.05\,\,{\rm g}\,{\rm cm}^{-3}$: SCvH EOS

• •

$0.3<\rho\leq 5.0\,\,{\rm g}\,{\rm cm}^{-3}$: EOS of Caillabet et al. (2011)

• •

$\rho>10.0\,\,{\rm g}\,{\rm cm}^{-3}$: CP98 EOS

Between these limits, a bicubic spline interpolation is performed which ensures continuity of the functions and their two first derivatives.

Whereas, as mentioned in §1, the calculations are performed over a vast density-temperature domain, namely $10^{-8}\leq\rho\leq 10^{6}\,\,{\rm g}\,{\rm cm}^{-3}$ and $10^{2}\leq T\leq 10^{8}$ K, several limits are identified.

(1) For T\,\hbox{\hbox{<}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,10^{3} K, hydrogen becomes solid over some pressure/density range. The melting line for H2 has been determined experimentally up to $T\simeq 1000$ K, $P\simeq 100$ GPa (Datchi et al. 2000, Deemyad & Silvera 2008, Eremets & Trojan 2009) and has been extrapolated up to about 300 GPa by the following fonctional form (Kechin 1995), based on QMD simulations (Bonev et al. 2004, Morales et al. 2010a, Caillabet et al. 2011222Note the typo in eqn. (41) of Caillabet et al. (2011), corrected in eqn.(1) here.):

[TABLE]

with $T_{0}=4.853$ K, $a=0.023$ GPa, $b=0.748$ and $c=0.0098$ GPa*-1* for H2. Note, however, that the turning point at P\,\hbox{\hbox{>}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,100 Mbar predicted by this form has not been confirmed unambiguously yet by experiments. The experimentally determined H2 melting curve is identified in Fig. 1 by the thick solid line and is continued by the dotted line at higher density/pressure according to the above functional form. At higher temperatures and pressures, hydrogen becomes fully dissociated and ionized, reaching the limit of the one component plasma model (OCP) whose melting line corresponds to $\Gamma_{m}=175$ (Potekhin & Chabrier 2000), where $\Gamma=(Ze)^{2}/(ak_{B}T)=2.25\times 10^{5}(Z^{2}/A)(\rho^{1/3}/T)$ is the usual plasma coupling parameter. This yields the following melting line:

[TABLE]

as identified in Fig. 1. Interestingly enough, extrapolating this line to lower temperatures and densities (long-dash) nicely joins the H2 melting line. Note, however, that this line is just indicated for illustrative purposes and can not be considered as a rigorously determined melting line at high pressure. Indeed, at very high densities/pressures, quantum diffraction effects between protons become significant and the classical OCP model becomes invalid (see below).

(2) At low temperature and high density, quantum (diffraction) effects between ions become important. In the CP98 model, these effects are treated within the $\hbar^{2}$ Wigner-Kirkwood expansion to second order. This yields the free energy quantum correction $f_{WK}=F_{WK}/Nk_{B}T=\eta^{2}/24$, where $\eta=\hbar\Omega_{p}/k_{B}T\approx 7.71\times 10^{3}\rho^{1/2}T^{-1}ZA^{-5/3}$ and $\Omega_{p}=(4\pi(Ze)^{2}n_{i}/M_{i})^{1/2}$ denotes the ion plasma frequency. For f_{WK}\,\hbox{\hbox{>}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,0.7, the CP98 model has been found to become of dubious validity and then the present EOS can not be used beyond this limit. This limit corresponds to

[TABLE]

indicated by the short-dash line in Fig. 1. Beyond this limit, the treatment of quantum effects requires fully quantum numerical calculations such as PIMC or PIMD. Such a quantum domain for hydrogen, however, does not concern any astrophysical body (see, e.g. Chabrier 1993).

The hydrogen QMD calculations in the intermediate density regime (see above) are based on Caillabet et al. (2011), which gather QMD simulations by Holst et al. (2008), coupled electron-ion Monte Carlo (CEIMC) calculations by Morales et al. (2010a), and PIMC calculations by Militzer & Ceperley (2000). These calculations have been supplemented by further QMD calculations for our present purpose. The excess free energy, $F_{ex}$, was fitted by a functional form similar to the one proposed in Chabrier & Potekhin (1998), which accurately recovers all appropriate limits. The accuracy of this analytical parameterization was verified by the fact that its temperature derivative properly recovers the excess internal energy, $U_{ex}$, obtained in the simulations. In the present calculations, however, we found out that, whereas the fit for $F_{ex}$ used in Caillabet et al. (2011) correctly recovers the H2 melting curve, it becomes less accurate away from these conditions. Therefore, in the present calculations, we have modified the $d(\rho)$ parameter of the fit given in eqn.(24) of Caillabet et al. (2011) in various density domains in order to recover the ab initio calculations of Morales et al. (2010a). The results will be illustrated in §2.3 below for H and in §4 for the H/He mixture.

As mentioned above, the EOS is calculated initially in a $T$-$\rho$ domain, appropriate to QMD or PIMC calculations, and then transformed into a $T$-$P$ one by bicubic interpolation procedures. In §2.3, we will make extensive comparisons between our results and available numerical results from ab-initio simulations for several thermodynamic quantities in order to verify the validity of these EOS calculations.

2.2. Comparison with experimental results

The validity of the EOS of hydrogen, or its isotope deuterium, can be first assessed by comparison with high-pressure Hugoniot experiments. As mentioned in the introduction, these latter include different techniques. The original discrepancies between these various data sets have been significantly reduced when using a revised EOS of quartz for the impedance-matching in the case of laser-compression experiments (Knudson & Desjarlais 2009), and all results now agree reasonably well to provide a robust compression Hugoniot curve up to about 200 GPa. The precision of these measurements has been improved recently with magnetically accelerated flyer plate experiments on deuterium, reaching a precision of $\sim$1.5%-$1.9\%$ in density along the Hugoniot, and carrying out reshock measurements from these Hugoniot states that provide off-Hugoniot data in the $\sim 100$-200 GPa and $\sim 5000$-15000 K regime (Knudson & Desjarlais 2017). Furthermore, experiments combining static and dynamic methods, generating laser-driven planar shock waves in precompressed samples of different initial densities have allowed to explore a larger domain off the principal Hugoniot, probing the EOS of hydrogen isotopes over an even larger pressure-temperature domain, directly probing the conditions in giant planet interiors (Loubeyre et al. 2012, Brygoo et al. 2015).

In order to compare our EOS calculations with experiments, we have calculated Hugoniot pressure-density and pressure-temperature curves for D2 and H2. Postshock conditions are calculated from mass, momentum and energy conservation across the shock by finding solutions from the tabulated EOS that satisfy the Hugoniot relation (Zel’dovich & Raizer 2002) for given initial conditions ($\rho_{0},P_{0},{\tilde{E}}_{0}$), where these quantities denote respectively the mass density, pressure and specific internal energy333Note that for the shock velocities under consideration, $U_{S}\approx 10$-50 km s*-1*, radiative effects in the energy balance are negligible.:

[TABLE]

Figure 2 compares the Hugoniot compression curve for deuterium obtained with our EOS with the most recent data of Knudson & Desjarlais (2017), which include also some of the aforementioned experimental results data for an initial state $\rho_{0}(D_{2})=\,0.167\,{\rm g}\,{\rm cm}^{-3}$ at $T_{0}=20$ K (the results originally obtained for a slightly larger initial density have been rescaled accordingly (see Knudson & Desjarlais (2017) for details)). To calculate the deuterium Hugoniot curve, we have rescaled the hydrogen EOS by a factor $2$ in density, but proper quantum corrections on the energy are taken into account. Since D2 and H2 have similar molar volumes at these conditions, the Hugoniot curves are nearly identical for these two elements. Our EOS is in excellent agreement with the data, as already noted in Caillabet et al. (2011), including with the most recent experiments. The maximum discrepancy occurs at $P=50$ GPa and amounts to $\sim 3\%$ on the density. In contrast, the SCvH EOS is less compressible in the low pressure domain and more compressible at higher pressures. Since the compression peak corresponds to the domain of molecular dissociation (energy goes into the breaking of internal levels and molecular bonds, yielding an increase of $\rho/\rho_{0}$), this behaviour reflects a well known shortcoming of this EOS, which underestimates H2 pressure dissociation. This stems essentially from the too stiff H-H and H2-H potentials used in the Saumon-Chabrier theory, which do not include the softening due to N-body interactions, in contrast to the case of the H2-H2 potential (see Saumon & Chabrier 1991). Indeed, high-pressure experiments at this time were not reaching high enough pressures to explore the dissociated regime, and thus could not provide experimental guidance to derive softened potentials for interacting atomic species. The inset in Fig. 2 displays the comparison between the present calculations and the PIMC simulations by Militzer & Ceperley (2000) at higher pressures.

Figure 3 compares our theoretical Hugoniots for H2 and D2 with the ones obtained for various pre-compressed initial states (Loubeyre et al. 2012, Brygoo et al. 2015). Initial states have pressures $P_{0}=0.16$ GPa, 0.3 GPa, 0.7 GPa and 1.5 GPa at 297 K, respectively. Again, the agreement between the data and the present EOS is excellent for all series of experiments. Also shown for comparion is a predictive Hugoniot calculated for an initial pressure $P_{0}=$6.0 GPa, as planed with future high-pressure experiments, as well as a typical Jupiter internal isentrope for a helium number fraction $x_{He}\simeq 0.08$ (mass fraction $Y\simeq 0.25$ (see §4)). As seen in the figure, conditions along this Hugoniot are very close to or intercept Jupiter’s internal density and temperature profiles (assuming an isentropic thermal structure), respectively, notably in the crucial $\sim$Mbar pressure ionization region, and thus directly probe Jupiter’s deep interior. The same aforementioned general behaviour of the SCvH EOS, i.e. underestimated molecular hydrogen pressure dissociation, is observed for all Hugoniots and is particlarly striking along the 6.0 GPa one.

Figure 4 portrays the temperature-pressure curves along the Hugoniots for the same sets of experiments. Interestingly enough, the difference between the present EOS, which includes ab initio simulations, and the semi-analytic SCvH EOS is much smaller than for the $P$-$\rho$ compression curves. We note, however, that molecular dissociation in the SCvH EOS not only occurs at too high pressures, as mentioned above, but takes place very abruptly, as shown by the kinks in the dotted curves, yielding cooler temperatures at given pressure as energy goes into molecular dissociation instead of raising $k_{B}T$. Experiments (Loubeyre et al. 2012), in contrast, have revealed that reflectivity, then electrical conduction, increases gradually along the Hugoniot above about $\sim$ 5000 K, before reaching a plateau, reflecting the dissociation and ionisation of molecular hydrogen H2 and suggesting that this process, under the conditions probed by present Hugoniot experiments, occurs continuously. This is in agreement with ab initio calculations, which show an increasing conductivity along the Hugoniot, but predict a discontinuous molecular-ionic transition around P$\simeq 100$ GPa at cooler temperatures, in the range $T\simeq 2000$-6000 K (Morales et al. 2010b, Mazzola et al. 2018).

Another important experimental constraint on the EOS comes from quasi-isentropic ramp compression of hydrogen or deuterium. Those experiments have direct astrophysical applications since, as will be discussed later, the interiors of low-mass objects are entirely convective such that their internal temperature profile follows an isentrope. Dynamic quasi-isentropic shock wave experiments using high explosives on deuterium have been carried out up to about 1500 GPa (15 Mbar) and densities about 4.5 $\,{\rm g}\,{\rm cm}^{-3}$, directly probing jovian planet deep interiors. While, in some cases, both the density $\rho$ and the pressure $P(\rho)$ were measured simultaneously (Boriskov et al. 2011), in other experiments only the densities were measured while the pressure was determined afterwards from a hydrocode with a model EOS (Fortov et al. 2007). An extension of these latter experiments was carried out by Mochalov et al. (2010) up $\rho=108\times\rho_{0}=4.3\,{\rm g}\,{\rm cm}^{-3}$, reaching an unprecedented experimental pressure for D2 of 1800 GPa (18 Mbar). In all cases, the temperatures were determined from a model EOS. As noted by Becker et al. (2013), however, the experimental points of Fortov et al. and Mochalov et al. were found not to lie on the same isentrope, questioning the validity of the results, at least of the model-dependent pressure determination from the measured density. Following Becker et al. (2013), we have calculated the isentropic compression path obtained with our and SCvH EOS’s, respectively, starting from Fortov et al. (2013) model-independent initial condition for D2, $\rho_{0}=0.04\,{\rm g}\,{\rm cm}^{-3}$, $T_{0}=283$ K. According to the present and SCvH EOSs, this corresponds to an entropy $S=9.9\,k_{B}$/atom. The result is displayed in Fig. 5. The empty symbols are the quoted experimental $\rho$-$P$ determinations while the solid symbols correspond to the values obtained by Becker et al. (2013, Table I) with their EOS for the above initial conditions. As seen in the figure, our isentrope agrees very well with the values obtained by Becker et al., with then all the experimental results lying on the same isentrope. We see in particular that the pressure rises continuously with the density along the isentrope, with no sign of discontinuity due to a first-order phase transition in this regime, as suggested by Fortov et al. (2007). Interestingly enough, we see that the SCvH EOS predicts larger pressures for a given density (by about $\sim 10\%$ at 2 $\,{\rm g}\,{\rm cm}^{-3}$), i.e. a significantly ($>20\%$) warmer isentrope.

2.3. Comparison with ab initio calculations

Figure 6 compares the specific internal energy, $\tilde{U}$, as a function of density for the present calculations with available PIMC (Militzer & Ceperley 2001 (MC01), Hu et al. 2011) and QMD (Soubiran & Militzer 2015, SM) simulations over the available temperature-density ranges. We also make comparisons with the QMD simulations of Becker et al. (2014) over a larger density domain. Note that some of these simulations encompass the domain of hydrogen dissociation and ionization, i.e. $\rho\sim 0.5$-$5\,{\rm g}\,{\rm cm}^{-3}$, $T\sim 3000$-50000 K. All results are rescaled to the zero of energy of the present EOS, which is the same as in SCvH, namely the ground state of the H2 molecule. In all the domain explored by MC01 PIMC simulations, we note the good agreement between all different calculations, including SCvH and these simulations. Clearly, these simulations do not probe a density regime where differences between the various EOS due to the treatment of hydrogen dissociation and ionization arise. The Hu et al. (2011) PIMC and Soubiran & Militzer (2015) simulations, in contrast, reach higher densities and enter the crucial dissociation/ionization regime. The agreement between the present calculations and these simulations is excellent. We notice, however, the surprising behaviour of the Hu et al. (2011) calculations at high density for the 15 kK, 30 kK and 62 kK isotherms. We need to stress here that, for these ($T,\rho$) conditions, the temperature is of the order of the electron Fermi temperature ($\theta=T/T_{F}\simeq 1$). Under such conditions, Monte Carlo samplings are known to be extremely inefficient and can lead to unreliable results. We also notice an energy shift (i) between the PIMC and Becker et al. (2014) energy on one hand and the present and SCvH ones on the other hand for the 30000 K isotherm for $\rho\lesssim 0.1\,{\rm g}\,{\rm cm}^{-3}$ and (ii) between Becker et al. and the present or SCvH calculations for the coolest isotherms even at very low densities. The shift at $T=30000$ K most likely stems from the underestimated H2 dissociation in the SCvH EOS and thus in the present one below $\sim 0.3\,{\rm g}\,{\rm cm}^{-3}$, due to the interpolation procedure (see §2.1). Recall that our QMD calculations only extend down to 0.2 $\,{\rm g}\,{\rm cm}^{-3}$. The maximum discrepancy, however, is about 15% around $\sim 0.05\,{\rm g}\,{\rm cm}^{-3}$ and becomes negligible below $\sim 0.01\,{\rm g}\,{\rm cm}^{-3}$. For this temperature, the PIMC simulations predict 57%/43% H+/H ionization fractions, with $x_{H_{2}}=0$ at $\rho=2.7\times 10^{-3}\,{\rm g}\,{\rm cm}^{-3}$ (see Table I of MC01) whereas SCvH predict 54%/46% H+/H, with $x_{H_{2}}=0$, quite a good agreement. The shift for the coolest isotherms, notably at very low density, between Becker et al. (2014) and the SCvH and present calculations is more surprising as at these densities thermal dissociation and ionization, when they take place, are well described by the Saha equation, a limit correctly recovered by SCvH. For $T\leq 2000$ K, for which H2 rotational levels are excited, but not vibrational ones, the SCvH correctly recovers the perfect gas limit, ${\tilde{U}}=(5/2)\times k_{B}T/\mu=1.03\times 10^{-2}\times T$ MJ kg*-1*, where $\mu=A_{H_{2}}\times m_{H}$, with $m_{H}=1.660\times 10^{-27}$ kg the atomic unit mass, which does not seem to be the case of the Becker et al. EOS. Note in passing that the spin dependence of the H2 molecule (ortho- and para-hydrogen) is correctly accounted for in the Saumon-Chabrier theory (see Saumon & Chabrier 1991).

All curves exhibit a sharp rise above $\sim 0.1$-5$\,{\rm g}\,{\rm cm}^{-3}$, depending on the temperature. This corresponds to the increasing (repulsive) interactions between hydrogen molecules and/or atoms and then to pressure ionization and the onset of the electron degeneracy contribution. Not surprisingly, then, the most noticeable differences between the present, SCvH and Becker et al. (2014) EOSs occur in this domain, a regime covered by QMD simulations in both the present and Becker et al. EOSs but described by a semi-analytical model in SCvH. Generally speaking, the SCvH EOS overestimates the internal energy compared with the simulations in this domain, except for the lowest temperatures where it first underestimates and then overestimates it. This reflects the now well identified shortcoming of the Saumon-Chabrier (SC) theory that overestimates the density domain of stability of molecular hydrogen, and then predicts a too abrupt ionization, globally underestimating the domain of hydrogen dissociation/ionisation. This behaviour was already noted on the Hugoniot experiments (§2.2), and arises essentially on one hand from the too stiff interatomic potential compared to the intermolecular one in the SC theory, as mentioned previously, and also from the fact that atomic and molecular ionization, i.e. the Stark effect, are underestimated in the theory (see Saumon & Chabrier 1991, 1992).

Figure 7 displays similar comparisons for the energy per atom (as in MC01) as a function of temperature between PIMC (Militzer & Ceperley 2000, 2001) and the present and SCvH EOSs along the available isochores. As mentioned above, within this density range, all types of calculations are in good agreement. We note a slight departure between SCvH and the simulations for the highest density in the temperature range $\sim 30000$-80000 K, as already noted by Militzer & Ceperley (2001). This again reflects the imperfect treatment of ionization in the SC model. The deviations, however, remain modest, within at most $\sim 6\%$. In contrast, the present calculations agree very well with the MC01 PIMC results.

Figure 8 portrays now similar comparisons with another set of first principle simulations, namely the Coupled-Electron-Ion Monte-Carlo (CEIMC) calculations of Morales et al. (2010a), probing the very pressure-dissociation/ionization regime, between $\rho=0.7\,\,{\rm g}\,{\rm cm}^{-3}$ and 2.4 $\,{\rm g}\,{\rm cm}^{-3}$. Remember that in this regime, both the present and Becker et al. (2014) EOS calculations rely on QMD simulations. We first note the excellent agreement between QMD-based and CEIMC results in the probed temperature-density range, which gives confidence that both types of methods can handle this crucial density regime. We also note the strong departure between the SCvH EOS and these results, by as much as 25%, for the same reasons as mentioned previously.

Figures 9 and 10 compare the pressure, more precisely $P/\rho$ to highlight non-ideal contributions, as a function of density of the present EOS with the same set of simulations. As seen in Fig. 9, we first note that, as for the internal energy, all EOS calculations agree very well with the simulations below about 0.1 $\,{\rm g}\,{\rm cm}^{-3}$. Above this value, the SCvH EOS significantly overestimates the pressure for T\,\hbox{\hbox{<}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,15 kK and underestimates it at larger temperatures, by as much as 25% around 1 $\,{\rm g}\,{\rm cm}^{-3}$, highlighting again the approximate treatment of pressure dissociation and ionization in the SC theory. In contrast, the present EOS is in excellent agreement with the Soubiran & Militzer (2015) simulations, for the available isotherms, and with the Hu et al. (2011) ones for $T>60000$ K. For these latter, we note the same spurious behaviour for $T=15000$ K and 30000 K, as for the energy, which confirms the dubious validity of these results in the partially degenerate domain ($T/T_{F}\sim 1$). Figure 10 confirms the excellent agreement between the present and Becker et al. (2014) EOS’s, as well as with the CEIMC (Hu et al. 2011) simulations at higher density, and the previously identified shortcomings of the SCvH EOS in this regime. For sake of completeness, we have carried out similar comparisons for the pressure with the calculations of Militzer & Ceperley (2001) for the various isochores given in these simulations. This is portrayed in Figure 11. As seen in Fig. 10, we note some small wiggles in both the present and Becker et al. EOS’s for the $T=1000$ K and 2000 K isotherms (bottom curves), which stem from the imperfect interpolation procedures in their construction.

As mentioned in the Introduction, a major improvement of the present EOS over previous calculations is that it provides the entropy. Indeed, Caillabet et al. (2011) derived a parameterization of the free energy $F$, yielding the entropy as $S=(F-U)/T$. However, as mentioned previously, the analytical fit of Caillabet et al. is valid over a rather limited temperature-density range, close to the hydrogen melting curve. In order to extend the validity of the free energy and entropy over a larger range, we have corrected the fit in various $T$-$\rho$ places to recover the results of Morales et al. (2010a) for H and Militzer & Hubbard (2013) for the H/He mixture (see §4). The comparisons between the present calculations and the fitting parameterization derived by Morales et al. (2010a) from their simulations for the free energy and the entropy are illustrated in Figures 12 and LABEL:fig13, respectively. Data points from QMD simulations by Soubiran & Militzer (2015) for the entropy are also shown in Fig. LABEL:fig13. The agreement between the present EOS and the Morales et al. results for $F$ is excellent, in contrast to the SCvH results. For the entropy, although not perfect, the agreement between the present calculations and the results fitted from the simulations is also much better than for the SCvH EOS, in particular for the coolest isotherms. The sudden rise of entropy above \,\hbox{\hbox{>}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,0.3\,{\rm g}\,{\rm cm}^{-3} reflects the pressure dissociation and ionization, leading notably to an increase in the number of particles. The most noticeable differences appear for the $T=6000$ K isotherm between the present (or SCvH) results and the Soubiran & Militzer (2015) ones in the range $\sim 0.3$-$0.6\,{\rm g}\,{\rm cm}^{-3}$. The discrepancy, however, remains modest, at less than $5\%$. For the coolest isotherms, $T\leq 3000$ K, the difference between the present calculations and the MC01 or SM ones rather stems from the onset of ion quantum effects (see Fig. 1), which are included throughout the $\hbar^{2}$ Wigner-Kirkwood correction in the present calculations but are not taken into account in the simulations.

Besides the specific internal energy $\tilde{U}$ and entropy $\tilde{S}$, the present EOS delivers all the necessary thermodynamic quantities. These include the specific heats at constant volume and pressure, $\tilde{C_{V}}$, $\tilde{C_{P}}$, from the relations:

[TABLE]

where $\chi_{T}=(\frac{\partial\log P}{\partial\log T})_{\rho}=-(\frac{\partial\log\rho}{\partial\log T})_{P}/(\frac{\partial\log\rho}{\partial\log P})_{T}$ and $\chi_{\rho}=(\frac{\partial\log P}{\partial\log\rho})_{T}$. Figure 14 compares these quantities for the present calculations and the SCvH ones for $T=$15000, 20000 and 30000 K. In the low density limit, we recover the value for molecular hydrogen with two rotational and two vibrational degrees of freedom (the rotation and vibration temperatures of H2 are $\theta_{rot}=85$ K and $\theta_{vib}=6120$ K, respectively), i.e. ${\tilde{C}}_{V}=(\frac{7}{2})\times 8.25\times 10^{-3}$ MJ kg*-1* K*-1* = 0.029 MJ kg*-1* K*-1*, and ${\tilde{C}}_{P}=\frac{9}{7}\times{\tilde{C}}_{V}$ = 0.038 MJ kg*-1* K*-1*. As molecular dissociation or ionization take place, the increase of the number of particles and the release of dissociation or ionization energy yield an increase of the specific heats, with maxima corresponding to the partial dissociation or ionization zones. Eventually, at high density, the specific heats decrease to reach the limit of a monatomic (ionized) gas, i.e. ${\tilde{C}}_{V}=\frac{3}{2}\times 8.25\times 10^{-3}$ MJ kg*-1* K*-1* = 0.012 MJ kg*-1* K*-1*, ${\tilde{C}}_{P}=\frac{5}{3}{\tilde{C}}_{V}$ for an ideal gas, potentially corrected by non-ideal contributions. As seen in the figure, these isotherms bracket the domain of hydrogen pressure dissociation and ionization which occurs between about $\sim 0.1$ and 3.0 $\,{\rm g}\,{\rm cm}^{-3}$. The inset clearly highlights the lack of H2 ionization at high density in the SCvH model.

Figure 15 displays the isothermal compressibility factor

[TABLE]

where $\kappa_{T}^{0}=m_{H}/(\rho k_{B}T)$ is the isothermal compressibility of a perfect monatomic hydrogen gas and $\mu=Am_{H}$ is the atomic weight, with $m_{H}=1.660\times 10^{-27}$ kg the atomic unit mass, for the present and SCvH calculations over the temperature and density domain characteristic of hydrogen pressure and thermal dissociation and ionization. The figure highlights the too large compressibility of the SCvH EOS in this domain, as noticed in §2.2 over the Hugoniot compression curves, due to the lack of H2 dissociation/ionization.

3. The Helium equation of state

3.1. Construction of the EOS model

The procedure for the helium EOS is exactly the same as for the hydrogen one, with the combination of different calculations. For $T\geq 1.0\times 10^{6}$ K, the plasma becomes fully ionized and we use the Chabrier & Potekhin (1998) EOS model. For $T<1.0\times 10^{6}$ K, the EOS is divided again into 3 density regimes. In the low density (atomic) one, we use the SCvH EOS for pure He. In the intermediate $T$-$\rho$ regime, the ab-initio calculations are based on original QMD simulations and will be described in details in a dedicated paper (Soubiran et al. 2019, in prep.). In the high density, fully ionized domain, we use the Chabrier-Potekhin (1998) EOS.

• •

$\rho\leq 0.1\,\,{\rm g}\,{\rm cm}^{-3}$: SCvH EOS

• •

$1.0<\rho\leq 100.0\,\,{\rm g}\,{\rm cm}^{-3}$: EOS of Soubiran et al. (2019), based on QMD calculations

• •

$\rho>100.0\,\,{\rm g}\,{\rm cm}^{-3}$: CP98 EOS

As for hydrogen, bicubic spline procedures are used to interpolate the thermodynamic quantities in the intermediate regime. The fact that we merge the QMD and CP98 calculations at $100\,{\rm g}\,{\rm cm}^{-3}$ is justified by the fact that QMD calculations (Soubiran et al. 2012), reanalyzing reflectivity measurements of dense fluid helium (Celliers et al. 2010) by including the effects of temperature on the helium gap, suggest that this latter closes at a density of about 10$\,{\rm g}\,{\rm cm}^{-3}$, in good agreement with previous semi-analytical models (Winisdoerffer & Chabrier 1995), implying that helium should be fully ionized above this density. The zero of energy for the helium EOS is the same as in SCvH, namely the zero of the isolated He atom.

As for hydrogen, although, for practical purposes, the tables are calculated over square $T$-$\rho$ and $T$-$P$ domains, part of these latter are meaningless, as they correspond to regions where either helium becomes solid (Loubeyre et al. 1993) or quantum diffraction effects for ions become dominant. The melting line for helium, determined by diamond anvil cell experiments, is well described by a simple Simon law (Datchi et al. 2000), even when extrapolated to Mbar pressures, where this expression is in good agreement with ab initio calculations (Lorenzen et al. 2009):

[TABLE]

where the pressure $P$ is in kbar ($=0.1$ GPa). This is indicated by the thick solid line in Fig. 16, while the OCP melting line (eqn.(2)) is shown by the long-dashed line. Interestingly enough, as for hydrogen, we see that this curve, when extended to low temperatures and densities, nicely merges with the experimental one. As for the onset of quantum diffraction effects, the characteristic parameter is the same as for hydrogen, f_{WK}\,\hbox{\hbox{>}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,0.7, and the condition given by eqn.(3) is indicated by the short-dashed line.

3.2. Comparison with ab initio calculations

Extensive comparisons between this pure He EOS and PIMC or existing QMD simulations will be presented in details in a forthcoming paper (Soubiran et al., 2019, in prep.). Meanwhile, Figures 17 and 18 compare the Hugoniot compression curves obtained with our He EOS with the recent data of Brygoo et al. (2015), for different precompressed initial conditions. The agreement between the present calculations and the data is very good, except for the two data points at lowest pressure of the $P_{0}=0.3$ and 0.5 GPa precompressed experiments, which seem to be surpringly stiff. As for hydrogen, we note that the SCvH model predicts a too abrupt ionization compared with both the experiments and the present calculations, which rather suggest a smoothly ongoing process. This may again point to a limitation of the so-called chemical semi-analytical model, based on the concept of pair potentials for the various species (at present He, He+ and He*++*), to describe N-body interactions, as the same abrupt ionization is found in the more sophisticated model of Winisdoerffer & Chabrier (2005) for helium.

4. The Hydrogen/Helium mixture equation of state

4.1. Calculation of the H/He EOS

The calculations of the EOS for the H/He mixture are carried out within the so-called ”additive volume law” (AVL), as in SCvH. This latter is based on the additivity of the extensive variables (volume, energy, entropy,…) at constant intensive variables (P,T). Although this method is formally exact for non-interacting, ideal mixtures, and excellent in the limit of fully ionized systems (Chabrier & Ashcroft 1990), it is no longer the case for interacting systems, i.e. between hydrogen and helium species in the present context, or in the domain of partial ionization. Nevertheless, we expect the correction to remain modest, of the order of a few per cent. Clearly, this is a limitation of the present EOS. Calculations for the interacting H/He mixture have been carried out recently with QMD simulations but have focused on a limited density-temperature domain characteristic of Jupiter internal adiabat (Morales et al. 2010b, Lorenzen et al. 2009, 2011, Militzer 2013, Mazzola et al. 2018, Schöttler & Redmer 2018) and only one of these calculations has calculated the entropy of the interacting mixture (Militzer & Hubbard 2013). Note also that all these simulations have been carried out with a rather small number of particles, of the order of 10 for the helium atoms. Indeed, for large numbers of particles, demixing can occur in the simulation box, preventing the calculation of the thermodynamic properties of the mixture (Lorenzen et al. 2009, Soubiran et al. 2013). Therefore, even though various schemes exist to correct for the finite size errors, it seems fair to say that the quantities derived from existing simulations still retain some degree of uncertainty.

As mentioned above, within the AVL, an extensive variable $W$ at given $(T,P)$ for the mixture reads (see, e.g. SCvH):

[TABLE]

where $X_{i}=M_{i}/(\sum_{i}M_{i})$ denotes the mass fraction of component $i$ (H or He in the present context). The density for the H/He mixture (which is an inverse specific volume) reads

[TABLE]

where $Y=M_{He}/(M_{H}+M_{He})$ denotes the helium mass fraction. For the specific entropy of the mixture, the ideal mixing entropy must be added to eqn.(8) in order to correctly recover the ideal gas limit, yielding

[TABLE]

For a mixture of $N=\sum_{i}N_{i}$ components $i$ of number fraction $x_{i}=N_{i}/N$ and atomic mass $A_{i}m_{H}$, with $m_{H}$ the atomic mass unit, the ideal mixing entropy reads

[TABLE]

where $\langle A\rangle=\sum_{i}x_{i}A_{i}$ and $k_{B}$ denotes the Boltzmann constant. It should be noted that in the above equation, we have neglected the contribution from the free electron entropy. Indeed, in contrast to semi-analytical so-called ”chemical models” such as SCvH, based on well defined chemical entities such as molecules, atoms and electrons, such an identification does not exist in QMD simulations, preventing the precise caracterization of a free electron density. Our approximation, however, is justified both in the regime of neutral hydrogen and helium, where there are no free electrons, and in the regime of full ionization, where the electrons become degenerate and thus have a negligible entropy. The approximation, however, fails in between these two regimes, i.e. in the regime of partial ionization.

As mentioned above, Militzer (2013) and Militzer & Hubbard (2013) have carried out QMD simulations for a H/He mixture over a significant temperature-density domain and have calculated the free energy $F$ by thermodynamic integration, which yields also the entropy $S=(U-F)/T$. These authors provide a fitting formula for $F$ over the range of their simulations, namely $\sim 0.2$-9.0$\,{\rm g}\,{\rm cm}^{-3}$ and 1000-80000 K, covering the domain of H and He pressure and temperature dissociation and ionization. The simulations were carried out for a helium number fraction $x_{He}=18/238=0.076$, corresponding to a helium mass fraction $Y=0.246$. Figures 19 and 20 display the comparison of the internal energy $E$ and the excess pressure $P/\rho$ as a function of density for several isotherms calculated by MH13 with the present and SCvH results. As noted in MH13, the SCvH EOS generally slightly overestimates the internal energy compared with the simulations over the probed density range. This is improved with the new EOS, although this latter predicts a lower internal energy than MH13 for the highest isotherms, most likely due to the temperature and density interpolation procedures between the QMD and CP98 calculations in this domain. We also have to recall that non-ideal H/He mixing effects are not included in the present calculations, based on the additive volume law. Indeed, it has been shown that in the ionized regime, this latter yields a lower energy than for the non-ideal mixture (Chabrier & Ashcroft 1990). The sharp increase in the energy in the density regime $\sim 0.5$-1.0 $\,{\rm g}\,{\rm cm}^{-3}$ in the displayed temperature regime stems from the onset of ionization, yielding a strong increase of the free electron energy contribution due to the Pauli principle (see e.g. Saumon & Chabrier 1992, Militzer & Hubbard 2013).

Figures 21 and 22 portray similar comparisons as a function of temperature for several isochores calculated by MH13. In Figure 21, we notice the abrupt increase of the SCvH energy w.r.t. both MH13 and the present calculations around $T\sim 30000$ K and $\rho\gtrsim 0.2\,{\rm g}\,{\rm cm}^{-3}$, i.e. in the regime of pressure dissociation/ionization, while the reverse is true for the pressure, with a cross over of the SCvH isochores with the present and MH13 ones in the regime $\sim$0.2-3.0 $\,{\rm g}\,{\rm cm}^{-3}$ around $T\sim 30000$ K ($\log T=4.5$) (see Fig. 22). This again illustrates the approximate treatment of this process in the semi-analytical SC model, as mentioned in the previous sections and as already noted by Militzer & Ceperley (2001) for pure hydrogen and Militzer & Hubbard (2013) for the mixture, highlighting the already mentioned lack of dissociation and too abrupt ionization in SCvH with increasing density. The underestimated degree of molecular dissociation and/or ionization in the SC model is also reflected by the increasing offset between SCvH and both MH13 and the present calculations for $T\lesssim 10^{4}$ K for both $U$ and $P/\rho$ in the density regime $\sim 0.75$ - 3.6$\,{\rm g}\,{\rm cm}^{-3}$. As pointed out by MH13, this discrepancy on the pressure as a function of temperature can have a significant impact on giant planet internal structures. In contrast, the agreement between the present EOS and MH13 in this crucial domain can be considered as satisfactory. At higher densities, when the system becomes dominantly ionized, all calculations agree. Generally speaking, the present EOS agrees well with the MH13 simulations, except possibly in the domain 4.5\,\hbox{\hbox{<}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,{\rm log}\,T\,\hbox{\hbox{<}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,5.0 for for $\rho\simeq 2.0$-6.0$\,{\rm g}\,{\rm cm}^{-3}$, as seen in the figures, with a maximum discrepancy of $\sim 8\%$ . Since this is within the domain of interpolation between the QMD-based simulations and the CP98 model in the present EOS calculations (see §2.1), the discrepancy is likely to be blamed upon this procedure.

As already mentioned, besides the pressure and the internal energy, the knowledge of the entropy is necessary to determine the thermal profile and the cooling rate of objects below about 0.6 ${\rm M}_{\odot}$. This domain encompasses low-mass stars, brown dwarfs and gaseous planets. Indeed, these objects are too cool for heat to be transported efficiently by radiation and not dense enough for electron conduction to be significant. Heat is thus transported by convection, yielding a nearly adiabatic internal profile. Deriving the entropy over a large enough temperature-density range to cover the evolution of these bodies is thus of prime importance for astrophysical applications as well as for isentropic high-pressure experiments aimed at characterising hydrogen and helium pressure ionization. So far, no such EOS has been derived.

Militzer (2013) and Militzer & Hubbard (2013) calculated the Helmholtz free energy from their QMD simulation data by performing a so-called thermodynamic integration technique (TDI). The advantage of this technique, where integration is performed over trajectories that are derived from a hybrid potential energy function between the one of a classical system and the one obtained with a Kohn-Sham functional, is that it does not require a prohibitively large number of $(T,\rho)$ simulation points. The other advantage of the TDI method is that it allows directly the determination of the ionic contributions to the entropy. Whereas in most cases this contribution represents essentially a measure of the total entropy of the system, this is no longer true when electronic excitations become important, i.e. once ionization takes place. In that case the electron contribution to the entropy must be taken into account in the TDI integration (see Militzer (2013) and references therein). It must also be kept in mind that the procedure becomes less straightforward in the molecular regime, where a rigorous classical reference system is more difficult to define, because of exchange reactions, leading to dissociation and recombination. Last but not least, in some (low)temperature-(high)density domain, corrections due to quantum effects in the motion of nuclei must also be taken into account in the DFT-MD results. Finally, finite-size effects due to the limited number of particles must be treated with extreme care to ensure they do not affect substantially the results. Computational calculation of the entropy of a system is thus a highly delicate task and is not free from uncertainties.

Figures 23 and 24 portray a comparison of the free energy $F$ per atom as a function of temperature and density, respectively, between the present calculations, the SCvH EOS and the MH13 simulations in the density-temperature range probed by these latter, using either their numerical data points or their polynomial fit within its domain of validity. As already noticed by MH13, the agreement for this quantity is much better than for the pressure and the internal energy, which are respectively the density and temperature derivatives of $F$. We note, however, the better agreement of the present calculations with the simulations compared with SCvH in the $T$-$\rho$ domain where ionization sets in.

Figures 25 and 26 show the same comparisons for the entropy. For the coolest isotherms ($T<5000$ K) and low densities ($\rho\lesssim 0.3\,{\rm g}\,{\rm cm}^{-3}$), i.e. in the molecular/atomic domain, all calculations agree quite well, showing that the SC model adequately handles this regime, even when interactions between H2 molecules or He atoms become significant. For higher temperatures and densities, the SC model starts to depart from both the MH13 and present results, first underestimating the entropy in the domain $0.2\lesssim\rho\lesssim 2\,{\rm g}\,{\rm cm}^{-3}$ and $5000\lesssim T\lesssim$ 10000 K and then showing an abrupt increase of the entropy in this density regime at higher temperatures. This corresponds to the very domain of pressure ionization and reflects the already mentioned inaccurate (and too abrupt) treatment of this process in the SC theory. In contrast, the agreement between the present calculations and the MH13 results can be considered as very satisfactory over the entire temperature-density range explored by the simulations. The sudden decrease in entropy for $T=5000$ K and $\rho\gtrsim 2\,{\rm g}\,{\rm cm}^{-3}$ in the fit derived from the MH13 simulations compared with both the present and SCvH results, as seen in both Figures 25 and 26, is rather surprising and might point to either the increasing contribution to the interactions between H and He species or an issue with the TDI procedure or the inferred fitting formula. Note, however, that quantum effects between ions become significant in this regime (see Fig. 1) and that either the present calculations treat them as a perturbation, with the Wigner-Kirkwood expansion (SCvH and present) or they are ignored (MH13). As seen in the figures, however, this region concerns a high density and low entropy ($S<6\,k_{b}/$e-) domain where there are no astrophysical objects.

As mentioned above, interiors of astrophysical bodies below $\sim 0.6\,{\rm M}_{\odot}$ are essentially convective and thus nearly adiabatic. Their internal profile is thus characterized by an isentrope for a given mass at a given age and their thermal evolution corresponds to a series of decreasing isentropes. Figure 27 portrays the temperature and pressure profiles of such isentropic structures for the present, SCvH and MH13 calculations for entropy values between 4 and 16 $k_{b}/$e-, i.e. $2.9\times 10^{-2}$ to $1.15\times 10^{-1}$ MJ kg*-1* K*-1* for the present H/He mixture ($Y=0.246$, i.e. $x_{He}=0.076$). For adiabats above $S\gtrsim 10\,k_{b}/$e-, we note that the SCvH adiabats underestimate both the temperature and the pressure in the density domain $\sim 0.1$-2.0 $\,{\rm g}\,{\rm cm}^{-3}$, i.e. the domain of ionization, compared with the two other types of calculations, both based on QMD calculations in this regime, which predict higher temperatures and pressures in this domain. Again, this reflects the inaccurate treatment of the ionization process (Stark effect and electronic excitations) in the SC model.

Interestingly enough, the behaviur reverses for cooler isentropes, with SCvH predicting higher temperatures and pressures than the two other models. In this regime, molecular hydrogen H2 is still present and the disagreement arises from the lack of a proper treatment of H2 pressure dissociation in the SCvH model. As seen in the figure, the present EOS agrees fairly well with the calculations of MH13, notably for the pressure. For the temperature, the present calculations start departing from MH13 for $S\lesssim 10\,k_{b}/$e-, predicting slightly warmer structures than these latter calculations in the pressure ionization regime. This reflects the increasing contribution of the H/He interactions, and thus of the non-ideal mixing entropy, in the mixture, not treated in the present calculations, for the cooler and denser domains, yielding eventually a H/He phase separation for the coolest isentropes (Lorenzen et al. 2009, Morales et al. 2010b, Militzer 2013, Mazzola et al. 2018, Schöttler & Redmer 2018).

Since we do not have the results from MH13 for the exact Jupiter isentrope ($T=166$ K, $P=1$ bar) for the cosmogonic helium mass fraction ($Y=0.275$), we can not make a comparison for the correct Jupiter adiabat. The $S=7\,k_{b}/$e- one, however, is close enough to Jupiter’s value to estimate the discrepancy between the various calculations under Jupiter-like conditions. As can be inferred from Fig. 25 and 26, we found out that, for this entropy value, the maximum discrepancy occurs at $P=100$ GPa, with MH13 giving a temperature of $T=4867$ K, the present EOS $T=5382$ K, i.e. a 10$\%$ difference, and SCvH $T=5623$ K, i.e. 15.5$\%$ difference. Conversely, for the corresponding density $\rho\simeq 0.75\,{\rm g}\,{\rm cm}^{-3}$, MH13 predicts a pressure $P=94$ GPa against 93 GPa for the present calculations ($\approx 1\%$ discrepancy) and 132 GPa for SCvH ($=37\%$ discrepancy). This is the very domain of pressure ionization, so these differences are not surprising and illustrate the better treatment of this process in the present calculations, based on QMD simulations in this domain, compared with SCvH. The remaining discrepancy with MH13 can thus have two origins. The first one is errors in the parameterization of the free energy in the present calculations. The second one is of course the missing treatment of H/He interactions in the present EOS calculations and thus the lack of non-ideal mixing entropy. Indeed, MH13 numerical simulations reveal an H/He phase separation in this regime (see Fig. 2 of Militzer 2013), even though other simulations reach a different conclusion for similar $T$-$P$ values (Schöttler & Redmer 2018), suggesting that H and He are still miscible under Jupiter internal adiabat conditions. Looking at figures 25, 26 and 27, we can infer the impact of the aforementioned discrepancies in $T$ and $P$ between the present, SCvH and MH13 calculations in the pressure ionization domain under Jupiter-like conditions. We found out that for $T=5000$ K and $\rho=0.75\,{\rm g}\,{\rm cm}^{-3}$ (and $Y=0.246$), we get $S=6.9$ and 6.63 $k_{B}$/e- for the present and SCvH EOS, respectively, against 7.07 $k_{B}$/e- for MH13, i.e. $\sim$2% and 6% differences, respectively. A precise quantification of these differences upon Jupiter’s internal properties requires deeper explorations, with exact models of Jupiter, to be conducted in forthcoming calculations.

Note that, in the above comparisons, we have not included the EOS recently derived by Miguel et al. (2016+erratum). Indeed, the entropy values given by their tables for various $T$ and $\rho$ conditions differ significantly from all the calculations displayed in this section. This points to a severe issue in these tables (see also Debras & Chabrier 2019, §2.2).

4.2. Thermodynamic quantities

As mentioned earlier, the present EOS delivers all the necessary thermodynamic quantities besides temperature $T$, pressure $P$, specific internal energy $\tilde{U}$ and specific entropy $\tilde{S}$. These include the specific heats at constant volume and pressure, $\tilde{C_{V}}$, $\tilde{C_{P}}$ (see eqn.(5)), and the adiabatic gradient, $\nabla_{ad}$, from the relation

[TABLE]

where $\chi_{T}=(\frac{\partial\log P}{\partial\log T})_{\rho}=-(\frac{\partial\log\rho}{\partial\log T})_{P}/(\frac{\partial\log\rho}{\partial\log P})_{T}$ and $\chi_{\rho}=(\frac{\partial\log P}{\partial\log\rho})_{T}$. Note that since, in our calculations, the ideal mixing entropy, $S_{mix}$, does not depend on $T$ or $P$ (see eq.(11)), the adiabatic gradient (as well as the other first derivative quantities) at given $(P,T$) for the mixture can easily be calculated from the linear interpolation

[TABLE]

which may happen to be less numerically noisy than calculating the second derivative of a spline. Other quantities include the thermal expansion coefficient $\alpha$, the adiabatic sound speed $C_{S}$, the isothermal and isentropic compressibilities $\kappa_{T},\kappa_{S}$, and the Grüneisen parameter $\gamma$, all easily derived from these relations:

[TABLE]

where $\mu=\langle A\rangle m_{H}$ and $\kappa_{T}^{0}=({\rho k_{B}T\over\mu})^{-1}$ denotes the value for a perfect gas of atomic mass $\mu$.

Figure 28 portrays the specific heats for an H/He mixture with solar helium abundance ($Y=0.275$). We recover the same features as in Fig. 14 with the obvious domains of H2 temperature and pressure dissociation and ionization.

The adiabatic temperature gradient is a quantity of prime importance in astrophysics because it is the quantity used in the Schwarzschild criterion to determine whether transport of energy occurs by convection or by microscopic diffusion processes (conduction or radiation). This quantity is displayed in Fig. 29 for the same cosmogonic H/He abundance ($Y=0.275$). The figure spans a $\rho$-$T$ range characteristic of pressure dissociation and ionization from 103 to 106 K and 0.05 to 6.0 $\,{\rm g}\,{\rm cm}^{-3}$. Some typical physical features can be identified in the figure. At low density, the low temperature limit corresponds to the domain of molecular hydrogen with excited rotational levels ($\theta_{rot}\leq T<\theta_{vib}$), $\nabla_{ad}=0.3$, while the high temperature domains corresponds to the perfect monatomic gas, $\nabla_{ad}=0.4$, potentially modified by non-ideal contributions. The two dips reflect the excitation of vibrational levels and the regions of hydrogen dissociation and H or He ionization, respectively. As density increases and dissociation/ionization take place, the two dips vanish and eventually the whole mixture becomes fully ionized for \rho\,\hbox{\hbox{>}\kern-8.00003pt\lower 4.30554pt\hbox{\sim}}\,6\,{\rm g}\,{\rm cm}^{-3}. The spikes around $\log T\approx 3.4$ for $\rho=0.108$ and 0.25 $\,{\rm g}\,{\rm cm}^{-3}$ reflect the onset of hydrogen crystallization. Being a combination of several second derivatives, the adiabatic gradient is very sensitive to the interpolation procedures in the calculations of the EOS table. This is reflected by the wiggles in the domain $4.7\lesssim\log T\lesssim 5.0$, which is the domain of interpolation between the QMD and CP98 calculations. This is the unfortunate consequence of the necessity to combine different calculations in order to construct large enough $(T,P,\rho)$ tables for astrophysical use.

5. Form of the EOS tables

As mentioned in the Introduction, QMD calculations for the EOS have been performed in the canonical ensemble, i.e. with ($T,\rho$) as independent variables. The additive volume law (AVL) procedure to calculate the thermodynamic quantities of the H/He mixture, however, requires the independent variables to be ($T,P$), which imposes transformations of the various quantities from the ($T,\rho$) ensemble into the ($T,P$) ensemble by spline interpolation procedures. Eventually, the ($T,P$) table for the H/He mixture was transformed back into a ($T,\rho$) one, as many astrophysical calculations use these latter quantities as input variables. Although the online H/He table corresponds to a solar (cosmogonic) helium abundance, $Y=0.275$, other mixtures can be easily obtained from the pure H and He tables with use of eqns.(8)-(11).

For reasons of practical interpolation, all tables have rectangular forms with the following limits:

[TABLE]

with grid spacings $\Delta\,\log T=0.05$, $\Delta\,\log P=0.05$, $\Delta\,\log\rho=0.05$, i.e. 121 isotherms, each with 441 values of $P$ or 281 values of $\rho$, and $T$ in K, $P$ in GPa, $\rho$ in $\,{\rm g}\,{\rm cm}^{-3}$. As mentioned in the previous sections, some parts of these tables are unphysical because they correspond to regions in the diagrams that are not handled by the present calculations. These regions concern essentially the domains of solid hydrogen and helium and regions where quantum diffraction effects on the nuclei can no longer be treated by a Wigner-Lirkwood expansion. They are identified in Fig. 1 and 16. Because of the rectangular format of the tables, values at very low density in the $(\log T,\log P$) table and at very low pressure in the $(\log T,\log\rho)$ tables also become unphysical and should not be considered.

Table I gives an example of the various quantities provided by the tables. All quantities are specific quantities, i.e. are given by unit mass, with $A_{H}=1.00794$, $A_{He}=4.00262$ and the atomic mass unit $m_{H}=1.66\times 10^{-27}$ kg. The main second derivatives are also provided. All necessary thermodynamic quantities can be derived from these derivatives, from eqns(5) and (13)-(18). Users, however, might prefer to use only values corresponding to the first derivatives of the Helmholtz free energy, namely $U$, $P$, $S$, and to carry out their own interpolation procedures to calculate second derivatives.

We also stress that the entropy of the spin of the nuclei $S_{nuc}^{id}={\rm Ln}(2s+1)$ $k_{B}$/proton, where $s$ is the spin of the nucleus, is not included in the calculations.

6. Conclusion

In this paper, we have presented new equations of state for pure fluid hydrogen and helium as well as for hydrogen/helium mixtures within the so-called additive volume law (AVL) approximation, i.e. simply taking into account the ideal mixing entropy contribution between the two species (H and He) to the thermodynamic quantities of the mixture. The calculations combine first-principle calculations, based on quantum molecular dynamics (MD-DFT) simulations, in the regime of pressure ionization, with semi-analytical calculations in the low density (molecular/atomic) and high density or temperature (fully ionized) regimes, and provide not only the pressure, internal energy and density but also the entropy and all necessary thermodynamic derivatives. The initial calculations are performed in the canonical ensemble, implying ($T$, $\rho$) as independent variables, and are transformed into ($T,P$) tables to be able to make use of the AVL for the mixture. Therefore, we provide tables in both sets of independent variables. The EOS tables cover a wide temperature-pressure-density domain which permits the calculations of the mechanical and thermal (cooling) structures of a large variety of astrophysical bodies, from solar-type stars to low-mass stars, brown dwarfs and (solar and extrasolar) gaseous planets down to Saturn-like masses. These calculations should supersede the previously widely used Saumon-Chabrier-vanHorn EOS in this domain. At higher densities and/or temperatures, the EOSs merge with those of Chabrier & Potekhin (1998) and Potekhin & Chabrier (2000) ones, devoted to the physics of compact, relativistic bodies such as white dwarfs and neutron stars.

These calculations are by no means without flaws and limitations. Flaws include unphysical numerical oscillations, notably in the calculations of second-derivative thermodynamic quantities, due to spline interpolations. For this reason, verifications of the Maxwell relations between thermodynamic derivatives would be meaningless, because they would undoubtedly be affected by the numerical interpolation procedures and thus have no real physical foundations (see, e.g. §8 and Figs. 18 and 19 of SCvH). Possible future improvements in these numerical treatments will be indicated in future versions of the online EOS tables. Note also that QMD simulations retain as well some degree of uncertainty, inherent to the exchange-correlation functional used in the calculations. Indeed, it has been shown that for liquid hydrogen, pressures obtained with the Perdew-Burke-Ernzerhof functional (Perdew et al. 1996) and by the van der Waals functional (Lee et al. 2010) functionals can differ by as much as $\sim 10$-20% for a given density in the present domain of interest (Morales et al. 2013; see also Pierleoni et al. 2016, Knudson & Desjarlais 2017, Mazzola et al. 2018).

The most challenging limitation of the present calculations is the use of the AVL in the treatment of the H/He mixture, which omits interactions between hydrogen and helium species. While relatively inconsequential for bodies with internal entropies larger than about 10 $k_{B}/e^{-}$ ($\sim 11\,k_{B}/at$ ), i.e. about $7\times 10^{-2}$ MJ kg*-1* K*-1*, as seen from the comparisons with Militzer & Hubbard (2013) (see Fig. 27), this contribution becomes important for cooler entropy values, which encompasses essentially all objects in the brown dwarf and planetary domain older than a few gigayears. Incorporating these non-ideal H/He contributions into the present calculations will be explored in the near future.

The H, He and H/He EOS tables for a solar mixture ($Y=0.275$) are available on the site: http://perso.ens-lyon.fr/gilles.chabrier/DirEOS

The authors are deeply indebted to Andreas Becker and Ronald Redmer for providing tables of their EOS calculations and to Burkhard Militzer for making available his own calculations, which enabled us to perform detailed comparisons between the different EOSs. We also thank Stéphanie Brygoo, Paul Loubeyre and Markus Knudson for sending us their data. This work has been partly supported by the Programme National de Planétologie (PNP).

REFERENCES

Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P., 2003, A&A, 402, 701

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G., 2015, A&A, 577, 42

Becker, A., Nettlelmann, N., Holst, B., & Redmer, R., 2013, Phys. Rev. B, 88, 45122

Becker, A., et al., 2014, Astrophys. J., Suppl. Ser. , 215, 21

Belov, S., et al., 2002, JETP Lett., 76, 443

Bonev, S., Schwegler, E., Ogitsu, T., & Galli, G., 2004, Nature, 431, 669

Boriskov, G.V., et al., 2003, Dokl. Phys., 48, 553

Boriskov, G.V., et al., 2011, CoPP, 51, 339

Brygoo, S. et al., 2015, J. Appl. Phys., 118, 195901

Caillabet, L., Mazevet, S., & Loubeyre, P., 2011, Phys. Rev. B, 83, 4101

Celliers, P., et al., 2010, Phys. Rev. Lett., 104, 184503

Chabrier, G., 1990, J. de Physique, 51, 1607

Chabrier, G., 1993, Astrophys. J. , 414, 695

Chabrier, G., & Potekhin, A., 1998, Phys. Rev. A, 58, 4941

Chabrier, G., & Ashcroft, N., 1990, Phys. Rev. A, 42, 2284

Chabrier, G. & Baraffe, I., 2000, ARA&A, 38, 337

Collins, G.W., et al., 1998, Science, 281, 1178

Datchi, F., Loubeyre, P., & Le Toullec, R., 2000, Phys. Rev. B, 61, 6535

Debras, F. & Chabrier, G., 2019, Astrophys. J. , in press

Deemyad, S., & Silvera, I., 2008, Phys. Rev. Lett., 100, 155701

Eremets, M., & Trojan, I., 2009, ZhETF 89, 198

Fortov, V., et al., 2007, Phys. Rev. Lett., 99, 185001

Hicks, D.G., et al., 2009, Phys. Rev. B79, 014112

Holst, B., Redmer, R., & Desjarlais, M., 2008, 77, 184201

Hu, S.X., Militzer, B., Goncharov, V., & Skupsky, S., 2011, Phys. Rev. B, 84, 224109

Kechin, V., J. Phys.: Condens. Matter, 7, 531

Knudson, M., 2004, Phys. Rev. B, 69, 4209

Knudson, M. & Desjarlais, M., 2009, Phys. Rev. Lett., 103, 225501

Knudson, M., & Desjarlais, M., 2017, Phys. Rev. Lett., 118, 035501

Lee, K., Murray, E., Kong, B., Lundqvist, B., & Langreth, D., 2010, Phys. Rev. B, 82, 081101

Loubeyre, P., Le Toullec, R. & Pinceaux, J-P., 1993, Phys. Rev. Lett., 70, 2106

Loubeyre, P. et al., 2012, Phys. Rev. B, 86, 144115

Lorenzen, W., Holst, B., & Redmer, R., 2009, Phys. Rev. Lett., 102, 115701

Lorenzen, W., Holst, B., & Redmer, R., 2011, Phys. Rev. B, 84, 235109

Mazzola, G., Helled, R., & Sorella, S., 2018, Phys. Rev. Lett., 120, 5701

Miguel, Y., Guillot, T., & Fayon, L., 2016, A&A, 596, 114

Militzer, B., 2009, Phys. Rev. B, 79, 5105

Militzer, B., 2013, Phys. Rev. B, 87, 4202

Militzer, B., & Ceperley, D., 2001, Phys. Rev. A, 63, 6404

Militzer, B., & Ceperley, D., 2000, Phys. Rev. Lett., 85, 1890

Militzer, B., & Hubbard, W., 2013, Astrophys. J. , 774, 148

Mochalov, M. A., 2010, JETP Lett., 2010, 92, 300

Morales, M., Pierleoni, C., & Ceperley, D., 2010a, Physical Review E., 81, 021202

Morales, M., Pierleoni, C., Schwegler, & Ceperley, D., 2010b, PNAS, 107, 12799

Morales, M., McMahon, J., Pierleoni, C., & Ceperley, D., 2013, Phys. Rev. Lett., 110, 065702

Pierleoni, C., Morales, M., Rillo, J., Holtzmann, M., & Ceperley, D., 2016, PNAS, 113, 4953

Perdew, J., Burke, K. & Ernzerhof, M., 1996, Phys. Rev. Lett., 77, 3865

Potekhin, A., & Chabrier, G., 2000, Phys. Rev. A, 62, 8554

Sano, T., et al., 2011, Phys. Rev. B83, 54117

Saumon, D. & Chabrier, G., 1991, Phys. Rev. A, 44, 5122

Saumon, D. & Chabrier, G., 1992, Phys. Rev. A, 46, 2084

Saumon, D., Chabrier,G., & van Horn, H., 1995, Astrophys. J., Suppl. Ser. , 99, 713

Schöttler, & Redmer, R., 2018, Phys. Rev. Lett., 120, 5703

Soubiran, F., Mazevet, S., Winisdoerffer, C., & Chabrier, G., 2012, Phys. Rev. B, 86, 115102

Soubiran, F., Mazevet, S., Winisdoerffer, C., & Chabrier, G., 2013, Phys. Rev. B, 87, 165114

Soubiran, F. & Militzer, B., 2013, Astrophys. J. , 806, 228

Wigner, E., & Huntington, H., 1935, J. Chem. Phys., 3, 764

Winisdoerffer, C., & Chabrier, G., 2005, Physical Review E., 71, 026402

Zel’dovich & Raizer, 2002, Physics of Shock Waves and High Temperature Hydrodynamic Phenomena, Dover