© The Authors 2025

1 Introduction

In the last 30 years, almost 6000 exoplanets have been detected¹, and atmospheric characterization has started to become possible for a subset of the known population. Due to observational constraints, to date, it has been possible to measure the atmospheres of gas giants, sub-Neptunes, and a few super-Earths spectroscopically. A few selected candidates of approximately Earth’s size have been studied (e.g., Kreidberg et al. 2019), with the Spitzer InfraRed Array Camera (IRAC) and, more recently, with the James Webb Space Telescope (JWST) (see, e.g., Zieba et al. 2023; Ducrot et al. 2024; Xue et al. 2024; Bello-Arufe et al. 2025). However, little is still known about the compositions of terrestrial exoplanet atmospheres or even of their existence. This gap can only be addressed through direct atmospheric characterization via transmission spectroscopy, emission spectroscopy, or reflected light observations (see, e.g., Seidler et al. 2024; Hammond et al. 2025), and near-future JWST observations will be dedicated to this effort.

Gas giants retain nebular primary atmospheres, and their composition is directly connected to their formation history, while rocky planets develop secondary atmospheres, which are acquired through various processes, including accretion of nebular gas (Sasaki & Nakazawa 1990) and devolatilization of impactors (Kuwahara & Sugita 2015), but mainly through outgassing during their initial magma state and subsequent volcanic degassing (see, e.g., Abe 2011; Carlson et al. 2014; Schaefer & Fegley 2017; Gaillard & Scaillet 2014). In this context, we cannot directly connect their atmospheres to the conditions at their formation (Brachmann et al. 2025), at least not without first understanding other intermediate mechanisms.

The geochemical outgassing processes responsible for the birth and evolution of secondary atmospheres occur most intensively in the early stages of planetary evolution, when the interior and surface are hot enough and are often in the magmaticmelt phase. If a planet experiences high tidal energy dissipation, outgassing processes can persist as long as the volatile budget of the interior remains unexhausted, with volcanic activity and volatile cycles playing an important role in the planet’s overall evolution over its lifetime. Rocky exoplanets that are estimated to exhibit high outgassing activity due to volcanism or those covered entirely or partially by magma oceans offer a rare opportunity to explore the conditions under which planets can develop secondary atmospheres. Rocky exoplanets with high surface and atmospheric temperatures offer enhanced observability in both transmission and emission spectroscopy. In transmission, elevated temperatures can partially offset the suppression of spectral signals caused by high mean molecular weight, increasing the scale height. In emission, higher temperatures enhance the thermal contrast with the host star. These effects are especially pronounced for close-in planets orbiting cool stars, such as M dwarfs, where the small stellar radii amplify transit signals and facilitate the characterization of even tenuous secondary atmospheres dominated by heavy molecules using current instrumentation. Empirical evidence is essential for unraveling the complex interplay between a rocky planet’s initial volatile inventory and the atmospheric loss and replenishment processes that shape its evolution. Observations of carefully selected candidates, combined with numerical simulations exploring a broad parameter space, will provide critical test cases for elucidating the observable manifestations of atmospheric processes. Typically, candidates with equilibrium temperatures up to 1000 K are preferred, as their close proximity to their host stars not only ensures high temperatures but also makes them likely to be in a tidally locked configuration, which may result in significant tidal heating (see, e.g., Seligman et al. 2024).

In this context, LP 791-18 d is an Earth-size planet orbiting the cool M6 dwarf LP 791-18 (Crossfield et al. 2019; Peterson et al. 2023) and is part of a coplanar system. This configuration presents a unique opportunity to study an exo-Earth alongside a sub-Neptune that may have retained its gaseous or volatile envelope. Given its potential for atmospheric characterization with JWST, this planet was selected as a case study for this work. Furthermore, having a sub-Neptune in the same system makes it possible to conduct comparative planetology. The gravitational interaction with the sub-Neptune prevents the complete circularization of LP 791-18 d’s orbit, resulting in continued tidal heating of the interior, possibly inducing high volcanic activity. This continuous tidal heating must be balanced by convective cooling through the surface, which can lead to significant outgassing activity. The estimation of the mantle equilibrium temperature from Peterson et al. (2023) for different melt-to-rock fractions ranges from 1630 to 1670 K. In the next section, we demonstrate that this calculation likely underestimates the magma temperature and that the interior thermal equilibrium is expected to reach higher temperatures by a few hundred kelvins. Observations with the JWST during Cycle 3 have been approved for LP 791-18 d (Benneke et al. 2024). The aim of this observation is to achieve a 4σ confidence level in distinguishing between a bare-rock scenario and a CO₂-rich atmosphere capable of efficient heat redistribution, even under conditions that might include Venus-like clouds. Observations are planned for five visits with the F1500W filter and the SUB256 sub-array, a 3 μm band-pass centered at 15 μm, which covers a strong absorption feature from CO₂. The detection of an atmosphere on LP 791-18 d would provide the first evidence that a temperate rocky planet orbiting an M dwarf can sustain one. All of the above make LP 791-18 d a benchmark candidate for the exploration of outgassing processes and the evolution of secondary atmospheres on terrestrial planets orbiting M dwarfs.

2 Approach

To assess the role of interior redox state on the atmospheric stability and observability of LP 791-18 d, we adopted a modeling framework that connects geochemical outgassing and atmospheric processes, and finally to emission spectra and observational diagnostics. In this section, we present a road-map of our approach and outline the steps followed in our theoretical framework.

A key reason for selecting LP 791-18 d as our case study is the ability to estimate its internal thermal state by balancing the energy flux with convective cooling. This calculation is presented in Section 3. Once the internal temperature is determined, it is assumed to correspond to the melt temperature, be it the temperature of surface lavas or the temperature of subsurface magma chambers. This temperature is used to minimize the multi-dimensional parameter space of planetary unknowns, namely the chemical inventory (parameterized by interior oxygen fugacity), surface pressure, melt temperature, and graphite activity. The outgassing rate is parametrized as described in Section 4. The volatile content of the planet is partitioned between the melt and the atmosphere through equilibration, achieved by minimizing the system’s Gibbs energy. Once the gases are released into the atmosphere, we apply full chemical kinetics, photochemistry coupled with radiative convective model, until a steady-state composition is reached. In Section 4, we describe the numerical tools used to perform simulations across a parameter grid defined by our main variables: oxygen fugacity, surface pressure, and graphite activity. Because this parameter space spans a wide range of atmospheric compositions – from highly reducing to highly oxidizing – and from tenuous to massive atmospheres, we aim for a self-consistent approach. For example, we used an analytical approximation for the heat redistribution efficiency and applied mixing length theory to estimate eddy diffusion coefficients. A sensitivity analysis on surface albedo is also performed to test the robustness and broader applicability of our modeling framework under a broad range of surface properties.

Once steady-state atmospheric compositions are obtained, we evaluated their survivability under the background stellar environment in Section 5 by estimating its escape regime and atmospheric mass loss rates. The next step in our analysis is the generation of mock spectra to evaluate observational detectability with JWST at selected photometric bands. We produced two categories of mock spectra. The first includes airless scenarios, representing LP 791-18 d after complete atmospheric escape. For this, we used Hapke reflectance and emittance theory, as described in Section 6. The second category includes the thermal emission spectra of atmospheric scenarios, following the methods in Section 7. We then analyzed all spectra by integrating them over JWST bandpasses and discuss the potential for distinguishing atmospheric scenarios from bare rocks. In addition, we performed parametric radiative transfer simulations that include the effects of photochemical hazes and clouds to assess their impact on our spectral predictions (Section 7.3). Figure 1 presents a concise mind map illustrating the interplay between the various physical mechanisms incorporated into our simulations and the resulting observable quantities.

In Section 8, we discuss the implications of our findings for upcoming JWST photometric observations of LP 791-18 d, and infer the possible current state of the planet. We conclude by outlining broader implications for the population of small, rocky exoplanets and suggest promising targets for future studies.

Fig. 1 Mind map describing the different parts of our methodology. The simulations are composed of three core pieces coupled together. Radiative convective model, chemical kinetics, and an equilibrium inter-phase multi-phase chemistry between the surface and the atmosphere. Post-numerical simulation analysis is done to produce transmission and emission mock spectra.

3 The outgassing nature of LP 791-18 d and its internal thermal state

The possible thermal state of the planet’s interior was investigated in Peterson et al. (2023) by balancing the competing contributions of interior tidal heating and convective cooling. Their calculations were performed using the publicly available melt module² (Barr et al. 2018; Dobos & Turner 2015; Moore 2003).

Their results indicated the existence of stable thermal equilibria at mantle temperatures ~1620–1650 K, i.e., above the presumed solidus temperature of the mantle. Consequently, the planet is expected to harbor a partially molten mantle that evades complete solidification by eccentricity-driven tidal heating balancing the escaping convective flux. We reproduce this family of interior thermal equilibria obtained by Peterson et al. (2023) in our Figure 2, represented by the blue squares. For these equilibria, the surface flux due to tidal heating is ~0.2–0.8 W m⁻², that is, one order of magnitude greater than the heat flux of the Earth (see, e.g., Davies & Davies 2010, which is mainly attributed to its radiogenic activity), and one order of magnitude less than that of Io (see, e.g., Lainey et al. 2009, which is primarily driven by the tidal action of Jupiter’s gravitational field).

These results are based on the assumption that the mantle, at all temperatures, responds to tidal forces as a viscoelastic solid (the specific rheology adopted in Peterson et al. (2023) to model viscoelasticity is that of a Maxwell solid; see Bagheri et al. 2022, for a recent review on solid tidal rheolgies). However, for mantle temperatures beyond the solidus, the increasing melt fraction causes the molten layer to exhibit fluid-like behavior under tidal stresses rather than as a solid. The tidal signature of this rheological transition in molten layers was recently explored by Farhat et al. (2025) by coupling the equations describing dynamical tides in fluid layers with those describing flexure in a viscoelastic solid. The resulting difference in the tidally generated heat flux between the two models, shown in Figure 2, is summarized as follows. In the framework of a solid mantle, the tidal heat flux is controlled by the mantle’s temperature-dependent viscosity: as the temperature increases, more melt is generated, the mantle’s viscosity decreases, and so does the generated heat. This explains the decaying solid tidal flux (black curves in Figure 2) as a function of temperature. In contrast, if the contribution of the fluid within the molten layer is taken into account, at temperatures beyond the solidus, the tidal response increases with increased temperature by virtue of the increased melt fraction.

Consequently, accounting for fluid tidal heating, interior thermal equilibria of LP 791-18 d would shift towards higher temperatures (~1680–1880 K), depending on the profile of convective cooling. The surface heat flux also increases to ~4.5–6 W m⁻², exceeding that which generates the vigorous volcanic activity on Io (2.24 ± 0.45 W m⁻², Lainey et al. 2009). Notably, these new equilibria, similar to those in the solid framework, are dependent on the actual mantle phase diagram, specifically on the radial profiles of the solidus and the critical melt fraction at which the rheological transition occurs. In the evident absence of such constraints, one can confidently argue that the equilibrium temperatures obtained in the solid framework can only serve as lower limits to the thermal state of the partially molten mantle, while the enhanced effect of fluid tides can significantly increase the equilibrium temperatures and the resulting surface flux of the planet. This is of utmost significance for atmospheric dynamics since the outgassed pressure scales exponentially with the melt’s temperature.

