Towards a self-consistent evaluation of gas dwarf scenarios for temperate sub-Neptunes
Frances E. Rigby Lorenzo Pica-Ciamarra Måns Holmberg Nikku Madhusudhan Savvas Constantinou Laura Schaefer Jie Deng Kanani K. M. Lee Julianne I. Moses



These authors contributed comparably to this work.
Institute of Astronomy, University of Cambridge,
Madingley Road, Cambridge CB3 0HA, UK


These authors contributed comparably to this work.
Institute of Astronomy, University of Cambridge,
Madingley Road, Cambridge CB3 0HA, UK


These authors contributed comparably to this work.
Institute of Astronomy, University of Cambridge, 
Madingley Road, Cambridge CB3 0HA, UK


These authors contributed comparably to this work.
Institute of Astronomy, University of Cambridge, 
Madingley Road, Cambridge CB3 0HA, UK


These authors contributed comparably to this work.
Institute of Astronomy, University of Cambridge, 
Madingley Road, Cambridge CB3 0HA, UK

Nikku Madhusudhan
nmadhu@ast.cam.ac.uk


Department of Geological Sciences, School of Earth, Energy, and Environmental Sciences,
 Stanford University, Stanford, CA 94305, USA


Department of Geosciences, Princeton University, Princeton, NJ 08544, USA


Department of Physics, United States Coast Guard Academy, New London, CT 06320, USA


Space Science Institute, 4765 Walnut Street, Suite B, Boulder, CO 80301, USA

§ ABSTRACT 
The recent JWST detections of carbon-bearing molecules in a habitable-zone sub-Neptune have opened a new era in the study of low-mass exoplanets. The sub-Neptune regime spans a wide diversity of planetary interiors and atmospheres not witnessed in the solar system, including mini-Neptunes, super-Earths, and water worlds. Recent works have investigated the possibility of gas dwarfs, with rocky interiors and thick H_2-rich atmospheres, to explain aspects of the sub-Neptune population, including the radius valley. Interactions between the H_2-rich envelope and a potential magma ocean may lead to observable atmospheric signatures. We report a coupled interior-atmosphere modelling framework for gas dwarfs to investigate the plausibility of magma oceans on such planets and their observable diagnostics. We find that the surface-atmosphere interactions and atmospheric composition are sensitive to a wide range of parameters, including the atmospheric and internal structure, mineral composition, volatile solubility and atmospheric chemistry. While magma oceans are typically associated with high-temperature rocky planets, we assess if such conditions may be admissible and observable for temperate sub-Neptunes. We find that a holistic modelling approach is required for this purpose and to avoid unphysical model solutions. We find using our model framework and considering the habitable-zone sub-Neptune K2-18 b as a case study that its observed atmospheric composition is incompatible with a magma ocean scenario. We identify key atmospheric molecular and elemental diagnostics, including the abundances of  CO_2, CO, NH_3 and, potentially, S-bearing species. Our study also underscores the need for fundamental material properties for accurate modelling of such planets. 




§ INTRODUCTION
 

Sub-Neptune planets, with radii 1 R_⊕≲ R_p≲ 4 R_⊕, have emerged as the new frontier of exoplanet science and constitute the most numerous class of planets detected to date (e.g., ). The nature of the sub-Neptune population remains debated, as their bulk densities can be explained by a number of degenerate interior compositions <cit.>. These include rocky planets with diverse atmospheric compositions, mini-Neptunes with volatile-rich interiors and deep H_2-rich atmospheres, and water worlds with substantial water mass fractions, including Hycean worlds <cit.>. 

The James Webb Space Telescope (JWST) is revolutionising our understanding of sub-Neptunes through high-precision atmospheric spectroscopy. JWST observations have led to confident detections and precise abundance constraints for CH_4 and CO_2 in the atmospheres of the habitable-zone sub-Neptune and candidate Hycean world <cit.> K2-18 b <cit.>, demonstrating the promise of JWST for detailed atmospheric characterisation. Furthermore, such observations are starting to be available for other temperate sub-Neptunes, including TOI-270 d <cit.> – where abundance constraints for CH_4 and CO_2 were also retrieved – and LHS 1140 b <cit.>. Such precise abundance measurements pave the way towards understanding the interactions between the planet's atmosphere and interior, including the presence and nature of an underlying surface, as well as the planetary formation processes that give rise to such planets.

One of the most distinct features of the sub-Neptune population is the radius valley, a bimodal distribution of sub-Neptune radii with a minimum around 1.8 R_⊕ <cit.>. Two competing hypotheses have been proposed to explain the origin of the radius valley. One explanation suggests that the valley is a consequence of differential atmospheric mass loss between planets of different masses. In this hypothesis, both populations would be composed of planets with predominantly rocky interiors. The more massive planets would retain their primary H_2-rich atmospheres, while the less massive ones would instead largely lose their envelope and hence have a smaller radius. We refer to the larger population, with rocky interiors and a deep H2-rich atmospheres, as gas dwarfs. The mechanism for the mass loss is debated, with the predictions of two hypotheses – photoevaporation <cit.> and core-powered mass loss <cit.> – both proposed to explain the observations <cit.>. A second explanation <cit.> suggests that the valley could instead be due to planets having different interior compositions. The smaller radius population would be rocky, as in the atmospheric mass-loss scenario, while the larger population would be composed of planets with water-rich interiors due to significant accumulation of icy planetesimals/pebbles during their formation and migration. Atmospheric observations of planets in the sub-Neptune range may be able to distinguish between these two scenarios <cit.>. 

While the gas dwarf hypothesis has garnered significant attention in the literature <cit.>, several open questions remain. Firstly, it is unclear whether it is possible for rocky cores to accrete a substantial H_2-rich envelope without significant accretion of other volatiles and ices <cit.>. Secondly, should such planets exist, would the atmosphere-interior interactions give rise to distinct atmospheric signatures? This might be expected if the rocky surface were to be molten, giving rise to a magma ocean scenario <cit.>. It is, however, not fully clear whether this scenario is possible, particularly for planets with a low equilibrium temperature. For these planets, only a subset of atmospheric structures, combining sufficient but not exceedingly high surface pressure and very high surface temperature, could result in magma at the base of the atmosphere. 


Several recent studies have explored the implications of a magma ocean on the atmosphere and interior compositions of diverse planets, both with terrestrial-like <cit.> and H_2-rich atmospheres <cit.>. These works identify several key factors, including temperature and oxygen fugacity at the bottom of the atmosphere, that influence the composition of the atmosphere, driven by thermochemical equilibrium at the gas-melt interface. For example, some notable atmospheric signatures of reduced conditions in a rocky interior include potential nitrogen depletion <cit.> and high CO/CO_2 ratio for H_2-rich atmospheres <cit.>. However, the interplay between the atmosphere, interior, and the corresponding surface-atmosphere interactions in sub-Neptunes is only beginning to be explored in a realistic manner <cit.>. 


In this work, we develop an integrated magma ocean framework for temperate, H_2-rich sub-Neptunes. Our framework, presented in Section <ref>, includes atmospheric and internal structure modelling, melt-gas interactions, and both equilibrium and disequilibrium processes in the atmosphere, resulting in spectroscopic predictions of atmospheric observables. We consider thermochemical equilibrium at the magma-atmosphere interface, and the solubility of volatile (H, C, N, O, S) bearing species in magma. We explore the extreme case of the habitable-zone sub-Neptune and Hycean candidate K2-18 b <cit.> to investigate the plausibility of a magma ocean <cit.> and, if present, its atmospheric signatures. In doing so, we first use our framework to perform a comparative assessment of previous works in this direction in Section <ref>, both on terrestrial-like atmospheres <cit.> and on H2-rich ones <cit.>, with a focus on the case study of the candidate Hycean world K2-18 b. We then present our model predictions in Section <ref>. Finally, we summarize our findings and discuss future work in Section <ref>, highlighting the need for physically consistent models, and new experimental and theoretical work to derive accurate fundamental material properties. 




§ METHODS
 
We develop an integrated modelling framework to evaluate gas dwarf scenarios for planets in the sub-Neptune regime. A schematic flowchart of the framework is shown in  Figure <ref>.
We start by considering the constraints that the observed bulk parameters (mass, radius, and hence density) and known atmospheric properties impose on the planet's atmospheric and internal structure. This enables us to infer the possible conditions at the surface-atmosphere boundary, and, by considering a relevant mineral phase diagram, assess whether such conditions can in principle lead to a magma ocean scenario. If they can, we proceed by modelling the chemistry at the magma-atmosphere interface, which is determined by equilibrium processes including the solubility of relevant volatiles in the silicate melt, providing us with the elemental abundances in the gas phase at the interface. 
These are then evolved to the rest of the atmosphere, assuming chemical equilibrium in the lower atmosphere, and non-equilibrium processes (photochemistry and vertical mixing) in the upper atmosphere. 
This allows us to compute the observable composition of the atmosphere, which can be compared with the molecular abundances retrieved through observations to finally assess the plausibility of a magma ocean scenario for the planet.
We now describe in detail each of the steps outlined above. 





 §.§ Atmospheric Structure and Composition


We begin by modelling the atmospheric temperature structure in a self-consistent manner. In order to do so, the atmospheric chemical composition needs to be assumed. This can be done either by assuming the elemental abundances and atmospheric chemistry, or by directly assuming the molecular mixing ratios in the atmosphere. Other parameters that need to be taken into account include the internal temperature T_int representing an internal heat flux, the incident irradiation, the stellar properties, the presence and characteristics of clouds/hazes in the planet's atmosphere, and the efficiency of day-night energy redistribution. The self-consistent calculation will yield a pressure-temperature (P-T) profile, which will be coupled to the internal structure model, as discussed in Section <ref>.

