Theoretical evidence of h-he demixing under jupiter and saturn conditions

ABSTRACT The immiscibility of hydrogen-helium mixture under the temperature and pressure conditions of planetary interiors is crucial for understanding the structures of gas giant planets

(e.g., Jupiter and Saturn). While the experimental probe at such extreme conditions is challenging, theoretical simulation is heavily relied in an effort to unravel the mixing behavior of

hydrogen and helium. Here we develop a method via a machine learning accelerated molecular dynamics simulation to quantify the physical separation of hydrogen and helium under the conditions

of planetary interiors. The immiscibility line achieved with the developed method yields substantially higher demixing temperatures at pressure above 1.5 Mbar than earlier theoretical data,

but matches better to the experimental estimate. Our results suggest a possibility that H-He demixing takes place in a large fraction of the interior radii of Jupiter and Saturn, i.e.,

27.5% in Jupiter and 48.3% in Saturn. This indication of an H-He immiscible layer hints at the formation of helium rain and offers a potential explanation for the decrease of helium in the

atmospheres of Jupiter and Saturn. SIMILAR CONTENT BEING VIEWED BY OTHERS EVIDENCE OF HYDROGEN−HELIUM IMMISCIBILITY AT JUPITER-INTERIOR CONDITIONS Article 26 May 2021 EVIDENCE FOR

SUPERCRITICAL BEHAVIOUR OF HIGH-PRESSURE LIQUID HYDROGEN Article 09 September 2020 UNDERSTANDING DENSE HYDROGEN AT PLANETARY CONDITIONS Article 01 September 2020 INTRODUCTION The

gravitational might of giant planets has played a key role in the formation of our solar system1. Jupiter and Saturn are the largest and most massive gas giants in the Sun’s planetary

system. Current models of Jupiter and Saturn suggest that the structures of both planets are similar in composition, both containing a visible cloud top, layers of gaseous hydrogen, liquid

hydrogen, and metallic hydrogen, and possibly a rocky core2,3. Helium is present in all three layers of hydrogen, albeit with different abundances. The ratio of helium mass density to the

sum of helium and hydrogen mass densities is 0.238 ± 0.054 in the atmosphere of Jupiter, and 0.06-0.085 in Saturn, both of which are conclusively lower than the estimated protosolar helium

mass fraction (0.275 ± 0.016). The observed helium reduction is thought to be caused by the demixing of hydrogen and helium that precipitate toward deeper layers in the planet’s interior.

The sinking helium, through the exchange of gravitational potential energy to thermal energy, is thought to be an additional energy source to power Saturn’s excess luminosity7. Thus, a

complete diagram of the solubility of helium in hydrogen at planetary _P_-_T_ conditions is highly required for an accurate modeling of Jupiter, Saturn, and other gas giants like them. A

recent experiment on giant planet modeling, through laser-driven shock compression of H-He mixtures, reveals a large region of H-He separation (~15% of the radial range) under Jovian

interior conditions8. However, previous first-principles simulations9,10,11,12,13 point toward much lower demixing temperatures and smaller immiscibility regions compared to the experiment.

In particular, the demixing temperatures predicted by different theoretical models fall on either side of the adiabatic lines of Jupiter and Saturn with a large discrepancy of the order 2000

 K. Such inconsistency would lead to completely different models for the planets’ internal structures. In addition to the range of separation, the separation intensity is also important for

planetary modeling – an ab initio calculation suggests excessive separation might cause Saturn’s cooling time to be longer than the lifespan of the Solar System14. In theory, the miscibility

range is usually determined based upon the Gibbs free energy of mixing (\(\Delta G\)), while different levels of theories (non-ideal entropy, choice of exchange-correlation functional) will

yield slightly different results. On the other hand, the physical separation process plays an important role in the demixing of hydrogen and helium, but it has not been accounted in the

previous studies. To this end, we will provide a new method to quantify the physical separation of H-He, and to determine the miscibility range directly from the separation process. We will

also present a nonequilibrium approach for an improved evaluation of \(\Delta G\) in large systems (~27,000 atoms), and use it to obtain the miscibility. Previous density functional theory

(DFT) evaluations of \(\Delta G\) are usually carried out in systems with a small number of atoms (several tens to hundreds) to ensure sufficient mixing at all helium fractions and _P_-_T_

conditions. This is because the conventional calculation method refers to equilibrium state, whereas a large system cannot represent the Gibbs free energy of a specific helium abundance once

it undergoes demixing during equilibration. The difficulty of reaching equilibrium in ab initio molecular dynamics (MD) simulation for large systems should also be noted. Aside from low

resolution, small systems tend to have non-negligible thermal fluctuations which may smear out the free energy difference between different configurations. One way to address this issue is

to construct simulations with a sufficiently large number of atoms to enable statistical sampling, while at the same time maintaining a first-principles quality. Effective efforts including

matter at extreme conditions were made by deep learning recently15,16,17. In this study, the MD simulation is substantially scaled up (to 27,000 atoms) using machine learning representation

of potential energy surfaces. We constructed a deep-learning potential (DLP)18 based on the strongly constrained and appropriately normed (SCAN) functional19 with van der Waals (vdW)

interaction from rVV1020. Details of the potential and its benchmarking are provided in Supplementary Information Note 2. RESULTS The separation behaviors of H-He mixture are captured in the

DLP accelerated MD simulations. However, a rigorous analysis of the MD trajectory must go beyond a visual observation to quantify the separation in the He-poor and He-rich zones. The

intuitive mean-field methods, i.e., calculating and averaging the thermodynamic properties in small subregions of the system, always yield unstable and non-convergent results. Since the

subregions are treated as individuals, correlations among them are neglected and fluctuations are smoothed out. This is more problematic when dealing with extremely skewed compositions under

real planetary conditions. Therefore, we designed an atomic miscibility analysis approach based on reweighted conditional probabilities (_P_cond) for ‘atoms in neighborhoods’ to quantify

the extent of separation in H-He mixtures. This method follows a similar statistical strategy to that of measuring the disorder in fluid mixture using conditional entropy21. Through

tabulating the probability at which particular types of atoms are presented in the neighborhood of each atom in the systems, we can circumvent the aforementioned limitations and obtain a

more precise numerical criterion for mixing and separation. We apply this method to the DFT-MD trajectories of 2048 atoms with helium abundances of _X_He = 0.073 or 0.2. In addition to

validating the accuracy of the DLP, it also demonstrates that the method can accurately differentiate helium abundance on a case-by-case basis, even when the separation levels are very close

(see Fig. S11 in Supplementary Information Note 2). The validation of this method and its comparison to the mean field method are presented in Supplementary Information Note 3. The main

concept of AMA-RCP method is to determine the _P_cond for each atom in the system based on the atomic composition of the _N__nb_ nearest neighbors. All atoms are categorized by their _P_cond

