# Understanding the Transport Properties of Methane–Hydrogen Mixtures under the Interior Conditions of Ice Giants

Argha Jyoti Roy, Mandy Bethkenhagen, Armin Bergermann, Ronald Redmer

► To cite this version:

Argha Jyoti Roy, Mandy Bethkenhagen, Armin Bergermann, Ronald Redmer. Understanding the Transport Properties of Methane–Hydrogen Mixtures under the Interior Conditions of Ice Giants. The Astrophysical Journal, 2025, 990 (2), pp.212. 10.3847/1538-4357/adf5be. hal-05380568

HAL Id: hal-05380568

https://hal.science/hal-05380568v1

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Understanding the Transport Properties of Methane–Hydrogen Mixtures under the Interior Conditions of Ice Giants

Argha Jyoti Roy<fnref n="1" />,<fnref n="2" /> ID, Mandy Bethkenhagen<fnref n="3" /> ID, Armin Bergermann<fnref n="1" />,<fnref n="4" /> ID, and Ronald Redmer<fnref n="1" /> ID

Received 2025 April 22; revised 2025 July 22; accepted 2025 July 28; published 2025 September 10

## Abstract

Inaccuracy in our knowledge of the transport properties of relevant mixtures under planetary interior conditions is a roadblock in predicting the observable properties of planets in our solar system and beyond. In this work, we investigate methane–hydrogen mixtures using data sets obtained from density functional theory calculations for the electronic structure, combined with molecular dynamics simulations for the ions for a wide range of pressure and temperature. Hydrogen concentration significantly affects the equation of state (EOS) but has little influence on transport properties. We provide an analytical expression to model thermal EOS and transport properties as a function of hydrogen content, with the maximum deviation observed at low P–T conditions. These insights are particularly relevant to improve the planetary models and enhance our ability to predict the properties of “ice” giants and beyond.

Unified Astronomy Thesaurus concepts: Uranus (1751); Neptune (1096); Methane (1042)

## 1. Introduction

Ice giants Uranus and Neptune are the least explored and are still relatively unknown planets in our solar system. Besides the Voyager 2 flyby, both planets still await dedicated orbiter missions similar to Juno (Jupiter) and Cassini-Huygens (Saturn) (see M. Dougherty et al. 2009; S. J. Bolton et al. 2017) respectively. Ideally, we need missions that include entry probes to gather essential data for a better understanding of their interior structure, chemical composition, thermal evolution, and magnetic field structure. The recent Planetary Science and Astrobiology Decadal Survey (National Academies of Sciences, Engineering, and Medicine 2022) has recommended such an orbiter and probe to in situ explore Uranus as NASA’s highest-priority mission (see K. E. Mandt 2023). In order to study the planetary properties, interior models are used that solve the structure equations (D. J. Stevenson et al. 1982; T. Guillot 1999; T. Guillot et al. 2015). Standard three-layer models are used as long as accurate gravity data from an orbiter mission are unavailable. Many such models have been proposed for the concerned planets; e.g., adiabatic cooling schemes using the equation of state (EOS) of “ice” materials at 0 K fail to explain the rate of cooling concerning the present-day luminosity (W. Hubbard 1978; W. B. Hubbard & J. J. MacFarlane 1980). Thermal evolution calculations (J. J. Fortney et al. 2011; N. Nettelmann et al. 2013; E. F. Linder et al. 2019; L. Scheibe et al. 2019), assuming adiabatic cooling of three-layer models, are used to explain the present-day luminosity of Uranus and Neptune.

With the different approaches to the problem, the basic input is EOS of the relevant materials, i.e., H–He (outer envelope), an icy mixture H-C-N-O (inner envelope), and rocks (core). For instance, R. Redmer et al. (2011), found excellent

agreement of their interior model based on ab initio EOS data of hydrogen, helium, and water mixture with the proposed thin shell geometry of dynamo models by S. Stanley & J. Bloxham (2004) and S. Stanley & J. Bloxham (2006).

The ice giants are predominantly composed of so-called “ices” and “rocks,” unlike Jupiter and Saturn, which are dominated by hydrogen and helium. The term “ices” subsumes compounds such as water (H₂O), ammonia (NH₃), and methane (CH₄), while “rocks” stand for silicate materials (MgSiO₃, SiO₂). Deep within the planetary interior, the chemistry and physical nature of these “ices” play an important role in their properties. For instance, water and ammonia transform into superionic phases (M. French et al. 2009; M. Bethkenhagen et al. 2013), and methane decomposes into diamond and sinks deeper into the planet due to gravity, leaving possibly a metallic hydrogen-rich phase (M. Ross 1981; L. R. Benedetti et al. 1999; B. Cheng et al. 2023). This could also be a source of heating in these planets, especially in Neptune, due to its higher luminosity than Uranus (M. Podolak & R. Helled 2012; L. Scheibe et al. 2019). Even in the atmospheres of Uranus and Neptune, due to the abundance of hydrocarbons, these two planets appear in a distinct bluish hue due to absorption of the red light (P. G. Irwin et al. 2022). While studying pure systems offers useful insights, it is essential to recognize that hydrogen, methane, and other hydrocarbons are typically found in mixtures in nature, and the interaction among the atoms can substantially affect the understanding of planet’s observable properties.

Recent experimental campaigns have provided critical insights into the behaviors of pure fluids under these extreme conditions, which are particularly relevant to understanding the interiors of ice giants. For example, G. Tabak et al. (2024) reached about 400 GPa on precompressed pure methane samples in diamond-anvil cells (DACs), where the experiments showed consistency with density functional theory-molecular dynamics (DFT–MD) simulations predicting a polymeric region. Additionally, laser-driven shock compression on complex hydrocarbons such as polystyrene by

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Institut für Physik, Universität Rostock, D-18051 Rostock, Germany; argha.roy@uni-rostock.de</footnote>

<footnote n="2">² Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstr. 400, 01328 Dresden, Germany</footnote>

<footnote n="3">³ LULI, CNRS, CEA, Sorbonne Université, École Polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau, France</footnote>

<footnote n="4">⁴ SLAC National Accelerator Laboratory, 2575 Sand Hill Road, MS 72 Menlo Park, CA 94025, USA</footnote>
D. Kraus et al. (2017) and Z. He et al. (2022) showed the formation of diamonds under conditions similar to Uranus and Neptune. Findings from M. Frost et al. (2024) showed diamond formation from statically compressed polystyrene above a temperature of 2500 K, at pressures ranging from 19 to 27 GPa.