In order to carry out the self-consistent modelling of the atmospheric structure, we use the GENESIS framework <cit.> adapted for sub-Neptunes (). GENESIS solves for radiative-convective equilibrium throughout the atmosphere, which is assumed to be plane-parallel, using the Rybicki scheme. It carries out line-by-line radiative transfer calculations through the Feautrier method <cit.> and the discontinuous finite element method <cit.>, while taking into account all of the parameters mentioned earlier in this section.  

For the atmospheric composition, we adopt uniform mixing ratios of molecular species based on the retrieved values at the terminator region of K2-18 b <cit.>. We use the median retrieved abundances for the one-offset case: logX_CH4 = -1.72 and logX_CO2 = -2.04. For H2O, we consider the 95% one-offset upper limit, logX_H2O = -3.01. We also assume the incident irradiation and stellar properties of K2-18 b, and uniform day-night energy redistribution. We then explore the remaining parameter space. In particular, we consider two end-member values for T_int, 25 K and 50 K, following <cit.> and <cit.>, and three values for a, the hazes' Rayleigh enhancement factor: 100, 1500 and 10000. We consider four combinations of these parameters, obtaining a cold case (designated C1, corresponding to T_int = 25 K, a = 10000), two canonical cases (both with a = 1500, designated C2 for T_int = 25 K and C3 for T_int = 50 K) and a hot case (C4, with T_int = 50 K and a = 100). These profiles are shown in Figure <ref>.

We place the upper boundary of the atmosphere at 10^-6 bar, and calculate the P-T profile self-consistently down to 10^3 bar, below the radiative-convective boundary. At higher pressures, we extrapolate the profile as an adiabat, using the H2/He equation of state (EOS), ρ=ρ(P,T), and adiabatic gradient from <cit.>. 
We note that, in principle, an appropriate P-T profile may be even colder than C1, considering the constraints on clouds/hazes at the terminator from observations of K2-18 b <cit.>.




 §.§ Internal Structure Modelling
 



We model planetary internal structures using the HyRIS framework, outlined in <cit.>. The model calculates the planet radius (R_p) from the planet mass (M_p), the mass fractions of the planet's components (x_i=M_i/M_p), and the corresponding EOS and P-T profile. HyRIS solves the equations for mass continuity and hydrostatic equilibrium using a fourth-order Runge-Kutta method, and solves for R_p using a bisection procedure. For the purpose of this study investigating magma-ocean scenarios, the internal structure model includes a H_2-rich envelope, a silicate mantle, and an iron core. 

The silicate mantle is described by EOSs valid for the liquid and solid phases – for simplicity, we adopt a separate EOS prescription on either side of a melting curve. The composition is nominally assumed to be peridotitic. The magma is described by an EOS for peridotitic melt compiled similarly to <cit.> by combining the densities of molten enstatite, forsterite, fayalite, anorthite and diopside, described by third-order Birch-Murnaghan/Mie-Gruneisen EOSs from <cit.>, weighted by their mass fractions. For the purpose of this initial study, we assume complete melting occurs at the liquidus, and hence do not include an EOS prescription for the partial melt between the solidus and liquidus curves. We use the peridotite liquidus from <cit.>, based on <cit.> – both the liquidus and solidus are shown in Figure <ref> <cit.>. The solid portion of the silicate mantle is described by the EOS of <cit.> for the high-pressure peridotite assemblage. At extreme mantle pressures beyond the pressure range of these experiments (107 GPa), we use the temperature-independent EOS of <cit.> for MgSiO_3 perovskite, originally derived at room temperature. The thermal effects for solid silicates at these pressures are small <cit.> with negligible effect on the internal structure. The iron core is described by the EOS of <cit.> for hexagonal close-packed Fe.

The temperature structure in the melt is assumed to be adiabatic. The adiabatic gradient is calculated using the specific heat for peridotite from <cit.> and the volume expansion coefficient that we calculate from the combined peridotite melt EOS. The adiabatic gradient in the upper portion of the solid mantle is calculated following <cit.>. Following previous studies <cit.>, the remaining solid portion of the interior is taken to be isothermal, as the EOSs used are temperature-independent <cit.>.  

The mass of the magma ocean follows from the adiabatic temperature profile in the melt, similarly to the calculation of water ocean depths by <cit.> and <cit.>. The melt adiabat and hence the magma base pressure are defined by the surface pressure and temperature. For a given interior composition and surface conditions, the mass of the melt can thus be calculated. We adapt HyRIS to automate the extraction of the relevant melt characteristics, similar to the methods for water oceans in <cit.>. The mass fraction of the melt is an important quantity for considerations of the available volatile reservoir, as discussed below. We note that the moderate increase of the magma ocean mass fraction that may result from partial melting is partly accounted for by our range of considered melt masses in Section <ref>.





 §.§ Melt-atmosphere Interface Chemistry
 

The atmospheric chemistry is constrained by the elemental composition at the bottom of the atmosphere, which is governed by the interactions at the magma ocean and atmosphere interface. At this boundary, we model the reactions and solubility of the gas species in thermochemical equilibrium. We include 82 H-C-N-O-S gas species and He, the set of which we denote 𝐗, and their equilibrium reactions, nominally excluding other effects such as condensation and exsolution. Of these volatile species, we consider the solubility in the melt of H_2 <cit.>, H_2O <cit.>, CO <cit.>, CO_2 <cit.>, CH_4 <cit.>, N_2 <cit.>, S_2 <cit.> and H2S <cit.>, as further motivated in Appendix <ref>. We note that the solubility of H2S is uncertain at high temperatures/pressures and may be higher if, for example, its solubility approaches that of S2. Furthermore, we remark that we are not considering the possible exsolution of FeS, which may affect the abundance of sulfur in the atmosphere. Likewise, the overall solubility of nitrogen is calculated here through N2, and may be higher if the solubility of NH3 is significant. The data on NH3 solubility in magma is currently limited and it is difficult to make any quantitative estimates of NH3 solubility. For the explicitly composition-dependent laws, we use the <cit.> Etna basalt melt composition. Similarly to <cit.>, we assume that the magma is well-stirred such that the equilibration at the surface sets the volatile abundance throughout the melt.

These solubility laws relate the partial pressures in the atmosphere to the concentrations of the volatiles in the melt. The amount of volatiles in the melt thus depends on the equilibrium chemistry, the solubility and the total mass of the melt, M_melt. For a given mass of the atmosphere and the melt, we have the following mass balance condition for each species i <cit.>,


    M_tot  w_i = M_atm  w_i, atm + M_melt  w_i, melt ,

where w_i is the total mass fraction of each species i.

To determine the chemical composition of the atmosphere and the melt, we solve the element conservation equations

    ε_j = ∑_i  ∈ 𝐗ν_ijn_i/n_⟨ H ⟩ ,

where n_i is the total amount of moles of species i, n_⟨ H ⟩ = n_H + 2n_H_2 + 2n_H_2O + … is the total amount of moles of hydrogen, ν_ij are the coefficients of the stoichiometric matrix, and ε_j is the elemental abundance of element j relative to hydrogen. Equation (<ref>) is coupled to Equation (<ref>) via n_i ∝ w_i / μ_i, where μ_i is the molar mass, which in turn is coupled to the law of mass action

    p_i/p^ = K_i ∏_j  ∈   E(  p_j/p^)^ν_ij ,

and the solubility laws, determining both w_i, atm and w_i, melt. Here, E is the set of all elements, p_i is the partial pressure of species i, K_i is the temperature-dependent equilibrium constant, and p^ is a standard pressure of 1 bar. For each gas species, we approximate the equilibrium constant as

    ln K(T) = a_0/T + a_1 ln T + b_0 + b_1 T + b_2 T^2 ,

using the coefficients provided by FASTCHEM <cit.>, mainly derived using thermochemical data from <cit.>.

Overall, Equation (<ref>) depends on the elemental partial pressures, with 6 unknowns, corresponding to the 6 elements considered. Nominally, we solve for these using the 5 equations in (<ref>) for all elements apart from H, together with

    P_s = ∑_i  ∈ 𝐗 p_i ,

to fix the total pressure. This treatment of oxygen yields a first-order estimate of the redox state as set by the atmosphere. Alternately, we consider oxygen fugacity (f_O2) as a free parameter, by determining p_O in Equation (<ref>) via f_O2 =  p_O2 = K_O2  p_O^2  / p^, allowing us to consider different redox conditions. In this framework, we assume ideal gas behavior such that fugacity and partial pressure are equivalent <cit.>. 

As a cross-check, we validate our new framework against a self-consistent atmosphere composition model <cit.> which uses the Gibbs energy minimization code IVTANTHERMO <cit.>. IVTANTHERMO uses a thermodynamic database based on <cit.>, which we modify to include the silicate-melt dissolved volatile species H_2, OH^-, O^2-, CO, CO_2, CH_4, N^3-, and S^2-. We calculate equilibrium between a total possible 366 gas species and 201 condensed species. For the dissolution reactions, we assume ΔC_p = 0 and that any temperature dependence in the equilibrium constant is due to the heat of the reaction. However, data is available only in limited temperature ranges for most dissolution reactions, so we assume a simple Henry's law solubility relation for all of the dissolved species except S^2-, OH^-, H_2, and CH_4. We also neglect non-ideality in both the gas phase and melt. Using IVTANTHERMO, we then compute self-consistent equilibrium between the gas phase and melt species as a function of pressure and temperature.

For this comparison, we use 50×solar bulk elemental abundances (not including He), P_s=10^4 bar, T_s=3000 K, and M_melt / M_atm = 0.20, which is given by the gas-to-melt mass ratio as calculated by IVTANTHERMO. We find that all major H-C-N-O-S gas species agree to within at most 0.35 dex (standard deviation of 0.1 dex), with the largest deviation coming from CO2. This deviation mostly stems from the oxygen fugacities being somewhat different between the two approaches, with IVTANTHERMO yielding a 0.35 dex lower value. Furthermore, we verify that we recover the atmospheric abundances given by FASTCHEM 2 <cit.> and GGCHEM <cit.> when setting M_melt = 0.




 §.§ Atmospheric Chemistry