Fig. 2 Modeled internal energy balance for LP 791-18 d. Plotted are convective (blue curves) and tidal (black and red curves) heat fluxes as a function of the mantle’s potential temperature. The isothermal mantle solidus is set at Tp = 1600 K. We follow the prescription of Peterson et al. (2023) and use the same parameters therein to compute: i) the convective flux assuming soft turbulence in the mantle (e.g., Solomatov 2015), whereby the four sets of curves cover a range of melt fraction coefficients (B = 10, 20, 30, and 40); and ii) the tidal heating flux of a Maxwell viscoelastic solid mantle. In addition, we compute the tidal heat flux allowing for the fluid contribution of melt in the mantle, following Farhat et al. (2025), using a conservative estimate of the dissipative timescale, σR = 10−3 s−1. Interior thermal equilibria, corresponding to the balance between tidal heating and convective cooling, are shown by the colored markers.

4 Numerical modeling

This work implements a self-consistent and systematic modeling framework for identifying spectra that can provide insight into the formation of LP 791-18 d’s atmosphere and serve as a proxy for rocky exoplanets in general. Unlike other state-of-the-art models of rocky exoplanet atmospheres to date, we model an outgassing-like planet as an integrated radiation-atmospheremelt physical system.

The atmospheric environment in an actively outgassing planetary state is expected to interact with the melt near the surface. The atmospheric composition is determined by the volatile components released through outgassing from the magma’s interior. At this hot lower boundary, multiphase chemical reactions are fast enough to assume equilibrium between the first layers of the atmosphere and the melt (see, e.g., Tian & Heng 2024, and the references therein). It is important to note that, in the volcanic degassing scenario, equating the atmospheric surface pressure to the outgassing pressure implicitly assumes that gas–melt equilibration always occurs at the surface, rather than deeper within the magma chamber. Once this hot gas is released into the cooler environment of the upper atmospheric layers, it can no longer be considered in equilibrium, and a proper chemical-kinetics treatment is required. We solved the chemical kinetics continuity equation as prescribed by Hu et al. (2012) and Tsai et al. (2017), including photochemistry in the upper layers. Hence, the chemical kinetics equation with eddy mixing, molecular diffusion, and production and loss terms for each species is mathematically described as: $$\frac{\mathrm{\partial }{C}_i}{\mathrm{\partial }t}={\mathrm{\Phi }}_i+P_i-L_i$$(1)

where the left-hand side of the equation represents the rate of change in the concentration of species i over time. The Φ_i on the right-hand side represents the vertical transport of species, which is given in more detail below. Finally, P_i and L_i are the chemical production and loss mechanisms of each species per altitude grid. In this work, we utilized the VULCAN model (see, Tsai et al. 2017, 2021) to self-consistently simulate atmospheric chemical kinetics, vertical mixing (including both molecular diffusion and eddy transport), and photochemistry. Specifically, VULCAN solves the full chemical kinetics continuity equation, accounting explicitly for diffusive vertical transport processes alongside chemical production and loss terms for each atmospheric species. VULCAN includes a database of over 700 neutral reactions. The photochemical calculations within VULCAN also include a built-in UV radiative transfer module with multiple scattering to simulate photodissociation processes accurately. Due to its versatility, this powerful model can simulate any solar or extrasolar atmosphere with relatively minor customization. Thermochemistry within the atmosphere strongly depends on temperature, which, in turn, is influenced by background radiation. To address this problem, we coupled radiative transfer with chemical kinetics. The governing equation can be written as: $$\frac{d{I}_{\nu }}{ds}=-{\alpha }_{\nu }{I}_{\nu }+{\eta }_{\nu }+{ϵ}_{\nu }-{\rho }_{\nu }{I}_{\nu }$$(2)

where $\frac{d{I}_{\nu }}{ds}$ represents the rate of change of intensity with penetration s within the atmosphere, I_ν is the radiation intensity at frequency ν, α_ν is the absorption coefficient, η_ν is the emission coefficient, ϵ_ν is the scattering coefficient, and ρ_ν is the reflection coefficient.

The radiative transfer equation is solved with the Helios³ self-consistent radiative transfer code described in Malik et al. (2017) and Malik et al. (2019). Opacities for atmospheric species were taken from the Data & Analysis Center for Exoplanets (DACE)⁴ database. The interaction between different computational schemes was established to study realistic couplings between interior and atmospheric environments. The coupling of the numerical codes is described in Drant et al. (2025).

4.1 Vertical transport

Vertical transport and mixing of atmospheric species is an important component of our framework. The balance between eddy mixing and molecular diffusion determines the location of the homopause layer. This, in turn, defines the spatial region where molecular diffusion acts such that it can modify the mean molecular weight at the base of the possible planetary wind, R_xuv (see Section 5 for details).

For the eddy diffusion coefficient, K_zz, we wanted to employ a self-consistent approximation that will be flexible across the wide range of atmospheric compositions simulated in this work. We based our approach on the Mixing Length Theory (MLT), assuming convection dominates vertical transport below the radiative region. For each layer, we compute eddy diffusion as: $$K_{zz}={\ell }^{2}N,$$(3)

where ℓ is the vertical mixing length and N is the Brunt–Väisälä frequency. The mixing length is assumed to be a fraction of the atmospheric scale height $$\ell =\alpha H,\text{−with−}H=\frac{k_{B}T}{\mu m_{H}g},$$(4)

where α is the mixing length parameter (in this work we assume α = 0.1 across all of our simulations, which was chosen by benchmarking with the background typical atmospheres for Venus, Earth and Mars), c_p is the local specific heat capacity at constant pressure for the atmospheric composition, T is the local temperature, and g is the gravitational acceleration. The Brunt–Väisälä frequency: $$N=\sqrt{\frac{g}{T}(\frac{dT}{dz}+\frac{g}{c_p})}.$$(5)

The specific heat capacity and mean molecular weight $\overline{\mu }$ are determined from the local atmospheric composition. This yields a self-consistent K_zz profile varying with altitude, temperature, and atmospheric composition.

For the molecular diffusion, we followed the methodology below. Since we solve self-consistently the photochemistry from the outgassing species from the lower boundary, we are agnostic to an a priori background dominant atmospheric gas around homopause; for example, Tsai et al. (2021) includes molecular diffusion of species to an a priori selection of background gas. In this work, we have to follow a full multi-component treatment and not assume any dominant background species. This is done by constructing a matrix of binary diffusion coefficients D_ij between all pairs of species as a function of altitude using the Chapman–Enskog formula, corrected by a species-specific Γ_ij factor (Chapman & Cowling 1970): $$D_{ij}=\frac{k_B{T}^{3/2}}{P{\sigma }_{ij}^2{\mathrm{\Omega }}_{ij}},$$(6)

where k_B is Boltzmann’s constant, T is the temperature, P is the pressure, σ_ij is the average collision diameter of species i and j, and Ω_ij is the corresponding collision integral. A correction factor Γ_ij, empirically determined or set to 1.0 by default, accounts for non-ideal interactions (Banks & Kockarts 1973a,b). The binary diffusion coefficients are then used to compute species-specific diffusion coefficients D_i using Wilke’s Rule $$D_i={\left(\sum _{j\ne i}\frac{X_j}{D_{ij}}\right)}^{-1},$$(7)

where X_j is the mole fraction of species j. The effective molecular diffusion coefficient D_eff of the mixture is then calculated as a weighted average $$D_{\mathrm{e}\mathrm{f}\mathrm{f}}=\sum _i{\varphi }_i{D}_i,$$(8)

where ϕ_i is the mixing ratio of species i. The collision diameters σ and collision integrals Ω are taken from standard tabulated values for common molecules (Carey 1988).

Finally, combining both coefficients in the vertical transport, we can introduce them in Equation (14) as $${\mathrm{\Phi }}_i=-C(K_{zz}\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }z}+D_{eff}[\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }z}+f_i(\frac{1}{H_i}-\frac{1}{H})]),$$(9)

with the Φ_i in cm⁻² s⁻¹, $f_i=\frac{C_i}{C}$ is the mixing ratio of species i, and $H_i=\frac{k_{B}T}{m_{i}g}$ is the scale height of species i.

4.2 Atmospheric heat redistribution

An important question in modeling exoplanet atmospheres with 1D models is how to represent the physical processes that govern atmospheric heat redistribution between the day and night hemispheres. A common approach is to modify the planet’s energy budget through a redistribution factor f, such that the dayside brightness temperature T_day is given by (Burrows 2014): $$T_{\mathrm{d}\mathrm{a}\mathrm{y}}^4={\left(\frac{R_{\star }}{d}\right)}^2(1-{\alpha }_B)f{T}_{\star }^4,$$(10)

where R_⋆ and T_⋆ are the stellar radius and temperature, d is the orbital distance, and α_B is the Bond albedo. The factor f encapsulates the efficiency of day–night heat redistribution, ranging from $f=\frac{2}{3}$ for no redistribution (dayside-only re-radiation) to $f=\frac{1}{4}$ for full, uniform redistribution.

In our framework, we explore a wide range of atmospheric compositions with mean molecular weights varying from approximately 2 to 30 amu. To remain physically consistent across this range, we approximate f using the analytic scaling derived by Koll (2022), which has been benchmarked against 3D general circulation models. This scaling combines weak-temperature-gradient theory and heat engine arguments to estimate the circulation strength and, consequently, the heat redistribution efficiency. The expression we adopt for f corresponds to Equation (10) in Koll (2022), and it smoothly interpolates between the thin- and thick-atmosphere regimes. It thus provides a flexible and physically motivated approximation of the redistribution factor across a wide variety of planetary scenarios.

4.3 Degassing

The above governing equations describe the atmosphere. At the lower boundary, at the surface-atmosphere interface, we bring the system into thermochemical equilibrium with an outgassing melt. The model of Tian & Heng (2024) is employed here to calculate outgassing at the lower boundary, taking oxygen fugacity, surface pressure, and graphite activity as inputs, and producing volume mixing ratios of outgassed species (e.g., CO, H₂O) as outputs.

As initial conditions, we prescribed oxygen fugacity, surface pressure, and graphite activity. The graphite activity (a_C) is a hypothetical parameter dictating the carbon content in the outgassing melt, and thus it indirectly affects the carbon content of the outgassed atmosphere. a_C values ranging from 10⁻⁷ to 10⁻³ are explored in this study, with the lower bound corresponding to carbon-poor outgassing and the upper bound representing relatively carbon-rich outgassing Tian & Heng (2024). The temperature is treated as constant for all simulations, with a fixed value of 1720 K, as explained in Section 3. The atmosphere–interior boundary condition described above sets only the starting chemical composition at the atmosphere’s lower boundary, assuming local thermochemical equilibrium at the gas–melt interface. It does not, however, directly impose global-scale equilibrium throughout the entire atmospheric column. The vertical atmospheric composition and thermal profiles are instead computed dynamically, driven by chemical kinetics, photochemistry, molecular diffusion, vertical mixing (eddy transport), and radiative–convective processes. Thus, although our lower-boundary conditions implicitly assume equilibrium gas–melt interactions (applicable both to global magma-ocean scenarios and localized volcanic degassing events), the resulting atmospheric structure naturally emerges from these self-consistent atmospheric processes. We assume the equilibrium occurs either at the surface in active volcanic vents or near the surface inside magma chambers just before degassing. Due to the high melt temperature assumed (1720 K), local chemical equilibrium is a reasonable assumption, and it does not imply a global-scale interface between melt and atmosphere; instead it merely means gas and melt are in equilibrium before the gas escapes melt to contribute to atmosphere growth (e.g. Bower et al. 2022, Shorttle et al. 2024, Tian & Heng 2024).

A constant surface gas emission flux is prescribed in each simulation until convergence to a steady-state atmosphere is reached. This lower boundary source term is included in the chemical continuity equation as a fixed influx of species, scaled by their mixing ratios in the outgassed composition. Our approach follows the simplified parameterization described by James & Hu (2018), in which the total surface outgassing flux f_j of a species j is expressed as: $$f_j=\dot{V}\cdot {\chi }_j,$$(11)

