Thermal conductivity of dissociating water—an ab initio study

The thermal conductivity of partially dissociated and ionised water is calculated in a large-scale study using density functional theory (DFT)-based molecular dynamics (MD) simulations. In doing so, the required heat current of the nuclei is calculated by mapping the effective particle interactions from the DFT-MD simulations onto classical pair potentials. It is demonstrated that experimental and theoretical thermal conductivity data for liquid heavy water and for ice VII are well reproduced with this efficient procedure. Moreover, the approach also allows for an illustrative interpretation of the characteristics of the thermal conductivity in the dense chemically reacting fluid. The thermodynamic conditions investigated here range from densities between 0.2 and 6 g cm−3 and temperatures between 600 and 50 000 K, which includes states highly relevant for understanding the interiors of water-rich planets like Uranus and Neptune and exoplanets of similar composition.


Introduction
Gaining accurate knowlegde of the thermodynamic and transport properties of water and water-rich molecular mixtures at high pressures and high temperatures is key to understanding the interior structure, thermal evolution, and magnetic fields of icy giant planets [1][2][3]. For example, a major unsolved puzzle is the existence of a thermal convection barrier deep inside the ice giant Uranus, which may explain its unusually low luminosity [4]. Furthermore, it is unknown to which degree the icy cores of Jupiter or Saturn are eroded and dissolved into their envelopes [5]. But quantitative planetary evolution models that include both convective and diffusive heat transport [6,4] require reliable thermal conductivity data, which are not yet available at the relevant conditions. Direct measurements of the thermal conductivity of materials at extreme pressures and temperatures are highly challenging since undesirable influences of convection, radiation, and heat flux through the confining device have to be inhibited [7]. Much of the progress has therefore been achieved on good conducting metals [8,9]. For example, the thermal conductivity of iron was recently measured up to pressures of 140 GPa and temperatures of 3000 K using laser-heated diamond anvil cells [10]. The most extreme conditions at which thermal conductivity of water was measured, so far, were only 3.5 GPa at temperatures less than 700 K [11].
This work describes an efficient procedure to calculate the thermal conductivity with ab initio computer simulation techniques. The method is then applied to generate a large data base for the thermal conductivity of water that spans a density range between 0.2 and 6 g cm −3 and temperatures between 600 and 50000 K.
The latter is tremendously more difficult to derive because a rigorous calculation of the nuclear heat current is complicated and has been achieved yet only for one-component fluids in ground-state DFT-MD simulations [17,18]. The reason is that the underlying Born-Oppenheimer approximation decouples the dynamics of electrons and ions. The electron-nucleus interactions then need to be properly reintroduced into the definition of the microscopic nuclear heat current so that it contains the dynamical information relevant for the time scale of nuclear motion. Furthermore, the numerical formalism inherently suffers from rather poor convergence behaviour of the time correlation functions in the Green-Kubo formula, which is another reason for the absence of any large-scale ab initio studies on nuclear thermal conductivity of fluid systems to date.