We carry out equilibrium and disequilibrium chemistry calculations to determine the atmospheric composition above the magma/rock surface. We use the VULCAN photochemical kinetics framework <cit.>, with the initial atmospheric chemistry obtained using the FASTCHEM equilibrium chemistry code <cit.>.

For equilibrium chemistry calculations, we consider thermochemical equilibrium involving H-C-N-O-S species as well as He, along with H_2O condensation. For calculations considering disequilibrium processes, we additionally include the effects of vertical mixing and photochemistry. We follow the K_zz parameterisation of <cit.>:

    K_zz / cm^2 s^-1 = 
    min(5.6 × 10^4/ (P/bar)^1/2, 10^10),     P ≤ 0.5 bar
    
        10^6,     P > 0.5 bar,

although we note that the K_zz in the troposphere could be higher (e.g., ∼ 10^7-10^8 cm^2 s^-1 in the deep convective region of the atmosphere) or lower (e.g., ∼10^4 cm^2 s^-1) in any radiative regions if moist convection is inhibited by the high molecular weight of water in the H_2-rich atmosphere (see ). Accordingly, we consider a wider range of K_zz values than our canonical treatment in Section <ref> and Appendix <ref>.

Additionally, we consider photochemical reactions including H-C-N-O-S species, using a nominal stellar spectrum from the HAZMAT spectral library <cit.> corresponding to a median 5 Gyr star of radius 0.45 R_⊙ following previous work <cit.>. We also specifically consider the condensation of H_2O to liquid and solid droplets, which fall at their terminal velocity, as described in <cit.>. We note that while the H-C-N-O chemistry has been extensively explored for sub-Neptunes in various studies <cit.>, the S chemistry has not been explored in significant detail and may be incomplete. Nevertheless, we include S using the VULCAN framework <cit.> for completeness.

With the above calculations we obtain the vertical mixing ratio profiles for a number of relevant chemical species in the atmosphere. The abundances of key species in the observable part of the atmosphere can then be compared against constraint retrieved from an atmospheric spectrum. 



 §.§ Spectral Characteristics



We use the results of the chemistry calculation described in Section <ref> to simulate how such an atmosphere would appear in transmission spectroscopy, including the spectral contributions of relevant species. For this, we use the forward model generating component of the VIRA retrieval framework <cit.>, which treats the planet’s terminator as a 1D atmosphere in hydrostatic equilibrium. We consider atmospheric opacity contributions from H_2O <cit.>, CH_4 <cit.>, NH_3 <cit.>, CO <cit.>, CO_2 <cit.>, C_2H_2 <cit.>, HCN <cit.>, H_2S <cit.> and SO_2 <cit.>. We do not include N_2 in the model, as it has no significant absorption features in the near-infrared and it is not present in significant enough quantities to affect the atmospheric mean molecular weight. We additionally consider atmospheric extinction arising from H_2-H_2 and H_2-He collision-induced absorption <cit.>, which provide the spectral baseline, as well as H_2 Rayleigh scattering. We simulate transmission spectra using the vertical mixing ratio profiles computed using VULCAN as described above, and the P-T profile appropriate to each case considered.  



§ RESULTS: COMPARISON WITH PREVIOUS WORK
 

We now apply the framework described in Section <ref> and compare with previous works on both terrestrial-like and sub-Neptune atmospheres.




 §.§ Terrestrial-like Atmospheres
 

Many previous studies have investigated surface-atmosphere interactions for magma oceans underneath terrestrial-like atmospheres <cit.>. Recent studies have explored the implications of diverse interiors of exoplanets for their atmospheric compositions. <cit.> proposed that, for reduced conditions, nitrogen is expected to be preferentially sequestered in the mantle, providing a valuable way to study the interior composition of such exoplanets. More recently, <cit.> investigated the primordial distribution of volatiles within the framework of melt-atmosphere interactions and discussed applications for Venus and Earth. For the early Earth, they find that reduced conditions, with oxygen fugacity two dex below the iron-wüstite (IW) buffer, f_O2≲IW -2, result in an atmosphere abundant in H2, CO and CH4 but depleted in CO2 and N2. On the other hand, for f_O2≳IW+2, CO2 becomes the main atmospheric component, with significant levels of SO2, N2 and H2O. In particular, the behaviour of nitrogen is a consequence of the high solubility of N2 as N^3- in silicate melt at reducing conditions <cit.>, via the following reaction

    1/2N2_ (gas) + 3/2O^2-_ (melt) ⇌ N^3-_ (melt) + 3/4O2_ (gas) .

As a result, the melt concentration of N^3- is proportional to f_N2^1/2 f_O2^-3/4, thus favouring low f_O2. In conclusion, these works predict that the abundance of atmospheric nitrogen may be used as a diagnostic for the redox state of a rocky planet's mantle.

As a benchmark, we compare our melt-atmosphere equilibrium chemistry framework with <cit.>. We use their case with a magma ocean mass of half the bulk silicate mantle at T = 1773 K and with volatile contents of 90, 102, 3.3, and 126 ppm-wt for C, H, N, and S, respectively. With this, we reproduce their atmospheric composition as shown in Figure <ref>. Compared to our nominal setup in Section <ref>, we added a constraint for the hydrogen abundance and solved for the resulting mass of the atmosphere, coupled to the surface pressure using <cit.>

    P_s = g M_atm/4 π R_p^2 ,

where g is the gravitational acceleration at R_p. For a like-to-like comparison, we added the condensation of graphite and used the same gas species (excluding Ar) and solubility laws as <cit.>. Overall, we find good agreement between both implementations, with the most deviation coming from N2 and CH4. We find that the N2 discrepancy comes from an inconsistency in the code by <cit.>, whereby they use a molar mass of 14 g/mol for N2 instead of 28 g/mol. The remaining discrepancy is likely a result of minor differences in the implementations of the different reactions. We find that by accounting for some of these differences we can better match the result by <cit.>, as shown in Figure <ref>. For this purpose, in addition to considering their adopted molar mass, we implemented the reactions CH4 + 2 O2⇌ 2 H2O + CO2 and H2O⇌ 0.5 O2 + H2 using the equilibrium constants by <cit.> to obtain the partial pressures of CH4 and H2O, instead of deriving these from the elemental partial pressures as described in Section <ref>. We also used the oxygen fugacity of the IW buffer from <cit.> instead of <cit.>.







 §.§ Sub-Neptunes with H_2-rich Atmospheres


Several recent studies have also explored magma-atmosphere interactions in sub-Neptunes with rocky interiors and H_2-rich atmospheres. <cit.> considered the impact of H2 solubility in silicate melts on the radius distribution of sub-Neptunes, addressing the radius cliff, a sharp decline in the abundance of planets with R_p ≳ 3 R_⊕. 
They find that the high solubility of H2 in magma, especially at high pressure, limits the maximum radius that can be attained by sub-Neptunes through accretion of atmospheric H2.
For a 10 M_⊕ core they find a limiting mass fraction of 1.5 wt% H2 in the atmosphere – corresponding to > 20 wt% H2 in the planet –, as any additional H2 would be stored almost exclusively in the interior. Looking at smaller planets (2 R_⊕≤ R_p≤ 3 R_⊕), <cit.> find that magma-atmosphere interactions would significantly affect the atmosphere's composition and mass. For example, a key insight is that the H2O/H2 ratio in the atmosphere reflects not only external water delivery, but also water production as a result of atmosphere-magma interactions. This would make the H2O/H2 a good diagnostic for atmospheric origin, as well as for magma composition. In particular, it is found to be proportional to the magma FeO content. 

Further investigations were carried out by <cit.>, <cit.> and, most recently, <cit.>.
Considering a surface temperature T_ s = 4500 K, 1 % to 14 % H mass fractions (of overall planet mass) and model parameters resulting in f_O2≲IW - 2, <cit.> find that the atmosphere is expected to be dominated across the explored parameter space by H2, SiO, CO, Mg and Na, followed by H2O, which should exceed CO2 and CH4 by two to three orders of magnitude. It should be noted they do not include N in their model. <cit.> and <cit.>, instead, consider total hydrogen pressures ranging between 10^-6 and 10^6 bar, temperatures between 1800 and 3500 K, and do not include any volatiles in their calculations, but also show that detectable absorption features of H2O and SiO should be expected. Additionally, the volatile-free investigation by <cit.> finds that silane (SiH4) should also be expected, dominating over SiO at P ≳ 0.1 bar for an isothermal T = 1000 K pressure-temperature (P-T) profile in the upper atmosphere. Most recently, <cit.> also investigated the outgassing mechanism for hybrid atmospheres in sub-Neptunes, but without considering solubilities in magma.




 §.§ End-member Scenario of K2-18 b
 

Some of the principles described above were recently applied to the habitable-zone sub-Neptune K2-18 b by <cit.>, hereafter S24. Similarly to <cit.> and <cit.>, S24 point to a high CO/CO2 ratio and, like <cit.> and <cit.>, a depletion in atmospheric N as signatures for the presence of a magma ocean and/or a reduced interior. It should be noted that the case of K2-18 b constitutes an end-member scenario. While most of the work on magma oceans has focused on very hot planets <cit.>, K2-18 b is a temperate sub-Neptune with equilibrium temperature T_eq = 272 K (assuming an albedo of 0.3), close to that of the Earth. Here, we assess the findings of S24 using the framework described in Section <ref> and Figure <ref>.