Several stable high-pressure compounds of methane-hydrogen mixture were studied in combination with X-ray diffraction experiments performed in DACs and ab initio simulation, as shown by U. Ranieri et al. (2022). Various mixtures under the same conditions have also been studied using DFT-MD techniques; however, the additive volume rule, applied at constant pressure and temperature, also known as linear mixing approximation (LMA), is the consistent one. The method is used for many different mixtures such as H-He (F. Soubiran & B. Militzer 2015), for gas giants with the inclusion of heavier species such as C, N, O, Si, and Fe (F. Soubiran & B. Militzer 2016), and for hydrocarbon C₂H₃ as a C and H mixture (J.-F. Danel & L. Kazandjian 2015). Similarly, M. Bethkenhagen et al. (2017) showed a maximum deviation of 4% in a mixture of CH₄, H₂O, and NH₃. Our previous work (A. J. Roy et al. 2024) testing the LMA also showed the same deviation of 4% in a binary mixture of CH₄ and H₂, as we found the clustering of carbon to occur at lower pressure in the presence of hydrogen than in pure methane and band gap closure is driven by metallic hydrogen at a lower temperature in comparison to pure methane (B. L. Sherman et al. 2012) and methane-water mixture (R. Hu et al. 2023).

Hence, we focus on extending the database of mixture and focus on transport properties. We use ab initio simulations and the Kubo-Greenwood formalism (R. Kubo 1957; D. Greenwood 1958) to understand better the EOS and transport properties of CH₄-H₂. Based on these results, we present an analytical model for the EOS and the transport properties as a function of hydrogen concentration. Our paper is organized as follows: In section 2, we explain the computational methods. Section 3 outlines the findings of this work. The EOS of the methane-hydrogen mixtures is discussed in Section 3.1, and we test the LMA with respect to the hydrogen molecular mass fraction. The transport properties of the pure substances and of the mixtures, as well as analytical model, are presented in Section 3.2. We discuss the behavior of the pure components and mixtures in the degenerate domain, i.e., the Wiedemann-Franz law in Section 3.3. In Section 4, we discuss the implications of our results, and Section 5 summarizes our findings.

## 2. Computational Details and Methods

The entire data set was obtained from DFT-MD simulations using the Vienna Ab Initio Simulation Package (VASP; G. Kresse & J. Hafner 1993, 1994; G. Kresse 1996; G. Kresse & J. Furthmüller 1996; J. Hafner 2008). The electronic and ionic dynamics are decoupled according to the Born-Oppenheimer approximation. The choice of the exchange-correlation (XC) functional plays a crucial role in determining the accuracy and computational costs of DFT calculations. As shown by A. Bergermann et al. (2024), computationally expensive XC functionals such as that of Heyd, Scuseria, and Ernzerhof (J. Heyd et al. 2003; see also J. Heyd & G. E. Scuseria 2004a, 2004b) shift the critical point of the first-order liquid-liquid phase transition in pure hydrogen by 22% in pressure and 150 K in temperature. However, the costs of such DFT-MD are

<table><caption>**Table 1**
Molecular Mass Fraction according to Equation (1) Related to the Molecule Number Combination</caption><thead><tr><th>System</th><th>X<sub>H2</sub></th><th>X<sub>CH4</sub></th></tr></thead><tbody><tr><td>Pure CH<sub>4</sub></td><td>0.000</td><td>1.000</td></tr><tr><td>[2:1] CH<sub>4</sub>-H<sub>2</sub></td><td>0.059</td><td>0.941</td></tr><tr><td>[2:2] CH<sub>4</sub>-H<sub>2</sub></td><td>0.116</td><td>0.884</td></tr><tr><td>[1:2] CH<sub>4</sub>-H<sub>2</sub></td><td>0.2008</td><td>0.7992</td></tr><tr><td>Pure H<sub>2</sub></td><td>1.000</td><td>0.000</td></tr></tbody></table>

almost 25 times and optical routines about 100 times (depending on the k-mesh grid) more than conventional generalized gradient approximation (GGA) calculations, and the cost increases further for complex molecules such as methane. Hence, for the evaluation of electronic properties, we used the GGA of Perdew, Burke, and Ernzerhof (PBE; J. P. Perdew et al. 1996) as the XC functional. The ionic trajectories are calculated using classical equations of motion and the Verlet algorithm, with the forces between electrons and ions obtained from the Hellman-Feynman theorem. This procedure combines DFT and classical MD simulations and allows us to treat systems at finite temperatures (P. Hohenberg & W. Kohn 1964; N. D. Mermin 1965; M. Weinert & J. Davenport 1992; R. M. Wentzcovitch et al. 1992). We set the plane wave cutoff energy to 1100 eV for every simulation and use projected augmented wave (PAW) pseudopotentials for hydrogen (PAW PBE H_h 06Feb2004) and carbon (PAW PBE C_h 06Feb2004; P. E. Blöchl 1994; G. Kresse & D. Joubert 1999). We control the temperature with the Nosé-Hoover thermostat (W. G. Hoover 1985; N. Shuichi 1991) with a Nosé-mass parameter corresponding to a temperature oscillation period of one per 40 time steps. The simulation box contains either two number combinations, i.e., 54 and 27 number of molecules of methane and hydrogen (2:1) and vice versa (1:2). We sampled the Brillouin zone with Baldereschi k-point (A. Baldereschi 1973). We choose a time step of 0.25 fs for every simulation. We continued the simulations for 8-12 ps after equilibration to evaluate thermodynamic and electrical properties. For pure H₂ and the 2:2 mixture, i.e., 54 molecules of methane and hydrogen, the data were obtained from the previous publication (A. J. Roy et al. 2024). For pure methane, we used the work of M. Bethkenhagen et al. (2017). All number permutations utilized in this work correspond to the molecular mass fractions (see Table 1) according to the following equation:

$$ \chi_i = \frac{n_i M_i}{\sum_i M_i}, \qquad (1) $$

where $n_i$ and $M_i$ are the number and molar mass of element i in the mixture, respectively.

The framework of linear response theory provides a means of calculating transport properties for our DFT-MD simulations. Though DFT, in principle, is a single-particle mean-field theory, the Kubo-Greenwood formalism works very well in tandem. In order to calculate the electronic properties for our DFT-MD simulations, we evaluate 40-100 snapshots for each simulation, depending on the system, using a Monkhorst-Pack k-point mesh of 4 × 4 × 4 to sample the Brillouin zone and perform static DFT runs to evaluate dynamic Onsager coefficients (see Equation (2)) in the framework of linear
<figure><img src="image_4.png" /><figcaption>Figure 1. EOS along an isotherm of T = 4000 K as a function of density for all the considered systems (color coded).</figcaption></figure>