Theoretical method and simulation procedure
To calculate the nuclear heat current in the FT-DFT-MD simulations with an efficient but yet sufficiently accurate procedure, we map all effective interactions between the nuclei on radial pair potentials. This allows us to employ an established expression for the heat current in the Green-Kubo formula for the thermal conductivity λ nuc [19][20][21]: In general, the heat current t J q ¢ ( ) can be written as sum of a kinetic part t J q,kin ¢ ( ) and an interaction part t J q,int ¢ ( ). The kinetic part is easily calculated from the particle velocities v i (t) via: where m i is the mass of particle i N , is the total particle number, and k T 5 2 B is the average kinetic enthalpy per particle. The interaction part for pairwise interacting particles can be written as: where F ij (t) and r ij (t) are force and distance between particles i and j, respectively. The interaction energy of a particle reads: The partial interaction enthalpy of particle i of species a is taken as: The above expression neglects nonideal mixing contributions to h F a , but simplifies the calculations [22,23]. Nonideal mixing effects on equation of state quantities were found to be very small in calculations of waterammonia-methane mixtures under similar conditions [24]. Nevertheless, avoiding this approximation is possible, and a comparison with results from an exact expression for λ nuc was made, see the appendix. The very good agreement shows the suitability of the approach described here for dense water plasmas.
The following ansatz for the potential between the nuclear species a and b is made: where r=r ij is the radius and E l l ab ab ab . The exponential term [25] in equation (6) serves to smoothly truncate the potential at half of the simulation box length L. It allows for the evaluation of all expressions that enter equation (3) within the minimum image convention [21] and for systematic convergence tests with the system size. Extensive tests have shown that there are universally suitable choices possible for some of the parameters, i.e. 0.1, 0.5 ab ab a g = = , and 4 ab d = -. To allow the potential(6) to reach its converged shape, which depends on density ñ and temperatureT, the cutoff parameter C ab needs to be set typically to numbers between 10 and 20. Potential and thermal conductivity both converge very similarly with C ab .
After the completion of an FT-DFT-MD simulation run, the coefficients c l ab are fitted by matching the force fields r V r F ij ij ij ij = - ( ) ( ) to the forces on the particles at each 20th (or 50th for very long simulations) time step, which fully determines the potentials (6) for that particular simulation run. This mathematical task can be formulated as a linear least-squares optimisation problem, and its exact solution is determined with means of linear algebra.
After the pair potentials are parametrised, equation (1) is evaluated using the particle trajectories of the same FT-DFT-MD simulation run. This is a very important point to take note of because it ensures that the correlation function of the heat current in equation (1) is evaluated using the ensemble (trajectories) á¼ñ generated in the FT-DFT-MD simulation 1 . The whole procedure is labelled here as DFT-MD+GKFF method 2 .
The formalism described above is now applied to various FT-DFT-MD simulations of warm dense water performed with the VASP code [26][27][28] and the PBE exchange correlation functional [29]. Most of the simulation parameters were set to values well-proven for calculating equation of state and diffusion coefficients [30][31][32][33]. The major modification here is the large extension in simulation time, so that now between 50000 and 250000 time steps (typical length of a time step: 0.3-0.5 fs) are simulated. This large increase is necessary to converge the time correlation function in equation (1) and reach a statistical uncertainty for the λ nuc data of about 10%-20%.

