Most exoplanets with a few Earth radii are more inflated than barerock planets with the same mass indicating a substantial volatile amount. Neither the origin of the volatiles nor the planets bulk composition can be constrained from the massradius relation alone and the spectral characterization of their atmospheres is needed to solve this degeneracy. Previous studies showed that chemical interaction between accreted volatile and possible molten rocky surface i.e. magma ocean can greatly affect the atmospheric composition. However a variety in the atmospheric compositions of such planets with different properties remains elusive. In this work we examine the dependence of atmospheric H O and C on planetary parameters atmospheric thickness planetary mass equilibrium temperature and magma properties such as redox state assuming nebula gas accretion on an Earthlike core using an atmospheremagma chemical equilibrium model. Consistent with previous work we show that atmospheric HSUB2SUBO fraction on a fully molten rocky interior with an Earthlike redox state is on the order of 10SUP2SUP10SUP1SUP regardless of other planetary parameters. Despite the solubility difference between H and Cbearing species CH increases only a few times above the nebula value except for atmospheric pressure 1000 bar and HSUB2SUBO fraction 10. This results in a negative OHCO trend and depleted CO below onetenth of the nebula gas value under an oxidized atmosphere which could provide a piece of evidence of rocky interior and endogenic water. We also highlight the importance of constraints on the highpressure material properties for interpreting the magmaatmospheric interaction of inflated planets.
['Exoplanet atmospheres', 'Exoplanet formation', 'Extrasolar rocky planets', 'Planetary interior', '487', '492', '511', '1248', 'Astrophysics - Earth and Planetary Astrophysics']
Role of Magma Oceans in Controlling Carbon and Oxygen of SubNeptune Atmospheres
{'Role of magma oceans in controlling carbon and oxygen of sub-Neptune atmospheres': '<!-- image -->  \n1 Department of Astronomical Science, School of Physical Sciences, Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  \n- 2 Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  \n3 Faculty of Science and Technology, Sophia University, Kioi-Cho 7-1, Chiyoda-ku, Tokyo, 102-8554 Japan  \n- 4 Department of Physics and Astronomy, University, College London Gower Street, WC1E 6BT London, United Kingdom', 'ABSTRACT': "Most exoplanets with a few Earth radii are more inflated than bare-rock planets with the same mass, indicating a substantial volatile amount. Neither the origin of the volatiles nor the planet's bulk composition can be constrained from the mass-radius relation alone, and the spectral characterization of their atmospheres is needed to solve this degeneracy. Previous studies showed that chemical interaction between accreted volatile and possible molten rocky surface (i.e., magma ocean) can greatly affects the atmospheric composition. However, a variety in the atmospheric compositions of such planets with different properties remains elusive. In this work, we examine the dependence of atmospheric H, O, and C on planetary parameters (atmospheric thickness, planetary mass, equilibrium temperature, and magma properties such as redox state) assuming nebula gas accretion on an Earth-like core, using an atmosphere-magma chemical equilibrium model. Consistent with previous work, we show that atmospheric H 2 O fraction on a fully molten rocky interior with an Earth-like redox state is on the order of 10 -2 -10 -1 regardless of other planetary parameters. Despite the solubility difference between H- and C-bearing species, C/H increases only a few times above the nebula value except for atmospheric pressure ≲ 1000 bar and H 2 O fraction ≳ 10%. This results in a negative O/H-C/O trend and depleted C/O below one-tenth of the nebula gas value under an oxidized atmosphere, which could provide a piece of evidence of rocky interior and endogenic water. We also highlight the importance of constraints on the high-pressure material properties for interpreting the magma-atmospheric interaction of inflated planets.  \nKeywords: Exoplanet atmospheres(487) - Exoplanet formation(492) - Extrasolar rocky planets(511) - Planetary interior(1248)", '1. INTRODUCTION': "Past exoplanet surveys have revealed that exoplanets within the radius range spanning Earth to Neptune-sized bodies are abundant (e.g., Howard et al. 2010, 2012; Fressin et al. 2013; Hsu et al. 2019). While some have a mass-radius relationship consistent with the bulk composition similar to that of Earth, others have larger radii than those of bare Earth-like planets with the same masses, suggesting that they contain thick layers of volatile species (e.g., Weiss & Marcy 2014; Rogers 2015; Owen 2019). The former tend to be called super-Earths, while the latter are often called sub-Neptunes, although the exact definition of these terms can depend on the context.  \nThe origin of these intermediate-sized planets close to their host stars has been the subject of active debate. One hypothesis posits that super-Earths and sub-Neptunes share a common origin characterized by a rocky core enveloped by an H-He atmosphere accreted from the protoplanetary disk and finds that such a model can explain the distinctive radius gap in observed exoplanets' radii at approximately 1 . 8 R ⊕ (e.g., Fulton et al. 2017; Van Eylen et al. 2018; Petigura et al. 2022). Previous works have discussed the possible origins of this radius gap, including atmospheric  \nescape, impact loss, and accretion processes (e.g., Owen 2019; Wyatt et al. 2020; Lee et al. 2022). Another hypothesis proposes that the planets with larger radii accreted a large amount of icy planetesimals beyond the so-called snow lines (e.g., Zeng et al. 2019; Mousis et al. 2020; Venturini et al. 2020; Zeng et al. 2021; Luque & Pall'e 2022; Burn et al. 2024). An examination of these proposed scenarios for the formation of such close-in small planets can be found in the review by Bean et al. (2021), along with pertinent references contained therein.  \nGiven that the mass-radius relationship alone cannot uniquely constrain the bulk composition (e.g., Rogers & Seager 2010; Valencia et al. 2013; Rogers et al. 2023), more clues to the origins of sub-Neptunes are expected to be obtained through atmospheric spectroscopy. This is particularly true as the James Webb Space Telescope ( JWST ) data become available with higher precision than past observations and start to constrain the atmospheric compositions (typically at the 0.01-1 bar levels) of sub-Neptunes/super-Earths (e.g., Moran et al. 2023 for GJ 486b; Madhusudhan et al. 2023 for K2-18b; May et al. 2023 for GJ 1132b; Hu et al. 2024 for 55 Cancri e; Benneke et al. 2024 for TOI-270 d).  \nWhen interpreting these atmospheric data, it is important to note that the atmospheric composition of these small planets could result from interactions between accreted volatile species and the planetary core. Therefore, studying how the accreted volatiles are altered and partitioned within the planet is critical (e.g., Hirschmann 2012; Grewal et al. 2020; Sossi et al. 2020; Deng et al. 2020; Gaillard et al. 2022). Several review papers (e.g., Gaillard et al. 2021; Wordsworth & Kreidberg 2022; Suer et al. 2023; Lichtenberg et al. 2023; Salvador et al. 2023; Lichtenberg & Miguel 2024) provide a comprehensive survey of the previous studies on the impacts of core-atmospheric chemistry on terrestrial planet atmospheres, starting with the solar system planets-Earth, Venus, and Mars-and more recently extending to larger rocky exoplanets such as super-Earths and sub-Neptunes.  \nOne of the new insights in the context of exoplanets is that the magma (and iron core) of sub-Neptunes serve as crucial volatile reservoirs, especially for hydrogen, through the dissolution of significant amounts of H 2 and H 2 O. If the total amount of volatiles is within a few percent, this role of the planetary core as a volatile reservoir significantly influences the planet's evolution and causes a decoupling between the amount of volatiles in the interior and the planetary radius (e.g., Chachan & Stevenson 2018; Kite et al. 2019, 2020; Schlichting & Young 2022). In addition to the radius-volatile content decoupling effect, the ingassing of certain volatile species affects the atmospheric composition by changing the elemental abundances of the atmosphere. For instance, Shorttle et al. (2024) used an atmosphere-magma model of K2-18b to demonstrate that the high solubility of nitrogen in magma in a reducing environment can deplete N from the atmosphere depending on the redox state. Based on this result, they argued that using N depletion as the sole indicator for the existence of a water ocean can be misleading.  \nAnother process through which the rocky interior affects the atmospheric composition is the redox reactions between the atmosphere and the magma. Notably, H 2 , the major atmospheric species of sub-Neptunes, can be oxidized by magma to produce substantial amounts of H 2 O(e.g., Sasaki 1990; Ikoma & Genda 2006). Kite et al. (2020) showed that for sub-Neptunes with relatively thin atmospheres, such alteration of atmospheric composition by magma challenges the interpretation of the atmospheric composition, as a certain H 2 Omixing ratio can be explained by different combinations of the accreted volatile composition and the redox state of magma. Kite & Schaefer (2021) proposed that water-rich atmospheres originated from the interaction between the primordial nebula-driven gas and magma would have the C/O ratio smaller than the value expected from the ice solids, which can be interpreted as the evidence of magmaatmosphere reaction. The C/O ratio change due to magma was also recently discussed by Tian & Heng (2024) using their magma-atmosphere H-C-N-O-S chemical reaction model. They argued that the C/O ratio is difficult to use as an indicator of the formation history of rocky planets due to its susceptibility to internal chemistry variations. Another parameter that can imply the presence of magma is the CO 2 /CO ratio. Shorttle et al. (2024) suggested that a small CO 2 /CO ratio in the deep atmosphere of K2-18b can indicate the presence of high-temperature magma that reacts with the atmosphere. However, the CO 2 /CO ratio at the altitudes observed through the transmission spectroscopy can differ from that in the deep atmosphere, and the former varies depending on the elemental abundance of the atmosphere and the efficiency of the vertical mixing rather than directly constraining the presence of magma.  \nWhile these previous works pointed out the importance of chemical interactions between the atmosphere and the planetary core in determining the final state, they did not explicitly examine how the final atmospheric composition depends on assumptions about planetary parameters. Several previous studies (e.g., Tian & Heng 2024; Shorttle et al. 2024; Maurice et al. 2024) have calculated the atmospheric composition of the H-C-N-O-(S) elemental system by assuming the final redox state of the magma. However, the final redox state of the system fundamentally depends on the relative amount of reactive magma compared to the atmosphere, the composition of the magma (iron content), and the iron speciation prior to the chemical reaction. These are parameters that are also relevant to planet formation scenarios.  \nAlthough these magma properties have been considered in some existing studies (e.g., Kite et al. 2020; Kite & Schaefer 2021), quantitative analyses for a wider range of elements (such as Carbon) are limited. Given the accumulating data from atmospheric spectroscopy of super-Earths and sub-Neptunes, it is also important to investigate the trends of the main atmospheric components, represented by O/H and C/O, as functions of various observable planetary parameters related to atmospheric structure (e.g., planetary radius, equilibrium temperature, mass).  \nIn this paper, we study such trends based on chemical-equilibrium calculations between the atmosphere and underlying magma, focusing on the scenario where sub-Neptunes form through the accretion of nebula gas onto rocky cores interior to the H 2 O snowline (Owen & Wu 2017; Modirrousta-Galian & Korenaga 2023). We show specific relations between O/H and C/O and their dependencies on magma properties (effective amount, Fe fraction, and initial Fe speciation), observable planetary parameters (planetary mass, equilibrium temperature, and radius), and an additional parameter for determining atmospheric structure (Radiative-Convective Boundary). Furthermore, we introduce the concept of planetary evolution using a toy model and discuss the impact of magma solidification on atmospheric composition. We explore how the O/H-C/O relation and their ranges could potentially help to distinguish the origins of sub-Neptunes from other major formation scenarios, such as the assembly of icy planetesimals. In addition to the numerical calculations with the non-ideal behavior of major species, analytic expressions will be presented for future use.  \nIn section 2, we introduce the structure, chemical reactions, and other assumptions of our sub-Neptune model. Section 3 presents the main results, showing the range of atmospheric compositions and their dependence on the planetary parameters. In section 4, we will address the implications of our findings on the formation scenario, potential spectral features, and the impact of different atmospheric structure models and atmospheric escape mechanisms on atmospheric composition. Section 5 summarizes our discussion.", '2. MODEL': 'Our study aims to investigate atmospheric compositions focusing on H-, O-, and C-bearing species under the formation scenario where the nebula gas accretes onto a rocky core (e.g., Owen & Wu 2017). The structure of our model follows that of Kite et al. (2020). However, we also introduce several modifications to their model. Specifically, we incorporate C-bearing species into the analysis and derive the surface conditions (temperature and pressure) from the atmospheric temperature profile and observable planetary parameters such as planetary mass, radius, and equilibrium temperature derived from stellar irradiation without using the pre-assumed value. We also introduce an additional buffer system, the Ferric-Ferrous equilibrium, alongside the Iron-Wustite buffer to expand the range of affordable redox states towards a more oxidized regime.  \nThe model components and assumptions are described in depth in this section.', '2.1. Vertical structure': "Our planetary model consists of three layers: an atmosphere, a silicate layer, and an innermost iron core, as shown in Figure 1.  \nWe assume a vertically homogeneous composition in the atmosphere when calculating the bulk atmospheric composition for simplicity. The elemental ratio of helium to hydrogen in the atmosphere (He/H) is fixed at 0.1, reflecting the nebula gas value (0.079, from Lodders 2003). This assumption implies that the partitioning of H 2 to the magma is comparable to that of He, which could be supported by the similar mean molecular radius between He and H 2 (e.g., Hirschmann et al. 2012; Suer et al. 2023). We note that there remains uncertainties in this assumption because the low-pressure experimental study shows smaller solubility of He than H 2 (e.g., Jambon et al. 1986). However, since the solubility of He does not significantly affect our results, we assume a fixed nebula gas-like atmospheric He/H ratio for simplicity. Only when we demonstrate the transmission spectra in Section 4.2 do we model the vertical profile of atmospheric molecules considering atmospheric chemistry and eddy diffusion (considered indirectly by assuming specific quenching points.)  \nAs the fiducial model, the atmospheric thermal profile is assumed to be composed of an isothermal stratosphere with equilibrium temperature, T eq , above the radiative-convective boundary (RCB) and the dry adiabatic troposphere between the RCB and the bottom of the atmosphere (BOA). The planetary radius, R p , is set at the pressure level of 0.1 bar, consistent with the altitude of the optical radius of planets (e.g., Kawashima & Ikoma 2018). The pressure of the RCB, P RCB , is treated as a variable and is either 0.1, 1, 10, or 100 bars. This simplified Temperature-Pressure profile allows for easy analytic estimates and helps us clarify the dependence on planetary parameters in Section 3.3.  \nFigure 1. The schematic figure of the vertical structure of the sub-Neptune model presented in this study. The black line within the atmosphere layer illustrates the assumed temperature-pressure relationship in our model. Chemical species involved in internal reactions are written in red within each structural layer. The RCB stands for the Radiative-Convective Boundary, while the BOA denotes the Bottom of the Atmosphere.  \n<!-- image -->  \nWe will also discuss the difference between our simplified thermal profile and a more realistic radiative-convective equilibrium profile in Section 4.4.  \nIn our model, the atmosphere is assumed to be in hydrostatic equilibrium and divided into 1000 layers, evenly spaced in altitude between the RCB and BOA, and 100 layers evenly spaced in pressure between the R p and the RCB. With this layered setup, we account for the altitude dependence of the gravitational acceleration g and the temperature dependence of the heat capacity C p when solving for structure assuming hydrostatic equilibrium. The lapse rate in the lower dry adiabatic troposphere is computed from g/C p , and for C p , we take into account the temperature dependence of the C p of H 2 , H 2 O, and He, which account for over 95% of the atmosphere in our results. We refer to the JANAF thermochemical table (Chase 1998) for the C p of the three species, which covers the temperature range from 100 K to 6000 K.  \nThe innermost iron core is assumed to take up 34% of the planetary mass, consistent with Earth's iron core mass fraction, and the remaining 66% of the planetary mass is attributed to the silicate layer. We ignore the atmospheric mass in this study considering the dominant mass of the planetary core. We calculate the radius of the planetary (silicate + iron) core using the mass-radius relationship equations of rocky cores from Fortney et al. (2007) (see their Figure 4 and Equation 8). Note that our study currently omits the chemical impact of the iron core. This assumption would align with scenarios where the segregation of the iron core occurred early, preventing the iron core from being involved in reactions with magma. The efficacy of iron core segregation depends on various factors, such as heat flow and iron droplet size (e.g., Lichtenberg 2021). In cases of sufficiently vigorous magma convection, it has been argued that the formation of the iron core may be delayed, allowing sufficient time for the iron core (and drifting iron droplets) to attain equilibrium with the surrounding magma and atmosphere (e.g., Lichtenberg 2021). The potential impact of the innermost iron core on volatile partitioning is described in Section 4.5.  \nThe silicate layer is stratified into an upper molten magma and a lower solidified rock layer. Under nominal conditions, magma convection is assumed vigorous enough to ensure homogeneous mixing and fully reactive magma. The dependence of atmospheric composition on the reactive amount of magma is discussed in Section 3.2.  \nIn the fiducial model, we assume the elemental composition of the magma to be the same as the Mid Ocean Ridge Basalt (MORB), with elemental ratios of Fe / Si = 0 . 17 and Mg / Si = 0 . 22 (Gale et al. 2013, and references therein). We note that this MORB-like composition differs from the peridotite-like Bulk Silicate Earth (BSE), which is characterized by a relatively higher magnesium fraction (Mg / Si ∼ 0 . 92, Javoy et al. 2010). This compositional difference can influence the equilibrium state by altering the abundance of reactive species (e.g., FeO) and the solubility of gases such as H 2 , H 2 O, and CO 2 (e.g., Hirschmann et al. 2012; Allison et al. 2022). We explore the dependence of atmospheric composition on the Fe/Si ratio in Section 3.2 and discuss the potential effects of varying gas solubility (H 2 , H 2 O, and CO 2 ) in Section 4.5.  \nThe lower solidified rock layer remains chemically inert and is assumed to be free of volatiles. We presume the molten fraction of the silicate layer to be 100% and introduce solidified magma layers into the model when discussing the magma solidification effect. Section 2.5 describes in detail how we account for magma solidification in this study.  \nImportantly, we calculate the chemical equilibrium state (Section 2.2 and Section 2.3 below) using the temperature and pressure conditions at the interface between the magma and the atmosphere. In other words, we assume that chemical reactions occur at the magma surface and the lowest atmospheric layer. In reality, the equilibrium state can be dependent on the pressure and the saturation amount of each volatile species that can vary with depth in the magma. Such effect will be the scope of the future work.", '2.2. Redox reactions': 'In our chemical equilibrium calculation, we consider H 2 , H 2 O, O 2 , CO, CO 2 , CH 4 , Fe, FeO, and Fe 2 O 3 as participating species. Note that He is one of the major gas components determining the atmospheric profile but does not participate in the chemical reactions. Its abundance remains fixed at the nebula gas value (He/H = 0.1), as mentioned in Section 2.1.  \nWe consider the redox equilibrium among these species as follows (see also Schlichting & Young 2022). In the atmosphere,  \n- · 2 H 2 + O 2 ⇆ 2 H 2 O,\n- · 2 CO + O 2 ⇆ 2 CO 2 ,\n- · CH 4 + H 2 O ⇆ CO + 3 H 2 .  \nIn the molten silicate layer, we consider the Iron-Wustite (IW) reaction (as assumed in Kite et al. 2020) and the Ferric-Ferrous reaction (in the most oxidized case) as the redox buffers:  \n- · 2 Fe + O 2 ⇆ 2 FeO\n- · 2 FeO + 0.5 O 2 ⇆ Fe 2 O 3  \nWe include the non-ideal behavior of atmospheric species by utilizing the fugacity ( f S, where S represents each species) instead of partial pressure in our equilibrium calculations. The equilibrium between the molten rock and the atmosphere indicates that the fugacity of O 2 ( f O 2 ) in these equations is equated. The fugacity of the species f S is determined by f S = ϕ S P S , where ϕ S is the fugacity coefficient and P S is the partial pressure. We calculate the fugacity coefficient of gases from the Equation of State (EoS) formulation by Zhang & Duan (2009).  \nFor reactions between fluids, the equilibrium constant is determined using the Gibbs free energy:  \nK p = exp ( -∆ G · RT ) (1)  \nwhere K p denotes the equilibrium constant, ∆ G · is the difference in standard Gibbs free energy of formation between the reactant(s) and the product(s), R is the gas constant, and T is the temperature. The standard Gibbs free energy of formation for volatile species is again adopted from the JANAF table (Chase 1998).  \nFor the reactions in the molten silicate layer, we include the change in volume of the reaction (∆ V ) term:  \nK p = exp ( -∆ G · + ∫ P P 0 ∆ V dP RT ) . (2)  \nThe ∆ G · of the IW buffer is obtained from Kowalski & Spencer (1995), which is validated for temperatures up to 3000 K. We refer to Komabayashi (2014) for the EoS of liquid Fe, validated up to 3300 K and 10 6 bar. For the EoS of liquid FeO, we utilize the data from Armstrong et al. (2019), which is validated under conditions up to 3000 K and 3 × 10 5 bar. We consider the non-ideal behavior of liquid FeO by introducing the activity coefficient of 1.5 (Holzheid et al. 1997; Kite et al. 2020). On the other hand, we assume that liquid Fe remains in a pure state within the magma, corresponding to an activity coefficient of 1.  \nIn the case of the Ferric-Ferrous equilibrium, we adopt the equilibrium constant equation proposed by Hirschmann (2022). This equation incorporates the EoS of liquid Fe from Komabayashi (2014) and Fe 2 O 3 from Armstrong et al.  \n(2019) and Deng et al. (2020). Additionally, it includes a term describing the dependence of the activity coefficients of FeO and Fe 2 O 3 on the silicate composition based on multiple experimental studies rather than assuming a fixed value. The suggested equation from Hirschmann (2022) has the form of  \nlog ( X FeO 1 . 5 X FeO ) = a log f O 2 -( ∆ G · + ∫ P P 0 ∆ V dP ) RT ln(10) -log ( γ FeO 1 . 5 γ FeO ) (3)  \nwhere a is a stoichiometric coefficient and γ is the activity coefficient, which is the function of temperature and silicate composition. In Hirschmann (2022), a has an empirical value of 0.1917. We note that this is different from the a = 0 . 25 of the ideal case, which we can expect from the reaction (FeO + 0.25 O 2 ⇆ FeO 1 . 5 ). The validation range of temperature and pressure for this equation extends up to 2200 K and 70000 bar, respectively.  \nThe f O 2 parameterization in magma depends on the metal saturation in the magma ocean (e.g., Itcovitz et al. 2022). We assume the [Fe], the molar fraction of iron in the metal phase, is always 1 (pure iron). The IW buffer has a limited range of f O 2 it can cover. Under the reduced condition (low f O 2 ), pure iron exists in the magma, and both the IW buffer and the Ferric-Ferrous equilibrium occur concurrently (Fe ←→ FeO ←→ Fe 2 O 3 ). As the magma becomes more oxidized, the Fe oxidizes to FeO, and FeO changes further to Fe 2 O 3 . Under the oxidized condition where the f O 2 level exceeds the range covered by the IW buffer, pure iron becomes depleted, and the ferric-ferrous equilibrium solely determines the f O 2 (FeO ←→ Fe 2 O 3 ).  \nFigure 2 illustrates the relationship between iron speciation and the oxygen fugacity of the system. Here, we parameterize the iron speciation of magma by the number ratio between the oxygen bound to iron and iron itself, denoted as N O(Fe) / N Fe . Depending on the oxidation state, N O(Fe) / N Fe ranges from 0 (all iron exists as pure Fe) to 1.5 (all iron exists as Fe 2 O 3 ). The figure clearly shows the O-rich magma corresponds to higher f O 2 . We utilize N O(Fe) / N Fe as the primary parameter to indicate the redox state of the magma hereafter.  \nWe note that the f O 2 shows higher sensitivity to the N O(Fe) / N Fe after the depletion of pure Fe due to the change in governing reactions and the specific parameterization, N O(Fe) / N Fe , that we use. From the chemical equations of IW buffer and Ferric-Ferrous equilibrium, the equilibrium constants for the two buffers are  \nK IW = ( γ FeO [FeO]) 2 f O 2 (4)  \nK FF = ( γ FeO 1 . 5 [FeO 1 . 5 ] γ FeO [FeO] ) 1 /a 1 f O 2 (5)  \nwhere a represents the stoichiometric coefficient of O 2 in the Ferric-Ferrous equilibrium, and K IW and K FF denote the equilibrium constants of the IW buffer and the Ferric-Ferrous equilibrium, respectively. From the two equilibrium constant equations, it is clear that the proportional relationship shifts from log( f O 2 ) ∝ log([FeO]) to log( f O 2 ) ∝ 1 a log( [FeO 1 . 5 ] [FeO] ) as the IW buffer is excluded from the reaction. Therefore, f O 2 under oxidized conditions rapidly changes as N O(Fe) / N Fe increases because the ratio between two iron species determines f O 2 .  \nWhile f O 2 also is related to the amount of O 2 , we do not track O 2 . f O 2 serves only as the variable that indicates the redox state of the system.  \nFigure 3 illustrates the temperature and pressure dependencies of equilibrium constants for the IW buffer, FerricFerrous equilibrium, H 2 -H 2 O, and CO-CO 2 reactions. All reactions share similar temperature trends, becoming more reduced as temperature increases, although the sensitivities vary depending on the reactions. The IW buffer and Ferric-Ferrous equilibrium show pressure dependence when the pressure exceeds ∼ 10 4 bars. It is worth mentioning that the H 2 -H 2 O and IW buffer exhibit similar temperature dependencies, suggesting a smaller temperature effect for the combined reaction (H 2 + FeO ⇆ H 2 O+ Fe) compared to each reaction.  \nAs shown in Figure 3, the temperature validation ranges of the IW and Ferric-Ferrous buffers are much narrower compared to the H 2 -H 2 O and CO-CO 2 reactions, extending only up to 3000 K and 2200 K. To discuss the results with higher surface temperature, we extrapolate the equilibrium constants of the two buffers as follows. We extrapolate the equilibrium constants up to 4000 K, based on the negligible difference (see Figure 3) between the two IW equilibrium constants from JANAF and Kowalski & Spencer (1995) and embracing arguments from Hirschmann (2022) regarding the thermodynamic model of Ferric-Ferrous equilibrium. Beyond 4000 K, we employ boundary values for the equilibrium constants of all redox reactions. In Appendix A, we confirmed that the current extrapolation setting does not yield notable differences from the boundary-value extrapolation for all equilibrium constants over 2200 K (see Appendix A for details).  \nFigure 2. The relation between the iron speciation and the oxygen fugacity. The x-axis is the number ratio between the oxygen bound to iron and iron itself, while the y-axis denotes the oxygen fugacity ( f O 2 ). The cyan (red) plot indicates the oxygen fugacity calculated from the IW buffer (Ferric-Ferrous equilibrium). The black dashed line represents the actual f O 2 of the system. The magenta-colored line shows the Fe 3+ / Fe ratio corresponding to each N O(Fe) /N Fe . The temperature and pressure are set at 2500 K and 1000 bar, respectively.  \n<!-- image -->', '2.3. Dissolution of volatile species into magma': "Here we describe the solubility law for each species used in this paper. We adopt the following law for the H 2 solubility in basaltic silicate, based on the experiments conducted at temperatures of 1673 to 1723 K and pressures ranging from 7000 to 30000 bar (Hirschmann et al. 2012):  \nX H 2 = f H 2 [bar] e -11 . 403 -0 . 76P[GPa] (6)  \nwhere X H 2 is the H 2 mole fraction, f H 2 is the H 2 fugacity (in bar), and P is the atmospheric pressure (in GPa). This can be rewritten as:  \nX H 2 ∼ 0 . 11 · s H 2 ( P ) ( P H 2 1 GPa ) (7)  \ns H 2 ( P ) ≡ ϕ H 2 (0 . 47) P [GPa] (8)  \nFor the H 2 O solubility to the basaltic magma, we refer to the equation from Schaefer et al. (2016). They applied a curve fit to the H 2 O solubility data for various silicate melts from Papale (1997) which has a temperature range of 900 to 1670 K and a pressure range between 100 and 10000 bar:  \nw H 2 O = 3 . 44 · 10 -8 ( P H 2 O [Pa]) 0 . 74 (9)  \nFigure 3. The temperature-pressure dependencies of the equilibrium constants for significant reactions. The x-axis represents temperature, while the y-axis indicates the equilibrium constant. Different colored lines represent the equilibrium constants of individual reactions, as specified in the legend. The solid line corresponds to the default pressure of 1 bar, while dashed and dotted lines denote the IW and Ferric-Ferrous buffer equilibrium constants at pressures of 10000 bar and 70000 bar (the upper limit of the validation range). The lines with faded colors denote temperatures beyond the validated range, where extrapolation is applied. The black solid line represents the equilibrium constant of the IW buffer at 1 bar, calculated using the JANAF table.  \n<!-- image -->  \nwhere w H 2 O is the weight fraction of H 2 O, and P H 2 O is the H 2 O partial pressure. This may be rewritten as:  \nX H 2 O =0 . 16 · s H 2 O ( P H 2 O 1 GPa ) 0 . 74 (10)  \ns H 2 O ≡ ( µ m µ H 2 O ) (11)  \nwhere µ m is the molar mass of magma. Note that the solubility of H 2 O does not linearly depend on the P H 2 O . It implies the H 2 O dissolution into silicate melt is not completely physical dissolution, and the chemical reaction between the magma and H 2 O fluid changes the solubility law. More specifically, H 2 O has two dissolution mechanisms in silicate (e.g., Kuramoto & Matsui 1996; Zhang 1999; Liu et al. 2005; Mitchell et al. 2017; Suer et al. 2023):  \n- · H 2 O fluid ⇆ H 2 O melt\n- · H 2 O melt + O 2 -melt ⇆ 2OH -melt  \nThe former reaction is the physical dissolution of H 2 O into the melt, which follows Henry's law, while the latter reaction establishes a proportionality between the dissolved H 2 O and the square root of partial pressure of H 2 O ( ∝ P 0 . 5 H 2 O ). Therefore, the solubility of H 2 O has an exponent between 0.5 and 1, with the exact value depending on the dissolved amount of OH -in the magma. A high concentration of hydroxide ions in magma makes the physical dissolution of H 2 O more significant. This transition is understood to occur at a boundary of a few weight percentages of dissolved H 2 O (e.g., Stolper 1982; Papale 1997; Mitchell et al. 2017; Suer et al. 2023).  \nThe solubility law for CO 2 is adopted from the fitting function given by Pan et al. (1991), which has a pressure range of 1000 to 15000 bar and a temperature range of 1473 to 1873 K. The original equation is :  \nw CO 2 [%] = 0 . 00119 + 4 . 81438 · 10 -5 · P CO 2 [bar] +5 . 019505 · 10 -9 · P 2 CO 2 -2 . 587138 · 10 -13 · P 3 CO 2 +5 . 96362 · 10 -18 · P 4 CO 2 -5 . 67816 · 10 -23 · P 5 CO 2 (12)  \nwhere w CO 2 is the weight fraction of CO 2 . However, the above fitting function may not be a good approximation at small P CO 2 because the mass fraction of CO 2 in silicate does not converge to zero, even at zero bar. Therefore, we use a re-fitted function that employs only a first-order term directly proportional to the partial pressure of CO 2 . Note that such a linear dependence on P CO 2 is expected for the following dissolution mechanism in an idealized situation (CO 2 , fluid + O 2 -melt ⇆ CO 2 -3 , melt ; Pan et al. 1991):  \nw CO 2 = 0 . 0054 P CO 2 [GPa] (13)  \nwhere w CO 2 is the weight fraction of CO 2 . This equation shows a below 5% difference from the original equation within the pressure range of 200 to 1700 bar. We rewrite this equation to:  \nX CO 2 = 0 . 0054 ( µ m µ CO 2 )( P CO 2 1 GPa ) . (14)  \nThe CO dissolution to MORB satisfies the following equation at temperatures between 1523 and 1873 K and pressures between 2000 and 30000 bar (Yoshioka et al. 2019):  \nw CO [%] = 10 -4 . 83+0 . 8 log f CO[bar] , (15)  \nor  \nX CO = 0 . 0002 · ( ϕ CO ) 0 . 8 · ( µ m µ CO )( P CO 1 GPa ) 0 . 8 . (16)  \nThe CH 4 solubility law refers to the solubility to haplo-basalt given in Ardia et al. (2013) with the experimental data with pressures within 7000 to 30000 bar and temperatures between 1673 and 1723 K. We omit the temperature dependence of CH 4 solubility due to their limited temperature range (e.g., Hirschmann et al. 2012). The derived solubility law is as follows:  \nX CH 4 = f CH 4 [GPa] e -7 . 63 -1 . 9P[GPa] , (17)  \nor  \nX CH 4 =0 . 0005 · s CH 4 ( P CH 4 1 GPa ) (18)  \ns CH 4 = ϕ CH 4 (0 . 15) P [GPa] (19)  \nAt this point, it is worth estimating the impact of dissolution into magma on the atmospheric composition by approximately calculating the fraction of dissolved species in the magma layer relative to those in the atmosphere at the BOA. For this, we use the following approximate relations for the total number of molecules in the atmosphere, N atm :  \nN atm ∼ 4 πR 2 core P BOA µ atm g (20)  \n∼ 4 πR 2 core P BOA µ atm R 2 core GM core (21)  \n∼ 4 πR 4 core P BOA GM p µ atm (22)  \nwhere P BOA is the pressure at BOA and µ atm is the mean molecular weight of the atmosphere. Here, we ignore the dependence of the gravity g on the elevation above the magma-atmosphere interface. Meanwhile, the total number of molecules in magma, N m , can be written as  \nN m = M core ζ µ m ∼ M p ζ µ m , (23)  \nwhere ζ is the mass fraction of the silicate layer that is molten and reactive. Denoting the mixing ratio of species S in the atmosphere is indicated by x S , the ratio of the species S dissolved in magma to those in the atmosphere, r S , can be expressed as follows:  \nr S = N m X S  \nN atm x S (24) = GM 2 p ζ 4 πR 4 core µ atm µ m X S x S P BOA (25)  \n∼ 42 ( M p 6 M ⊕ ) 2 ( R p 1 . 6 R ⊕ ) -4 ( ζ 0 . 66 )( µ atm 0 . 1 µ m ) X S P BOA ,S [GPa] , (26)  \nwhere P BOA , S is the partial pressure of the species S at BOA. Substitution of Equations (16) and (18) indicates that r CH 4 and r CO range from about 0.01 to 0.1 if the atmospheric pressure is lower than a few GPa. This implies that the reduced C-bearing species (CO and CH 4 ) will remain in the atmosphere in most cases (neglecting uptake by the iron core).  \nThe solubility of CO 2 could marginally affect the result:  \nr CO2 =0 . 23 ( M p 6 M ⊕ ) 2 ( R p 1 . 6 R ⊕ ) -4 ( ζ 0 . 66 )( µ atm 0 . 1 µ CO 2 ) . (27)  \nThus, only in the case where the atmosphere is highly oxidized ( µ atm ≳ 0 . 1 µ CO 2 ), and the substantial fraction of silicate is almost fully molten and well-mixed, i.e., ζ ∼ 0 . 66, atmospheric CO 2 is somewhat affected by the dissolution. We do not set an upper limit on the dissolved mole fraction of volatiles in magma. Within our parameter range, the maximum mole fraction of volatiles in the magma reaches about 29%, corresponding to a mass fraction of about 4.0%. Given the small mass fraction of the volatile species in magma, we neglect the density decrease and compositional changes in magma due to dissolved volatiles, as well as their effects on chemical reactions and the mass-radius relation of the core.", '2.4. Numerical procedure': 'We calculate the equilibrium atmospheric compositions of the atmosphere-magma system in terms of five species (H 2 , H 2 O, CH 4 , CO, CO 2 ). The elemental composition of the total accreted volatile species matches the solar system composition from Lodders (2003); C/H= 2 . 9 × 10 -4 and C/O= 0 . 5. This implicitly assumes that nebula gas accretion took place interior to the snowline of H 2 O and C-bearing species, which would be consistent with the assumption of the rocky interior (e.g., Bean et al. 2021). The five species both in the atmosphere and the magma (2 × 5 = 10 variables) and the amount of O bounded to Fe in the rock (1 variable) can be determined by the three redox equations, the five dissolution laws, and the total amount of H, O, and C in the atmosphere-magma system (three constraints); see also Schlichting & Young (2022).  \nThe full set of equations is inherently non-linear. The non-linearity would disappear when we assume that (1) the atmospheric amount of CO, CO 2 , and CH 4 is negligible, and (2) the non-linear solubility law of CO is negligible, both of which are valid in the case studied here where the C/H in the accreted volatile is small. However, for future applications to general cases, we solve the equations using a scipy Nelder-Mead root-finding algorithm (Virtanen et al. 2020) in the following way.  \n- 1. We begin with estimating the radius of the planetary core R core and the number of molecules of the magma layer N m , using the mass-radius relationship of bare rock by Fortney et al. (2007) with the given mass fraction of the silicate layer (default: 0.66) and the externally given planetary mass, M p . Note that we ignore the mass of the atmosphere since the mass fraction of the atmosphere for sub-Neptunes considered here is at most 0.46%. The difference between the planetary radius R p and R core gives the thickness of the atmosphere, ∆ R .  \n- 2. We then prescribe a value for the atmospheric H 2 O fraction x H 2 O . The much smaller amount of C compared to H ( ∼ 10 -4 ) and the fixed atmospheric He H ratio allow us to derive the atmospheric H 2 fraction x H 2 and the He fraction x He from x H 2 O . The temperature and pressure at the BOA are either specified directly (Section 3.1) or calculated using the adiabatic lapse rate for the given atmospheric composition and planetary parameters (Section 3.2).\n- 3. The oxygen fugacity ( f O 2 ) and rock composition (the number of pure Fe in magma, [FeO], [Fe 2 O 3 ]) are determined based on the redox equilibrium given the temperature and pressure at BOA as well as the x H 2 O x H 2 .\n- 4. From the composition of the accreted volatiles (with the abundance ratio of H, O, and C being the solar values), surface temperature ( T BOA ), surface pressure ( P BOA ), f O 2 , atmospheric He/H ratio, N m , and rock composition, we determine the amounts of H 2 , H 2 O, CO, CO 2 , He, and CH 4 in the atmosphere and magma under the constraints on the O/H and C/H of the accreted gas. We employ a root-finding algorithm on the atmospheric C/H ratio due to the non-linearity. If the atmospheric C/H ratio is given, atmospheric amounts of C-bearing species ( x CO , x CO 2 , and x CH 4 ) can be derived from x H 2 , x H 2 O , and surface conditions ( T BOA , P BOA , and f O 2 ). Then, solubility laws and N m give the dissolved amounts of all volatiles, and planetary amounts of all volatiles naturally emerge. The root-finding algorithm finds the atmospheric C/H ratio that makes the planetary volatile composition consistent with the assumed nebula composition: C/H= 2 . 9 × 10 -4 . To be specific, the difference between the calculated planetary C/H and the assumed nebula C/H is minimized until it becomes smaller than 10 -5 .\n- 5. The equilibrium calculation gives the total amount of O ( N O ). The difference between N O and the amount of O from accreted gas is interpreted as the amount of N O(Fe) before the chemical reaction, which is also related to the redox state of magma before the gas accretion.  \nBy repeating this process with varying x H2O , we obtain the relation among x H2O , the mixing ratios of C-bearing species, and the rock composition (both before and after the gas accretion).', '2.5. Consideration of magma solidification': 'When considering magma solidification, we assume a scenario where the atmosphere first reacts with the initially hot, fully molten homogeneous silicate core, and the chemical equilibrium changes as the magma solidifies from the bottom over time. In this study, we simplify the situation by considering only two discrete stages rather than the continuous solidification and the change in equilibrium state: (1) the initial fully molten magma reacting with the atmosphere and (2) the completion of solidification (i.e., the desired solidified fraction has been reached). We also assume that there is no exchange of material with the outside of the magma-atmosphere system and that the surface temperature remains constant. The concrete calculation method is as follows:  \n- 1. We calculate the chemical equilibrium state of a planet with fully molten homogeneous magma by assuming specific values for the atmospheric H 2 O fraction x H 2 O , T BOA , and P BOA . From this result, we determine the planetary amounts of H, C, He, and O (Section 2.4).\n- 2. We find the new equilibrium state under the updated molten fraction. As mentioned in Section 2.1, we assume that no volatiles can exist in solidified magma, so the amounts of H, C, and He in the magma and atmosphere remain the same as in the initial state. However, some of the O is excluded from the reaction in the form of FeO and Fe 2 O 3 contained in the solidified magma. The excluded amount can be calculated from the iron speciation of the magma obtained in step 1. Again, a root-finding algorithm is applied to determine the P BOA , atmospheric C/H ratio, atmospheric He/H ratio, and x H 2 O that satisfy the given amounts of H, C, He, and O.', '2.6. Parameter range': 'In order for some of our assumptions to be valid, we consider the results only when the following criteria are met:  \n- · T BOA > 1800 K: This ensures that the silicate layer is at least partially molten.\n- · P BOA < 30000 bar: This range is determined based on the pressure range of the solubility experimental study of the most dominant gas species in our results, H 2 . With our fiducial set of planetary parameters ( M p = 6 M ⊕ ,  \nFigure 4. The relation between the mole fraction of each species and N before O(Fe) / N Fe . The variable for the x-axis, N before O(Fe) / N Fe , represents the redox state of the silicate before the reaction. The upper panel corresponds to the atmospheric composition, while the lower panel denotes the species present in the magma. The surface temperature and pressure are fixed at 3000 K and 10000 bar, respectively.  \n<!-- image -->  \nT eq = 750 K) and assuming an H 2 -rich atmosphere ( x H 2 O = 0 . 01), 30000 bar corresponds to the atmospheric mass fraction of 0 . 5% and the atmospheric thickness of ∆ R ∼ 0 . 8 R ⊕ , corresponding to R p = 2 . 4 R ⊕ . This provides an idea of the upper limit of the planetary radius discussed in this paper. Due to the lack of available data, we do not explicitly discuss larger planets such as K2-18 b, although we expect qualitatively similar trends.', '3.1. Dependence on surface conditions and magma redox state': "This section presents the relation between atmospheric composition and magma redox state before the reaction (as represented by the iron speciation of the magma), with a specific focus on the influence of two fundamental parameters in determining the equilibrium state: surface pressure and temperature. Here, the rocky layer is assumed to be fully molten.  \nBefore discussing elemental ratios, we first introduce how volatile speciation depends on the magma redox state. Figure 4 shows the volatile speciation in the planetary model as a function of the magma redox state before reactions, N before O(Fe) / N Fe . The upper and lower panels represent the atmosphere and magma, respectively. Surface pressure ( P BOA ) and temperature ( T BOA ) are fixed at 10000 bar and 3000 K. It is clear that as the magma becomes more oxidized, the fractions of H 2 O, CO, and CO 2 increase relative to H 2 and CH 4 . Additionally, the larger solubility of H-bearing species, particularly H 2 O, results in their relative mole fractions in magma being substantially larger than those of C-bearing species.  \nFigure 5 shows the atmospheric O/H and C/H as a function of N before O(Fe) / N Fe . The upper and lower panels represent the dependence on surface pressure ( P BOA ) and temperature ( T BOA ), respectively.  \nThe atmospheric O/H ratio is primarily governed by the redox state of the silicate before the reaction, with its dependence on BOA pressure becoming non-negligible as pressure increases. When the pressure exceeds approximately  \nFigure 5. The relation between the atmospheric H-C-O ratio and N before O(Fe) / N Fe , with its dependency on surface pressure and temperature. The solid line corresponds to the O/H ratio, while the dashed line represents the C/H ratio. We note that log(O / H) and log(C / H) of accreted volatiles are assumed to be about -3 . 23 and -3 . 53. (Upper) We show four surface pressure cases: 1000, 3000, 10000, and 30000 bar. The surface temperature is fixed at 3000 K. (Lower) Three surface temperature cases are shown: 2000, 3000, and 4000 K. The surface pressure is fixed at 10000 bar.  \n<!-- image -->  \n30000 bar, the atmospheric O/H ratio decreases to about half of what it is at 1000 bar under the condition where N before O(Fe) / N Fe ∼ 1. Overall, the atmospheric H 2 O fraction ( x H 2 O ∼ 2 O H ) ranges from a few percentages to a few tens of percentages (for thin atmospheres with P BOA ≲ 3000 bar) unless the magma before the reaction is highly reduced, where N before O(Fe) / N Fe ≲ 0 . 1. On the other hand, the atmospheric C/H ratio appears relatively insensitive to N before O(Fe) / N Fe and surface pressure. The atmospheric C/H ratio maintains the nebula gas composition except in the lowest-pressure ( < 1000 bar) and most oxidized atmospheres with an H 2 O mixing ratio exceeding a few percent, leading to an enrichment of C/H several dozen times. The impact of temperature on both O/H and C/H ratios appears to be relatively limited, at least for temperatures below 4000 K.  \nThis can be understood from an approximate analytic expression for O/H and C/H, described below. We first focus on the O/H ratio. We begin by writing a simple mass conservation of O:  \nN before O(Fe) = N after O(Fe) + N atm x H 2 O + N m X H 2 O , (28)  \nwhere N atm and N m represent the total number of molecules in the atmosphere and magma, respectively. However, due to the large solubility of H 2 O in magma, the second term on the right-hand is substantially smaller than the third term if ζ 0 . 66 ≳ 0 . 5 under our fiducial planetary parameter combination (See Equation 37). Therefore, for simplicity, the second term is omitted below.  \nUsing Equation (10) and assuming that the Iron-Wustite buffer dominates the reactions after the reaction ([FeO] after ≫ [Fe 2 O 3 ] after ),  \nN before O(Fe) = N m { [FeO] after +0 . 16 ( µ m µ H 2 O )( P H 2 O 10000 bar ) 0 . 74 } (29)  \nTherefore, Equation (29) becomes  \nN before O(Fe) = N m ( a 1 x H 2 O + a 2 x 0 . 74 H 2 O ) (32)  \nN before O(Fe) N m = N before O(Fe) N Fe f Fe = a 1 x H 2 O + a 2 x 0 . 74 H 2 O (33)  \na 1 ≡ 1 . 2 ϕ H 2 O γ FeO ϕ H 2 K (34)  \na 2 =0 . 54 ( P BOA 10000 bar ) 0 . 74 , (35)  \nwhere f Fe is the mole fraction of Fe in magma. We note that 1 x H 2 ∼ 1 . 2 under the condition that H 2 and He dominate the atmospheric composition, and the He/H ratio is 0.1. The first term makes the relation N before O(Fe) ∝ x H 2 O . The pressure and temperature dependence of the part of the proportionality constant, K ' = ϕ H 2 O γ FeO ϕ H 2 K , is shown in Figure 6. While K ' increases with increasing temperature and decreases with increasing pressure, K ' does not vary significantly. On the other hand, the second term leads to N before O(Fe) ∝ x 0 . 74 H 2 O . The proportionality constant for this term is influenced solely by surface pressure, with higher pressure resulting in a larger constant. These assessments suggest that the second term becomes dominant for massive ( P BOA ≳ 10000 bar) and reduced atmospheres.  \nAlthough these two terms have opposite trends in response to pressure change (i.e., an increase in pressure leads to a decrease in the first term and an increase in the second term), their sum simply increases as the pressure increases for a given x H 2 O , as indicated by the upper panel of Figure 5. This is because the pressure dependence of K ' is minimal under low pressure, where the first term dominates. This result is consistent with the intuitive understanding that as the atmospheric amount increases, x H 2 O decreases due to the limited oxidizing power of magma.  \nFigure 6. The temperature- and pressure-dependence of the K ' = ϕ H 2 O γ FeO ϕ H 2 K .  \n<!-- image -->  \nwhere µ m µ H O ∼ 3 . 4 from the assumed magma composition.  \n2 Denoting the equilibrium constant of the following reaction: H 2 + FeO ⇆ H 2 O+ Fe, we may write  \nK = f H 2 O[Fe] γ FeO f H 2 [FeO] (30)  \n= ϕ H 2 O x H 2 O γ FeO ϕ H 2 x H 2 [FeO] . (31)  \nThe temperature dependence of K ' results in the trend that the same N O(Fe) corresponds to a more reduced atmosphere under higher temperatures. However, due to the limited variation in the value of K ' mentioned above, the overall dependence is minor, as depicted in the bottom panel of Figure 5.  \nFigure 5 shows that the slope of the O/H ratioN before O(Fe) relationship changes towards the most oxidized end at low pressure (1000 and 3000 bar). This abrupt change corresponds to the depletion of pure iron in magma, where the reaction transitions from Fe ←→ FeO ←→ Fe 2 O 3 to FeO ←→ Fe 2 O 3 after the depletion of pure iron. This corresponds to the change in the slope of the relation between f O 2 and N O(Fe) , shown in Figure 2. When the atmospheric pressure is higher than ∼ 3000 bar, the reducing power of the atmosphere is enough to produce pure iron, even if the rock before the reaction is fully oxidized. Therefore, a sudden change in the slope of the O/H ratioN before O(Fe) relationship is not observed.  \nNext, we will discuss the dependence of C/H on planetary parameters. As expected, the atmospheric C/H ratio is always larger than the input (total) C/H value (10 -3 . 53 ) due to the dissolution of H-bearing species into magma. Combining Equation (26) and the solubility laws, the ratio of the dissolved H 2 and H 2 Oto the atmospheric counterparts are:  \nr H 2 =4 . 6 ( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( ζ 0 . 66 )( µ atm 0 . 1 µ m ) ϕ H 2 (0 . 47) P BOA [GPa] (36)  \nr H 2 O =22 ( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( ζ 0 . 66 ) ( x H 2 O 10 -2 ) -0 . 26 ( µ atm 0 . 1 µ H 2 O )( P BOA 10000 bar ) -0 . 26 . (37)  \nIn Equation (36), the term ϕ H 2 (0 . 47) P BOA [GPa] is close to 1 when the pressure is less than a few thousand bars. At high pressures, it exhibits temperature dependence and ranges from ∼ 1 . 5 (2000 K) to ∼ 0 . 4 (4000 K) at 30000 bar. However, we note that the range of ϕ H 2 (0 . 47) P BOA [GPa] is relatively narrow and does not affect the result significantly. If we ignore the He fraction in atmosphere, the conservation law of H can be expressed as follows:  \nN H = { (1 + r H 2 )(1 -x H 2 O ) + (1 + r H 2 O ) x H 2 O } N atm , H (38)  \nwhere N atm , H is the total number of H in the atmosphere. As previously discussed in Section 2.3, the solubility of C-bearing species, except for the highly oxidized atmosphere case, is negligible. Ignoring the solubility of C-bearing species yields:  \n( C H ) atm ∼ N C N atm , H (39)  \n∼{ (1 + r H 2 )(1 -x H 2 O ) + (1 + r H 2 O ) x H 2 O } ( C H ) total . (40)  \nWhen x H 2 O ≪ 1, this reduces to  \n( C H ) atm ∼ (1 + r H 2 ) ( C H ) total . (41)  \nWhen combined with Equation (36) and µ m ∼ 0 . 06 kg / mol, it becomes apparent that the dissolution of H 2 can increase ( C H ) atm approximately 1.6 to 8.0 times, depending on the surface conditions (see also Chachan & Stevenson 2018).  \nOn the other hand, if the second term is dominant,  \n( C H ) atm ∼ r H 2 O x H 2 O ( C H ) total (42)  \n∼ 0 . 22 ( M p 6 M ⊕ ) 2 ( R core 1 . 6 M ⊕ ) -4 ( ζ 0 . 66 )( P BOA 10000 bar ) -0 . 26 ( µ atm 0 . 1 µ H 2 O ) ( x H 2 O 10 -2 ) 0 . 74 ( C H ) total . (43)  \n( C H ) atm increases as x 0 . 74 H 2 O where the power-law index comes from the solubility law of H 2 O.  \nThe condition when the second term becomes dominant can be estimated by taking the ratio between the two terms on the right-hand side of Equation (40), approximating that r H 2 ≫ 1 and r H 2 O ≫ 1:  \n2ndterm 1st term = (1 + r H 2 O ) x H 2 O (1 + r H 2 )(1 -x H 2 O )  \n∼ 16 1 0 . 46 P (GPa) ϕ H 2 ( x H 2 O 1 -x H 2 O )( P BOA 10000 bar ) ( x H 2 O 10 -2 ) -0 . 26  \n∼ 0 . 16 1 0 . 46 P (GPa) ϕ ( P BOA 10000 bar ) ( x H 2 O 10 -2 ) 0 . 74 .  \n(44) -0 . 26 (45) H 2 -0 . 26 (46)  \nIf we let 0 . 46 P (GPa) ϕ H 2 ∼ 1, the second term becomes dominant when  \n( P BOA 10000 bar ) -0 . 26 ( x H 2 O 10 -2 ) 0 . 74 ≳ 60 . (47)  \nTherefore, the C/H ratio can increase significantly, by more than several tens of times, in a low-pressure, oxidizing atmosphere case. This effect is illustrated by the red line in the upper panel of Figure 5, which represents the 1000 bar surface pressure case, showing the highly enriched C/H ratio within the oxidized region that the Ferric-Ferrous governs.  \nEquations (36), (41), and (43) suggest that the atmospheric C/O ratio would decrease as the atmosphere becomes more H 2 O-rich, as the increase of the atmospheric C/H ratio is not as significant as that of H 2 O. If H 2 dissolution dominates the change in atmospheric C/H ratio and we assume the amounts of C-bearing species and He is negligible, the atmospheric C/O ratio can be approximated as follows, based on Equation (41):  \n( C O ) atm ∼ { 1 + 4 . 6 ϕ H 2 (0 . 47) P BOA [GPa] ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( µ atm 0 . 1 µ m ) } ( O H ) -1 atm ( C H ) total ∼ { 200 + 920 ϕ H 2 (0 . 47) P BOA [GPa] ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( µ atm 0 . 1 µ m ) } ( x H 2 O 10 -2 ) -1 ( O H ) accreted ( C O ) accreted ∼ { 0 . 12 + 0 . 53 ϕ H 2 (0 . 47) P BOA [GPa] ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( µ atm 0 . 1 µ m ) } ( x H 2 O 10 -2 ) -1 ( C O ) accreted (48)  \nwhere ( O H ) accreted is the O/H ratio of accreted volatiles, which has a value of 5 . 8 × 10 -4 in this study. It is evident that the C/O ratio follows a simple inverse proportional relation with x H 2 O . If x H 2 O exceeds a few percent, the C/O ratio diminishes to 10% of that of the assumed composition of accreted gas.  \nOn the other hand, if the H 2 O dissolution governs the change in atmospheric C/H, the atmospheric C/O follows the relation as follows:  \n( C O ) atm ∼ 0 . 22 ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( P BOA 10000 bar ) -0 . 26 ( µ atm 0 . 1 µ H 2 O ) ( x H 2 O 10 -2 ) 0 . 74 ( O H ) -1 atm ( C H ) total ∼ 44 ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( P BOA 10000 bar ) -0 . 26 ( µ atm 0 . 1 µ H 2 O ) ( x H 2 O 10 -2 ) -0 . 26 ( O H ) accreted ( C O ) accreted ∼ 0 . 026 ( ζ 0 . 66 )( M p 6 M ⊕ ) 2 ( R core 1 . 6 R ⊕ ) -4 ( P BOA 10000 bar ) -0 . 26 ( µ atm 0 . 1 µ H 2 O ) ( x H 2 O 10 -2 ) -0 . 26 ( C O ) accreted . (49)  \nEquation (49) demonstrates that when the dissolution of H 2 O plays a dominant role in C/H enrichment, the relationship between C/O depletion and x H 2 O is milder ( ∝ x -1 H 2 O → x -0 . 26 H 2 O ).  \nFigure 7 shows the O/H-C/O relationship and its dependence on surface pressure. As discussed earlier using the two C/O equations, the decrease in C/O begins with a monotonic trend, becomes gentler as x H 2 O increases, and finally reaches ∼ 10% of the nebula gas value under an oxidized atmosphere with x H 2 O ≥ a few percentages. We note that the C/O ratio decreases with increasing pressure due to the lower r H2 and r H2O , as Equation (36) and Equation (37) show.  \nFigure 7. The relationship between atmospheric O/H and C/O ratios and their dependence on surface pressure. Different colors represent different surface pressures: 1000, 3000, 10000, and 30000 bar. The blue vertical and horizontal lines denote the solar system values of O/H and C/O. A dotted horizontal line marks 10% of the solar C/O.  \n<!-- image -->", '3.2. Dependence on Fe/Si and size of magma': 'In this section, we discuss two additional magma properties that affect the atmospheric composition beyond the magma redox state (iron speciation) already covered in Section 3.1: the magma composition (Fe/Si ratio) and the reactive amount of magma (reactive mass fraction to the total magma mass, ζ 0 . 66 ). Our model assumes fiducial parameters of a Fe/Si ratio of 0.17 (MORB-like) and ζ 0 . 66 of 1. The ranges considered in this section vary from 10 -1 to 10 0 . 5 times the MORB value for Fe/Si and from 10 -2 to 1 for ζ 0 . 66 . We assume that the T BOA is 3000 K and the P BOA is 10000 bar.  \nFigure 8 illustrates the relationship between the atmospheric O/H and C/H ratios and the magma iron speciation (redox state), focusing on the dependence on magma properties. The upper panel shows that iron-poor magma corresponds to a smaller O/H ratio for the same iron speciation, which is expected due to the limited amount of oxidized Fe species. This result is also consistent with the atmospheric composition trends with respect to different FeO content (e.g., Kite et al. 2020).  \nThe lower panel demonstrates that a smaller effective magma fraction results in a lower O/H ratio when ζ 0 . 66 < 0 . 1. This can be understood based on the discussion in Section 3.1. As long as we can ignore the second term of Equation (28), the ζ -dependence of O/H is minor. When ζ becomes smaller, however, the second term can no longer be ignored, and Equation (33) becomes  \nN before O(Fe) N Fe f Fe = a 1 x H 2 O + N atm N m x H 2 O + a 2 x 0 . 74 H 2 O . (50)  \nNow, for a given N before O(Fe) / N Fe , the smaller N m is compensated by a lower x H2O .  \nNote that the atmospheric O/H-C/O relation remains consistent across different Fe/Si cases. However, a smaller reactive magma amount causes the same O/H ratio to correspond to a lower C/O ratio, as can also be inferred from Equation (48) and Equation (49).  \nFigure 8. The relation between the atmospheric H-C-O ratio and N before O(Fe) / N Fe , with its dependency on two magma properties, Fe/Si ratio and size of reactive magma ( ζ 0 . 66 ). (Upper) We show three Fe/Si cases: 10 -1 , 1, and 10 0 . 5 times the MORB value (0.17). The reactive magma fraction is fixed at 1. (Lower) Three reactive magma fraction cases are shown: 10 -2 , 10 -1 , 1. The Fe/Si is fixed at the MORB value. The surface temperature and pressure are fixed at 3000 K and 10000 bar, respectively.  \n<!-- image -->', '3.3. Dependence on atmospheric parameters': 'In our model, the two key factors that control the atmospheric composition in contact with magma ocean, P BOA and T BOA , are related to the observable parameters R p , T eq , and M p through the assumption on P RCB and Earth-like rocky core. In this section, we show the atmospheric composition as a function of these parameters. The fiducial parameters are M p = 6 M ⊕ , T eq = 750 K, and P RCB = 10 bar. Instead of R p , we use the atmospheric thickness, ∆ R = R p -R core , to represent the planetary radius. The fiducial value for ∆ R is 0 . 4 R ⊕ . Note that under the assumption of an Earth-like rocky core, the core radius is found to be R core = 1 . 6 R ⊕ , and thus ∆ R = 0 . 4 R ⊕ corresponds to R p = 2 . 0 R ⊕ . We change one of these parameters at a time to clarify the dependence of each parameter. The ranges of parameters considered are 0.1 to 0.8 R ⊕ (∆R), 550 to 1150 K ( T eq ), 0.1 to 100 bar ( P RCB ), and 2 to 8 M ⊕ ( M p ).  \nFigure 9 illustrates the ranges of P BOA and T BOA under an oxidized atmosphere and their dependencies on the planetary parameters. A thicker atmosphere (larger ∆ R ), a colder upper atmosphere (smaller T eq ), and a heavier planetary mass (larger M p ) correspond to larger P BOA and T BOA . On the other hand, the effects of P RCB are different for pressure and temperature: the larger P RCB corresponds to higher P BOA and lower T BOA .  \nWe show the relation between the atmospheric composition (O/H and C/H) and the magma redox state and its dependence on the planetary parameters in Figure 10. It demonstrates that atmospheres with smaller ∆ R , M p , and P RCB and larger T eq result in larger O/H and correspondingly elevated C/H ratios, consistent with the expectations from the P BOA result shown in Figure 9. From these investigations, we find that the effect of the planetary parameters, except for ∆ R and P RCB , is almost negligible. However, O/H is sensitive to small ∆ R and small P RCB since both can make a small atmospheric pressure.  \nFigure 9. The relationship among P BOA (red contour), T BOA (black contour), ∆ R , T eq , P RCB , and M planet in our model. The atmosphere is composed of H 2 (82%), He (17%), and H 2 O (1%). In the left panel, P RCB and M planet are fixed at 10 bar and 6 M ⊕ , respectively. In the right panel, ∆ R and T eq are fixed at 0 . 4 R ⊕ and 750 K, respectively.  \n<!-- image -->  \nIn this section, we demonstrate the effect of magma solidification on the atmospheric composition following the methodology described in Section 2.5. We assume that the P BOA before the solidification is 5000 bar and the T BOA remains constant at 3000 K. Note that fixing the BOA temperature does not cause significant error in estimating the atmospheric composition (see Section 3.1).  \nFigure 11 shows the case with ζ 0 . 66 = 10 -1 , 10 -2 , and 10 -3 scenarios. It is clear that the x H2O increases under the same N before O(Fe) / N Fe due to the extraction of dissolved H 2 O from the solidified rock layer to the upper molten magma and atmosphere. At the same time, the C/H ratio decreases due to the decrease of ζ . These two effects imply that the C/O ratio can decrease further with the magma solidification effect.', '4.1. Implications for the internal structure and formation scenarios': 'We have examined the range of atmospheric composition expected for the scenario where the rocky planetary core accreted the nebula gas. Can such atmospheric composition be a useful clue to the formation pathways of subNeptunes? To address this question, we investigate the expected volatile composition of sub-Neptunes under an alternative formation scenario: the assembly of icy planetesimals outside of the snow line (referred to as ice-rich sub-Neptunes hereafter) and compare it to the case of magma-bearing sub-Neptunes.  \nThe composition of ice-rich sub-Neptunes would be determined from the elemental composition of the solid components in the protoplanetary disk (e.g., Oberg et al. 2011). The elemental abundance ratio of gaseous and solid components within the protoplanetary disk is expected to vary as a function of the orbital radius relative to the sublimation lines of volatile species, although there are ongoing debates on the underlying mechanisms and chemistry of the protoplanetary disk (e.g., Oberg et al. 2023). For example, the C/O ratio of solids (gas) beyond the H 2 O snow line would be smaller (larger) than that interior of the H 2 O snow line. The exact C/O value in solids at a given location would depend on the ratio of carbon in volatile form to refractory form. Under the typical scenario where about 75%  \nFigure 10. The relation between the atmospheric H-C-O ratio and N before O(Fe) / N Fe , with respect to four planetary parameters: (∆ R , T eq , P RCB , and M planet ). The x-axis denotes N before O(Fe) / N Fe , while the y-axis represents atmospheric composition. Solid lines represent O/H ratio, while dashed lines represent C/H ratio. Horizontal solid and dashed lines indicate log(O / H) and log(C / H) of accreted volatiles, which are approximately -3 . 23 and -3 . 53, respectively. (Upper left) Four different ∆ R cases: 0.1, 0.4, 0.6, and 0.8 ( R ⊕ ). (Upper right) Four P RCB cases: 0.1, 1, 10, and 100 bar. (Lower left) Four T eq results: 550, 750, 950, and 1150 K. (Lower right) Four M planet lines: 2, 4, 6, and 8 M ⊕ . Cyan lines serve as the fiducial values: ∆ R = 0 . 4 R ⊕ , P RCB = 10 bar, T eq = 750 K, and M planet = 6 M ⊕ . Data for cases where T BOA and P BOA are out of the parameter range of this study are not shown.  \n<!-- image -->  \nof the carbon exists in volatile form and a rock-to-water ratio of approximately 1, the C/O in solids beyond the H 2 O snow line would be a few tens of percent of the nebula value. If the planetesimals beyond the CO 2 and CO snow lines also accrete, the planetary C/O would further increase. In these cases, we anticipate that the oxidized atmosphere in contact with magma ocean, which exhibits a depleted C/O ratio well below 10% of nebula value, would be distinct from the expectation for the ice-rich sub-Neptune scenario.  \nHowever, uncertainties exist in the previously described simplistic scenario. In particular, the fraction of refractory carbon could be smaller due to the irreversible destruction processes (e.g., Li et al. 2021), which would reduce C/O in ice-rich sub-Neptune scenarios to an unknown extent. This uncertainty might overlap with the C/O expected in the magma scenario considered in this paper. Further studies on the plausible range of the volatile compositions of planetesimals beyond the H 2 O snow line and the formation pathways of small planets (e.g., Kimura & Ikoma 2022) are needed. These studies will help to clarify further the atmospheric compositions that can distinguish the formation pathways of sub-Neptunes, and this is left for future work.', '4.2. Implications for observation': 'In this section, we discuss the spectral features of the atmospheric scenarios we find. Using the open-source code TauREx3 (Al-Refaie et al. 2021), we simulate the transmission spectra for two representative atmospheric scenarios from our main results (see Figure 7): (1) the reduced case with log 10 O / H = -3 . 5 and log 10 C / O = 0 . 3, and (2) the oxidized case with log 10 O / H = -1 . 5 and log 10 C / O = -1 . 3.  \nThe atmospheric temperature-pressure structure is calculated with the nominal assumption: ∆ R = 0 . 4 R ⊕ , T eq = 750 K, P RCB = 10 bar, and M p = 6 M ⊕ . To compute the molecular composition of the atmosphere, we assume that molecular composition is constant above the quenching point, where the vertical mixing timescale overwhelms  \nFigure 11. The relation between the atmospheric H-C-O ratio and silicate redox state (before the reaction), with its dependency on the molten fraction of silicate core. The variable for the x-axis is N before O(Fe) / N Fe , which denotes the redox state of the silicate before the reaction. The y-axis corresponds to the elemental ratio on a logarithmic scale. The solid line corresponds to the O/H ratio, while the dashed line represents the C/H ratio. log(O / H) and log(C / H) of accreted volatiles are approximately -3 . 23 and -3 . 53, respectively. The lines of different colors represent various mass fractions of magma ( ζ ), with ζ = 6 . 6 × 10 -1 depicted in black, ζ = 6 . 6 × 10 -2 in blue, and ζ = 6 . 6 × 10 -3 in magenta.  \n<!-- image -->  \nthe chemical timescale (e.g., Hu & Seager 2014). We approximately represent this quench point by a single pressure. Noting that, in reality, the quenching point can depend on the molecules and the background atmospheric composition, we explore the dependence on the quenching pressure by varying it between 100 and 10 bar. We also tested the case where the quenching pressure is located at 1 bar, but it does not differ significantly from the 10 bar case. In all of these cases, we find that O and C in the atmosphere exist predominantly as H 2 O and CH 4 respectively (see Figure 11 of Hu & Seager 2014). Therefore, the mixing ratio of H 2 O and CH 4 is simply ∼ 2 log(O / H) and ∼ log(C / H), respectively.  \nFigure 12 illustrates the transmission spectra with a spectral resolution of 100 for the two atmospheric scenarios we mentioned above. The left and right columns correspond to two different quenching points. We assumed that the parent star is an M-type star with a temperature of 3600 K and a stellar radius of R ⋆ = 0 . 42 R ⊙ . If one wants to consider the transit spectra for the planets accompanied by other stars, the transit depth can be scaled by R -2 ⋆ . All spectra are primarily shaped by the absorption bands of H 2 Oand CH 4 . The dominant features for the reduced scenario come from CH 4 (1.65 µ m, 2.4 µ m, 3.3 µ m, and 6.15 µ m), while the H 2 O features are also prominent (1.9 µ m and 2.87 µ m) in the oxidized case. Although the absorption bands of H 2 O and CH 4 tend to overlap at shorter wavelengths, they diverge at longer wavelengths, in particular beyond 2 µ m, enabling the detection of both species with high S/N observations. Notably, the oxidized atmosphere exhibits a weak signal of CO 2 at 4.3 µ m in cases where the quenching point is as deep as 100 bar, providing additional constraints on O/H. Through the strong spectral features of CH 4 and H 2 O, we anticipate that finding the atmospheric H-C-O ratio from the spectral observation would be straightforward. Additionally, the mean molecular weight, another observable parameter (e.g., Gao et al. 2023; Benneke et al. 2024), increases from 2.4 g/mol to 3.2 g/mol as the O/H ratio shifts from reduced to oxidized states.  \nNote that the signal strengths of the molecule spectral peaks in our demonstration (Figure 12) could be overestimated compared to transmission spectra shown in previous studies (e.g., Kawashima & Ikoma 2019). This discrepancy arises because we did not consider photolysis caused by stellar UV irradiation, which can significantly reduce the abundance of  \nFigure 12. Simulated transmission spectra of two representative atmospheric cases from our results (reduced, C-enriched one and oxidized, C-poor one) with varying quenching pressures. The top row corresponds to the reduced atmosphere (atmospheric O/H ratio of 10 -3 . 5 and C/O ratio of 10 0 . 3 ), while the bottom row illustrates the oxidized atmosphere (atmospheric O/H ratio of 10 -1 . 5 and C/O ratio of 10 -1 . 3 ). Different colored lines represent the contributions of different species: H 2 O (blue), CO (cyan), CO 2 (green), and CH 4 (magenta). The left column shows the case with a quenching point at 100 bar, while the right column denotes the shallower quenching point at 10 bar.  \n<!-- image -->  \nmolecules in low-pressure regions. Such inhomogeneity is prominent at altitudes where the timescale of photochemistry is shorter than the mixing timescale; for instance, Kawashima & Ikoma (2019) shows that H 2 O and CH 4 can be fully photo-dissociated at pressures below 10 -4 . 5 bar under intense UV irradiation (10 5 times that of GJ 1214). In our simulations, for the reduced and carbon-enriched case (upper left panel of Figure 12), we confirmed that complete photolysis under 10 -4 . 5 bar would set the maximum transit depth to be ∼ 2 . 2 × 10 -3 , reducing the spectral signal of CH 4 and H 2 O. Haze, which could photochemically form from CH 4 , and refractory clouds that have not been considered here could further weaken the signals (e.g., Ma et al. 2023). A detailed investigation of the impact of photochemistry and photochemical haze/clouds on transmission spectra will be our future study.', '4.3. Comparison to the C-enrichment by atmospheric escape': "While we focused on the effect of a magma ocean on the atmospheric composition, another well-known factor that alters (or enriches) the atmospheric composition is atmospheric loss (e.g., Watson et al. 1981). When the planetary upper atmosphere receives the intense XUV, H 2 -rich atmospheres of sub-Neptunes undergo hydrodynamic escape, where individual atmospheric species escape to space at a rate that depends on the mass of the species. In this section, we estimate the change in atmospheric composition solely through such an atmospheric escape process and compare its effect against the effect of interaction with magma. Through this comparison, we aim to find the parameter space where the effect of magma can be identified.  \nWe mainly follow the methodology applied in Schaefer et al. (2016). The lightest atom, hydrogen, escapes from the atmosphere and drags heavier atoms (such as C and O) away. The escaping rate of heavier element ( ϕ 2 ) and the escaping rate of lighter element ( ϕ 1 ) satisfy the following equation (e.g., Hunten et al. 1987):  \nϕ 2 =    ϕ 1 X 2 X 1 µ c -µ 2 µ c -µ 1 , if µ 2 < µ c 0 if µ 2 > µ c , (51)  \nwhere the ( X 2 and X 1 ) are the molar concentration, and ( µ 2 and µ 1 ) are the atmic mass of each species. The µ c is the crossover mass, determining the mass criterion for drag to occur. The crossover mass follows the equation:  \nµ c = µ 1 + k B Tϕ 1 b 12 gX 1 m p , (52)  \nwhere k B is the Boltzmann constant, m p is the proton mass, and b 12 is the binary diffusion coefficient between two species. We refer to the neutral-neutral atom pair binary diffusion coefficient from Zahnle et al. (1990) and Yoshida et al. (2022), given by:  \nb ij = 1 . 96 × 10 6 T 0 . 5 µ 0 . 5 ij , (53)  \nwhere µ ij represents the reduced mass between species i and j. The escape flux without the heavier element (reference flux) is assumed to follow an energy-limited escape, satisfying the following equation (e.g., Chassefi'ere 1996; Schaefer et al. 2016):  \nϕ ref = ϵF XUV r p 4 Gµ 1 m p M p , (54)  \nwhere F XUV is the stellar XUV flux, r p is the planetary radius, M p is the planetary mass, and ϵ is the heating efficiency. The three escape rates ( ϕ ref , ϕ 1 , ϕ 2 ) are related to each other: ϕ ref = µ 1 ϕ 1 + µ 2 ϕ 2 .  \nThe escape rate and crossover mass equations suggest that the escape rate of heavier elements is smaller than that of lighter elements, and the difference in escape rates increases as the escape rate of the lighter element decreases. Therefore, we anticipate that atmospheric escape will lead to an increase in the atmospheric O/H and C/H ratios and a slight decrease in the atmospheric C/O ratio. The selective escape effect is expected to be more pronounced under weaker XUV radiation and smaller atmospheric amounts, provided that the XUV radiation remains sufficiently strong to cause atmospheric escape.  \nUsing the above equations, we calculate the enrichment of O/H and C/O under the following assumptions. First, we assume that the planet has evolved for 5 Gyr and finally has an atmosphere of the corresponding thickness and equilibrium temperature. Second, we consider an early M-type host star (M1-M2 type) with its XUV flux following the early M-type stellar evolution model proposed by Peacock et al. (2020). Third, we neglect the evolution of the planetary orbital radius and planetary radius. The orbital radius is translated to the planetary equilibrium temperature so as to match the parameterization of our magma-atmosphere model, with the assumption that the planetary bond albedo is 0.3 (e.g., Demory 2014). Finally, we maintain the heating efficiency ( ϵ ) at a constant value of 0.3, following the assumption from Schaefer et al. (2016).  \nThe solid lines in Figure 13 show the atmospheric O/H and C/H enrichment of 6M ⊕ sub-Neptunes due to the atmospheric escape as a function of atmospheric thickness and equilibrium temperature. The dashed lines represent the enrichment due to interaction with the magma ocean (without the effect of atmospheric loss), where the N before O(Fe) / N Fe ratio is set to 1.0. Data with ∆ R = 0 . 1 R ⊕ are excluded from the figure because magma cannot exist if N before O(Fe) / N Fe = 1 under such thin atmospheres due to low BOA temperature. The black region in Figure 13 indicates the region where the surface pressure exceeds 3 GPa, which is beyond the pressure range validated by experiments. We excluded this non-validated region from our calculations. The enrichment due to atmospheric escape is most pronounced when T eq ∼ 450 K, as a compromise between the amount of atmospheric loss (higher at higher equilibrium temperature) and the efficiency of the selective escape (higher at lower equilibrium temperature). As expected, thinner atmospheres exhibit a stronger atmospheric escape effect.  \nWe find that in the broad parameter space shown here, the effect of a magma ocean on O/H enrichment is dominant over that of atmospheric escape. Additionally, while magma does not effectively enrich atmospheric C/H, atmospheric escape enriches C/H to the same level as O/H because of the similar mass of C and O. Therefore, large O/H with small C/O (relative to the nebula values) would be safely attributed to the effect of the interaction between the atmosphere and magma, rather than the atmospheric escape. Escape can be accelerated by ionization (e.g., Dos Santos et al. 2023) compared to the neutral-neutral pair, in which case the enrichment of C/H and O/H due to escape is likely to be smaller. However, the effect on the C/O ratio is limited.", '4.4. Dependence on the atmospheric structure: Comparison to the radiative-convective profile': 'Our fiducial models assumed a simplified atmospheric T-P profile consisting of an isothermal stratosphere and adiabatic troposphere, which simplifies the analytic interpretation. The temperature profiles of real planets can  \nFigure 13. The enrichment in atmospheric O/H (left panel) and C/H ratio (right panel) relative to the initial nebula value due to species-dependent atmospheric escape (solid contour) and due to the interaction with magma (dashed contour). The x-axis denotes the planetary equilibrium temperature while the y-axis represents the atmospheric thickness. The black region indicates the parameter space where the BOA pressure is larger than the validated range ( ≲ 3 GPa) for the calculation of magma-bearing cases.  \n<!-- image -->  \ndeviate from this simple profile, and the relation between the observable parameters and the condition at the BOA can be different from our nominal models. To evaluate this effect, we additionally calculated temperature profiles in one-dimensional radiative-convective equilibrium using the open-source code HELIOS (Malik et al. 2017, 2019). In this subsection, we present the calculation set-up and the results, showing that the results with our nominal assumptions exhibit similar planetary parameter dependencies as discussed in Section 3.3 but with a lower O/H-C/H ratio and stronger C/O depletion.  \nThe calculation of radiative-convective equilibrium profiles with HELIOS requires (1) the assumption for the atmospheric composition, (2) the opacities of considered atmospheric species, and (3) the adiabtatic T-P profile for the convective regions. For (1), we consider H 2 , He, H 2 O, Na, and K. We ignore the presence of C-bearing species because they hardly affect the opacity of the atmosphere and the mean molecular weight. We should emphasize that the T-P profiles can be affected by the species not considered in our chemistry model. In particular, alkali metals such as Na and K have significant opacities at the optical wavelengths where H 2 , He, and H 2 O are largely transparent and thus substantially alter the T-P profile. The extensive survey of the range of the T-P profile will be the scope of future papers. In this paper, we simply set the ratio of He/H, Na/H, and K/H to be the values of the solar-composition atmosphere in chemical equilibrium at the given temperature, while the H 2 O mixing ratio is varied.  \nAs for (2), the k-coefficient tables up to 10 5 bar were calculated with HELIOS-K (Grimm & Heng 2015; Grimm et al. 2021). The opacity tables for the collision-induced absorption, H 2 -H 2 and H 2 -He, were created with HITRAN (Gordon et al. 2022), and that for H 2 O was based on the ExoMOL database (Polyansky et al. 2018). We extrapolate the opacity tables for Na and K from the tables provided in the default dataset up to 1000 bar, assuming that the opacity is proportional to the pressure.  \nRegarding (3), the adiabatic T-P profile, specifically d log T/d log P , was computed in a way consistent with our default model taking account of the temperature dependence of C p for H 2 , H 2 O, and He. The T-P profile deeper than 10 5 bar is simply extrapolated using the adiabatic lapse rate at the BOA.  \nOther major scalar input parameters include the internal temperature, which is assumed to be 50 K, similar to the value of Neptune (Pearl & Conrath 1991), the BOA pressure, which is set to 10 5 bar, and the stellar parameters where the radius is 0.45 R ⊙ and the effective temperature is 3500 K. The surface gravity in the radiative transfer routine is fixed to 20 m / s 2 , close to but smaller than the surface gravity at the BOA for our nominal case (23.4 m / s 2 ); the  \nFigure 14. The relationship among P BOA (red contour), T BOA (black contour), ∆ R , and T eq in the cases of radiative-convective equilibrium atmospheric profile. The radiative-convective equilibrium profile is further divided into the profile without Na and K (left) and including them (right). The x-axis denotes the planetary irradiation (equilibrium) temperature, while the y-axis represents the atmospheric thickness. The atmospheric composition is assumed to be H 2 (82%), He (17%), H 2 O(1%), Na (solar metallicity), and K (solar metallicity).  \n<!-- image -->  \nvariation of surface gravity in the range of 10-23.4 m / s 2 (corresponding to 0 < ∆ R < 0 . 8 R ⊕ ) has only minor effects in the T-P profile.  \nAfter computing the radiative-convective T-P profile, we use these profiles to find the surface pressure and temperature given ∆ R , now taking account of the altitude dependence of surface gravity. Figure 14 shows the relation among ∆ R , T eq , P BOA , and T BOA based on the radiative-convective calculation of HELIOS (see Figure 9 for comparison). Here, T eq is estimated from the irradation flux F irr by T eq ∼ ( F irr / (4 σ SB )) 1 / 4 , where σ SB is the Stefan-Boltzmann coefficient.  \nNotably, the HELIOS results indicate a much smaller T BOA than T BOA from the T-P profile of our model if effective opacity sources like Na and K are excluded. On the other hand, the inclusion of Na and K maintains relatively higher T BOA , which is slightly lower than or comparable to the T BOA in the case of our model assumption. Consequently, the range of planetary parameters (atmospheric thickness and equilibrium temperature) that ensures the presence of magma would be much narrower compared to the result with our nominal model setting if Na and K cannot contribute as effective opacity sources. However, it is essential to note that the effect of different values for T BOA on the trend and range of atmospheric composition is limited, as discussed in Section 3.1.  \nConversely, the P BOA of the HELIOS result is larger than that of our model. Consequently, atmospheric composition from the HELIOS profile is expected to have lower O/H and C/H ratios under the same N before O(Fe) / N Fe , as illustrated in Figure 5. In addition, this increased pressure effect will further lower the C/O ratio with a given O/H ratio, reinforcing our result that magma leads to a depleted atmospheric C/O ratio (see Figure 7).', '4.5. Limitations': 'There are mechanisms not covered in our model that are expected to suppress the C/O ratio further by either decreasing the amount of C or increasing the amount of O. The precipitation of C-bearing species into graphite or diamond can influence carbon partitioning, potentially reducing the atmospheric carbon content (e.g., Hirschmann 2012 and Maurice et al. 2024). Mechanisms that increase the amount of O include the precipitation of iron droplets. Contrary  \nto our assumption that iron droplets are suspended in magma, weaker magma convection would precipitate them, resulting in an oxidized magma. The partitioning of silicon between magma and the iron core (2Fe metal +SiO 2 ←→ 2FeO + Si metal ) may also introduce additional FeO into the mantle, leading to a more oxidized state than would otherwise be expected from the same initial rock composition. These mechanisms only reinforce our conclusion that magma can deplete the atmospheric C/O ratio.  \nIn addition to the oxidation of magma through silicon partitioning mentioned above, the role of the innermost iron core as a reservoir for hydrogen, carbon, and oxygen through efficient partitioning is widely recognized (e.g., Davies et al. 2018; Fischer et al. 2020; Li et al. 2020; Schlichting & Young 2022; Suer et al. 2023; Horn et al. 2023; Kim et al. 2023). Given that the iron core constitutes approximately 34% of the total planetary mass, volatile partitioning to the iron core can significantly impact atmospheric composition. However, these partitioning coefficients are sensitive to pressure, temperature, and silicate composition, emphasizing the importance of understanding the core formation process for accurate characterization. Therefore, a detailed discussion of the impact of the iron core on atmospheric composition remains in our future scope.  \nThere are admittedly uncertainties in our toy model of the solidification process. In particular, approximating the solidification as a two-step process likely overestimates the final atmospheric O, as the gradual oxidization of the magma-atmosphere system during solidification would leave more oxidized iron in the solidified rock. Also, the trapping of H 2 O in the solidified layer (e.g., Salvador et al. 2023) and the separation of magma into shallow surface and deep basal magma layers due to mid-crystallization (e.g., Caracas et al. 2019), which are ignored here, would suppress the oxidation of the magma-atmosphere system. While our conclusion that the solidification further oxidizes the atmosphere is likely to be valid, the quantitative assessment of the effect of the above-mentioned processes requires further modeling.  \nThe exclusion of rocky vapor in our model introduces uncertainties in the results, particularly under high temperatures. For example, rocky vapor can influence atmospheric composition through reactions such as the formation of silane: SiO + 3H 2 ←→ SiH 4 +H 2 O (e.g., Misener et al. 2023). Furthermore, the atmospheric structure can be affected by rocky vapor (e.g., Misener & Schlichting 2022). Investigating the impact of rocky vapor on the internal chemical reaction system and upper atmospheric composition will be our future scope.  \nWe anticipate that future geochemical studies will provide important implications, especially for the parameter ranges not covered in this study due to a lack of validated data. For instance, the dependence of solubility on melt composition has not been extensively investigated, except for a few species like H 2 , CO 2 , and H 2 O (e.g., Hirschmann et al. 2012; Allison et al. 2022). This is why we do not include a comprehensive analysis of melt composition (except for Fe content) in this study. If we were to infer the effect of melt composition on H 2 solubility (which has been discussed for three cases: Andesite, Basalt, and Peridotite; see Hirschmann et al. 2012), peridotite melt could reduce H 2 solubility by about a factor of three, further depleting atmospheric C/H and C/O ratios. The melt composition effects on the other two gases (H 2 O and CO 2 ) are likely limited since the dependence of H 2 O is less than factor of two for f H 2 O ≤ 2000 bar (e.g., Allison et al. 2022) and CO 2 is the most minor dissolving species in our model (see also Figure 4). The dependence of the solubility on the temperature is not well constrained from data either. Given that many sub-Neptunes likely have high-temperature and high-pressure surface conditions, the solubility data under extreme temperature-pressure environments will be invaluable for future modeling to properly interpret the sub-Neptune data.', '5. CONCLUSION': 'Motivated by the bulk composition degeneracy of sub-Neptunes, this study aimed to investigate the range of atmospheric H-C-O composition in magma-bearing sub-Neptunes formed through the nebula gas accretion on a dry, rocky core. We clarify the trend in the atmospheric composition and their dependence on magma properties (effective amount, Fe fraction, and initial Fe speciation), observable planetary parameters (planetary mass, equilibrium temperature, and radius), and an additional parameter for determining atmospheric structure (Radiative-Convective Boundary), using a chemical equilibrium model between the atmosphere and magma core. Furthermore, We discuss the effect of magma solidification on atmospheric composition. Our key findings are summarized as follows:  \n- 1. The atmospheric O/H ratio (or equivalently x H 2 O ) after the reaction between the accreted volatile and magma is a function of the redox state of the magma (which we parameterize by the number of iron-bound oxygen before the reaction, N before O(Fe) / N Fe ) and the surface pressure (P BOA ). The atmospheric H 2 O fraction is proportional to N before O(Fe) / N Fe at low pressure ( P BOA ≲ 10000 bar) and ( N before O(Fe) /N Fe ) 1 0 . 74 P -1 BOA at high pressure ( P BOA ≳ 10000 bar). The atmospheric  \n- H 2 O fraction ranges from a few percent to a few tens of percent (only for thin atmospheres, P BOA << 10000 bar) unless the initial rock is substantially reduced, with more than about 90% of Fe present in the form of pure iron.\n- 2. Despite the higher solubility of H 2 and H 2 O than C-bearing species, the C/H ratio experiences a modest increase, typically two to three times higher than the accreted (nebula) value. Only in highly oxidized thin-atmosphere scenarios does the C/H ratio exhibit enrichment of several tens of times.\n- 3. These two trends can be summarized in a simple relation between the atmospheric C/O ratio and atmospheric O/H ratio (or equivalently x H 2 O ) under the existence of magma. The atmospheric C/O ratio is proportional to x -1 H 2 O under a reduced atmosphere and x -0 . 26 H 2 O under an oxidized atmosphere. Thus, the atmospheric C/O ratio decreases to below 10% of the nebula value under an oxidized atmosphere with larger than a few percent of H 2 O fraction.\n- 4. A higher magma Fe fraction leads to a more oxidized atmosphere under constant N before O(Fe) / N Fe , as it increases the oxygen available for atmospheric-magma reactions. However, the Fe fraction does not alter the O/H-C/O relationship itself, maintaining a similar range of C/O ratios regardless of the Fe fraction.\n- 5. A smaller amount of the reactive magma causes the atmospheric composition to approach that of the accreted volatiles. A significant shift in the O/H ratio trend is observed when the reactive magma fraction falls below ∼ 10%, assuming M p = 6 M ⊕ , x H 2 O = 0 . 01, and P BOA = 10000 bar. Moreover, the quantity of reactive magma influences the O/H-C/O trend, where a smaller volume of the reactive magma leads to a smaller C/O ratio for a given O/H ratio.\n- 6. A thicker atmosphere with low equilibrium temperature, high radiative-convective boundary pressure, and heavier planetary core results in lower atmospheric O/H and C/H ratios due to the larger surface pressure, again not changing the O/H-C/O relation.\n- 7. Solidification of the magma further increases the atmospheric O/H ratio and decreases the atmospheric C/O\n- ratio, potentially by one or two orders of magnitude, through the extraction of H 2 O to the magma-atmosphere system.\n- Our results suggest that the atmospheric O/H - C/O ratio relation and the resulting depletion of the atmospheric C/O ratio to below 10% of the nebula value under an oxidized atmosphere can be evidence of magma-atmospheric interaction as well as no substantial additional supply of C. Further constraints on the volatile composition of icy solids from which sub-Neptunes could have formed would provide clearer insights into the distinction between icy sub-Neptunes from those formed through the nebula gas accretion onto the rocky core.\n- The authors deeply thank E. S. Kite for fruitful comments and H. Kuwahara and T. Tsuchiya for the discussions on 1\n- the material properties at high pressure. YI was supported by JSPS KAKENHI grant 22K14090. 2', 'A. THE DIFFERENCES IN RESULTS BASED ON THE EXTRAPOLATION METHOD FOR THE EQUILIBRIUM CONSTANT.': 'In this section, we present the results obtained with a different extrapolation method for equilibrium constants. As explained in section 2.2, we extrapolate the equilibrium coefficient equations for the IW reaction and Ferric-Ferrous equilibrium up to 4000 K even though the validation range for Ferric-Ferrous equilibrium ends at 2200 K. Figure 15 shows the results where we simply use the equilibrium constants at the boundary value of 2200 K if the temperature of the system exceeds 2200 K. The two results (Figure 5 and Figure 15) show that the boundary value extrapolation leads to a slightly more oxidized atmospheric composition, consistent with expectations from Figure 5. However, the difference is minute and does not affect the main results of this study.', 'REFERENCES': "```\nAl-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37, doi: 10.3847/1538-4357/ac0252 Allison, C. M., Roggensack, K., & Clarke, A. B. 2022, Contributions to Mineralogy and Petrology, 177, 40, doi: 10.1007/s00410-022-01903-y Ardia, P., Hirschmann, M. M., Withers, A. C., & Stanley, B. D. 2013, GeoCoA, 114, 52, doi: 10.1016/j.gca.2013.03.028 Armstrong, K., Frost, D. J., McCammon, C. A., Rubie, D. C., & Boffa Ballaran, T. 2019, Science, 365, 903,\n```  \n```\ndoi: 10.1126/science.aax8376\n```  \nFigure 15. The same figure as Figure 5, but generated using results obtained with a different extrapolation method. Note the negligible difference from Figure 5.  \n<!-- image -->  \nBean, J. L., Raymond, S. N., & Owen, J. E. 2021, Journal of Geophysical Research (Planets), 126, e06639,  \ndoi: 10.1029/2020JE006639 Benneke, B., Roy, P.-A., Coulombe, L.-P., et al. 2024, arXiv e-prints, arXiv:2403.03325, doi: 10.48550/arXiv.2403.03325 Burn, R., Mordasini, C., Mishra, L., et al. 2024, Nature Astronomy, 8, 463, doi: 10.1038/s41550-023-02183-7 Caracas, R., Hirose, K., Nomura, R., & Ballmer, M. D. 2019, Earth and Planetary Science Letters, 516, 202, doi: 10.1016/j.epsl.2019.03.031 Chachan, Y., & Stevenson, D. J. 2018, ApJ, 854, 21, doi: 10.3847/1538-4357/aaa459 Chase, M. 1998, NIST-JANAF Thermochemical Tables, 4th Edition (American Institute of Physics, -1) Chassefi'ere, E. 1996, J. Geophys. Res., 101, 26039, doi: 10.1029/96JE01951 Davies, C. J., Pozzo, M., Gubbins, D., & Alf'e, D. 2018, Geophys. Res. Lett., 45, 6042, doi: 10.1029/2018GL077758 Demory, B.-O. 2014, ApJL, 789, L20, doi: 10.1088/2041-8205/789/1/L20 Deng, J., Du, Z., Karki, B. B., Ghosh, D. B., & Lee, K.  \nK. M. 2020, Nature Communications, 11, 2007,  \ndoi: 10.1038/s41467-020-15757-0  \nDos Santos, L. A., Garc'ıa Mu˜noz, A., Sing, D. K., et al. 2023, AJ, 166, 89, doi: 10.3847/1538-3881/ace445  \nFischer, R. A., Cottrell, E., Hauri, E., Lee, K. K. M., & Le Voyer, M. 2020, Proceedings of the National Academy of Science, 117, 8743, doi: 10.1073/pnas.1919930117 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661, doi: 10.1086/512120 Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81, doi: 10.1088/0004-637X/766/2/81 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb Gaillard, F., Bouhifd, M. A., Furi, E., et al. 2021, SSRv, 217, 22, doi: 10.1007/s11214-021-00802-1 Gaillard, F., Bernadou, F., Roskosz, M., et al. 2022, Earth and Planetary Science Letters, 577, 117255, doi: 10.1016/j.epsl.2021.117255 Gale, A., Dalton, C. A., Langmuir, C. H., Su, Y., & Schilling, J.-G. 2013, Geochemistry, Geophysics, Geosystems, 14, 489, doi: 10.1029/2012GC004334 Gao, P., Piette, A. A. A., Steinrueck, M. E., et al. 2023, ApJ, 951, 96, doi: 10.3847/1538-4357/acd16f Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, JQSRT, 277, 107949, doi: 10.1016/j.jqsrt.2021.107949"}