Material Properties of Saturn’s Interior from Ab Initio Simulations

We investigate the material properties of a mixture of hydrogen, helium, and oxygen representative of Saturn’s interior at pressure–temperature conditions of a recent Saturn model (see Mankovich & Fortney) with molecular dynamics simulations based on density functional theory. Their model considers the demixing of hydrogen and helium and predicts a He-rich layer above a diluted core. We calculate the thermodynamic and transport properties and discuss the impact on Saturn’s evolution and interior structure. We find a significant impact of the He-rich layer on the specific heat capacity, speed of sound, viscosity, diffusion coefficients, thermal and electrical conductivity, Lorenz number, and magnetic and thermal diffusivities.


Introduction
Saturn is the second-largest planet in our solar system.Composed mainly of hydrogen and helium, Saturn belongs to the class of gas giants such as Jupiter.Exploration of Saturn with space missions has increased our knowledge about its interior structure, evolution, and magnetic field enormously over the past decades.Pioneer 11 (1979) and Voyager 1 (1981) and 2 (1981) visited Saturn only briefly in flybys (Hubbard & Stevenson 1984).The Cassini spacecraft, the first to orbit Saturn (2004Saturn ( -2017)), performed a long-term in situ exploration of the planet, including its rings and moons (Ingersoll 2020).Cassini also carried the Huygens probe that landed on Saturn's moon Titan in 2005.
Among other findings, the Cassini mission drastically improved our knowledge of the Saturnian magnetosphere (Simon et al. 2015), rotation rate, and composition.However, many important issues are still not well understood (Baines et al. 2019).For instance, the enrichment of heavy elements compared to the Sun and Jupiter, as well as their distribution throughout the planet are under debate.This problem is closely related to the size and composition of the core.A key question is the region of phase separation of helium from hydrogen inside the planet and how this process impacted the evolution of Saturn.
Helium rain (Stevenson 1975), the formation of He-rich droplets and the subsequent descent of the droplets, might explain the depletion of helium in the outer atmosphere and act as an additional heat source due to the release of gravitational energy (Stevenson & Salpeter 1977a;Fortney et al. 2007;Püstow et al. 2016).However, in order to form helium droplets, the H-He mixture needs to undergo a demixing process which heavily depends on the free enthalpy of mixing.
Schöttler & Redmer (2018, hereafter SR18) calculated the demixing of H-He mixtures at conditions relevant to giant planets.Their work improved on the pioneering studies of Lorenzen et al. (2009) and Morales et al. (2009).The exact location of the demixing region in P-T space is under lively debate after predictions based on shock-wave compression experiments with precompressed targets have been reported (Brygoo et al. 2021).Mankovich & Fortney (2020, hereafter MF20) developed a Saturn evolutionary model based on a shifted SR18 demixing diagram, assuming efficient helium rain such that the mixingratio profile of helium to hydrogen instantaneously follows the thermodynamic equilibrium (i.e., saturation) curve given by the Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
demixing diagram.As a result, a helium gradient in the P  100 GPa region grows as the planet cools, the helium being increasingly sequestered at depth as a function of time, as described previously by a number of publications (Stevenson & Salpeter 1977b;Hubbard et al. 1999;Fortney & Hubbard 2003;Püstow et al. 2016).MF20 were the first to apply the SR18 demixing diagram, a major implication of which is that the helium gradient would quickly extend down to Saturn's heavy element core (or center, in a model lacking a core).After that time, excess He-rich droplets raining from above would accumulate above the core-mantle boundary, forming a ∼17 M ⊕ ocean of nearly pure helium by the present day.While the details would be affected by a more realistic accounting of Saturn's heavy element distribution (e.g., as informed by ring seismology; Mankovich & Fuller 2021), the MF20 model raises important questions: How does a deep layer with only trace amounts of hydrogen affect the profiles for the conductivity, pressure, temperature, and density?And ultimately, what are the implications for the all-important dynamo generation process in Saturn?
Earlier studies on the material and transport properties of Jupiter by French et al. (2012) and brown dwarfs by Becker et al. (2018) provided valuable data for improved planetary models.However, to this date, similar data for Saturn are still missing.The main purpose of this work is to provide a strong basis of material properties for future dynamic models of Saturn's interior, especially for magnetic field simulations (Wicht & Tilgner 2010;Wicht et al. 2018;Yan & Stanley 2021;Yadav et al. 2022).So far, the origin of Saturn's observed strongly dipolar and axisymmetric magnetic field is not well understood (Smith et al. 1980;Schubert & Soderlund 2011).Its apparent geometry is in conflict with Cowling's theorem (Cowling 1933), which states that a stationary magnetic field cannot have the same symmetry axis as the rotation.An appropriate database of material properties based on Saturn's differentiated interior structure is a necessary prerequisite to resolve this question with future magnetohydrodynamic simulations (Stevenson 1982).
Section 2 outlines details of the MF20 Saturn model that this work is based on.We used the model's mass fractions of hydrogen, helium, and heavier elements (we chose oxygen as being representative of the heavier elements) and calculated thermophysical properties along the corresponding P-T path.Details of the ab initio simulation method are described in Section 3. Results for the thermal expansion coefficient α, isothermal compressibility κ T , heat capacities c P and c V , sound velocity c s , Grüneisen parameter γ, self-diffusion coefficients D H , D He , and D O , dynamic and kinetic shear viscosities η and ν, thermal conductivity λ, thermal diffusivity κ, electrical conductivity σ, magnetic diffusivity β, and Lorenz number L are shown in Section 4. We discuss characteristic changes in these quantities along the radius on the basis of the different layers present in Saturn.The pressures and temperatures of our simulations range from 1.22-1359 GPa and 1966-9051 K, corresponding to radius coordinates of 0.95-0.25 R S .We summarize our results in Section 5.

