© The Authors 2025

1 Introduction

The NASA Kepler mission dedicated to the search for exoplanets (Borucki et al. 2010) revealed that objects with radii smaller than Neptune are the most abundant in the galaxy, making up more than three-fourths of the entire population (Borucki et al. 2011; Batalha et al. 2013). The distribution of exoplanet radii measured with Kepler exhibits a gap at ≈ 1.5 R_⊕ (Fulton et al. 2017), referred to in the literature as the “radius gap,” that marks the transition between the categories of super-Earth (R_p < 1.5 R_⊕) and sub-Neptune (R_p > 1.5 R_⊕). The latter objects are likely composed of a rocky core surrounded by a H₂-He primitive envelope that survived erosion during the early active phase of the star, whereas the former lost this light primitive atmosphere through extreme ultraviolet (XUV) stellar heating (i.e., photo-evaporation) and/or the transfer of energy from a cooling interior (i.e., core-powered mass loss) (Owen & Wu 2013; Lopez & Fortney 2013; Fulton & Petigura 2018; Loyd et al. 2020; Van Eylen et al. 2018; MacDonald 2019; Gupta & Schlichting 2019).

The dominant mass-loss mechanism shaping the radius valley remains debated in the literature; both photo-evaporation and core-powered mass-loss can play a crucial role at different stages in the planet’s evolution (Owen & Schlichting 2024). With this primitive envelope lost, planets with sizes below the radius gap could host secondary atmospheres produced via outgassing, similarly to the rocky planets in the inner Solar System (Gaillard & Scaillet 2014; Gillmann et al. 2024).

For rocky exoplanets below the radius gap, a broad range of bulk densities have been reported, suggesting variations from water worlds to Earth-like silicate-rich objects to Fe-rich super-Mercuries (Weiss & Marcy 2014; Luque & Pallé 2022). The variety of metallicities observed in the host stars (Buchhave et al. 2012; Brewer & Fischer 2016) further supports the expected diversity in interior composition. Planetary mass and radius (M-R) measurements have provided crucial first insights into the internal structure, suggesting that a large fraction of rocky exoplanets have a massive core (Dorn et al. 2015; Liu & Ni 2023; Unterborn et al. 2023). Although the core to mantle fraction can be well constrained with M-R measurements, the degeneracy on the mantle composition is substantial (Seager et al. 2007) and can only be solved if the relative abundances of refractories from the star are known (Dorn et al. 2015; Unterborn et al. 2016). This approach assumes a comparable composition between star and planet, a reasonable and validated starting hypothesis (Thiabaud et al. 2015), with some exceptions where planet formation mechanisms could modify the correlation (Adibekyan et al. 2021; Schulze et al. 2021).

The rocky interior properties, including volatile budget and redox state, control the composition of the surrounding secondary atmosphere (Gaillard et al. 2021). The speciation of gas species is primarily controlled by the interior redox state quantified with the oxygen fugacity (fo₂), a key parameter reflecting the relative abundance of ferric iron relative to the total iron content (Fe³⁺/ΣFe) in the melt, with high fo₂ favoring the oxidation of ferrous iron in the buffer reaction FeO + $\frac{1}{4}$ O₂ ⟷ FeO_3/2 or the oxidation of metallic iron in the buffer reaction Fe + $\frac{1}{2}$ O₂ ⟷ FeO (Frost & McCammon 2008; Burgisser et al. 2015; Sossi et al. 2020; Deng et al. 2020; Ortenzi et al. 2020; Liggins et al. 2022; Guimond et al. 2023). This latter buffer reaction, called the iron-wustite (IW) buffer, is often used to express the fo₂ relative to a reference (e.g., Tian & Heng 2024). The melt oxygen fugacity is controlled by pressure (Deng et al. 2020) as well as (volatile) composition (Young et al. 2023).

For ultra-short-period (USP) exoplanets, the high surface temperature forced by incoming stellar radiation should favor the formation of a sustained magma ocean in these so-called “lava worlds.” In this scenario, the atmosphere is directly at equilibrium with the upper layers of the magma ocean, and the characterization of the gas composition can provide crucial information on the interior properties (Gaillard et al. 2021, 2022). Several models have been developed to describe the composition of the atmosphere at equilibrium with the magma ocean accounting for the solubility of volatile species in silicate melts (Sossi et al. 2020; Gaillard et al. 2022; Liggins et al. 2022; Bower et al. 2022; Wolf et al. 2023). A large diversity of atmospheric compositions can be produced depending on the fo₂, from CO- to H₂- to CO₂-dominated atmospheres (Gaillard et al. 2022; Tian & Heng 2024). However, steam atmospheres are less likely given the high solubility of water in silicate melts (Bower et al. 2022; Sossi et al. 2023). Based on the strong correlation between the fo₂ of the magma and the atmospheric composition, it has been suggested that specific relative abundances of volatiles in the atmosphere can help constrain the oxygen fugacity of the interior, with CO₂/CO being the most robust candidate (Tian & Heng 2024; Shorttle et al. 2024). Beyond volatiles, the atmospheres of lava worlds can contain significant amounts of vaporized refractories including silicate and magnesium oxides, which can also help constrain the interior redox state (Wolf et al. 2023; Seidler et al. 2024).

For secondary atmospheres of cooler planets with a crystallized crust, the mass and composition of the atmosphere result from a balance between volcanic degassing and atmospheric escape (Ortenzi et al. 2020; Liggins et al. 2020). The degassing rate is mainly controlled by planet mass and internal heating, the latter being controlled via tidal or radiogenic processes that allow partial melting in the upper mantle (Dorn et al. 2018; Quick et al. 2020; Oosterloo et al. 2021). For M_p < 2−3 M_⊕, the degassing rate increases with planet mass, making secondary atmospheres around massive super-Earths more likely (Dorn et al. 2018) unless the convection efficiency is reduced by a higher core fraction (Oosterloo et al. 2021). In simplified models, the degassing rate is derived using the mantle melt fraction, rock-melt partition coefficients, and the approximated mantle convection rate (Ortenzi et al. 2020; Liggins et al. 2022). The process of degassing is thus subjected to too many unknown variables, although the outgassed composition is mainly controlled by the volatile budget and the oxygen fugacity (fo₂) in a manner similar to a magma ocean scenario (Gaillard & Scaillet 2014; Liggins et al. 2020, 2022; Tian & Heng 2024). Although the planetary size is a key parameter controlling the interior fo₂, variations in the mantle mineralogy can also modify fo₂ by orders of magnitude (Guimond et al. 2023). In other words, the large diversity in planetary mass-radius and composition (refractories) observed in the exoplanet population will likely produce rocky exoplanet atmospheres with a variety of compositions. Contrary to a magma ocean scenario, a difference in temperature is expected between the melt and the surface of the planet. This process where the warm degassed volatiles cool in a new atmospheric environment, which we refer to as atmospheric cooling, can modify the composition of the atmosphere and thus affect observations (Liggins et al. 2023). Other processes affecting the composition of the upper atmosphere probed during spectroscopic observations, such as escape and photochemistry, can modify the degassed composition (Liggins et al. 2020; Krissansen-Totton & Fortney 2022) and thus bias our objective, which is focused on inferring properties of the rocky interior.