where $\dot{V}$ is the global magma production rate (in m³ s⁻¹), and χ_j is the mole fraction of species j released per unit volume of erupted magma. This expression is used to fix the lower boundary condition for the chemical continuity equations and is implemented in our framework as described in Tsai et al. (2017). Although $\dot{V}$ is treated as constant during the simulations, it can be physically motivated by planetary parameters. The internal energy release from tidal heating provides an estimate of the total energy available to drive mantle melting. Assuming that a fraction η of this heat contributes to melt production, and using the latent heat of fusion L and the mantle density ρ, the magma production rate can be estimated as: $$\dot{V}=\frac{\eta \dot{Q}}{\rho L}.$$(12)

For LP 791-18 d, $\dot{Q}$ and related parameters are taken from tidal dissipation models provided in Section 3. We estimate a magma production rate of $\dot{V}$ ≈ 20–27 km³ yr⁻¹. This value is comparable to the volcanic production rates observed on Earth and Venus, and lies within the plausible range considered in exo-planetary volcanic outgassing studies (e.g., Gaillard & Scaillet 2014; James & Hu 2018).

Once equilibrium is established between the melt and the atmosphere, numerical simulations are executed until a steady state is reached, accounting for thermochemistry and photo-chemistry.

The stellar input flux was taken from the MUSCLES⁵ catalog. MUSCLES is a spectral survey of M and K exoplanet host stars covering wavelengths from X-ray and ultraviolet to infrared (Behr et al. 2023, Wilson et al. 2025).

Each simulation reached statistical steady-state conditions after the change in parameters reached a threshold, where the minimum limit for convergence was set to ΔU ≈ 10⁻⁶, which was computed at each time step as: $$\mathrm{\Delta }U=\sqrt{\frac{1}{N_{var}}\sum _{u=1}^{N_{\text{var−}}}\left[\frac{1}{N_{grid}}\sum _{grid=1}^{N_{\text{grid}}}\frac{{(U_u^{n+1}-U_u^n)}^2}{{\left(U_u^n\right)}^2}\right]}$$(13)

where N_var is the number of variables, N_grid is the number of grid cells, and U is the u^th variable vector at time n. Examples of the converged atmospheric volume mixing ratios and thermal structure are shown in Figures 3 and 4, respectively.

In Table 1, the simulation grid is presented in the oxygen fugacity and surface pressure space for a_c = 10⁻⁶. We choose fO₂ − IW values ranging from −6 to 5.0, where IW represents the pressure – and temperature-dependent iron-wüstite buffer for oxygen fugacity Tian & Heng (2024), and log f_O2 – IW is the offset from this reference buffer in log₁₀ units. We choose surface pressures ranging from 0.1 bar to 100 bar. In Table 1, we present the mean molecular weight (in amu) of the total atmospheric mass for each steady-state solution, where $$\overline{\mu }=\sum _i{x}_i{\mu }_i$$

accounts for all atmospheric layers. We note a smooth logarithmic dependence of the mean molecular weight, $\overline{\mu }\sim \log \left(f{O}_2\right)$, transitioning from values near those of nebular compositions in highly reduced atmospheres to values resembling those typical of rocky solar system planets. In Figure 5, the mean molecular weight is presented for a combination of parameters. Oxygen fugacity versus surface pressure and graphite activity a_c. Additionally, the graphite activity as a function of surface pressure. We visualize those parameters to allow for intercomparison of their effect on the outgassed mean molecular weight. Surface pressure affects the outgassed mean molecular weight only in oxidized interiors, where [fO₂ − IW] ≳ 1, which generally decreases with pressure. In contrast, changes in graphite activity affect the mean molecular weight across all redox conditions except for highly reduced atmospheres with graphite activity ≲10⁻⁵. In the following section, we highlight the importance of mean molecular weight in determining atmospheric stability. All examined atmospheric scenarios demonstrate that, given sufficient interior volatile inventories, the planet can generate a wide range of secondary atmospheres, in contrast to the solar system’s dichotomy of high-mean-molecular-weight atmospheres for rocky planets and low-mean-molecular-weight atmospheres for gas giants. Once we obtain an atmosphere in steady-state conditions, we generate synthetic spectra for emission viewing geometries. Section 7 provides details on this process. Finally, we convert our synthetic spectra into observational predictions for remote sensing observations with JWST to evaluate the feasibility of future detections under various atmospheric scenarios and to identify which spectral features could theoretically be detected using emission spectroscopy. For this transformation, we utilize the PANDEXO⁶ version software, developed to optimize JWST and Hubble Space Telescope (HST) observations, as described in Batalha et al. (2017).

Fig. 3 Converged simulations for a surface pressure of 1 bar. Volume mixing ratios (VMR) of important greenhouse gases that significantly contribute to emission and transmission spectroscopy are shown. Three representative oxygen fugacity cases are presented to illustrate the transition from a reduced to an oxidized atmosphere.

Fig. 4 Converged simulations for a surface pressure of 1 bar. The thermal structure is shown for all simulated scenarios within the 1 bar surface pressure grid. We note the thermal inversion in the upper layers of the simulated atmosphere, caused by infrared absorption due to high water vapor abundances.

Fig. 5 Mean molecular weight (μ) of the outgassed atmosphere (in amu) as a function of varying surface pressure, oxygen fugacity (fO2), and graphite activity (aC). Each panel isolates the effect of one parameter by holding it fixed while allowing the other two to vary, as indicated in each panel. The value of μ shown corresponds to mean molecular weight of the volatile partition at a given degassing rate. This figure illustrates how redox state, pressure, and carbon content influence the dominant molecular species and therefore control the overall atmospheric composition.

4.4 Surface albedo

For our baseline simulations, we adopted a surface albedo of 0.1, consistent with volcanic material expected for LP 791-18 d, given its anticipated high volcanic activity. To assess the sensitivity of our results to this assumption, we performed a series of simulations with albedo values ranging from 0.01 to 0.6 for both highly reduced and highly oxidized atmospheres. We explored how variations in surface albedo affect the atmospheric thermal structure, considering this broad range of plausible values representative of rocky exoplanet surfaces. Figure 6 shows the resulting temperature–pressure profiles for both redox scenarios, assuming a surface pressure of 1 bar and a fixed graphite activity of a_c = 10⁻⁶. The greatest sensitivity to albedo occurs in the atmospheric layers closest to the surface. In particular, surface-near temperatures vary from approximately 700 to 800 K for the oxidized case and from 550 to 600 K for the reduced one as albedo increases. However, at higher altitudes, the temperature profiles converge regardless of surface reflectivity.

While changes in Bond albedo directly affect the net stellar energy absorbed by the planet, the resulting impact on the atmospheric temperature profile is modulated by the presence of greenhouse gases. In particular, absorption and re-radiation throughout the atmosphere buffer the thermal response to increased stellar flux. In our simulations, although lowering the albedo from 0.6 to 0.01 results in a ~25% increase in equilibrium temperature (as expected from T_eq ∝ (1 − A_b)^0.25), the temperature difference is more pronounced at deeper atmospheric layers. The upper atmosphere, where radiative cooling and shortwave absorption dominate, remains comparatively less sensitive to surface albedo. This is because the greenhouse effect redistributes the incoming energy vertically, and the deposition of stellar radiation does not shift substantially upward with decreasing albedo. The upper atmospheric structure appears relatively stable across the albedo range explored here.

As a result, the impact on the thermal emission spectrum is found to be negligible—at least within the spectral bands considered in this study. This is because the photospheric pressure of those bands, P_photosphere ~ g/κ, lies at relatively low pressures (≲0.1 bar), above the region where albedo makes a significant influence.

Fig. 6 Top: vertical temperature profile for a reduced and oxidized atmosphere at surface pressure of 1 bar and ac = 10−6. We show simulations for surface albedo ranging from 0.01 to 0.6. Bottom: detail from the top plot, focusing on the layers close to the surface.

5 Atmospheric stability

All planetary atmospheres undergo escape, losing species to space, with mass loss rates varying according to planetary parameters such as gravity, atmospheric composition, and temperature structure, as well as the nature of the interaction with the stellar environment. Generally, we can distinguish two main escape regimes. The first is the planetary wind regime, which can sculpt the atmosphere and define its structure, or even its very existence. This regime involves a hydrodynamic flow where gas escapes in bulk to space (Watson et al. 1981; Tian et al. 2005; Owen 2019; Gronoff et al. 2020). The second regime occurs through Jeans escape and non-thermal processes, where species escape due to their ballistic velocities exceeding the escape velocity above the exobase (the non-collisional upper layer of the atmosphere) in a direction opposite to the planet (see, e.g., Shizgal & Arkos 1996). In the latter case, the thermospheric part of the atmosphere, the region where the energetic portion of the stellar spectrum is deposited, can be considered hydrostatic. Here, the incoming energy is balanced by thermal conduction to the lower layers of the atmosphere, where molecular emissions radiate the energy back to space. In the non-hydrostatic scenario, however, the atmosphere balances excess incoming energy mainly through PdV work, leading to a planetary wind regime similar to Parker’s wind mechanism for stellar winds (Parker 1958). In our simulations, we examine the stability of these escape scenarios across the parameter space grid by applying the following methodology. We first determine the steady-state composition and thermal structure of the background atmosphere down to an atmospheric pressure of 10⁻⁸ bar. We then numerically define the R_xuv [m] radius, which corresponds to the atmospheric layer where the optical depth of hard XUV radiation reaches unity, τ_xuv ≈ 1, and assume that this is the layer where XUV energy is deposited. Next, we assume that the system is in an isothermal steady-state wind, a valid assumption if radiative processes occur on timescales shorter than dynamical processes, allowing the gas to regulate at a constant temperature. The spherically symmetric flow is described by Lamers & Cassinelli (1999) and follows the equation: $$\dot{M}=4\pi r^2\rho (r)v(r).$$(14)

Is the continuity equation of the flow, and the momentum equation: $$v\frac{dv}{dr}+\frac{1}{\rho }\frac{dp}{dr}+\frac{GM}{r^2}=0.$$(15)

The energy equation is to be T(r) = T₀ for an isothermal flow, which gives us the following solutions for the so-called critical point where the sonic velocity of the flow is achieved. $$v_s=\sqrt{\frac{k{T}_0}{\mu m_H}}.$$(16)

μm_H is the mean molecular weight in amu. The sonic radius is found by: $$r_s=\frac{G{M}_{\mathrm{p}}}{2v_s^2}.$$(17)

The velocity profile of the wind is then given by: $$\frac{v(r)}{v_s}\exp [-\frac{v^2(r)}{2v_s^2}]={\left(\frac{r_s}{r}\right)}^2\exp (-\frac{2r_s}{r}+\frac{3}{2}).$$(18)

By combining (14) and (18) we find the density profile: $$\frac{n(r)}{n_s}=\exp [\frac{2r_s}{r}-\frac{3}{2}-\frac{v^2(r)}{2v_s^2}].$$(19)

To achieve closure for the previous set of equations and solve for the temperature T₀, we assume an energy-limited constant flow following (e.g., Sanz-Forcada et al. 2011): $$\dot{M}=\frac{\eta \pi F_{\mathrm{X}\mathrm{U}\mathrm{V}}{R}_{xuv}^3}{G{M}_p}$$(20)

where $\dot{M}$ is the atmospheric escape rate [kg/s], η is the efficiency of the escape process, F_XUV is the XUV flux incident on the planet, R_p is the planetary radius, M_p is the planetary mass, and G is the gravitational constant. The efficiency of the escape process, which can range from 0.01 to 0.6 depending on atmospheric composition and irradiation conditions (e.g., Tian et al. 2005; Shematovich et al. 2014; Owen 2019; Watson et al. 1981; Yelle 2004; Salz et al. 2016), is assumed to be η = 0.1. This conservative value is appropriate for secondary atmospheres, typically enriched in heavier molecules and atoms, where efficient radiative cooling processes limit the energy available for escape (Johnstone et al. 2018; Tian et al. 2005).

