Heat and charge transport in H2O at ice-giant conditions from ab initio molecular dynamics simulations

The impact of the inner structure and thermal history of planets on their observable features, such as luminosity or magnetic field, crucially depends on the poorly known heat and charge transport properties of their internal layers. The thermal and electric conductivities of different phases of water (liquid, solid, and super-ionic) occurring in the interior of ice giant planets, such as Uranus or Neptune, are evaluated from equilibrium ab initio molecular dynamics, leveraging recent progresses in the theory and data analysis of transport in extended systems. The implications of our findings on the evolution models of the ice giants are briefly discussed.


SUPPLEMENTARY NOTE 1. AB INITIO SIMULATIONS
Ab initio molecular dynamics simulations were run using the plane-wave pseudo-potential method and the Car-Parrinello (CP) Lagrangian formalism [1], as implemented in the cp.x component of the Quantum ESPRESSO suite of computer codes [2,3]. We adopted norm-conserving pseudopotentials [4] from Ref. 5, which are accurate also at the present pT-conditions [6], and the PBE energy functional [7]. We employed a plane-wave energy cutoff of E cut = 110 Ry for wave-functions and four times as much for the charge density. The fictitious electronic mass was set to 25 electronic masses and the integration time step to ≈ 0.0484 fs. These choices ensure adequate adiabatic decoupling of the electron fictitious dynamics from the ionic motion and conservation of the extended total energy in the relevant pT range considered here. All our simulations have been performed neglecting nuclear quantum effects (NQEs), as these are expected to become negligible as temperature increases. Even though NQEs may be counter-enhanced by high pressure [8], their account is beyond the scope of the present work.
Flux time series were collected from N V E simulations for at least 60 ps, after an initial N P T equilibration of a few ps at the target pT conditions, followed by ∼ 4 ps of further N V E equilibration. N P T simulations were run with a chain of three Nosé-Hoover thermostats [9] with a frequency of 60 THz, while pressure was controlled via a Parrinello-Rahman barostat [10]. In order to minimize the effects of Pulay stresses, the kinetic energy functional was modified as: [11] with E 0 = 100 Ry, A = 170 Ry, and σ = 15 Ry.
Our simulations were run on samples of 128 atoms for ice-X, SI-BCC, and PDL water, and of 108 atoms for SI-FCC. In order to estimate the magnitude of finite-size effects arising in the AIMD simulation of the thermal conductivity, we have performed extensive MD simulations of PDL water samples of 128, 256, 512, and 1024 H 2 O units, using a deep-neural-network potential trained on our ab initio trajectory [12], at the same conditions (T ≈ 1970 K and p ≈ 33 GPa). These results indicate that finite-size corrections to the thermal conductivity computed for our sample of 128 H 2 O units should be smaller than 5%, in qualitative agreement with systematic tests performed in Ref. 13.

SUPPLEMENTARY NOTE 2. FLUX TIME SERIES
The energy and charge fluxes, J E (t) and J Z (t), were sampled every 30 or 40 time steps of the original CP simulation (∆ t = 1.45 fs or 1.94 fs), depending on the specific simulation. This sampling rate is high enough to make the flux power spectra vanish at the corresponding Nyqvist frequency, f Ny = 1/(2∆ t ), thus avoiding any aliasing effects. We verified that the energy band gap is far larger than k B T all along our simulations, thus excluding any possible non adiabatic effects in heat and charge transport at the pT conditions of interest here. The time series of the energy gap occurring in two different SI phases of water are reported in Supplementary Figure 1.
The energy flux was estimated using the expression derived in Ref. 14. The charge flux was estimated using two different expressions: i) the charge flux, J Z , estimated from the classical expression derived from the formal oxidation numbers q H = +1 and q O = −2, according to the topological theory of adiabatic charge transport introduced in Ref. 15, and ii) the standard quantum-mechanical one, J Z , that is the sum of the adiabatic electron-charge flux, J el , and is the number of H2O stoichiometric units employed in the simulation. Temperature and pressures are reported together with their standard deviation. The last two simulations concern the defected structure, where one H2O was removed. The temperatures and pressures are reported with their standard deviations. The uncertainty on the diffusivity is estimated via a 6 block analysis with ≈ 10-ps segments.
the classical current of the atomic cores (Z H = +1 and Z O = +6). The two expressions read: The adiabatic electron charge current, J el , can be computed as: where the sum runs over the occupied Kohn-Sham states, {φ v }, and the orbitals are the projections over the empty-state manifold of the action of the position operator over the v-th occupied orbital, and of its adiabatic time derivative [3]. The operatorsP v andP c = 1 −P v are projectors over the occupied-and empty-states manifolds, respectively. Equations (4) are well defined in periodic boundary conditions and have been computed from standard density-functional perturbation theory [16].