with a weight assigned for each category. With the categorized and reweighted results (Fig. 1a), \({x}_{1}\) and \({x}_{2}\) can be calculated directly,

$$\left\{\begin{array}{c}{x}_{1}=\frac{({N{{{\rm{He}}}}}_{{{{\rm{He}}}}-{{{\rm{poor}}}}})}{({N{{{\rm{He}}}}}_{{{{\rm{He}}}}-{{{\rm{poor}}}}})+({N{{{\rm{H}}}}}_{{{{\rm{He}}}}-{{{\rm{poor}}}}})}\\

{x}_{2}=\frac{({N{{{\rm{He}}}}}_{{{{\rm{He}}}}-{{{\rm{rich}}}}})}{({N{{{\rm{He}}}}}_{{{{\rm{He}}}}-{{{\rm{rich}}}}})+({N{{{\rm{H}}}}}_{{{{\rm{He}}}}-{{{\rm{rich}}}}})}\end{array}\right..$$

(1) The value ∆_x_ = _x_2 – _x_1 represents the abundance difference for helium in He-rich and He-poor zones, which is used as a measure of immiscibility (See Methods for the definition of

immiscibility using \(\Delta x\)). Figure 1b shows three representative examples with _X_He = 0.073 for illustrating immiscible, weakly immiscible and miscible systems determined using this

method. In the miscible system (12,000 K and 2 Mbar), the filtered hydrogen and helium atoms are few in number and gathered, indicating that phase separation does not occur, but rather the

corresponding regions exhibit significant fluctuations. In contrast, the immiscible system filtered with the same _P_cond shows the prevalence of phase separation in the simulation box. With

the AMA-RCP, we conduct a rigorous analysis of mixing and separation behaviors in the MD trajectories of large H-He mixtures calculated using LAMMPS22 package with DLP (see “Methods”

section for simulation details). Long-time trajectories obtained after the systems have reached equilibrium are used as inputs to obtain the degree of immiscibility for each system. Multiple

helium abundances are used, including _X_He = 0.073 and 0.089 calculated with 27,000-atom simulation boxes and _X_He = 0.2, 0.357, 0.5, 0.643, 0.8, and 0.91 with 9604-atom simulation boxes.

The MD simulations are carried out at various temperatures and pressures to construct the miscibility diagram, for which a total of 529 MD trajectories are generated. To compare with

previously published results9,13, we also calculated the miscibility lines at 4 and 10 Mbar using thermodynamic integration for calculating free energy \(\Delta G\) with 64-atom systems as

the equation, $$\Delta G\left({X}_{{{\rm{He}}}}\right)=G\left({X}_{{{\rm{He}}}}\right)-{X}_{{{\rm{He}}}}G\left(1\right)-\left(1-{X}_{{{\rm{He}}}}\right)G(0).$$ (2) The values of _x_1 and

_x_2 can be determined by applying a common tangent construction to \(\Delta G\). The calculated \(\Delta G\) and corresponding comparison with results of the previous work13 are presented

in the “Methods” section. In Fig. 2, we present the Δ_x_ contour line diagram across various compositions of H-He mixture calculated at different pressures. The measure of immiscibility,

Δ_x_, ranges from 0.07 to 1. It is worth mentioning that Δ_x_ = 0.07 remains in a separate state, where the color part only represents the degree of separation, not the separation boundary.

At 4 and 10 Mbar, the miscibility lines (black dots) are calculated using thermodynamic integration on 64-atom systems to compare with the previous studies. Our results show a good

consistence with the results obtained using vdW-DF and non-ideal entropy13, with the miscibility lines shifting slightly toward higher temperature at the He-rich zone. This trend is more

apparent in the comparison to results obtained using the PBE and ideal entropy9. These shifts are the consequence of non-ideal effects in entropy, which lowers the demixing temperatures at

the He-poor zone13 and captures the weak proton pairs at the He-rich zone leading to a higher demixing temperature12. Furthermore, the \(\Delta x\) profiles reveal a rather complex landscape

of H-He immiscibility driven by multiple factors, including temperature, pressure, and atomic compositions. The boundary of each Δ_x_ domain is shown as a function of _X_He, which is

smoothed using a 5th order polynomial regression. These boundaries show a general trend of humping up in the middle. This suggests that the increase of temperature at a fixed pressure will

lead to H-He mixing but the temperature required in the intermediate _X_He region is higher than that for the two ends. Moreover, the temperature-induced mixing is not symmetric with respect

to _X_He – it has different gradients in He-poor and He-rich regions. As the pressure increases, the miscibility gap widens, and the immiscibility intensifies. However, the changes

gradually slow down after the system has reached 7 Mbar; the miscibility diagram obtained at 10 Mbar sees much smaller changes from the former. It is observed that the immiscibility of

hydrogen and helium is highly sensitive to helium abundance in planetary conditions. Once demixing is initiated even by a very weak miscibility gap, the accumulation of helium will

accelerate the process. When the slow sedimentation of helium reaches a sufficient concentration, it might trigger stronger separation. Figure 3a shows the miscibility diagram with current

and previously published results. The calculated immiscibility line at the protosolar helium abundance (_X_He = 0.089) is shifted to higher temperatures compared with previous theoretical

results11,12,13, and it has the closest match to the experimental estimate8. Moreover, different helium abundances used in previous calculations and experiment should be noted. The four

non-ideal adiabats are calculated using thermodynamic integration with the entropy of _s_J = 7.613 _k_B atom−1 for Jupiter23,24 and _s_S = 7.357 _k_B atom−1 for Saturn24,25, and two other

entropy values, 7.822 _k_B atom−1 and 8.030 _k_B atom−1. The calculated demixing temperature (black and white dots) increases sharply upon increasing the pressure from 0.5 Mbar to 2 Mbar,

coinciding with hydrogen metallization, and then flats out to a slow slope. This trend is the same as shown by DFT calculations12,13 and experiment8. According to the miscibility diagram and

adiabats, Jupiter would enter the miscibility gap at a pressure of 1.09 (±0.03) Mbar and a temperature of 5440 (±50) K. For Saturn, the values are 0.93 (±0.03) Mbar and 4720 ( ± 50) K. The

demixing region in Jupiter ceases at about 10,000 K when pressure is 10 Mbar, the highest pressure considered in this study, where Saturn is still in miscibility gap. The miscibility gap

inferred from structural analysis spans a very wide thermodynamic range, yet the intensity of immiscibility is relatively weak. This implies that the sedimentation of helium could be a

protracted and subtle process. In Fig. 3a, the immiscibility line calculated for _X_He = 0.089 corresponds to the protosolar helium abundance. When the pressure exceeds 1 Mbar, this

immiscibility line is notably higher than that of _X_He = 0.073, which is coherent to the positive slope of Δ_x_ in He-poor region (Fig. 2). Based on the isentropes of various surface

temperatures at protosolar helium abundance14, the calculated immiscibility line for _X_He = 0.089 suggests that the onset of H-He phase separation in Saturn is about 0.746 Gyr (Fig. S15).

Due to the early onset and gradual accumulation, Saturn might experience a significant phase separation after sufficient quantitative change, whereas Jupiter would require a much longer

time. The miscibility diagram calculated at protosolar helium abundance is provided in Fig. S15. To assess whether different functionals fundamentally influence the conclusions, we also

trained a DLP based on vdW-DF and conducted identical simulations under protosolar helium abundance, as indicated by the pink squares in Fig. 3a. It can be observed that changing the

functional still yields similar conclusions. The miscibility diagram calculated at protosolar helium abundance by vdW-DF based DLP is provided in Fig. S16. To get an energy perspective, we

calculated the \(\Delta G\) using the systems of the same size (27,000 atoms) in a zoom-in helium abundance region at five pressures (0.5, 2, 4, 7, and 10 Mbar) and temperature of 5000 K

(Fig. 3b). The zoom-in region includes planetary helium abundance (_X_He = 0.073, plus symbol) and 8 adjacent points (circles). The calculation employs the Jarzynski equality26 to a virtual

integrable system27, which enables the calculation of the Gibbs free energy in nonequilibrium states prior to demixing in large systems (see “Methods” section for calculation of \(\Delta

G\)). At 0.5 Mbar, all points on the \(\Delta G\) curve have positive curvature (second derivative), suggesting that the state is stable without demixing. At 2 Mbar, small negative

curvatures appear in Δ_G_, signaling the occurrence of demixing. The demixing region is determined by common tangent construction, which is separated into unstable and metastable regions

(red area in Fig. 3b). At 4 Mbar, although the \(\Delta G\) curve still exhibits a downward trend, the negative curvature becomes quite evident – the phase separation is prevalent. At 7 

Mbar, the \(\Delta G\) curve starts to rise, and due to the constraint\(\,\Delta G\left(0\right)\,=\,\Delta G\left(1\right)\,=\,0\) (Eq. 2), there will be at least two minima where the

system is stable. At 10 Mbar, \(\Delta G\) curve rises to positive values, indicating the presence of strong separation in the system. These results correlate the change in curvature of

\(\Delta G\) to the demixing. Phase separation initially emerges in localized regions of negative curvature, and progressively proceeds as the line shape evolves. Combining the calculated

miscibility diagram with a planet model3, and assume the interior of the planet is adiabatic, we derived the nowadays internal structures of Jupiter and Saturn (Fig. 3c). In Jupiter, the

liquid layer of hydrogen, where helium droplets condensate and rain down, is estimated to be in between 0.572 and 0.847 of the radius. The calculated starting point for this layer is very

close to the experimental estimate (0.84)8 while it has a deeper ending point than the composition-corrected experimental estimate (0.68)8. This discrepancy is due to the different behaviors

of demixing temperature predicted by theory and experiment (Fig. 3a). While the demixing temperature remains constant above 2 Mbar in experiment8, it has a slow but continuous increase in

theory. In Saturn, the separation layer is estimated to be in between 0.194 and 0.677 of the radius, indicates that nearly half of Saturn’s radial distance would have helium rain. While the

outer molecular envelope may still be well approximated as adiabatic, the deep interior is expected to depart from an adiabat due to helium rain and associated double-diffusive convection

effects. Modeling these non-adiabatic regions is crucial for accurately describing the H-He demixing in Jupiter and Saturn. Incorporating these complexities into the model can lead to

temperature shifts in the predicted hydrogen-helium separation region5,28. This work offers more precise parametric conditions for intricate planetary modeling. The ab initio H-He phase

diagram delineates the _P-T_ conditions at which helium becomes immiscible with hydrogen. It is crucial for determining the equilibrium helium abundance, influencing the depletion of the

outer molecular envelope as rainout proceeds. The helium gradient, driven by the increasing helium concentration inward influences the profile of the double-diffusive convection region. The

energy released from helium differentiation scales with the amount of helium that rains out, which is controlled by the phase diagram demarcating immiscible versus miscible regions. Through

our large-scale calculations and structural analysis, we can also deduce the density ratios from the degree of demixing. This parameter can be used to parameterize the efficiency of heat

transport by double-diffusive convection in the region where helium is raining out. We hope that our presented data can be incorporated into future planetary physics research to enhance the

evolutionary modeling of planets. DISCUSSION In this work, we present direct observations suggesting the H-He immiscibility from calculations of ab initio quality. The structural analysis on

machine learning accelerated large scale MD simulations establishes a renewed miscibility diagram of hydrogen and helium under the planetary conditions of Jupiter and Saturn. The

immiscibility of hydrogen and helium is shown to be strongly dependent on temperature, pressure, and atomic composition of the mixtures. The miscibility gap is recalculated using the Δ_G_

with the systems of the same size, which reinforces the structural analysis. Based on the results, we propose a new hypothesis for the mechanism of H-He separation in Jupiter and Saturn –

the subtle increase in helium abundance through gradual helium accumulation in the early stage of planet formation potentially initiates the H-He phase separation. This process is similar to

the formation of rain droplets in gathered clouds. The slightly larger immiscibility and lower adiabatic temperature within the interior of Saturn make it more likely to undergo this

process. In contrast, in Jupiter the speed of this pre-sedimentation is lower, along with its higher adiabatic temperature decreases the likelihood of further separation. METHODS DLP

POTENTIAL TRAINING DLP is constructed based on SCAN+rVV10. This scheme is chosen for its ‘best-of-both-worlds’ feature – SCAN captures the short and intermediate-range interactions and rVV10

describes the long-range vdW interaction29. The combination of SCAN+rVV10 has been shown to have accuracy better than 1 kcal/mol for benchmarking the MP2 + ∆CCSD(T) results and be able to

reproduce subtle features of the potential energy surface30. In addressing gigapascal pressures, SCAN+rVV10 has been shown to improve the agreement with experimental data compared with PBE

and SCAN on liquid-liquid phase transition in nitrogen-oxygen mixtures, thanks to the revised noncovalent interactions of short- to long-range van der Waals interaction31. A similar better

match appears in the works of Tantardini32 and Anh33. Regarding metallic materials, SCAN+rVV10 has been demonstrated to align well with experimental results34, outperforming PBE, PBEsol, and

SCAN. Notably, the addition of rVV10 to SCAN yields a highly accurate method for diversely bonded systems. We use the projector augmented wave35 (PAW) approach for electron-ion interaction,

specifically the hard PAW pseudopotentials for hydrogen (H_h, 06Feb2004), and helium (two valence electrons, He 05Jan2001) as provided with VASP. The DFT convergence tests are provided in

Supplementary Information Note 1. We use Deep Potential Generator (DP-GEN)36 to cover temperatures from 1000 K to 13,000 K and pressures from 0.5 Mbar to 10 Mbar. 28 helium proportions

evenly distributed from 0 to 1 were used as training set. A total of 106,817 configurations in the training set, and trained for 12,000,000 steps. EXPLAINATION OF RE-WEIGHTING PROCEDURE In

Fig. 1, after assigning a \({P}_{{{{\rm{cond}}}}}\) value to each atom, we obtained atoms in their respective enriched regions (iv, \({P}_{{{{\rm{cond}}}}} \, > \, 0.5\)) and atoms in

another element-enriched regions (iv, \({P}_{{{{\rm{cond}}}}} \, < \, 0.5\)). This way, we preliminary classified all atoms; however, atoms belonging to the same category may not have

equal contributions to their respective regions. For instance, in the 6000 K case in iii, a helium atom with \({P}_{{{{\rm{cond}}}}}=0.9\) contributes differently than a helium atom with

\({P}_{{{{\rm{cond}}}}}=0.55\) (although both values are larger than 0.5). While the latter also appears in the helium-rich region, it behaves more like an interface atom. Therefore, it is

necessary to further consider their weights when counting. To address this issue, we re-weighted atoms in the self-enrichment zone by the corresponding frequency of \({P}_{{{{\rm{cond}}}}}\)

values. And because of the causal links between two types of atoms within the same zone, the weights of atoms in the alternative element-enriched area should also be defined by the relevant

\({P}_{{{{\rm{cond}}}}}\) values frequency of the enriched type. As a result, the weights of all atoms exhibit central symmetry, as depicted by the gray shadows in iii. The effect of the

re-weight procedure is visibly noticeable in iv, the weights of atoms at the interface decrease, while the weights of the genuine “shapers” are intentionally retained. In Fig. S13, the

values of \(\Delta x\) corresponding to helium abundances ranging from pure hydrogen to pure helium under the condition of complete mixing are provided, which should be 0 under the

assumption of perfect ideal conditions for infinite systems. However, in practical applications involving finite systems, where the dividing lines of \({X}_{{{\rm{He}}}}\) and

\(1-{X}_{{{\rm{He}}}}\) are not sharp but have a Gaussian width, this introduces a certain level of uncertainty to the values of \({x}_{1}\) and \({x}_{2}\). As a result, achieving a fully

mixed \(\Delta x\) value of 0 is not feasible. And it can be observed that the order closer to the ideal scenario is (27,000 atoms, re-weight) > (27,000 atoms) > (9604 atoms,

re-weight) > (9604 atoms). To enhance accuracy, it is necessary to employ as many atoms as possible and subsequently reweight the results of conditional probabilities. IMMISCIBILITY

DEFINITION Through re-weighted conditional probability classification, we gain the \(\Delta x{=x}_{2}-{x}_{1}\) as a measure of immiscibility. Under the assumption of perfect ideal

conditions for infinite systems, ∆_x_ = 0 means fully mixed and ∆_x_ = 1 means fully separated. However, in practical applications involving finite systems, when the atom number of the

system is once fixed, the fully mixed \(\Delta x\) value is typically correlated with \({X}_{{{\rm{He}}}}\) of the system. Thus, the individualized immiscibility limits are considered for

different \({X}_{{{\rm{He}}}}\) in this paper. The distribution of \(\Delta x\) corresponds to normal distribution through the Kolmogorov-Smirnov test37, as shown in Fig. S13. 20,000 frames

for each concerned composition were built randomly to better construct normal distributions of the mixed system. At the same time, trajectories from MD simulation were randomly shuffled to

achieve a mixed system. From the analysis of these coordinates, the mean value of \(\Delta x\) shows to be at largest at _X_He = 0.5 and accelerates approaching 0 with the system’s

composition tends to be extreme. To explain this, we have to first assume a non-ideal mixed system with helium fraction \({X}_{{{\rm{He}}}}\) which is misjudged separated because of random

fluctuation. The fake _x_1 and _x_1 will be very close and are symmetric concerning the \({X}_{{{\rm{He}}}}\): $${x}_{1}{+{{\rm{x}}}}_{2}={2X}_{{{\rm{He}}}},$$ (3) the definition of \(\Delta

x\) is: $$\Delta {x}_{{{{\rm{mix}}}}}={x}_{2}{-x}_{1},$$ (4) with Eq. 1 in the main text, and introduce fluctuation variable _δ_ for non-ideal mixed system, it can be derived that: $$\Delta

{x}_{{{{\rm{mix}}}}}

=2{x}_{2}-2{X}_{{{{\rm{He}}}}}=\frac{2}{1+\frac{N({{{{\rm{H}}}}}_{{{{\rm{He}}}}-{{{\rm{rich}}}}})}{N({{{{\rm{He}}}}}_{{{{\rm{He}}}}-{{{\rm{rich}}}}})}}-2{X}_{{{{\rm{He}}}}}\\

\,=\,\left\{\begin{array}{c} \displaystyle \frac{2}{1+\frac{1-{X}_{{{{\rm{He}}}}}}{{X}_{{{{\rm{He}}}}}}}-2{X}_{{{{\rm{He}}}}}=0,\,{{{\rm{ideal}}}}\,{{{\rm{mix}}}} \hfill \\ \displaystyle

\frac{2}{1+\frac{1-{X}_{{{{\rm{He}}}}}}{{X}_{{{{\rm{He}}}}}}\delta }-2{X}_{{{{\rm{He}}}}},\,{{{\rm{non}}}}-{{{\rm{ideal}}}}\,{{{\rm{mix}}}}\end{array}\right..$$ (5) In ideal mixed system,

the He-rich region completely overlaps with the He-poor region, thus \(\frac{({N{{{\rm{H}}}}}_{{{\rm{He}}}-{{{\rm{rich}}}}})}{N({{{\rm{He}}}}_{{{\rm{He}}}-{{{\rm{rich}}}}})}\) is equal to

\(\frac{1-{X}_{{{\rm{He}}}}}{{X}_{{{\rm{He}}}}}\). For non-ideal mixed systems, a fluctuation variable _δ_ has to be added to thought. The dashed function fit lines are from Eq. 5 with

proper _δ_ to fit random and shuffled trajectories. The relation between \(\overline{\Delta {x}_{{{\rm{mix}}}}}\) and \({X}_{{{\rm{He}}}}\) of mixed system indicates different demixing

definitions for systems with different compositions. The inset of Fig. S13 shows the distribution of \(\Delta {x}_{{{\rm{mix}}}}\) of 27,000 atoms. Because \(\Delta {x}_{{{\rm{mix}}}}\) of

mixed system abides by normal distribution, the three-sigma rule38 is used to determine the mixed interval. μ is the mean of the distribution, and σ is its standard deviation. _μ_ + 3_σ_

values are marked in the figure as the right boundaries of intervals, in which 99.865% of \(\Delta {x}_{{{\rm{mix}}}}\) will lie in. Now, with a \(\Delta x\) calculated from a single frame

or a \(\overline{\Delta x}\) from a trajectory containing many frames, we can compare it with _μ_ + 3_σ_ value of \(\Delta {x}_{{{\rm{mix}}}}\). If the value falls outside of the 3_σ_

intervals, indicating that the probability of the system being mixed is only 0.135%. For a \(\Delta x\) from a single frame, it risks with these small probabilities, but for a

\(\overline{\Delta x}\) from a long trajectory of a wide time scale containing many frames, this risk can be further reduced. Thus, the _μ_ + 3_σ_ values of \(\Delta {x}_{{{\rm{mix}}}}\) for

every \({X}_{{{\rm{He}}}}\) are used as the immiscibility threshold. The inset of Fig. S13 presents an application example of immiscibility definition at helium abundance of Jupiter and

Saturn along 2 Mbar. With the decrease of the temperature, the distributions of \(\Delta x\) gradually move away from the top mixed one. When the mean value of the distribution is greater

than _μ_ + 3_σ_, the system is determined to be in the miscibility gap, and so is Jupiter and Saturn’s case at this thermodynamic condition. The _μ_ + 3_σ_ values for systems with 9604 and

27,000 atoms at certain helium abundances are shown in the Table. S1. For binary systems, the results are symmetric around a helium abundance of 0.5, so we only list the values for _X_He 

< 0.5. MD SIMULATIONS DETAILS To analyze through structures of trajectories, two sets of MD simulations were performed. The first set of calculations was carried out on 9604-atom

simulation boxes at compositions of _X_He = 0.2, 0.357, 0.5, 0.643, 0.8, and 0.9. The second set of calculations was performed on 27,000-atom simulation boxes with helium abundances of 0.073

(representing the current helium abundance in the atmospheres of Jupiter and Saturn) and 0.089 (representing protosolar helium abundance, which is the theoretical initial helium abundance

of Jupiter and Saturn). For each \({X}_{{{\rm{He}}}}\), calculations were performed at pressures of 0.5, 1, 1.5, 2, 4, 7, and 10 Mbar, and necessary temperatures ranging from 1000 to 13,000 

K. A timestep of 0.2 fs was used and 100 ps MD simulation was performed to guarantee a long enough trajectory for full sample after thermodynamic equilibrium was reached. \(\DELTA G\)

CALCULATION OF LARGE SYSTEMS When calculating the \(\Delta G\), it is crucial to ensure that the term \(G({X}_{{{\rm{He}}}})\) in Eq. 2 refers to a fully mixed system, as it cannot represent

the Gibbs free energy of a specific helium abundance \({X}_{{{\rm{He}}}}\) once the system undergoes demixing. Unlike previous work using small enough systems to guarantee the fully mixed

state in simulation, we apply the Jarzynski equality to a virtual integrable system27 to calculate the Gibbs free energy difference of volume change: $${{{{\rm{e}}}}}^{{-\beta}

({G}_{{{{\rm{B}}}}}-{G}_{{{{\rm{A}}}}})}={({V}_{{{{\rm{B}}}}}/{V}_{{{{\rm{A}}}}})}^{N} {\big\langle {{\rm{e}}}^{{\beta}

({U}_{{{{\rm{A}}}}}-{U}_{{{{\rm{B}}}}})}\big\rangle}_{{{\rm{A}}},{{\bf{r}}}}.$$ (6) Here \(\beta=1/{k}_{{{\rm{B}}}}T\) is the inverse temperature, \({U}_{{{{\rm{A}}}}}\) and

\({U}_{{{{\rm{B}}}}}\) are the interaction energies of the initial state A and final state B, and the corresponding Gibbs free energies are represented by \({G}_{{{\rm{A}}}}\) and

\({G}_{{{\rm{B}}}}\). \({{{\bf{r}}}}\,=\,(x,\,y,\,z)\) is the coordinates of all particles at initial state A. Since the transition from state A to state B only involves a change in volume

while keeping the relative positions of all atoms unchanged, the miscibility of state B remains the same as state A. Therefore, we only need to locate a completely mixed point on the

miscibility diagram to be the initial state A. The calculation process for determining \(\Delta G\) in large systems can be roughly divided into two steps: 1) calculating the absolute free

energy of completely mixed state A using thermodynamic integration, and 2) repeatedly sampling the process from state A to state B using Eq. 6. After considering structure analysis results