Saturn Model
The MF20 reference Saturn structure model is a solar age model from an evolutionary sequence based on the SR18 demixing diagram and the Militzer (2013) equation of state (EOS) for hydrogen and helium, blended with the Saumon, Chabrier, and van Horn (Saumon et al. 1995) EOS for pure helium to enable the modeling of arbitrary mixtures (following Miguel et al. 2016).The model adopts an additive offset of +540 K to all SR18 demixing curves, a step that MF20 took to guarantee consistency with Jupiter's known atmospheric helium depletion at the solar age (von Zahn et al. 1998).
The MF20 models considered only simplistic heavy element distributions consisting of a pure water core and a uniform water mixture in the hydrogen-/helium-dominated envelope.In the meantime, evidence from ring seismology has lent strong support for the existence of large-scale gradients in Saturn's heavy element abundance (e.g., Mankovich & Fuller 2021).No model yet exists that addresses both realistic heavy element distributions and detailed expectations from a hydrogen/ helium demixing diagram.
The present model has a pure water core of 15 M ⊕ within R = 0.25 R S , surrounded by a helium-dominated layer of 17 M ⊕ .This layer has Y ≈ 0.95 as determined self-consistently from the He-rich side of the SR18 demixing curves.Atop this layer is a helium gradient region extending from P ≈ 500 GPa to P = 200 GPa, outside of which the molecular envelope is convective and well mixed with a helium mass fraction of Y = 7%.Outside of its isothermal core, the model's temperature is assumed to be adiabatically stratified.Figure 1 displays the model's pressure and temperature conditions.For further details, we refer the reader to Mankovich & Fortney (2020) (particularly their Figures 2, 3, 10, and 16).

Ab Initio Simulations
Our simulations employ a combination of classical molecular dynamics (MD) for the ionic system coupled to a quantum mechanical description of the electrons within density functional theory (DFT; Hohenberg & Kohn 1964;Kohn & Sham 1965;Mermin 1965).We use the implementation in the Vienna Ab Initio Simulation Package (VASP; Kresse & Hafner 1993, 1994;Kresse & Furthmüller 1996, 1996).
In order to be consistent with calculations of the H-He demixing diagram of SR18, which was used in the calculation of the Saturn model (Mankovich & Fortney 2020), we employed the vdW-DF1 XC functional (Dion et al. 2004), projector augmented wave pseudopotentials PAW_PBE H 15Jun2001, PAW_PBE He 05Jan2001, and PAW_PBE O_h 06Feb2004 (Blöchl 1994), and a plane wave cutoff of 900 eV.Electronic temperatures enter the DFT formalism in the form of electronic Kohn-Sham orbitals, which are occupied according to the Fermi-Dirac distribution.We adjust the number of orbitals such that the total number of orbitals for each condition is 10%-20% above the last orbital for which the occupation is >5 × 10 −6 in the DFT-MD simulation.A Brillouin zone sampling with the Baldereschi mean value point (Baldereschi 1973) suffices for convergence of the DFT-MD simulations.
We control the ionic temperature with a Nosé-Hoover thermostat (Nosé 1984;Hoover 1985;Nosé 1991;Bylander & Kleinman 1992).The MD time step is between 0.2 and 0.3 fs in order to adequately resolve the vibrations of the hydrogen molecules.
We constructed 22 simulation boxes with 102-320 atoms, mostly hydrogen and helium, to reproduce the MF20 model's mass fractions for all radius coordinates.Table 1 lists the respective number of atoms in the simulation boxes.We chose oxygen to represent the heavy elements of the MF20 model.Figure 2 displays the particle compositions in terms of mass fractions in dependence of the radius.In each of the boxes, the initial configuration of the hydrogen atoms was molecular.We then ran DFT-MD simulations for at least 100,000 time steps.After about 20,000-100,000 time steps, depending on the condition, the simulation reached thermodynamic equilibrium.We then continued the simulation as long as necessary to yield converged values for the thermodynamic and transport properties, typically for another 100,000 time steps.In order to investigate the convergence of the self-diffusion coefficients for hydrogen, D H , and helium, D He , regarding the number of atoms, we performed additional calculations at R S = 0.90 and 0.55.Simulations without the oxygen atom reproduce the D H and D He of calculations that included the oxygen atom.We therefore found no apparent effect of the oxygen atom on D H and D He .Smaller simulations with 50, 100, 150, 200, and 250 hydrogen atoms and one to five helium atoms showed that 150 hydrogen atoms reproduce the result for D H obtained by a simulation with 300 hydrogen atoms within 15%.For helium, due to the larger uncertainties from the comparatively low number of atoms, a simulation with one helium atom reproduces the result of a simulation with six helium atoms.With the identical procedure for the viscosity η, the results of calculations with 100 hydrogen and two helium atoms reproduce our data for a total of 307 atoms within 12%.
Electronic transport properties were derived within the Kubo-Greenwood formalism (Holst et al. 2011;Kubo 1957;Greenwood 1958) using the HSE functional (Heyd et al. 2003(Heyd et al. , 2006)), see Section 4.3.Evaluation of their convergence behavior regarding the number of ionic snapshots, Brillouin zone sampling, number of atoms in the DFT-MDs, and number of electronic bands ensured that corresponding uncertainties are well below 5%.
Densities and temperatures were chosen such that the MF20 model's pressure and temperature conditions are reproduced  within 1%. Figure 1 shows the corresponding pressuretemperature conditions and the equilibrium pressures of our simulations, and Figure 2 the associated particle compositions of the simulation boxes.