Results and discussion
At first, we perform a special benchmark simulation in order to make a direct comparison with the sole example for rigorous calculation of the thermal conductivity of water made with DFT-MD so far [17]. In that work, Marcolongo et al derived a special variant of density functional perturbation theory [34,35] that can describe the heat flow contribution of bound electrons via the energetic polarisation induced by nuclei in motion. Combining this formalism with a ground-state DFT-MD run enabled them to calculate the nuclear thermal conductivity for electronically insulating systems, which is adequate for liquid water. However, this elaborate procedure enhances the computing time of a DFT-MD simulation run by significant factor, similarly as occurs in analogous calculations of ionic conductivities [36,37].
Our benchmark simulation is performed with liquid heavy water (64 molecules) at the same density of 1.11 g cm −3 and same temperature of 385 K as reported by Marcolongo et al [17]. The resulting integrated thermal conductivity from the DFT-MD+GKFF approach is displayed in the upper panel of figure 1. After about 100fs the short-time correlations have decayed and the final value of 0.7±0.1W Km −1 is reached. This number is in excellent agreement with the result from Marcolongo et al [17] of 0.74±0.12W Km −1 and also with the experimental value of about 0.68W Km −1 [38].
As second example we compute the thermal conductivity of ice VII at room temperature and at 4 densities using proton-ordered 3 crystal structures of 54 and 128 molecules. The experimental thermal conductivity of ice VII, which was measured up to pressures of 22GPa [39,40], is reproduced very well, see lower panel of figure 1. The same plot also shows that the DFT-MD+GKFF approach produces converged values for the thermal conductivity already with 54 molecules in the simulation box. Success at this crucial convergence test indicates that the mapping of the long-ranging forces between charged ice nuclei on short-range pair potentials is adequate for calculating the heat current. This can be understood through the structure of equation (3), which requires the force fields only to reproduce well the fluctuations of virial and energy terms around average values. These fluctuations are mainly determined by the effective interactions on short-range (particle collisions). Accurate reproduction of the total average pressure or Ewald energy, which are more influenced by the longrange character of the Coulomb forces, is not required.
In general, the interatomic potentials(6) look unspectacular for the molecular and fully dissociated fluid states, i.e. they show characteristic repulsive and/or attractive properties. Their behaviour in the partially 1 An obvious alternative, i.e. performing a subsequent classical MD simulation with the derived potentials, would lead to a different ensemble because the radial force fields cannot reproduce the DFT forces exactly. This would significantly worsen the results because such an error would then also propagate into the other relevant dynamic quantities (positions, velocities). 2 GKFF stands for Green-Kubo (with) force fields. The prefix FT is omitted for simplification, albeit all simulations were made with thermally excited electrons. 3 Proton order was assumed for simplicity and to enable unambiguous convergence checks with increasing particle number. The thermal conductivity of solids is mainly due to the acoustic phonons, so that the influence of proton disorder is probably small. dissociated region is worth explicit discussion, see figure 2 (upper panel), which shows an example at a density of 1 g cm −3 and a temperature of 6000 K. This state represents a reacting dense oxyhydrogen gas mixture, where very short-lived H 2 and O 2 molecules are present as indicated by the pair correlation functions in figure 2 (lower panel). The most prominent attractive well occurs in V OH at the typical length of an intramolecular O-H bond. In addition, both V HH and V OO have shallow valleys at the distances of the diatomic bond lengths of H 2 and O 2 molecules, respectively. These valleys can be reached only by particles that overcome significant potential barriers even higher than the valley depths. All this illustrates the endothermic character of the dissociation reaction of water molecules in this supercritical fluid state.
The main effort of this work is directed towards a general understanding of the thermal conductivity of fluid water across large density and temperature ranges. For such purpose, we first evaluate equation (1) by plugging in the kinetic part of the heat current, equation (2), only. The respective results, λ kin , are shown as dots connected with dashed lines in figure 3. The values rise systematically with temperature but they are relatively independent of the density, especially below 3000 K, where water is largely molecular [30,32]. This behaviour is somewhat surprising and it suggests that the change in correlations between the water molecules with increasing density does not lead to drastic effects on distribution of their velocities. The averaged values of λ kin for the  temperatures 600, 1000, and 2000 K are in good concordance with the known low-density limiting values of the thermal conductivity for molecular water gas [42,43].
Several interesting changes occur when λ nuc is calculated with the complete heat current that includes both kinetic and interaction part. For molecular water (below T3000 K), systematic deviations from λ kin arise as the thermal conductivity acquires a pronounced increase with the density. At densities larger than 1 g cm −3 , this leads to a merging of the λ nuc curves for 600 K T 3000 K with no significant effect of temperature. An absent or weak temperature dependence of λ nuc is common for many liquids [44]. Liquid water at few 100°Celsius is a known exception because its thermal conductivity increases, and this can be understood quantitatively through an excess contribution from the dissociation of hydrogen bonds [45,46]. It is still debated if and at which density and temperature this rise in thermal conductivity eventually levels off [11,46,47]. Our results indicate that this happens somewhat below 600 K, which was observed earlier in classical MD simulations [47]. The only experiment made at these conditions, however, reported of a T nuc l~behaviour that persists up to more than 600 K [11], but this unusual scaling law 4 has not yet been independently confirmed. On the other hand, the accuracy of FT-DFT-MD in the description of hydrogen bonds might not yet be sufficient to fully understand the temperature dependence of the thermal conductivity in liquid water, either due to the exchange correlation functional [48] or the classical treatment of the ionic motion. Despite the very good agreement with the measurements in the examples shown in figure 1, our four 600 K data points for λ nuc are lower by 20%-30% than recommended reference values [42]. Still, the density dependence of our thermal conductivites above 1 g cm −3 is very well in line with that measured by Abramsonet al [11], see figure 3. This indicates that the increase in strength of the repulsive molecular interactions under compression is captured well in the DFT-MD+GKFF approach.
A prominent increase in λ nuc emerges at temperatures greater than 3000 K, where water dissociates into atoms or smaller molecules. Figure 3 shows this especially at 6000 and 10000 K, where λ nuc does not reach the limiting values of λ kin at the lowest density considered here. The reason is the large contribution due to the dissociation reactions of the molecules that enters the interaction part, t J q,int ¢ ( ), through the shape of the potentials(6), as described above on the basis of figure 2. This contribution is analoguous to that from hydrogen-bond dissociation in the liquid, but much larger in number because the intramolecular bonds have greater dissociation energies. Therefore, this effect is more commonly included in chemical models of dilute dissociating gases as an excess term [49,50]. At T16000 K, water is fully dissociated and λ nuc is dominated by the kinetic part t J q,kin ¢ ( ), so that it approaches λ kin again at low densities.
For practical purposes, the generated data for the nuclear thermal conductivity are fitted to an analytic expression, which we recommend to use for T1000 K:  [42,43]. Temperatures are indicated via the colour code in the legend. The blue/ cyan stripe is the curve from Abramsonet al [11] scaled to T=600 K. 4 Although a T nuc l~scaling law exists for gases of hard spheres (derived within the Chapman-Enskog kinetic theory) [19], this is mere coincidence because underlying assumptions of that theory do not hold for liquid water.
, and all other coefficients are given in table 1. Figure 4 shows the representation of this fit together with the data points. The fit function is purely empirical and designed to become solely temperature-dependent in the low-density limit.
The nuclear thermal conductivity from equation (7) is intended to be added to the electronic contribution λ el from [16]. Figure 4 also includes λ el for comparison. The nuclear contribution is naturally dominant at low temperatures, and it can remain comparable to λ el even up to T≈20 000 K. This is very different from the electrical conductivity, which becomes dominated by electronic transport already near 5000 K [30,32]. Temperatures in the interiors of Uranus and Neptune are typically less than 6000 K, so that the thermal conductivity of water inside both planets is caused mainly by the interactions between the nuclei. This should also apply to possibly occurring superionic phases [51], which have similar proton diffusion coefficients as the fluid but strongly reduced electronic conductivity [30,32,52].