We briefly note that in addition to gas dwarf and Hycean world scenarios, a mini-Neptune scenario with a thick H_2-rich atmosphere has also been proposed for K2-18 b <cit.>. <cit.> conduct photochemical modelling of mini-Neptune cases for K2-18 b, suggesting a plausible solution. However, as noted in <cit.>, the calculated abundances are unable to match the retrieved abundances <cit.>. In particular, the mixing ratios of CO and NH_3 are too large compared to the retrieved abundances, and so is the CO/CO_2 ratio.



  §.§.§ Consistency with bulk parameters



At the outset, it is important to ensure that any assumption about the internal structure is consistent with the planetary bulk parameters. Previous studies have shown that the bulk parameters of K2-18 b allow a degenerate set of solutions between a mini-Neptune, a Hycean world, or a rocky world with a thick H_2-rich atmosphere, i.e. a gas dwarf <cit.>. Considering the present gas dwarf scenario, a purely rocky interior would require a minimum H_2-rich envelope mass fraction of ∼1% <cit.>, as discussed below.

The model grid of S24 contains four values of mantle mass fraction relative to the total planet mass (0.001, 0.01, 0.1, and 1) and five values of the hydrogen mass fraction relative to the mantle mass (1, 10, 100, 1000 and 10000 ppm). Firstly, all the cases with a mantle mass fraction of 1 violate mass balance, as the sum of the mantle and atmospheric masses would exceed the total planet mass. Secondly, for the gas dwarf scenario, as noted above, the bulk density of K2-18 b requires an H_2-rich atmosphere with a minimum mass fraction of ∼1%. In the S24 model grid, there is only one model which has an atmospheric mass fraction of 1%, and it corresponds to a mantle mass fraction of 1, as noted above. It follows that all the remaining cases, with H_2 mass fraction below 1%, are incompatible with the planet's bulk density. 

In order to estimate the allowed atmospheric mass fractions for K2-18 b in the gas dwarf scenario, we consider four possible interior compositions, illustrated in Table <ref>: f_silicate=100 %, Earth-like (f_silicate=67 %), Mercury-like (f_silicate=30 %) and f_silicate=5 %, where f_silicate is the mass fraction of the interior (i.e. excluding the envelope) in the silicate mantle. We include f_silicate=5 % as an end-member case, close to the upper limit for the allowed envelope mass fraction. Similarly, the extreme pure-silicate interior case is included as an end-member, yielding the lower limit on the allowed envelope mass fraction for a gas dwarf scenario. We adopt the median planetary mass M_p=8.63 M_⊕ <cit.> and radius R_p=2.61 R_⊕ <cit.> of K2-18 b. The allowed envelope mass also depends on the choice of P-T profile, with hotter profiles leading to lower envelope masses for a given interior composition, as shown in Table <ref>. 

Considering the four self-consistent P-T profiles described in Section <ref>, we find that an envelope mass fraction x_env≥ 1.34 % is required for consistency with the bulk parameters. This limit corresponds to the extreme case of a 100 % silicate interior, adopting C4 for the envelope P-T profile. For a like-to-like comparison with the S24 model grid, we also consider their P-T profile, which is the profile from <cit.> log-linearly extrapolated to higher pressures. For this profile, we find envelope mass fractions of x_env≥ 0.90% are required, again corresponding to the extreme 100 % silicate interior case. Overall, we find that all the models in the model grid of S24 are incompatible with mass balance and/or the bulk density of the planet considered. We demonstrate a self-consistent approach of accounting for the observed bulk parameters of K2-18 b in such calculations in Section <ref>.




  §.§.§ Feasibility of a magma ocean
 

As described in Section <ref> and shown in Figure <ref>, given an interior composition, the choice of P-T profile affects the resulting envelope mass fraction. This, in turn, determines the surface pressure and temperature and the liquid/solid phase of the rocky surface underneath. Therefore, it is important to consider a physically motivated P-T profile in the envelope. As mentioned above, S24 consider the P-T profile from <cit.> at low pressures (P ≤ 4 bar) and perform a log-linear extrapolation to the deep atmosphere (P ≳ 10^5 bar). The resulting temperature gradient can be significantly different from other self-consistent model P-T profiles for the H_2-rich envelope <cit.>; an example is shown in Figure <ref>.

We also note, however, that the actual surface temperature at the magma-atmosphere interface (T_ s) used in S24 appears to be a free parameter rather than self-consistently determined from their P-T profile. The T_ s ranges between 1500 K and 3000 K, but the corresponding pressure is not clear, considering their assertion that the maximum surface pressure allowed by the model is 10^8 bar. This pressure also appears to be inconsistent with their maximum envelope mass fraction of 1%. Across the range of rocky compositions we consider, the maximum pressure reached is ∼5-7×10^5 bar for envelope mass fractions ∼5-7% depending on the P-T profile as shown in Table <ref> and Figure <ref>. 

Nevertheless, in order to establish the feasibility of achieving melt conditions in the S24 model, we consider the five highest envelope mass fractions used in S24. We adopt their mantle mass fraction of 1 and the corresponding five H2 mass fractions in their model grid, with a maximum of 1%. We then use these envelope mass fractions and the S24 P-T profile to determine the corresponding expected surface pressures and temperatures, independent of satisfying the planetary bulk properties. These model points are shown in Figure <ref> along with the liquidus and solidus curves for peridotite <cit.>. We find that only two of these five cases result in a magma surface in our framework. Finally, since we considered only the five highest envelope mass fractions of S24, it follows that all of the other models would also be unlikely to result in melt. We further note that for the two cases that result in a magma surface in S24, the magma mass fraction they consider is equal to the planet mass. However, based on the temperature structures shown in Figure <ref>, we find that the maximum magma mass fraction across the different interior compositions is ∼13%, potentially somewhat higher as a result of partial melting, but not 100%. 





  §.§.§ Magma-atmosphere interactions


If the plausibility of a magma ocean is established, the melt-atmosphere interaction must be considered to determine its effect on the atmospheric composition. As described in Section <ref>, the gas phase composition depends on the pressure and temperature at the interface, the elemental abundances, the amount of magma available, the solubilities of the chemical species, and the chemical properties of the melt. 

For the case of K2-18 b, S24 consider oxygen fugacity as a free parameter and assess the abundances of several H-C-N-O species in the lower atmosphere following melt-atmosphere interactions. They determine the atmospheric composition by considering three reactions, CO2 + 2 H2⇌CH4 + O2, 2CO2⇌ 2 CO + O2, and 2H2O⇌ 2 H2 + O2, in thermochemical equilibrium, and solubilities of CH_4, N_2, CO_2, and H_2O in the magma. However, we note that these reactions do not encompass all the prominent H-C-N-O molecules at the considered conditions. In particular, NH_3 is expected to be the dominant N-bearing species at the base of the atmosphere. 
By not including NH_3 and its equilibrium with N_2 and H_2, S24 may be overestimating the nitrogen depletion in the atmosphere, given that all of the nitrogen is assumed to be in N_2, which is very soluble in magma at reducing conditions, as we show in Figure <ref> in Appendix <ref>.

In our framework, described in Section <ref>, we find that nitrogen depletion in the atmosphere increases by several orders of magnitude by not including NH_3. Ultimately, this highlights the importance of the completeness of the reactions and solubilities considered. Finally, we note that it is also possible to not have significant N depletion even in the presence of a molten surface depending on the pressure and temperature, as shown in Table <ref>.



  §.§.§ Atmospheric composition and observables

The properties at the surface determine the composition in the upper layers of the atmosphere, and hence its observable characteristics. These are strongly influenced by model assumptions on elemental abundances. S24 allow the C/H ratio to vary between 0.01 × solar and 100 × solar, while keeping the N/H ratio fixed to solar, i.e., N/H = 6.76 × 10^-5 by number. This itself limits the NH3 log-mixing ratio to at most logX_NH3∼ -4, close to the upper bound of -4.46 found by <cit.>, and biases the model by construction to allow for up to 100 × more (or down to 100 × fewer) C-based molecules than N-based ones. The dependence of the S24 model outcomes on the choice of C/H values is not reported. It should be noted that a 100× enhancement or depletion of C/H without any change in N may be difficult to reconcile with potential formation mechanisms.

We note two further points regarding the abundance of C- and N-bearing species predicted by S24. Firstly, S24 appear to indicate that the total abundances of carbon in their models reach up to 3.8 wt% of the planet mass. It is, however, unclear how this may be compatible with their assumptions of a C/H ratio of at most 100 × solar  and an H mass fraction ≤ 1 %, given they adopt the <cit.> value for (C/H)_⊙, i.e., 3.2 × 10^-3 by mass. Secondly, as argued in Section <ref>, only the largest atmospheric mass fractions S24 consider can potentially lead to a magma ocean. At the resulting surface pressures, however, their model predicts a log-mixing ratio for CO2 of log X_CO2≲ -3. This is at the lowest end, if not outside, of the 1 σ confidence interval presented in <cit.>. Furthermore, the CO abundance or the CO2/CO ratio are not reported in S24, making it difficult to assess the validity of the chemical estimates. 

Finally, S24 argue that the model spectra from their model ensemble provide a qualitatively reasonable match to the data. Even if the model spectra were taken at face value, the lack of a reported goodness-of-fit metric precludes a reliable assessment of the match to data. More generally, a limited grid of forward models is insufficient to robustly explore the full model space taking into account all the degeneracies involved in an atmospheric spectral model and to obtain a statistically robust fit to the data; that is the purpose of atmospheric retrievals <cit.>. A more reliable approach in the present context is to compare the model-predicted chemical abundances with the abundance constraints obtained from robust atmospheric retrievals of the observed spectra. As discussed above, the cases of S24 with the highest surface pressure, i.e. those that may allow a magma surface, still predict lower CO_2 abundances than those retrieved for K2-18 b <cit.>. The CO and H_2O abundances are not reported in S24, which prevents a clear assessment of the agreement between the chemical predictions and the retrieved abundances. 



§ RESULTS: A CASE STUDY OF K2-18 B
 