response theory (B. Holst et al. 2011):

$$L_n(\omega) = \frac{2\pi(-1)^n}{V\omega} \sum_{k\alpha\beta} \mathbf{Re}(\langle k\alpha|\mathbf{v}|k\beta \rangle) \\ \otimes \langle k\beta|\mathbf{v}|k\alpha \rangle)(f_{k\alpha} - f_{k\beta}) \\ \times \left(\frac{\epsilon_{k\alpha} + \epsilon_{k\beta}}{2} - h_e\right)^n \\ \delta(\epsilon_{k\alpha} - \epsilon_{k\beta} - \hbar\omega), \qquad (2)$$

where $\omega$ is the frequency, $\hbar$ is the reduced Planck constant, $V$ is the volume of the simulation box, $|k\alpha\rangle$ are the Bloch states— where, $k$ is the reciprocal wavevector and $\alpha, \beta$ are the band indices, having eigen energy $\epsilon_{k\alpha}$, $\mathbf{v}$ is the velocity operator, $f_{k\mu}$ is the Fermi distribution function, and $h_e$ is the enthalpy per electron. The matrix elements $\langle k\alpha|\mathbf{v}|k\beta \rangle$ are calculated using optical routines of VASP. The delta function is replaced by a Gaussian broadening, which is required to circumvent the discrete spectrum of eigenvalues arising from the periodic boundary constraining the simulation box (B. Holst et al. 2008, 2011; D. Knyazev & P. Levashov 2013; M. French et al. 2022).

The macroscopic transport coefficients are determined by the zero-frequency value of the Onsager coefficients given by Equation (2). The electrical conductivity follows from

$$\sigma = e^2 \lim_{\omega \to 0} L_0(\omega), \qquad (3)$$

and the thermal conductivity is determined via

$$\lambda = \frac{1}{T} \lim_{\omega \to 0} (L_2(\omega) - L_1(\omega)L_0^{-1}(\omega)L_1(\omega)). \qquad (4)$$

### 3. Results

#### 3.1. EOS

In Figure 1, we present a comprehensive analysis of the EOS of pure and mixed methane–hydrogen systems at

$T=4000$ K. For the isotherms of $T=3000$ K and $6000$ K, we provide the EOS in Figure 6 in the Appendix. Quantities such as pressure $P$ and temperature $T$ are calculated by averaging over the simulation after equilibration. Temperature is determined through the motion of ions, while pressure is assessed according to the relation

$$P_T = \rho k_b T + P_e, \qquad (5)$$

where $P_e$ is the virial pressure of electrons, $\rho$ the ion number density, and $k_b$ is the Boltzmann constant. Our findings reveal a significant difference in the slope $\frac{\Delta P}{\Delta \rho}$, with hydrogen exhibiting a steeper slope. As the volume of the simulation box decreases, the dissociation of the molecules starts to increase, leading to a sharp rise, particularly for hydrogen. Dissociation creates a more homogeneous electronic charge density, intensifying the forces between electron clouds. Importantly, we have chosen our conditions to avoid any discontinuous first-order phase transition in pure hydrogen (W. Lorenzen et al. 2010; M. Zaghoo et al. 2016), although we cross the 50% dissociated line at $T=3000$ K, as noted by I. Tamblyn & S. A. Bonev (2010). The structural transformation of the pure components significantly influences electronic arrangements in the mixture. However, it is unclear whether the hydrogen molecules in the mixture or the mixture of C–H altogether undergo a first-order phase transition. Notably, D. Murayama et al. (2023) demonstrated that the strong covalent nature of C–C remains intact at densities up to $3 \text{ g cm}^{-3}$ at low temperatures, while methane starts to become unstable at around 4000 K.

In the case of pure methane, the rise of pressure with increasing density is at a slower rate compared to hydrogen, as methane remains molecular. However, as density increases further, methane dissociates, releasing hydrogen. When we increase the number of hydrogen molecules in the mixture to a molecular mass fraction of $\chi_{H_2} = 0.2008$, we see a sharp rise
<table><caption>**Table 2**
Fitted Parameters Obtained from Equation (7) for Different Hydrogen Mass Fractions at Selected Temperatures</caption><thead><tr><th rowspan="2">X<sub>H2</sub> / T [K]</th><th colspan="3">A [GPa]</th><th colspan="3">B [GPa]<sup>1/c</sup></th><th colspan="3">c</th></tr><tr><th>3000</th><th>4000</th><th>6000</th><th>3000</th><th>4000</th><th>6000</th><th>3000</th><th>4000</th><th>6000</th></tr></thead><tbody><tr><td>0.059</td><td>556.086</td><td>67.965</td><td>60.5415</td><td>-555.514</td><td>67.102</td><td>-31.589</td><td>0.6</td><td>2.36</td><td>2.16</td></tr><tr><td>0.116</td><td>157.324</td><td>63.5897</td><td>80.0283</td><td>-179.596</td><td>-41.995</td><td>-34.809</td><td>1.3</td><td>2.6</td><td>2.16</td></tr><tr><td>0.2008</td><td>282.103</td><td>153.063</td><td>243.737</td><td>152.890</td><td>-76.5746</td><td>-128.436</td><td>2.9</td><td>2.0</td><td>1.00</td></tr></tbody></table>

in the pressure. In the mixture with a molecular mass fraction of $X_{H_2}$ = 0.059, where the methane molecule is twice the hydrogen molecule, we see a decremental effect on pressure. Notably, the pressure for $X_{H_2}$ = 0.2008 mixture is 100 GPa higher than the pressure of the $X_{H_2}$ = 0.059 mixture at the same density, and this difference grows with the increasing density. For the mixture with $X_{H_2}$ = 0.116, the isotherm sits almost in the middle of the other isotherms of either mixture. The sharp changes in the slopes as a function of hydrogen concentration are attributed to additional hydrogen and the hydrogen from the methane. When either component dominates in molecular mass fraction, the overall motion of ions slightly dominates over the other.