The atmospheric composition and relative abundances of volatile species with different redox states is first set by geochemical outgassing but modified by atmospheric chemistry (photochemistry and thermochemistry). In the present study, we focus on secondary atmospheres of evolved planetary bodies with a crystallized crust and assess the biases caused by photochemistry and atmospheric cooling as well as their influence on future atmospheric characterization in order to constrain the interior fo₂. In Sect. 2, we present the models we used and the simulations performed for this study considering geochemical outgassing with self-consistent calculations of atmospheric chemistry and radiative transfer. In Sect. 3, we describe the effect of atmospheric cooling (thermochemistry) on relative molecular abundances used as tracers of the interior fo₂. In Sect. 4, we assess the effect of atmospheric cooling, photochemistry, and vertical mixing on the relative abundances for a 1D case study of the super-Earth GJ 1132b. In Sect. 5, we discuss the implications for future observations with the James Webb Space Telescope (JWST).

2 Model description

2.1 Outgassing model: C–H–O gas at equilibrium with carbon-bearing melt

The outgassing model used for this study is described in detail in Tian & Heng (2024) following an approach initially introduced by French (1966). The model assumes chemical equilibrium between a C-H-O fluid phase and a carbon-bearing melt. At the pressures considered in this study (below 100 bars), ideal gas may be assumed (Tian & Heng 2024) which simplifies the determination of the gas partial pressures to an analytical calculation. Only four parameters are needed to determine the speciation of the gas at equilibrium with the carbon-bearing melt: oxygen fugacity, activity of graphite (A_c), melt temperature (T_m) and degassing pressure (P_d). The chemical system can be described by a set of four independent reactions: $$\mathrm{C}+{\mathrm{O}}_2⟷{\mathrm{C}\mathrm{O}}_2$$(R1) $$\mathrm{C}+\frac{1}{2}{\mathrm{O}}_2⟷\mathrm{C}\mathrm{O}$$(R2) $$\mathrm{C}+2{\mathrm{H}}_2⟷{\mathrm{C}\mathrm{H}}_4$$(R3) $${\mathrm{H}}_2+\frac{1}{2}{\mathrm{O}}_2⟷{\mathrm{H}}_2\mathrm{O}$$(R4)

The IW buffer is expressed following Ballhaus et al. (1991) based on the experimental data of O’Neill (1987). The current model does not consider solubility, and it thus does not require volatile inventories as an input parameter unlike other models applying mass balance (e.g., Gaillard & Scaillet 2014 or Liggins et al. 2020). In practice, it means that the total gas pressure is set in our model and not calculated considering the solubility of each gaseous species in the melt. We emphasize that the lack of solubility in the model does not bias the correlation between CO₂/CO and fo₂ reflected with reactions (R1) and (R2). CO₂/CO remains directly correlated to fo₂ following X_CO2/X_CO = fo₂^0.5/K_eq, with X the molecular abundance and K_eq the equilibrium constant of the reaction CO + 0.5 O₂ ⟷ CO₂ (i.e., (R1)–(R2)). Although the model of Tian & Heng (2024) can consider N- and S-bearing species, these will not be considered in the present work since our approach targets molecular abundances that can directly inform on the interior redox state such as CO₂/CO. Speciation between N₂ and NH₃ is controlled by the amount of H₂ in the system which will be difficult to constrain with atmospheric characterization. In addition, the solubility of ammonia in the melt, which is not considered here, can be important to explain the abundances found in the atmosphere (Shorttle et al. 2024). Molecular abundances of sulfur-bearing compounds are also highly sensitive to solubility, even at low pressures.

In our calculations, T_m is set to 1600 K, a typical melt temperature used in previous studies (Gaillard & Scaillet 2014; Gaillard et al. 2022). Degassing is assumed to occur under atmospheric pressure (i.e., P_d = P_s) similar to the approach of Gaillard & Scaillet (2014).

2.2 OD calculations of atmospheric cooling

Rocky exoplanets currently observed with JWST exhibit different masses and radii (Luque & Pallé 2022) with likely different interior compositions (suggested by stellar observations) that would create different atmospheric compositions. To avoid the large number of unknowns required to compute outgassing and escape fluxes, we consider a state where both processes balance using a parametrized surface pressure. In other words, the outgassing model is used to set the initial composition of the atmosphere at a given pressure, we do not consider a continuous volcanic flux at the surface. Boundary conditions including surface deposition and escape are not considered as we focus on the effect of atmospheric chemistry. Most rocky exoplanets observed with JWST present equilibrium temperatures between 400 K and 1500 K. If these objects maintained an atmosphere, the volatiles outgassed from the interior will cool to the new temperature found in the atmosphere. This atmospheric cooling mechanism will modify the stable gas abundances according to the new pressure-temperature conditions found in the atmosphere (Liggins et al. 2023).

After atmospheric cooling, the molecular abundances might change significantly compared to the prediction of the outgassing model. To assess the change in molecular abundances of a volcanic mixture subjected to atmospheric cooling, we combine the geochemical outgassing calculations of Tian & Heng (2024) with atmospheric thermochemical calculations performed using FastChem and VULCAN. FastChem, described in detail in Stock et al. (2018) and Stock et al. (2022), uses a semi-analytical approach to determine the molecular abundances at chemical equilibrium for a given pressure, temperature and set of elemental abundances. VULCAN, described in Tsai et al. (2017) and Tsai et al. (2021), is a 1D chemical kinetics code used to study equilibrium chemistry and/or photochemistry. VULCAN is used to understand the change in molecular abundances using a detailed analysis of the chemical network and timescales. We use FastChem in parallel to compare with VULCAN and ensure that the kinetic simulations successfully reached chemical equilibrium. This step is crucial, as we consider a large range of pressures and temperatures that will result in various timescales of thermochemistry. The C−H−O chemical network was validated in Tsai et al. (2017) and Tsai et al. (2018). We updated the chemical kinetic constants taken from the NIST chemical kinetics database¹. A few reactions were added following the work of Liggins et al. (2023). Most kinetic constants used are compiled in the reviews of Tsang & Hampson (1986), Tsang (1987) and Baulch et al. (1992) as they were validated in a broad range of gas temperature ranging from 300 to 2500 K. For a few reactions, we were compelled to use kinetic constants derived from theoretical calculations when experimental data were not available in a broad range of temperature. The database counts 428 chemical reactions (forward and reverse) involving only C−H−O species.

The aim of this first 0D analysis is to assess the reliability of CO₂/CO as an atmospheric tracer for the interior fo₂ despite the expected change in temperature between the melt and the surface of the evolved crystallized rocky planet. We explore the parameter space of the model now including the surface temperature T_s in addition to the parameters of the outgassing model. We perform random sampling using log-uniform and uniform laws to cover the entire parameter space within the range presented in the Table 1.

2.3 1D calculations combining radiative transfer and atmospheric chemistry

The 0D calculations provide crucial insights into the main chemical mechanisms resulting from atmospheric cooling. In a 1D atmosphere, however, the atmospheric temperature and composition are strongly correlated. The abundances of atmospheric species are affected by the atmospheric temperature changing throughout the column atmosphere. The composition and more specifically the greenhouse properties of the atmospheric mixture also controls the temperature profile. This complexity compels us to combine different calculations using a climate-chemistry model. For that purpose, we coupled the radiative-convective model HELIOS (Malik et al. 2017, 2019a) with the chemistry model VULCAN to perform self-consistent calculations including both radiative transfer and chemical kinetics (photochemistry and thermochemistry). The geochemical calculations of Tian & Heng (2024) are used to set the initial atmospheric composition. It can also be used as a lower boundary condition in VULCAN using a parameterized volcanic flux. As stated previously, we refrain from using this continuous outgassing feature in this work given the lack of constraints regarding the flux. Given the broad range of temperature encountered in our case study, competition between mixing, radiative and chemical timescales is expected which motivates the need to use a self-consistent model. In our atmospheric model, the temperature profile is changed by HELIOS during the VULCAN simulation using the radiative timescale as an iteration time step. The radiative timescale is updated throughout the simulation to modify the time step if needed, it is expressed as follows (Malik et al. 2017): $$t_{rad}=\frac{C_{p}P}{\sigma g{T}^3},$$(1)