After having established the consistency of our results with <cit.>, and having discussed the S24 findings for K2-18 b, we proceed to apply our framework ex novo. We do so for K2-18 b in the present section, starting, as outlined in Figure <ref>, with internal and atmospheric structure modelling that ensures consistency with the known bulk parameters. 
Through considering magma-atmosphere interactions, equilibrium chemistry in the lower atmosphere and non-equilibrium processes in the upper atmosphere, we make predictions for the observable composition and spectral signatures of a sub-Neptune magma world.



 §.§ Atmospheric Structure

As discussed in Section <ref>, the dayside atmospheric structure is calculated self-consistently from the atmospheric constraints retrieved in the one-offset case of <cit.>: the median  logX_CH4 = -1.74, logX_CO2 = -2.04, and the 2σ upper bound logX_H2O = -3.01. The P-T profile depends on a wide range of parameters, not all of which are observationally well-constrained: these include the internal temperature T_int, the properties of clouds/hazes if present, and the efficiency of day-night heat redistribution. A detailed exploration of the temperature profiles in deep H2-rich sub-Neptune atmospheres has been carried out before, in <cit.>. Here, we assume uniform day-night heat redistribution, and consider four cases for the P-T profiles, varying the internal temperature T_int and the Rayleigh enhancement factor (a) for the hazes: C1, corresponding to T_int = 25 K, a = 10000; C2 and C3, both with a = 1500, with T_int = 25 K and T_int = 50 K, respectively; C4, with T_int = 50 K and a = 100. We note that, in principle, even colder profiles are plausible, given the clouds/haze properties retrieved from observations <cit.>. All the P-T profiles are shown in Figure <ref>.


 


 §.§ Internal Structure




For each of these profiles, we obtain the permitted H_2-rich envelope mass fraction (x_env) and corresponding surface conditions (P_s, T_s) based on the bulk properties of the planet, as discussed in Sections <ref> and <ref> and shown in Table <ref>. We vary the interior composition from f_silicate=5% to f_silicate=100%, adopting the median M_p=8.63 M_⊕ <cit.> and R_p=2.61 R_⊕ <cit.>. We note that the pure silicate and 95% iron (f_silicate=5%) interior cases are unrealistic end-member interior compositions, but we consider them for completeness. We adopt P_0=0.05 bar as the outer boundary condition for the internal structure modelling, corresponding to the pressure at R_p, based on <cit.>.

In Figure <ref> we show the P-T profiles considered, along with the surface conditions (black circles) and adiabatic profiles in the melt for our nominal C2 and C3 scenarios, which we further discuss below. The results for all P-T profiles are given in Table <ref>.

The presence and amount of magma depend on the adopted P-T profile. We start by considering one of the colder profiles, C2. For an Earth-like interior, we find x_env=3.76%, with surface conditions P_s=2.00×10^5 bar and T_s=3084 K. The melt mass fraction (x_melt) in this case is 1.81%. For a Mercury-like interior, i.e. with higher Fe content, we find x_env=5.29%, with surface conditions P_s=3.79×10^5 bar and T_s=3290 K. Based on our assumption of the liquidus as the melt curve, we class this as having 0% melt in Table <ref>. In reality, these surface conditions lie between the liquidus and solidus, which would lead to a partially molten surface. This is also the case for the f_silicate=5% interior, with x_env=6.60%, with surface conditions P_s=6.27×10^5 bar and T_s=3461 K. On the other hand, for the extreme case of a pure silicate interior, we find a melt mass fraction of 3.16%, for x_env=2.62%, P_s=1.12×10^5 bar and T_s=2819 K. 


We next consider the higher-temperature P-T profile C3, which permits solutions with a magma ocean surface for all the interior compositions considered. For each interior composition, the permitted envelope mass fraction, and hence the surface pressure, is lower for this hotter P-T profile. For an Earth-like interior, we find x_env=2.52%, with surface conditions P_s=1.36×10^5 bar and T_s=3870 K. The melt mass fraction in this case is 10.16%. For a Mercury-like interior, we find x_env=3.83%, with surface conditions P_s=2.83×10^5 bar and T_s=4200 K. The corresponding x_melt is 5.78%. For the extreme case of a pure silicate interior, we find a lower x_env=1.62%, with P_s=6.97×10^4 bar and T_s=3512 K. The melt mass fraction in this case is larger, at 11.91%. For the other extreme of f_silicate=5%, we obtain x_env=5.06%, for P_s=5.05×10^5 bar and T_s=4503 K, with x_melt=2.62%. We note that including modelling of partial melt would somewhat increase the melt mass fraction in all cases.

As shown in Table 1, the envelope mass fractions and surface conditions we find for profiles C1 and C4 are very similar to C2 and C3 respectively. This is despite the differences in envelope temperature structure resulting from differing haze properties; the difference between C2 and C3 is primarily due to the differing T_int.







 §.§ Volatile Abundances at the Interface



At the surface-atmosphere interface, the interactions between the gas phase equilibrium reactions and solubility of the gases in the magma, if any is present, drive the elemental abundances in the atmosphere. We consider the four P-T profiles presented in Table <ref> and assume 50×solar metallicity, using solar abundances by <cit.>. The assumed metallicity is approximately based on the median retrieved CH4 abundance for K2-18 b <cit.>. Across all considered cases, we find that the dominant H-C-N-O-S gas species at the surface are H2, H2O, CH4, NH3, and H2S. The resulting atmospheric elemental abundances from these scenarios are shown in Table <ref>. As expected, the atmosphere is highly reduced, with oxygen fugacities varying between IW-8.8 and IW-4.9 <cit.> among the 12 cases with magma. We note that although our calculations of the oxygen fugacities agree to within 0.35 dex with the self-consistent IVTANTHERMO code at P_s=10^4 bar and T_s=3000 K, as described in Section <ref>, the redox state at higher pressures/temperatures is not well understood. Future work is needed to better understand the redox state of gas dwarfs at these conditions. 

Overall, we find that H2O and molecules containing N and S are the most dominant volatile species in the magma ocean, with high surface pressures strongly favouring the solubility of N2. As such, for a given interior composition, we find that cooler P-T profiles, resulting in higher P_s, act to increase the depletion of nitrogen in the atmosphere - until the temperature is too low to support a molten surface. The dependence of N depletion on P_s is stronger than that on the melt fraction. Therefore, a hotter temperature profile does not necessarily result in higher N depletion. In terms of the internal structure, we find that the interior needs to be more iron-rich than Earth's interior to result in nitrogen depletion larger than ∼2 dex. 





Whilst we find that nitrogen can be depleted under certain conditions, in line with previous works investigating the solubility of nitrogen in reduced interiors <cit.>, we do not reproduce the six orders of magnitude depletion found by S24. Additionally, we also identify sulfur as a potential atmospheric tracer of a magma ocean; however, the depletion is less than that of nitrogen. Finally, we find that the solubility of H_2, CO, CO_2, and CH_4 is less prominent at the considered conditions and does not drive the abundances of these species far from chemical equilibrium expectations without a magma ocean. However, we note that, as further detailed in Appendix <ref>, many molecular species lack solubility data at the extreme conditions considered here. Hence, further work is needed to improve our knowledge of the solubility of prominent volatiles in silicate melts. 


In Figure <ref>, we show the mixing ratios of the major C-H-O-N-S species in the lower atmosphere and the corresponding elemental abundances for a range of oxygen fugacities using the C3 profile and a Mercury-like interior (f_silicate=30 %). This represents the case with the strongest nitrogen depletion, excluding the extreme 5% silicate interior cases, with atmospheric N/H being ∼2.5 dex lower than the assumed metallicity of 50×solar. We also see the onset of sulfur depletion in the atmosphere due to the solubility of S2 at very reducing conditions (∼IW-6 in this case). On the other hand, the carbon abundance remains unchanged, as mentioned above. We also highlight the potential effect of partial melt by doubling the melt mass fraction, shown by the dotted line in Figure <ref>, leading to an approximately linear increase in the depletion of nitrogen. 





 §.§ Atmospheric Chemistry




We now use the elemental abundances obtained above to determine the atmospheric composition above the surface, using equilibrium and non-equilibrium calculations. From across all the models shown in Table <ref>, we focus on two realistic cases, one with and one without melt. For the molten case, we consider the C3 profile with 30% silicate fraction, which gives a significant N depletion. For the case with no melt we consider the C2 profile with 30% silicate which has no N depletion. For each case, we set the atmospheric elemental budget to that obtained in Section <ref> and reported in Table <ref>. As expected from the model set-up, the no-melt scenario results in all elemental abundances being identical to those of a 50×solar metallicity gas. 

Across all cases considered, the primary O, C, N, and S reservoirs are H_2O, CH_4, NH_3 and H_2S over most of the atmosphere, as indicated by the dashed lines in Figure <ref>. This is seen for a pressure range spanning over 10 orders of magnitude and a temperature profile ranging between ∼260-2700 K. 




The mixing ratio profiles obtained for the no-melt case are shown on the left-hand-side of Figure <ref>. In both the equilibrium and disequilibrium cases, the abundance of H_2O in the upper atmosphere is significantly depleted by a cold trap below the ∼1 bar pressure level. While CO and CO_2 are absent from the photosphere in the equilibrium case, they are present in the disequilibrium case, arising from photochemical processes. However, their abundance is significantly hindered by the limited availability of O, with the main carrier H_2O being depleted by the cold trap. The abundance of CO_2 is lower than that of CO throughout the atmosphere.

Compared to the retrieved atmospheric composition of K2-18 b <cit.>, shown as errorbars and arrows in Figure <ref>, the computed CH_4 abundance is consistent with the retrieved constraint. However, there is a substantial difference of ∼8 dex between the computed abundance of CO_2 with the measured value across the observable pressure range. Additionally, the retrieved upper limits for H_2O and CO are consistent with the computed amounts. Lastly, the computed value of NH_3 is higher than, and therefore inconsistent with the retrieved upper limit.