of this work, as well as multiple previous calculations and experiments, we have selected the temperature of 5000 K and pressure of 0.5 Mbar as the condition for the completely mixed state

A, which can be observed from Fig. 2. Moreover, at a temperature of 5000 K, there will be various degrees of phase separation occurring with changes in pressure. This provides an opportunity

for more comprehensive validation. To compare systems with different particle compositions, it is necessary to calculate the absolute Gibbs free energy. We calculate the absolute Gibbs free

energy of 11 components (including pure hydrogen, pure helium, and components near the abundance of planetary helium) at the initial state of 0.5 Mbar and 5000 K through thermodynamic

integration (TI). As a result, we obtain the \(\Delta G\) curve representing the completely mixed initial state. Then we apply the Jarzynski equality to a virtual integrable system27 to

calculate Gibbs free energy difference of volume change and repeat this sampling process to calculate the ensemble average on the right-hand side of Eq. 6. Subsequently, several sets of

\(\Delta G\) calculations were performed near the planetary helium abundance along the 5000 K isotherm in the _P-T_ miscibility diagram. By taking the logarithm of both sides of Eq. 6 and

considering values per atom, we obtain: $$\frac{-({G}_{{{{\rm{B}}}}} \! -{ \! G}_{{{{\rm{A}}}}})}{N}\approx \frac{{V}_{{{{\rm{B}}}}}}{\beta {V}_{{{{\rm{A}}}}}}+\frac{{\langle

{U}_{{{{\rm{A}}}}} \! -{ \! U}_{{{{\rm{B}}}}}\rangle }_{{{{\rm{A}}}},{{{\bf{r}}}}}}{N}.$$ (7) Thus, the error of \(\Delta G\) mainly comes from two parts: \(\frac{{V}_{{{\rm{B}}}}}{{\beta

V}_{{{\rm{A}}}}}\) and \(\frac{{\langle {U}_{{{{\rm{A}}}}}-{U}_{{{{\rm{B}}}}}\rangle }_{{{{\rm{A}}}},{{{\bf{r}}}}}}{N}\). Generally, our calculations are based on a system with 27,000 atoms,