Summary
We have presented an efficient approach for calculating the nuclear thermal conductivity with DFT-MD simulations via the determination of effective radial force fields directly from the simulations. One can expect the method to perform generally well for anharmonic crystalline, glassy and superionic solids, dense atomic gases, liquids, and plasmas. Molecular and chemically reacting systems like water can be treated as well if their structure is not too complex. It remains to be explored up to which accuracy liquids consisting of much larger polyatomic molecules can also be described. For such cases, the DFT-MD+GKFF approach could be further expanded and improved straight-forwardly, for example, by using many-body potentials [53] or tight-binding coefficients [54] to describe non-radial interactions between the particles.  Compared to more rigorous techniques for thermal conductivity calculations with DFT-MD that employ different variants of polarisation theory [17,18], the DFT-MD+GKFF method has three main advantages. First, it does not enhance the computing time of DFT-MD runs at all because λ nuc is calculated within a fast postprocessing procedure. Second, the partial enthalpy per particle can be calculated directly from the force fields, either approximately via equation (5) or with exact expressions similar to those in [22]. Third, expressing the heat current as time derivative of polarisation of the electronic system is formally possible only when the material is an electronic insulator. Since the DFT-MD+GKFF circumvents the use of polarisation theory, it is certainly first choice for applications to semimetals and materials with thermally excited electrons.
Finally, we have made novel predictions for the nuclear thermal conductivity of fluid water across a large density-temperature region. Deeper insight into its characteristic changes with density and temperature was achieved by discussing the effective interactions between the particles and their influence on the microscopic heat current. Some of our qualitative predictions, e.g. the existence of a large, mostly temperature-independent region in the dense molecular (or weakly dissociated) fluid will require experimental verification. Apart from that, a fit formula is given for convenient implementation of our results in applications. s = å · , which allows for the derivation of linear phenomenological transport equations for the currents j l generated by forces X l that weakly perturb a thermodynamic equilibrium state [20]: where L ik are Onsager coefficients, which are scalar quantities in an isotropic system. Several different choices for j l and X l are possible, and they lead to different phenomenological equations that, in principle, describe the same physics but are more or less convenient to handle in certain cases [20, 55].