We estimate the XUV flux reaching LP 791-18 d based on the spectral type and the age of the star as F_XUV ≈ 0.1 [W/m²/s] following: $$\frac{L_{XUV}}{L_{\star }}=\{\begin{array}{ll}{10}^{-3.6},& t<t_{\mathrm{s}\mathrm{a}\mathrm{t}}\\ {10}^{-3.6}{\left(\frac{t}{t_{\mathrm{s}\mathrm{a}\mathrm{t}}}\right)}^{-\alpha },& t\ge t_{\mathrm{s}\mathrm{a}\mathrm{t}}.\end{array}$$(21)

For L_⋆ ≈ L_⊙ × (M_⋆/M_⊙)⁴, t_sat = 100 My, the saturation interval of high-energy flux (Vilhu & Walter 1987, Wright et al. 2011), t ≈ 500 Myr, the estimated age of LP 791-18 (Reiners & Basri 2009, Newton et al. 2018), and α = 0.86 for the full XUV luminosity (see King & Wheatley 2021, their Figure 1), these estimations align with previous studies. We compare our results with Ramirez & Kaltenegger (2014), Johnstone et al. (2015), Johnstone et al. (2021), and Richey-Yowell et al. (2023) to extract F_XUV based on our target’s characteristics. The exobase, R_exo, is approximated by setting l/H(r) ≈ 1, where $l=\frac{1}{n{\sigma }_{\text{col}}}$ is the mean free path, with σ_col = 10⁻¹⁵ cm² as the collision cross-section (Chubb et al. 2024) and H(r) as the local scale height. The stability criterion is then defined at the supersonic point of the flow. In scenarios where the flow becomes supersonic below the exobase, we consider a planetary wind regime with mass loss rates given by Eq. (20). For cases where the critical point is reached in the collisionless environment, the wind cannot be formed, and we assume a stable thermosphere. An example of this analysis is shown in Figure 7 for a scenario with a surface pressure of 1 bar and oxygen fugacity of 0. We plot the velocity profile of an isothermal wind, where the temperature is solved from the system of equations presented in this section. The implications of these calculations for the long-term evolution of the planet’s atmosphere are discussed in Section 8. In scenarios classified as stable (i.e., not experiencing hydrodynamic escape), atmospheric loss occurs primarily through Jeans escape and non-thermal processes. We estimate these to contribute at most ~1–3 kgs⁻¹. At such low escape rates, and assuming a typical volatile inventory for a terrestrial planet, the atmosphere could persist for timescales exceeding ~1 Gyr. The upper panel of Figure 8 shows the steady-state mean molecular weight of the total atmospheric mass across the explored parameter space. The lower panel displays the mean molecular weight near the R_xuw layer, along with the region of parameter space where atmospheric stability is achieved.

Fig. 7 Wind velocity profile as a function of radius for a surface pressure of 1 bar and fO2 − IW = 0. The blue solid line represents the wind velocity profile as a function of radius, the green dashed line indicates the exobase estimation, and the orange dashed line shows the escape velocity, which decreases as ~ r−2. The red dashed line marks the critical radius where the sonic velocity of the flow is reached, while the black solid line represents the Hill radius, beyond which material is no longer gravitationally bound to the planet. In this example, a planetary wind regime is established since rs ≤ Rexo.

6 Thin atmosphere to airless case

Following the atmospheric stability analysis in the previous section, we must consider the possibility of volatile saturation in LP 791-18 d due to continuous hydrodynamic escape throughout the planet’s lifespan. This requires accounting for the scenario of an airless rocky planet or one with a very thin atmosphere (e.g., < 0.01 bar), where there is no significant heat redistribution. In such cases, the observational manifestation of a so-called bare-rock planet depends on the emissivity and reflectance of the surface material. At high surface temperatures (≈1600 K), the planetary surface becomes molten and the near-infrared spectra are largely dominated by blackbody radiation. This behavior is supported from lava flow Earth observations, including laboratory experiments and in situ measurements in Hawaii (e.g., Flynn & Mouginis-Mark 1992, Lombardo et al. 2020). However, at lower temperatures, the crust interacts with electromagnetic radiation through a combination of reflection, absorption, scattering, and emission, which depends on the geometry of the interaction, the physical characteristics of the surface, the radiation wavelength, and micro-crystal formation of the material (Hapke 1993; Shkuratov & Helfenstein 2001; Nelson et al. 2000).

To simulate this interaction and generate synthetic spectra, the optical properties of the surface (e.g., single-scattering albedo) are required. Even for Solar System objects, determining unique photometric parameters from disk-integrated observations alone is challenging, even when wide-phase-angle coverage is available (Domingue et al. 1991). Thus, we must approach the problem inversely, particularly for exoplanet observations.

Since there is no explicit guiding parameter for exploring the parameter space, unlike in atmospheric composition studies where oxygen fugacity serves as a reference, we must infer composition based on our knowledge of Earth, Solar System volcanic activity, and other analogs. Therefore, we assume a surface material and use laboratory measurements of each material to extract its optical properties.

Reflected light depends on the geometry of the system, i.e., phase angle, latitude, and longitude of the surface area, but also on the optical properties of the surface material. Since reflected light and thermal emission may overlap, a unified treatment is required. To address this, we numerically generate a mesh-grid disk where each surface unit is assigned a set of angles, namely, the incident angle i, the emergent angle e, and the phase angle g. The geometric transformation relations between the star-planet system and the system-observer are given by: $${\mu }_0\equiv \cos i=\cos \theta \cos (\alpha -\varphi )$$(22) $$\mu \equiv \cos e=\cos \theta \cos \varphi$$(23) $$g=\alpha .$$(24)

The total incoming light to the observer from the whole disc can then be calculated following Shkuratov & Helfenstein (2001) and Hu et al. (2012) as: $$\frac{F_p}{F_{\ast }}=\frac{1}{F_{\text{inc−}}}{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}{I}_p(\theta ,\varphi ){\cos }^2\theta \cos \varphi d\theta d\varphi \times {\left(\frac{R_p}{D_p}\right)}^2.$$(25)

By default, the stellar coplanar radiation originates from the direction (θ = 0, ϕ = α). The incident radiation at each surface element is given by: $$F_{\mathrm{i}\mathrm{n}\mathrm{c}}=\pi B_{\lambda }\left[T_{\ast }\right]{\left(\frac{R_{\ast }}{D_p}\right)}^2.$$(26)

With B representing the wavelength-dependent stellar blackbody function, R_* the stellar radius, and D_p the star-planet distance. The total incoming radiation at each surface element is a combination of reflected, scattered, and thermally emitted radiation and is expressed in Equation (25) as the I_p function: $$I_p(\theta ,\varphi )=I_s(\theta ,\varphi )+I_t(\theta ,\varphi )$$(27)

where: $$I_s(\theta ,\varphi )=\frac{F_{\mathrm{i}\mathrm{n}\mathrm{c}}{\mu }_0}{\pi }{r}_c({\mu }_0,\mu ,g)$$(28) $$I_t(\theta ,\varphi )=ϵ(\mu )B_{\lambda }[T(\theta ,\varphi )].$$(29)

In Equations (28) and (29), two critical parameters are introduced: r_c, the bidirectional reflectance coefficient or “radiance coefficient” of solid material. This coefficient represents the brightness of a surface relative to the brightness of a Lambert surface under identical illumination (a Lambertian sphere has r_c = 1), and it depends on the direction of both incident and scattered light. Additionally, it is influenced by the chemical composition and crystal structure of the material. Secondly, ϵ, the emissivity, quantifies the deviation of a surface material at temperature T_s from behaving as a perfect Planck emitter, given by B_λ[T_s].

We employed Hapke’s theory of bidirectional reflectance and directional–hemispherical reflectance (Hapke 1981, 2002) to derive the optical and photometric properties of each assumed surface composition scenario. This analytical framework accounts for single and multiple scattering in particulate surfaces and relates emissivity to reflectance through Kirchhoff’s law. We assume isotropic scatterers for multiple scattering and a simplified, cosine-based phase function for single scattering. The key optical parameter in Hapke’s model is the single-scattering albedo, which depends on the relative contributions of absorption and scattering by individual surface particles. The mean directional–hemispherical reflectance is used to compute surface emissivity for energy balance calculations in the absence of an atmosphere. For LP 791-18d, utilizing the system parameters given in Peterson et al. (2023), the secondary eclipse lasts about ≈2° in orbital angles. We then have to account for the opposition surge effect – an increase in reflectance at low phase angles – through a wavelength-dependent enhancement term, S(λ), informed by the microscopic texture and porosity of the surface material (Shkuratov & Starukhina 1999; Helfenstein 1988; Gkouvelis 2025). This approach allows for a consistent treatment of both photometric and thermal properties across the spectrum. The resulting emissivity values are then used to solve for the surface temperature under radiative equilibrium, with incident flux absorbed and re-emitted according to the surface’s wavelength-dependent properties.

Space weathering effects (see, e.g., Hu et al. 2012) will not be considered in the reflected properties of the surface material, since it has already been demonstrated by Peterson et al. (2023) and this work that the interior is highly convective, making crust replenishment likely to occur on geological timescales.

We consider various surface material scenarios to investigate the magnitude of the bare-rock spectrum and compare it with the synthetic spectra of the atmospheric scenarios. Each surface material scenario has its own composition, where the total sum must equal one (i.e., 100%). We simulate metal-rich, ultramafic, feldspathic, basaltic, and wüstite scenarios, which cover the spectrum from a nearly featureless thermal emission to one with strong spectral signatures—e.g., metal-rich and ultramafic, respectively. The general shape and strength of spectral features is in agreement with the works of Hu et al. (2012), Lyu et al. (2024) and Hammond et al. (2025) for each surface scenario. The composition of five basic scenarios is summarized in Table 2.

For each of these five scenarios, we use laboratory measurements from the United States Geological Survey Spectral Library (USGS) (Milliken et al. 2021, Milliken 2020, Clark et al. 2007), which are typically provided for an incidence angle of i = 30°, an emission angle of e = 0°, and a phase angle of g = 30°. We extract the optical parameters required for each model using the Hapke bidirectional model. To solve the inverse problem and retrieve the desired parameters, we apply the Bayesian algorithm described in Appendix B of Gkouvelis et al. (2024), which is based on χ² minimization by evaluating the Jacobians $\frac{\mathrm{\partial }F(x)}{\mathrm{\partial }x}$ of the forward model, which in this case is the Hapke model. The retrieved optical properties for all scenarios are the single-scattering albedo ω, the wavelength-dependent b asymmetry factor which defines the phase function scattering parameter and the photometric parameters relevant to the opposition effect.

The secondary eclipse depth of the system is evaluated using Equation (25) at a phase angle of g = 0°, which is standard practice in observations (e.g., when the planet approaches eclipse, see Benneke et al. 2024) as it maximizes the contrast. The energy balance equation is evaluated numerically, with convergence solutions achieved at angular resolutions finer than ≈7° for both μ and ϕ. We present the synthetic spectra of LP 791-18 d as a bare-rock planet for the five surface scenarios, shown in Figure 9. Additionally, we overplot the isothermal brightness temperature functions T_b(λ) as indicators of the planet’s thermal emissivity (Seager & Lissauer 2010). These brightness temperature functions serve as a visual guide to understanding how different wavelengths trace temperature variations for each scenario. However, it is important to note that these spectra represent disk-averaged values. In practice, this can be determined by solving for the brightness temperature T_b(λ) using: $$F_p=\pi B_{\lambda }\left(T_b\right){\left(\frac{R_p}{D}\right)}^2.$$(30)

For featureless surface scenarios the brightness temperature function T_b(λ) varies by only a few kelvins across the spectral range due to the smooth wavelength dependence of surface emissivity. In contrast, for surfaces with strong spectral features, T_b(λ) can vary by several tens of kelvins, reflecting the stronger wavelength-dependent modulation of the emergent thermal flux.

Fig. 8 Mean molecular weight across the simulation parameter grid space. Top: Total atmospheric mass mean molecular weight (in amu), calculated once the simulations reach steady state. Bottom: Mean molecular weight at Rxuv, the region where the energetic portion of the stellar spectrum is deposited, forming the base of the wind. The thick dashed black line separates the parameter space between the wind regime and stable thermospheres.

7 Observational manifestations

In this section, we focus on observable manifestations that can be used to characterize the planet and derive physical properties, such as atmospheric chemical composition and its correlation with oxygen fugacity, which serves as an indicator of the interior state. The forthcoming JWST observations of LP 791-18 d will target only the CO₂ 15 μm strong absorption band. Therefore, we first analyzed this band based on our simulated spectra. We demonstrate that, based on our simulations, the objectives proposed in Benneke et al. (2024) can be achieved for a limited subset of scenarios within our simulation grid space (see also the discussion section). We then propose alternative strategies for characterizing LP 791-18 d, which may also serve as a general approach for studying similar Earth-size exoplanets.

7.1 Emission

Thermal emission originates directly from the planet and can be detected through its thermal contrast with the host star, expressed as: $$\frac{F_p(\lambda )}{F_{\ast }(\lambda )}=\frac{f_p(\lambda )}{f_{\ast }(\lambda )}\frac{R_p^2}{R_{\ast }^2}\mathrm{\Phi }(\alpha )+{\left(\frac{R_p}{a}\right)}^2{A}_g\mathrm{\Phi }(\alpha ).$$(31)

Where Φ(α) is the orbital phase angle of the planet. It is defined as α = 0 at primary transit, when the observer views the planet’s night side directly, and α = 1 at secondary eclipse, when the planet is behind the host star. The second term represents the contribution from scattered light; more details on its derivation can be found in Seager & Lissauer (2010). Here, α denotes the semi-major axis of the star-planet system, and A_g is the geometric albedo. As discussed in previous sections, our numerical simulations include detailed radiative transfer calculations. As a result, we obtain the outgoing radiation from the planet in 73 spectral bins, which serves as f_p(λ), while the incoming radiation at the top of the atmosphere, f_*(λ), is one of our model inputs. For the latter, we used the stellar spectrum of LP 791-18 from the MUSCLES catalog (Behr et al. 2023; Wilson et al. 2025). This allows us to directly compute the thermal emission spectrum for each of our simulations. Our theoretical framework solves the radiative transfer using the correlated-k method at relatively low spectral resolution (R ~ 70) for computational efficiency. To assess the impact of spectral resolution on the computed emission spectra and the resulting photometric observables, we conducted a numerical experiment to test the robustness of our approach. For a representative subset of simulations, we generated high-resolution spectra using the line-by-line mode of petitRADTRANS⁷ (Mollière et al. 2019, 2020; Alei et al. 2022; Nasedkin et al. 2024) at a resolving power of R = 60 000, which were then convolved to R = 200. We also computed spectra using the correlated-k method at R = 1000, likewise convolved to R = 200. We found that the synthetic photometry derived from these higher-resolution spectra was in excellent agreement with that obtained using the correlated-k method as well as the original HELIOS output. This confirms that the spectral resolution of the HELIOS spectra is sufficient for producing reliable photometric predictions when integrated over the broad MIRI filter bandpasses (R ~ 10–20). A comparison of the results from HELIOS and petitRADTRANS is shown in Fig. 10.

The photometric results presented in this work are based on synthetic spectra generated with petitRADTRANS using the correlated-k method at intermediate resolution (R = 1000), subsequently convolved to R = 200. The HELIOS code was used exclusively to compute the radiative–convective equilibrium structure of the atmosphere and achieve steady-state temperature profiles. These profiles were then passed to petitRADTRANS to calculate self-consistent emission spectra. While petitRADTRANS was used here for its flexibility and higher spectral resolution capabilities, either tool could be used to produce synthetic spectra, and the choice does not impact the scientific conclusions drawn in this work. An example is presented in Figure 11 for different simulated scenarios. We also overplot the isothermal brightness temperature, T_b(λ), as an indicator of the temperature traced at various atmospheric layers at each wavelength (Seager & Lissauer 2010). This follows the relation $P_{\text{photosphere−}}\approx \frac{g}{\kappa }$. In the plotted example, the atmospheric layers are enclosed within a temperature range of approximately 320–500 K.

Due to variations in atmospheric composition, we observe differences in the thermal emission spectra across the oxygen fugacity parameter space. The primary spectral features of interest include: 1) A strong absorption band around 15 μm, which is caused by CO₂ absorption and is therefore more prominent in oxidized atmospheres. 2)An enhanced continuum feature appears around 9 μm, arising from the absence of CH₄ and H₂O opacity. Under more oxidizing conditions, the depletion of these hydrogen-bearing species opens a spectral window, allowing radiation at this wavelength to escape from deeper and hotter atmospheric layers. The strength of the emission feature increases rapidly with oxygen fugacity, owing to the logarithmic nature of the f_O2 scale. As f_O2 crosses a critical threshold (e.g., IW +3), the atmospheric composition undergoes a sharp transition, with the near-complete loss of CH₄ and H₂O. This transition enhances thermal emission at 9 μm, making the feature especially prominent in that regime.