where P is the pressure, C_p is the heat capacity, σ is the Stefan-Boltzmann constant, g is the gravity and T the atmospheric temperature.

The mean heat capacity of the atmosphere is derived using the abundances of the main atmospheric molecules (CO, CO₂, CH₄, H₂ and H₂O) and their C_p at the given surface temperature. C_p(T) of each species is calculated using the thermodynamic data (polynomial constants) from Holland & Powell (1998) also used in the geochemical outgassing calculations (Tian & Heng 2024). The change in temperature is read by VULCAN to update the kinetic constants and reaction rates. Temperature-dependent photo-absorption cross sections and molecular diffusion rates are also changed according to the new temperature conditions.

We include a transit spectra simulator to read the outputs of the atmospheric model. We use the forward model of HELIOS-T first described in Fisher & Heng (2018) which we modified to apply generally for any atmospheric composition and consider temperatures and mixing ratios varying with pressure. HELIOS provides the secondary eclipse emission spectrum of the planet along with the contribution function of different spectral bins to assess the pressure range probed during spectroscopic observations.

Pressure is expected as a critical parameter influencing the dominant chemical regime shaping the atmospheric composition between thermochemistry and photochemistry. Since pressure itself affects both the chemistry and the radiative balance, it is essential to consider both atmospheric chemistry and radiative transfer in the model. The structure of the model is illustrated in Fig. 1. We explore the parameter space of the atmospheric model guided by a specific question, namely, how relative molecular abundances can help constrain the interior fo₂ from observations. We want to clarify the role of photochemistry and thermochemistry. If deterministic mechanisms are at play, our aim is to identify how we can correct for them in future applications with a retrieval framework.

We focus on a case study of GJ 1132b as recent observations suggest the presence of an atmosphere (Swain et al. 2021; May et al. 2023). In addition, the equilibrium temperature of this exoplanet is close to that of several other super-Earths which ensures that the results and underlying mechanisms identified can be generalized to other objects. We focus on three different end-members of oxygen fugacity and atmospheric surface pressure to evaluate the influence of both parameters. The different end-member cases are summarized in Table 2.

In the outgassing calculation, melt temperature is kept fixed at 1600 K following Gaillard & Scaillet (2014). The process of atmospheric cooling is mainly sensitive to the difference between melt temperature and atmospheric temperature, hence our choice to only vary atmospheric temperature. Activity of graphite is set to ensure an atmosphere rich in carbon (see Table 2). If we consider degassing from a carbon-poor melt, only H₂ and H₂O would present high abundances that could be detected limiting our ability to infer a quantitative redox state. In the atmospheric model, photochemical haze formation is not included, the surface albedo does not impact our simulations as surface emissivity is set accordingly (albedo = 1 – emissivity). Stellar temperature, radius and mass are set for GJ 1132 (Berta-Thompson et al. 2015) and used to produce the blackbody spectrum in HELIOS. The planet’s surface gravity is also taken from Berta-Thompson et al. (2015). Surface pressure is set according to our three different scenarios (Table 2). The upper boundary pressure (top of atmosphere) is set to 10⁻⁷ bar to ensure the accuracy of the UV optical depth and vertical distribution of actinic flux in the 1D atmosphere for the photochemical calculations. The stellar flux of GJ 1132 is subjected to negative values. We thus used the stellar flux of GJ 1214, which presents a similar temperature and radius, with data available as part of the MUSCLES survey (France et al. 2016). The climate model and transit spectra generator includes the opacities of CO (Rothman et al. 2010; Li et al. 2015), CO₂ (Rothman et al. 2010), CH₄ (Rothman et al. 2010; Hargreaves et al. 2020), H₂O (Polyansky et al. 2018) along with collision-induced absorption (CIA) by H₂-H₂ (Abel et al. 2011). Other CIAs are not included given the lack of data in a broad spectral and temperature range (Chubb et al. 2024). Rayleigh scattering is considered in the atmospheric model (VULCAN, HELIOS, HELIOS-T) for the different species following the formulation in MacDonald & Lewis (2022) using the King factor, refractive index and reference number density with most data taken from Sneep & Ubachs (2005). The water refractive index is pressure-temperature dependent, the data is taken from Murphy (1977) and Schiebener et al. (1990). Vertical mixing described by the eddy diffusion coefficient K_zz is one of the main sources of uncertainty in the chemical model, especially since the temperature is not known a priori. The role of K_zz is assessed and discussed in Sects. 4.1 and 4.2. Molecular diffusion is considered in VULCAN (Tsai et al. 2021) with diffusion coefficients taken from Marrero & Mason (1972) following the formalism in Hu et al. (2012). Photo-absorption and -dissociation cross sections were compiled from the Leiden Observatory database² in the original study of Tsai et al. (2021). Temperature-dependent cross sections are available for CO₂ (150–800 K; Venot et al. 2018) and water (ExoMol³; 300–573 K), although not in the entire range of temperature covered with our simulations (≈400–1600 K). More temperature-dependent cross-section data will be needed in the future to address their influence properly (Chubb et al. 2024). The database of chemical reactions is similar to the one used for the 0D simulations (Sect. 2.2) with the addition of 40 photolysis reactions to consider photochemistry.

3 Influence of atmospheric cooling on CO₂/CO and CH₄ production

3.1 Relative molecular abundances: Inferring fo₂ and temperature using CO, CO₂, and CH₄

The relative molecular abundances obtained with the 0D atmospheric cooling simulations are shown in Fig. 2. We remind the reader that these simulations assume thermochemical equilibrium at the planet’s surface temperature set by T_s, lower than the melt temperature. Even after atmospheric cooling, CO₂/CO and fo₂ remain strongly correlated with a Pearson coefficient at 0.88. The lower correlation compared to the outgassing prediction is explained by a decrease of CO₂/CO as T_s increases for a given fo₂ (top right panel, Fig. 2). The correlation between CO₂/CO and fo₂ decreases as the atmospheric temperature departs from the melt temperature, it is thus caused by atmospheric cooling. For T_s between 500 and 1600 K, CO₂/CO increases by up to 3–4 orders of magnitude for a given fo₂. A constraint on the temperature is therefore essential to correct biases induced by atmospheric cooling.

Figure 2 shows that a single CO₂/CO can correspond to a broad range of atmospheric temperature. In other words, this molecular pair cannot be used in observations to infer atmospheric temperature. CO–CH₄ conversion is, however, very sensitive to atmospheric temperature (top right panel, Fig. 2). The methane abundance can thus be used to constrain the atmospheric temperature and improve the retrieved fo₂. Using the abundances of CO, CO₂ and CH₄, one can in theory retrieve accurately the oxygen fugacity and atmospheric temperature (top panels, Fig. 2). The feasibility of this latter statement lies in the absolute abundances of each species and their detectability in the spectral range of JWST. In practice, high abundances of methane are found only under reducing conditions. Inferring temperature would be more difficult for oxidized atmospheres. In theory, H₂/H₂O provides a similar information on fo₂ as CO₂/CO (bottom left panel, Fig. 2). Although the H₂ abundance cannot be retrieved from observations, it suggests that low oxygen fugacity objects would create thicker atmospheres which would affect transit observations through the scale height as discussed in Ortenzi et al. (2020).

3.2 Chemical network and timescales