A.1. Thermal conductivity in variant A
For a two-component system in centre-of-mass frame, where j j 1 2 = -, equation (A.1) immediately leads to the transport equations: The thermal conductivity λ is usually measured in a stationary state, where j 1 =j 2 =0, that is reached after applying a temperature gradient and after the development of a concentration gradient [20]. Fourier's law for heat conduction is then derived by eliminating the gradient term containing the chemical potentials: With the above choice of currents and forces, we can derive the following microscopic expression for the thermal conductivity using linear response theory [19][20][21]: This expression can be directly evaluated in MD simulations, but its structure requires the calculation of three correlation functions. Both terms in the large bracket are of equal importance and they are subtracted from each other, which is a rather unfavourable structure for error propagation. The microscopic mass current density of species 1 reads: where N a is the particle number of species a and m i is the mass of particle i. For a system with pairwise interacting particles, the heat current density is defined as: where the modified heat current density reads: Here the partial enthalpies appear as additional quantities. The phenomenological equations of our two-component system then read (using again j 1 =−j 2 ): In these equations, the driving forces are now temperature gradient and gradients in chemical potential at constant T and p. The thermal conductivity λ in the stationary state acquires the same form as in variant A and reads: A major advantage in variant B is that L qq ¢ now contains the main contribution to λ. The small reduction from the second term is due to thermodiffusion and amounts only to very few percent at maximum [20], and we can therefore neglect it. The microscopic equation for L B qq l » ¢ from linear response theory then requires only the calculation of a single correlation function: A difficulty in this variant is the determination of the partial enthalpies, whose direct calculation with MD simulations is relatively involved [22,56], especially in simulations at constant volume. However, it is possible to approximate the thermodynamic derivative with the mean enthalpy per particle of the same species: This approximation neglects nonideal mixing contributions to the partial enthalpies and has performed very well in thermal conductivity calculations for Lennard-Jones mixtures [23,57].

A.3. Comparative validation of both variants against each other
In brief, we made two small and reasonable approximations in variant B, whereas variant A is a formally exact procedure. Moreover, note that systematic differences between both variants may come into effect only in a multi-component fluid system with inter-species diffusion, i.e. neither if water is fully molecular (where j j j , which can only be fulfilled if both currents are zero) nor in a solid state (where j 1 =0 and j 2 =0). Now we calculate both λ A and λ B for three representative isotherms and compare the results in figure A1. Distinguishment between kinetic (superscript: kin) and total (kinetic plus interaction, superscript: nuc) thermal conductivities is also made, just as described in the main text.
The chosen isotherms are 1000 K (fully molecular water up to to least 1.25 g cm −3 , see [58]), 6000 K (chemically reacting water at low densities, dense plasma at high densities), and 20000 K (fully dissociated water at low densities, dense plasma at high densities). The most noticeable difference between λ A and λ B is the much larger scattering in the values of λ A compared to that of λ B . The reason is that equation (A.6) contains a difference of two equally important terms in the bracket, each of which contributes significantly to the statistical uncertainty. This effect seems to be especially strong when the interaction contribution to the heat current is small (or set to zero as in A B , kin l ). There are no significant differences between the corresponding data sets of λ A and λ B .
In conclusion, the well-controllable approximations made in variant B do not introduce significant systematic errors into the resulting values for λ B . In contrast, although the λ A values are derived from formally exact equations, they are spoiled by much larger statistical uncertainties. Since the high computational cost of DFT-MD simulations limits the possibility of reducing the statistical error to arbitrary smallness, variant B is clearly the preferred path to follow in thermal conductivity calculations with the DFT-MD+GKFF method. This will hold even more strongly for systems that contain more than two nuclear species because λ A would then be given by a more complex arrangement of correlation functions, whereas the expression (A.15) can remain the same, regardless of the number of species. All of the main results in this work are derived within variant B, i.e. λ nuc =λ B .
Of course, it would be generally desirable to calculate the partial enthalpies (A.11) exactly, e.g. similar as in [22,56]. This may be required for systems where nonideal mixing effects in the thermodynamic functions are not as small as those in warm dense water and other planetary ices [24].