In Figure 12, we show the variation of the 15 μm spectral band (where the MIRI F1500W filter’s effective wavelength is λ = 14.79 μm) integrated over ±1.5 μm to match the MIRI F1500W filter bandwidth (Greene et al. 2023, Hammond et al. 2025). We plot the 15 μm spectral band for each of our simulations, including the airless planet scenarios, using the same scale (right-hand section of Figure 12). We find that bare-rock eclipse depths are comparable to those of reduced atmospheric scenarios, with mean differences of approximately 30–40 ppm relative to oxidized atmospheres. To assess the detectability of these scenarios with JWST’s F1500W filter, we applied our synthetic spectra to the PANDEXO software (Batalha et al. 2017) to estimate uncertainties for each simulation. Using the observational characteristics described in Benneke et al. (2024), we find error bars on the order of ≈ ± 30 ppm, depending on the specific case within our simulation grid.

The estimated level of uncertainty suggests that distinguishing an atmospheric detection from an airless case is possible, but with low confidence, and only in the case of a thin and highly oxidized atmosphere. We discuss this issue in more detail in the final section. In the following parts of this section, we explore alternative approaches beyond solely relying on 15 μm photometric observations.

Firstly, we analyzed the secondary eclipse depth for spectral bands adjacent to the 15 μm band. This observational approach was recently adopted by Ducrot et al. (2024), who conducted follow-up observations of TRAPPIST-1 b in the 12.8 μm band using MIRI’s F1280W filter. We present the ratio of the F1280W to F1500W filters for our simulated scenarios as a function of oxygen fugacity in Figure 13 (upper plot). Furthermore, from the full list of MIRI photometric bands, we identify two bands that correlate with the 15 μm band: The 9 μm band, where CH₄ exhibits strong absorption. The 20 μm band, which remains largely insensitive to oxygen fugacity or surface pressure variations (at least within our numerical experiments).

Across all filter indices, variations of up to 50% are observed in the oxygen fugacity space. Additionally, increasing surface pressure suppresses the proposed filter index ratios, leading to conclusions similar to those of Hammond et al. (2025): the degeneracy between thick atmospheres and bare-rock planets remains unresolved, even when analyzing multiple photometric bands. Nevertheless, examining multiple bands provides more information than a single-band approach. A clear correlation with the redox state emerges across all photometric band indices, suggesting that, in principle, such observations could offer insights into the planet’s interior redox state. Seeing Figure 13, the ratio of any of those three photometric pairs could indicate a reduced or oxidized interior. For example, a low ratio value between F1500W with F1200W or F1000W, on average, is an indicator of a reduced interior while F1500W with F2100W is the opposite. However, for LP 791-18 d, this capability remains beyond the reach of current instrumentation.

Fig. 9 Airless body scenarios. Secondary eclipse depth spectra generated for various surface material scenarios without considering space weathering effects. Brightness temperature curves are overplotted as a comparative scale to highlight spectral features.

Fig. 10 Thermal emission spectrum for a simulation with surface pressure of 1 bar, redox state of fO2 − IW = 2, and ac = 10−6. Top: comparison between the line-by-line radiative transfer method at R = 60 000 and the correlated-k method at R = 1000, both convolved to R = 200. Bottom: comparison of convolved spectra from both line-by-line and correlated-k methods generated using petitRADTRANS, along with the output spectrum from our framework. All three approaches show negligible differences when predicting photometric fluxes in the JWST MIRI broad-band filters (see text for more details).

Fig. 11 Thermal emission spectra of 1bar and 100bar surface pressure simulations across various oxygen fugacities at ac = 10−6. The CO2 15 μm absorption band becomes prominent for oxidized atmospheres but it is subdued for thick atmospheres.

7.2 A JWST color-color diagram

A natural continuation of the analysis conducted so far, both in this manuscript and in the literature, is to introduce an additional photometric band and explore its potential. With three photometric bands, we can plot a color–color diagram and look for clustering of the possible scenarios. In Figure 14, we extract photometric secondary eclipse depths from our simulations for the characteristics of MIRI filters F1000W, F1500W, and F2100W, which correspond to 10 μm, 15 μm, and 21 μm, respectively. To estimate the observational uncertainties, we used the PANDEXO simulator (Batalha et al. 2017), adopting the same configuration proposed for LP 791-18 d observations by Benneke et al. (2024), namely five visits using the F1500W filter with the SUB256 sub-array. Identical integration times were assumed for the F1000W and F2100W bands to derive corresponding error bars. We observe clustering of the bare-rock scenarios in the upper-left corner of the figure, while the more reduced scenarios nearly overlap with the bare-rock color-color space. As the interior redox state approaches more oxidized conditions, the population gradually separates. The stable atmospheres are clearly separated from the bare-rock scenarios in the color-color diagram. Both surface pressure and graphite activity effects become negligible in the diagram, meaning that, for example, a separation of stable atmospheres with 1 bar is nearly indistinguishable from 100 bar or graphite activity values ranging from 10⁻⁷ to 10⁻³. To visualize the density of scenarios in color–color space, we compute a two-dimensional Kernel Density Estimation (KDE) using a Gaussian kernel with bandwidth determined via Scott’s Rule. The shaded contours in Fig. 14 correspond to the 2σ (95%) confidence region of the KDE distribution for each population. The two point clouds are for bare rocks and the stable atmospheres for the other.

7.3 Impact of potential photochemical hazes and clouds

Although our current modeling framework does not explicitly include photochemical haze formation, we acknowledge that such aerosols may play a non-negligible role in shaping the atmospheric observables of LP 791-18 d, particularly given the planet’s M-dwarf host star and potential for UV-driven chemistry. In analogy to Titan’s atmosphere, photochemical production of complex organic molecules could lead to high-altitude haze layers, especially in reducing or hydrocarbon-rich atmospheres. These hazes can increase the planetary albedo, reducing the dayside brightness temperature, and can also flatten or obscure infrared molecular features in emission spectra. To investigate the impact of photochemical hazes on the emission spectra of LP 791-18 d, we employ a two-step simplified parametric approach that captures their scattering and extinction. First, we scale the Rayleigh scattering cross sections of gas species using the haze_factor (H_f) parameter in petitRADTRANS (Mollière et al. 2019). This scaling enhances the native Rayleigh opacity of the gas mixture, simulating the wavelength-dependent scattering introduced by small-particle hazes. H_f multiplies the Rayleigh scattering opacity (in cm² g⁻¹) of gas species. That is, if the baseline Rayleigh opacity is denoted by κ_Rayleigh, the total scattering opacity becomes: $${\kappa }_{\text{total−}}=H_f\times {\kappa }_{\text{Rayleigh}}.$$(32)