Phase Transitions and Molecular Structure
We expect that several phase transitions occur in the H-He-O mixture between the outer layers and the core: The outermost layer will be molecular (H 2 ) and nonmetallic, followed by a continuous transition zone of hydrogen dissociation and simultaneous increase in conductivity deeper in the planet.At even lower radius coordinates, the system will be a metallic liquid.The He-rich layer at 0.35 R S 0.25 is expected to have a drastically lower conductivity than the metallic region as the insulator-to-metal transition in helium occurs at pressures of the order of tens of terapascals (Monserrat et al. 2014;Preising & Redmer 2020).
The density of states (DOS) is a useful quantity to characterize the electronic system.Within DFT-MD, the DOS is computed from a histogram of the Kohn-Sham eigenenergies in thermodynamic equilibrium (Preising & Redmer 2020).When the electronic DOS of a H-He-dominated system has values of zero in the vicinity of the chemical potential μ e , the system has a finite band gap and is therefore not metallic.
Figure 3 shows the DOS (upper panel) and the radial distribution function g H−H (r) (lower panel) for the radius coordinates R S = 0.80, 0.75, 0.37, and 0.35, obtained from calculations with the HSE functional.
From the outer layer to the deep interior, the H-He-O mixture undergoes several dissociation and ionization transitions.Between a radius coordinate of 0.80 and 0.70, the molecular hydrogen continuously transitions to a dissociated state, which is consistent with earlier findings regarding the location of the critical point of this transition (Lorenzen et al. 2010;Knudson et al. 2015Knudson et al. , 2018;;Mazzola et al. 2018).As can be seen in Figure 3, the band gap closes between 0.80 and 0.75 R S .The DOS is not exactly zero at said radius coordinates, which is typical of a bad metal with a pseudogap.The first g H−H (r) peak decreases concomitantly, indicating that hydrogen is dissociating in this region.The band gap keeps closing at lower radius coordinates down to the onset of the He-rich layer.Helium is an insulator under the conditions of the He-rich layer (Soubiran et al. 2012;Preising et al. 2018;Preising & Redmer 2020), in which the system's band gap opens again.
The correlations between the hydrogen atoms weaken upon dissociation; therefore, the peak of the H-H pair distribution g H−H (r) decreases with decreasing R S coordinates.In the Herich layer; however, the significantly higher number of helium atoms decreases the phase space available for hydrogen, and H 2 molecules form again.However, their bond length in the Herich layer (0.63 Å) is lower than for the P-T conditions in the outer layer (0.725 Å), again due to the reduced phase space.This finding is consistent with earlier studies of H-He mixtures (Lorenzen et al. 2010;Schöttler & Redmer 2018;Vorberger et al. 2007).