In the geochemical outgassing model, no information on the kinetics is provided as equilibrium between the gas and the melt is assumed (Tian & Heng 2024; French 1966). By calculating the Gibbs free energy and equilibrium constant of each independent reaction (R1)–(R4), one can derive the partial pressures of the main stable species involved in the network (French 1966; Tian & Heng 2024). For the atmospheric calculations, one could use a similar approach where the activity of graphite and oxygen fugacity are replaced by the elemental abundances of the gas mixture (Heng & Tsai 2016; Heng et al. 2016). With chemical kinetics, one understands that the partial pressures of the stable molecules and the timescales of chemical equilibrium are controlled by reactions with the radical products of thermal dissociation, for example, H₂O ⟶ OH + H. The abundances of these radicals is changed with temperature. Since we previously showed that atmospheric cooling can strongly modify CO₂/CO, we can look at the radicals and reactions to understand why. Figure 3 shows the number density ratios OH/H and CO₂/CO as a function of oxygen fugacity and temperature for the different Monte Carlo simulations performed, covering the parameter space described in Table 1. The results teach us that OH/H is correlated to both fo₂ and atmospheric temperature. For a given temperature, the fo₂ corresponds to a specific O/H or OH/H. We now need to look at the chemical network to understand why the change in radical density also affects molecular abundances such as CO₂/CO.

For each simulation, we have a specific atmospheric/degassing pressure, atmospheric temperature and elemental abundances predicted by the outgassing calculations. Using kinetic simulations performed with VULCAN, we evaluate the different reaction rates and find the fastest chemical pathway between two molecules (e.g., CO and CH₄) using a Dijsktra algorithm (Dijkstra 1959; Tsai et al. 2018). This analysis gives us the main chemical pathway controlling a specific conversion at a given pressure-temperature and for a given set of elemental abundances. Within this fastest pathway, we identify the rate limiting step (RLS) or slowest reaction which bottlenecks the chemical pathway and controls the timescale. In most scenarios, the chemical timescale can approximated as $$t_{eq,X}=\frac{n_X}{V_{RLS,X}},$$(2)

where n is the number density (molecules cm⁻³) of species X at chemical equilibrium and V_RLS is the slowest reaction rate in the fastest pathway (molecules cm⁻³ s⁻¹).

The outgassing model predicts a broad range of elemental abundances, although C/O is usually below 1 given the difficulty in producing CH₄-rich atmospheres. O/H varies significantly depending on fo₂ giving very high metallicities poorly studied as previous work focused on hot Jupiters or hot Neptunes and the effect of C/O (Moses et al. 2011, 2013a,b; Venot et al. 2015; Tsai et al. 2018). Our network analysis shows that CO–CO₂ and H₂− H₂O conversions are always controlled by the following scheme in the entire parameter space (pressure, temperature, elemental abundances): $$\mathrm{C}\mathrm{O}+\mathrm{O}\mathrm{H}⟷{\mathrm{C}\mathrm{O}}_2+\mathrm{H}$$(R5) $${\mathrm{H}}_2\mathrm{O}+\mathrm{H}⟷{\mathrm{H}}_2+\mathrm{O}\mathrm{H}$$(R6)

net: CO+H₂O ⟷ CO₂ + H₂

From this pathway, one can understand why variations of CO₂/CO and OH/H are directly correlated. First, the increase in CO₂/CO with fo₂ at a given temperature is caused by the increase in OH/H favoring the oxidation of CO. An increase in fo₂ corresponds to an increase in O/H and therefore OH/H. However, OH/H also varies with temperature for a given O/H, as thermal dissociation rates of water and molecular hydrogen vary differently with temperature (Fig. 3). H₂/H₂O is controlled by OH/H and therefore decreases with increasing fo₂ and increasing temperature. This analysis explains the trends of relative abundances as a function of fo₂ and temperature in Fig. 2 (bottom panels). CO₂/CO on the other hand increases with increasing fo₂ but decreases with increasing temperature. This result can seem counter-intuitive but is explained by the balance in elemental abundances. The decrease in OH/H and H₂O/H₂ with increasing temperature at a given fo₂ leads to an increase of CO₂/CO to conserve elemental abundances and especially O/H. The trends of relative molecular abundances for CO₂, CO, H₂ and H₂O are completely explained following the change of the radical density ratio OH/H with temperature and fo₂.