The right-hand-side of Figure <ref> shows the case with a molten surface. This configuration results in very similar abundances for O- and C-carrying molecules as the no-melt case. This includes the significant depletion of H_2O due to a cold trap, the limited production of CO and CO_2, and CO being more abundant than CO_2.  The main difference from the no-melt case is the notable depletion of NH_3, due to N dissolving in the magma. Specifically NH_3 and N_2 are at much lower mixing ratios than in the no-melt case, by ∼2 dex. Compared to constraints from observations, CH_4, H_2O, CO and in this case NH_3 as well are consistent with the retrieved constraints and upper limits. However, the resulting CO_2 abundance is still substantially lower than the observed abundance.

In summary, we find that even for the case with significant melt, corresponding to our hotter P-T profile with a high T_ int, the NH_3 abundance is close to the observed 95% upper limit, while the CO_2 abundance is still significantly discrepant from the observed value and lower than CO. Therefore, we find that the retrieved atmospheric composition of K2-18 b <cit.> is inconsistent with a magma ocean scenario, or more generally with a deep H_2-rich atmosphere with or without melt, for this planet. In principle, the absence of a cold trap could lead to higher H_2O abundance in the troposphere, which in turn could lead to higher CO_2 abundance. However, such a scenario would also give rise to a significant amount of H_2O and CO, which are presently not detected.



 §.§ Sensitivity to Atmospheric Parameters
 


We also explore other values for the three key atmospheric parameters that may influence the observable composition: the metallicity, the eddy diffusion coefficient K_zz, and the internal temperature T_int. We consider the two cases shown in Figure <ref> as the canonical cases corresponding to the two P-T profiles (C2 and C3). Both cases assume a median metallicity of 50×solar and K_zz of 10^6 cm^2s^-1 in the deep atmosphere. It may be argued that a higher metallicity could result in higher CO_2 abundances than the canonical cases and better match the observed abundances. Similarly, a broader range of K_zz may also influence the abundances. Therefore, for each of the two canonical cases, we investigate models with different values for the metallicity and K_zz. We consider metallicities of 100×solar and 300×solar, representing cases with significantly higher metallicities beyond the median retrieved value of ∼50×solar. For K_zz, we explore two end-member scenarios of 10^4 cm^2s^-1 and 10^8 cm^2s^-1. Based on <cit.> and <cit.>, for our canonical cases we considered values of 25 K and 50 K for T_int. We additionally consider the effect of using a P-T profile with a higher T_int of 60 K, as has been considered by <cit.>. As found for our canonical cases, we find that the observed CO and CO_2 abundances remain unexplained by these models with different values of metallicity, K_zz and T_int. These results are discussed in full in Appendix <ref>.





 §.§ Spectral Characteristics


We use the atmospheric compositions computed in Section <ref> to examine the spectral signatures of CH_4, NH_3, CO and CO_2, which have been previously identified as key diagnostics of the presence of a magma surface. Using the VIRA retrieval framework's <cit.> capability of considering non-uniform vertical mixing ratios, we directly use the atmospheric composition profiles computed using the VULCAN <cit.> non-equilibrium code described above and shown in Figure <ref>. We specifically consider the melt case discussed above, to evaluate the spectral implications for the presence of a magma layer. For all cases, we consider parametric grey cloud and Rayleigh-like haze properties corresponding to the median retrieved constraints of <cit.>, to facilitate a qualitative comparison with the observations. Specifically, we set log(a) = 10^7.31, γ = -11.67, P_c = 10^-0.55 and ϕ_c = 0.63.

The resulting spectral contributions and transmission spectrum are shown in Figure <ref>. As can be seen in the top panel, CH_4 has prominent spectral features throughout the 1-5 μm wavelength range, while CO_2 and CO give rise to absorption features at ∼4.3 and ∼4.7 μm respectively. NH_3 shows a spectral feature at ∼3 μm. Due to the depletion of atmospheric nitrogen arising from its dissolution in the magma, the NH_3 spectral contribution is relatively weak and not detected in the present data. Without such a depletion, i.e. with nitrogen at a solar elemental abundance ratio, NH_3 would have prominent spectral features across the wavelength range of comparable strength to CH_4. While CO is more abundant than CO_2 in the observable atmosphere, as described in Section <ref>, the low absolute abundances of both molecules give rise to comparably weak spectral features.

The resulting transmission spectrum provides a reasonable match for the NIRISS observations of K2-18 b at shorter wavelengths due to the strong CH4 features. However, the spectrum does not fit the prominent CO_2 absorption feature seen in the NIRSpec G395H data. Moreover, the spectral contribution of CO is also minimal. Together, the two molecules are present at abundances that do not provide a good fit to the data in the 4-5 μm range. 

Overall, we find that the gas dwarf scenario of a thick H_2-rich atmosphere of K2-18 b in equilibrium with a magma ocean at depth is not consistent with the existing JWST observations. In particular, irrespective of the NH_3 depletion, the models predict a low CO_2 abundance and CO > CO_2 which are inconsistent with the retrieved abundances. Future studies need to investigate if other effects may contribute to the observed composition. For example, similar to that  discussed in <cit.>, in order for the detected abundance of CO2 to be compatible with a deep H_2-rich atmosphere, an unphysically low C/O ratio of ∼ 0.02–0.06, together with a moderate C/H ratio (∼30-50× solar) and vertical quenching may be required. However, such an atmosphere could also lead to significant CO abundances that may not be consistent with the observations, and the deep atmosphere would have more H_2O than H_2.



§ SUMMARY AND DISCUSSION
 

In this study, we report an integrated framework to investigate the plausibility of magma oceans on temperate gas dwarfs, and their potential atmospheric  signatures. Our framework models the various components of a planet, and their interplay. Specifically, it includes atmospheric and internal structure modelling, magma-atmosphere chemical interactions, and equilibrium as well as disequilibrium (photochemistry and vertical mixing) processes in the atmosphere. Considering all these coupled factors, it predicts the observable abundances of molecular species in the atmosphere and the expected spectral features. 

We apply our framework to perform a comparative assessment of previous works, validating our modelling of magma-atmosphere interactions against <cit.> and assessing the model predictions of <cit.>, S24, for a temperate sub-Neptune. Our findings highlight the importance of considering physically plausible models, set up in a holistic framework. In particular, we note that the use of stand-alone magma-atmosphere interaction models, which do not consider the complex interplay of interior and atmospheric factors, can lead to erroneous results. 



 §.§ Summary


Magma oceans are normally expected for rocky planets with high equilibrium temperatures. In the present work, we have tested the limits of this scenario by exploring whether K2-18 b, a habitable-zone sub-Neptune, can host a magma ocean, as previously suggested by S24, and what the observable signatures could be. We summarise our key findings as follows: 


    
  * An integrated framework is essential to obtain physically plausible and self-consistent results for modelling sub-Neptune gas dwarfs. Our framework includes an atmosphere and interior structure model, including phase diagrams and equations of state of appropriate silicates; thermochemical equilibrium calculations for the silicates-atmosphere interface and lower atmosphere; and disequilibrium processes throughout the atmosphere.
    
    
  * The melt fraction admissible in a gas dwarf depends on atmospheric and interior properties, specifically the interior composition and the atmospheric P-T profile. The P-T profile, in turn, depends strongly on the internal temperature T_int, as well as on the presence and properties of clouds/hazes and on the molecular absorbers present in the atmosphere. For a gas dwarf scenario assuming the bulk parameters of K2-18 b, we find that, with an Earth-like interior composition, maximal melt mass fractions of ∼10% are possible, and may increase somewhat if partial melting is considered. 

    
  * A planet's bulk parameters and temperature structure place both upper and lower limits on the envelope mass fraction, assuming a gas dwarf scenario. For the K2-18 b models considered in this work, these limits are ∼1% and 7% of the planet mass, corresponding to a pure silicate and a 95% iron interior, respectively. The envelope mass fraction affects the surface pressure at the rock-atmosphere boundary, which, in turn, affects the potential melt conditions. 
    
    
  * We find using our framework that the current chemical constraints for K2-18 b are inconsistent with a magma ocean scenario or any gas dwarf scenario, contrary to S24. Firstly, the high observed abundance of CO2 along with low H2O is inconsistent with the chemical expectations for the gas dwarf scenario. Secondly, we find CO to be higher than CO2 by over 1 dex which is also inconsistent with the observations. We find this to be the case with or without a magma ocean, and relatively  independent of the uncertainties in magma-atmosphere interactions at the extremely reduced conditions as described in Appendix <ref>. Finally, we find that N depletion in the atmosphere depends on a wide range of atmospheric and interior parameters, and can range between no depletion and ∼2.5 dex for a realistic model space, given available solubility data. 
    
    
  * Overall, we find that key atmospheric signatures for identifying a gas dwarf include the CO and CO2 abundances, and, if melt is present, possible nitrogen depletion, consistent with some previous studies (cf. Section <ref>). In particular, we expect that CO/CO2> 1 if no H2O is observed (as a result, e.g., of condensation), or, in the presence of H2O, CO/CO2 ≲ 1, due to photolysis of H2O making more oxygen available for the formation of CO2. Furthermore, we find that N depletion is more sensitive to the surface pressure than to the amount of melt present, provided this is non-zero. Thus, the presence of a magma ocean does not ensure a significant N depletion in the atmosphere.  

    
    
  * Our models predict significant H2S for a deep H2-rich atmosphere scenario. Hence, a lack of H2S may be indicative of a shallow atmosphere. However, we note that there are significant uncertainties in the behaviour of S-bearing species in silicate melts at such extremely reducing conditions. Therefore, more robust data for such conditions is needed in order for this signature to be used with a higher degree of confidence. We also note that there is uncertainty in the sulfur photochemical network for such planetary conditions.

    
  * As discussed below, a number of important unknowns remain. In particular, as discussed in Appendix <ref>, the solubility of NH3 in magma remains poorly understood, especially at extremely reducing conditions, as is also the case for H2S at high pressures and temperatures.

   




 §.§ Future Work