SUPPLEMENTARY NOTE 3. CEPSTRAL ANALYSIS
Transport coefficients were estimated by cepstral analysis of the flux time series, as explained in Refs. 17 and 18, using the thermocepstrum code [19]. This technique allows one to obtain accurate transport coefficients, as well as a quantitative estimate of their statistical accuracy, depending on two parameters: an effective Nyqvist frequency, f * , used to limit the analysis to a properly defined low-frequency portion of the spectrum, and the number P * of inverse Fourier ("cepstral ") coefficients of the logarithm of the spectrum in the low-frequency region thus defined. The estimated conductivities depend very little on f * , whereas an optimal value of P * can be estimated in most cases by statistical model-selection techniques, as explained in Ref. 17. In our applications P * is determined via the so-called Akaike Information Criterion [17,20] Supplementary Figure 2 displays the sample power spectra of J el , J Z and J Z , scaled in units such that the zerofrequency value is the electrical conductivity in S/cm. The faint lines are obtained by applying a moving average filter [21] with a window of 0.1 THz to the raw periodograms; the solid smooth lines are obtained by applying a cepstral filter [17] with the parameters reported in Supplementary Table 1 conductivity, even if their time series values differ, as is clear from the inset, in agreement with Ref. [15], where this finding was proved to be the consequence of gauge-invariance of transport coefficients and topological quantization of adiabatic charge transport in electronically gapped materials.
As stated in the main text, a 3-variate analysis on thermal conductivity adopting the adiabatic electron flux J el as an additional flux reduces the statistical uncertainty without affecting the estimate of the thermal conductivity, as shown in Supplementary Figure 3 (green). The 3-variate analysis produces a flatter spectrum featuring a reduced total power, corresponding to a smaller number of cepstral coefficients, P * , thus ensuring a smaller statistical error.
A comprehensive summary of our results, along with the values of f * and P * used for cepstral analysis, is given in Supplementary Table 1.

SUPPLEMENTARY NOTE 4. DEFECTED STRUCTURE
We also performed two further simulations-one in solid ice X and one in the SI phase with BCC oxygen structurefor defected structures where two H and one O atoms were removed, to preserve the stoichiometry of the ideal systems, resulting in a defect relative concentration of 1/128 ≈ 0.008. In the SI phase, vacancy hopping in the oxygen lattice is observed in the form of jumps of one or more than one O atoms within the same hopping event. From the mean square displacement of oxygen atoms with respect to their center of mass, we estimated the O self-diffusion coefficient as D O ≈ (8.5 ± 0.6) × 10 −3Å2 /ps. The linear growth of the oxygen mean square displacement, together with the related statistical uncertainty from a 6-block analysis of 10ps segments, is reported in Supplementary Figure 4. On the contrary, in the defected ice X structure below the SI transition we observed no oxygen self-diffusion.
Even if such defect concentration is far larger than what it is expected to be at the pT-conditions of Uranus and Neptune from extrapolation of experimental data [22], it is still insufficient to strongly affect the thermal conductivity, when compared to the ideal-structure value. In fact, we obtained values for κ which are compatible with those of the ideal structures, namely κ ≈ 17 and κ ≈ 12 W/(Km) for the solid and the SI phases, respectively.