Thermodynamic Properties
We require additional EOS points to determine derivatives of pressure P, internal energy u, temperature T, and density ρ, which are needed for the calculation of thermodynamic material properties via Equations (1)-( 6).Therefore, we calculated four additional EOS points for each P-T point shown in Table 2; two on a corresponding isotherm with densities of ±4%, and two on a corresponding isochore with temperatures of ±4%.We consider the quantum effects of the ionic motion on the caloric EOS by calculating respective contributions via the vibrational power spectrum according to Berens et al. (1983) and Bethkenhagen et al. (2013).From those additional EOS points, we compute the partial derivatives required for the calculation of the thermal expansion coefficient α, the isothermal compressibility κ T , the heat capacities at constant volume c V and constant pressure c P , the sound velocity c s , and the Grüneisen parameter γ: The reported uncertainties are based on a statistical analysis of the time variations of P, u, and T during the DFT-MD simulations.The density ρ is an input parameter, and therefore not subject to uncertainties.We compare all quantities with values from a study of material properties under Jovian conditions by French et al. (2012).The following Figures 4-13 Figure 4 displays our results for the thermal expansion coefficient α, Equation (1), and the isothermal compressibility κ T , Equation (2).
Our values for the isothermal compressibility κ T decrease monotonically; the results for Saturn are systematically larger than those for Jupiter.The He-rich layer in the Saturn model results in a steeper decline of κ T compared to values at larger radii, which is a result of the dominating helium properties compared to those of hydrogen at higher radii.The dissociation of hydrogen molecules around R = 0.92 R J leads to a dip in the Jovian α values.Saturn is smaller and less dense than Jupiter.Therefore, dissociation takes place deeper in the planet and over a wider range of radius coordinates, starting from R < 0.80 R S .The thermal expansion coefficient displays a significant drop to negative values between R = 0.717 and 0.733 R S , which is due to the volume reduction of metallic hydrogen compared to molecular hydrogen, and therefore a consequence of pressure ionization.Aside from the region at which hydrogen dissociates, 0.783 R S 0.7, Saturn's values for α are systematically larger than those of Jupiter.Saturn's higher values for α and κ T account for the lower electron degeneracy in Saturn compared to the matter inside Jupiter.
Figure 5 shows our DFT-MD results for the specific heat capacities c V , Equation (3), and c P , Equation (4).Both curves for Saturn exhibit a pronounced peak around the dissociation transition of hydrogen (R ≈ 0.75 R S ), which accounts for the energy required to break molecular bonds.The corresponding slopes in Jupiter (French et al. 2012) show the same behavior at R ≈ 0.92 R J .At smaller radius coordinates, c V and c P decrease and exhibit a jump down to ∼6 J (gK) −1 at the helium shell boundary.The jump originates from helium's higher mass compared to that of hydrogen: helium has a lower specific heat capacity per mass unit than hydrogen, with a corresponding impact on the He-rich layer.
The sound velocity c s follows from κ T , c V , and c P according to Equation (5), resulting in the values displayed in Figure 6.Aside from variations near the radius coordinates at which hydrogen dissociates, our results for the sound velocity in Saturn increase with mass density.In the He-rich layer, the higher mass density of helium has a greater impact than the nearly constant ratio of c P /c V or the lower values of κ T , and c s decreases accordingly.Our results for the sound velocity in Saturn are systematically smaller than those for Jupiter (French et al. 2012), which reflects Saturn's lower density compared to Jupiter.The profiles for the sound velocity represent valuable information for the further development of seismology for both planets.They could help to better constrain their interior structure and composition (Gaulme et al. 2015).Similar to c s , the Grüneisen parameter γ follows from the quantities discussed above, namely, α, c V , and κ T according to Equation (6).Our results are shown in Figure 7.We additionally compare to values of an ideal Fermi gas in the deep interior and an ideal H-He mixture without vibrational excitations.
In the dissociation region of molecular hydrogen, the DFT-MD results for γ in both Saturn and Jupiter show a pronounced minimum.In Saturn, the minimum value is lower than γ = −0.09while for Jupiter a value of γ = 0.24 has been reported (French et al. 2012).Grüneisen parameters in Jupiter reproduce ideal gas values in limiting cases: the value of an ideal Fermi gas, γ = 2/3, at R J 0.5, and the value of a classical H 2 -He gas without vibrational excitations of hydrogen, γ(Y = 0.238) ≈ 0.487, in the outermost layer.Our results for Saturn right above the He-rich layer approach the ideal Fermi gas value.The He-rich layer itself exhibits slightly higher values.Our results for γ in the outer layer reproduce the value for an ideal mixture with the Saturn model's H-He concentration in the vibrational ground state, γ(Y = 0.07) ≈ 0.575, within the uncertainties.

