Theoretical evidence of H-He demixing under Jupiter and Saturn conditions

Domain: earth_science
Published: 2024-10-02
Editorial outcome: major_revision

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.

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–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 (Δ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 of H-He, which is a non-equilibrium process, is also a key factor 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 non-equilibrium approach for an improved evaluation of ΔG in large systems (~27,000 atoms), and use it to obtain the miscibility.

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 (Pcond) 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 XHe = 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 Pcond for each atom in the system based on the atomic composition of the N nearest neighbors. All atoms are categorized by their Pcond with a weight assigned for each category. With the categorized and reweighted results (Fig. 1a), x1 and x2 can be calculated directly, as the relation of physical separation will yield slightly different results on the other: x1 = (NHe,poor) / (NHe,poor + NH,poor); x2 = (NHe,rich) / (NHe,rich + NH,rich). (1)

The value Δx = x2 − x1 represent 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 Δx). Figure 1b shows three representative examples with XHe = 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 Pcond 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 XHe = 0.073 and 0.089 calculated with 27,000-atom simulation boxes and XHe = 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 ΔG with 64-atom systems as the equation, ΔG(XHe) = G(XHe) − XHeG(1) − (1 − XHe)G(0). (2)

The values of x1 and x2 can be determined by applying a common tangent construction to ΔG. The calculated Δ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 Δ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 XHe, 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 XHe region is higher than that for the two ends. Moreover, the temperature-induced mixing is not symmetric with respect to XHe – 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 (XHe = 0.089) is shifted to higher temperatures compared with previous theoretical results11–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 = 7.613 kB atom−1 for Jupiter23,24 and s = 7.357 kB atom−1 for Saturn24,25, and two other entropy values, 7.822 kB atom−1 and 8.030 kB 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 XHe = 0.089 corresponds to the protosolar helium abundance. When the pressure exceeds 1 Mbar, this immiscibility line is notably higher than that of XHe = 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 XHe = 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 Δ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 (XHe = 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 ΔG). At 0.5 Mbar, all points on the Δ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 ΔG curve still exhibits a downward trend, the negative curvature becomes quite evident – the phase separation is prevalent. At 7 Mbar, ΔG curve starts to rise, and due to the constraint ΔG(0) = ΔG(1) = 0 (Eq. 2), there will be at least two minima where the system is stable. At 10 Mbar, ΔG curve rises to positive values, indicating the presence of strong separation in the system. These results correlate the change in curvature of Δ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 rain out 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.

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, 06 Feb 2004), and helium (two valence electrons, He05 Jan 2001) 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 Pcond value to each atom, we obtained atoms in their respective enriched regions (iv, Pcond > 0.5) and atoms in another element-enriched regions (iv, Pcond < 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 Pcond = 0.9 contributes differently than a helium atom with Pcond = 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 Pcond 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 areas should also be defined by the relevant Pcond 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 Δ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 XHe and 1 − XHe are not sharp but have a Gaussian width, this introduces a certain level of uncertainty to the values of x1 and x2. As a result, achieving a fully mixed Δ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 Δx = x2 − x1 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 Δx value is typically correlated with XHe of the system. Thus, the individualized immiscibility limits are considered for different XHe in this paper. The distribution of Δ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 Δx shows to be at largest at XHe = 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 XHe which is misjudged separated because of random fluctuation. The fake x1 and x2 will be very close and are symmetric concerning the XHe: x1 + x2 = 2XHe, (3)

When calculating the Δx in the main text, and introduce fluctuation variable δ for non-ideal mixed system, it can be derived that: Δxmix = 2 |x2 − XHe| = |1 + N(He,rich) − N(He,poor) / (N(He,rich) + N(He,poor))| − 2XHe = 0, ideal mix (5)

In ideal mixed system, the He-rich region completely overlaps with the He-poor region, thus N(He,rich) is equal to 1 − XHe. 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 Δx and XHe of mixed system indicates different demixing definitions for systems with different compositions.

The inset of Fig. S13 shows the distribution of Δx of 27,000 atoms. Because Δx 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 Δx will lie in. Now, with a Δx calculated from a single frame or a Δx from a trajectory containing many frames, we can compare it with μ + 3σ value of Δx. If the value falls outside of the 3σ intervals, indicating that the probability of the system being mixed is only 0.135%. For a Δx from a single frame, it risks with these small probabilities, but for a Δx from a long trajectory of a wide time scale containing many frames, this risk can be further reduced. Thus, the μ + 3σ values of Δx for every XHe 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 Δ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 XHe < 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 XHe = 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 XHe, 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.

ΔG calculation of large systems
When calculating the ΔG, it is crucial to ensure that the term G(XHe) in Eq. 2 refers to a fully mixed system, as it cannot represent the Gibbs free energy of a specific helium abundance XHe 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: β(G(B) − G(A)) = −ln < e−β(U(B) − U(A)) > A,r (6)

Here β = 1/kBT is the inverse temperature, and UA, UB are the interaction energies of the initial state A and final state B, and the corresponding <...>A,r is the ensemble average over the initial state A. Since the relative positions of all atoms remain 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 ΔG in large systems can be roughly divided into two steps: 1) calculating the free energy of a completely mixed state A, and 2) calculating the free energy difference of volume change using the Jarzynski equality. In this work, 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 Δ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 Δ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: −(GB − GA) = 1/β ln < e−β(UB − UA) > A,r (7)

In this way, we have performed 10−7 VASP instructions for obtaining the <UB − UA>A,r, and we performed convergency test on sample size to ensure the accuracy. See Fig. S14.

ΔG calculation by thermodynamic integration
Thermodynamic integration calculates differences in free energies to get the target state (described by potential U1) free energy integrated from a reference state (described by potential U0). Then define a transition potential with two switching functions f(λ), g(λ):
U = f(λ)U0 + g(λ)U1, (8)

f(λ) and g(λ) satisfy f(0) = 1, f(1) = 0, g(0) = 0 and g(1) = 1. With Eq. 8, the free energy F is a function of N, V, T, λ:
F(λ) = −kBT ln Q(N, V, T, λ) (9)

The difference between the target state and reference state is:
F1 − F0 = ∫01 dλ ∂F(λ)/∂λ (10)

And it can be derived that ∂F(λ)/∂λ = <∂U/∂λ>λ through:
∂F/∂λ = −kBT ∂lnQ/∂λ = −kBT ∂lnZ/∂λ = −kBT ∂/∂λ ln ∫ dNr e−βU(r1,...,rN,λ) = −kBT ∫ dNr (−β ∂U/∂λ) e−βU(r1,...,rN,λ) / ∫ dNr e−βU(r1,...,rN,λ) = <∂U/∂λ>λ (11)