however, when relaxing to obtain the \({V}_{B}\), we employ larger system sizes. 729,000 atoms systems are used for 2 Mbar while 216,000 atoms systems are used for 4 Mbar, 7 Mbar, and 10 

Mbar. These larger sizes are 27 times and 8 times the target size, respectively (The smaller the pressure, the larger the volume with fluctuation, and the greater the number of atoms

required to maintain consistent accuracy). In this way, \(\frac{{{{\rm{\sigma }}}}(\overline{{V}_{{{{\rm{B}}}}}})}{\beta {V}_{{{{\rm{A}}}}}}\) is of the order 10−7 eV. As for repeated

sampling of \(\frac{{\langle {U}_{{{{\rm{A}}}}}-{U}_{{{{\rm{B}}}}}\rangle }_{{{{\rm{A}}}},{{{\bf{r}}}}}}{N}\), we have performed convergency test on sample size to ensure the accuracy. See

Fig. S14. \(\DELTA G\) CALCULATION BY THERMODYNAMIC INTEGRATION Thermodynamic integration calculates differences in free energies to get the target state (described by potential \({U}_{1}\))

free energy integrated from a reference state (described by potential \({U}_{0}\)). Then define a transition potential with two switching functions _f_(λ), _g_(λ): $$U=f({{{\rm{\lambda

}}}}){U}_{0}+g({{{\rm{\lambda }}}}){U}_{1},$$ (8) _f_(λ) and _g_(λ) satisfy _f_(0) = 1\(,\) _f_(1) = 0\(,\) _g_(0) = 0 and _g_(0) = 1. With Eq. 8, the free energy _F_ is a function of _N_,