Self-diffusion Coefficients
We calculate the self-diffusion coefficients D x for each species x = H, He, and O, via a Green-Kubo equation (Allen & Tildesley 1989): x where D x depends on the number of atoms N x for each species x, the velocity vector v i , and an integral over the time t.Equation ( 7) is evaluated from the ionic trajectories of the DFT-MD simulations.Our results are depicted in Figure 8.The diffusion coefficient of hydrogen decreases with density in Saturn's outermost layer.It reaches a minimum at the onset    of dissociation around R ≈ 0.75 R S and then increases again as more and more molecules dissociate and become more mobile.The higher helium concentration in the He-rich layer decreases the mobility of the hydrogen atoms, which becomes visible at 0.45 R S and cumulates in a significant decrease of D H to a nearly constant value of 0.06 mm 2 s −1 in the He-rich layer.The slope of D H in Saturn is consistent with its slope in Jupiter, where the dissociation of hydrogen takes place at R ≈ 0.92 R J and the influence of the He-rich layer is absent.Jupiter's higher temperatures result in higher particle velocities and therefore higher D H .
While the values for the diffusion coefficient of hydrogen follow the slope of D H in Jupiter at lower absolute values, the behavior of Saturn's D He is more complex: In the He-rich layer, D He ≈ D H ≈ 0.06 mm 2 s −1 .Due to the comparatively low number of two to 10 helium atoms, D He outside the He-rich layer is subject to large fluctuations, which can be seen in the slope shown in Figure 8.
This holds true for D O as well, as the simulation boxes only contain one to three oxygen atoms.Therefore, those values represent only rough estimates for the tracer diffusion coefficient of oxygen in the H-He mixture.Soubiran et al. (2017) performed DFT-MD simulations of a mixture of hydrogen, helium, and silicon dioxide.They calculated diffusion coefficients for hydrogen and helium at 5000 K, which corresponds to a pressure of ≈222 GPa and a radius coordinate of ≈0.54 in the MF20 model.We interpolated their diffusion coefficients for hydrogen and helium on the 5000 K isotherm to match the corresponding pressure of the MF20 model.Their lower hydrogen diffusivity originates from a greater amount of helium and heavy elements in their simulations.Their result for helium agrees with ours.Wilson (2015) calculated diffusion coefficients for hydrogen, helium, carbon, silicon, and iron with DFT-MD.However, their lowest temperature under consideration was 10,000 K. Therefore, we do not compare their results to our data which range from around 2000 K to around 9000 K.
From the self-diffusion coefficients D H and D He , we can estimate the interdiffusion coefficient D H,He via Darken's relation (Krishna & van Baten 2005;Liu et al. 2013).Together with the thermal diffusivity and viscosity, D H,He can influence the convection behavior of the fluid, which might have a semiconvective or diffusive portion in certain parts of the planet (Stevenson & Salpeter 1977a;Rosenblum et al. 2011;Leconte & Chabrier 2012;Nettelmann et al. 2015;French & Nettelmann 2019).

Viscosity
We evaluate the shear viscosity η via a Green-Kubo formula employing autocorrelation functions of off-diagonal elements of the pressure tensor p ij in the framework of linear response theory (Kubo 1957;Allen & Tildesley 1989;Alfè & Gillan 1998;French & Nettelmann 2019).We extract the pressure tensors from our DFT-MD simulations in thermodynamic equilibrium.
Aside from p ij , the dynamic shear viscosity η depends on the volume V of the simulation box, the Boltzmann constant k B , the temperature T, and a sum over the directions xy, yz, and zx.Our conservative estimate of uncertainty for our results for the viscosities is 15%, except for the highest radius coordinate, where it is 50%.Figure 9 displays our results for the dynamic shear viscosity.
The DFT-MD results for the dynamic shear viscosity η in Jupiter increase with decreasing radius.Our results for Saturn follow the same slope and are systematically lower.In Saturn's He-rich layer, however, the density increase leads to values greater than those in Jupiter.
In a DFT-MD study of hydrogen, helium, and silicon dioxide under planetary interior conditions, Soubiran et al. (2017) also calculated the viscosity of their mixture at 5000 K.As mentioned in the discussion of Figure 8, we interpolated their data to a pressure of ≈222 GPa, which corresponds to a radius coordinate of ≈0.54 in the MF20 model.Their predicted viscosity is roughly 25% higher than our result at this radius coordinate, which is consistent with the greater amount of heavy elements in their simulations.
We additionally compare values extrapolated from lowdensity experimental data on pure hydrogen (Stiel & Thodos 1963).The extrapolated data feature significantly lower pressures than this study.We therefore compare the data at their highest pressure ≈130 kPa and the temperatures of the MF20 model between 0.522 and 0.98 R S .Our result at 0.95 R S agrees well with the extrapolated data.At lower radius coordinates the mixture becomes more correlated, resulting in increased viscosities, which is why our results differ from their low-density values.
Figure 10 shows our results for the kinematic shear viscosity ν = η/ρ.DFT-MD results for the kinematic shear viscosity in Jupiter are approximately constant and increase only slightly with radius.Our results for Saturn's He-rich layer exhibit a significantly higher ν than the predictions for Jupiter at the same relative radius coordinates.Our predictions for the dynamic shear viscosity in the layer of the dissociated fluid in Saturn coincide with Jovian values in the region directly above the He-rich layer.Their increase with radius is more pronounced than the predictions for Jupiter.Aside from our results for the He-rich layer, ν is almost constant within the  (French et al. 2012) in orange, compared to viscosity data for pure hydrogen from fitting to experimental data in blue (Stiel & Thodos 1963), and interpolated data from DFT-MD results for a similar mixture with a higher content of heavy elements (Soubiran et al. 2017) as light blue diamond.
uncertainties.Within the region that is most relevant for the operation of a dynamo, the approximation of a constant ν, which is a common assumption in dynamo modeling (Wicht & Tilgner 2010), is justified for Saturn according to our calculations.
Together with the higher density, the higher viscosity in the He-rich layer slows convection and might contribute to a stabilization of the layer within the planet.