In order to aid accurate modelling of potential gas dwarf magma ocean planets, further developments are needed in three areas: (1) solubility laws for volatiles at extremely high pressures and temperatures and very reducing conditions, (2) equations of state (EOS) of silicates at the conditions relevant to temperate sub-Neptunes, and (3) complete reactions lists for all relevant atmospheric species.

As discussed in Appendix <ref>, there is a pressing need for further experimental data and/or ab initio simulations on the solubility of volatile species in silicate melt at the physical and chemical conditions that we have shown in this study to be relevant to the magma-atmosphere interface on sub-Neptunes. This includes high pressure and temperature, and low oxygen fugacity. In particular, the availability of NH3 solubility laws at these conditions would allow more precise prescriptions than assuming its solubility to be negligible, avoiding the resulting likely overestimation of the abundance of N-bearing species in the atmosphere. 
In general, present laws are expected to give an order-of-magnitude estimate of the solubility at the conditions explored in this study; future work is needed to improve the solubility data.

Furthermore, once more accurate and precise solubility laws become available, the non-ideality of gas behaviour at the high pressures relevant at the interface may become a notable source of error if ignored, and will thus need to be appropriately treated <cit.>. We also note that, as a result of the lack of knowledge on the solubility of volatiles in the melt, the phase of the melt itself is not well-constrained. In particular, it is possible that some of the models considered here fall in a regime where there is no surface, and the atmosphere and magma become a single continuous phase at some lower pressure. This would happen if the volatiles were completely miscible in the melt, as is the case for water above a few GPa <cit.>. It is however not known whether this behaviour applies to H2-dominated atmospheres such as the one considered here. Furthermore, even if complete miscibility is not achieved, it is possible that the presence of volatiles in the magma may lead to a change in its EOS, which has not been accounted for here, where we have instead assumed a volatile-free melt for the internal structure calculations.  

There is also scope for future work on the internal structure modelling, including the melt. This includes implementing the partial melting that would occur due to the magma's heterogeneous nature between the solidus and liquidus, as shown in Figures <ref> and <ref>. This is expected to result in a larger fraction of the mantle being at least partially melted than when considering the fully melted region alone, hence further depleting the atmosphere of the most soluble species. This effect is however in part addressed in this work, by considering the impact of a doubled melt mass fraction, as shown in Figure <ref>. Furthermore, future work will include more detailed prescriptions for the mantle, including alternate mineral compositions, and a fully temperature-dependent EOS for the solid portion. 

Overall, JWST provides a promising avenue for atmospheric characterisation of sub-Neptune exoplanets. The high quality of the observations means that concomitant advances need to be made in theoretical models to maximise the scientific return from the data. In this work, we have outlined an end-to-end framework for gas dwarf sub-Neptunes to enable an evaluation of this scenario given high precision JWST data, and highlight the need for more accurate inputs for these models. Such advancements in both observations and theory promise a new era in the characterisation of low-mass exoplanets with JWST.








Acknowledgements: We thank the reviewer for their helpful comments on the manuscript. N.M., L.P.C. and F.R. acknowledge support from the Science & Technologies Facilities Council (STFC) towards the PhD studies of L.P.C. and F.R. (UKRI grants 2886925 and 2605554). N.M. and M.H. acknowledge support from STFC Center for Doctoral Training (CDT) in Data Intensive Science at the University of Cambridge (grant No. ST/P006787/1), and the MERAC Foundation, Switzerland, towards the doctoral studies of M.H. N.M. and S.C. acknowledge support from the  UK Research and Innovation (UKRI) Frontier Grant (EP/X025179/1, PI: N. Madhusudhan). J.D. acknowledges support from Grant EAR‐ 2242946 of National Science Foundation. K.K.M.L acknowledges support from the US Coast Guard Academy Faculty Research Fellowship. J.M. acknowledges NASA grant 80NSSC23K0281. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (<www.csd3.cam.ac.uk>), provided by Dell EMC and Intel using Tier-2 funding from STFC (capital grant EP/P020259/1), and DiRAC funding from STFC (<www.dirac.ac.uk>).




aasjournal






§ AVAILABILITY OF SOLUBILITY LAWS
 
We discuss here the availability of silicate melt solubility laws for the volatile species of interest, at the chemical and physical conditions relevant for magma oceans on temperate sub-Neptunes in the gas dwarf scenario. These findings motivate our choices for the solubility laws adopted in this work. We compile a bibliography of the solubility laws consulted for the preparation of this work in Table <ref>, and show a selection of them in Figure <ref>. For most composition-dependent laws, we adopt a basalt composition for the melt – specifically, when a law explicitly depends on melt composition parameters, we set these to the values corresponding to the Mt Etna basalt from <cit.>. This choice is due to the wide availability of solubility laws for basaltic melt and because of the association of basalt with peridotite, which we assume to be the mantle composition. 



 §.§ Nitrogen Species

The solubility of N2 has been explored for a wide range of parameters (e.g., ), at pressures up to 14.8 GPa and temperatures up to 2800 K <cit.>. By compiling the available data at P ≤ 8.2 GPa and adding their own measurements, <cit.> proposed the solubility law which we use in our calculations.
This law, however, does not appear to extrapolate well at higher pressures and moderately reduced conditions (f_O2∼IW - 2). As warned by <cit.>, experimental data indicate the solubility seems to reach a plateau, while the law predicts solubility to monotonically increase with pressure. A direct comparison with <cit.>'s 10 GPa and 14.8 GPa data points reveals indeed a true solubility ∼ 1 order of magnitude lower than predicted using <cit.>'s law for pure nitrogen vapor. At the extremely reduced conditions explored here, the plateau effect is expected to already be significant at lower pressures <cit.>.

It is also noteworthy that <cit.> - whose data was included in the <cit.> dataset - find some indication of a decrease in the physical solubility of N2 already at P=8 GPa. We note that physical solubility is expected to be the dominant solubility mechanism at the oxidized conditions explored by , as opposed to the chemical solubility relevant at reduced conditions <cit.>.
Nevertheless, as the relevant quantity is not the total pressure, but rather the nitrogen partial pressure, we believe that the <cit.> law can still be a reasonably good approximation even at the reducing, high-pressure conditions that apply at the magma-atmosphere interface, given that the expected N2 mixing ratio in the atmosphere is ≲ 10^-4 in the present models.

The lack of data or simulations for the solubility of NH3 in silicate melt leads us to neglect it, with the caveat that this will lead to our calculations setting only an upper limit on the abundance of N-bearing species in the atmosphere.


 §.§ Carbon Species

Of the three prominent carbon species (CO2, CO, CH4), CO2 is by far the one for which the most complete experimental data on solubility in magma is available (e.g., ). 
Considering this wide dataset, for the case of T = 2273 K and a bulk silicate Earth (BSE) melt composition, <cit.> find that the solubility of CO2 is well-approximated by Henry's Law, which they fit to the data. <cit.>'s law is in excellent agreement with high pressure (P ≥ 8 GPa) molecular dynamics simulations by <cit.> for the T = 2273 K and rhyolite case, and in good agreement with the corresponding mid-ocean ridge basalt (MORB) case. Interestingly, the agreement is slightly worse with the kimberlite melt case, where instead the melt is closest to Suer's BSE composition. 

At lower pressures, the agreement with the simulations is worse, but still within a factor of order unity. In any case, the agreement between <cit.>'s law and <cit.>'s simulations is always satisfactory, which also highlights the weak dependence of the solubility of CO2 on the melt composition, particularly at P ≤ 8 GPa <cit.>. 
Despite the fact that the <cit.> law is intended for lower temperatures than those relevant in this study, due to the lack of more appropriate alternatives we adopt it in our calculations.


For CO, there is a lack of solubility data at high pressure and temperature. Solubility laws are provided by, e.g., <cit.> and <cit.> (for both MORB and rhyolite melts), both of whom carried out experiments at P ∼ 1 GPa and T ∼ 1500 ^∘C. 
The lack of data may be explained by the fact that exploring the solubility of CO at high pressures is especially complicated, because the 2CO = C + CO2 reaction gets skewed to the right as pressure grows, making an initially pure CO vapour spontaneously become mostly CO2 at P ≳ 1 GPa <cit.>. 

An alternative prescription, used by <cit.> as informed by <cit.>, is to instead set the solubility of CO to be one third of that of CO2. This method, taking <cit.>'s BSE law for the CO2 solubility, yields a CO solubility significantly higher than any of the other laws mentioned so far. This might be due to <cit.>'s law being tested for different temperatures and melt compositions. This, however, seems unlikely: on the one hand, the solubility of CO is only weakly dependent on temperature <cit.>; on the other hand, <cit.>'s two laws for MORB and rhyolite yield results less than an order of magnitude apart. This indicates a comparatively weak dependence of solubility on melt composition, while <cit.>'s prescription applied to <cit.>'s law results in a solubility between 1 and 2 orders of magnitude higher, depending on the pressure.
Ultimately, we use the <cit.> MORB law in our calculations, due to it being more recent and calibrated at higher pressures than <cit.>.


The data on CH4 seems to be even sparser: the only solubility law we encountered in the literature – which we use here – is the one in <cit.>, resulting from experiments at 0.7 ≤ P ≤ 3 GPa and 1400 ≤ T ≤ 1450 ^∘C, and the Henry's Law fit to their data by <cit.> - the law by <cit.>, indeed, follows Henry's law for total pressures P ≲ 10^4 bar (at T=3000 K, regardless of the CH4 partial pressure).




 §.§ Other Volatiles

The other major volatiles of note are expected to be sulfur species, water, and H2, as well as He, which, however, as a noble gas, has no impact on the atmospheric chemistry.