The robust bonding between the carbon atoms results in larger molecules that occupy volume and make the pure methane hard to compress as a whole. This effect is observed in the all the mixtures, with formation of cluster of carbon at a lower pressure due to added hydrogen molecules (A. J. Roy et al. 2024); this may influence carbon cluster structures formed within these temperature and pressure conditions inside a planet such that they will be amorphous before reaching the deep interior.

Modeling the interior of planets requires not only the knowledge of the EOS of relevant materials but also an understanding of how it changes due to varying composition. Therefore, we provide a means of modeling the thermal EOS of the mixture of methane and hydrogen as a function of hydrogen molecular mass fraction (X), incorporating nonideal effects. Our model consists of the ideal contributions, i.e., the linear mixing of pressure and ideal configuration entropy, each weighted by the molecular mass fraction. Here, $P_i$ represents the pressure of the pure components methane and hydrogen, respectively, and R is the universal gas constant. To encapsulate the nonideal effect, an additional term $\Delta P$ is introduced:

$$P_{\text{fit}}(\mathcal{X}, \rho, T) = \sum_{i} \mathcal{X}_{i}P_{i}(\rho, T) + RT\rho\sum_{i} \mathcal{X}_{i}\log\mathcal{X}_{i} + \Delta P, \quad (6)$$

$$\Delta P(\mathcal{X}, \rho) = \mathcal{X}_{H_2}(1 - \mathcal{X}_{H_2})(A + B\rho^c). \quad (7)$$

The excess pressure $\Delta P$ is modeled using the polynomial expansion technique given by O. Redlich & A. T. Kister (1948) to map the excess thermodynamic quantity in binary mixtures. However, we modify the expression by replacing the summation with a polynomial fit to incorporate the increasing difference between the pressure from DFT-MD and ideal mixing, as given in Equation (7). We interpolate the pressure of pure components on to the density axis of the respective mixture, and we utilize least-square fitting to determine the best-fit values for the coefficients A, B, and c. The fitted parameters are given in Table 2 for each isotherm. In Figure 2, we present deviation of the thermal EOS from the DFT-MD for all the mixtures using Equations (6) and (7). For the lowest

and highest concentration of hydrogen, the deviation does not exceed more than 5%; however, for $X_{H_2}$ = 0.116, the deviation exceeds more than 8% at $1.00 \text{ g cm}^{-3}$, coinciding with dissociation of the hydrogen molecule. We note that the deviation tends to be high in regions of interpolation between isochores and/or isotherms, indicating that the fluctuation of cubic splines significantly contributes to the error.

Alongside the model for thermal EOS, planetary interior modelers rely on mixing the EOS of pure materials linearly at constant pressure and temperature, i.e., LMA. Based on the work of W. C. Demarcus (1958) and P. J. E. Peebles (1964), we extract density and specific internal energy as given in the following equations:

$$\frac{1}{\rho_{\text{LM}}(P, T)} = \sum_{j=1}^{N} \frac{\mathcal{X}_j}{\rho_j(P, T)}, \quad (8)$$

and, for the specific internal energy,

$$u_{\text{LM}}(P, T) = \sum_{j=1}^{N} \mathcal{X}_j u_j(P, T), \quad (9)$$

where, $\rho_j$, $u_j$, and $\mathcal{X}_j$ are the density, specific internal energy, and mass fraction, respectively, of individual components j of the mixture comprised of N components. The deviations from real mixtures are calculated using the following equations:

$$\Delta\rho(P, T) = \frac{\rho_{\text{LM}}(P, T)}{\rho_{\text{real}}(P, T)} - 1, \quad (10)$$

$$\Delta u(P, T) = u_{\text{LM}}(P, T) - u_{\text{real}}(P, T). \quad (11)$$

In Figure 7 (see the Appendix), we present the deviation in density and specific internal energy between the real mixture (DFT-MD) and LMA for molecular mass fraction of $X_{H_2}$ = 0.059 and 0.2008. The maximum deviation in density and specific energy occurs at 3000 K, with the deviations decreasing as the pressure and temperature increase. A similar trend in the deviations is observed as the temperature increases for $X_{H_2}$ = 0.2008. However, no such trend can be established with the pressure increase, with the maximum deviation of $3.8 \text{ kJ g}^{-1}$ occurring at 3000 K. The overall deviations in both quantities do not exceed 4%, holding the LMA good within the considered P-T space for the respective molecular mass fraction, $X_{H_2}$.

## 3.2. Transport Properties

Figure 3 and Figure 8 (see the Appendix) depict the electrical properties of the pure and mixed systems along the isotherms of T = 4000 K. We see a systematic increase of the electrical conductivity with increasing pressure. The system changes from a molecular to an atomic phase with a typical metallic conductivity. For pure methane, the electrical
<figure><img src="image_5.png" /><figcaption>Figure 2. Deviations in pressure from the real CH$_{4}$–H$_{2}$ mixtures $\chi_{H_2} = 0.059$ (left), $\chi_{H_2} = 0.116$ (middle), $\chi_{H_2} = 0.2008$ (right), and the predicted from analytical model Equation (6). Black dots denote the DFT–MD points. Between the isotherms and isochores, we interpolated the deviations cubically.</figcaption></figure>

conductivity is substantially lower than for pure hydrogen, even at higher pressures. One would expect a similar conductivity since every methane molecule yields four hydrogen atoms upon dissociation. However, carbon has a large scattering cross section and the formation of clusters, i.e., polymers, which restricts the transport of the free electrons. Doubling the number of hydrogen molecules increases the overall electrical conductivity (star symbol in Figure 3) by a factor of about 1.5 with respect to pure methane; however, it remains about 2.5 times less than that of pure H$_{2}$. A similar trend is observed for thermal conductivity, where heat transfer becomes less efficient with the presence of methane molecules that provide thermal insulation. Such an effect has been previously reported by W. Lorenzen et al. (2011) and M. French et al. (2012), where the electrical and thermal conductivity decreased systematically with increasing helium concentration. In the methane–hydrogen mixtures, the slower rise of the steepness of the electrical and thermal conductivity curve could be attributed to two different mechanisms, metallization of hydrogen occurs due to the closing of the band gap, with a high degree of sensitivity to pressure increase, which is not observed for carbon atoms. This is due to the electronic structure of carbon, which is more resilient to changes in pressure due to the substantial overlap of orbitals among carbon atoms.

While mixing rules are well established and tested to obtain EOS data of mixtures from the EOS from pure species, few models were developed to depict optical properties. It is incredibly challenging to obtain such mixing rules based on data for pure components. Approximate methods such as those given by C. Starrett et al. (2012, 2020), W. Johnson et al. (2006), and M. Y. Kuchiev & W. R. Johnson (2008) include only information of the pure species and lack the input on cross-species correlations. The recent work of A. J. White et al. (2024) developed analytical models for mixing optical properties of CH mixtures, where the mixing rule describing low-frequency excitations does not match the high-frequency regime. Most of these studies are done in conditions relevant to inertial fusion confinement experiments and not necessarily for planetary conditions, where bonding and intermolecular interactions are important. We propose the following fit formulas for the electrical and thermal conductivity along the isotherms, where x is the pressure axis:

$$ \lambda_{\text{fit}}, \sigma_{\text{fit}}^{(4000-6000)\text{K}} = A x^m e^{-B/x^n} + C x^d + E^k. \quad (12) $$

For the isotherm at T = 3000 K, we replace the power function in Equation (12) by a shifted sigmoid function as lower temperatures yield very low values for electrical conductivity, as shown in Equation (13). For thermal conductivity, Equation (12) provides a very good fit to the DFT data.

$$ \sigma_{\text{fit}}^{3000\text{K}} = A x^m e^{-B/x^n} + \frac{C}{1 + e^{-d(x-j)}} + E^k. \quad (13) $$

We use analytical fits to our data points to construct a model to describe the transport properties of mixtures, as shown in Figure 4 (e.g., we chose pure CH$_{4}$ and refer to Figure 9 in the Appendix for pure H$_{2}$) and the fitted parameters are provided in Table 3 and Table 4. Our fit reaches a minimal electrical and thermal conductivity value, with the first part of the fit equation being dominant, while the second part (power function or sigmoid function) dominates at higher pressure, with some deviations due to the analytical fit. The maximum deviation from the DFT points is observed to be about 9% and 12.5% in electrical and thermal conductivity, respectively. We observe that the maximum deviation for thermal conductivity occurs at around 75 GPa and 3000 K. For the electrical conductivity, it occurs at a lower pressure and temperature region, as well as at 25.15 GPa and 4000 K, as shown in Figures 4 and 9. Finding a trend of electronic properties along with temperature would require more isotherms, which we have not evaluated in this study. We provide $\sigma$ and $\lambda$ of pure species H$_{2}$ and CH$_{4}$ along isotherms, which were then used to model the transport properties of the mixtures.

Here we use two models, (i) a model proposed by A. Bar-ysevich & S. Cherkas (2011), Equation (14), and (ii) the reciprocal model according to Equation (15):

$$ \varphi_{\text{mix}} = \varphi_{\text{H}_2} \left( 1 + \frac{1 - \chi_{\text{H}_2}}{\frac{\chi_{\text{H}_2}}{n} - \frac{\varphi_{\text{H}_2}}{\varphi_{\text{H}_2} - \varphi_{\text{CH}_4}}} \right), \quad (14) $$

$$ \frac{1}{\varphi_{\text{eff}}} = \sum_{i} \frac{\chi_{i}}{\varphi_{i}}, \quad (15) $$

where $\varphi$ labels the electrical ($\sigma$) and thermal conductivity ($\lambda$). The numerical model (i) was introduced to model the electrical conductivity of copper obtained from wire
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Electrical conductivity mixing models as function of pressure along an isotherm of 4000 K. Black broken line: (i) Barysevich model according to Equation (14) with $n = -2.8$. Green broken line: (ii) reciprocal mixing Equation (15). DFT-MD data are color-coded symbols; circles denote pure H₂, and stars denote pure CH₄. Diamond, triangles, and squares denote the various mixed systems.</figcaption></figure>

explosions, including experimental and magnetohydrodynamic simulation data. The transport theory for heterogeneous systems allows modeling the conductivity also in the two-phase region, i.e., vapor and liquid (see Equation (14)). In our work, we used the same equation to model the electrical and thermal conductivity with the tuning parameter $n_σ = n_λ = -2.8$, for which our model predicts the trend better. Similarly, we also use the “resistance in series” model with the appropriate number fraction. Both models (i) and (ii)

show a good agreement with the electrical and thermal conductivity data obtained from Equations (3) and (4) for pure methane and hydrogen (see Figure 3, Appendix Figure 8). Despite the origin of Equation (14) requiring a two-phase region, its predictions agree well with the trend of the simulation data, whereas the “resistance in series” model deviates substantially from DFT data as the pressure increases, predicting a decrease in the conductivity stronger than Equation (14).
<figure><img src="image_7.png" /><figcaption>(a) Electrical conductivity</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Thermal conductivity</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Deviations in electrical conductivity</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Deviations in thermal conductivity</figcaption></figure>

**Figure 4.** (a), (b) Electrical and thermal conductivity of pure CH₄ as a function of pressure along the isotherms. The hollow circles represents the points obtained from the fit Equations (12) and (13), and hollow stars denote DFT simulation data using Equation (2). The deviations between the predicted values of the fits using Equations (12) and (13) and the DFT data for electrical and thermal conductivity are shown in (c) and (d), respectively. Between the isotherms and isochores, we interpolated the deviations cubically.

### 3.3. Lorenz Number

A well-known relation between the electrical and thermal conductivity is the Lorenz number $L$:

$$L = \frac{\lambda}{T\sigma} \left(\frac{e}{k_b}\right)^2. \qquad (16)$$

For simple metals, i.e., for the degenerate electron gas, this relation yields a constant $L_0 = \pi^2/3$ (in units of $(e/k_b)^2$) known as the Wiedemann-Franz law. We calculate the Lorenz number to test the validity of this relation for the mixtures and conditions, as considered here. The results are shown in Figure 5. Pure hydrogen at high pressures shows no temperature dependence, typical for metallic systems, satisfying the Wiedemann-Franz law. However, below 100 GPa, we see deviations from the value $L_0$ in the molecular region with nonmetallic conductivity. For pure methane (star symbol), we also see deviations from $L_0$ at low pressure, irrespective of temperature. $L_0$ is only valid at higher pressures. At lower

pressures, temperature has a greater effect on disrupting the bonds between carbon and hydrogen.

The temperature dependence is quite strong in the case of CH₄-H₂ mixtures. We see a significant deviation from the value of $L_0$ along the isotherm $T=3000$ K. The calculated Lorenz number for $\chi_{H_2} = 0.2008$ is lower than $L_0$ even at higher pressures. This effect could be tied to the stable clustering of carbon atoms, enhancing electron scattering with the mobile hydrogen atoms, which reduces the electronic thermal conductivity. With increasing temperature, this effect starts to disappear. At $T=6000$ K, we find $L_0$ valid for all mixtures. Similarly, we encounter the same deviation for the $\chi_{H_2} = 0.059$ and $\chi_{H_2} = 0.116$, although not as strong as for the 1:2 mixture.

## 4. Discussion

For pure methane, we find the direct current (dc) value of electrical conductivity to be approximately 2.63 orders of magnitude larger than that of water at 4000 K and $2 \text{ g cm}^{-3}$
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Lorenz number of pure components and methane–hydrogen mixtures as a function of pressure for different temperatures. Symbols ○, ∇, □, ◊, ⋆ denote pure hydrogen, mixtures of methane and hydrogen for different mass fraction $\chi_{H_2}$, and pure methane, respectively. Brown, yellow, and blue color symbols denote 3000, 4000, and 6000 K, respectively. The dashed gray line shows the limiting value of $L_0 = \pi^2/3$, known as the Wiedemann–Franz law.</figcaption></figure>

(M. French et al. 2010); however, the difference in orders of magnitude drops to 1.4 at 6000 K for the same density. Adding hydrogen to pure methane affects the dc values of electrical properties, as hydrogen atoms upon dissociation act as charge carriers. Using Mott’s criterion for the minimum metallic electrical conductivity (N. Mott 2004; W. Lorenzen et al. 2011) predicts the metallization in H–He mixtures to occur around 100 GPa at 4000 K for a hydrogen concentration of 66%. We find similar metallization pressures at the same temperature for a substantially lower hydrogen molecular mass

fraction of 5% in a CH$_4$–H$_2$ mixture and a further drop in metallization pressure of about 20 GPa when the hydrogen concentration is increased to 20%. We also see a rise in the electrical and thermal conductivity values with increasing hydrogen concentration in the CH$_4$–H$_2$ mixtures; however, the dc conductivity values are significantly less than what we observe for pure hydrogen, despite the dominating number of hydrogen atoms. It can have a substantial effect on the modeling of planetary interiors as lower values of thermal conductivity may reduce the convective nature of heat
transport, which some models rely on. The situation becomes
more complex when other high-pressure phases of carbon
(such as diamond) and superionic phases of water and
ammonia are taken into consideration. The nonmetallic nature
of some of the phases can enhance phonon scattering, further
reducing the thermal conductivity and the efficacy of
convective heat transport. These key differences may mediate
into variations in the age and radius of the planets, as changes
in conductivity gradients may or may not allow for facilitating
multipolar magnetic fields with shallower density stratification
(L. D. Duarte et al. 2013; K. Soderlund & S. Stanley 2020) or
allow more advanced interior models that incorporate the
chemical properties of planetary matter.

**5. Conclusion**

We treated mixtures of methane and hydrogen for different
molecular mass fractions of hydrogen and studied the effect on
the EOS and transport properties. We focused on a temperature
range of 3000-6000 K and hydrogen molecular mass fractions
of 0, 0.059, 0.116, 0.2008, and 1.000. We present an analytical
expression to model the thermal EOS obtained from DFT-MD
simulations. Our model incorporates nonideal effect as a
function of density and molecular mass fraction. We find our
model to be applicable within the chosen temperature range. In
Figure 2, we find a maximum deviation of 8.3% for the
mixture with a molecular mass fraction of $\chi_{H_2}$ = 0.116. For
the lowest and highest concentration of hydrogen in the
mixtures, the deviations are below 5%. At higher density and
temperature, the model produces pressure with a deviation of
less than 2%, despite the transient nature of carbon bonds.
Experiments reported the formation of diamonds above
2000 K, which we do not see in our simulation box and that
could well be attributed to finite size effect. However, we can
infer the diffusion coefficients (A. J. Roy et al. 2024) of carbon
atoms. We suspect there may be demixing between carbon and
hydrogen, driven by the metallic nature of hydrogen, and the
two regions are separated from each other due to the large
scattering cross section of carbon and polymeric bond
formation of carbon. Calculations for the transport coefficients
show a slight increment in electrical and thermal conductivity
with an increase in the hydrogen molecular mass fraction. At

higher pressure, the insulator–metal transition is driven by the
delocalization of the wave function of hydrogen (D. Muray-
ama et al. 2023) and the dissociation of methane. We use an
analytical fit for all systems’ electrical and thermal conductiv-
ity, and we find good agreement with the DFT–MD data, and
we observe a maximum deviation at low pressure and
temperature conditions for both the quantities. Our fit allows
us to extrapolate to minimal pressure values, where DFT
calculations become expensive.

Future DFT-MD simulations for relevant materials will
advance the modeling of ice giants’ interiors and dynamo
models and aid experimental campaigns. Additionally, our
findings will contribute to understanding planets like Uranus
and Neptune and their properties, which will, in turn, highlight
the ever-increasing number of Neptune-like and super-Earth
planets.

**Acknowledgments**

We thank Nadine Nettelmann, Martin Preising, and Uwe
Kleinschmidt for their helpful discussions. This work was
supported by the North German Supercomputing Alliance
(HRN), the ITMZ of the University of Rostock, and German
Israel Foundation for Scientific Research and Development
(GIF) within the Research grant No. I-1486-303.7/2019. We
also thank Helmholtz-Zentrum Dresden-Rossendorf e.V.
(HZDR) for support.

**Appendix**

In Figure 6 we plot the EOS of pure and the mixtures for
isotherms T= 3000 K and 6000 K. Figure 7 shows the
deviation between DFT-MD EOS and LMA for molecular
hydrogen concentration of 5.9% and 20%, respectively. In
Figure 8 we show the dc values of thermal conductivity, the
influence of hydrogen molecules and the models as given in
Equations (14) and (15) and finally, Figure 9 shows the
deviation of fitted equations according to Equations (12) and
(13) from the simulation data, with Tables 3 and 4 holding the
values of fitted parameters.
<figure><img src="image_12.png" /><figcaption>**Figure 6.** EOS isotherms of T = 3000 K and 6000 K as a function of density for all the concerned systems. Symbols ◦, ∇, □, ◊, ⋆ denote pure hydrogen, mixtures of methane and hydrogen for different mass fraction $\chi_{H_2}$, and pure methane, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Deviation of densities (light blue; in percent) and specific internal energy (kJ g⁻¹; teal) from the real mixture having a molecular mass fraction of $\chi_{H_2}$ = 0.059 (top) and 0.2008 (bottom) according to Equations (10) and (11), respectively.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Thermal conductivity mixing models as function of pressure along an isotherm of 4000 K: Black broken lines: (i) Barysevich model Equation (14) ($n = -2.8$). Green broken lines: (ii) reciprocal mixing Equation (15). DFT-MD data are color-coded symbols, circles denote pure H₂, and stars denote CH₄. Diamonds, triangles, and squares denote the various mixed systems.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>(a) Electrical conductivity</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Thermal conductivity</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Deviation in electrical conductivity</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Deviation in thermal conductivity</figcaption></figure>

**Figure 9.** (a), (b) Electrical and thermal conductivity of pure H₂ as a function of pressure along the isotherms. The hollow circles represent the points obtained from the fit Equations (12) and (13), and hollow stars are DFT simulation data using Equation (2). The deviations between the predicted values of the fits using Equations (12) and (13) vs. DFT data for electrical and thermal conductivity are shown in (c) and (d), respectively. Between the isotherms and isochores, we interpolated the deviations cubically.

<table><caption>**Table 3**
Parameters Obtained from Fitting the Analytical Equation for Electrical Conductivity Equations (12) and (13)</caption><thead><tr><th rowspan="2">Parameters [Units]/ T [K]</th><th colspan="3">H<sub>2</sub></th><th colspan="3">CH<sub>4</sub></th></tr><tr><th>3000</th><th>4000</th><th>6000</th><th>3000</th><th>4000</th><th>6000</th></tr></thead><tbody><tr><td>A [MS/(GPa)<sup>n</sup>]</td><td>0.071980</td><td>4.8654 &times; 10<sup>-8</sup></td><td>-0.019929</td><td>0.2336</td><td>0.033782</td><td>0.03241</td></tr><tr><td>m</td><td>0.60</td><td>2.00</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td></tr><tr><td>B [(GPa)<sup>n</sup>]</td><td>2025.2</td><td>-3478.6</td><td>6.062 &times; 10<sup>5</sup></td><td>2.6302 &times; 10<sup>4</sup></td><td>3.3782 &times; 10<sup>3</sup></td><td>2.170 &times; 10<sup>3</sup></td></tr><tr><td>n</td><td>1.85</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td></tr><tr><td>C</td><td>-0.167</td><td>0.102</td><td>0.0464</td><td>-0.128</td><td>0.0142</td><td>0.0307</td></tr><tr><td>d</td><td>0.47415</td><td>0.55615</td><td>0.6540</td><td>1.000</td><td>0.5482</td><td>0.4498</td></tr><tr><td>E [MS/k]</td><td>-0.13283</td><td>-0.90996</td><td>-0.6557</td><td>0.70779</td><td>-0.68997</td><td>-3.740 &times; 10<sup>-3</sup></td></tr><tr><td>k</td><td>5.00</td><td>5.00</td><td>3.00</td><td>8.00</td><td>7.00</td><td>7.00</td></tr><tr><td>j</td><td>50.00</td><td>...</td><td>...</td><td>80.00</td><td>...</td><td>...</td></tr></tbody></table>

Note. m, n, d, k are dimensionless quantities.
<table><caption>**Table 4**
Parameters Obtained from Fitting the Analytical Equation for Thermal Conductivity Equation (12)</caption><thead><tr><th rowspan="2">Parameters [Units]/ T [K]</th><th colspan="3">H<sub>2</sub></th><th colspan="3">CH<sub>4</sub></th></tr><tr><th>3000</th><th>4000</th><th>6000</th><th>3000</th><th>4000</th><th>6000</th></tr></thead><tbody><tr><td>A [w/(m·K·(GPa)<sup>n</sup>)]</td><td>6.733 × 10<sup>-6</sup></td><td>2.1067</td><td>-2.324 × 10<sup>-5</sup></td><td>1.4559</td><td>-0.010493</td><td>-0.01819</td></tr><tr><td>m</td><td>2.00</td><td>0.90</td><td>2.00</td><td>0.900</td><td>1.80</td><td>1.20</td></tr><tr><td>B [(GPa)<sup>n</sup>]</td><td>-3.4034 × 10<sup>3</sup></td><td>1.1612 × 10<sup>3</sup></td><td>8.773 × 10<sup>4</sup></td><td>1.3833 × 10<sup>4</sup></td><td>1.0211 × 10<sup>-4</sup></td><td>-2.051 × 10<sup>3</sup></td></tr><tr><td>n</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>-2.00</td><td>2.00</td></tr><tr><td>C [w/(m·K·(GPa)<sup>d</sup>)]</td><td>20.5</td><td>-0.0487</td><td>8.27</td><td>-1.45</td><td>11.2</td><td>4.72</td></tr><tr><td>d</td><td>0.43756</td><td>1.0427</td><td>0.6217</td><td>0.89033</td><td>0.41380</td><td>0.6128</td></tr><tr><td>E [w/(m·K)<sup>1/2</sup>]</td><td>-2.5557</td><td>-1.7889</td><td>-2.155</td><td>1.6748</td><td>-38.199</td><td>0.01312</td></tr><tr><td>k</td><td>5.00</td><td>4.00</td><td>5.00</td><td>8.00</td><td>1.00</td><td>8.00</td></tr></tbody></table>

Note. m, n, d, k are dimensionless quantities.

## ORCID iDs

Argha Jyoti Roy https://orcid.org/0009-0006-4798-650X

Mandy Bethkenhagen https://orcid.org/0000-0002-1838-2129

Armin Bergermann https://orcid.org/0000-0002-2849-1791

Ronald Redmer https://orcid.org/0000-0003-3440-863X

## References

* Baldereschi, A. 1973, PhRvB, 7, 5212

* Barysevich, A., & Cherkas, S. 2011, PhPl, 18, 052703

* Benedetti, L. R., Nguyen, J. H., Caldwell, W. A., et al. 1999, Sci, 286, 100

* Bergermann, A., Kleindienst, L., & Redmer, R. 2024, JChPh, 161, 234303

* Bethkenhagen, M., French, M., & Redmer, R. 2013, JChPh, 138, 234504

* Bethkenhagen, M., Meyer, E. R., Hamel, S., et al. 2017, ApJ, 848, 67

* Blöchl, P. E. 1994, PhRvB, 50, 17953

* Bolton, S. J., Lunine, J., Stevenson, D., et al. 2017, SSRv, 213, 5

* Cheng, B., Hamel, S., & Bethkenhagen, M. 2023, NatCo, 14, 1104

* Danel, J.-F., & Kazandjian, L. 2015, PhRvE, 91, 013103

* Demarcus, W. C. 1958, AJ, 63, 2

* Dougherty, M., Esposito, L., & Krimigis, S. M. 2009, Saturn from Cassini-Huygens (Berlin: Springer)

* Duarte, L. D., Gastine, T., & Wicht, J. 2013, PEPI, 222, 22

* Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011, ApJ, 729, 32

* French, M., Becker, A., Lorenzen, W., et al. 2012, ApJS, 202, 5

* French, M., Mattsson, T. R., Nettelmann, N., & Redmer, R. 2009, PhRvB, 79, 054107

* French, M., Mattsson, T. R., & Redmer, R. 2010, PhRvB, 82, 174108

* French, M., Röpke, G., Schörner, M., et al. 2022, PhRvE, 105, 065204

* Frost, M., McWilliams, R. S., Bykova, E., et al. 2024, NatAs, 8, 174

* Greenwood, D. 1958, PPS, 71, 585

* Guillot, T. 1999, Sci, 286, 72

* Guillot, T., Gautier, D., & Schubert, G. 2015, Treatise on Geophysics (2nd ed.; Amsterdam: Elsevier)

* Hafner, J. 2008, JCoCh, 29, 2044

* He, Z., Rödel, M., Lügtert, J., et al. 2022, SciA, 8, eabo0617

* Heyd, J., & Scuseria, G. E. 2004a, JChPh, 121, 1187

* Heyd, J., & Scuseria, G. E. 2004b, JChPh, 120, 7274

* Heyd, J., Scuseria, G. E., & Ernzerhof, M. 2003, JChPh, 118, 8207

* Hohenberg, P., & Kohn, W. 1964, PhRv, 136, B864

* Holst, B., French, M., & Redmer, R. 2011, PhRvB, 83, 235120

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

* Hoover, W. G. 1985, PhRvA, 31, 1695

* Hu, R., Li, Y., Wang, C., Hu, Y., et al. 2023, PhRvB, 107, 174108

* Hubbard, W. 1978, Icar, 35, 177

* Hubbard, W. B., & MacFarlane, J. J. 1980, JGRB, 85, 225

* Irwin, P. G., Teanby, N. A., Fletcher, L. N., et al. 2022, JGRE, 127, e2022JE007189

* Johnson, W., Guet, C., & Bertsch, G. 2006, JQSRT, 99, 327

* Knyazev, D., & Levashov, P. 2013, ComMS, 79, 817

* Kraus, D., Vorberger, J., Pak, A., et al. 2017, NatAs, 1, 606

* Kresse, G. 1996, PhRvB, 54, 11169

* Kresse, G., & Furthmüller, J. 1996, ComMS, 6, 15

* Kresse, G., & Hafner, J. 1993, PhRvB, 47, 558

* Kresse, G., & Hafner, J. 1994, PhRvB, 49, 14251

* Kresse, G., & Joubert, D. 1999, PhRvB, 59, 1758

* Kubo, R. 1957, JPSJ, 12, 570

* Kuchiev, M. Y., & Johnson, W. R. 2008, PhRvE, 78, 026401

* Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85

* Lorenzen, W., Holst, B., & Redmer, R. 2010, PhRvB, 82, 195107

* Lorenzen, W., Holst, B., & Redmer, R. 2011, PhRvB, 84, 235109

* Mandt, K. E. 2023, Sci, 379, 640

* Mermin, N. D. 1965, PhRv, 137, A1441

* Mott, N. 2004, Metal-Insulator Transitions (2nd ed.; Boca Raton, FL: CRC Press)

* Murayama, D., Ohmura, S., Kodama, R., & Ozaki, N. 2023, JAP, 134, 095902

* National Academies of Sciences, Engineering, and Medicine 2022, Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023–2032 (Washington, DC: The National Academies Press)

* Nettelmann, N., Helled, R., Fortney, J., & Redmer, R. 2013, P&SS, 77, 143

* Peebles, P. J. E. 1964, ApJ, 140, 328

* Perdew, J. P., Burke, K., & Ernzerhof, M. 1996, PhRvL, 77, 3865

* Podolak, M., & Helled, R. 2012, ApJL, 759, L32

* Ranieri, U., Conway, L. J., Donnelly, M.-E., et al. 2022, PhRvL, 128, 215702

* Redlich, O., & Kister, A. T. 1948, Ind. Eng. Chem. Res., 40, 345

* Redmer, R., Mattsson, T. R., Nettelmann, N., & French, M. 2011, Icar, 211, 798

* Ross, M. 1981, Natur, 292, 435

* Roy, A. J., Bergermann, A., Bethkenhagen, M., & Redmer, R. 2024, PCCP, 26, 14374

* Scheibe, L., Nettelmann, N., & Redmer, R. 2019, A&A, 632, A70

* Sherman, B. L., Wilson, H. F., Weeraratne, D., & Militzer, B. 2012, PhRvB, 86, 224113

* Shuichi, N. 1991, PThPS, 103, 1

* Soderlund, K., & Stanley, S. 2020, RSPTA, 378, 20190479

* Soubiran, F., & Militzer, B. 2015, ApJ, 806, 228

* Soubiran, F., & Militzer, B. 2016, ApJ, 829, 14

* Stanley, S., & Bloxham, J. 2004, Natur, 428, 151

* Stanley, S., & Bloxham, J. 2006, Icar, 184, 556

* Starrett, C., Clérouin, J., Recoules, V., et al. 2012, PhPl, 19, 102709

* Starrett, C., Shaffer, N., Saumon, D., et al. 2020, HEDP, 36, 100752

* Stevenson, D. J., et al. 1982, AREPS, 10, 257

* Tabak, G., Rygg, J., Millot, M., et al. 2024, PhRvB, 109, 064102

* Tamblyn, I., & Bonev, S. A. 2010, PhRvL, 104, 065702

* Weinert, M., & Davenport, J. 1992, PhRvB, 45, 13709

* Wentzcovitch, R. M., Martins, J. L., & Allen, P. B. 1992, PhRvB, 45, 11372

* White, A. J., Craven, G. T., Sharma, V., & Collins, L. A. 2024, PhPl, 31, 042706

* Zaghoo, M., Salamat, A., & Silvera, I. F. 2016, PhRvB, 93, 155128