This parameter allows a simplified treatment of enhanced atmospheric scattering due to photochemical hazes and scales with the assumed aerosol content as ~n_aerosols × σ_scat. For example, H_f = 1: standard Rayleigh scattering from gaseous species only. H_f > 1: mimics additional scattering by sub-micron aerosol particles. H_f ≫ 1 represents optically thick hazes with significant scattering at short wavelengths. This parameterization serves as a computationally efficient proxy for the enhanced scattering from small photochemical haze particles, such as complex hydrocarbons (e.g., tholins), which can form particles with radii ≲0.5 μm and introduce strong wavelength-dependent opacity in the UV to near-IR range. Additionally, we incorporate a toy model for wavelength-dependent absorption, representing Mielike opacity to mimic additional broadband absorption from larger aerosols, including condensate clouds, defined as: $${\kappa }_{\text{haze−}}(\lambda )={\kappa }_0{\left(\frac{\lambda }{1\mu \mathrm{m}}\right)}^{\beta }$$(33)

where κ₀ is the opacity at 1 μm which can vary from 0.1 to 10 cm² g⁻¹ depending on composition. β is the power-law index, typically β = −4 for small haze particles (Lecavelier Des Etangs et al. (2008)) and size α ≪ λ where Rayleigh scattering dominates. For larger particles, α ~ λ, for which Mie scattering dominates β factor tends to values closer to −1. Typical haze particles are in the range of 0.01–0.3 μm (Lavvas et al. (2008), Lavvas et al. (2010), Lavvas et al. (2019), Gao et al. (2021)). Nevertheless, larger particles can be produced from coagulation and growth (see, e.g., Figure 8 from Gao et al. 2021).

This absorption-like attenuation is applied to the synthetic spectra as an exponential suppression exp(−τ_haze) to explore the impact of Mie-type extinction. Together, these methods allow us to approximate a broad range of possible haze and condensate cloud effects in a computationally efficient way, without explicitly modeling microphysical processes.

In Figure 15a, we show the effect of H_f across four orders of magnitude, with solid lines. We have performed emission spectrum calculations across a range of haze enhancement factors, from H_f = 1 (clear atmosphere) up to 10⁴ (opaque hazy atmosphere), to assess their impact on the observable planet-to-star flux contrast (see Figure 15). These values encompass plausible ranges for hazy atmospheres, as observed in solar system bodies like Titan and in photochemically active exoplanets (Arney et al. 2016; Gao et al. 2021). The haze is assumed to be well-mixed and vertically distributed throughout the atmospheric pressure range. While this approach does not model haze microphysics explicitly, it provides a computationally efficient means to explore the sensitivity of thermal emission features to scattering and extinction by photochemical aerosols. Our simulations indicate that significant haze enhancement (H_f ≳ 10³) is required to appreciably mute emission features, particularly at near-infrared wavelengths. At H_f ~ 10⁴ the atmosphere is saturated and the short wavelengths (≲4 μm) are flattened.

However, their effect is negligible at higher wavelengths (>5 μm). This indicates that the important factor for our targeted spectral range is β. For this reason, we fixed H_f in an intermediate value (10³) and simulated various scenarios for β ranging from −4 to 0. Where in the limit of 0, the gray clouds are considered at various pressure positions (e.g., 1, 10, 100 mbar). The gray cloud extinction is wavelength-independent. These decks represent fully optically thick cloud layers that block emission from deeper layers. In this approach, the spectrum is truncated at the cloud-top level, simulating the limiting case of high condensate number density and large optical depth. This is shown in Fig. 16, where spectral features become increasingly muted as the cloud-top pressure decreases.

These treatments provide a unified and flexible framework to explore aerosol effects, bridging the continuum between photo-chemical hazes and condensate clouds, as well as from optically thin to thick conditions. However, the formation efficiency, composition, and optical properties of such hazes remain highly uncertain for closed-in tidally locked rocky exoplanets, particularly under the temperature regime of LP 791-18 d. A more detailed treatment is deferred to future studies incorporating coupled photochemistry and microphysics models.

In Figure 17 we reproduce again the spectral indexes shown in Figure 13 only for surface pressure 1 bar as a function of oxygen fugacity and a_c = 10⁻⁶ but now we introduce the effects of gray cloud decks at various vertical layers. Those effects can flatten the oxygen fugacity trend across the three JWST ratio indices we have experiment and produce degeneracies in the separation with a bare rock scenarios for the same indices. Finally, for each value of β and each gray cloud deck pressure considered, we simulated emission spectra across the entire parameter space, as shown in Figure 14, and calculated the average trajectory in the color–color space for each scenario affected by photochemical haze or clouds. We present these curves, which encompass the full range of possible effects from haze or clouds, in Figure 18. The shadowed regions are calculated in the same manner as Figure 14, With pink for bare rock clustering, gray the combined population of β −4 to −0.4 and green the combined gray cloud population at all pressure levels. Finally, with a red point we overplot the possible detection, including error bars, representing the expected observational uncertainty based on propagated errors from the JWST/MIRI secondary eclipse setup proposed by Benneke et al. (2024). This serves as a benchmark for evaluating the extent to which different scenarios are observationally distinguishable given current instrumental capabilities.

Fig. 12 Left: 15 μm band secondary eclipse depth as a function of fO2 across the surface pressure grid space with main graphite activity ac = 10−6. Dashed circles indicate stable atmospheres, as approximated using the methodology described in Section 5. Symbols in black indicate the stable atmospheres of all surface pressure scenarios for the main graphite activity value, while all other colors indicate the variability in graphite content for a 100 bar atmosphere. Right: 15 μm band for the bare-rock with various surface material scenarios, plotted on the same vertical scale as the left panel.

Fig. 13 Photometric band ratios for various surface pressures as a function of oxygen fugacity at ac = 10−6. We overplot the extreme values of graphite activity, ac, 10−7 and 10−3 for 100bar surface pressure. Top: ratio of the 12.8 μm to 15 μm spectral bands. Middle: ratio of the 15 μm to 21 μm. Bottom: ratio of the 10 μm to 15 μm. Variation in surface pressure is indicated with different colored shapes.

7.4 Sensitivity of observables to key system parameters

A central goal of this study was to understand how the physical and chemical properties of LP 791-18 d influence its observable spectral features, particularly in the context of thermal emission spectroscopy and the photometric capabilities of JWST. Among the various atmospheric characteristics, the most important drivers of emission signatures are the chemical composition and the thermal structure of the atmosphere. The chemical composition, encapsulated by the mean molecular weight, strongly affects atmospheric opacity and therefore the depth and shape of spectral features. In Figure 5, the dependence sensitivity of the mean molecular weight on three key interior-related parameters is shown: oxygen fugacity, graphite activity, and surface pressure. Our analysis shows that oxygen fugacity is the dominant factor, producing large compositional transitions between reduced and oxidized regimes. In contrast, graphite activity and surface pressure introduce subtler variations in molecular weight. However, atmospheric thickness, quantified by surface pressure, emerges as a crucial parameter for spectral observability. As illustrated in Figures 12 and 13, increasing surface pressure pushes the photosphere to higher altitudes, where temperatures are lower. This results in compressed and muted spectral features, making thick atmospheres appear more spectrally similar to bare-rock scenarios under certain conditions. Lastly, for atmospheric stability and escape, the most critical parameters are the mean molecular weight near the base of the wind and the stellar XUV flux, which together control whether the atmosphere can persist over geological timescales.

The thermal structure is shaped by several factors, but we find that surface albedo plays only a secondary role. As shown in Section 4.4, varying the surface albedo mainly affects the lower atmospheric layers, with negligible influence on the temperature profile at higher altitudes, where most of the thermal emission originates. Therefore, surface albedo does not significantly alter the spectral features observed in the mid-infrared.

Fig. 14 Proposed color-color diagram that breaks the degeneracy between bare-rock and thick oxidized stable atmospheres. The combination and names of the three employed JWST/MIRI photometric bands are shown on the axes. A set of bare-rock surface material scenarios is represented by colored stars (⋆). The various shapes indicate values across our grid of surface pressures, with colors representing the oxygen fugacity values. Blue represents reduced conditions and red represents oxidized conditions. Finally, we overplot different graphite activity (ac) values with different symbols (see legend: †, +, ×, and *) at random surface pressures. The shadowed areas represent the KDE probability density of the bare-rock and stable-atmosphere clusters.

8 Conclusions and discussion

There is a substantial body of literature on the theory of secondary atmospheres, curated in Earth and planetary sciences over decades (see, e.g., Gaillard & Scaillet 2014). However, we need more empirical evidence on the early phases of secondary atmosphere generation, along with a diverse sample of the star-planet population. Our current observing capabilities limit our ability to fully characterize the abiotic baseline of a planet, which is crucial for assessing the potential for false-positive biosignature signals.

An important goal, as community, is to identify selected candidates, such as LP 791-18 d, that are in a phase similar to early-Earth, where plain outgassing, a process responsible for secondary atmospheres, is amenable to observations from contemporary instrumentation in the upper layers of the atmosphere, where physical processes are out of chemical equilibrium and cooler. We can make a connection to the interior redox state from those atmospheric observations. In this work, we have shown that the bulk mean molecular weight (μ) of the atmosphere is controlled by oxygen fugacity rather than bulk planetary metallicity. Variations in redox conditions determine the dominant molecular species and their partitioning between oxidized and reduced forms, which in turn shape the overall atmospheric composition and μ. We have demonstrated that a variety of atmospheres can be generated from different interior redox states, assuming that the surface pressure of LP 791-18 d varies from 0.1 to 100 bars. More simulations are needed for similar candidates, in combination with observations, to help disentangle the value of oxygen fugacity. A few promising candidates have been recently identified. For example, Bello-Arufe et al. (2025) present evidence for volcanic activity on the sub-Earth planet L98-59b. A theoretical analysis of a larger sample of similar objects, combined with future observations of atmospheric characterization, could provide insightful information on the planet’s chemical inventory relative to stellar abundances. This is a crucial connection when estimating the abiotic baseline of a rocky planet and comparing it with its observable signal to detect potential deviations that might be associated with extraterrestrial life.

Fig. 15 Effects of photochemical haze aerosols. (a) Reduced atmospheric scenario with ac = 10−6 and surface pressure 1 bar is simulated across the haze enhancement factor Hf and β. (b) Same atmosphere as (a), but the full spectral range is presented and across β variances. (c) Oxidized atmosphere with ac = 10−6 and surface pressure 1 bar is presented to show the effects on the 15 μm CO2 absorption band.

8.1 The bare rock – thick atmosphere degeneracy

The efforts to characterize small rocky exoplanets have been initiated through JWST observations targeting a handful of such worlds. The observation strategies have primarily focused on secondary eclipse photometry in the 15 μm band, where CO₂ shows strong absorption. The rationale behind this observational approach is that comparing the 15 μm band with models of bare-rock and CO₂-rich atmosphere scenarios was expected to reveal a distinction in secondary eclipse depths (see, e.g., Zieba et al. 2023; Xue et al. 2024). In the recent work by Ducrot et al. (2024), they performed observations with one additional photometric MIRI-JWST filter, namely F1280W in the 12.8 μm band, and concluded that a likely bare-rock scenario exists for TRAPPIST-1 b. Hammond et al. (2025) argued that photometric observations in the 15 μm band alone are not sufficient to unambiguously distinguish between a bare-rock planet and an atmosphere, even with the combination of the two MIRI filters F1280W and F1500W. This is because thick atmospheres can still produce thermal emission spectra that resemble nearly featureless blackbodies at these wavelengths. In particular, they significantly suppress the 15 μm CO₂ absorption band. As a result, distinguishing such atmospheres from bare-rock emission becomes challenging with current observational capabilities. Instead, they proposed using F1500W phase curve observations to detect a possible heat redistribution signal to the night side.