We also evaluated CO–CH₄ conversion for the different simulations. We identified the following scheme: $$\mathrm{C}\mathrm{O}+\mathrm{H}+\mathrm{M}⟷\mathrm{H}\mathrm{C}\mathrm{O}+\mathrm{M}$$(R7) $$\mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{H}}_2⟷{\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{H}$$(R8) $${\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{H}+\mathrm{M}⟷{\mathrm{C}\mathrm{H}}_3\mathrm{O}+\mathrm{M}$$(R9a) $${\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{H}+\mathrm{M}⟷{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}+\mathrm{M}$$(R9b) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}+{\mathrm{H}}_2⟷{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{H}$$(R10a) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}+{\mathrm{H}}_2\mathrm{O}⟷{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{O}\mathrm{H}$$(R10b) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{M}⟷{\mathrm{C}\mathrm{H}}_3+\mathrm{O}\mathrm{H}+\mathrm{M}$$(R11a) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{H}⟷{\mathrm{C}\mathrm{H}}_3+{\mathrm{H}}_2\mathrm{O}$$(R11b) $${\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}+\mathrm{H}⟷{\mathrm{C}\mathrm{H}}_3+\mathrm{O}\mathrm{H}$$(R12) $${\mathrm{C}\mathrm{H}}_3+{\mathrm{H}}_2⟷{\mathrm{C}\mathrm{H}}_4+\mathrm{H}$$(R13a) $${\mathrm{C}\mathrm{H}}_3+{\mathrm{H}}_2\mathrm{O}⟷{\mathrm{C}\mathrm{H}}_4+\mathrm{O}\mathrm{H}$$(R13b) $$\mathrm{O}\mathrm{H}+{\mathrm{H}}_2⟷{\mathrm{H}}_2\mathrm{O}+\mathrm{H}$$(R14) $$\mathrm{H}+\mathrm{H}+\mathrm{M}⟷{\mathrm{H}}_2+\mathrm{M}$$(R15)

net: CO + 3 H₂ ⟷ CH₄ + H₂O

This general scheme is valid for the different simulations performed, although some steps can be associated with different reactions depending on the temperature, pressure and elemental abundances. The oxygenated species produced during CO–CH₄ conversion are, however, similar. As our outgassing calculations predict C/O≤1, the CO–CH₄ scheme involves primarily oxygenated species such as methanol and formaldehyde. For C/O≥1, the fastest pathway would also involve saturated and unsaturated reduced hydrocarbons (Moses et al. 2011, 2013a; Venot et al. 2015; Tsai et al. 2018).

The reaction converting CH₃O−CH₃ OH involves either H₂ and H (R10a) or H₂O and OH (R10b) depending on the oxygen content in the atmosphere and more specifically depending on O/H. The reaction (R10a) dominates for O/H < 0.5, whereas (R10b) dominates for higher O/H. A similar result is found for the CH₃-CH₄ conversion. For high O/H, oxidation of methane dominates through reaction (R13b), whereas hydrogenation dominates for lower O/H (R13a). The CH₃-CH₄ reaction might be changed with H₂S when sulfur is included (Liggins et al. 2023), although this depends on the sulfur content in the mixture.

The reaction converting CH₃-CH₃OH or CH₃-CH₂OH depends on the temperature-pressure conditions. Using the analysis of the reaction rates in the VULCAN simulations, we find that (R11a) and (R12) are competing at temperatures above ≈900 K. Reaction (R11a) dominates at high pressure as methanol formation/dissociation is pressure-dependent. Reaction (R11a) is indeed a three-body reaction which becomes inefficient at low pressures. The two-body reaction (R12) becomes more efficient than (R11a) at low pressures. The fastest pathway is thus going directly from CH₂OH to CH₃ without methanol at low pressure. The pressure value marking the transition from (R11a) to (R12) depends on the temperature. This pressure value increases with temperature from about 10⁻² bar at 900 K to 10 bar at 1600 K in agreement with the analysis of Tsai et al. (2018). At temperatures below 900 K, (R11a) and (R11b) are competing. For a high OH/H, (R11a) would dominate as it involves OH compared to (R11b) involving H. The transition between (R11a) and (R11b) for the limiting reaction depends on both temperature and O/H since both parameters affect OH/H.

For every simulation, either (R11a), (R11b) or (R12) is identified as the limiting reaction (RLS) for the thermochemical production of methane. The CO–CH₄ conversion requires the dissociation of three bonds binding C and O together. The timescale of this conversion is limited by the dissociation of the last bond in methanol or CH₂OH. These results are in agreement with previous analyses of CO–CH₄ conversion for exoplanet atmospheres (Moses et al. 2011; Tsai et al. 2018; Liggins et al. 2023). We refer the reader to previous work by Tsai et al. (2018) for a more detailed analysis of RLS changes in various conditions of pressure, temperature and elemental abundances.

Our results show that the change in OH/H controlled by oxygen fugacity and temperature has a significant impact on the CO–CH₄ and CO–CO₂ conversions. Differences are thus expected in terms of timescales. Figure 4 presents the thermochemical timescales of CO derived for each simulation as a function of temperature, pressure and fo₂. In general, a broad range of thermochemical timescale is possible for a single temperature depending on the fo₂ and pressure (Fig. 4).

It was previously shown that the timescale of CO can be derived by identifying the RLS of the CO–CH₄ conversion for low metallicity atmospheres representative of hot Jupiters (Tsai et al. 2018). For high O/H, that statement is no longer valid, as the oxidation of CO with the radical OH becomes more efficient than its hydrogenation that initiates the CO–CH₄ conversion. The timescale of CO can be approximated using (R5) as RLS in Eq. (2) when the rate of (R5) is higher than the rate of (R7). We note that the timescale of CO at low temperatures and high fo₂ joins the timescales trend obtained for low fo₂ (Fig. 4, bottom panel) and is thus controlled by CO–CH₄ conversion. This stems from the decrease in OH/H at a lower temperature and the resulting less efficient CO oxidation relative to hydrogenation. For low O/H, the timescale of CO is controlled by CO–CH₄ conversion and thus by the RLS in the given pressure-temperature conditions.

4 Case study using a climate-chemistry model: How pressure, temperature, mixing, and UV photons impact the atmospheric tracers of fo₂

The trends in relative molecular abundances and timescales observed inFigs. 2–4 with atmospheric cooling (thermochemistry alone) are caused by the changes in OH/H with both temperature and fo₂. It is thus important to constrain atmospheric temperature from observations to improve the retrieval of fo₂. Generally speaking, the thermochemical timescales are shorter at increasing temperatures, pressures and fo₂. Below 700−800 K, the timescale can be very long depending on O/H and pressure (Fig. 4). At these temperatures, equilibrium chemistry might not dominate over photochemistry and both must be considered carefully. The generalization of the timescales for photochemistry is difficult as it depends on several parameters including atmospheric composition, temperature, pressure, orbital distance and the distribution of stellar flux across the 1D column atmosphere. In this section we perform 1D modeling including climate and chemical calculations (see Table 2) to understand the combined role of photochemistry and thermochemistry, and their influence on the molecular abundances of C–H–O species.

4.1 Production of methane in volcanic atmospheres: A chemistry-climate feedback

Previous outgassing calculations showed that methane-rich degassing only occurs at low melt temperature and high degassing pressure (Wogan et al. 2020; Tian & Heng 2024). However, methane can be produced by atmospheric cooling (thermochemistry) using the degassed CO and H₂ (Liggins et al. 2023). CH₄ is a strong greenhouse gas that can significantly affect the temperature profile. For our simulation with a surface pressure of 0.01 bar, we found CH₄ buildup by atmospheric cooling to be inefficient given the low surface temperature around 650 K. CH₄ is produced efficiently in the 1 and 100 bar simulations however. Figure 5 shows the evolution of temperature and CO–CH₄ mixing ratios during the process of atmospheric cooling for the two end-member fo₂ (IW-3 and IW+3) with a surface pressure of 1 bar. We focus on this 1-bar surface pressure simulation to describe the mechanism forming methane. In the low fo₂ scenario (IW-3), the initial atmosphere presents a surface temperature much cooler than the high fo₂ scenario. An oxidized atmosphere made of CO₂−H₂O indeed possesses a higher IR opacity compared to the reduced CO–H₂ atmosphere. The 900 K initial surface temperature of the low fo₂ scenario is explained by the high abundance of water in the atmosphere with a smaller contribution of H₂−H₂ CIA. The surface temperature of the low fo₂ atmosphere, however, increases by 300 K in only a few years, whereas the high f₂ surface temperature is stable over time (top panels, Fig. 5). Both simulations converge to a similar final surface temperature, although warming is driven by CO₂ and H₂O in the high fo₂ case, and by H₂O and CH₄ in the low fo₂ case.

Figure 5 (left panels) shows that the increase in surface temperature in the low fo₂ scenario is correlated with the increase in CH₄ and decrease in CO mixing ratios. CO–CH₄ conversion indeed favors methane at a temperature of 900 K, and the chemical timescale is still short (Fig. 4), thus allowing fast methane buildup. At 900 K, chemical equilibrium would, however, predict a methane mixing ratio above the CO mixing ratio (see dashed black curve, Fig. 5, middle left panel), whereas our simulation shows that the methane abundance stabilizes around 5.10⁻³. This is explained by the radiative properties of methane limiting its own buildup as the abundance of CH₄ is very sensitive to the temperature. The increase of temperature caused by the increase in methane mixing ratio indeed changes the CO–CH₄ conversion making CO more stable at the final temperature of ≈1200 K. In other words, the strong IR opacity of methane prevents its buildup above percentage level. This statement is true as the radiative timescale is here shorter than the chemical timescale leading to a climate feedback evolving progressively with the chemistry. In the high fo₂ scenario, the higher O/H favors oxidation of CO (R5) making methane buildup less efficient (right panels, Fig. 5). Liggins et al. (2023) previously showed that high methane abundances can be obtained in low fo₂ conditions by atmospheric cooling, although their simulations do not consider this climate feedback which likely lead to biases in the final abundances. In both fo₂ scenarios, equilibrium chemistry would predict high methane abundances in the upper atmosphere, at pressures below 0.1 bar, given the lower atmospheric temperature (middle panels, black dashed curves, Fig. 5). Vertical mixing, however, quenches the methane abundance at higher pressures making the observed abundances dominated by disequilibrium chemistry. In the case of an oxidized rocky interior, for instance, the abundances of CO and CH₄ in the lower atmosphere are at chemical equilibrium (Fig. 5, middle & bottom right panels), whereas their abundances in the upper atmosphere are homogenized, a process referred to as mixing-driven quenching. This occurs as mixing become more efficient (shorter timescale) than chemical equilibrium in the upper atmosphere following the decrease in temperature and pressure (Moses 2014; Tsai et al. 2017). The quench point refers to the altitude (or pressure) at which the timescale of mixing becomes equal to that of chemical equilibrium. In other words, one must know the pressure-temperature conditions in the lower atmosphere close to the surface to interpret the methane abundance quenched by mixing in the upper atmosphere.

In the high surface pressure scenario (100 bar) of the low fo₂ end-member simulation, methane is more stable and can build up to higher abundances reaching percentage levels (not shown). In general, the methane buildup and lower ratio CO/CH₄ suggest the presence of a thick atmosphere leading to high temperatures at the surface. High-pressure/High-temperature chemistry dominates the observables as thermochemistry and mixing-driven quenching is faster than photochemistry in these conditions. K_zz does not affect the final outcome of the feedback in terms of temperature and mixing ratio. It however affects the evolution over time as the mixing and radiative timescales are competing. For low K_zz values, the CH₄ buildup starting in the lower atmosphere at higher pressures is not transported upwards fast which affects the radiative feedback. In the 100 bar scenario, the chemical timescale is shorter than the radiative timescale making the model converge fast, with fewer iterations. The dominant chemical mechanism changes for lower pressures where photochemistry dominates the evolution, this is discussed in Sect. 4.2.

If one would change the planet’s orbital distance, the resulting effect on the surface temperature would also modify the outcome of the climate-chemistry feedback induced by methane. This process is therefore sensitive to the atmospheric redox state and surface temperature. This latter parameter is itself controlled by the abundances of greenhouse gases, the surface pressure and the orbital distance.

4.2 Atmospheric CO₂/CO as a tracer of the interior fo₂: Effect of thermochemistry and photochemistry

As discussed in Sect. 3, CO₂/CO decreases with increasing temperature driven by the change in OH/H. In a 1D column atmosphere, the CO–CO₂ abundances will change throughout the atmosphere given the decrease of temperature with altitude. The vertical abundances of the main species involved in the net reaction CO + H₂O ⟷ CO₂ + H₂ are shown in Fig. 6 for the two end-member fo₂ (IW-3 and IW+3) simulations considering a 1-bar atmosphere. Since CO and CO₂ quench at a lower pressure than CH₄, the change of temperature in the lower atmosphere is important and will control their abundance in the upper atmosphere (i.e., in the region probed during spectroscopic observations). For both fo₂ scenarios, we observed that CO₂/CO increases with altitude caused by the decrease in temperature (Fig. 6). K_zz controls the quench point and thus the CO₂/CO in the upper atmosphere.

In our simulations, the CO₂/CO predicted by outgassing can change in the atmosphere by a factor of ≈ 2−3 depending on the temperature and mixing rate (see Fig. 6). In theory, CO₂/CO could change by several orders of magnitude (Fig. 2) as discussed in Sect. 3 but the long chemical timescales obtained for the low temperature simulations make thermochemistry inefficient compared to photochemistry. This change of CO₂/CO by a factor of 2−3 is reasonable and suggest that the molecular abundances in the atmosphere still encode crucial information on the redox state of the interior.

Although thermochemistry dominates the vertical abundances of the species in the 1-bar and 100-bar simulations given the high surface pressure and temperature, photochemistry dominates the evolution in the 0.01-bar simulation. The temporal evolution of the temperature profile and vertical mixing ratios of CO₂, CO and H₂O are shown in Fig. 7 for a 0.01-bar atmosphere produced around GJ 1132b from a reduced rocky interior (IW-3). Analysis of the chemical network reveals an evolution similar to thermochemistry with CO + H₂O ⟷ CO₂ + H₂ as a net reaction: $${\mathrm{H}}_2\mathrm{O}+\mathrm{h}\nu ⟶\mathrm{O}\mathrm{H}+\mathrm{H}$$(R16) $$\mathrm{C}\mathrm{O}+\mathrm{O}\mathrm{H}⟶{\mathrm{C}\mathrm{O}}_2+\mathrm{H}$$(R17) $$\mathrm{H}+\mathrm{H}+\mathrm{M}⟶{\mathrm{H}}_2+\mathrm{M}$$(R18)

net: CO + H₂O ⟶ CO₂ + H₂

The change here is that radicals OH and H are produced by photo-dissociation in the upper atmosphere instead of thermo-dissociation in the lower atmosphere. In the low fo₂ scenario (IW-3), the abundance of CO exceeds that of water. The high abundance of water, however, leads to an efficient production of the OH radical which oxidizes and destroys CO. Over long periods of time (tens to hundreds of thousands of years), the loss of water and CO becomes significant with the oxidation of CO leading to an increase of CO₂ abundances. The photolysis of CO₂ reproduces CO but this process is less efficient than the water photolysis and oxidation of CO. As the system evolves over time, water mixing ratio decreases and CO will stabilize to an abundance resulting from a balance between destruction by oxidation and production by CO₂ photolysis. Overall, CO₂/CO increases by ≈1 order of magnitude in tens of thousands of years (Fig. 7). The temperature profile changes over time but only in the upper atmosphere (Fig. 7, top left panel) since the initial stratospheric inversion induced by water (Malik et al. 2019b) is removed as a result of its photochemical loss. This photochemical mechanism suggests that photochemistry can introduce a bias in CO₂/CO, although the long timescale of this process would make it inefficient if the planet is still volcanically active.

4.3 Overall influence of pressure and fo₂ on atmospheric abundances

The chemical regime controlling CO₂/CO, i.e. photochemistry or thermochemistry, is mainly influenced by oxygen fugacity itself (through varying O/H), orbital distance and surface pressure. All of these parameters influence the surface temperature and thus the kinetic competition between photochemistry and thermochemistry. Since the stellar flux and orbital distance are well constrained for the different rocky exoplanets observed, the main parameters one should assess carefully to interpret relative molecular abundances are surface pressure and fo₂. Mixing can be important, although it should affect the quenched relative abundances by a factor of 2−3 at the most (Fig. 6). If the atmosphere is thin (low surface pressure and temperature) and continuously replenished by geochemical outgassing, the CO₂/CO will be preserved given the long timescales of the photochemical mechanism (see Sect. 4.2).

For each scenario, we calculated with HELIOS the contribution function at different spectral bins to assess the pressure and temperature probed in emission. The contribution function quantifies the weight of each atmospheric layer in the emission of radiation escaping at the top of the atmosphere. In spectral bins with high molecular opacity, the radiation escaping at the top of the atmosphere mainly originates from the low-pressure layers of the upper atmosphere. In a spectral range with no clear molecular opacity, deeper layers are probed. The contribution function at different spectral bins corresponding to transparent windows and absorption peaks of CO₂ and CH₄ mainly is shown in Fig. 8 (top right panel) for the high fo₂ (IW+3) 1-bar scenario. We clearly observe that the deep atmosphere is probed at 3.9 μm, between the absorption peaks of CO₂ and H₂O.

Figure 8 (top left panel) shows the emission spectra (expressed as ratio of flux emitted by planet relative to that emitted by the star) for the two different fo₂ end-member scenarios assuming a 1-bar atmosphere. For the high fo₂ scenario (IW+3), emission between 3.2 and 4.2 μm is high as we are probing the deeper atmospheric layers in agreement with the information provided by the contribution function (Fig. 8, top right panel). For the low fo₂ scenario, we see a significant difference in this spectral range as methane is absorbing around 3.4 μm. In other words, the difference between the two fo₂ end-members in emission is mainly reflected between 3.2 and 4.2 μm where the presence or absence of methane leads to a change in the temperature/pressure probed. In the case where methane is not present, i.e. high fo₂, the temperature of the deep atmosphere could be inferred from these spectroscopic observations. We remind the reader that the absence of methane can also be explained by a thin atmosphere with low surface temperature in a low fo₂ scenario.

In the region between 10 and 15 μm previously used to assess the presence of a thick CO₂ atmosphere around Trappist-1c (Zieba et al. 2023) and −1b (Ducrot et al. 2024), we observe a change between the low fo₂ and high fo₂ scenarios. In both cases, the 15-μm CO₂ band is observed but the ratio of signal between 12 and 15 μm is different and traces the abundances of water. This ratio will likely depend strongly on H₂O/CO₂ and therefore does not directly reflect fo₂.

Figure 8 (bottom panels) shows the simulated transmission spectra of GJ 1132b for a 1-bar atmosphere and two different fo₂ (IW-3 and IW+3). In transmission, variations between low fo₂ and high fo₂ scenarios mainly lie in the abundance of methane and carbon monoxide as shown in Fig. 8 (bottom panels). In the low fo₂ case, the CO–H₂ atmosphere is seen with the clear signatures of methane at 2.3 and 3.4 μm and CO dominating the contribution at 4.7 μm. CO₂ is still abundant and detected in the low fo₂ case using its 4.4 μm feature. In the high fo₂ scenario (bottom right panel, Fig. 8), the high abundances of CO₂ and H₂O dominate the spectrum and contribute significantly at 4.7 μm given the lower opacity and abundance of CO. Retrieving accurately high CO₂/CO might be challenging as the data would probably be fitted efficiently without considering the contribution from CO. Future work should assess in detail the limitations in inferring CO₂/CO for rocky exoplanet atmospheres using a retrieval framework. Only with this detailed retrieval analysis, one could evaluate how accurate the CO₂/CO value can be.

In our analysis, we observed that the probed pressure at the peak of the opacity bands of CO₂, CO, CH₄, and H₂O is below 0.1 bar (top right panel, Fig. 8). The probed mixing ratio thus corresponds to the quenched value controlled by vertical mixing. Our analysis suggest that the effect of K_zz has negligible impact on the transit spectra as the variation in mixing ratios are small for the dominant atmospheric species and usually within 50% for the secondary species (CO–H₂ in high fo₂ case and H₂O−CO₂ in low fo₂ case). In Fig. 9, we review the quenched abundances of the different simulations across the P_s−fo₂ space and summarize the changes in CO₂/CO, caused by atmospheric cooling for a 1-bar and 100-bar atmosphere and photochemistry for a 0.01-bar atmosphere, compared to the outgassing prediction. We also assess the changes in CO₂/CH₄ as the distinct signatures of methane in transit spectra are easier to resolve in comparison to CO (Fig. 8, bottom panels).

In all the scenarios, CO₂/CO is increased (Fig. 9, top panel) by atmospheric chemistry compared to the outgassing prediction following the chemical mechanisms described in Sects. 3.2 and 4.2. With thermochemistry, the change in CO₂/CO remains below one order of magnitude as the high abundances of greenhouse gases in volcanic atmospheres (CH₄, CO₂ and H₂O) tend to produce high surface temperatures with a relatively low difference with the melt temperature. Photochemistry can modify CO₂/CO for thin atmospheres, although this process occurs over timescales likely longer than the replenishment by outgassing. These two mechanisms changing CO₂/CO in the atmosphere are summarized in Fig. 9 (bottom panel). Following the outgassing calculations, CO₂/CO changes by one order of magnitude when fo₂ varies by two orders of magnitude (see Fig. 9, top panel). The effect of atmospheric chemistry introduces a limit into how accurate one can be in inferring fo₂ if the temperature of the quench point cannot be constrained. We estimate this bias around 1–2 orders of magnitude.

Figure 9 (top panel) shows that CO₂/CH₄ decreases following the thermochemical buildup of methane induced by atmospheric cooling for thick atmospheres (1-bar and 100-bar cases) and described in Sect. 4.1. At fo₂ around IW and below, CO₂ and CH₄ coexist at high abundances. CO₂/CH₄ can help constrain atmospheric fo₂ if the CO abundance cannot be retrieved accurately. We emphasize that unlike CO₂/CO and fo₂, the correlation between CO₂/CH₄ and fo₂ is not linear (see Fig. 9, top panel) as it also depends on the hydrogen budget; however, it provides a crucial insight into the redox state. Above IW, CO₂/CH₄ is very high, and methane will likely not be observed given its low stability in an atmosphere with high O/H (see Fig. 8, bottom right panel). The high abundance of water and CO₂ and the absence of methane might be the best indicator of high fo₂ scenarios, although constraining the exact fo₂ value is limited by the detection efficiency of CO. At low fo₂, the co-existence of CO, CH₄, and CO₂ could be constrained and provide crucial information on the interior fo₂ and the atmospheric temperature. CO/CH₄ can help constrain the atmospheric temperature (Fig. 2, top-right panel), whereas CO₂/CO or CO₂/CH₄ can be used to constrain fo₂ (Fig. 9, top panel). Figure 9 also shows that the difference in CO₂/CO and CO₂/CH₄ between the 1-bar and 100-bar scenarios is small, below the accuracy expected with JWST. This tells us that these relative molecular abundances are not clearly correlated to the surface pressure. Surface pressure is, however, one of the main parameters affecting the outcome of atmospheric chemistry as it also controls temperature. A change in pressure from 1 bar to 100 bars does not significantly change the methane abundances. Methane is more stable at high pressure, but it is also less stable at a high temperature. Although a thick atmosphere has the ideal pressure conditions for methane to be stable, it also creates warmer conditions that will tend to de-stabilize methane. In the end, the methane abundance is not changed significantly compared to a 1-bar atmosphere. Since retrieving P_s from observations is difficult (e.g., Heng & Kitzmann 2017), retrieval framework and fo₂ inference should consider this parameter carefully.

5 Model limitations and discussions

5.1 Limitations of the model

The geochemical outgassing calculations following the approach of French (1966) provide a simple solution for the partial pressure of volatiles in the C−H−O fluid system, although it does not consider their solubility in a silicate melt. In practice, one needs to extend the simple system in French (1966) to include gas-melt equilibria, for example for water (dissolved as OH⁻) and CO₂ (dissolved as carbonate ${\mathrm{C}\mathrm{O}}_3^{2-}$), following the formalism in Iacono-Marziano et al. (2012), Gaillard & Scaillet (2014), Burgisser et al. (2015) or Liggins et al. (2020). Ascension of magma and surface degassing is typically simulated by solving the set of gas-gas and gas-melt equilibria equations at decreasing pressure until soluble species ex-solve from the melt and the surface pressure is reached. This approach is used in the D-compress and EVolve models of Burgisser et al. (2015) and Liggins et al. (2022) respectively, both used for applications in planetary and exo-planetary sciences.

The added complexity in these calculations leads to additional model parameters including mantle melt fraction, partition coefficients of the species and volatile budget, all of which are largely unconstrained for exoplanets. Our approach focuses on the main parameter affecting gas speciation, oxygen fugacity, using simple calculations that can be easily included in the analysis of JWST data to assess the error on the retrieved fo₂. The solubility of water could, however, affect the O/H budget of the atmosphere and change the outcome of thermochemistry and photochemistry. Water is indeed favored at higher pressures in the gas phase (Tian & Heng 2024), although water solubility becomes important above 1 bar (Gaillard & Scaillet 2014; Bower et al. 2022). Since the change of CO₂/CO caused by atmospheric cooling and photochemistry is directly correlated to O/H and the abundance of water, solubility must be considered carefully when applying this approach to warm rocky exoplanets. For instance, the increase in CO₂/CO for thin atmospheres shown in Fig. 9 (top panel) and driven by photochemistry would change in an atmosphere depleted in water vapor. In a similar way, atmospheric escape can change the O/H atmospheric budget over time. Given the conclusions of this study, both melt solubility and atmospheric escape should be considered in future work to assess in greater details their influence on the variations of O/H compared to fo₂. This complexity increases the difficulty in inferring the interior fo₂ from observations especially since the surface pressure (controlling solubility) and the escape history are not constrained from observations.

Given the high solubility of N in the melt at low fo₂ (Bernadou et al. 2021) and its significant impact on the atmospheric abundances of NH₃ (Shorttle et al. 2024), we leave these calculations for a future study where solubility will be included in the calculations. N-rich volcanic atmospheres are expected (Liggins et al. 2022, 2023), one might foresee a similar effect of O/H on the N₂-NH₃ conversion with N₂ favored at high O/H. The effect of NH₃ on the radiative-chemical feedback discussed for CH₄ (Fig. 5, Sect. 4.1) should be assessed given its high IR opacity and its thermochemical stability also strongly dependent on temperature.

5.2 Implications for JWST observations

For rocky exoplanets, recent observations of 55 Cancri e (Tsiaras et al. 2016; Hu et al. 2024), GJ 486b (Moran et al. 2023) or GJ 1132b (Swain et al. 2021; May et al. 2023) suggest the presence of an atmosphere where the main gas components can be identified. Although the observations of Trappist-1c (Zieba et al. 2023), Trappist-1b (Ducrot et al. 2024) and LHS 475b (Lustig-Yaeger et al. 2023) point to a thin atmosphere or no atmosphere, the data quality suggest that a thick atmosphere should be detected if present. Even in emission, one should expect to resolve CO₂ and H₂O bands efficiently if enough transits were accumulated (Piette et al. 2022). The current observations are limited by the spectral range and the presence of cool stellar spots leading to detections of methane and water bands that might not originate from the planet’s atmosphere (May et al. 2023; Moran et al. 2023). Despite these limitations, interpretation of molecular detections in rocky exoplanet atmospheres is starting including reasoning around the interior redox state (Heng 2023; May et al. 2023; Hu et al. 2024). The variations in the secondary eclipse observations of 55 Cancri e could be caused by transient degassing unbalanced with atmospheric escape given the very short orbital period of the planet (Heng 2023). Recent observations suggest an atmosphere rich in CO₂ and/or CO as revealed by the clear absorption between 4 and 5 μm (Hu et al. 2024). The large error bars on the data, however, do not allow for clear constraints on CO₂/CO and fo₂ (Hu et al. 2024). Our theoretical simulations do not indeed predict a clear distinction between high fo₂ and low fo₂ between 4 and 5 μm (Fig. 8, top left panel). Additional measurements of 55 cancri e around the 15-μm CO₂ feature with MIRI might help constrain the main atmospheric species and the ratio CO₂/CO. For GJ 1132b, the detection of H₂O and CH₄, if not caused by stellar spots, limits the interpretation of the redox state as CO₂ and CO abundances are not well constrained (May et al. 2023). The recent data suggest a H₂O-dominated atmosphere which could be explained by a thin atmosphere produced from an oxidized interior. H₂O-dominated atmospheres are unlikely at high pressures because of the high solubility in melts (Bower et al. 2022; Sossi et al. 2023) and recent emission observations of GJ 1132b rule out the presence of a thick atmosphere (Xue et al. 2024).

Our work is also relevant for sub-Neptune planets as recent observations suggest that a magma ocean could be present beneath the thick H₂ envelope with degassing explaining a high metallicity (Shorttle et al. 2024; Benneke et al. 2024). In addition, the high atmospheric pressure scale height of these objects make them ideal targets with JWST. The high abundances of CH₄ and CO₂ inferred from JWST observations of K2-18b, although compatible with biotic methane emissions in an ocean planet with a shallow H₂ envelope and habitable surface temperature (Madhusudhan et al. 2023), could also be explained by a thick H₂ atmosphere with high metallicity around 100 Xsolar (Wogan et al. 2024). The co-existence of CO₂-CH₄ at high abundances in the atmospheres of K2-18b may be explained by a deep magma ocean with reducing conditions at equilibrium with the thick H₂ atmosphere (Shorttle et al. 2024). Following our study, the high abundances of methane in the atmosphere of K2-18b point to a low OH/H, although one would need to constrain the atmospheric temperature to retrieve the fo₂ with greater accuracy. At these high metallicities and low oxidation states, abundances of CO higher than CO₂ are expected (Shorttle et al. 2024) but not suggested by the retrieval analysis of Madhusudhan et al. (2023). The difference between the habitable shallow atmosphere and warm deep atmosphere scenarios lies in the CO abundance relative to CO₂ and CH₄. As described in Sect. 4.3, CO could be unconstrained given its low IR opacity and its main band overlapping with that of CO₂. CO₂/CH₄ can, however, be seen as a tracer of fo₂ as well (see Sect. 4.3). More quantitatively, Madhusudhan et al. (2023) and Benneke et al. (2024) suggested CO₂/CH₄ ≈ 1 for both K2−18b and TOI-270d respectively which would roughly correspond to a fo₂ around IW-3 (see Fig. 9) if we consider atmospheric cooling. A more complex model considering equations of state is needed for sub-Neptunes to accurately consider the high pressure conditions at the boundary between atmosphere and rocky body. The nature of these sub-Neptune planets is still unclear, although the composition of the atmosphere may be hybrid, made of primitive H₂-He mixed with heavier gases from volcanic degassing (Benneke et al. 2024). One would need to assess the CO₂/CH₄ and if possible CO₂/CO for several sub-Neptune planets with different equilibrium temperatures to confirm the deep magma ocean hypothesis and help unveil the nature of these mysterious objects.

6 Conclusions

Using relative abundances of the simple molecules (CO, CO₂., CH₄, and H₂O) that are starting to be constrained from observations of rocky exoplanet atmospheres with JWST, the implications for the redox state of the rocky interior can start to be discussed. The relative molecular abundance CO₂/CO provides crucial information given the direct correlation with oxygen fugacity. This relative abundance can, however, be modified in the atmosphere by thermochemical cooling (referred to as atmospheric cooling) and photochemistry. Our numerical simulations revealed the following:

• Thermochemical cooling leads to an increase of CO₂/CO by up to one order of magnitude, which would bias the inferred fo₂ by two orders of magnitude. This process is strongly dependent on the temperature difference between the melt and the quench point in the atmosphere. Constraints on the atmospheric temperature can help correct the bias on the retrieved fo₂;

• Photochemistry leads to an increase of CO₂/CO over long timescales and is driven by the oxidation of CO via the photolysis of water vapor. The bias on CO₂/CO would depend on the competition between outgassing flux and photochemistry;

• The transition between these two chemical regimes (i.e., photochemistry-driven and thermochemistry-driven) is not associated with a single value of temperature. It depends not only on surface pressure and temperature but also on the redox state of the atmosphere (via O/H). Thick atmospheres (high surface pressure and temperature) with high fo₂ tend to be driven by thermochemical cooling, whereas thin atmospheres with low fo₂ tend to evolve controlled by photochemistry.

Although CO₂/CO is probably the most reliable tracer of fo₂, the weak detectability of CO in spectroscopic observations can limit our ability to infer fo₂. CO₂/CH₄ can be used as an alternative tracer, as methane is produced via atmospheric cooling from the outgassed CO and H₂.

Acknowledgements

T.D. acknowledges support by the “ADI 2021” project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. T.D. and N.C. thank the European Research Council for funding via the ERC OxyPlanets project (grant agreement No. 101053033).

References

---

All Tables

All Figures