Thermal Conductivity
We calculate the thermal conductivity as the sum of ionic and electronic contributions: Here, the generalized heat currents ¢ ( ) J t q and the Green-Kubo formula (Equation ( 10)) are evaluated as described by French (2019).Our estimate for the uncertainties in λ i is 15%.
For the electronic contribution, we calculate the frequencydependent Onsager coefficients L mn (ω) from (Holst et al. 2011): which contain the frequency ω, the charge q = −e, the electron mass m e , a sum over the Bloch states |kμ〉 with the occupation f kμ and the eigenenergies E kμ , the enthalpy per electron h e , and the reduced Planck's constant ÿ.
The matrix elements n 6 that take into account the nonlocal contributions from the Hartree-Fock exchange and from PAW potentials as implemented in the optical routines of VASP (Gajdoš et al. 2006;French & Redmer 2017).
We extract 30 ionic configurations from each of the DFT-MD simulations and then perform static DFT calculations with a Monkhorst-Pack (Monkhorst & Pack 1976) k-point sampling of 2 × 2 × 2. We use the HSE (Heyd et al. 2003(Heyd et al. , 2006) ) exchange-correlation functional with a Hartree-Fock screening parameter of 0.2 Å −1 for a more accurate description of the band gap, and by extension, all electronic transport properties.Our estimated uncertainties for λ e vary from 5%-50%, depending on the shape of the function we evaluate according to Equation (11).
The relation between the thermal diffusivity κ and the thermal conductivity is shows our results for the thermal conductivity λ and the thermal diffusivity κ.
In the molecular fluid of the outer layer, ionic contributions to the thermal conductivity dominate in our calculations.Upon the dissociation of molecular hydrogen, heat conduction from the electronic system becomes more pronounced and finally dominates in the dissociated fluid until the onset of the He-rich layer.Here, hydrogen molecules form in an atomic helium environment, which leads to only a few conducting electrons, and thus, a dominating λ i .The thermal conductivity in Jupiter is significantly larger than in Saturn for R 0.90 R Obj due to higher temperatures and pressures in Jupiter's interior.
The slopes for the thermal diffusivities κ from DFT-MD simulations are similar for Jupiter and Saturn.In the He-rich layer, however, κ is significantly smaller in Saturn than the thermal diffusivity in Jupiter due to the direct proportionality with the He layer's lower thermal conductivity λ, higher density, and lower heat capacity c P at constant pressure (see Figure 4).
The He-rich layer exhibits a higher density, higher viscosity (see Section 4.3.2),and lower thermal conductivity than the layers above it.The corresponding reduction in heat flow from the hot core to the H-He envelope might have a significant  impact on the cooling of the planet.To the best of our knowledge, this has never been accounted for in evolution studies for Saturn.