_V_, _T_, \(\lambda\): $$F(\lambda )=-{k}_{{{{\rm{B}}}}}T\,{{\mathrm{ln}}}\,Q(N,V,T,\lambda ).$$ (9) The difference between the target state and reference state is: $${F}_{1}-{F}_{0}={\int

}_{0}^{1}d{\lambda }\frac{\partial F({\lambda })}{\partial {\lambda }},$$ (10) And it can be derived that \(\frac{\partial F(\lambda )}{\partial \lambda }={\langle \frac{\partial U}{\partial

\lambda }\rangle }_{\lambda }\) through: $$\frac{\partial F(\lambda )}{\partial \lambda } =\, -\frac{{k}_{{{{\rm{B}}}}}T}{Q}\frac{\partial Q(\lambda )}{\partial \lambda

}\,=-\frac{{k}_{{{{\rm{B}}}}}T}{Z}\frac{\partial Z(\lambda )}{\partial \lambda }\,=-\frac{{k}_{{{{\rm{B}}}}}T}{Z}\frac{\partial }{\partial \lambda }\int

{{{{\rm{d}}}}}^{N}{{{\bf{r}}}}{{{{\rm{e}}}}}^{-\beta U({{{{\bf{r}}}}}_{{{{\bf{1}}}}}{{{,}}}{{\ldots }}{{{,}}}{{{{\bf{r}}}}}_{{{{\boldsymbol{N}}}}},\lambda )}\\ =\,

