Ab initio based equation of state of dense water for planetary and exoplanetary modeling

As a first step toward a multi-phase equation of state for dense water, we develop a temperature-dependent equation of state for dense water covering the liquid and plasma regimes and extending to the super-ionic and gas regimes. This equation of state covers the complete range of conditions encountered in planetary modeling. We use first principles quantum molecular dynamics simulations and its Thomas-Fermi extension to reach the highest pressures encountered in giant planets several times the size of Jupiter. Using these results, as well as the data available at lower pressures, we obtain a parametrization of the Helmholtz free energy adjusted over this extended temperature and pressure domain. The parametrization ignores the entropy and density jumps at phase boundaries but we show that it is sufficiently accurate to model interior properties of most planets and exoplanets. We produce an equation of state given in analytical form that is readily usable in planetary modeling codes and dynamical simulations (a fortran implementation can be found at http://www.ioffe.ru/astro/H2O/). The EOS produced is valid for the entire density range relevant to planetary modeling, for densities where quantum effects for the ions can be neglected, and for temperatures below 50,000K. We use this equation of state to calculate the mass-radius relationship of exoplanets up to 5,000M_Earth, explore temperature effects in ocean and wet Earth-like planets, and quantify the influence of the water EOS for the core on the gravitational moments of Jupiter.


Introduction
With the advent of a new generation of space and ground based instruments, the constraints of the interior structure of planets and exoplanets have been continuously improving. This is for example the case with Jupiter where the Juno space mission (Bolton et al. 2017) is currently measuring gravitational moments to an unprecedented accuracy, or the various transit and radial velocity programs such as HARPS or Kepler that provide, when combined, density measurements for more than 600 exoplanets (Exoplanet Team 2018). These continuously improving observational constraints on the inner structure of planets and exoplanets call for a proportional effort on the modeling side to achieve a better understanding on the nature of these objects. The modeling of planetary interiors directly relies on the thermodynamic properties of matter at the extreme temperature and pressure conditions encountered within a planet. These can reach several hundreds of Mbars (1 Mbar=100 GPa) and up to ⋆ Current address: Lycée Jean Dautet, La Rochelle, France 500 000 K for the brown dwarfs and giant planets several times the size of Jupiter (see, e.g., Baraffe et al. 2010, for review).
Great progress has been made over the past ten years at understanding this extreme state of matter, which is not directly accessible to laboratory experiments, by using first principles or ab initio simulations based on density functional theory (Benuzzi-Mounaix et al. 2014). This computational intensive approach, that can be validated on the limited densitytemperature range accessible to shock or high-pressure experiments, provides a fully quantum mechanical description for the electronic structure of this state of matter without adjustable parameters. With computational resources greatly increasing, this approach provides the most reliable means to calculate the properties of matter in the thermodynamical range most relevant to planetary modeling, extending from the experimentally accessible thermodynamic conditions to the ones where analytical and semi-analytical approaches become valid. This method has been recently applied to provide comprehensive equations of state (EOS) for the two most abundant elements, hydrogen and helium (Caillabet et al. 2011;Becker et al. 2014;Militzer 2013), which A&A proofs: manuscript no. h2o21 brought about a renewed understanding of the internal structure of Jupiter (Nettelmann et al. 2012;Militzer et al. 2016;Wahl et al. 2017;Guillot et al. 2018). A similar effort is underway for water, and ab initio simulations are now probing the physical properties of water at conditions encountered in planetary interiors.
Following the pioneering work of Cavazzoni et al. (1999), Mattsson & Desjarlais (2006) and French et al. (2009) calculated the properties of the superionic phase for dense water at conditions encountered within Uranus and Neptune. For water, the superionic phase is defined as oxygen atoms locked into either a body-centered cubic (BCC) or face-centered cubic (FCC) crystalline structure with the hydrogen atoms diffusing like in a liquid. This particular phase, which appears at pressures and temperatures above the regular solid ice phases, provides electrical conductivities compatible with the unusual magnetic fields observed for these objects (Redmer et al. 2011). Subsequent works attempted to identify the stable solid phase underlying the superionic region of the phase diagram (Wilson et al. 2013;French et al. 2016) and investigated the miscibility of water in a H-He dense plasma anticipated near Jupiter's core (Soubiran & Militzer 2015). While the debate is ongoing regarding the precise localization and nature of the superionic phase for dense water (Millot et al. 2018), there is still no comprehensive EOS of dense water available for planetary modeling.
As the focus in exoplanetary science is now turning to characterize the Earth-like to Neptune-like continuum, there is a great need for an EOS for water that spans thermodynamic conditions ranging from the atmosphere of an Earth-like planet to the core of a giant planet or brown dwarf several times the size of Jupiter. In the first part of the manuscript, we expand on the work of French et al. (2009) and apply ab initio molecular dynamics simulations and its high-pressure hightemperature Thomas Fermi limit to calculate the properties of water up to a density of 100 g cm −3 and reach conditions encountered in these massive objects. We supplement this data set by the free-energy parametrization developed by the International Association for Properties of Water and Steam (IAPWS) 1 (Wagner & Pruß 2002) that provides an accurate account of the behavior of water in the vapor and liquid phases at pressures below 1 GPa. Using these data sets, we built a wide-range EOS that covers the complete thermodynamical state relevant for planetary modeling. We approximate this EOS by an analytical fit of the free energy, whose derivatives simultaneously provide the fits to pressure and internal energy in agreement with the data, and provides an estimation for the total entropy.
We apply this EOS to probe the effect of temperature on the standard mass-radius diagram used to identify exoplanets by considering wet Earth-like, ocean planets and pure water planets. Finally, we use this EOS for dense water to calculate the gravitational moments of Jupiter currently measured by the Juno probe (Guillot et al. 2018).

Ab initio simulations
The EOS developed in the present work is based on ab initio molecular dynamics simulations for densities ρ between 1 g cm −3 and 50 g cm −3 and temperatures T between 1000 K and 50 000 K, complemented with the free-energy parametrization of Wagner & Pruß (2002) at lower densities and temperatures. At ρ > 50 g cm −3 , we used the Thomas-Fermi molecular dynamics (TFMD) simulations (Lambert et al. 2006;Mazevet et al. 2007).

Computational details
To complement the data obtained previously for water using ab initio molecular dynamics simulations (French et al. 2009;Wilson et al. 2013;French et al. 2016), we carried out simulations using the ABINIT (Gonze et al. 2009) electronic structure package. This consists in treating the electrons quantum mechanically using finite temperature density functional theory (DFT) while propagating the ions classically on the resulting Born-Oppenheimer surface by solving the Newton equations of motion. We used the generalized gradient approximation (GGA) formulation of the DFT with the parametrization of the exchange-correlation functional provided by Perdew et al. (1996) We used two sets of pseudopotentials to cover the density range from 1 g cm −3 to 50 g cm −3 . For densities up to 5 g cm −3 , we used two projector augmented wave (PAW) pseudopotentials generated by Jollet et al. (2014). These pseudopotentials are designed to reproduce accurately the all-electrons results obtained for the individual atomic species. This provides a warranty that the use of the pseudopotential does not cause any important spurious effect. For the two atomic species considered here, hydrogen and oxygen, this consists in a cutoff radius of respectively 0.7a 0 and 1.2a 0 , where a 0 = 2 /(m e e 2 ) is the Bohr radius, and an oxygen pseudopotential with the 1s state treated as a core state. To reach densities above 7 g cm −3 , we use the ATOM-PAW (Holzwarth et al. 2001) package to generate pseudopotentials with cutoff radius of r paw = 0.4 a 0 and r paw = 0.6 a 0 for, respectively, the hydrogen and oxygen atomic species. We further find that the oxygen 1s state needs to be included as a valence state to reach the highest density treated, 50 g cm −3 . The accuracy of the two pseudopotentials produced was inferred by directly comparing the cold curves obtained for the individual atomic species in the FCC phase with the corresponding allelectrons calculations (Jollet et al. 2014). This significant reduction in the cutoff radius requires an increase of the plane wave cutoff from 30 to 100 Ha to reach convergence in pressure and energy below 1%.
The convergence tests performed regarding the number of particles in the simulation cell and the k-point grids in momentum space confirm the results reported by French et al. (2009French et al. ( , 2016. We paid particular attention to the influence of the superionic phase and performed calculations using both the FCC and BCC crystallographic structures Wilson et al. 2013 pointed out that a superionic phase where the oxygen ions remain in an FCC rather than BCC structure may be more stable at intermediate temperature. For the EOS points, we used 54 atoms in the BCC superionic phase. For the FCC superionic phase, we used 108 atoms for densities below and up to 15 g cm −3 while we found sufficient to use 32 atoms at the highest densities. For both phases, we performed the simulations at the Γ-point and integrated the equations of motion with a timestep of 5 a.u. (1 a.u. = 0.024 fs). We attribute this slight difference with the simulation parameters reported by French et al. (2016) to the level of accuracy required in their calculations to evaluate the thermodynamic potentials in the superionic FCC and BCC phases. TF-MD this work SI-BCC this work Liquid this work Wagner 2002 Ice VII+X French 2009Liquid French 2009SI-BCC-French 2009 Fig. 1. The phase diagram of dense water obtained using the ab initio and Thomas-Fermi simulations. Each symbol represents a simulation point. The phase state is indicated by a colored symbol according to the legend. Previous ab initio points obtained by French et al. (2009) as well as representative points of Wagner & Pruß (2002) are also shown.

Ab initio simulation results
Figure 1 displays the pressure P and temperature T values at which ab initio simulations were performed. We have also included sample points from the IAPWS free energy formulation (Wagner & Pruß 2002) for completeness. For the internal energy and pressure, we found a good agreement between our calculations and the previous results of French et al. (2009). We have thus directly included these points in our ab initio set. Figure 1 also shows that a superionic phase remains stable up to the highest pressures for temperatures up to 16000 K.
The superionic phase is identified in our molecular dynamics simulations by looking at the mean square displacement of the hydrogen and oxygen ions as a function of time. Figure 1 shows that at ρ > 15 g cm −3 the superionic phase remains stable up to T = 16000 K when considering either the BCC or FCC structures. With the temperature grid used here, this suggests that the superionic-plasma phase boundaries for both the BCC and FCC structures vary slowly in this pressure range and are both located between 16000 K and 20000 K. We also point that the simulations performed here do not allow us to identify the superionic phase that is the most stable in this thermodynamical regime. They do not indicate either whether another superionic phase may exist in this thermodynamical range.
Here, we do not explore further the exact determination of either the superionic-plasma boundary nor the nature of the superionic state. The results previously obtained at low pressures indicated that this issue has little consequences for the EOS (French et al. 2016). To confirm that this remains the case for the entire density range considered here, we show in Table 1 the results obtained for the internal energy and pressure at representative densities and considering both the BCC and FCC superionic states. We see in Table 1 that the pressure and internal energy values agree to within 1.5%. We thus started our simulations with the BCC superionic lattice at ρ > 20 g cm −3 , as this enables smaller simulation cells. We tested using larger simulation cells that convergence is reached for pressure and internal energy up to 50 g cm −3 .

Thomas-Fermi extension
Beyond 50 g cm −3 , we switch from full ab initio simulations to Thomas Fermi molecular dynamics simulations (TFMD) (Lambert et al. 2006). This consists in using the Thomas-Fermi approximation to describe the electrons while propagating the ions on the resulting Born-Oppenheimer surface. In this framework, the kinetic energy operator in the electronic Hamiltonian is replaced by a functional of the density (Martins 2004). This greatly simplifies the calculation as a plane wave basis is no longer needed and the electronic density is obtained by just solving the Poisson equation. This represents the natural highdensity limit to the DFT and hence to the ab initio simulations. We show in Fig. 2 the relative difference for the pressure P and internal energy U between the full ab initio and the TFMD simulations at densities between 34 and 50 g cm −3 . For both quantities, we see that the difference between the two methods is rapidly reducing as the density increases to be well under 1% at 50 g cm −3 . This result clearly shows that the high-density limit is reached and justifies the use of the Thomas-Fermi approximation beyond this density.
We further point out that the Thomas-Fermi approximation requires the use of a regularization potential. We make the choice of using the pseudopotential determined at ρ = 75 g cm −3 and T = 2000 K throughout the entire range of interest. This approximation introduces an uncertainty on the internal energy of a few percent as the regularization formally breaks the transferability of the pseudopotential in density and temperature (Lambert et al. 2006). The internal energy obtained by the TFMD method is adjusted to the ab initio one at ρ = 50 g cm −3 and T = 6000 K. While the simulations were all started in the FCC phase, we make the choice of not recording the stability of the superionic phase as quantum effects for the protons may start to play a nonnegligible role in its stability (French et al. 2016). The effect on the pressure and internal energy is a higher order effect.

Analytical fit of the Helmholtz free energy
The full ab initio simulation results presented in the previous section are used to construct a functional form of the Helmholtz A&A proofs: manuscript no. h2o21 ρ (g cm −3 ) T (K) P BCC (Mbar) ∆P (Mbar) ∆P/P U BCC (eV/atom) ∆U (eV/atom) |∆U/U| free energy valid over the entire density-temperature domain relevant to planetary modeling. Such an analytical fit of the free energy provides a convenient means to combine various data for use in simulations of interior structure or evolution of planets. For water, this includes the ab initio results presented above, valid down to about a few GPa, and the low-energy IAPWS formulation (Wagner & Pruß 2002) that provides the EOS of water at P < 1 GPa in the vapor and liquid phases as constrained by experimental measurements.

Formulation
The parametrization of the Helmholtz free energy is expressed as where each term has its own physical meaning. The first term is the translational (ideal molecular gas) contribution. Here, is the thermal wavelength of a molecule.
In the second and third terms of Eq. (1), are the analytical expressions for the excess free energy in the moderate-density liquid regime (i.e., at ρ 1 g cm −3 and 300 K T 2000 K) and at high densities (ρ ≫ 1 g cm −3 ), respectively, which, together with F tran , provide the fit to the pressure as a function of density through the thermodynamic relation and w(ρ, T ) = 1 1 + (ρ/2.5 g cm −3 + T/3509 K) 4 is an interpolating function, which varies from 0 to 1 and ensures fitting the pressure as a function of density in the entire ρ − T domain considered. In Eq. (5), F e (N e , T, V) is the Helmholtz free energy of the ideal nonrelativistic Fermi gas of N e = n e V electrons at temperature T , and Z * is an effective charge number, which is expressed as an analytical fitting function so as to adjust the pressure derived through Eq. (6) to the pressure data from the ab initio calculations. Explicitly, where is the effective (ideal Fermi gas) electron pressure, λ e = (2π 2 /m e k B T ) 1/2 is the electron thermal wavelength, is the effective electron chemical potential, and In Eq. (9), is the standard Fermi integral, and X ν (I) in Eq. (10) is the inverse function. For both I ν (X) and X ν (I) we use the Padé-type approximations of Antia (1993). In Eq. (11), Γ e and r s are the usual electron Coulomb parameter and density parameter, respectively: Γ e = e 2 /(a e k B T ) in the CGS system, and r s = a e /a 0 , where a e = ( 4 3 πn e ) −1/3 is the electron-sphere radius, and n e = 10 n H 2 O is the total number density of all electrons (free and bound). In Eq. (4), a vdW and b vdW are the van der Waals constants (respectively, 5.524 × 10 12 erg cm 3 mol −2 and 30.413 cm 3 mol −1 , Grigoriev & Meilikhov 1997). The first line in Eq. (4) reproduces the van der Waals EOS, which is sufficiently accurate at ρ ≪ 1 g cm −3 , and the second line adjusts the EOS at ρ ∼ 1 g cm −3 .
The fourth term in Eq.
(1) reads where b 1 = 3 × 10 −13 erg, b 2 = 1.35 × 10 −13 erg, b 3 = 2.43 × 10 −13 erg, and τ = T/T crit = T/647 K. It is derived from the fitting correction to the residual internal energy, This correction does not affect pressure but improves the fit to the internal energy, through the thermodynamic relation Note that we measure the total internal energy from its minimum at the ground state of the molecular phase, so that U > 0 at any ρ and T . This definition is the same as in Wagner & Pruß (2002). It differs from the definition adopted in Table 1 and Fig. 2 by the ground-state energy value of 11.14 MJ/g. It also differs by a constant of 77 kJ/g from the internal energy given in French et al. (2009);French & Redmer (2015); Soubiran & Militzer (2015). The last term in Eq. (1), −S 0 T is an additional correction, which affects neither P nor U, but shifts the entropy by constant S 0 . We find that the value S 0 = 9.8k B N at provides the best fit (within ±0.3k B N at ) to the results presented for S by Soubiran & Militzer (2015) at ρ ≈ (1 − 2.5) g cm −3 and T = (1000 − 6000) K. However, entropy evaluation from our present fit should be used with caution especially when crossing the boundaries of different phases, where one can expect a discontinuity. This may lead to a value of S 0 that differs from one phase to the other. The present analytical fit describes the EOS of liquid water at ρ 1 g cm −3 and T 2000 K, as well as plasma at 1 g cm −3 ρ 10 2 g cm −3 and 10 3 K T ×10 5 K. While not including the super-ionic phase as a different phase, the single-phase approximation used here provides a satisfactory description of the thermodynamical properties in this superionic regime. It has, however, a limited applicability for the ice VII and ice X phases that occurs at T 2000 K in the range (0.02 − 0.5) Mbar P 3 Mbar (Petrenko & Whitworth 1999). We also point out that quantum effects for the ions, which could be relevant at the highest densities and for low temperatures, are not included in the current parametrization. To build a fully multi-phase EOS for water, one can supplement our fit by the parametrizations constructed specifically for the ice and superionic phases (French & Redmer 2015;French et al. 2016). This will be the topic of further work. Finally, we also point out that our analytical fit is less reliable in the domain of thermal ionization and dissociation of molecular water, where ρ ≪ 1 g cm −3 and T ≫ 10 3 K. This regime is indeed poorly constrained by either the ab initio simulations or the IAPWS parametrization.

Validation of the analytical fit
We first verify the ability of the analytical fit at reproducing both the results of theoretical calculations and the IAPWS free energy model. In Figs. 3 and 4, we compare the behavior of pressure and internal energy obtained with the input data at respectively low and high densities and temperatures.
As we are primarily interested in planetary interiors, the accurate description of the liquid-vapor transition below the critical point located at T crit = 647 K and P crit = 22.064 MPa is beyond the scope of this study. Furthermore, these conditions are tied to an accurate modeling of the atmosphere of the planet that do not directly involve an EOS as the one developed here. Figures 3a-b suggest that without atmospheric treatment, interior models should consider the liquid state for surface conditions when the surface temperature is below the critical point. We will expand further on this point in the following sections. At the lowest densities, we see in figure 3a that the pressure turns negative along the 300 K and 600 K isotherms for densities below 1 g cm −3 . Figure 3b indicates that this translates into a minimum for the internal energy. This corresponds to the crossing of the liquid-vapor phase boundary and a region where the pressures are formally negative, which in fact corresponds to phase coexistence. For instance, for the 300 K isotherm the analytical fit gives a region of negative pressures expanding from 0.9 g cm −3 to 0.1 g cm −3 . Note that the IAPWS data (Wagner & Pruß Fig. 4. Comparison of the data with the analytical fitted isotherms at high densities for the pressure normalized to the ideal atomic gas value, P/n at k B T (a), and for the internal energy U (b). Symbols show the data: the IAPWS table (Wagner & Pruß 2002, straight crosses); results of ab initio calculations by Soubiran & Militzer (2015) (empty diamonds); ab initio calculations by French et al. (2009) (oblique crosses for liquid, inverted triangles for ice X, filled diamonds for the superionic [SI] phase, asterisks for the plasma phase); and the results of our present ab initio calculations (filled triangles for the SI phase and filled dots for the plasma phase), supplemented with our TFMD calculations at ρ > 50 g cm −3 (empty circles). Solid lines represent the present fit; the dotted line represents the fit of French & Redmer (2015) for ice X at T = 1000 K.
2002) indicate a wider phase-transition region, with the low density ∼ 10 −3 g cm −3 on this isotherm. The agreement improves for the 600 K isotherm. We see that the analytical fit reproduces the overall behavior of the pressure across the liquid-vapor boundary. However, it extends this boundary to a higher temperature, giving the critical point at 683 K and 0.331 g cm −3 (to be compared with the experimental values of 647 K and 0.322 g cm −3 ).
Figures 3a-b also indicate that the agreement with the ab initio data at higher temperatures is satisfactory up to 6000 K. We note that the ab initio results and the free energy model predictions are not in perfect agreement at 1000 K. As the ab initio method becomes less reliable as density decreases, mainly due to the deficiency of density functional theory in under-dense regime, the ab initio results fail to match the IAPWS formulation at low density. Our analytical fit eliminates this mismatch by interpolation between the low-density IAPWS and high-density ab initio data. Figures 4a-b show the data and fitted isotherms for the pressure and internal energy across the entire density range at higher temperatures, 1000 K ≤ T ≤ 50 000 K. Figure 4a displays the pressure normalized to the atomic ideal gas contribution n at k B T . As noted above, the fit does not perfectly reproduce the ab initio data for the ice X phase along the 1000 K isotherm at 2.5 g cm −3 ≤ ρ ≤ 4 g cm −3 . However, it satisfactorily reproduces the thermodynamical properties, despite the fact that the system crosses a number of various phases in the ρ−T domain displayed in the figures.
We also point out that the data represented in Figs. 3 and 4 by empty symbols (the ice phase, the data by Soubiran & Militzer 2015, and the Thomas-Fermi results) have not been explicitly included to construct the analytical fit. We see that the data by Soubiran & Militzer (2015) is in good agreement with the prediction of the analytical fit. Furthermore, the good agreement found with the Thomas-Fermi results up to 100 g cm −3 shows the validity of the analytical fit up to this high density. Fig. 5 summarizes the applicability of the current analytical EOS for planetary modeling. We display the accuracy of the analytical fit as a color code corresponding to the residual difference between the predicted and input P and U values. On the same figure, we also show representative interior profiles for Jupiter, Saturn, and Neptune that display a significant amount of water in their interior. This shows that the analytical fit is accurate for modeling these objects. For comparison, the profile of a 9M J planet is also shown. As pointed out before, the single-phase approximation used for this analytical fit ignores the discontinuities due to the phase changes between the different phase states. For the phase transition between the liquid state and ice X, the discrepancies increase to tens percent. Otherwise we see that the fit remains a reasonable approximation for both the pressure or internal energy throughout the relevant thermodynamical domain.

Comparison with experimental data and previous EOSs
We now turn to compare the predictions of the analytical fit developed in Sect. 3 with existing experimental data from both static and dynamical experiments. We compare these predictions with EOSs commonly used in planetary interior models. In Fig. 6, we compare the predictions of the analytical fit developed with the static high pressure data obtained using diamond anvil cells (Hemley et al. 1987;Sugimura et al. 2008). We also show in Fig. 6 the 1000 K isotherm and its corresponding ab initio data to illustrate further the approximation made by not accounting explicitly for the solid ice phases and extending the liquid throughout the solid phases. We see that the analytical fit thus misses the jump between the liquid and ice X phases at ρ ∼ 2.5 g cm −3 , as we have already seen in Fig. 4. The T = 300 K isotherm behaves similarly. We see that both the IAPWS free energy formulation and our analytical fit, being  Baraffe et al. 2008Baraffe et al. , 2010. continued from the low-density region at T = 300 K, overestimate the pressure (underestimate the density) in the ice phases at higher densities. At P = 10 GPa, the resulting density is about 7% lower than in the ice phases. This should be compared with the predictions of Seager et al. (2007) who used T = 0 K DFT calculations for the various ice phases to construct their EOS. By neglecting the effect of temperature and the phonon contribution, they overestimate the density of dense water by about 3% at 10 GPa. We also note that this difference with the experimental data tends to decrease as the pressure increases for both EOSs. Therefore, the most significant difference resulting from neglecting the ice phases that we can anticipate for interior structure calculations could be for the pressure profile of the outer layer of a planet, if it had T ∼ 10 3 K at ρ ∼ 3 g cm −3 . From Fig. 5 we see that it is not the case for the giant planets, for which the temperature is much higher and thus the temperature profile passes well above this phase jump.
We also point out some differences between the two ab initio calculations beyond 4 g cm −3 . We attribute this difference to the use of different functionals in the Thomas-Fermi approximation. We will investigate further its impact on the interior structure models in the following sections. Figure 7 shows a comparison between the predictions of our analytical fit deduced from ab initio simulations and the measured experimental data along the principal shock Hugoniot line. For a given initial state, the Rankine-Hugoniot relation determines the final states allowed by conservation of energy and where subscript 0 indicates the initial state. The principal Hugoniot corresponds to a single shock, obtained with initial state at rest at normal conditions. Figure 7 shows that shock experiments probe a range of pressures almost an order of magnitude higher than when using diamond anvil cells (Fig. 6). This also corresponds to a significant increase in temperature. The temperature reaches about 5 000 K around 100 GPa for a shock, while it remains near 300 K in diamond anvil cell experiments. With the increase of the shock pressure beyond 500 GPa, the temperature exceeds 50 000 K. Since the input ab initio data that underlie our fit have been obtained for T ≤ 50 000 K, the accuracy of the fit at higher temperatures is not guaranteed (the corresponding part of the Hugoniot line is drawn by long dashes in Fig. 7). Figure 7 shows a good agreement at low pressures with early data obtained using explosion (Volkov et al. 1980) and gas-gun techniques Lyzenga et al. 1982). The unique measurement of the water EOS at P > 1000 GPa, published long ago by Podurets et al. (1972), is satisfactorily described by the above-mentioned continuation of our fit beyond the range where it has been constrained by the data. The first laser shock data obtained in the 100 to 1000 GPa range by Celliers et al. (2004) are much softer than the analytical fit. This experimental data set is in good agreement with SESAME 7150 predictions. In contrast, we see that the analytical fit agrees very nicely with the more recent experimental data of Knudson et al. (2012) obtained using the Z-pinch techniques, as well as with the latest results of Kimura et al. (2015). This confirms that earlier laser shock experiments likely suffer from systematic errors which could be caused by the standard used in the impedance matching method (Knudson & Desjarlais 2009). This experimental set will thus not be considered to validate the behavior of water at high pressures and temperatures. This also highlights that the SESAME 7150 predictions, that are in good agreement  Volkov et al. (1980);; Lyzenga et al. (1982); Celliers et al. (2004); Knudson et al. (2012), and Kimura et al. (2015). For comparison, the principal Hugoniot lines predicted by the SESAME (Lyon & Johnson 1992) and ANEOS (Thomson & Lauson 1972) models are shown by the dotted line and short-dashed line, respectively. The inset shows the principal Hugoniot line in the ρ − T plane calculated from the fit and the experimental data points from Lyzenga et al. (1982) and Kimura et al. (2015).
with the data of Celliers et al. (2004), should be ruled out for planetary modeling. The SESAME 7150 model is not considered a reliable EOS nowadays, but it is still in use for planetary modeling (Miguel et al. 2016). This is also the case for the ANEOS model (Thomson & Lauson 1972). Figure 7 shows that ANEOS predictions depart from the data of Knudson et al. (2012) at ρ > 2.5 g cm −3 significantly outside the experimental errorbars. Figure 8 shows a comparison between the predictions of our analytical fit and the EOS measurements obtained using doubleshock experiments (Volkov et al. 1980;Knudson et al. 2012), where the initial shock wave is reflected from a surface of a standard material (aluminum or quartz). In this case, the initial state in Eq. (17) lies on the principal Hugoniot line and represents the initial condition for the secondary Hugoniot. Since the latter initial condition is measured with some uncertainties, the position of the secondary Hugoniot is not firmly defined. We therefore show the 1σ limits for each secondary Hugoniot that arise from these uncertainties in the initial state. Figure 9 shows a comparison between the predictions of our analytical fit with two EOS measurements in experiments using laser-induced shock in statically compressed water. We find a better agreement with the latest dataset of Kimura et al. (2015) compared to the ones of Lee et al. (2006). We show these for completeness as the scatter in the experimental data cannot further constrain the validity of the analytical fit.  Volkov et al. (1980) and of Knudson et al. (2012), respectively. Fig. 9. Comparison between the calculations based on our analytical fit and the experimental shock data for initially precompressed water. (a) Points measured by Lee et al. (2006); (b) points obtained by Knudson et al. (2012). The inset in each panel shows the initial pressures and densities measured.
Overall, we find that our present analytical fit provides a satisfactory description of both the static and dynamical experimental results available to date for dense water. We now turn to illustrate the applications of the EOS developed for the interior structure of different classes of planets.

Implication for planetary interiors
The internal structure calculations are performed by solving the standard hydrostatic, mass, and energy conservation equations (e.g., Schwarzschild 1958;Kippenhahn et al. 2012) where P is the pressure, ρ the density, g = Gm/r 2 the gravity, G is the gravitational constant, m is the mass enclosed within a sphere of radius r, and ∇ T = d ln T/d ln P depends on the mechanism of energy transport. If the medium is stable against convection, then where T eff is the effective surface temperature and K is the effective opacity. If the transport of energy is dominated by convection, then in the simplest (Schwarzschild) approximation is the adiabatic gradient. The interior structure of a planet of a given mass, M, is obtained by integrating inward the set of equations (18)-(20), starting from a boundary condition defined by a fixed surface temperature and pressure. This is formally given by either in-situ measurements for planets of the solar system or full atmospheric calculations for exoplanets. As we are primarily interested in testing our EOS for dense water, we will devise strategies to overcome this difficulty for each type of planets that we will consider.

Water planets
In order to test the accuracy of our EOS, we first turn to the calculation of the inner structure and mass-radius relationship of a planet entirely made of water. This is a purely academic exercise disconnected from the outcomes of planetary formation but the resulting mass-radius relationship remains a well used benchmark to classify exoplanets (Seager et al. 2007). It also allows us to decipher the influence of temperature on both the inner structure and mass-radius relationship. Figure 10 shows the effect of temperature in isothermal interior structure models for planet of respectively 0.5 and 5 M Earth . In this situation, the temperature throughout the planet is constant and equal to the surface temperature, T surf . The outer boundary conditions are chosen as 1 bar for temperature above the critical point and at the liquid-vapor boundary below. This corresponds to different surface densities and close to 1 g cm −3 for the two models using a surface temperature below the critical value. Figure 10 shows that within this model, the effect of temperature is twofold. As the EOS developed fully accounts for temperature throughout the density range covered by the conditions existing in planetary interiors, it affects the extent of the outer edge boundary as well as the compressibility deeper within the planet. Figure 10 shows that both these effects are important for low mass planets where we notice a radius increasing by a factor of two when the surface temperature varies from 300 K to 2000 K. This increase in radius comes from a significant expansion of the low density outer edge beyond 1.5R Earth that follows the rapid expansion of the supercritical liquid at high temperatures and low surface gravity. It also comes from the significant temperature dependence of the water EOS for densities below 2 g cm −3 previously pointed out in Fig. 6. For a planet of 0.5M Earth , the maximum pressure reached at the center of the planet is 0.18 Mbar for a temperature of 2000 K and 0.27 Mbar for a temperature of 300 K. Figure 10 shows that the effect of temperature decreases significantly when the planet size increases. This comes from both a smaller expansion of the outer edge that comes from a larger surface gravity and a reduced temperature dependence of the EOS as the density increases. Figure 10 shows that the density profile of the planet is dominated by densities between 2 g cm −3 and 4 g cm −3 for a planet of 5M Earth . Figure 6 shows that the temperature dependence of the EOS is almost negligible in this density range. The effect of temperature is thus reduced to a small expansion of the outer edge that leads to an increase of 10% of the radius. This suggests that the effect of temperature becomes negligible as the mass of the planet increases. Figures 11a and 11b give the mass-radius relationship obtained with the current EOS and calculated for the complete range of planets detected to date, and a zoom for planets below 10 M Earth (Exoplanet Team 2018). As anticipated above, we see that the temperature dependence decreases as the mass of the planet increases. Figure 11b shows that the temperature effect is rather important for planets smaller than 10M Earth . We also see in figure 11a that the temperature dependence for pure water planets can be neglected for planets larger than 15M Earth when considering isothermal temperature profile. We also find that the mass-radius relationship obtained with the EOS developed in the current work is consistent with the calculations of Seager et al. (2007) for the entire range of planets detected. Fig-1 10 100 Fig. 11. Mass-radius relationship for pure water planets in an isothermal model. The surface temperature and pressure are indicated in the figure and compared to the result of Seager et al. (2007). (a) over the entire range of planets currently detected. (b) Same as (a) for planets with mass less than 10M Earth . Isothermal models for pure silicates (MgSiO 3 ), iron (Fe) (Bouchet et al. 2013;Mazevet et al. 2015) ures 11a and 11b show that the radius obtained with the EOS developed here and for an isothermal model with surface temperature of 250 K is 3% higher than the results of Seager et al. for planets less than 10M Earth and 2% lower above this planetary mass. This latter result comes from the difference in the behavior of the two EOS at high densities already mentioned in the previous section. For low mass planets, where inspection of Fig. 7 indicates the largest difference between the EOS developed here with both the experimental data and the zero temperature EOS of Seager et al., we also find a rather satisfactory agreement. This indicates that neglecting the ice phases has a minimal impact on the resulting mass-radius relationship even for low mass planets. For benchmarking purposes, we also show in Fig. 11b (Bouchet et al. 2013;Mazevet et al. 2015). The benchmark calculations for wet super-Earth are from Thomas & Madhusudhan (2016). the current EOS by considering a surface pressure of 10 GPa (0.1 Mbar). Finally, we notice in Fig. 11b that several planets detected fall within the temperature dependent range of the isothermal pure water model. This suggests that temperature-dependent effect on the mass-radius relationship needs to be included when considering the interior structure to identify the nature of these objects. We will pursue further this suggestion using more realistic interior structure models and by considering wet super-Earths and ocean planets.

Wet super-Earth planets or ocean planets
Wet super-Earth or ocean planets are objects that have no equivalent in the solar system. They have been introduced to interpret the continuum of planets detected between pure water and Earthlike object. Earth-like planets consist of objects of varying mass but following the Earth composition (33% Fe, 66% MgSiO 3 by mass). Wet Super-Earths consist of an Earth-like core made of 33% iron and 66% silicates and a significant fraction of water outside the core (Valencia et al. 2006;Thomas & Madhusudhan 2016). To test the validity of our EOS at describing these objects, we show in figure 13 the mass-radius relationship obtained by considering isothermal models and wet super-Earth planets constituted of 50% water. The iron and silicates EOSs used are temperature dependent and, similarly to the water EOS developed here, are described by free energy functional forms adjusted directly to ab initio calculations (Bouchet et al. 2013;Mazevet et al. 2015). Figure 12 shows the mass-radius relationship for isothermal models, obtained using two different surface temperatures. For a surface temperature of T surf = 500 K, we obtained a rather satisfactory agreement with the previous calculations of  Fig. 13. Temperature dependence of the mass-radius relationship for adiabatic models of wet super-Earth planets containing 50% water as solid lines. The surface temperatures and pressures are indicated in the figure. The dots represent detected planets in this range of planetary mass. The isothermal models for pure water planets (dashed) and wet super Earth (purple) with the corresponding surface temperatures are also displayed in the figure. These correspond to the data shown in, respectively, Fig. 11 and Fig. 12. Thomas & Madhusudhan (2016). We point out that the outer boundary condition used in the current work follows a somewhat different prescription. As mentioned before, we do not consider the vapor phase for surface temperature below the critical point and use the liquid-vapor boundary at 1 bar as outer boundary conditions. Figure 12 shows that the temperature dependence of the mass-radius relationship is rather significant and can not be neglected when assessing the interior structure of these objects, in agreement with the conclusions of Baraffe et al. (2008). Indeed, we see in Fig. 12 that several planets lie between the boundaries delimited by the two surface temperatures. As for the case of the pure water planets, the temperature dependence decreases as the mass of the planet increases and becomes negligible beyond 15M Earth . The wet Earth-like model that consists by 50% of water is at a midpoint between the dry Earth-like planets and pure water planets. By extension of the temperature dependence shown for planets composed by 50% of water, we can anticipate that this also needs to be taken into account to deduce the composition of objects such as Kepler-18b or 55Cnc that hold a significant fraction of water. In this case, the amount of water deduced will directly depend on the surface temperature considered. Along the same line of reasoning, we can also anticipate that the overlap between the pure water and wet Earth-like planet with a surface temperature of T surf = 2000 K shown below 2M Earth for a composition of 50% water will expand to planets with larger mass when the fraction of water is increased. To make a more quantitative statement on this issue requires to go beyond the internal structure calculations and to account for an accurate modeling of the atmosphere as well as the amount of light received from the host star. Both these topics are beyond the scope of this study that aims at validating the EOS for dense water developed here. We will instead turn to the evaluation of the effect of temperature on the mass-radius relationship beyond the simple isothermal model where the temperature is kept constant throughout the planet.
The EOS developed in the current work provides the total entropy thanks to the parametrization of the Helmholtz free energy using the ab initio results. This allows us to calculate more realistic interior models by considering a water layer undergoing nearly adiabatic convection. The uncertainty in defining S 0 pointed out in Sect. 3 does not affect these results as long as the convective layer remains within a single thermodynamic phase. The previous isothermal models and the adiabatic ones bracket the maximum impact of the water EOS upon the mechanical structure of the body. We see in Fig. 13 that the effect of temperature on the mass-radius relationship is almost two times more important when considering the water layer as adiabatic rather than isothermal. When considering a surface pressure of 1 bar, Fig. 13 shows that the radius obtained up to 10M Earth is significantly bigger than the pure water case at zero temperature. This effect is the most spectacular for surface temperature above the critical temperature. It remains limited when the surface temperature is below the critical point. We also point out that the surface pressure tends to reduce this effect. Figure 13 shows that the radius decreases by 20% at 1200 K when the surface pressure increases from 1 bar to 100 bar. This suggests that neglecting the effect of temperature when identifying the internal structure of exoplanets tends to over-estimate the overall amount of water, especially for objects close to their parent star and receiving a significant amount of light. This needs to be tempered by the probable escape of the atmosphere under these particular conditions and will need to be assessed on a case by case basis with a proper treatment of the atmosphere.

Core of giant planets
We now turn to the last situation where a temperature-dependent EOS for dense water is important for planetary modeling, the core of giant planets. Following the well accepted core-accretion scenario (Pollack et al. 1996), a significant amount of water is expected in the core of giant planets. This potentially significant amount of water stems from the likely composition of the initial core that triggered the accretion of a large fraction of hydrogen and helium. The exact amount of water as well as the size of the core for a planet like Jupiter is a matter of debate and one of the main scientific goal for the Juno mission (Bolton et al. 2017). The probe is currently measuring Jupiter's high order gravitational moments for Jupiter to decipher the amount of metallic element present in the interior as well as their distribution throughout the planet (Wahl et al. 2017). This goal requires an accurate modeling of the EOS for all the elements potentially constituting the planet, including water. Figure 14 shows the dependence of Jupiter first two gravitational moments, J 2 and J 4 , on the dense water EOS used in the modeling of the interior. Jupiter's interior is obtained by solving the standard hydrostatic equilibrium equations (18)-(20) and by considering the planet as composed of an H-He envelope and a pure water core (note that we assume a fully adiabatic interior profile). The hydrogen and helium EOSs used for the envelope are also based on ab initio simulation results (Caillabet et al. 2011;Soubiran 2012). The gravitational moments are calculated using the theory of figures to the third order (Zharkov & Trubitsyn 1978). Figure 14 shows the values of the first two gravitational moments obtained using two different water EOSs, namely the present and the widely used ANEOS ones, to describe the core and by considering two different helium concentrations. The size  of the core is varied to the values indicated in the graph. The first mass fraction of helium, Y He = 0.2384, corresponds to the Galileo measurements, while a slightly higher mass fraction, Y He = 0.26, allows us to reproduce the values of the two gravitational moments J 2 and J 4 . Within a two layers model of Jupiter, the first case reproduces the observed radius of the planet as well as the value of J 2 measured but, as shown in figure 14, misses the J 4 moment by slightly more than 1%. Conversely, in the case where the mass fraction of helium is increased to Y He = 0.26, the observed radius is underestimated by slightly more than 1% while matching the first two gravitational moments. This shortcoming of the two layer model for Jupiter, namely a central core surrounded by a homogeneous gaseous envelope, has been well documented elsewhere (Miguel et al. 2016). We point out here that this issue is not resolved by using different water EOSs for the core. We see in figure 14 that using different water EOSs leads to different predictions regarding the size of the core. These two predictions in the size of the core differ by close to 20%, and do not depend on the helium concentration. We note that this also translates into different average and maximum densities reached in the central core. The EOS developed here predicts a pure water core density more than 25% higher than in the case of the widely used ANEOS one. The highest density reached in the former case is close to 12.45 g cm −3 , while it remains close to 9.67 g cm −3 in the latter case. This difference is comparable with the variation of the core density when a pure ice core is replaced by a pure rock core, as pointed out by Guillot (1999). In our case, this can be traced back to significant discrepancies between the two EOSs, which predict different compressibilities for pressures and temperatures relevant to Jupiter's inner core. These conditions correspond to pressures above 40 Mbar and temperatures between 15 000 K and 20 000 K depending on the hydrogen-helium EOSs used Miguel et al. (2016). In this thermodynamical region, the two water EOSs predict pressures that differ by more than 10%. This shows that the current EOS is useful to calculate accurately the exact size of the core that is potentially present at the center of giant planets.

Summary
In summary, we developed an equation of state for water applicable to the full range of thermodynamic conditions relevant to planetary modeling. This encompasses the range from outer layers of wet super-Earth to the core of giant planets. This equation of state includes an evaluation of the entropy (within an arbitrary constant value S 0 , Sect. 3, which can differ from one phase to the other but has no impact on all other thermodynamic quantities) thanks to the parametrization of the Helmholtz freeenergy functional form on the ab initio results. Using this EOS, we show that the temperature dependence of the EOS needs to be accounted for when analyzing the composition and nature of exoplanets using the standard mass-radius relationship, a point already stressed by Baraffe et al. (2008). We also show that an accurate description of the thermodynamic properties of dense water is required to deduce the mass of the core of giant planets from gravitational moments as it is currently measured for Jupiter by the Juno mission. As discussed in Sect. 3, the EOS developed here is not very accurate for the solid ice phases. This does not noticeably affect the planetary gross properties, such as the mass-radius relationship, but may produce an error of several percent in calculations of temperature profiles. To include explicitly all the ice phases as well as treating the super-ionic phase will be the topic of further development of the present EOS.