Static Electrical Conductivity
The static electrical (DC) conductivity σ follows from the static limit of the Onsager coefficient L 11 (ω): We estimate the uncertainties, depending on the low-frequency behavior of L 11 (ω), to be 5% in the metallic region between 0.733 and 0.37 R S , 20% at 0.75 R S , and 50% at all other radius coordinates.Figure 12 displays our results for σ.
Our results for σ have a slope similar to that of the results for Jupiter.The DC conductivity increases from the outer layer with increasing density.Due to a finer grid regarding radius coordinates, our results resolve a noticeable change in slope in the region of hydrogen dissociation.Our results predict an increase of σ in Saturn almost up to the values calculated for Jupiter (French et al. 2012).The electronically insulating helium in Saturn's He-rich layer results in a significant drop in σ by several orders of magnitude.The DC conductivity in the He-rich layer itself increases with density.
Results for the electronic contribution to the thermal conductivity (see Section 4.3.3)and the DC conductivity demonstrate that the He-rich environment deep in the planet impairs electronic transport, i.e., both electronic heat and charge conduction.Therefore, even if the He-rich layer was vigorously convective, it would barely contribute to the generation of Saturn's planetary magnetic field.Liu et al. (2008Liu et al. ( , 2014) ) calculated interior models for Jupiter and Saturn based on the EOS of Saumon et al. (Saumon et al. 1995) and then inferred DC conductivities from a semiconductor model for hydrogen.Their earlier model (Liu et al. 2008) predicts DC conductivities in the molecular region that are systematically lower than our DFT-MD results but reproduces the slope of our results.Their later model (Liu et al. 2014) extends to smaller radius coordinates and features two distinct changes in slope; one around 0.85 R S , which results from the onset of hydrogen metallization, and another around 0.6 R S marking the full hydrogen metallization.Our ab initio simulations predict a comparable slope in the molecular region, and steeper slopes in the region of gradual hydrogen metallization as well as the fully dissociated regime.The main difference between our results and the Liu model is the location of the onset of hydrogen metallization at ≈0.775 R S (Liu model: ≈0.85 R S ) and the location of full metallization at ≈0.7 R S (Liu model: ≈0.615R S ).
The magnetic diffusivity β follows from σ according to b sm = ( ) 1 , 1 5 0 where μ 0 is the vacuum permittivity.Table 3 lists our results.Notice that our DFT-MD results for σ(R) (as well as other material properties reported in this work) along the P-T path for both planets are consistent with the corresponding EOS data.Other conductivity models in this parameter space are only valid for limiting cases.For instance, a multicomponent plasma model for the calculation of the composition and conductivity in the atmosphere of gas giants, in particular of hot Jupiters, has been developed recently (Kumar et al. 2021).
The widely used Lee-More conductivity model (Lee & More 1984) is based on the relaxation time approximation and can only be applied to the fully ionized plasma in the deeper interior.Interpolation formulae for the electrical conductivity of strongly coupled, fully ionized plasmas based on correlation function approaches were derived for the same region (Ichimaru & Tanaka 1985;Röpke & Redmer 1989).
The Lorenz number L is a relation between electrical and thermal conductivities (Lorenz 1872): In a completely degenerate electron gas, where only electrons are responsible for electrical and thermal conduction, like in metals and dense plasmas, L = π 2 /3 according to the Wiedemann-Franz law (Franz & Wiedemann 1853).When ions (or atoms and molecules) substantially contribute to the thermal conductivity, in weakly ionized gases or fluids, L can be orders of magnitude larger.Figure 13 displays our results for L.
Our results for L in Saturn's He-rich layer are at least an order of magnitude larger than π 2 /3 due to the contributions of  nonionized helium in the He-rich layer.In regions with smaller helium content, our values for L are closer to those in Jupiter.They reproduce the Wiedemann-Franz law within the uncertainties between the He-rich layer and the onset of hydrogen dissociation.However, L has large positive values in the lower He-rich and outer molecular layer, which is typical for the behavior of a semiconductor or a bad metal (Holst et al. 2011).
In light of the drastically lower electrical conductivity we find between R = 0.25 and 0.35 R S as a result of helium sedimentation, it is clear that the processes in the region of Saturn's dynamo generation depend intimately on the H-He phase diagram.We encourage future magnetohydrodynamic studies to consider the effects of a deep insulating layer in which properties of helium dominate the properties of the H-He mixture due to the higher helium concentration.In particular, the apparent contradiction of Saturn's magnetic field with Cowling's theorem (axisymmetric with the rotation axis) has to be resolved by magnetohydrodynamic modeling, which can now be based on our ab initio database for material properties.Table 3 compiles our results for transport properties of matter in Saturn's interior.

Conclusions
We calculated thermodynamic and transport properties along the P-T path of a recent evolutionary Saturn model of MF20 in a range between 0.25 and 0.95 R S .The MF20 model is based on ab initio data for the EOS (Becker et al. 2014) and the demixing diagram of hydrogen and helium (Schöttler & Redmer 2018).Both the EOS and the demixing diagram employed MD simulations based on a self-consistent quantum description of the electronic structure within DFT-MD.
Our findings significantly improve and enhance our knowledge about material properties (thermodynamic and transport) in Saturn.They provide robust results through the consistent use of DFT-MD throughout the generation of data for the MF20 model as well as for our calculation of material properties.
Compared to earlier results for Jupiter (French et al. 2012), the lower temperatures in the MF20 model lead to more pronounced implications of the insulator-to-metal phase transition on material and transport properties in the hydrogen-dominated outer layers of Saturn.We predict a significant impact of the metallization of hydrogen and the MF20 model's He-rich layer, especially on the thermal and electrical conductivity, as well as considerable consequences regarding the cooling behavior of Saturn.
Our results can serve as input for future models of Saturn's interior dynamics, especially hydrodynamic simulations for its magnetic field.Furthermore, they will help infer constraints for the interior structure and composition from future Saturn seismology data.