-\frac{{k}_{{{{\rm{B}}}}}T}{Z}\int {d}^{N}{{{\bf{r}}}}(-\beta \frac{\partial U}{\partial \lambda }){{{{\rm{e}}}}}^{-\beta U({{{{\bf{r}}}}}_{{{{\bf{1}}}}}{{{,}}}{{\ldots

}}{{{,}}}{{{{\bf{r}}}}}_{{{{\boldsymbol{N}}}}},\lambda )}=\frac{\int {d}^{N}{{{\bf{r}}}}(\frac{\partial U}{\partial \lambda }){{{{\rm{e}}}}}^{-{{{\rm{\beta

}}}}U({{{{\bf{r}}}}}_{{{{\bf{1}}}}}{{{,}}}{{\ldots }}{{{,}}}{{{{\bf{r}}}}}_{{{{\boldsymbol{N}}}}},\lambda )}}{\int {d}^{N}{{{\bf{r}}}}{{{{\rm{e}}}}}^{-\beta

U({{{{\bf{r}}}}}_{{{{\bf{1}}}}}{{{,}}}{{\ldots }}{{{,}}}{{{{\bf{r}}}}}_{{{{\boldsymbol{N}}}}},\lambda )}}={\left\langle \frac{\partial U}{\partial \lambda }\right\rangle }_{\lambda }.$$ (11)

Replace the corresponding term in Eq. 10: $${F}_{1}-{F}_{0}={\int }_{0}^{1}d\lambda {\left \langle \frac{\partial U}{\partial \lambda }\right \rangle }_{\lambda }.$$ (12) From Eq. 12, the

difference of free energies can be obtained by integrating transition potential _U_ defined by Eq. 8. In this work, the switching functions were defined as _f_(_λ_) = 1 – _λ_ and _g_(_λ_) =

_λ_. Ideal gas was taken as the reference state. An Ideal gas is chosen as the reference state, from which we integrate 47 coupling constants10 to the target state with DLP potential. After

a relaxation of at least 50 ps of a certain composition to thermodynamic equilibrium, all MD simulations for integration have been run for 2 ps with a timestep of 0.2 fs. 7 temperatures and

64 compositions are considered for the two concerned pressures. We also calculate two \(\Delta G\) lines with vdW-DF and SCAN+rVV10 using VASP39 code at 10 Mbar and 5000 K to benchmark the

comparison. With free energy _F_ of a certain fully mixed system with helium fraction \({X}_{{{\rm{He}}}}\), the delta Gibbs free energy can be calculated by Eq. 2, where

\(G({X}_{{{\rm{He}}}})=F({X}_{{{\rm{He}}}})+{PV}\). Then, by applying a common double tangent construction to ∆_G_ we can determine _x_1 and _x_2. According to Fig. S17, the negative

curvature of ∆_G_ gradually disappear with the increase of temperature. DATA AVAILABILITY The relevant data for this research is accessible on Zenodo40. Source data are provided with this

paper. CODE AVAILABILITY LAMMPS and DeepMD are free and open source codes available at https://lammps.sandia.gov and http://www.deepmd.org, respectively. VASP is a commercial code available

from https://www.vasp.at. Detailed instructions for obtaining and using these codes can be found on their respective websites. The AMA-RCP code is provided on Zenodo40. REFERENCES * Guillot,

T., Gautier, D., Spohn, T. & Schubert, J. _Treatise of Geophysics_. Vol. 10. (_Planets and Moons_, Elsevier, 2007). * Stevenson, D. J. Metallic helium in massive planets. _Proc. Natl

Acad. Sci. USA_ 105, 11035–11036 (2008). Article  ADS  CAS  PubMed Central  Google Scholar  * Helled, R. _The interiors of Jupiter and Saturn_. Vol. 175 (Oxford Univ. Press, 2018). * von

Zahn, U., Hunten, D. M. & Lehmacher, G. Helium in Jupiter’s atmosphere: results from the Galileo probe Helium Interferometer Experiment. _J. Geophys. Res. Planets_ 103, 22815–22829

(1998). Article  ADS  Google Scholar  * Mankovich, C. R. & Fortney, J. J. Evidence for a dichotomy in the interior structures of Jupiter and Saturn from helium phase separation.