Chemical equilibrium calculations indicate that, at conditions relevant at the interface, sulfur will be mostly in H2S, with little S2. 
For S2 we use the law by <cit.>.
It should be noted that this law is calibrated only with data collected at atmospheric pressure, relatively low temperature (T ≤ 1673 K), and not very reducing conditions (ΔIW≥ -1). As such, its extrapolation to the extreme conditions explored in this paper should be considered only as a zeroth-order estimate of the true solubility. No significant high-pressure/high-temperature data to compare with <cit.>'s predictions were found either, the <cit.> high-pressure data being for a carbonate-silicate melt.

For H2S, we found two laws in the literature, by <cit.> – for rhyolite – and by <cit.> – for basaltic melts. The former is calibrated for 1073 ≤ T ≤ 1273 K and P = 2 × 10^3 bar, while the latter for 1323 ≤ T ≤ 1473 K and 250 ≤ P ≤  2 × 10^3 bar. The two laws differ significantly in their temperature dependence: <cit.> find that the solubility of H2S moderately increases with increasing temperature, while the law by <cit.> indicates an extremely strong and negative temperature dependence. Furthermore, <cit.> include a dependence on the mole fraction of FeO in the magma, while the law by <cit.> only depends on thermodynamic parameters. 
However, when extrapolated to high temperature (T ∼ 3000 K) and pressure (P ∼ 10^5 bar), both laws predict negligibly small solubility for H2S at the expected mixing ratios (shown in Figure <ref>). Hence, we do not expect the results of our investigation to be noticeably impacted by the choice of one law over the other. In this investigation, we chose to use the law by <cit.>.

For H2, the law most used in the literature we reviewed is by <cit.>, who carry out experiments at 0.7 ≤ P ≤ 3 GPa and 1400 ≤ T ≤ 1500 ^∘C, and give two expressions, for basaltic and andesitic melt. Their law for basaltic melt is in excellent agreement with that given by <cit.> for BSE melt up to P ∼ 1 GPa, and so is, to a slightly lesser extent, their andesitic melt law. At higher pressures, however, they diverge, with <cit.> predicting Henrian behaviour to arbitrary pressure, while <cit.>'s laws predict a decline in solubility as pressure increases, consistently with their experimental results. 

As the laws given in <cit.> have a robust high-pressure experimental background, we use those, and here, specifically, the basaltic melt case.
 

For H2O, there is a great deal of experimental data on the solubility in silicate melts <cit.>, a complete review of which is beyond the scope of this work. We focus here on two solubility laws: <cit.>, the most recent law available, and <cit.>, which is the one we choose to implement in our study. 
<cit.> provide two slightly different estimates depending on the value of the molar absorption coefficient ϵ_3550, each depending linearly on the square roots of the atmospheric fugacities of both water and molecular hydrogen. These are the result of experiments carried out at very low pressure (P = 1 atm) and high temperature (T=2173 K). 

Higher-pressure experiments are carried out in <cit.>, who also propose a solubility law, calibrated upon a vast but low-temperature experimental database (10^2 ≤ P ≤ 10^4 bar, 1100 ≤ T ≤ 1400 ^∘C), which is in rough agreement with that of <cit.> for an H2-rich envelope. 

The fact that <cit.>'s law depends on a linear combination of the square roots of the fugacities of both H2 and H2O, however, risks breaking element conservation for oxygen: indeed, it would predict some dissolved O in the magma even if no O is present – in any species – in the initial atmospheric composition. This effect is expected to be particularly relevant at the very reduced conditions explored here, where the abundance of O is expected to be low. 

We thus consider extrapolating the law of <cit.> to higher temperatures to be a more accurate prescription than extrapolating that by <cit.> to high pressures, and hence do so here, assuming an Etna basalt composition for the melt. This choice is also consistent with that in <cit.>.  


 §.§ Summary

Data on solubility in silicate melt are available, at some conditions, for several species of interest, with the one exception being NH3, for which we were unable to find any solubility laws. We list the bibliography on solubility laws and/or data points we have explored for this study in Table <ref>, and we show a selection of them in Figure <ref>. In general, the scenario explored in this study, relevant for magma oceans on temperate gas dwarfs, is extreme in a threefold way: it leads to high temperatures (T ≳ 2500 K), high pressures (P ≳ 10^5 bar), and very reduced melts compared to Earth (ΔIW≲ - 5). There is no data, for any species, at such conditions in all three ways. Only for N2 data at both very high temperature and pressure exists, but that is for relatively oxidised conditions <cit.>. High-pressure (P ≥ 10^5 bar) simulations exist for CO2, but only at T ≤ 2273. For S2, for a co-existing fluid phase, high-pressure data only exists at low temperature, and only for carbonate-silicate melt <cit.>. All other species seem to lack high-pressure data.
 
Exploring this region of the parameter space, either experimentally or through simulations, will be crucial for improving our understanding of potential magma oceans in sub-Neptunes, and our ability to lift observational degeneracies with other possible internal structures.







§ SENSITIVITY TO ATMOSPHERIC PARAMETERS
 

As described in Section <ref>, we explore a range of values for three key atmospheric parameters that could influence the observable composition: the metallicity, the eddy diffusion coefficient K_zz, and the internal temperature T_int. Our canonical cases, shown in Figure <ref>, correspond to P-T profiles C2 and C3 with median metallicity of 50×solar, with and without elemental depletion respectively, and K_zz of 10^6 cm^2s^-1 in the deep atmosphere. We investigate if a higher metallicity, a broader range of K_zz values and/or a higher T_int could better match the observed abundances than our canonical cases, for example with higher CO_2 abundance. We therefore consider models with higher metallicities of 100×solar and 300×solar, and two end-member scenarios of 10^4 cm^2s^-1 and 10^8 cm^2s^-1 for K_zz in the deep, convective region. We also consider the effect of using a higher value of T_int of 60 K, as previously considered by <cit.>. Disequilibrium effects due to photochemistry and vertical mixing are included in all cases discussed here.


We start with investigating departures from the canonical C2 case, as shown in Figure <ref>. We first fix the K_zz profile to that used in the canonical case and vary the metallicity as described above. The resulting vertical mixing ratio profiles are shown in Figure <ref>, along with those for 50× metallicity from Figure <ref> for comparison. For both the 100× and 300× solar cases, the abundance of CO_2 remains lower than that of CO throughout the atmosphere, as for the 50× solar case. Similarly, the CO_2 and NH_3 abundances are inconsistent with the retrieved values in the photosphere, between ∼0.01-10 mbar, in all cases. Additionally, the CH_4 abundance for 300× solar metallicity is higher than the retrieved abundance. 

Next we consider a range of K_zz values in the deep atmosphere, using the C2 P-T profile. We vary K_zz at P>0.5 bar from 10^4 to 10^8 cm^2s^-1, with our canonical value at 10^6 cm^2s^-1. The metallicity remains fixed at the canonical value of 50× solar. As shown in Figure <ref>, both the higher and lower K_zz values negligibly affect the computed mixing ratios at observable pressures. Increasing K_zz shifts the quench point to higher (deeper) pressures, as shown by the mixing ratio profile for CO_2 in the right-hand panel of Figure <ref>.

We now consider the hotter P-T profile case with NH_3 depletion due to magma; this is the C3 profile with 30% silicates, as discussed above. The higher metallicities of 100× and 300× solar are implemented by proportionately enhancing the canonical elemental abundances for this case. These originally corresponded to 50× solar, hence we increase the relevant elemental abundances in Table <ref> by factors of 2 and 6, respectively. The results are shown in Figure <ref>. As for the C2 profile, the predicted CO_2 abundance remains significantly below the retrieved value in both cases, with the CO mixing ratio exceeding that of CO_2 throughout the atmosphere. The CH_4 abundance for 300× solar metallicity is additionally too high compared to the retrieved abundance.

As an end-member case, we consider each of the C2 and C3 profiles discussed above and adopt our extreme values of 300×solar metallicity and K_zz=10^8 cm^2s^-1 in the deep atmosphere. The resulting vertical mixing ratio profiles are shown in Figure <ref> along with the canonical cases. These end-member cases are similarly unable to match the retrieved CO_2 abundance constraints. A higher K_zz would further increase the abundances of both CO and CO_2, however CO remains more abundant than CO_2.

Thus far we have considered values of 25 K and 50 K for T_int, corresponding the C2 and C3 profiles, respectively. Lastly, we explore the effect of increasing T_int to a higher value of 60 K for completeness, as has been considered by other works for K2-18 b <cit.>. We adopt the P-T profile of <cit.> with 100× solar metallicity, extrapolated to higher pressures (1000 bar) using an adiabat. We consider two cases: 1) 100× solar metallicity with depletion (i.e. twice the C3 30% silicates abundances from Table <ref>) and our canonical K_zz treatment, and 2) a high K_zz=10^8 cm^2s^-1 and a high metallicity of 300× solar (i.e. 6× the C3 30% silicates abundances). With the canonical K_zz, we find that the computed CO abundance exceeds the retrieved upper limit, while the computed CO_2 abundance remains significantly lower than the retrieved abundance. The retrieved CH_4 abundance and NH_3 upper limits can be explained by this model. For the high K_zz and high metallicity case, the computed CO abundance similarly exceeds the retrieved abundance. In this case the retrieved CO_2 abundance can be explained by the model. However, the computed CH_4 abundance exceeds the retrieved value. Due to the higher temperatures in this P-T profile, the H_2O abundance exceeds the retrieved value for both cases of metallicity and K_zz considered.

Overall, we have explored a wide parameter space for the atmospheric chemistry, considering a range of values for K_zz, metallicity and T_int. In this exploration, we do not find a case resulting in CO_2 > CO that would satisfy the retrieved atmospheric abundance constraints for K2-18 b <cit.>.