In this work, we chose LP 791-18 d as a case study planet and performed comprehensive climate simulations, which were computationally feasible in 1D modeling. These simulations include outgassing at the lower boundary, thermochemistry, and photo-chemistry coupled with a radiative-convective model to solve for the thermal structure. This was possible because, based on tidal heat flux calculations for LP 791-18 d, we could constrain its interior melt temperature, which is a critical parameter for our outgassing model. The output of secondary atmospheres was tested for stability, utilizing the criterion of hydrodynamic escape, which was modeled using the isothermal wind approximation. The isothermal approximation was used as a diagnostic tool to explore atmospheric stability across the parameter space defined by surface pressure, oxygen fugacity, and graphite activity. Including radiative cooling through atomic and molecular line emission would alter the thermospheric temperature and pressure structure, reducing the efficiency of hydrodynamic escape. Then, the mass-loss rates presented here should be regarded as qualitative upper bounds. This underscores the value of future models incorporating a more detailed treatment of upper atmospheric physics. A more accurate assessment of atmospheric stability could be achieved by coupling our 1D lower atmosphere model with an upper atmosphere module that solves the full energy and momentum equations (see, e.g., Yelle 2004; Tian et al. 2005; Johnstone et al. 2018). Nonetheless, the energy-limited approximation employed here already accounts for unresolved processes via the efficiency factor η, which encapsulates the complexity of upper-atmospheric energy redistribution in a time-independent and agnostic manner (Owen 2019). To avoid overcomplicating the numerical framework, we adopted a simplified approach that still captures the peak temperature of the thermal profile, located near the base of the wind region (Lampón et al. 2020). If the conditions for wind formation are not met, the atmosphere is considered stable. In such cases, Jeans and non-thermal escape processes can shape the chemical inventory of the atmosphere over geological timescales, but they are unlikely to significantly reduce the total atmospheric mass. The sensitivity of atmospheric stability to the planet’s interior redox state is illustrated in several figures. In the bottom panel of Figure 8, we show the mean molecular weight at R_xuv, where most of the XUV energy is deposited, and overplot the stable thermosphere solutions. Stability in these models arises from two main factors: the intensity of incoming XUV radiation and the mean molecular weight of the background atmosphere at the wind base. In our case, the XUV energy is constant since we study the current scenarios in time. Atmospheric stability is generally achieved for scenarios with an oxidized redox state in the interior of LP 791-18 d. This stability occurs for oxygen fugacity values higher than ≈3, and for surface pressures approaching 100 bar, stability is achieved at ≈2. For varying graphite content, stability can be achieved starting from oxygen fugacity ≈2 for all surface pressures. Our solutions show smaller, stable regions, as seen in Figure 8. Although our results highlight these stable regions (as shown in Figure 8), higher-resolution simulations in the redox–pressure parameter space are needed to better delineate their boundaries.

For all the simulations that achieved steady-state solutions, we generated synthetic spectra for the secondary eclipse depth (e.g., Figure 11 for 1 bar mock spectra as a function of interior redox state). In addition, we considered a variety of bare-rock scenarios to cover cases ranging from featureless spectra (e.g., a metal-rich surface) to cases with strong features (e.g., ultramafic). We then examined the effectiveness of various observational strategies. Specifically, for LP 791-18 d, the F1500W filter is the only one utilized in JWST observations (Benneke et al. 2024). In Figure 12, we present the secondary eclipse depths as a function of the interior redox state and surface material composition for the 15 μm band. In the same figure, we indicate which generated secondary atmospheres can maintain a stable thermosphere, considering the planet’s space environment. Given the estimated age of the LP 791-18 system (~500 Myr) and our calculated atmospheric mass-loss rates for unstable scenarios (~10⁵−10⁸ kg s⁻¹), we find that continuous volcanic outgassing alone is likely insufficient to sustain the atmosphere over geological timescales. To assess this, we estimate the total cumulative atmospheric loss over the planet’s lifetime as follows: $$\begin{array}{rl}\mathrm{\Delta }M=\dot{M}\times t\approx & \left({10}^5\text{−to−}{10}^8\mathrm{k}\mathrm{g}{\mathrm{s}}^{-1}\right)\\ & \times (500\times {10}^6\mathrm{y}\mathrm{r})\approx (1.6\times {10}^{21}\text{−to−}1.6\times {10}^{24})\mathrm{k}\mathrm{g}.\end{array}$$(34)

At the lower end of this range, the total mass lost approaches, or even exceeds, typical estimates for the initial volatile inventory of rocky planets. The initial volatile budget of terrestrial planets is poorly constrained but is generally estimated to be in the range of 10⁻⁵ to 10⁻² of the planetary mass (Elkins-Tanton & Seager 2008; Zahnle et al. 2007; Schaefer & Fegley 2010). For Earth, this corresponds to a bulk silicate inventory of approximately (1.8–3.0) × 10²² kg of volatiles, primarily in the form of H₂O, CO₂, and other light elements (Hirschmann 2006; Dasgupta & Hirschmann 2010). For LP 791-18 d, a conservative volatile reservoir on the order of 10⁻³ M_p would yield ~6 × 10²¹ kg of volatiles. This estimate suggests that atmospheric escape at mass-loss rates of 10⁷−10⁸ kg s⁻¹ would exhaust such a reservoir in less than ~200 Myr. Even for lower escape rates (10⁵−10⁶ kg s⁻¹), the volatile budget would be significantly eroded within 1–2 Gyr. Furthermore, the XUV activity of the host star is expected to be substantially higher in its early evolutionary stages. Therefore, the long-term persistence of an atmosphere on LP 791-18 d requires either (1) formation with an unusually large volatile inventory (e.g., from a carbon-rich or water-rich disk), or (2) a redox state that favors the buildup of a dense, high-mean-molecular-weight atmosphere capable of resisting hydrodynamic escape. Under these assumptions, we conclude that reduced interior scenarios (low f_O2) are unlikely to retain secondary atmospheres over gigayear timescales, especially if the outgassed composition is dominated by light molecules. This reinforces our finding that only oxidized atmospheric compositions (with f_O2 − IW ≳ 2) produce atmospheres that are stable against hydrodynamic loss.

Nevertheless, our stability criterion is based on approximations. While we cannot determine with precision the exact redox threshold below which atmospheres undergo continuous hydrodynamic escape throughout their lifetimes, our results suggest that reduced atmospheres with f_O2 − IW ≲ −1 are likely to be dominated by very light species, both in total atmospheric mass and at R_xuv. In the absence of other mitigating mechanisms, such atmospheres are unlikely to suppress the planetary wind regime, except through a significant decline in the stellar XUV flux over time. For redox states between 0 ≲ fO₂ − IW ≲ 2, even though our current framework predicts planetary wind regimes, the atmospheric composition is rich in heavier elements like carbon and oxygen. Thus, due to molecular diffusion in the upper atmosphere, which preferentially allows light elements to escape, we expect the bulk atmospheric composition to become progressively enriched in heavier species over geological timescales. This could potentially stop the hydrodynamic escape flow until degassing from the lower boundary brings the atmospheric composition to a state where escape is initiated once more. In the absence of hydrodynamic escape, the remaining mass loss mechanisms, Jeans escape and non-thermal processes, can still cause significant compositional evolution across Gyr timescales. These processes preferentially remove lighter species (e.g., H₂, He, CH₄), enriching the upper atmosphere in heavier molecules and atoms (e.g., CO₂, O₂, N₂). This selective loss increases the mean molecular weight near the exobase, which in turn reduces escape efficiency and may help stabilize the atmosphere over time (Zahnle & Catling 2017). As a result, the present-day atmospheric composition may differ from the initially outgassed one, especially in intermediate redox states scenarios, where both light and heavy species are outgassed in the atmosphere. This evolutionary effect should be considered when interpreting spectral features and mean molecular weight as indirect diagnostics of the planet’s interior redox state.

For all the above mock spectra, we utilized PANDEXO (Batalha et al. 2017) to estimate the observational uncertainties, given the characteristics described in Benneke et al. (2024). The uncertainties are estimated to be on the order of 30-45 ppm, depending on the scenario. From Figure 12, we can expect to distinguish an atmospheric detection if the atmosphere is highly oxidized and thin or highly oxidized with low graphite content a_c ≲ 10⁻⁴. Nevertheless, it would be difficult to make a bold claim even for a thin atmosphere (≈0.1 bar) given the estimated uncertainties. Following the work of Ducrot et al. (2024), we investigated the combination of two photometric measurements from the available MIRI filter list. After performing numerical experiments, we show the best combinations in Figure 13. We display the ratios of the filters as a function of the interior redox state. For the combination of filters F1280W-F1500W, thin atmospheres show differences of up to 25% in secondary eclipse depth, which is at the limit of detection, considering the uncertainties and the thermal contrast estimated for LP 791-18 d. For thick atmospheres, the ratio approaches unity, confirming the degeneracy observed in Hammond et al. (2025). If we consider alternative pairs of filters, F2100W-F1500W and F1000W-F1500W, we obtain similar results to the first pair. However, we can achieve better contrast in the filter fraction and can more effectively separate at least the thin atmospheres. For thick atmospheres, the signal is suppressed. Additionally, degeneracies can be produced even for thin atmospheres if we add effects of thick photochemical haze or gray cloud deck at low pressure levels as it is shown in Figure 17. Our simulations show that we cannot unambiguously distinguish between an atmosphere and bare rock with observations of one or even two photometric bands.

Finally, we propose a novel approach, presented in Figure 14, which utilizes three MIRI photometric filter observations and plots the estimated results as a color–color diagram. In this diagram, all the bare-rock scenarios cluster in the upper-left corner, while the stable, highly oxidized atmospheres are located in the bottom-right corner. Surprisingly, when comparing these three photometric bands, the dependence on surface pressure and graphite activity is eliminated, since most of the plotted points follow a smooth slope with small dispersion. Due to the non-linearity of the redox state and the strength of individual spectral bands, the intermediate region of the color-color points is mixed with both reduced and oxidized scenarios, while the extremes are well-separated, as expected. Similarly, the clustering of reduced atmospheres is observed, as they do not exhibit strong spectral features in contrast to oxidized atmospheres, as shown in Figure 11. A real observation would likely result in a measurement that lies in one of the two shadowed areas. With the estimated uncertainties of ±30–45 ppm, it is expected to be able to separate an atmosphere from an airless case for LP 791-18 d, regardless of its surface pressure or graphite content. We performed a parametric study of photochemical hazes and both wavelength-dependent semi-gray and gray clouds across all simulations presented in this work and examined their effects on our proposed color–color diagram, shown in Figure 18. Our results show that, while variations in the properties of photochemical hazes, semi-gray, and gray clouds can modify the trends followed by different atmospheric scenarios in the diagram, the region corresponding to highly oxidized, stable atmospheres remains clearly distinguishable from that of bare rock surfaces. Based on the uncertainty estimation and its propagation for constructing the color-color diagram, we overplotted a representative synthetic detection, as shown in Figure 18. In order for the proposed color–color system to be observationally useful in the JWST era, the typical uncertainty required is on the order of ~50 ppm, which is currently achievable only for highly oxidized atmospheres. Finally, we emphasize that the proposed framework provides a foundation for prioritizing future observations, particularly as instrumental capabilities and retrieval techniques advance.

Fig. 16 Gray cloud deck effects on the emission spectra. We color-coded the various pressure levels of the position of the gray cloud layer (see legend in inset).

