The recent JWST detections of carbonbearing molecules in a habitablezone subNeptune have opened a new era in the study of lowmass exoplanets. The subNeptune regime spans a wide diversity of planetary interiors and atmospheres not witnessed in the solar system including miniNeptunes superEarths and water worlds. Recent works have investigated the possibility of gas dwarfs with rocky interiors and thick HSUB2SUBrich atmospheres to explain aspects of the subNeptune population including the radius valley. Interactions between the HSUB2SUBrich envelope and a potential magma ocean may lead to observable atmospheric signatures. We report a coupled interioratmosphere modeling framework for gas dwarfs to investigate the plausibility of magma oceans on such planets and their observable diagnostics. We find that the surfaceatmosphere 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 hightemperature rocky planets we assess if such conditions may be admissible and observable for temperate subNeptunes. We find that a holistic modeling approach is required for this purpose and to avoid unphysical model solutions. Using our model framework we consider the habitablezone subNeptune K218 b as a case study and find that its observed atmospheric composition is incompatible with a magma ocean scenario. We identify key atmospheric molecular and elemental diagnostics including the abundances of COSUB2SUB CO NHSUB3SUB and potentially Sbearing species. Our study also underscores the need for fundamental material properties for accurate modeling of such planets.
['Exoplanets', 'Exoplanet atmospheres', 'Planetary interior', '498', '487', '1248', 'Astrophysics - Earth and Planetary Astrophysics']
Toward a Selfconsistent Evaluation of Gas Dwarf Scenarios for Temperate SubNeptunes
{'Towards a self-consistent evaluation of gas dwarf scenarios for temperate sub-Neptunes': 'Frances E. Rigby, 1, ∗ Lorenzo Pica-Ciamarra, 1, ∗ M˚ans Holmberg, 1, ∗ Nikku Madhusudhan, 1, ∗ Savvas Constantinou, 1, ∗ Laura Schaefer, 2 Jie Deng, 3 Kanani K. M. Lee, 4 and Julianne I. Moses 5  \n1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK 2 Department of Geological Sciences, School of Earth, Energy, and Environmental Sciences, Stanford University, Stanford, CA 94305, USA  \n3 Department of Geosciences, Princeton University, Princeton, NJ 08544, USA  \n- 4 Department of Physics, United States Coast Guard Academy, New London, CT 06320, USA  \n5 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.', '1. 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., Fressin et al. 2013; Fulton & Petigura 2018). The nature of the sub-Neptune population remains debated, as their bulk densities can be explained by a number of degenerate interior compositions (e.g. Rogers et al. 2011; Valencia et al. 2013). These include rocky planets with diverse atmospheric compositions, mini-Neptunes with volatile-rich interiors and deep H 2 -rich atmospheres, and water worlds with sub-  \nCorresponding author: Nikku Madhusudhan  \nnmadhu@ast.cam.ac.uk  \n- ∗ These authors contributed comparably to this work.  \nstantial water mass fractions, including Hycean worlds (e.g., Rogers et al. 2011; Valencia et al. 2013; Dorn et al. 2017; Zeng et al. 2019; Madhusudhan et al. 2020, 2021; Rigby & Madhusudhan 2024).  \nThe 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 (Madhusudhan et al. 2021) K2-18 b (Madhusudhan et al. 2023b), 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 (Holmberg & Madhusudhan 2024; Benneke et al. 2024) where abundance constraints for CH 4 and CO 2 were also  \nretrieved - and LHS 1140 b (Doyon 2024; Damiano et al. 2024). 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.  \nOne 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 ⊕ (Fulton et al. 2017; Fulton & Petigura 2018; Cloutier & Menou 2020). 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 H 2 -rich atmospheres, as gas dwarfs. The mechanism for the mass loss is debated, with the predictions of two hypotheses - photoevaporation (e.g., Lopez & Fortney 2013; Jin et al. 2014; Owen & Wu 2017; Jin & Mordasini 2018) and core-powered mass loss (e.g., Ginzburg et al. 2016, 2018; Gupta & Schlichting 2019, 2020) - both proposed to explain the observations (Rogers et al. 2021). A second explanation (e.g., Zeng et al. 2019; Venturini et al. 2020; Izidoro et al. 2021) 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 (e.g., Kite et al. 2019, 2020; Daviau & Lee 2021; Gaillard et al. 2022; Schlichting & Young 2022; Charnoz et al. 2023; Misener et al. 2023; Falco et al. 2024).  \nWhile the gas dwarf hypothesis has garnered significant attention in the literature (e.g., Lopez & Fortney 2013; Jin et al. 2014; Ginzburg et al. 2016; Owen & Wu 2017; Jin & Mordasini 2018; Ginzburg et al. 2018; Gupta & Schlichting 2019; Kite et al. 2019; Gupta & Schlichting 2020; Kite et al. 2020; Bean et al. 2021; Schlichting & Young 2022; Charnoz et al. 2023), 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 (Fortney et al. 2013; Venturini et al. 2024). Secondly,  \nshould 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 (Schaefer et al. 2016; Schaefer & Fegley 2017; Kite et al. 2019, 2020; Daviau & Lee 2021; Gaillard et al. 2022; Schlichting & Young 2022; Misener et al. 2023; Charnoz et al. 2023; Falco et al. 2024; Shorttle et al. 2024; Tian & Heng 2024). 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.  \nSeveral recent studies have explored the implications of a magma ocean on the atmosphere and interior compositions of diverse planets, both with terrestrial-like (Schaefer et al. 2016; Schaefer & Fegley 2017; Daviau & Lee 2021; Gaillard et al. 2022; Tian & Heng 2024) and H 2 -rich atmospheres (Kite et al. 2019, 2020; Schlichting & Young 2022; Misener et al. 2023; Charnoz et al. 2023; Falco et al. 2024; Shorttle et al. 2024; Tian & Heng 2024). 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 (e.g., Daviau & Lee 2021; Dasgupta et al. 2022; Suer et al. 2023; Shorttle et al. 2024) and high CO/CO 2 ratio for H 2 -rich atmospheres (Gaillard et al. 2022; Schlichting & Young 2022). 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 (e.g., Kite et al. 2020; Schlichting & Young 2022).  \nIn this work, we develop an integrated magma ocean framework for temperate, H 2 -rich sub-Neptunes. Our framework, presented in Section 2, 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 habitablezone sub-Neptune and Hycean candidate K2-18 b (Madhusudhan et al. 2020) to investigate the plausibility of a magma ocean (e.g. Kite et al. 2020; Shorttle et al. 2024) and, if present, its atmospheric signatures. In doing  \nso, we first use our framework to perform a comparative assessment of previous works in this direction in Section 3, both on terrestrial-like atmospheres (Gaillard et al. 2022) and on H 2 -rich ones (Kite et al. 2019, 2020; Schlichting & Young 2022; Charnoz et al. 2023; Misener et al. 2023; Falco et al. 2024; Shorttle et al. 2024; Tian & Heng 2024), with a focus on the case study of the candidate Hycean world K2-18 b. We then present our model predictions in Section 4. Finally, we summarize our findings and discuss future work in Section 5, highlighting the need for physically consistent models, and new experimental and theoretical work to derive accurate fundamental material properties.", '2. METHODS': "We develop an integrated modelling framework to evaluate gas dwarf scenarios for planets in the subNeptune regime. A schematic flowchart of the framework is shown in Figure 1. 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.", '2.1. 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,  \nFigure 1. Flowchart showing an integrated modelling framework to assess gas dwarf scenarios for planets in the sub-Neptune regime.  \n<!-- image -->  \nthe stellar properties, the presence and characteristics of clouds/hazes in the planet's atmosphere, and the efficiency of day-night energy redistribution. The selfconsistent calculation will yield a pressure-temperature ( P -T ) profile, which will be coupled to the internal structure model, as discussed in Section 2.2.  \nIn order to carry out the self-consistent modelling of the atmospheric structure, we use the GENESIS framework (Gandhi & Madhusudhan 2017) adapted for sub-Neptunes (Madhusudhan et al. 2020; Piette & Madhusudhan 2020; Madhusudhan et al. 2021, 2023a). 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 (Hubeny 2017) and the discontinuous finite element method (Castor et al. 1992), while taking  \ninto account all of the parameters mentioned earlier in this section.  \nFor the atmospheric composition, we adopt uniform mixing ratios of molecular species based on the retrieved values at the terminator region of K2-18 b (Madhusudhan et al. 2023b). We use the median retrieved abundances for the one-offset case: log X CH 4 = -1 . 72 and log X CO 2 = -2 . 04. For H 2 O, we consider the 95% oneoffset upper limit, log X H 2 O = -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 Madhusudhan et al. (2020) and Valencia et al. (2013), 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 5.  \nWe 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 H 2 /He equation of state (EOS), ρ = ρ ( P, T ), and adiabatic gradient from Chabrier et al. (2019). 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 (Madhusudhan et al. 2023b).", '2.2. Internal Structure Modelling': "We model planetary internal structures using the HyRIS framework, outlined in Rigby & Madhusudhan (2024). 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.  \nThe 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 Monteux et al. (2016) by combin-  \nFigure 2. Cross-section of the internal structure of a potential gas dwarf, including the H 2 -rich envelope, silicate mantle and iron core.  \n<!-- image -->  \ning the densities of molten enstatite, forsterite, fayalite, anorthite and diopside, described by third-order BirchMurnaghan/Mie-Gruneisen EOSs from Thomas & Asimow (2013), 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 Monteux et al. (2016), based on Fiquet et al. (2010) - both the liquidus and solidus are shown in Figure 3 (Fiquet et al. 2010; Monteux et al. 2016). The solid portion of the silicate mantle is described by the EOS of Lee et al. (2004) 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 Seager et al. (2007) for MgSiO 3 perovskite, originally derived at room temperature. The thermal effects for solid silicates at these pressures are small (Seager et al. 2007) with negligible effect on the internal structure. The iron core is described by the EOS of Seager et al. (2007) for hexagonal close-packed Fe.  \nThe temperature structure in the melt is assumed to be adiabatic. The adiabatic gradient is calculated using the specific heat for peridotite from Monteux et al. (2016) 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 Lee et al. (2004). Following previous studies (e.g., Rogers et al. 2011; Nixon & Madhusudhan 2021; Rigby & Madhusudhan 2024), the remaining solid  \nportion of the interior is taken to be isothermal, as the EOSs used are temperature-independent (Seager et al. 2007).  \nThe mass of the magma ocean follows from the adiabatic temperature profile in the melt, similarly to the calculation of water ocean depths by Nixon & Madhusudhan (2021) and Rigby & Madhusudhan (2024). 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 Rigby & Madhusudhan (2024). 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 4.3.", '2.3. 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 X , 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 (Hirschmann et al. 2012, basalt case), H 2 O (IaconoMarziano et al. 2012), CO (Yoshioka et al. 2019, MORB case), CO 2 (Suer et al. 2023), CH 4 (Ardia et al. 2013), N 2 (Dasgupta et al. 2022), S 2 (Gaillard et al. 2022) and H 2 S (Clemente et al. 2004), as further motivated in Appendix A. We note that the solubility of H 2 S is uncertain at high temperatures/pressures and may be higher if, for example, its solubility approaches that of S 2 . 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 N 2 , and may be higher if the solubility of NH 3 is significant. The data on NH 3 solubility in magma is currently limited and it is difficult to make any quantitative estimates of NH 3 solubility. For the explicitly composition-dependent laws, we use the Iacono-Marziano et al. (2012) Etna basalt melt composition. Similarly to Kite et al. (2020), we assume that the magma is well-stirred such that the  \nequilibration at the surface sets the volatile abundance throughout the melt.  \nThese 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 (similar to Gaillard et al. 2022),  \nM tot w i = M atm w i, atm + M melt w i, melt , (1)  \nwhere w i is the total mass fraction of each species i .  \nTo determine the chemical composition of the atmosphere and the melt, we solve the element conservation equations  \nε j = ∑ i ∈ X ν ij n i n ⟨ H ⟩ , (2)  \nwhere n i is the total amount of moles of species i , n ⟨ H ⟩ = n H + 2 n H 2 + 2 n H 2 O + . . . 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 (2) is coupled to Equation (1) via n i ∝ w i /µ i , where µ i is the molar mass, which in turn is coupled to the law of mass action  \np i p · -= K i ∏ j ∈ E ( p j p · -) ν ij , (3)  \nand 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  \nln K ( T ) = a 0 T + a 1 ln T + b 0 + b 1 T + b 2 T 2 , (4)  \nusing the coefficients provided by FASTCHEM (Stock et al. 2018, 2022), mainly derived using thermochemical data from Chase (1998).  \nOverall, Equation (3) 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 (2) for all elements apart from H, together with  \nP s = ∑ i ∈ X p i , (5)  \nto 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 O 2 ) as a free parameter, by determining p O in Equation (3) via f O 2 = p O 2 = K O 2 p 2 O /p · -, allowing us to  \nFigure 3. Left: Pressure in the H 2 -rich envelope against the envelope mass fraction above this pressure level. Solid and dashed lines show these assuming C2 and the profile used by S24 in the envelope respectively. The black squares indicate expected surface conditions for different envelope mass fractions, independent of satisfying the bulk properties of the planet. Right: The nominal pressure-temperature profiles generated for this work and S24 shown against the liquidus and solidus for peridotite (Fiquet et al. 2010; Monteux et al. 2016). The solid and liquid phases are shaded. The black squares again show the corresponding surface conditions expected for the different envelope mass fractions.  \n<!-- image -->  \nconsider different redox conditions. In this framework, we assume ideal gas behavior such that fugacity and partial pressure are equivalent (e.g. Bower et al. 2022; Schlichting & Young 2022).  \nAs a cross-check, we validate our new framework against a self-consistent atmosphere composition model (e.g. Schaefer & Fegley 2017) which uses the Gibbs energy minimization code IVTANTHERMO (Belov et al. 1999). IVTANTHERMO uses a thermodynamic database based on Gurvich & Veyts (1990), 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.  \nFor this comparison, we use 50 × solar bulk elemental abundances (not including He), P s = 10 4 bar,  \nT 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 CO 2 . 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 (Stock et al. 2022) and GGCHEM (Woitke et al. 2018) when setting M melt = 0.", '2.4. 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 (Tsai et al. 2021), with the initial atmospheric chemistry obtained using the FASTCHEM equilibrium chemistry code (Stock et al. 2018).  \nFor equilibrium chemistry calculations, we consider thermochemical equilibrium involving H-C-N-OS species as well as He, along with H 2 O condensation. For calculations considering disequilibrium processes, we additionally include the effects of vertical mixing and photochemistry. We follow the K zz parameterisation of  \nMadhusudhan et al. (2023a):  \nK 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 , (6)  \nalthough 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 Leconte et al. 2024). Accordingly, we consider a wider range of K zz values than our canonical treatment in Section 4.5 and Appendix B.  \nAdditionally, we consider photochemical reactions including H-C-N-O-S species, using a nominal stellar spectrum from the HAZMAT spectral library (Peacock et al. 2020) corresponding to a median 5 Gyr star of radius 0.45 R ⊙ following previous work (Madhusudhan et al. 2023a). We also specifically consider the condensation of H 2 O to liquid and solid droplets, which fall at their terminal velocity, as described in Tsai et al. (2021). We note that while the H-C-N-O chemistry has been extensively explored for sub-Neptunes in various studies (e.g. Yu et al. 2021; Hu et al. 2021; Tsai et al. 2021; Madhusudhan et al. 2023a), the S chemistry has not been explored in significant detail and may be incomplete. Nevertheless, we include S using the VULCAN framework (Tsai et al. 2021) for completeness.  \nWith 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.', '2.5. Spectral Characteristics': "We use the results of the chemistry calculation described in Section 2.4 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 (Constantinou & Madhusudhan 2024), which treats the planet's terminator as a 1D atmosphere in hydrostatic equilibrium. We consider atmospheric opacity contributions from H 2 O (Barber et al. 2006; Rothman et al. 2010), CH 4 (Yurchenko & Tennyson 2014), NH 3 (Yurchenko et al. 2011), CO (Li et al. 2015), CO 2 (Tashkun et al. 2015), C 2 H 2 (Chubb et al. 2020), HCN (Barber et al. 2014), H 2 S (Azzam et al. 2016; Chubb et al. 2018) and SO 2 (Underwood et al. 2016). We do not include N 2 in the model, as it has no significant absorption features  \nin 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 (Borysow et al. 1988; Orton et al. 2007; Abel et al. 2011; Richard et al. 2012), 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.", '3. RESULTS: COMPARISON WITH PREVIOUS WORK': 'We now apply the framework described in Section 2 and compare with previous works on both terrestrial-like and sub-Neptune atmospheres.', '3.1. Terrestrial-like Atmospheres': "Many previous studies have investigated surfaceatmosphere interactions for magma oceans underneath terrestrial-like atmospheres (e.g., Matsui & Abe 1986; Elkins-Tanton 2008; Hamano et al. 2013; Lebrun et al. 2013; Wordsworth 2016; Kite & Schaefer 2021; Lichtenberg et al. 2021; Bower et al. 2022; Gaillard et al. 2022). Recent studies have explored the implications of diverse interiors of exoplanets for their atmospheric compositions. Daviau & Lee (2021) 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, Gaillard et al. (2022) 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-wustite (IW) buffer, f O 2 ≲ IW -2, result in an atmosphere abundant in H 2 , CO and CH 4 but depleted in CO 2 and N 2 . On the other hand, for f O 2 ≳ IW+2, CO 2 becomes the main atmospheric component, with significant levels of SO 2 , N 2 and H 2 O. In particular, the behaviour of nitrogen is a consequence of the high solubility of N 2 as N 3in silicate melt at reducing conditions (e.g., Libourel et al. 2003; Dasgupta et al. 2022), via the following reaction  \n1 2 N 2 (gas) + 3 2 O 2(melt) ⇌ N 3(melt) + 3 4 O 2 (gas) . (7)  \nAs a result, the melt concentration of N 3is proportional to f 1 / 2 N 2 f -3 / 4 O 2 , thus favouring low f O 2 . 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.  \nAs a benchmark, we compare our melt-atmosphere equilibrium chemistry framework with Gaillard et al. (2022). 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 4. Compared to our nominal setup in Section 2.3, we added a constraint for the hydrogen abundance and solved for the resulting mass of the atmosphere, coupled to the surface pressure using (Gaillard et al. 2022)  \nP s = gM atm 4 πR 2 p , (8)  \nwhere 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 Gaillard et al. (2022). Overall, we find good agreement between both implementations, with the most deviation coming from N 2 and CH 4 . We find that the N 2 discrepancy comes from an inconsistency in the code by Gaillard et al. (2022), whereby they use a molar mass of 14 g/mol for N 2 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 Gaillard et al. (2022), as shown in Figure 4. For this purpose, in addition to considering their adopted molar mass, we implemented the reactions CH 4 + 2O 2 ⇌ 2H 2 O + CO 2 and H 2 O ⇌ 0 . 5O 2 +H 2 using the equilibrium constants by Gaillard et al. (2022) to obtain the partial pressures of CH 4 and H 2 O, instead of deriving these from the elemental partial pressures as described in Section 2.3. We also used the oxygen fugacity of the IW buffer from Gaillard et al. (2022) instead of Hirschmann (2021).", '3.2. Sub-Neptunes with H 2 -rich Atmospheres': "Several recent studies have also explored magmaatmosphere interactions in sub-Neptunes with rocky interiors and H 2 -rich atmospheres. Kite et al. (2019) considered the impact of H 2 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 ≳ 3R ⊕ . They find that the high solubility of H 2 in magma, especially at high pressure, limits the maximum radius that can be attained by sub-Neptunes through accretion of atmospheric H 2 . For a 10 M ⊕ core they find a limiting mass fraction of 1 . 5 wt% H 2 in the atmosphere - corresponding to > 20 wt% H 2 in the planet -, as any additional H 2 would be stored almost exclusively in the interior. Looking at smaller planets  \nFigure 4. Comparison with Gaillard et al. (2022), showing the atmospheric composition as a function of oxygen fugacity. The solid lines show the partial pressures from our melt-atmosphere framework, described in Section 2.3. The transparent lines are computed using the code from Gaillard et al. (2022), and the dashed lines are computed with our framework modified to approximate the results from Gaillard et al. (2022). The dotted line at the top illustrates the total surface pressure.  \n<!-- image -->  \n(2 R ⊕ ≤ R p ≤ 3 R ⊕ ), Kite et al. (2020) find that magma-atmosphere interactions would significantly affect the atmosphere's composition and mass. For example, a key insight is that the H 2 O/H 2 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 H 2 O/H 2 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.  \nFurther investigations were carried out by Schlichting & Young (2022), Charnoz et al. (2023) and, most recently, Falco et al. (2024). Considering a surface temperature T s = 4500 K, 1% to 14% H mass fractions (of overall planet mass) and model parameters resulting in f O 2 ≲ IW -2, Schlichting & Young (2022) find that the atmosphere is expected to be dominated across the explored parameter space by H 2 , SiO, CO, Mg and Na, followed by H 2 O, which should exceed CO 2 and CH 4 by two to three orders of magnitude. It should be noted they do not include N in their model. Charnoz et al. (2023) and Falco et al. (2024), 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 H 2 O and SiO should be expected. Additionally, the volatile-free investigation by Misener et al. (2023) finds that silane  \n(SiH 4 ) should also be expected, dominating over SiO at P ≳ 0 . 1 bar for an isothermal T = 1000 K pressuretemperature ( P -T ) profile in the upper atmosphere. Most recently, Tian & Heng (2024) also investigated the outgassing mechanism for hybrid atmospheres in sub-Neptunes, but without considering solubilities in magma.", '3.3. 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 Shorttle et al. (2024), hereafter S24. Similarly to Schlichting & Young (2022) and Gaillard et al. (2022), S24 point to a high CO/CO 2 ratio and, like Daviau & Lee (2021) and Gaillard et al. (2022), 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 (e.g., Kite et al. 2016; Schaefer et al. 2016; Kite et al. 2020; Gaillard et al. 2022; Charnoz et al. 2023; Misener et al. 2023; Falco et al. 2024), 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 2 and Figure 1.  \nWe 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 (e.g. Hu et al. 2021; Wogan et al. 2024). Wogan et al. (2024) conduct photochemical modelling of miniNeptune cases for K2-18 b, suggesting a plausible solution. However, as noted in Glein (2024), the calculated abundances are unable to match the retrieved abundances (Madhusudhan et al. 2023b). 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.', '3.3.1. 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 (Madhusudhan et al. 2020, 2021; Rigby & Madhusudhan 2024). Considering the present gas dwarf scenario, a purely rocky interior would require a minimum H 2 -rich envelope mass fraction of ∼ 1% (Madhusudhan et al. 2020), as discussed below.  \nThe model grid of S24 contains four values of mantle mass fraction relative to the total planet mass (0.001,  \n0.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.  \nIn 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 1: 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 ⊕ (Cloutier et al. 2019) and radius R p = 2 . 61 R ⊕ (Benneke et al. 2019) 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 1.  \nConsidering the four self-consistent P -T profiles described in Section 2.1, 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 Benneke et al. (2019) 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 4.", '3.3.2. Feasibility of a magma ocean': 'As described in Section 2.2 and shown in Figure 3, given an interior composition, the choice of P -T profile  \naffects 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 Benneke et al. (2019) at low pressures ( P ≤ 4 bar) and perform a loglinear 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 (e.g., Hu 2021; Madhusudhan et al. 2023a; Leconte et al. 2024); an example is shown in Figure 3.  \nWe 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 selfconsistently 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 1 and Figure 3.  \nNevertheless, 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 H 2 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 3 along with the liquidus and solidus curves for peridotite (Fiquet et al. 2010; Monteux et al. 2016). 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 5, 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%.  \nIf 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 2.3, 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.  \nFor 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, CO 2 +2H 2 ⇌ CH 4 +O 2 , 2CO 2 ⇌ 2CO+O 2 , and 2H 2 O ⇌ 2H 2 +O 2 , in thermochemical equilibrium, and solubilities of CH 4 , N 2 , CO 2 , and H 2 O 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 9 in Appendix A.  \nIn our framework, described in Section 2.3, 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 1.', '3.3.4. 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 NH 3 log-mixing ratio to at most log X NH 3 ∼ -4, close to the upper bound of -4 . 46 found by Madhusudhan et al. (2023b), 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  \nTable 1. H/He envelope mass fraction, resulting surface temperature, pressure, melt mass fraction constrained by the median values of the K2-18 b bulk parameters ( R p = 2 . 61 R ⊕ , M p = 8 . 63 M ⊕ ), and atmospheric elemental abundances. We use four interior compositions, where 5% and 100% silicate are unrealistic extreme cases included for completeness, and the four self-consistent P -T profiles generated for this work. We note that the cases with 0 % may include a region of partial melt. The bulk elemental abundances in the melt and atmosphere combined are set to 50 × solar.  \nN may be difficult to reconcile with potential formation mechanisms.  \nWe 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 Asplund et al. (2009) value for (C/H) ⊙ , i.e., 3 . 2 × 10 -3 by mass. Secondly, as argued in Section 3.3.2, 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 logmixing ratio for CO 2 of log X CO 2 ≲ -3. This is at the lowest end, if not outside, of the 1 σ confidence interval presented in Madhusudhan et al. (2023b). Furthermore, the CO abundance or the CO 2 /CO ratio are not reported in S24, making it difficult to assess the validity of the chemical estimates.  \nFinally, 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 at-  \nmospheric spectral model and to obtain a statistically robust fit to the data; that is the purpose of atmospheric retrievals (Madhusudhan 2018). 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 (Madhusudhan et al. 2023b). The CO and H 2 O abundances are not reported in S24, which prevents a clear assessment of the agreement between the chemical predictions and the retrieved abundances.', '4. RESULTS: A CASE STUDY OF K2-18 b': 'After having established the consistency of our results with Gaillard et al. (2022), 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 1, 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 subNeptune magma world.  \nFigure 5. Atmospheric pressure-temperature profiles shown against the adopted phase boundary for the silicate mantle, with the dashed lines corresponding to the Monteux et al. (2016) liquidus and solidus for peridotite. The black circles indicate the surface conditions for the cases discussed and the coloured lines show the adiabatic temperature structure in the melt, adopting the liquidus as the melt-solid transition.  \n<!-- image -->', '4.1. Atmospheric Structure': 'As discussed in Section 2.1, the dayside atmospheric structure is calculated self-consistently from the atmospheric constraints retrieved in the one-offset case of Madhusudhan et al. (2023b): the median log X CH 4 = -1 . 74, log X CO 2 = -2 . 04, and the 2 σ upper bound log X H 2 O = -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 H 2 -rich sub-Neptune atmospheres has been carried out before, in Piette & Madhusudhan (2020). 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 (Madhusudhan et al. 2023b). All the P -T profiles are shown in Figure 5.', '4.2. 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  \nof the planet, as discussed in Sections 2.2 and 3.3.1 and shown in Table 1. We vary the interior composition from f silicate = 5% to f silicate = 100%, adopting the median M p = 8 . 63 M ⊕ (Cloutier et al. 2019) and R p = 2 . 61 R ⊕ (Benneke et al. 2019). 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 Madhusudhan et al. (2020).  \nIn Figure 5 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 1.  \nThe 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 1. 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.  \nWe 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  \nwould somewhat increase the melt mass fraction in all cases.  \nAs 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 .', '4.3. 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 1 and assume 50 × solar metallicity, using solar abundances by Asplund et al. (2021). The assumed metallicity is approximately based on the median retrieved CH 4 abundance for K218 b (Madhusudhan et al. 2023b). Across all considered cases, we find that the dominant H-C-N-O-S gas species at the surface are H 2 , H 2 O, CH 4 , NH 3 , and H 2 S. The resulting atmospheric elemental abundances from these scenarios are shown in Table 1. As expected, the atmosphere is highly reduced, with oxygen fugacities varying between IW-8.8 and IW-4.9 (using the oxygen fugacity of the IW buffer by Hirschmann 2021) 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 2.3, 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.  \nOverall, we find that H 2 O 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 N 2 . 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.  \nWhilst we find that nitrogen can be depleted under certain conditions, in line with previous works investigating the solubility of nitrogen in reduced interiors (Daviau & Lee 2021; Dasgupta et al. 2022; Suer et al. 2023; Shorttle et al. 2024), we do not reproduce the six  \norders 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 A, 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.  \nIn Figure 6, 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 S 2 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 6, leading to an approximately linear increase in the depletion of nitrogen.", '4.4. 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 1, 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 4.3 and reported in Table 1. 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.  \nAcross all cases considered, the primary O, C, N, and S reservoirs are H 2 O, CH 4 , NH 3 and H 2 S over most of the atmosphere, as indicated by the dashed lines in Figure 7. This is seen for a pressure range spanning over 10 orders of magnitude and a temperature profile ranging between ∼ 260-2700 K.  \n<!-- image -->  \nFigure 6. Atmospheric composition at the melt-atmosphere interface as a function of oxygen fugacity, at T s = 4200 K, P s = 2 . 83 × 10 5 bar, M melt = 0 . 0578 M p , M atm = 0 . 0383 M p and 50 × solar metallicity - corresponding to the C3 profile with a Mercury-like interior composition ( f silicate = 30%) in Table 1. Left: Atmospheric mixing ratios of major H-C-O-N-S species. The solid and dashed lines show the abundances with and without solubility, respectively. Right: Atmospheric elemental abundances normalised to hydrogen. The dotted line corresponds to a case with twice the melt mass fraction to highlight the potential effect of partial melt. The dashed horizontal line shows 50 × solar. The grey region centred at IW-7 ( ± 0 . 5), shown in both figures, corresponds to the approximate oxygen fugacity obtained with a total oxygen budget of 50 × solar.  \n<!-- image -->  \nFigure 7. Vertical mixing ratio profiles for several H-C-O-N-S molecular species. X i denotes the volume mixing ratio of a species i . Solid and dashed lines denote profiles computed with and without disequilibrium effects. Horizontal data points and arrows denote the mixing ratio constraints and 95% upper estimates retrieved by Madhusudhan et al. (2023b). The gray shaded region denotes the pressure range typically probed by transmission spectroscopy (e.g. Constantinou & Madhusudhan 2024). Left: Mixing ratio profiles corresponding to the C2 P -T profile and 30% interior silicate fraction case shown in Table 1. This corresponds to a 50 × solar elemental abundances. Right: Profiles computed for the C3 P -T profile and 30% interior silicate fraction case shown in Table 1. N is depleted due to dissolution in the magma surface.  \n<!-- image -->  \nFigure 8. Top: Spectral contributions arising from CH 4 , NH 3 , CO and CO 2 in the transmission spectrum of K2-18 b. The atmospheric abundance of each molecule corresponds to the right-hand side plot of Figure 7, generated with the C3 P -T profile, and includes disequilibrium effects. The dashed grey line corresponds to the spectral contribution of NH 3 if it were not depleted by dissolution in the magma. Bottom: The resulting transmission spectrum from all four species' spectral contributions. Blue and orange errorbars are JWST NIRISS and NIRSpec G395H observations of K2-18 b, which include the -41 ppm retrieved offset reported by Madhusudhan et al. (2023b). It can be seen that the magma ocean scenario does not result in sufficient CO 2 to explain the observations at ∼ 4.3 µ m. We emphasise that the present comparison to the data is solely for illustration. A robust comparison necessitates considering the constraints obtained from a detailed retrieval analysis, as done in Figure 7.  \n<!-- image -->  \nThe mixing ratio profiles obtained for the no-melt case are shown on the left-hand-side of Figure 7. In both the equilibrium and disequilibrium cases, the abundance of H 2 O 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 2 O being depleted by the cold trap. The abundance of CO 2 is lower than that of CO throughout the atmosphere.  \nCompared to the retrieved atmospheric composition of K2-18 b (Madhusudhan et al. 2023b), shown as errorbars and arrows in Figure 7, 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 2 O 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.  \nThe right-hand-side of Figure 7 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 2 O 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 2 O, 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.  \nIn 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 (Madhusudhan et al. 2023b) 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 2 O abundance in the troposphere, which in turn could lead to higher CO 2 abundance. However, such a scenario would also give  \nrise to a significant amount of H 2 O and CO, which are presently not detected.", '4.5. 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 7 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 2 s -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 2 s -1 and 10 8 cm 2 s -1 . Based on Madhusudhan et al. (2020) and Valencia et al. (2013), 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 Hu (2021). 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 B.', '4.6. Spectral Characteristics': "We use the atmospheric compositions computed in Section 4.4 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 (Constantinou & Madhusudhan 2024) capability of considering non-uniform vertical mixing ratios, we directly use the atmospheric composition profiles computed using the VULCAN (Tsai et al. 2021) non-equilibrium code described above and shown in Figure 7. 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 Madhusudhan et al. (2023b), 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.  \nThe resulting spectral contributions and transmission spectrum are shown in Figure 8. As can be seen in the  \ntop 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 4.4, the low absolute abundances of both molecules give rise to comparably weak spectral features.  \nThe resulting transmission spectrum provides a reasonable match for the NIRISS observations of K2-18 b at shorter wavelengths due to the strong CH 4 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.  \nOverall, 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 Madhusudhan et al. (2023b), in order for the detected abundance of CO 2 to be compatible with a deep H 2 -rich atmosphere, an unphysically low C/O ratio of ∼ 0.020.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 2 O than H 2 .", '5. 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, magmaatmosphere chemical interactions, and equilibrium as well as disequilibrium (photochemistry and vertical mixing) processes in the atmosphere. Considering all these  \ncoupled factors, it predicts the observable abundances of molecular species in the atmosphere and the expected spectral features.  \nWe apply our framework to perform a comparative assessment of previous works, validating our modelling of magma-atmosphere interactions against Gaillard et al. (2022) and assessing the model predictions of Shorttle et al. (2024), 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 magmaatmosphere interaction models, which do not consider the complex interplay of interior and atmospheric factors, can lead to erroneous results.', '5.1. 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:  \n- · 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 silicatesatmosphere interface and lower atmosphere; and disequilibrium processes throughout the atmosphere.\n- · 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.\n- · Aplanet'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,  \nthese 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.  \n- · 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 CO 2 along with low H 2 O is inconsistent with the chemical expectations for the gas dwarf scenario. Secondly, we find CO to be higher than CO 2 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 magmaatmosphere interactions at the extremely reduced conditions as described in Appendix A. 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.\n- · Overall, we find that key atmospheric signatures for identifying a gas dwarf include the CO and CO 2 abundances, and, if melt is present, possible nitrogen depletion, consistent with some previous studies (cf. Section 3). In particular, we expect that CO/CO 2 > 1 if no H 2 O is observed (as a result, e.g., of condensation), or, in the presence of H 2 O, CO/CO 2 ≲ 1, due to photolysis of H 2 O making more oxygen available for the formation of CO 2 . 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.\n- · Our models predict significant H 2 S for a deep H 2 -rich atmosphere scenario. Hence, a lack of H 2 S 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.\n- · As discussed below, a number of important unknowns remain. In particular, as discussed in Appendix A, the solubility of NH 3 in magma remains poorly understood, especially at extremely reducing conditions, as is also the case for H 2 S at high pressures and temperatures.", '5.2. 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 subNeptunes, and (3) complete reactions lists for all relevant atmospheric species.  \nAs discussed in Appendix A, 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 magmaatmosphere interface on sub-Neptunes. This includes high pressure and temperature, and low oxygen fugacity. In particular, the availability of NH 3 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.  \nFurthermore, 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 (Kite et al. 2019; Schlichting & Young 2022). 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 (Ni et al. 2017). It is however not known whether this behaviour applies to H 2 -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  \nwe have instead assumed a volatile-free melt for the internal structure calculations.  \nThere 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 3 and 5. 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 6. Furthermore, future work will include more detailed prescriptions for the mantle, including alternate mineral compositions, and a fully temperaturedependent EOS for the solid portion.  \nOverall, 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  \na new era in the characterisation of low-mass exoplanets with JWST.  \nAcknowledgements: 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).", 'REFERENCES': 'Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2011, Journal of Physical Chemistry A, 115, 6805, doi: 10.1021/jp109441f  \nArdia, P., Hirschmann, M. M., Withers, A. C., & Stanley, B. D. 2013, Geochimica et Cosmochimica Acta, 114, 52, doi: 10.1016/j.gca.2013.03.028  \nArmstrong, L. S., Hirschmann, M. M., Stanley, B. D., Falksen, E. G., & Jacobsen, S. D. 2015, Geochimica et Cosmochimica Acta, 171, 283, doi: 10.1016/j.gca.2015.07.007  \nAsplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141, doi: 10.1051/0004-6361/202140445  \nAsplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of Astronomy and Astrophysics, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222  \nAzzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, MNRAS, 460, 4063,  \ndoi: 10.1093/mnras/stw1133  \nBarber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828, doi: 10.1093/mnras/stt2011  \nBarber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087, doi: 10.1111/j.1365-2966.2006.10184.x  \nBean, J. L., Raymond, S. N., & Owen, J. E. 2021, Journal of Geophysical Research: Planets, 126, e2020JE006639, doi: 10.1029/2020JE006639  \nBelov, G. V., Iorish, V. S., & Yungman, V. S. 1999, Calphad, 23, 173,  \ndoi: https://doi.org/10.1016/S0364-5916(99)00023-1  \nBenneke, B., Wong, I., Piaulet, C., et al. 2019, The Astrophysical Journal Letters, 887, L14,  \ndoi: 10.3847/2041-8213/ab59dc  \nBenneke, B., Roy, P.-A., Coulombe, L.-P., et al. 2024, arXiv e-prints, arXiv:2403.03325,  \ndoi: 10.48550/arXiv.2403.03325  \nBernadou, F., Gaillard, F., Furi, E., Marrocchi, Y., & Slodczyk, A. 2021, Chemical Geology, 573, 120192, doi: 10.1016/j.chemgeo.2021.120192  \nBorysow, J., Frommhold, L., & Birnbaum, G. 1988, ApJ, 326, 509, doi: 10.1086/166112', 'A. 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 2, and show a selection of them in Figure 9. For most compositiondependent 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 Iacono-Marziano et al. (2012). 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.', 'A.1. Nitrogen Species': "The solubility of N 2 has been explored for a wide range of parameters (e.g., Libourel et al. 2003; Miyazaki et al. 2004; Roskosz et al. 2013; Mallik et al. 2018; Boulliung et al. 2020; Bernadou et al. 2021; Gao et al. 2022), at pressures up to 14.8 GPa and temperatures up to 2800 K (Roskosz et al. 2013). By compiling the available data at P ≤ 8 . 2 GPa and adding their own measurements, Dasgupta et al. (2022) 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 O 2 ∼ IW -2). As warned by Dasgupta et al. (2022), experimental data indicate the solubility seems to reach a plateau, while the law predicts solubility to monotonically increase with pressure. A direct comparison with Roskosz et al. (2013)'s 10 GPa and 14 . 8 GPa data points reveals indeed a true solubility ∼ 1 order of magnitude lower than predicted using Dasgupta et al. (2022)'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 (Dasgupta et al. 2022).  \nIt is also noteworthy that Gao et al. (2022) - whose data was included in the Dasgupta et al. (2022) dataset find some indication of a decrease in the physical solubility of N 2 already at P = 8 GPa. We note that physical solubility is expected to be the dominant solubility mechanism at the oxidized conditions explored by Gao et al. 2022, as opposed to the chemical solubility relevant at reduced conditions (Libourel et al. 2003). Nevertheless, as the relevant quantity is not the total pressure, but rather the nitrogen partial pressure, we believe that the Dasgupta et al. (2022) 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 N 2 mixing ratio in the atmosphere is ≲ 10 -4 in the present models.  \nThe lack of data or simulations for the solubility of NH 3 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.", 'A.2. Carbon Species': "Of the three prominent carbon species (CO 2 , CO, CH 4 ), CO 2 is by far the one for which the most complete experimental data on solubility in magma is available (e.g., Pan et al. 1991; Dixon et al. 1995; Papale et al. 2006; Iacono-Marziano et al. 2012). Considering this wide dataset, for the case of T = 2273 K and a bulk silicate Earth (BSE) melt composition, Suer et al. (2023) find that the solubility of CO 2 is well-approximated by Henry's Law, which they fit to the data. Suer et al. (2023)'s law is in excellent agreement with high pressure ( P ≥ 8 GPa) molecular dynamics simulations by Guillot & Sator (2011) 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.  \nAt lower pressures, the agreement with the simulations is worse, but still within a factor of order unity. In any case, the agreement between Suer et al. (2023)'s law and Guillot & Sator (2011)'s simulations is always satisfactory, which also highlights the weak dependence of the solubility of CO 2 on the melt composition, particularly at P ≤ 8 GPa (Guillot & Sator 2011). Despite the fact that the Suer et al. (2023) 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.  \nFor CO, there is a lack of solubility data at high pressure and temperature. Solubility laws are provided by, e.g., Armstrong et al. (2015) and Yoshioka et al. (2019) (for both MORB and rhyolite melts), both of whom carried out  \nexperiments 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 + CO 2 reaction gets skewed to the right as pressure grows, making an initially pure CO vapour spontaneously become mostly CO 2 at P ≳ 1 GPa (Yoshioka et al. 2019).  \nAn alternative prescription, used by Schlichting & Young (2022) as informed by Hirschmann (2016), is to instead set the solubility of CO to be one third of that of CO 2 . This method, taking Suer et al. (2023)'s BSE law for the CO 2 solubility, yields a CO solubility significantly higher than any of the other laws mentioned so far. This might be due to Suer et al. (2023)'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 (Yoshioka et al. 2019); on the other hand, Yoshioka et al. (2019)'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 Schlichting & Young (2022)'s prescription applied to Suer et al. (2023)'s law results in a solubility between 1 and 2 orders of magnitude higher, depending on the pressure. Ultimately, we use the Yoshioka et al. (2019) MORB law in our calculations, due to it being more recent and calibrated at higher pressures than Armstrong et al. (2015).  \nThe data on CH 4 seems to be even sparser: the only solubility law we encountered in the literature - which we use here - is the one in Ardia et al. (2013), resulting from experiments at 0 . 7 ≤ P ≤ 3 GPa and 1400 ≤ T ≤ 1450 · C, and the Henry's Law fit to their data by Lichtenberg et al. (2021) - the law by Ardia et al. (2013), indeed, follows Henry's law for total pressures P ≲ 10 4 bar (at T = 3000 K, regardless of the CH 4 partial pressure).", 'A.3. Other Volatiles': "The other major volatiles of note are expected to be sulfur species, water, and H 2 , as well as He, which, however, as a noble gas, has no impact on the atmospheric chemistry.  \nChemical equilibrium calculations indicate that, at conditions relevant at the interface, sulfur will be mostly in H 2 S, with little S 2 . For S 2 we use the law by Gaillard et al. (2022). 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 Gaillard et al. (2022)'s predictions were found either, the Woodland et al. (2019) high-pressure data being for a carbonate -silicate melt.  \nFor H 2 S, we found two laws in the literature, by Clemente et al. (2004) - for rhyolite - and by Lesne et al. (2015) - 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: Clemente et al. (2004) find that the solubility of H 2 S moderately increases with increasing temperature, while the law by Lesne et al. (2015) indicates an extremely strong and negative temperature dependence. Furthermore, Lesne et al. (2015) include a dependence on the mole fraction of FeO in the magma, while the law by Clemente et al. (2004) 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 H 2 S at the expected mixing ratios (shown in Figure 6). 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 Clemente et al. (2004).  \nFor H 2 , the law most used in the literature we reviewed is by Hirschmann et al. (2012), 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 Suer et al. (2023) 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 Suer et al. (2023) predicting Henrian behaviour to arbitrary pressure, while Hirschmann et al. (2012)'s laws predict a decline in solubility as pressure increases, consistently with their experimental results.  \nAs the laws given in Hirschmann et al. (2012) have a robust high-pressure experimental background, we use those, and here, specifically, the basaltic melt case.  \nFor H 2 O, there is a great deal of experimental data on the solubility in silicate melts (e.g., Stolper 1982; Silver et al. 1990; Moore et al. 1995; Dixon et al. 1995; Papale et al. 2006; Iacono-Marziano et al. 2012; Sossi et al. 2023), a  \ncomplete review of which is beyond the scope of this work. We focus here on two solubility laws: Sossi et al. (2023), the most recent law available, and Iacono-Marziano et al. (2012), which is the one we choose to implement in our study. Sossi et al. (2023) 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.  \nThese 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 Iacono-Marziano et al. (2012), 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 Sossi et al. (2023) for an H 2 -rich envelope.  \nThe fact that Sossi et al. (2023)'s law depends on a linear combination of the square roots of the fugacities of both H 2 and H 2 O, 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.  \nWe thus consider extrapolating the law of Iacono-Marziano et al. (2012) to higher temperatures to be a more accurate prescription than extrapolating that by Sossi et al. (2023) to high pressures, and hence do so here, assuming an Etna basalt composition for the melt. This choice is also consistent with that in Gaillard et al. (2022).", 'A.4. Summary': 'Data on solubility in silicate melt are available, at some conditions, for several species of interest, with the one exception being NH 3 , 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 2, and we show a selection of them in Figure 9. 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 N 2 data at both very high temperature and pressure exists, but that is for relatively oxidised conditions (Roskosz et al. 2013). High-pressure ( P ≥ 10 5 bar) simulations exist for CO 2 , but only at T ≤ 2273. For S 2 , for a co-existing fluid phase, high-pressure data only exists at low temperature, and only for carbonate-silicate melt (Woodland et al. 2019). All other species seem to lack high-pressure data.  \nExploring 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.', 'B. SENSITIVITY TO ATMOSPHERIC PARAMETERS': "As described in Section 4.5, 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 7, 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 2 s -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 2 s -1 and 10 8 cm 2 s -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 Hu (2021). Disequilibrium effects due to photochemistry and vertical mixing are included in all cases discussed here.  \nWe start with investigating departures from the canonical C2 case, as shown in Figure 7. 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 10, along with those for 50 × metallicity from Figure 7 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.  \nNext 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 2 s -1 , with our canonical value at 10 6 cm 2 s -1 . The metallicity remains fixed at the canonical value of 50 × solar. As shown in Figure 11, both the higher and lower K zz values negligibly affect the computed mixing  \nFigure 9. Behaviour of selected solubility laws in silicate melt for some prominent molecules. Solid blue lines indicate the laws used in this study. All temperature-dependent laws are shown here for T = 3000 K; the ∆IW-dependent N 2 and S 2 laws are shown for ∆IW = -6 . 4; all laws with free composition-dependent parameters are shown for the Mt. Etna basalt composition (Iacono-Marziano et al. 2012). For H 2 O, we label as (a) the law given in Sossi et al. (2023) for ϵ 3550 = 6 . 3 m 2 /mol, and as (b) that for ϵ 3550 = 5 . 1 m 2 /mol. The mixing ratio of each species, informed by Section 4.3, is indicated in the title of the respective subplot; the x-axis indicates total pressure. The gases are treated ideally, i.e. we take the fugacity of each species' to be equal to its partial pressure. We do not plot the CH 4 solubility, for which Ardia et al. (2013) - whose law we use - predicts a much lower value than the other carbon-bearing species, nor the solubility of H 2 S, for which both Clemente et al. (2004) and Lesne et al. (2015) predict much lower solubility than for S 2  \n<!-- image -->  \nTable 2. Sources of solubility data and laws considered in this study. The (total) pressure and temperature ranges indicated are those corresponding to the experiments carried out in the respective studies, or, if the studies calibrate a solubility law based on data from previous works, the range spanned by those.  \n- a : Molecular dynamics simulation.\n- b : Gaillard et al. (2022) state that their law is calibrated against data obtained with gas at a pressure of one atmosphere. However, they refer to Zajacz et al. (2013), whose experiments were carried out at 200 MPa.\n- c : For carbonate-silicate melt.  \nratios 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 11.  \nWe 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 1 by factors of 2 and 6, respectively. The results are shown in Figure 12. 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.  \nAs 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 2 s -1 in the deep atmosphere. The resulting vertical mixing ratio profiles are shown in Figure 13 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 .  \nFigure 10. Effect of enhanced metallicity on the vertical mixing ratio profiles for several H-C-O-N-S molecular species with the C2 P -T profile. Dotted lines show the profile from the left-hand side of Figure 7, for C2 with 30% silicates, equivalent to 50 × solar elemental abundances. Left: solid lines indicate the corresponding profiles for 100 × solar metallicity. Right: solid lines indicate the corresponding profiles for 300 × solar metallicity.  \n<!-- image -->  \nThus 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 (e.g. Hu 2021). We adopt the P -T profile of Hu (2021) 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 1) and our canonical K zz treatment, and 2) a high K zz = 10 8 cm 2 s -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 2 O abundance exceeds the retrieved value for both cases of metallicity and K zz considered.  \nOverall, 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 (Madhusudhan et al. 2023b).  \nFigure 11. Effect of varying K zz on the vertical mixing ratio profiles for several H-C-O-N-S molecular species for the C2 P -T profile. Dotted lines show the profile from the left-hand side of Figure 7, for C2 with 30% silicates, equivalent to 50 × solar elemental abundances with our canonical treatment of K zz , with a value of 10 6 cm 2 s -1 in the deep atmosphere. Left: solid lines indicate the corresponding profiles with a lower K zz = 10 4 cm 2 s -1 in the deep atmosphere. Right: solid lines indicate the corresponding profiles with a higher K zz = 10 8 cm 2 s -1 in the deep atmosphere.  \n<!-- image -->  \nFigure 12. Effect of enhanced metallicity on the vertical mixing ratio profiles for several H-C-O-N-S molecular species with the C3 P -T profile. Dotted lines show the profile from the right-hand side of Figure 7, for C3 with 30% silicates, with N depletion due to the presence of magma. Left: solid lines indicate the corresponding profiles for 100 × solar metallicity, i.e. 2 × the respective elemental abundances given in Table 1. Right: solid lines indicate the corresponding profiles for 300 × solar metallicity, i.e. 6 × the respective elemental abundances in Table 1.  \n<!-- image -->  \nFigure 13. Effect of high metallicity and high K zz on the vertical mixing ratio profiles for several H-C-O-N-S molecular species. Dotted lines show the canonical profiles from Figure 7. Left: solid lines indicate the profiles with the C2 P -T profile with 300 × solar abundance and K zz = 10 8 cm 2 s -1 in the deep atmosphere. Right: solid lines indicate the profiles with the C3 P -T profile with 6 × the elemental abundances from Table 1, i.e. equivalent to 300 × solar abundance, and K zz = 10 8 cm 2 s -1 in the deep atmosphere.  \n<!-- image -->  \nFigure 14. Effect of higher T int on the vertical mixing ratio profiles for several H-C-O-N-S molecular species. We adopt the P -T profile from Hu (2021) with T int of 60 K and 100 × solar metallicity, extrapolated to 1000 bar. Dotted lines show the canonical profiles using the hotter C3 profile from Figure 7. Left: solid lines indicate profiles assuming 100 × solar metallicity including N depletion, adopting K zz = 10 6 cm 2 s -1 in the deep atmosphere, for the Hu (2021) P -T profile. Right: solid lines indicate profiles assuming 300 × solar metallicity including N depletion, adopting K zz = 10 8 cm 2 s -1 in the deep atmosphere, for the Hu (2021) P -T profile.  \n<!-- image -->"}