_Astrophys. J._ 889, 51 (2020). Article  ADS  CAS  Google Scholar  * Bahcall, J. N., Pinsonneault, M. H. & Wasserburg, G. J. Solar models with helium and heavy-element diffusion. _Rev.

Mod. Phys._ 67, 781–808 (1995). Article  ADS  CAS  Google Scholar  * Fortney, J. J. Looking into the giant planets. _Science_ 305, 1414–1415 (2004). Article  CAS  PubMed  Google Scholar  *

Brygoo, S. et al. Evidence of hydrogen−helium immiscibility at Jupiter-interior conditions. _Nature_ 593, 517–521 (2021). Article  ADS  CAS  PubMed  Google Scholar  * Lorenzen, W., Holst, B.

& Redmer, R. Demixing of hydrogen and helium at megabar pressures. _Phys. Rev. Lett._ 102, 115701 (2009). Article  ADS  PubMed  Google Scholar  * Morales, M. A. et al. Phase separation

in hydrogen—helium mixtures at Mbar pressures. _Proc. Natl Acad. Sci. USA_ 106, 1324–1329 (2009). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Lorenzen, W., Holst, B. &

Redmer, R. Metallization in hydrogen-helium mixtures. _Phys. Rev. B_ 84, 235109 (2011). * Morales, M. A., Hamel, S., Caspersen, K. & Schwegler, E. Hydrogen-helium demixing from first

principles: From diamond anvil cells to planetary interiors. _Phys. Rev. B_ 87, 174105 (2013). Article  ADS  Google Scholar  * Schöttler, M. & Redmer, R. Ab initio calculation of the

miscibility diagram for hydrogen-helium mixtures. _Phys. Rev. Lett._ 120, 115703 (2018). Article  ADS  PubMed  Google Scholar  * Püstow, R., Nettelmann, N., Lorenzen, W. & Redmer, R.

H/He demixing and the cooling behavior of Saturn. _Icarus_ 267, 323–333 (2016). Article  ADS  Google Scholar  * Zeng, Q. et al. Ab initio validation on the connection between atomistic and

hydrodynamic description to unravel the ion dynamics of warm dense matter. _Phys. Rev. Res._ 3, 033116 (2021). Article  CAS  Google Scholar  * Zhang, L., Wang, H., Car, R. & E, W. Phase

diagram of a deep potential water model. _Phys. Rev. Lett._ 126, 236001 (2021). Article  ADS  CAS  PubMed  Google Scholar  * Niu, H., Yang, Y. I. & Parrinello, M. Temperature dependence

of homogeneous nucleation in ice. _Phys. Rev. Lett._ 122, 245501 (2019). Article  ADS  CAS  PubMed  Google Scholar  * Wang, H., Zhang, L., Han, J. & E, W. DeePMD-kit: A deep learning

package for many-body potential energy representation and molecular dynamics. _Comput. Phys. Commun._ 228, 178–184 (2018). Article  ADS  CAS  Google Scholar  * Sun, J., Ruzsinszky, A. &

Perdew, J. P. Strongly constrained and appropriately normed semilocal density functional. _Phys. Rev. Lett._ 115, 036402 (2015). Article  ADS  PubMed  Google Scholar  * Sabatini, R., Gorni,

T. & de Gironcoli, S. Nonlocal van der Waals density functional made simple and efficient. _Phys. Rev. B_ 87, 041108 (2013). Article  ADS  Google Scholar  * Brandani, G. B. et al.

Quantifying disorder through conditional entropy: an application to fluid mixing. _PLoS ONE_ 8, e65617 (2013). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Plimpton, S. Fast

parallel algorithms for short-range molecular dynamics. _J. Comput. Phys._ 117, 1–19 (1995). Article  ADS  CAS  Google Scholar  * Militzer, B. Equation of state calculations of

hydrogen-helium mixtures in solar and extrasolar giant planets. _Phys. Rev. B_ 87, 014202 (2013). * Schöttler, M. & Redmer, R. Simulations of H–He mixtures using the van der Waals

density functional. _J. Plasma Phys_. 84, 755840401 (2018). * Hubbard, W. B. & Militzer, B. A preliminary Jupiter model. _Astrophys. J._ 820, 80 (2016). Article  ADS  Google Scholar  *

Jarzynski, C. Nonequilibrium equality for free energy differences. _Phys. Rev. Lett._ 78, 2690–2693 (1997). Article  ADS  CAS  Google Scholar  * Zhu, L. & Wang, J. Calculating the free

energy difference by applying the Jarzynski equality to a virtual integrable system. _Phys. Rev. E_ 98, 022117 (2018). Article  ADS  MathSciNet  CAS  PubMed  Google Scholar  * Mankovich, C.,

Fortney, J. J. & Moore, K. L. Bayesian evolution models for jupiter with helium rain and double-diffusive convection. _Astrophys. J._ 832, 113 (2016). Article  ADS  Google Scholar  *

Peng, H., Yang, Z.-H., Perdew, J. P. & Sun, J. Versatile van der Waals density functional based on a meta-generalized gradient approximation. _Phys. Rev. X_ 6, 041005 (2016). Google

Scholar  * Chen, J., Xie, W., Li, K., Zhang, S. & Sun, Y.-Y. Benchmarking PBE+D3 and SCAN+rVV10 methods using potential energy surfaces generated with MP2+ΔCCSD(T) calculation. _Chin.

Phys. B_ 29, 013102 (2020). Article  ADS  CAS  Google Scholar  * Lan, Y.-S. et al. First-order liquid-liquid phase transition in nitrogen-oxygen mixtures. _Phys. Rev. B_ 103, 144105 (2021).

Article  ADS  CAS  Google Scholar  * Tantardini, C., Jalolov, F. N. & Kvashnin, A. G. Crystal structure evolution of fluorine under high pressure. _J. Phys. Chem. C_ 126, 11358–11364

(2022). Article  CAS  Google Scholar  * Anh, L. T., Wada, M., Fukui, H., Kawatsu, T. & Iitaka, T. First-principles calculations of the epsilon phase of solid oxygen. _Sci. Rep._ 9, 8731

(2019). Article  ADS  PubMed  PubMed Central  Google Scholar  * Patra, A., Bates, J. E., Sun, J. & Perdew, J. P. Properties of real metallic surfaces: effects of density functional

semilocality and van der Waals nonlocality. _Proc. Natl Acad. Sci. USA_ 114, E9188–E9196 (2017). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Kresse, G. & Joubert, D.

From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 59, 1758–1775 (1999). Article  ADS  CAS  Google Scholar  * Zhang, Y. et al. DP-GEN: a concurrent

learning platform for the generation of reliable deep learning based potential energy models. _Comput. Phys. Commun._ 253, 107206 (2020). Article  MathSciNet  CAS  Google Scholar  * Massey,