Fig. 17 Gray cloud effects in the spectral indexes introduced in Figure 13. With red squares we show the clear atmospheric scenario and with colored lines the gray cloud deck at various pressure levels. With the shadowed gray area, we show the 1σ variation of all the bare rock scenarios simulated in this work. The example shown is the 1 bar surface pressure at various oxidation states for ac = 10−6. With this figure we emphasize that gray clouds introduce degeneracies between bare rocks and atmospheres even for highly oxidized atmospheres.

Fig. 18 Proposed color-color diagram including the effects of cloud and hazes. The black solid line gives the mean value of the atmospheric scenarios shown in Figure 14 and the colored stars and pink shaded area are the same as in Figure 14. With dotted colored lines we show the effects of β factor from −4 to −0.4 and with gray shadowed region their combined KDE for the stable atmospheres. With dashed colored curves the limit of β → 0 is shown which represents the effect of gray clouds at various pressure layers and the green shaded area is the KDE of the combined gray cloud stable atmospheres scenarios. For reference, a red point with error bars is overplotted, representing the expected observational uncertainty based on propagated errors from the JWST/MIRI secondary eclipse setup proposed by Benneke et al. (2024). This serves as a benchmark for evaluating the extent to which different scenarios are observationally distinguishable given current instrumental capabilities.

8.2 The importance of stable thermospheres in habitability and detectability of rocky exoplanets

As we discussed above, stable thermospheres retain the atmospheric mass against bulk escape and play a crucial role in the habitability of rocky planets. This stability primarily results from the combination of two physical parameters: stellar incoming radiation and the mean molecular weight. While these parameters help stabilize the atmosphere, they pose challenges for detectability. Due to cooler atmospheres with high mean molecular weight, detecting such atmospheres becomes extremely challenging, even for future-generation detectors. In transmission, as mentioned above, the spectral feature is proportional to the atmospheric scale height parameter, given by $\delta D\approx \frac{2R_{p}H}{R_{\ast }^2}$ where H depends on the temperature and R_p is the radius where the optical depth of the spectral wavelength is ≈1. It is thus evident that alternative approaches must be considered that can enhance the detectability of small rocky planets.

When comparing the photoabsorption cross-sections of various species across different spectral regions, UV absorption (particularly in the FUV–NUV range) is typically many orders of magnitude stronger than in the optical to near-infrared (see, e.g., Heays et al. 2017; Venot et al. 2018; Malik et al. 2019). The higher-energy photons of the stellar spectrum are absorbed in the upper atmosphere (e.g., from the thermosphere to the exosphere), such as the R_xuv concept used in Section 5. Therefore, it is evident that transit radius in the FUV-NUV-XUV wavelengths are optimizing our probabilities for characterization of small rocky exoplanets with transmission spectroscopy.

Furthermore, stable upper atmospheres have one more characteristic that makes them appealing targets for future observations: their thermal structure is easier to model than the lower atmosphere. This is because the dominant factors in the energy balance equation are the incoming UV energy as the heating mechanism and thermal conduction as the cooling mechanism (Bougher & Roble 1991; Gkouvelis et al. 2021). The thermospheric thermal structure has the morphology of an inversion profile, which leads to an isothermal layer at the exosphere (Banks & Kockarts 1973a; Brecht et al. 2022). Typical temperatures in this region range from a few hundred to a few thousand kelvins, which further adds to their observability. These components benefit future UV transmission spectroscopic observations when coupled with Bayesian inference frameworks to characterize small rocky exoplanets. While no detector is nowadays capable of performing these observations, an effort is ongoing to include a UV spectropolarimeter on the future Habitable Worlds Observatory mission (see, Muslimov et al. 2024).

Acknowledgements

L.G. acknowledges the “Severo Ochoa excellence” visitor program for partially funding this study. Author F.J.P. acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S MICIU/AEI/10.13039/501100011033 and Ministerio de Ciencia e Innovación through the project PID2022-137241NB-C43.

References

All Tables

All Figures

	Fig. 1 Mind map describing the different parts of our methodology. The simulations are composed of three core pieces coupled together. Radiative convective model, chemical kinetics, and an equilibrium inter-phase multi-phase chemistry between the surface and the atmosphere. Post-numerical simulation analysis is done to produce transmission and emission mock spectra.
In the text

Fig. 2 Modeled internal energy balance for LP 791-18 d. Plotted are convective (blue curves) and tidal (black and red curves) heat fluxes as a function of the mantle’s potential temperature. The isothermal mantle solidus is set at Tp = 1600 K. We follow the prescription of Peterson et al. (2023) and use the same parameters therein to compute: i) the convective flux assuming soft turbulence in the mantle (e.g., Solomatov 2015), whereby the four sets of curves cover a range of melt fraction coefficients (B = 10, 20, 30, and 40); and ii) the tidal heating flux of a Maxwell viscoelastic solid mantle. In addition, we compute the tidal heat flux allowing for the fluid contribution of melt in the mantle, following Farhat et al. (2025), using a conservative estimate of the dissipative timescale, σR = 10−3 s−1. Interior thermal equilibria, corresponding to the balance between tidal heating and convective cooling, are shown by the colored markers.

In the text

	Fig. 3 Converged simulations for a surface pressure of 1 bar. Volume mixing ratios (VMR) of important greenhouse gases that significantly contribute to emission and transmission spectroscopy are shown. Three representative oxygen fugacity cases are presented to illustrate the transition from a reduced to an oxidized atmosphere.
In the text

	Fig. 4 Converged simulations for a surface pressure of 1 bar. The thermal structure is shown for all simulated scenarios within the 1 bar surface pressure grid. We note the thermal inversion in the upper layers of the simulated atmosphere, caused by infrared absorption due to high water vapor abundances.
In the text

Fig. 5 Mean molecular weight (μ) of the outgassed atmosphere (in amu) as a function of varying surface pressure, oxygen fugacity (fO2), and graphite activity (aC). Each panel isolates the effect of one parameter by holding it fixed while allowing the other two to vary, as indicated in each panel. The value of μ shown corresponds to mean molecular weight of the volatile partition at a given degassing rate. This figure illustrates how redox state, pressure, and carbon content influence the dominant molecular species and therefore control the overall atmospheric composition.

In the text

	Fig. 6 Top: vertical temperature profile for a reduced and oxidized atmosphere at surface pressure of 1 bar and ac = 10−6. We show simulations for surface albedo ranging from 0.01 to 0.6. Bottom: detail from the top plot, focusing on the layers close to the surface.
In the text

Fig. 7 Wind velocity profile as a function of radius for a surface pressure of 1 bar and fO2 − IW = 0. The blue solid line represents the wind velocity profile as a function of radius, the green dashed line indicates the exobase estimation, and the orange dashed line shows the escape velocity, which decreases as ~ r−2. The red dashed line marks the critical radius where the sonic velocity of the flow is reached, while the black solid line represents the Hill radius, beyond which material is no longer gravitationally bound to the planet. In this example, a planetary wind regime is established since rs ≤ Rexo.

In the text

Fig. 8 Mean molecular weight across the simulation parameter grid space. Top: Total atmospheric mass mean molecular weight (in amu), calculated once the simulations reach steady state. Bottom: Mean molecular weight at Rxuv, the region where the energetic portion of the stellar spectrum is deposited, forming the base of the wind. The thick dashed black line separates the parameter space between the wind regime and stable thermospheres.

In the text

	Fig. 9 Airless body scenarios. Secondary eclipse depth spectra generated for various surface material scenarios without considering space weathering effects. Brightness temperature curves are overplotted as a comparative scale to highlight spectral features.
In the text

Fig. 10 Thermal emission spectrum for a simulation with surface pressure of 1 bar, redox state of fO2 − IW = 2, and ac = 10−6. Top: comparison between the line-by-line radiative transfer method at R = 60 000 and the correlated-k method at R = 1000, both convolved to R = 200. Bottom: comparison of convolved spectra from both line-by-line and correlated-k methods generated using petitRADTRANS, along with the output spectrum from our framework. All three approaches show negligible differences when predicting photometric fluxes in the JWST MIRI broad-band filters (see text for more details).

In the text

	Fig. 11 Thermal emission spectra of 1bar and 100bar surface pressure simulations across various oxygen fugacities at ac = 10−6. The CO2 15 μm absorption band becomes prominent for oxidized atmospheres but it is subdued for thick atmospheres.
In the text

Fig. 12 Left: 15 μm band secondary eclipse depth as a function of fO2 across the surface pressure grid space with main graphite activity ac = 10−6. Dashed circles indicate stable atmospheres, as approximated using the methodology described in Section 5. Symbols in black indicate the stable atmospheres of all surface pressure scenarios for the main graphite activity value, while all other colors indicate the variability in graphite content for a 100 bar atmosphere. Right: 15 μm band for the bare-rock with various surface material scenarios, plotted on the same vertical scale as the left panel.

In the text

Fig. 13 Photometric band ratios for various surface pressures as a function of oxygen fugacity at ac = 10−6. We overplot the extreme values of graphite activity, ac, 10−7 and 10−3 for 100bar surface pressure. Top: ratio of the 12.8 μm to 15 μm spectral bands. Middle: ratio of the 15 μm to 21 μm. Bottom: ratio of the 10 μm to 15 μm. Variation in surface pressure is indicated with different colored shapes.

In the text

Fig. 14 Proposed color-color diagram that breaks the degeneracy between bare-rock and thick oxidized stable atmospheres. The combination and names of the three employed JWST/MIRI photometric bands are shown on the axes. A set of bare-rock surface material scenarios is represented by colored stars (⋆). The various shapes indicate values across our grid of surface pressures, with colors representing the oxygen fugacity values. Blue represents reduced conditions and red represents oxidized conditions. Finally, we overplot different graphite activity (ac) values with different symbols (see legend: †, +, ×, and *) at random surface pressures. The shadowed areas represent the KDE probability density of the bare-rock and stable-atmosphere clusters.

In the text

Fig. 15 Effects of photochemical haze aerosols. (a) Reduced atmospheric scenario with ac = 10−6 and surface pressure 1 bar is simulated across the haze enhancement factor Hf and β. (b) Same atmosphere as (a), but the full spectral range is presented and across β variances. (c) Oxidized atmosphere with ac = 10−6 and surface pressure 1 bar is presented to show the effects on the 15 μm CO2 absorption band.

In the text

	Fig. 16 Gray cloud deck effects on the emission spectra. We color-coded the various pressure levels of the position of the gray cloud layer (see legend in inset).
In the text

Fig. 17 Gray cloud effects in the spectral indexes introduced in Figure 13. With red squares we show the clear atmospheric scenario and with colored lines the gray cloud deck at various pressure levels. With the shadowed gray area, we show the 1σ variation of all the bare rock scenarios simulated in this work. The example shown is the 1 bar surface pressure at various oxidation states for ac = 10−6. With this figure we emphasize that gray clouds introduce degeneracies between bare rocks and atmospheres even for highly oxidized atmospheres.

In the text

Fig. 18 Proposed color-color diagram including the effects of cloud and hazes. The black solid line gives the mean value of the atmospheric scenarios shown in Figure 14 and the colored stars and pink shaded area are the same as in Figure 14. With dotted colored lines we show the effects of β factor from −4 to −0.4 and with gray shadowed region their combined KDE for the stable atmospheres. With dashed colored curves the limit of β → 0 is shown which represents the effect of gray clouds at various pressure layers and the green shaded area is the KDE of the combined gray cloud stable atmospheres scenarios. For reference, a red point with error bars is overplotted, representing the expected observational uncertainty based on propagated errors from the JWST/MIRI secondary eclipse setup proposed by Benneke et al. (2024). This serves as a benchmark for evaluating the extent to which different scenarios are observationally distinguishable given current instrumental capabilities.

In the text