Figure 1 .
Figure 1.Temperature T over pressure P of the MF20 Saturn model (alternating red and yellow lines), our simulations (black circles), calculations of the DC conductivity by Knudson et al. (2018) performed with DFT-MD and the vdW-DF1 functional (orange circles), dissociation of hydrogen from DFT-MD simulations with the vdW-DF1 functional by Knudson et al. (2015) (orange line) and from quantum Monte Carlo (QMC) calculations by Mazzola et al. (2018) (blue line).Different colors and line styles of the MF20 model denote (from lowest to highest pressures): layer with molecular, dissociating, and fully dissociated H 2 , He-rich layer, and isothermal core, respectively.

Figure 3 .
Figure 3. DOS (upper panel) over energy and radial distribution function of hydrogen g H−H (r) over the interatomic distance r (lower panel) at color-coded radius coordinates.

Figure 5 .
Figure5.Specific heat capacities at constant volume c V (dashed lines) and constant pressure c P (solid lines) in Saturn (black) and Jupiter(French et al. 2012) (orange).

Figure 7 .
Figure 7. Grüneisen parameter γ in Saturn (black) and Jupiter (French et al. 2012) (orange).Blue-dashed lines denote values of an ideal Fermi gas (γ = 2/3) and a classical ideal gas in the vibrational ground state with the H-He concentration predicted by the Saturn model (γ = 0.575).The red-dashed line (γ = 0.487) indicates values for a H-He ideal gas with Jovian composition.

Figure 8 .
Figure8.Self-diffusion coefficients D from DFT-MD for Saturn (solid lines) and Jupiter(French et al. 2012) (orange-dashed lines).Upper panel: D of hydrogen (black), helium (orange), and oxygen (red).The two lower panels show D H and D He , respectively.Blue diamonds denote interpolated data from DFT-MD results for a similar mixture with a higher content of heavy elements(Soubiran et al. 2017).

Figure 9 .
Figure9.Dynamic shear viscosity η from DFT-MD for Saturn in black and Jupiter(French et al. 2012) in orange, compared to viscosity data for pure hydrogen from fitting to experimental data in blue(Stiel & Thodos 1963), and interpolated data from DFT-MD results for a similar mixture with a higher content of heavy elements(Soubiran et al. 2017) as light blue diamond.

Figure 10 .
Figure10.Kinematic shear viscosity ν from DFT-MD for Saturn in black and Jupiter(French et al. 2012) in orange.Figure 11.Upper panel: results for the thermal conductivity λ in Saturn (black) and Jupiter (French et al. 2012) (orange).Saturn: blue diamonds connected by dashed-dotted lines represent electronic contributions to λ and light blue circles connected by dashed lines represent ionic contributions to λ. Lower panel: thermal diffusivity κ in Saturn and Jupiter.

Figure 11 .
Figure10.Kinematic shear viscosity ν from DFT-MD for Saturn in black and Jupiter(French et al. 2012) in orange.Figure 11.Upper panel: results for the thermal conductivity λ in Saturn (black) and Jupiter (French et al. 2012) (orange).Saturn: blue diamonds connected by dashed-dotted lines represent electronic contributions to λ and light blue circles connected by dashed lines represent ionic contributions to λ. Lower panel: thermal diffusivity κ in Saturn and Jupiter.

Figure 12 .
Figure12.DFT-MD results for the DC conductivity σ for Saturn (black) and Jupiter(French et al. 2012) (orange).Data shown in blue and light blue represent the results of earlier studies, obtained from calculations in the chemical picture, byLiu et al. (2008Liu et al. ( , 2014)).

Figure 13 .
Figure 13.Lorenz number L for Saturn (black) and Jupiter (French et al. 2012) (orange) from DFT-MD compared to the Wiedemann-Franz law (gray-dashed line).The inset shows the results for the outer layer.

Table 1
Simulation Box CompositionsFigure 2. Mass fraction over radius coordinate.Solid lines represent the MF20 model's mass fractions, open circles those of our simulations.The colors denote the elements (hydrogen in black, helium in orange, and oxygen in red).

Table 2
Thermodynamic Material Properties in Saturn

Table 3
Transport Properties in Saturn