Replace the corresponding term in Eq. 10:
F1 − F0 = ∫01 dλ <∂U/∂λ>λ (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 coupling constants 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 Δ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 XHe, the delta Gibbs free energy can be calculated by Eq. 2, where G(XHe) = F(XHe) + PV. Then, by applying a common tangent construction to the ΔG curve, we can determine x1 and x2. According to the curvature of ΔG, the unstable and metastable regions are determined and disappear with the increase of temperature.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Overall, the topic of the paper is timely and relevant. As discussed by the authors,\nunderstanding the phase diagram of H-He mixtures is important for modelling the evolution\nand internal structure of Jupiter and Saturn. The work by Chang et al. presents another step\nforward in the investigation of H-He separation in planetary conditions.\nWhile the work is interesting, details regarding the calculations are clearly missing. In\naddition, the results should not be taken as “direct evidence” as claimed by the authors.\nThe results are obtained thanks to a very large-scale molecular dynamics (MD) simulation\nenabled by a machine learning potential (MLP) based on deep neural networks and trained\non density functional theory (DFT) data. The MLP is generated using a software package,\nand it seems that its functions are somehow “a black box”. The MLP's accuracy is validated\nagainst smaller-scale full DFT simulations under various conditions such as pressure,\ntemperature, and H/He fraction. The majority of the main text is dedicated to post-\nprocessing the MLP trajectory to identify demixing phenomena. However, as discussed\nbelow, the weak spots of the methods come before that.\nFinally, the connection to planetary science would require modelling of the planetary\ninterior/evolution using planetary models.\nThe detailed comments are given below:\n— Calculations:\n- While the results appear reasonable, unfortunately, the data presented in the paper are\ninsufficient to claim \"direct evidence.\" This is due to several reasons. First of all, it is\nimportant to note that the energy difference between competing phases is on the order of 1\nmilliHartree per atom (27 milli-electronvolts per atom). This energy difference is typically\nsmaller than the variations resulting from different exchange correlation functionals (such\nas PBE, vdW, SCAN) or higher levels of theory like QMC. The qualitative outcomes are\nacceptable, but not the quantitative phase boundary, which is crucial for planetary\nmodeling, as the phase boundary in pressure and temperature closely parallels the planet's\nadiabatic profile.\n- The authors justify their choice of the SCAN+rVV10 functional by referencing a relatively\nniche paper [Ref 22] where benchmarks against quantum chemistry data were feasible.\nUnfortunately, these benchmarks (an amino acid and an RNA base) do not adequately\nrepresent the behavior of the functional at gigapascal pressures and metallic materials.\nWhile the choice may seem reasonable, other options like PBE or HSE could have been\nequally valid. The quality of the MLP is directly influenced by the reference theory chosen\nfor fitting.\n- Another major concern revolves around the MLP itself. The authors present validation\nplots in the Supplementary note, which are not fully informative, and only demonstrate that\nthe MLP is not significantly wrong. The root mean square error (RMSE) is not explicitly\nstated. The heat map of the relative RMSE error of the forces indicates an error of about\n20%. In other words, the ML \"fit\" has an error of approximately 20%. This margin of error\naligns with a recent MLP study on dense hydrogen (Cheng et al., Nature 585 (7824), 217-\n220, 2020 not cited in this paper). The fact that an MLP with a 20% error in forces can still\nyield a sensible equation of state and radial function distribution is noteworthy and likely\nrelies on the self-cancellation of errors. This is not fully understood. However, while the\noverall equation of states remain consistent across a broad pressure range (as shown in Fig\nS3, although the image resolution is poor), identifying phase transitions is far more\nchallenging. Cheng et al. observed unexpected results compared to standard DFT, and this\nremains a topic of debate. The situation in this study is even more complex due to the\nmixture. To efficiently capture the nucleation process, the MLP must be exposed to rare\nevent transition configurations, such as phase separation, in the training set, where the MLP\nneeds to learn the energetics of an interface. The satisfactory behavior demonstrated in Figs\nS4-S5 indicates that the MLP can describe a fully mixed liquid of a small to medium-sized\nMD, where helium is dispersed in hydrogen.\nOverall, there is still a big accumulation of errors in the whole procedure (DFT training,\nsampling of structures and MLP fit) that prevents a conclusive answer.\n- Connection to Jupiter and Saturn\n- The authors present a new phase diagram for H-He demising and argue that their\ncalculations provide “direct evidence” for helium rain in Jupiter and Saturn. However, in\noder to apply the results to Jupiter and Saturn proper modelling of the interior and the\nthermal evolution are required\n(see Mankovitz & Fortney, 2020, Mankovitz et al. 2016 and references therein).\nIn this work all the authors do is to simply present the adiabatic of Jupiter and Saturn and\nclaim a “direct evidence” for helium rain in the planets. This is clearly oversimplified and in\norder to make such a statement planetary simulations are required. I acknowledge that the\nauthors are not planetary scientists and focus on the phase diagram, but in that case they\nneed to be fair and clearly state that their result must be followed by planetary modelling.\nIn addition, it is now clear that adiabatic models are inappropriate to describe the interiors\nof Jupiter and Saturn as both planets are expected to have non-adiabatic regions. The\nminimum the authors could do is to plot different published P-T profiles and see how it\naffects their predictions. In any case, it is clear that the authors can’t claim that they present\na “direct evidence of helium rain”.", "reviewer_2": "n5. the data underlying Figure 1 should be made available so that they can be used for\napplication to Jupiter and Saturn modeling.\n6. Which functional was used to produce the black dots?\n7. Figure 2 should include the two experimental points for xHe=0.11.\n8. I would like to see the curves", "reviewer_with_the_reweighted_conditional_propability": "From Figure 1 I understand\nP_cond, what I do not understand why reweighting is needed and nor the meaning of the\nantisymmetric grey shadows. More explanation is needed, here or in the methods section.\nThe caption of Extended Data Figure 1 suggests that reweighting is needed to infer the\nresults that would have been obtained with 27.000 particles from simulations done with\n9600 particles. Is that so? please explain in the main text.\n12. Figure S1: Inlay shows a \"Distribution of Errors\". This label suggests that the true value is\nwhere DLP agrees with the DFT-based predictions. If this is intended it should be\nmade clear. Also, why do potential and force have a zero-transition, eg what does that\nphysically mean? In particular for the force the relative and the absolute deviations are\nlargest around\nf ~ 0, so this point deserves an explanation.\n13. Figure S2: It should be mentioned in the caption that the panels are for different He-\nabundances ranging from pure hydrogen (top left) to pure helium (bottom right). The color\ncode legend is missing. Values near the contour lines suggest that in H-rich systems, the\nforce mismatch is up 30% at low T and 1 Mbar. Is that correct?\n14. Figure S3: What is ''MD data'' ? Maybe used the common name DFT-MD consistently\n15. Figure S4, panel P over yHe: what do the small and the large filled circles stand for?\nThe large filled circles are beyond the\nxrange and yrange, or is this a 3D projection? Please explain in the figure caption."}, "author_response": "Response to Reviewers:\n----------------------------------------------------\nReply to Reviewer #1:\n----------------------------------------------------\nOverall, the topic of the paper is timely and relevant. As discussed by the authors,\nunderstanding the phase diagram of H-He mixtures is important for modelling the\nevolution and internal structure of Jupiter and Saturn. The work by Chang et al. presents\nanother step forward in the investigation of H-He separation in planetary conditions.\nWhile the work is interesting, details regarding the calculations are clearly missing.\nIn addition, the results should not be taken as “direct evidence” as claimed by the authors.\nThe results are obtained thanks to a very large-scale molecular dynamics (MD)\nsimulation enabled by a machine learning potential (MLP) based on deep neural\nnetworks and trained on density functional theory (DFT) data. The MLP is generated\nusing a software package, and it seems that its functions are somehow “a black box”.\nThe MLP's accuracy is validated against smaller-scale full DFT simulations under\nvarious conditions such as pressure, temperature, and H/He fraction. The majority of\nthe main text is dedicated to post-processing the MLP trajectory to identify demixing\nphenomena. However, as discussed below, the weak spots of the methods come before\nthat. Finally, the connection to planetary science would require modelling of the\nplanetary interior/evolution using planetary models.\nResponse:\nWe thank the reviewer for their insightful comments on our paper. We are encouraged\nby recognizing our work as another step forward to investigate H-He separation in\nplanetary conditions. Our work aims to achieve higher resolution intervals for H-He\ndemixing through large-scale and high-precision molecular dynamics (MD)\nsimulations. These simulations are enabled by employing a machine learning potential\n(MLP) trained on density functional theory (DFT) data. Initially, we used the term\n‘direct evidence’ to suggest that we could observe the degree of demixing in the H-He\nsystem directly through large-scale MD simulations and structural analysis. This\napproach offers a more straightforward visualization compared to calculating the free\nenergy by adjusting the H-He ratio in a small system. However, following discussions\nwith all authors, we acknowledged that the term \"direct evidence\" was inaccurate and\nmisleading, so we replaced it in the title with “Theoretical evidence of H-He demixing\nunder Jupiter and Saturn conditions”.\nWe appreciate the reviewers’ concern regarding the perception of machine learning\nmodels as a ‘black box’ to some extent and the ongoing challenge of achieving fully\nexplainable models. Nevertheless, the machine learning potential we constructed is\ndesigned to ensure mathematical smoothness and physical symmetry and adhere to\nconservation laws (energy, momentum, and angular momentum). This design strategy\nis to avoid the generation of non-physical results. As long as the fitting and\ngeneralization accuracy is high enough, the MD simulation with MLP gives physically\nmeaning results. Other potential construction methods, such as those with analytic\nforms, have limited expressive capability and often cannot accurately describe\nbehaviors across diverse temperature and pressure ranges. A typical example is the\nphase diagram prediction of water. The TIP4P/2005 model, renowned for its accuracy\nin modeling water, performs well in predicting the phase diagram at low and\nintermediate pressures. However, at higher pressures, significant deviations from\nexperimental observations are observed, particularly in the boundary lines separating\nthe ice VIII, VII, and VI phases. Furthermore, the TIP4P/2005 model predicts a first-\norder transition from ice VII to a plastic phase, which has not been experimentally\nconfirmed to date [Phys. Rev. Lett. 126, 236001 (2021)]. Numerical potentials are\nsignificantly less effective to neural networks in addressing the problems posed by the\n‘curse of dimensionality’. It will result in larger errors in the description of phonon\ndispersion, which is a property derived from higher-order derivatives of the potential\nenergy surface (PES). The modified embedded-atom-method (MEAM), a numerical\npotential commonly employed for modeling metal systems, exhibits good agreement\nwith DFT results in the low-frequency range of phonon dispersion, but significant\ndiscrepancies arise in the high-frequency range [J. Chem. Phys. 159, 204702 (2023)].\nThe superior fitting ability of machine learning for such problems has been rigorously\ndemonstrated in mathematical studies on this topic [IEEE Transactions on Information\nTheory, 39, 930-945 (1993); Communications in Mathematical Sciences, 17, 1407-\n1425 (2019); Constr. Approx. 55, 369-406 (2022)].\nThe MLP is constructed through rigorous physical principles and considerations,\nensuring its adherence to the fundamental laws and symmetries. Firstly, a local\ncoordinate frame is established for every atom and its neighboring atoms within a cutoff\nradius, thereby preserving the translational, rotational, and permutational\nsymmetries of the system [Phys. Rev. Lett. 120, 143001 (2018)]. Secondly, the\npotential energy associated with each atomic configuration is the sum of all \"atomic\nenergies\" within a cutoff radius. The environmental dependence, representing the\nmany-body character of interactions, is complex and nonlinear. The neural network can\ncapture the analytical dependency on the coordinates of atoms within the environment\nconcerning the composition of the sequence of mappings associated with the individual\nhidden layers, which collectively model the energy contributions of atoms within the\nspecified cutoff radius. The additive nature of potential energy naturally preserves the\nextensive character of the potential energy. Hence, if all local configurations are\nadequately represented in the training data, the model's accuracy is ensured for an\nexpansive system. The attainment of uniform accuracy across all local environments is\nfacilitated by a thorough exploration within the relevant thermodynamic regime,\nemploying the DP-GEN methodology. This approach collects the most crucial\nconfigurations for enhancing the model's accuracy, thus enables it to achieve a user-\nspecified level of accuracy with the most concise training dataset. Further details\nregarding the DP-GEN scheme can be found in [Computer Physics Communications\n253, 107206 (2020)].\nWe have conducted intensive DFT and MLP energy and force verifications on 28-atom\nsystems under various temperatures, pressures, and H/He fractions. Structure and\ndynamics (RDF & MSD) are validated using 216-atom and 1728-atom systems, albeit\nless densely sampled. These findings are detailed in Supplementary Information Note\n1. Given the inherent extensive character of MLP, which ensures consistent accuracy\nacross larger systems, we believe that the comprehensive testing conducted on smaller\nsystems at different scales demonstrates the reliability of our MLP. Furthermore,\nrecognizing that our primary focus revolves around utilizing the MLP trajectory to\ncalculate and identify demixing phenomena, we acknowledge the critical importance\nof ensuring the accuracy of our structural analysis. Thus, we conducted additional\nvalidations using a 2,048-atom system to benchmark the accuracy of our\nstructural analysis, which is essential for ensuring the integrity and reliability of\nour results.\nWe are grateful for the reviewers’ insights that prompted us to revisit the planetary\nmodel. Initially, our approach for identifying the formation of helium rain relied on\nobserving whether the planetary adiabat crossed a separation boundary in the phase\ndiagram. However, this approach seemed to provide only a rough approximation. We\nrecognize that the process of helium rain formation is more complex. As the hydrogen-\nhelium mixture cools within the planetary interior, it reaches temperatures and pressures\nwhere helium becomes immiscible and separates from the hydrogen phase. This leads\nto helium rainout and the formation of a helium-rich interior surrounded by a helium-\ndepleted outer envelope. In the region with a stabilizing helium gradient, double-\ndiffusive convection can occur, leading to a super-adiabatic temperature gradient and\nless efficient heat transport compared to adiabatic convection. This double-diffusive\nregion may consist of layered semi-convective zones with varying degrees of super-\nadiabaticity rather than a single adiabatic profile. Therefore, while the outer molecular\nenvelope may still be adequately approximated as adiabatic, the deep interior is\nexpected to deviate from an adiabat due to helium rain and associated double-diffusive\nconvection/layering effects. Modeling these non-adiabatic regions is crucial for\naccurately describing the H-He demixing in Jupiter and Saturn. Incorporating these\ncomplexities into the model may result in temperature adjustments in the anticipated\nhydrogen-helium separation region. However, a thorough understanding of the\nplanetary model requires specialized knowledge in planetary physics and advanced\ntechniques such as the MESA stellar evolution code, Bayesian statistics, and Markov\nChain Monte Carlo sampling. This topic appears to exceed the scope of this article. Our\nobjective is to provide foundational phase diagram data to inform future modeling\nefforts. An accurate phase diagram is needed to predict at what stage of cooling this\nseparation begins in the planets' interiors. Our calculations offer detailed insights,\nincluding the extent of H-He demixing, which could offer new perspectives for\ncalibrating the extent of helium depletion or the location of the stabilizing gradient. Our\nprimary focus is on providing data support for planetary modeling rather than engaging\ndirectly in the modeling process. We will undertake a comprehensive adjustment of our\nmanuscript to ensure that our findings are appropriately positioned.\nMajor comments:\nComment 1:\nWhile the results appear reasonable, unfortunately, the data presented in the paper are\ninsufficient to claim \"direct evidence.\" This is due to several reasons. First of all, it is\nimportant to note that the energy difference between competing phases is on the order\nof 1 milliHartree per atom (27 milli-electronvolts per atom). This energy difference is\ntypically smaller than the variations resulting from different exchange correlation\nfunctionals (such as PBE, vdW, SCAN) or higher levels of theory like QMC. The\nqualitative outcomes are acceptable, but not the quantitative phase boundary, which is\ncrucial for planetary modeling, as the phase boundary in pressure and temperature\nclosely parallels the planet's adiabatic profile.\nResponse:\nWe agree with the reviewer that the energy difference between different exchange-\ncorrelation functionals is critical in this particular interval; thus, after careful\nexamination, we choose to train a DLP based on SCAN+rVV10 functional (this will be\nfurther explained in Comment 2). We believe it is important to emphasize that our study\nrepresents the current state-of-the-art in terms of accuracy, particularly within the\ncontext of large-scale systems and considering size effects through long-timescale (100\nps) and large-system (27,000) molecular dynamics.\nTo further address the reasonable concern raised by the reviewer, we also retrained the\npotential function based on vdW-DF and constructed a phase separation boundary for\nprotosolar helium abundance, as illustrated in the figure below.It can be observed that\nthe results of the structural analysis of large systems using MLP based on vdW-\nDF reveal a larger phase separation region compared to DFT employing the same\nfunctional. Moreover, both the vdW-DF-based MLP and the SCAN+rVV10-based\nMLP exhibit similar phase separation landscapes that extend substantially beyond\nplanetary adiabatic lines. This consistency reflects the credibility of our\ncomputational findings.\nAdded Figure:\nExtended Data Fig. 4. Miscibility diagram at helium abundance of\nof vdW-DF-based DLP. The solid black line with pink squares represent the\n𝑿𝑯𝒆 = 𝟎.𝟎𝟖𝟗\nhydrogen-helium immiscibility lines of vdW-DF-based DLP. The color of the colored\narea represents the value of , and 0.0181 on the color bar is the minimum value\nthat can be considered as the sign of immiscibility for helium abundance of protosolar.\n∆𝑥\nOthers are the same as Fig. 2 in the main text.\nHowever, selecting different functionals can indeed impact the temperature of phase\nseparation, and we have not attempted higher levels of theory such as QMC. Therefore,\nwe have decided to revise the statement in the manuscript to a more focused discussion\nof the results rather than diverging into other aspects.\nThe modification of the picture:\n1) Fig. 3\nAdd X = 0.089 demixing points of vdW-DF-based DLP in Fig .3a;\nHe\nAdd pressure axis in Fig .3c.\n2) Extended Data Fig. 4\nAdd Extended Data Fig. 4 displayed above.\nThe modification of the article:\n1) Modified the title of the paper to be “Theoretical evidence of H-He demixing under\nJupiter and Saturn conditions”\n2) Modified lines 26-27\n“Our results revise the structures of Jupiter and Saturn where H-He demixing takes\nplace in a large fraction of the interior radii.” to\n“Our results suggest a possibility that H-He demixing takes place in a large fraction\nof the interior radii of Jupiter and Saturn”.\n3) Modified lines 28-29\n“This direct evidence of an H-He immiscible layer supports the formation of helium\nrain and explains the helium reduction in atmosphere of Jupiter and Saturn” to\n“This indication of an H-He immiscible layer hints at the formation of helium rain\nand offers a potential explanation for the decrease of helium in the atmospheres of\nJupiter and Saturn.”\n4) Added “To assess whether different functionals fundamentally influence the\nconclusions, we also trained a DLP based on vdW-DF and conducted identical\nsimulations under protosolar helium abundance, as indicated by the pink squares in\nFig. 3a. It can be observed that changing the functional still yields similar\nconclusions. The miscibility diagram calculated at protosolar helium abundance is\nprovided in Extended Data Fig. 4.” at line 210.\n5) Modified lines 218-219\n“Combining the calculated miscibility diagram with a planet model3, we derived\nthe internal structures of Jupiter and Saturn (Fig. 3c)” to\n“Combining the calculated miscibility diagram with a planet model3, and assume\nthe interior of the planet is adiabatic, we derived the nowadays internal structures\nof Jupiter and Saturn (Fig. 3c).”\n6) Modified line 227\n“In summary, we give direct evidence for the H-He immiscibility from first-\nprinciples” to\n“In summary, we present direct observations suggesting the H-He immiscibility\nfrom first-principles.”\nComment 2:\nThe authors justify their choice of the SCAN+rVV10 functional by referencing a\nrelatively niche paper [Ref 22] where benchmarks against quantum chemistry data were\nfeasible. Unfortunately, these benchmarks (an amino acid and an RNA base) do not\nadequately represent the behavior of the functional at gigapascal pressures and metallic\nmaterials. While the choice may seem reasonable, other options like PBE or HSE could\nhave been equally valid. The quality of the MLP is directly influenced by the reference\ntheory chosen for fitting.\nResponse:\nWe thank the reviewer for pointing out the inadequacies in our explanation regarding\nthe selection of the potential function, particularly in terms of the functional’s behavior\nunder gigapascal pressures and with metallic materials.\nRegarding gigapascal pressures, SCAN+rVV10 has been shown to have improved\nagreement with experimental data compared to PBE and SCAN on liquid-liquid phase\ntransition in nitrogen-oxygen mixtures, attributed to its revised noncovalent interactions\nspanning short- to long-range van der Waals interactions [Yang-Shun Lan, Phys. Rev.\nB 103, 144105 (2021)].\nSimilar improvements have been noted in the following works [Tantardini, C., J. Phys.\nChem. C 126, 11358–11364 (2022), Anh, L. T., Sci. Rep. 9, 8731 (2019)]. Regarding\nmetallic materials, SCAN+rVV10 has been demonstrated to align well with\nexperimental results [Patra A., Proc. Natl. Acad. Sci. 114, E9188–E9196 (2017)],\noutperforming PBE, PBEsol, and SCAN. Notably, the addition of rVV10 to SCAN\nyields a highly accurate method for systems with diverse bonds.\nWe understand the reviewer’s point about considering alternatives such as PBE or HSE,\nrecognizing that different functionals may provide valid options. Therefore, in addition\nto the previously mentioned effective practice of SCAN+rVV10, our selection also\nconsiders the specific characteristics of the system under investigation. In a\nbenchmarking study of density functionals for hydrogen-helium mixtures based on\nQMC [Phys. Rev. B 93, 035121 (2016)], TPSS, BLYP, and vdW-DF were identified as\nthe most accurate functionals by most metrics. Although this study predates the\nmaturity of the SCAN+rVV10 method, we can infer that metaGGA (such as TPSS and\nSCAN) should be a better option. The liquid-liquid phase transition (LLPT) is believed\nto play an important role in H-He demixing. Geng [Geng, H. Y. Phys. Rev. B 100,\n134109 (2019)] notes that dense hydrogen near dissociation is best represented between\nPBE and vdW-DF, with PBE underestimating the dissociation pressure and vdW-DF\noverstabilizing the molecule.\nHowever, we conducted comparative analyses and found that the LLPT line of\nSCAN+rVV10 lies right in the middle and is consistent with QMC. In conclusion,\nwe believe that the SCAN+rVV10 functional offers unique advantages and is well-\nsuited for the current research objectives.\nWe have incorporated all mentioned references into our revised manuscript to provide\na more comprehensive rationale for our potential function selection.\nThe modification of the article:\n1) Removed lines 76-80 “This scheme is chosen for its ‘best-of-both-worlds’ feature\n– SCAN captures the short and intermediate-range interactions and rVV10\ndescribes the long-range vdW interaction21.The combination of SCAN+rVV10 has\nbeen shown to have accuracy better than 1 kcal/mol for benchmarking the\nMP2+∆CCSD(T) results, and be able to reproduce subtle features of the potential\nenergy surface22.”\n2) Add a new section in Methods “DLP potential training”.\n“DLP is constructed based on SCAN+rVV10. This scheme is chosen for its ‘best-\nof-both-worlds’ feature – SCAN captures the short and intermediate-range\ninteractions and rVV10 describes the long-range vdW interaction37. The\ncombination of SCAN+rVV10 has been shown to have accuracy better than 1\nkcal/mol for benchmarking the MP2+∆CCSD(T) results, and be able to reproduce\nsubtle features of the potential energy surface38. In addressing gigapascal pressures,\nSCAN+rVV10 has been shown to improve the agreement with experimental data\ncompared with PBE and SCAN on liquid-liquid phase transition in nitrogen-oxygen\nmixtures, thanks to the revised noncovalent interactions of short- to long-range van\nder Waals interaction39. Similar better match appears in the works of Tantardini40\nand Anh41. Regarding metallic materials, SCAN+rVV10 has been demonstrated to\nalign well with experimental results42, outperforming PBE, PBEsol, and SCAN.\nNotably, the addition of rVV10 to SCAN yields a highly accurate method for\ndiversely bonded systems.\nWe use Deep Potential Generator (DP-GEN)43 to cover temperatures from 1,000 K\nto 13,000 K and pressures from 0.5 Mbar to 10 Mbar. 28 helium proportions evenly\ndistributed from 0 to 1 were used as training set. A total of 106,817 configurations\nin the training set, and trained for 12,000,000 steps.”\n3) Removed Supplementary Information lines 3-5 “The deep-learning potential (DLP)\nis constructed based on SCAN+rVV10. We use Deep Potential Generator (DP-\nGEN)1 to cover temperatures from 1,000 K to 13,000 K and pressures from 0.5\nMbar to 10 Mbar. 28 helium proportions evenly distributed from 0 to 1 were used\nas training set and another”.\n4) Modified Supplementary Information lines 4-5\n“15 helium proportions evenly distributed from 0 to 1 were used as testing set” to\n“In addition to the 28 helium proportions used during training, we evenly selected\nan additional 15 helium proportions evenly distributed from 0 to 1 were used as\ntesting set.”\n5) Added references\n37. Peng, H., Yang, Z.-H., Perdew, J. P. & Sun, J. Versatile van der Waals Density\nFunctional Based on a Meta-Generalized Gradient Approximation. Phys. Rev. X 6,\n041005 (2016).\n38. Chen, J., Xie, W., Li, K., Zhang, S. & Sun, Y.-Y. Benchmarking PBE+D3 and\nSCAN+rVV10 methods using potential energy surfaces generated with\nMP2+ΔCCSD(T) calculation. Chin. Phys. B 29, 013102 (2020).\n39. Lan, Y.-S. et al. First-order liquid-liquid phase transition in nitrogen-oxygen\nmixtures. Phys. Rev. B 103, 144105 (2021).\n40. Tantardini, C., Jalolov, F. N. & Kvashnin, A. G. Crystal Structure Evolution of\nFluorine under High Pressure. J. Phys. Chem. C 126, 11358–11364 (2022).\n41. Anh, L. T., Wada, M., Fukui, H., Kawatsu, T. & Iitaka, T. First-principles\ncalculations of the epsilon phase of solid oxygen. Sci. Rep. 9, 8731 (2019).\n42. Patra, A., Bates, J. E., Sun, J. & Perdew, J. P. Properties of real metallic\nsurfaces: Effects of density functional semilocality and van der Waals nonlocality.\nProc. Natl. Acad. Sci. 114, E9188–E9196 (2017).\n43. Zhang, Y. et al. DP-GEN: A concurrent learning platform for the generation of\nreliable deep learning based potential energy models. Comput. Phys. Commun. 253,\n107206 (2020).\nComment 3:\nAnother major concern revolves around the MLP itself. The authors present validation\nplots in the Supplementary note, which are not fully informative, and only demonstrate\nthat the MLP is not significantly wrong. The root mean square error (RMSE) is not\nexplicitly stated. The heat map of the relative RMSE error of the forces indicates an\nerror of about 20%. In other words, the ML \"fit\" has an error of approximately 20%.\nThis margin of error aligns with a recent MLP study on dense hydrogen (Cheng et al.,\nNature 585 (7824), 217-220, 2020 not cited in this paper). The fact that an MLP with a\n20% error in forces can still yield a sensible equation of state and radial function\ndistribution is noteworthy and likely relies on the self-cancellation of errors. This is not\nfully understood. However, while the overall equation of statesremain consistent across\na broad pressure range (as shown in Fig S3, although the image resolution is poor),\nidentifying phase transitions is far more challenging. Cheng et al. observed unexpected\nresults compared to standard DFT, and this remains a topic of debate. The situation in\nthis study is even more complex due to the mixture. To efficiently capture the nucleation\nprocess, the MLP must be exposed to rare event transition configurations, such as phase\nseparation, in the training set, where the MLP needs to learn the energetics of an\ninterface. The satisfactory behavior demonstrated in Figs S4-S5 indicates that the MLP\ncan describe a fully mixed liquid of a small to medium-sized MD, where helium is\ndispersed in hydrogen. Overall, there is still a big accumulation of errors in the whole\nprocedure (DFT training, sampling of structures and MLP fit) that prevents a conclusive\nanswer.\nResponse:\nThe reviewer’s concerns about the MLP’s performance in the study, particularly the\nvalidation plots provided in the Supplementary note and the challenges in fully\nassessing the MLP’s accuracy, are duly noted.\nThere are two critical factors to consider in the accuracy of the MLP: energy and force.\nWe consider the precision of energy fitting to be particularly sensitive when it comes to\nthermodynamic boundaries separating phases. Therefore, in the training of potential\nfunctions, we specifically focus on the accuracy of energy fitting. The graphs\npresented below illustrate the distribution of RMSE of energies and RMSE over\nstandard deviation (std) of energies. It can be observed that across a broad range of\ntemperature and pressure, the fitting accuracy of energies is consistently high.\nAdded Figure:\nFig. S4. The energy RMSE. The top right corner of each subplot displays their\nrespective helium abundance. The red points represent thermodynamic conditions\nsampled as testing sets.\nFig. S5. The ratio of energy RMSE between DLP and DFT over energy std of DFT.\nThe top right corner of each subplot displays their respective helium abundance. The\nred points represent thermodynamic conditions sampled as testing sets.\nAs for the fitting accuracy of force, it is noteworthy that Fig. S2 represents the\n“RMSE of forces between MLP and DFT over the standard deviation (std) of DFT\nforces”, rather than the “relative RMSE”. This measure, also dimensionless, focuses\non the comparison with fluctuations of real values. For an MLP, a force RMSE below\nthe force std (RMSE/std < 1) can be considered as having some fitting capability.\nAchieving 10% to 20% indicates a high level of precision. From a statistical perspective,\nthe MLP we developed has achieved an R-squared value (1- (RMSE/std)^2) of over 96%\nfor force fitting, demonstrating accurate model fitting. Physically, the current fitting\nerror of MLP can be approximated as increasing the system's force fluctuations. We\nhave statistically analyzed the force distribution within a single frame, as shown in the\ngraph below.\nIt can be observed that forces tend to follow a Gaussian distribution within the statistical\naverage of long-term large systems. This implies that for certain property predictions,\nconducting a sufficiently long statistical average leads to the cancellation of force\nfluctuation effects, as noted by the reviewer. Therefore, we conducted a series of tests\non RDF, EOS, and diffusion coefficients and found that these properties exhibited good\nconsistency with DFT (Fig. S7, Fig.S8, Fig.S9). This demonstrates the robustness of\nour approach to handle reasonable levels of noise or errors in force calculations. Our\nprevious work has confirmed a similar perspective. The errors in force calculations have\na large convergence tolerance and can be considered noise in the Langevin dynamics\nfor constant-temperature molecular dynamics simulations. Even with introducing\nGaussian noise, the results obtained from quantum Langevin MD remain consistent\nwith those obtained from highly accurate methods, and the dynamical properties of the\nsystems are reproduced accurately [Europhys. Lett. 88, 20001 (2009); Phys. Rev. Lett.\n109, 175701 (2012)].\nThis idea of approximating accurate results due to fluctuation self-cancellation is also\nutilized in DFT calculations using the CPMD method. Despite some offset between\nforces and potential energy surfaces, the overall behavior oscillates around the true\nvalues, indicating reliability and applicability in predicting various properties such as\nelectronic structures. In addressing the hydrogen-helium phase separation issue under\nplanetary conditions, we relaxed every system for a minimum of 100 ps and conducted\nlong-term statistical averaging. Therefore, when the fitting error of the potential\nfunction is covered by the thermodynamic fluctuations of the system and remains\nwithin a relatively low range, we can average out this error through the self-cancellation\neffect of fluctuations.\nWe acknowledge the reviewer’s concerns regarding the accuracy of MLP functions in\nphase transition processes, especially given the reviewer’s reference to Cheng et al.’s\nprevious work on the hydrogen liquid-liquid phase transition, where unexpected results\nwere encountered despite achieving a certain level of fitting accuracy of the potential\nfunction compared to DFT. Following Karasiev's commentary on Cheng et al.’s\nwork in response to their publication, as outlined in the literature [Karasiev, V. V.\nNature 600, E12–E14 (2021)], we directly validated the results of phase separation\nusing DFT calculations of a large system containing 2,048 atoms.\nAdded text:\nWe uniformly selected four temperature-density points within the system described in\nFig. 3 of the main text, all of which are situated within the intervals of relatively higher\nrelative errors as shown in Fig. S2. We conducted DFT-MD simulations for a\nconsiderably long time scale (2 ps) with a system comprising 2,048 atoms (Certainly,\nthe outcomes derived from such systems are only suitable for comparative analysis. For\na genuine structural analysis, 2,048 atoms remain insufficient). For each frame of the\ntrajectory, structural analysis, as described in the main text, was performed to compare\nwith the results obtained from the MLP model, thereby assessing whether the accuracy\nof the potential function meets the requirements for addressing the core issues discussed\nin this study. As depicted in the figure below, the comparison reveals that, in\nconjunction with the validation results of the potential function's energy, forces,\nstructures, and dynamics obtained previously, the MLP results are deemed reliable. It\nshould be noted that the energy and force RMSE of DLP and DFT for the following\nfour cases are listed in the table below, which are consistent with those of the smaller\nsystems, thus the model is valid in the MD simulations of 27K atoms.\nP [Mbar] T [K] X RMSE_energy[eV/atom] RMSE_force[eV/angstrom]\nHe\n1 3000 0.073 0.005 0.300\n1.5 4000 0.073 0.010 0.303\n7 4000 0.073 0.010 0.297\n2 6000 0.2 0.011 0.337\nWe fully agree with the reviewer's perspective regarding the importance of\nincorporating rare event transition configurations, especially those related to phase\nseparation, into the training dataset. Previous machine learning studies on the liquid-\nliquid phase transition of phosphorus have highlighted the importance of collecting\nlow-probability yet highly significant transition state arrangements [M. Yang, et al.,\nPhys. Rev. Lett. 127, 080603 (2021)]. Consequently, we employed the existing MLP\npotential function to generate 1000 phase-separated structures, each consisting of 216\natoms, uniformly across 10 pressure points ranging from 1 to 10 megabars. These\nstructures underwent first-principles single-point calculations and were incorporated\ninto the training dataset. The resulting potential functions obtained in this manner were\nalso directly compared with structural analysis (the “DLP(with interface)-MD” line).\nFrom the comparison results, it seems that high-precision potential functions can be\nobtained regardless of whether phase-separated structures are included in the process\nof phase transition. We believe that a possible reason for this observation is that for\nliquid-liquid phase transitions (like phosphorus or hydrogen), rare phase separation\noccurs at the potential barrier on the free energy surface. Collecting the structures at the\npotential barrier requires other methods, such as enhanced sampling via recursively\nreconstructing the free energy surface. However, for the H-He phase separation, phase\nseparation occurs at the potential well on the free energy surface, and most phase\nboundaries are not pure hydrogen and helium, but rather a physical mixture of varying\ndegrees of hydrogen-rich and helium-rich regions. Therefore, they can be equivalently\ntreated as various mixtures of components. Thus, if the sampling of each component is\nsufficiently performed, reasonable potential functions can be obtained.\nWe thank the reviewer for pointing out the missing aspect of force RMSE and the poor\nresolution of Fig. S3. We ensure that these concerns are fully addressed in the revised\nmanuscript.\nThe modification of the picture:\n1) Fig. S1\nAdd RMSE value of energy and force;\nDisplay data density through color coding.\n2) Fig. S2\nAdd Fig. S2 displayed below.\nFig. S2. The force RMSE. The top right corner of each subplot displays their\nrespective helium abundance. The red points represent thermodynamic conditions\nsampled as testing sets.\n3) Fig. S4\nAdd Fig. S4 displayed above.\n4) Fig. S5\nAdd Fig. S5 displayed above.\n5) Fig. S6\nEnhanced the resolution;\nDisplay overlapping portions separately;\nModified title.\nFig. S6. Comparison of the equations of state from DLP-MD data and DFT-\nMD data. a, Energy as a function of pressure for three temperatures and five He\nfractions x. b, Volume as a function of pressure for the same conditions as a. Solid\nand dashed lines are fits.\n6) Fig. S9\nAdd Fig. S9 displayed below.\nFig. S9. Comparison of helium abundance differentiation achieved by\nreweighted conditional probability between DLP and DFT with a system size\nof 2,048 atoms.\nThe modification of the article:\n1) Added “We apply this method to the DFT-MD trajectories of 2,048 atoms with\nhelium abundances of X = 0.073 or 0.2. In addition to validating the accuracy of\nHe\nthe DLP, it also demonstrates that the method can accurately differentiate helium\nabundance on a case-by-case basis, even when the separation levels are very close\n(see Fig .S9 in Supplementary Information Note 1)” in line 92.\n2) Added paragraphs indicated in the preceding text in Supplementary Information\nline 34.\nComment 4:\nThe authors present a new phase diagram for H-He demising and argue that their\ncalculations provide “direct evidence” for helium rain in Jupiter and Saturn. However,\nin oder to apply the results to Jupiter and Saturn proper modelling of the interior and\nthe thermal evolution are required (see Mankovitz & Fortney, 2020, Mankovitz et al.\n2016 and references therein). In this work all the authors do is to simply present the\nadiabatic of Jupiter and Saturn and claim a “direct evidence” for helium rain in the\nplanets. This is clearly oversimplified and in order to make such a statement planetary\nsimulations are required. I acknowledge that the authors are not planetary scientists and\nfocus on the phase diagram, but in that case they need to be fair and clearly state that\ntheir result must be followed by planetary modelling. In addition, it is now clear that\nadiabatic models are inappropriate to describe the interiors of Jupiter and Saturn as both\nplanets are expected to have non-adiabatic regions. The minimum the authors could do\nis to plot different published P-T profiles and see how it affects their predictions. In any\ncase, it is clear that the authors can’t claim that they present a “direct evidence of helium\nrain”.\nResponse:\nWe thank the reviewer for recognizing our background as non-planetary scientists. We\nacknowledge that our expertise in planetary modeling is limited. Our primary objective\nis to present a theoretical phase diagram that enhances the understanding of planetary\ninternal material structures and evolutionary states. To ensure the reliability of our\nconclusions, we primarily present in the main text the findings pertaining to the\nobserved helium abundance in Jupiter and Saturn, while also noting the results for the\nprotosolar helium abundance for comparison with previous publications. This can be\nseen as an illustration of the conditions spanning from the early stages of planet\nformation to the present day.\nWe have also considered the reviewer's suggestion and have now included adiabats\ncorresponding to different surface temperatures within the planetary interiors. In our\nprevious miscibility diagram at protosolar helium abundance, only the current planetary\nadiabatic lines were shown. We believe that by incorporating adiabatic lines of planets\nwith varying surface temperatures, comparisons can be made more conveniently.\nDetermining the formation of helium rain solely based on whether the planetary adiabat\ncrosses the separation boundary in the phase diagram would be an oversimplification.\nAccurately assessing the formation of helium rain requires a combination of thermal\nevolution models, planetary physical models, and various factors such as strong\nsuperadiabaticity due to double-diffusive convection. This endeavor requires a range of\nparameters, deep knowledge of planetary physics, and skillful application of\ntechnologies like the MESA stellar evolution code or Bayesian statistics, which are\nbeyond the scope of our current work. Our primary focus is to provide parametric\nconditions for planetary modeling. For instance, the ab initio H-He phase diagram\nspecifies the pressure-temperature conditions under which helium becomes immiscible\nand separates from the hydrogen phase. An accurate phase diagram is crucial for\npredicting the stage of cooling at which this separation initiates in the planets’ interiors.\nThe shape of the immiscibility region in the phase diagram determines the equilibrium\nhelium abundance as a function of pressure and temperature, influencing the depletion\nof the outer molecular envelope as rainout proceeds. Since the helium gradient arises\nfrom the increasing helium concentration inwards following the phase diagram, the\nprofile of this gradient depends on the diagram's shape, impacting the double-diffusive\nconvection region. The energy released from helium differentiation scales with the\namount of helium that rains out, which is controlled by the phase diagram demarcating\nimmiscible versus miscible regions. Through our large-scale calculations and structural\nanalysis, we can also deduce the density ratios from the degree of demixing. This\nparameter can be used to parameterize the efficiency of heat transport by double-\ndiffusive convection in the region where helium is raining out. We also hope that our\npresented data can be integrated into future research on planetary physics to advance\nthe evolution model. Given that our study did not involve dynamic planetary modeling,\nwe have removed phrases like \"direct evidence of helium rain\" and adopted a more\nrestrained expression, as per the reviewers' constructive comments.\nThe modification of the picture:\n1) Extended Data Fig. 3\nAdd adiabats corresponding to different surface temperatures;\nAdd X = 0.089 demixing points of vdW-DF-based DLP.\nHe\nExtended Data Fig. 3.Miscibility diagram at helium abundance of\nof vdW-DF-based DLP. The solid black line with pink squares represent the\n𝑿𝑯𝒆 = 𝟎.𝟎𝟖𝟗\nhydrogen-helium immiscibility lines of vdW-DF-based DLP. The color of the colored\narea represents the value of ∆x, and 0.0181 on the color bar is the minimum value\nthat can be considered as the sign of immiscibility for helium abundance of protosolar.\nThe dash black lines represent the isentropes of various planetary surface temperature of\nprotosolar helium concentration, in which the isentrope corresponding to a planetary\nsurface temperature of 200K reflects a Saturn age of 0.746 Gyr. Others are the same as\nFig. 2 in the main text.\nThe modification of the article:\n1) Added “While the outer molecular envelope may still be well approximated as\nadiabatic, the deep interior is expected to depart from an adiabat due to helium rain\nand associated double-diffusive convection effects. Modeling these non-adiabatic\nregions is crucial for accurately describing the H-He demixing in Jupiter and Saturn.\nIncorporating these complexities into the model can lead to temperature shifts in\nthe predicted hydrogen-helium separation region35,36. This work offers more precise\nparametric conditions for intricate planetary modeling. The ab initio H-He phase\ndiagram delineates the P-T conditions at which helium becomes immiscible with\nhydrogen. It is crucial for determining the equilibrium helium abundance,\ninfluencing the depletion of the outer molecular envelope as rainout proceeds. The\nhelium gradient, driven by the increasing helium concentration inward influences\nthe profile of the double-diffusive convection region. The energy released from\nhelium differentiation scales with the amount of helium that rains out, which is\ncontrolled by the phase diagram demarcating immiscible versus miscible regions.\nThrough our large-scale calculations and structural analysis, we can also deduce the\ndensity ratios from the degree of demixing. This parameter can be used to\nparameterize the efficiency of heat transport by double-diffusive convection in the\nregion where helium is raining out. We hope that our presented data can be\nincorporated into future planetary physics research to enhance the evolutionary\nmodeling of planets.” in line 245.\n2) Added references\n35. Mankovich, C., Fortney, J. J. & Moore, K. L. BAYESIAN EVOLUTION\nMODELS FOR JUPITER WITH HELIUM RAIN AND DOUBLE-DIFFUSIVE\nCONVECTION. Astrophys. J. 832, 113 (2016).\n36. Mankovich, C. R. & Fortney, J. J. Evidence for a Dichotomy in the Interior\nStructures of Jupiter and Saturn from Helium Phase Separation. Astrophys. J. 889,\n51 (2020).\n----------------------------------------------------\nReply to Reviewer #2:\n----------------------------------------------------\nThis paper addresses the interior structure and evolution of cool gas giant planet like\nJupiter and Saturn through an attempt to improve the understanding of hydrogen-\nhelium mixtures at megabar pressures and moderate temperatures of several 1000 K, as\nthey are predicted to occur in the interiors of Jupiter and Saturn.\nSuch work is essential for understanding the interior structure of the two planets Jupiter\nand Saturn and important to several areas in planetary science. For instance, the mass\nof heavy elements in these two gas giants, together with the about 30 Earth masses of\nheavy elments in Uranus and Neptune, sets the minimum solar nebular mass, which in\nturn is an important constraint for the formation of the solar system. Determining the\nheavy elemet content of Jupiter and Saturn has therefore been the goal of previous\nstudies such as gravity field measurements by spacecraft and interior modeling.\nHowever, our understanding of Jupiter and Saturn including their heavy element\ncontent is still mostly limited by the uncertainty of H-He mixtures, both in regard of the\nequation of state and the immiscibility behavior. Furthermore, these two planets are\noften referred to as benchmark planets for extrasolar planets, but as long as even the\nrather basic property of heavy element content remains poorly constrained, also interior\nstructure assumptions for extrasolar giant planets remain weakly guided. Therefore, this\nwork addresses a central topic of solar system science with relations to exoplanets.\nThis work significantly exceeds the current standard of constructing H/He phase\ndiagrams. Previous approaches picked an exchange-correlational functional for the\ncomputation of the pressure and internal energy of different miscible H-He systems,\ncalibrated the resulting entropy of mixing, and used these values to compute the\npressures, entropy, energy, and finally the Gibbs free energy of mixted systems.\nDeviations from the apparent energy minimum were interpreted as sign of demixing.\nSuch approach was chosen due to the limitations in the computational cost for\nsimulating many-particle systems. In this work, systems with 9600 to 27000 particles\nare simulated compared to the previously order-of-100 particles. This improvement is\nachieved through the application of machine-learning to infer the interaction\npotential in the many-particle system. Furthermore, this study offers not only an alternative\nH/He phase diagram but presents a novel method to quantity the likelihood of demixing in\nform of statistical distribution of the neighborhod population around each atom.\nMain finding of this work are\n(i) An immiscibility gap that stretches over a large portion of the Jupiter (abstract\n'Our results revise the structures of Jupiter and Saturn where H-He demixing\ntakes place in a large fraction of the interior radii, i.e., 27.5% in Jupiter'). For\nJupiter, this would be more than twice as thick as the previously considered ~1-\n(ii) The demixing curve for about solar H/He-ratio as relevant to the planets is found\nto fall in between the previous theoretical curves on one side and the\nexperimental data-based curve on the othe side. Previous theoretical and\nexperimental curves are vastly inconsistent.\n(iii) It is demonstrated that immiscibile many-particle stystems can be directly\nsimulated if advantage is made of a machine-learning method to infer the\neffective potential in the system. This is a major progress over previous methods,\nwhich simulated mixed systems to infer immiscibility.\n(iv) A novel method is presented to analyze the mixing or demixing behavior in\nmulti-particle systems that can be applied to other materials in the future.\nThe presented work is excellent and a major progress to the field. But the paper requires\nmore clarifications to support the major conclusions. I expect that after the points listed\nbelow have been addressed in a sligtly revised version I would recommend its\npublication in a highly-ranked but more spezialized journal. At present I can neighter\ncommend it to Nature Communications nor reject it because of two major concerns I\nhave. To address them will require further work by the authors.\nResponse:\nWe thank the reviewer for their thorough evaluation of our manuscript and their further\nelaboration on the research background. We are encouraged by the reviewer’s\nrecognition of the significance of this work.\nWe are very interested in the behavior of matter under planetary conditions. We conduct\nresearch using machine learning and anticipate obtaining data from ultra-large atomic\nsystems, which can provide support for modeling real planetary systems. When facing\nthe intriguing question of whether helium rain exists within Jupiter and Saturn, we\nencounter a controversy arising from the computational challenges and the proximity\nof planetary adiabats and phase separation boundaries in the P-T diagram. However,\nthis issue, which can be directly addressed through structural analysis, is ideally suited\nfor investigations at the super-sized atomic scale. In this context, the significance of\nmachine learning potentials extends beyond mere speed enhancement; it offers an\nalternative, more intuitive perspective for analyzing the problem, as highlighted by the\nreviewer.\nWe thank the reviewer for summarizing the main findings of this work. We believe that\nthese results can provide important data support for modeling Saturn and Jupiter and\nserve as a reference for similar studies. In the revised version of the manuscript, we\nhave further emphasized these findings and provided additional analysis to strengthen\nthe credibility and significance of these conclusions. Regarding the concerns raised, we\nacknowledge the need for further clarification to support our major conclusions. We\nwill carefully consider the reviewer's suggestions to ensure that our research meets the\nhighest standards. Once again, we appreciate the reviewer's thorough review and\nvaluable feedback.\nMajor comments:\nComment 1: the DLP potential\nCurrent evidence is insufficient that the main conclusions (i,ii) are due to the choice of\nthe \"best-of-both-worlds\" features of the SCAN+rVV10 functional in conjunction with\nthe new method (iv) of analyzing immiscibility. The other way round, evidence must\nbe provided that the obtained location of the P-T immiscibility curve (conclusion i-ii)\nis not due to the up to 10-30% strong deviations of the Deep-Learning Potential (DLP)\nfrom the SCAN functional-based DFT method. The benchmarking cases in the\nSupplmenetal Figure S2 demonstrate this large deviation. Therefore, it is mandatory\nthat the uncertainty on the H/He-phase diagram from this deviation be quantified. The\nauthors can design their own evidence but here are my three suggetions to select from:\neither simulate the large (27.000 particles) system with a modified DLP-functional that\nbetter matches the SCAN-functional in terms of the atomic force and derive the\nxHe=0.089 demixing curve, or 2nd, simulate the small 64-particles systems with\nmodified SCAN-functional that resemble the DLP-based forces to obtain the demixing\ncurve, or, to display the atomic forces of PBE, vdw-DF1, and SCAN-functional and\notherwise apply the same method to compute the demixing curve to see what magnitude\nin force-difference between the established potentials leads to what magnitude of\nchange in the demixing curve.\nResponse:\nWe appreciate the reviewer for pointing out the issues in verifying DLP potential errors,\nwhich may have made it difficult to convince the readers. We highly value this feedback\nand have attempted to elaborate and demonstrate the following from multiple\nperspectives.\nRegarding the error issue illustrated in Fig S2, we have to point out that the\ndimensionless ratio is “RMSE of forces between DLP and DFT over the standard\ndeviation (std) of DFT forces”, which focuses on the comparison with fluctuations\nof real values. For a DLP, a force RMSE below the force std (RMSE/std < 1) can be\nconsidered as having some fitting capability. Achieving 10% to 20% indicates a high\nlevel of precision. From a statistical perspective, our MLP exhibits a coefficient of\ndetermination R-squared (1- (RMSE/std)^2) over 96%, which signifies accurate model\nfitting for force. Therefore, from Fig S2, it can be observed that the DLP generally\nexhibits good accuracy. We have made numerous attempts and conducted training on\nsegmented intervals. After many trials, it seems that we may have reached the limit of\npotential function accuracy for this range. Thus, we opt to provide evidence, as pointed\nout by the reviewer, addressing the core issue of whether the DLP accuracy is adequate\nfor the purposes of this study. We directly validated the phase separation results using\nDFT calculations of a large system containing 2,048 atoms.\nAdded text:\nWe uniformly selected four temperature-density points within the system described in\nFig. 3, all of which are situated within the intervals of relatively higher relative errors\nas shown in previous Fig. S2. We conducted DFT-MD simulations for a considerably\nlong time scale (2 ps) with a system comprising 2,048 atoms. For each frame of the\ntrajectory, structural analysis, as described in the main text, was performed to compare\nwith the results obtained from the DLP model, thereby assessing whether the accuracy\nof the potential function meets the requirements for addressing the core issues discussed\nin this study. As depicted in the figure below, the comparison reveals that, in\nconjunction with the validation results of the potential function's energy, forces,\nstructures, and dynamics obtained previously, the DLP results are deemed reliable. It\nshould be noted that the energy and force RMSE of DLP and DFT for the following\nfour cases are listed in the table below, which are consistent with those of the smaller\nsystems, thus the model is valid in the MD simulations of 27K atoms.\nP [Mbar] T [K] X RMSE_energy[eV/atom] RMSE_force[eV/angstrom]\nHe\n1 3000 0.073 0.005 0.300\n1.5 4000 0.073 0.010 0.303\n7 4000 0.073 0.010 0.297\n2 6000 0.2 0.011 0.337\nAdded Figure:\nFig. S9. Comparison of helium abundance differentiation achieved by\nreweighted conditional probability between DLP and DFT with a system size\nof 2,048 atoms.\nWe believe the current fitting error of potential functions can be approximated as\nincreasing the system's force fluctuations. We have statistically analyzed the force\ndistribution within a single frame, as shown in the graph below. The Gaussian\ndistribution of forces might cause the effects of fluctuations canceling out in a\nsufficiently long statistical average. Our previous work has validated a similar\nperspective. The errors in force calculations exhibit a large convergence tolerance and\ncan be treated as noise in the Langevin dynamics for constant-temperature molecular\ndynamics simulations. Even when introducing Gaussian noise, the results obtained\nfrom quantum Langevin MD remain consistent with those obtained from highly\naccurate methods, and the dynamical properties of the systems are accurately\nreproduced [Europhys. Lett. 88, 20001 (2009); Phys. Rev. Lett. 109, 175701 (2012)].\nThis demonstrates our MLP’s robustness in handling reasonable levels noise or errors\nin force calculations while preserving the fidelity of the simulated dynamical properties.\nIn addressing the hydrogen-helium phase separation issue under planetary conditions,\nwe relaxed every system for at least 100 ps and conducted long-term statistical\naveraging. Therefore, when the fitting error of the potential function is covered by the\nsystem’s thermodynamic fluctuations and remains within a relatively low range, we can\naverage out this error through the self-cancellation effect of fluctuations.\nWe consider the precision of energy fitting to be particularly sensitive when it comes to\nthermodynamic boundaries separating phases. Therefore, in training potential functions,\nwe prioritize on the accuracy of energy fitting. The graphs presented below illustrate\nthe distribution of RMSE of energies and RMSE over std of energies. These graphs\ndemonstrate consistently high fitting accuracy of energies across a broad range of\ntemperature and pressure.\nAdded Figure:\nFig. S4. The energy RMSE. The top right corner of each subplot displays their\nrespective helium abundance. The red points represent thermodynamic conditions\nsampled as testing sets.\nFig. S5. The ratio of energy RMSE between DLP and DFT over energy std of DFT.\nThe top right corner of each subplot displays their respective helium abundance. The\nred points represent thermodynamic conditions sampled as testing sets.\nWe agree with the reviewer that the simultaneous variations in both the functional and\nmethodological aspects would raise confusion and it is imperative to clarify this issue.\nWe also greatly appreciate the reviewer’s insightful identification of three feasible\napproaches. We trained the potential function based on vdw-DF within the pressure\nrange of 0.5 to 4 Mbar and conducted 27,000-atom calculations for protosolar helium\nabundance, as shown in the figure below. It can be observed that at the\nlevel, our results closely align with Schottler's results (violet line and squares),\n∆𝒙 = 𝟎.𝟏\nindicating methodological consistency at lower resolutions. In cases where the free\nenergy methods for small systems cannot provide further resolution, methods\napplicable to larger systems can reveal the weak separation of hydrogen and helium.\nAt the same time, it can be seen that the phase separation intervals obtained by the\nvdW-DF-based DLP (pink squares and black solid line) closely approximate those\nof the SCAN+rVV10-based DLP (white circles with black dashed line).\nAdded Figure:\nExtended Data Fig. 4. Miscibility diagram at helium abundance of\nof vdW-DF-based DLP. The solid black line with pink squares represents the\n𝑿𝑯𝒆 = 𝟎.𝟎𝟖𝟗\nhydrogen-helium immiscibility line of vdW-DF-based DLP. The color of the colored\narea represents the value of , and 0.0181 on the color bar is the minimum value\nthat can be considered as the sign of immiscibility for helium abundance of protosolar.\n∆𝑥\nOthers are the same as Fig. 2 in the main text.\nThe modification of the picture:\n1) Fig. 3\nAdd X = 0.089 demixing points of vdW-DF-based DLP in Fig .3a;\nHe\nAdd pressure axis in Fig .3c.\n2) Extended Data Fig. 4\nAdd Extended Data Fig. 4 displayed above.\n3) Fig. S4\nAdd Fig. S4 displayed above.\n4) Fig. S5\nAdd Fig. S5 displayed above.\n5) Fig. S9\nAdd Fig. S9 displayed below.\nThe modification of the article:\n1) Added “We apply this method to the DFT-MD trajectories of 2,048 atoms with\nhelium abundances of X = 0.073 or 0.2. In addition to validating the accuracy of\nHe\nthe DLP, it also demonstrates that the method can accurately differentiate helium\nabundance on a case-by-case basis, even when the separation levels are very close\n(see Fig .S9 in Supplementary Information Note 1)” in line 92.\n2) Added “To assess whether different functionals fundamentally influence the\nconclusions, we also trained a DLP based on vdW-DF and conducted identical\nsimulations under protosolar helium abundance, as indicated by the pink squares in\nFig. 3a. It can be observed that changing the functional still yields similar\nconclusions. The miscibility diagram calculated at protosolar helium abundance is\nprovided in Extended Data Fig. 4.” at line 210.\n3) Added paragraphs indicated in the preceding text in Supplementary Information\nline 34.\nComment 2: the homogeneous adiabat\nConclusion (i) seems to be based on the overlap region of a homogeneous Jupiter\nadiabat with the demixing curve in P-T space for solar composition as shown in Figure\n3a. While such comparisons are common in the literature, this is insufficient to\ndetermine the extent of the demixing region in the planet. In the planet, the local He\nabundance will follow the phase curve, and leave it at some P-T if the remaining He-\nabundance becomes too low (Jupiter) to stay on the phase boundary. While the high-T\nlocation of the new demixing curve suggests that this region can extend deeper, only a\nJupiter adiabat that adopts locally the He concentration on the phase boundary and\nleaves the demixing regions due to He mass conservation would yield a reliable\nestimate. If such a planet model has not been constructed conclusion (i) should be\nphased more qualitatively as a possible outcome rather than a request for 'revision' of\nJupiter models.\nResponse:\nWe agree with the cautious approach suggested by the reviewer and acknowledge the\nlimitations in our expertise regarding planetary modeling. In response to the reviewer’s\nsuggestion, conclusion (i) has been revised to adopt a more qualitative approach rather\nthan demanding revisions to planetary models. It is emphasized that the results provided\nrepresent a possible outcome rather than a definitive prediction.\nFurthermore, additional adiabats of different planetary surface temperatures are\ncompared in the Extended Data Fig. 3, to strengthen the credibility of our conclusions.\nThis will better elucidate the limitations of our study and ensure the scientific rigor of\nour conclusions.\nThe modification of the picture:\n3) Extended Data Fig. 3:\nAdd adiabats corresponding to different surface temperatures;\nAdd X = 0.089 demixing points of vdW-DF-based DLP.:\nHe\nExtended Data Fig. 3. Miscibility diagram at helium abundance of\nof SCAN+rVV10-based DLP. The solid line with black circles represent the\n𝑿𝑯𝒆 = 𝟎.𝟎𝟖𝟗\nhydrogen-helium immiscibility lines of SCAN+rVV10-based DLP . The color of the\ncolored area represents the value of , and 0.0181 on the color bar is the minimum\nvalue that can be considered as the sign of immiscibility for helium abundance of\n∆𝑥\nprotosolar. The dash black lines represent the isentropes of various planetary surface\ntemperature of protosolar helium concentration. Others are the same as Fig. 2 in the\nmain text.\nThe modification of the article:\n1) Modified the title of the paper to be “Theoretical evidence of H-He demixing under\nJupiter and Saturn conditions”\n2) Modified lines 26-27\n“Our results revise the structures of Jupiter and Saturn where H-He demixing takes\nplace in a large fraction of the interior radii.” to\n“Our results suggest a possibility that H-He demixing takes place in a large fraction\nof the interior radii of Jupiter and Saturn”.\n3) Modified lines 28-29\n“This direct evidence of an H-He immiscible layer supports the formation of helium\nrain and explains the helium reduction in atmosphere of Jupiter and Saturn” to\n“This indication of an H-He immiscible layer hints at the formation of helium rain\nand offers a potential explanation for the decrease of helium in the atmospheres of\nJupiter and Saturn.”\n4) Modified lines 218-219\n“Combining the calculated miscibility diagram with a planet model3, we derived\nthe internal structures of Jupiter and Saturn (Fig. 3c)” to\n“Combining the calculated miscibility diagram with a planet model3, and assume\nthe interior of the planet is adiabatic, we derived the nowadays internal structures\nof Jupiter and Saturn (Fig. 3c).”\n5) Modified line 227\n“In summary, we give direct evidence for the H-He immiscibility from first-\nprinciples” to\n“In summary, we present direct observations suggesting the H-He immiscibility\nfrom first-principles.”\n6) Added “While the outer molecular envelope may still be well approximated as\nadiabatic, the deep interior is expected to depart from an adiabat due to helium rain\nand associated double-diffusive convection effects. Modeling these non-adiabatic\nregions is crucial for accurately describing the H-He demixing in Jupiter and Saturn.\nIncorporating these complexities into the model can lead to temperature shifts in\nthe predicted hydrogen-helium separation region35,36. This work offers more precise\nparametric conditions for intricate planetary modeling. The ab initio H-He phase\ndiagram delineates the P-T conditions at which helium becomes immiscible with\nhydrogen. It is crucial for determining the equilibrium helium abundance,\ninfluencing the depletion of the outer molecular envelope as rainout proceeds. The\nhelium gradient, driven by the increasing helium concentration inward influences\nthe profile of the double-diffusive convection region. The energy released from\nhelium differentiation scales with the amount of helium that rains out, which is\ncontrolled by the phase diagram demarcating immiscible versus miscible regions.\nThrough our large-scale calculations and structural analysis, we can also deduce the\ndensity ratios from the degree of demixing. This parameter can be used to\nparameterize the efficiency of heat transport by double-diffusive convection in the\nregion where helium is raining out. We hope that our presented data can be\nincorporated into future planetary physics research to enhance the evolutionary\nmodeling of planets.” in line 245.\n7) Added references\n35. Mankovich"}]

Decision letter