F. J. The Kolmogorov-Smirnov test for goodness of fit. _J. Am. Stat. Assoc._ 46, 68–78 (1951). Article  Google Scholar  * Grafarend, E. W. _Linear and Nonlinear Models: Fixed Effects, Random

Effects, and Mixed Models_ (Berlin; New York: Walter de Gruyter, 2006). * Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a

plane-wave basis set. _Phys. Rev. B_ 54, 11169–11186 (1996). Article  ADS  CAS  Google Scholar  * Chang, X. et al. Theoretical evidence of H-He demixing under Jupiter and Saturn

conditions.H-He-Immiscibility (v1.0.1). _Zenodo_ https://doi.org/10.5281/zenodo.13358797 (2024). * Nettelmann, N., Püstow, R. & Redmer, R. Saturn layered structure and homogeneous

evolution models with different EOSs. _Icarus_ 225, 548–557 (2013). Article  ADS  CAS  Google Scholar  * Celliers, P. M. et al. Insulator-metal transition in dense fluid deuterium. _Science_

361, 677–682 (2018). Article  ADS  CAS  PubMed  Google Scholar  * Lu, B., Kang, D., Wang, D., Gao, T. & Dai, J. Towards the same line of liquid–liquid phase transition of dense hydrogen

from various theoretical predictions*. _Chin. Phys. Lett._ 36, 103102 (2019). Article  ADS  CAS  Google Scholar  * Pierleoni, C., Morales, M. A., Rillo, G., Holzmann, M. & Ceperley, D.

M. Liquid–liquid phase transition in hydrogen by coupled electron–ion Monte Carlo simulations. _Proc. Natl Acad. Sci. USA_ 113, 4953–4957 (2016). Article  ADS  CAS  PubMed  PubMed Central 

Google Scholar  * Knudson, M. D. et al. Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium. _Science_ 348, 1455–1460 (2015). Article  ADS  CAS  PubMed 

Google Scholar  * Zaghoo, M. & Silvera, I. F. Conductivity and dissociation in liquid metallic hydrogen and implications for planetary interiors. _Proc. Natl Acad. Sci. USA_ 114,

11873–11877 (2017). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Ohta, K. et al. Phase boundary of hot dense fluid hydrogen. _Sci. Rep._ 5, 16560 (2015). Article  ADS  CAS 

PubMed  PubMed Central  Google Scholar  Download references ACKNOWLEDGEMENTS This work was supported by the National Key R&D Program of China under Grant No. 2017YFA0403200, the National

Natural Science Foundation of China under Grant Nos., No. 12047561, and No. 12104507, the NSAF under Grant No. U1830206, the Science and Technology Innovation Program of Hunan Province

under Grant No. 2021RC4026. We acknowledge Jianmin Yuan, Yexin Feng, and Rei Cao for their helpful discussions. Jupiter image (PIA02852) and Saturn image (PIA08415) in Fig. 3 are provided by

NASA/JPL-Caltech. AUTHOR INFORMATION Author notes * These authors contributed equally: Xiaoju Chang, Bo Chen. AUTHORS AND AFFILIATIONS * College of Science, National University of Defense

Technology, Changsha, China Xiaoju Chang, Bo Chen, Qiyu Zeng, Kaiguo Chen, Qunchao Tong, Xiaoxiang Yu, Dongdong Kang, Shen Zhang, Fangyu Guo, Yong Hou, Zengxiu Zhao & Jiayu Dai * Hunan

Key Laboratory of Extreme Matter and Applications, National University of Defense Technology, Changsha, China Xiaoju Chang, Bo Chen, Qiyu Zeng, Kaiguo Chen, Qunchao Tong, Xiaoxiang Yu, 

Dongdong Kang, Shen Zhang, Fangyu Guo, Yong Hou, Zengxiu Zhao & Jiayu Dai * Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing, P.

R. China Han Wang * Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Saskatchewan, Canada Yansun Yao * State Key Lab of Superhard Materials and

International Center for Computational Method and Software, College of Physics, Jilin University, Changchun, China Yanming Ma * International Center of Future Science, Jilin University,

Changchun, China Yanming Ma Authors * Xiaoju Chang View author publications You can also search for this author inPubMed Google Scholar * Bo Chen View author publications You can also search

for this author inPubMed Google Scholar * Qiyu Zeng View author publications You can also search for this author inPubMed Google Scholar * Han Wang View author publications You can also

search for this author inPubMed Google Scholar * Kaiguo Chen View author publications You can also search for this author inPubMed Google Scholar * Qunchao Tong View author publications You

can also search for this author inPubMed Google Scholar * Xiaoxiang Yu View author publications You can also search for this author inPubMed Google Scholar * Dongdong Kang View author

publications You can also search for this author inPubMed Google Scholar * Shen Zhang View author publications You can also search for this author inPubMed Google Scholar * Fangyu Guo View

author publications You can also search for this author inPubMed Google Scholar * Yong Hou View author publications You can also search for this author inPubMed Google Scholar * Zengxiu Zhao

View author publications You can also search for this author inPubMed Google Scholar * Yansun Yao View author publications You can also search for this author inPubMed Google Scholar *

Yanming Ma View author publications You can also search for this author inPubMed Google Scholar * Jiayu Dai View author publications You can also search for this author inPubMed Google

Scholar CONTRIBUTIONS J.D. designed the project. X.C., B.C. Y.Y. Y.M. and J.D. suggested the specific scientific problem and the general idea on methodology, X.C. and J.D. performed the MD

simulations and analyzed data, K.C. and H.W. contributed to the AMA-RCP. X.C., B.C., Q.Z., H.W., K.C., Q.T., X.Y., D.K., S.Z., F.G., Y.H., Z.Z, Y.Y, Y.M, and J. D. interpreted the results,

X.C., Y.Y., and J.D. wrote the paper, and Q.Z., H.W., K.C., Q.T., X.Y. D.K., Y.H., Z.Z, Y.M, X.Y., S.Z., and D.K. edited the manuscript before submission. CORRESPONDING AUTHORS

Correspondence to Yansun Yao, Yanming Ma or Jiayu Dai. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing interests. PEER REVIEW PEER REVIEW INFORMATION _Nature

Communications_ thanks Nadine Nettelmann and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available. ADDITIONAL INFORMATION

PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION

PEER REVIEW FILE SOURCE DATA SOURCE DATA RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International

License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the

source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived

from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line

to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will

need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. Reprints and permissions ABOUT THIS

ARTICLE CITE THIS ARTICLE Chang, X., Chen, B., Zeng, Q. _et al._ Theoretical evidence of H-He demixing under Jupiter and Saturn conditions. _Nat Commun_ 15, 8543 (2024).

https://doi.org/10.1038/s41467-024-52868-4 Download citation * Received: 10 October 2023 * Accepted: 19 September 2024 * Published: 02 October 2024 * DOI:

https://doi.org/10.1038/s41467-024-52868-4 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative