© The Authors 2025

1 Introduction

Among the thousands of exoplanets discovered so far, only a few tens of them are temperate, with a size and/or a mass compatible with a rocky composition. Their small size makes them difficult to detect and even more difficult to characterize. Hu et al. (2024) report the potential detection of a secondary atmosphere on the hot planet 55 Cancri e. A sulfur-rich atmosphere may have been detected on L 98-59 d (Banerjee et al. 2024; Gressier et al. 2024); however, no definitive detection of a secondary atmosphere has been made to date on a temperate Earth-sized rocky exoplanet. The James Webb Space Telescope (JWST) has the requisite sensitivity to probe the atmospheres of a few specific temperate rocky planets, mostly through transit spectroscopy and occultations (i.e., secondary eclipses). Those specific targets are exoplanets orbiting the coolest and smallest close-by stars, which are called red dwarf stars. Red dwarf stars are favorable targets (Triaud 2021; Nutzman & Charbonneau 2008; Sullivan et al. 2015; Barclay et al. 2018), because (1) their small size makes it easier to detect and characterize small transiting planets, (2) their lower luminosities result in more frequent planetary transits and occultations for the same stellar irradiation, (3) they are the most common objects in the Galaxy, and (4) temperate planets around them are more likely to be short-period rocky planets, allowing us to accumulate many transits and occultations or to perform phase curve observations.

With its seven transiting temperate rocky planets around the coolest host star known to date, the TRAPPIST-1 system is the most studied system of this kind (Turbet et al. 2020a). Thanks to very precise measurements of radii (Ducrot et al. 2020) and masses (Agol et al. 2021), TRAPPIST-1 planets are the most well-known rocky planets outside our Solar System. Since its discovery in 2016–2017 (Gillon et al. 2016, 2017), a large number of studies has assessed the possibility of atmospheres on these planets (Turbet et al. 2020a), and their detectability with the JWST (Lustig-Yaeger et al. 2019; Fauchez et al. 2019). The host star, with its high flux in the extreme ultraviolet and X-rays (XUV) and its long and hot pre-main sequence phase, is expected to result in significant atmospheric escape (Bolmont et al. 2017; Bourrier et al. 2017). However, depending on the formation and evolution pathways of the planets, they could have retained a part of their primary atmosphere or formed a secondary atmosphere (see Moore & Cowan (2020); Kite & Barnett (2020) for the evolution of rocky planets and potential secondary atmospheres, Coleman et al. (2019); Turbet et al. (2020a) for the particular case of TRAPPIST-1).

Considering their transiting nature combined with the infrared brightness (K = 10.3) and the Jupiter-like size of their host star (∼12% the radius of Sun), they have been designated as the most promising candidates for the first detailed study of rocky Earth-sized exoplanets with JWST (Gillon et al. 2020). Since its launch, all the TRAPPIST-1 planets have been observed in transit (JWST programs GTO 1201, GTO 1331, GO 1981, GO 2420, GO 2589, GO 6456). Unfortunately, the activity of the star can mimic the signatures from the planet’s atmosphere in transit (Lim et al. 2023). The host star’s surface presents inhomogeneities that can cause artificial spectral features in the transit, due to unocculted spots on the star rather than an atmosphere on the planet (Rackham et al. 2018). One way out of this problem is to observe the planets in emission rather than in transmission spectroscopy. In this case the surface inhomogeneities of the star have a much weaker impact, although not inexistent (Fauchez et al. 2025). The planets in the TRAPPIST-1 system that are hot enough to be observed with this method by the JWST are TRAPPIST-1 b and c. The two planets’ emissions have been observed through their secondary eclipses in the 15 µm MIRI filter (JWST programs GTO 1177 and GO 2304) and in the 12.8 µm filter (GTO 1279 and GO 5191).

The 15 µm eclipse observations of TRAPPIST-1 b already ruled out – at first sight – a thick CO₂ atmosphere (Greene et al. 2023), but they still left space for thinner atmospheres with different compositions (Ih et al. 2023). Taking into account the second photometric point, at 12.8 µm, Ducrot et al. (2025) pointed out that secondary eclipses can also theoretically be compatible with a CO₂-dominated atmosphere rich in hazes. A way to break this degeneracy inherent to dayside emission is to observe phase curves in order to get the nightside thermal emission of the planet (Hammond et al. 2025).

More recently, a phase curve of the system has been observed and carefully selected to be centered on the eclipses of TRAPPIST-1b and c (GO 3077). In this work, we present 1D radiative-convective and 3D numerical climate simulations of various types of atmosphere of TRAPPIST-1 b and compare them to all the available thermal emission observations (eclipses and phase curve). Although this work is an in-depth study of a single planet (the only moderately irradiated rocky planet for which we have such a large amount of observation data), its approach and results are meant to be broadly applicable to other temperate rocky exoplanets for which we have emission observations. This includes TRAPPIST-1c, but also LHS-1140 c, LHS-1478 b, GJ-3473 b, GJ-357 b, HD-260655 b, L-98-59 c, LTT-3780 b, TOI 1468 b, and TOI-270 b (the nine planets of the Hot Rocks Survey, GO 3730); LP 791-18 d (GO 6457); and Gl 486 b (GO 1743). Moreover, a 500-hour Director’s Discretionary Time program on JWST (Rocky Worlds DDT) will provide more mid-infrared observations of secondary eclipses of small rocky exoplanets, in order to probe their ability to maintain an atmosphere (Fig. 1).

Fig. 1 Exoplanets under consideration for the Rocky Worlds DDT. Planets above the cosmic shoreline are expected to have lost their atmosphere. The circled exoplanets have already been observed in emission during JWST’s first three cycles. The exoplanets for which the observations concluded there is no atmosphere are in red. We also added the Solar System terrestrial planets. The XUV irradiation is computed from Zahnle & Catling (2017).

2 Methods

Starting from the secondary eclipse observations at 15 and 12.8 µm (Greene et al. 2023; Ducrot et al. 2025), we attempted to identify all the possible families of atmosphere of TRAPPIST-1 b that could account for the high dayside emission flux measured in secondary eclipses. For this, we used a suite of numerical climate models, ranging from (1) 1D radiative-convective simulations designed for a fast exploration of the parameter space to (2) 3D global climate calculations to more reliably simulate the thermal structure and temperature contrasts across the planet, which play a major role in the planet’s emission signal. These simulations were then post-processed with radiative transfer tools to generate synthetic eclipse spectra and thermal phase curves, for comparison to the observations.

In practice, we first carried out dozens of numerical climate simulations with the 1D model exo_k (see Section 2.2) for a wide range of atmospheric compositions, and computed secondary eclipse spectra (see Section 2.3) for a large range of plausible atmospheric redistribution factors. For atmospheric compositions that yield good match with eclipse data, we used a 3D climate model to produce much more realistic simulations (see Section 2.1) with self-consistent heat redistribution, verified or not the good match with the eclipse data, and calculated thermal phase curves in the 15 µm MIRI filter. The list of all the 3D simulations used in this work is given in Table 1.

The models and tools used in our study are detailed in the subsections below.

2.1 3D simulations with the generic PCM

For the 3D atmospheric simulations of TRAPPIST-1 b, we used the Generic Planetary Climate Model (Generic PCM), a 3D numerical climate model which has now been extensively applied to a large diversity of exoplanets (Wordsworth et al. 2011; Charnay et al. 2015; Bolmont et al. 2016; Turbet 2018; Charnay et al. 2021; Turbet et al. 2023; Teinturier et al. 2024).

2.1.1 Model parameters

For the stellar, planetary and orbital parameters, we used the results of Agol et al. (2021). We assumed (unless otherwise stated) a circular, 0° obliquity, tidally locked orbit. Table 2 summarizes all the parameters used to set up our simulations. Sensitivity tests for some of these parameters are detailed in Section 3.

2.1.2 Model resolutions

For all our simulations, we used a standard horizontal resolution of 72 × 46 in longitude × latitude. For the vertical resolution, our standard simulations use 40 pressure layers, but for some simulations we had to change the number and distribution of pressure layers (suppressing layers for the simulation with atmospheric collapse; adding and refining the upper layers grid for the simulation with a thermal inversion). Given that we simulated a broad range of atmospheric compositions in the Generic PCM, with various radiative and dynamical timescales involved, we adapted the time step (ranging from ~0.5 s to ~60 s) for all the simulations.

2.1.3 Radiative transfer

Regarding the radiative transfer calculations, the model is based on the correlated-k approach (Fu & Liou 1992). For the correlated-k tables, we use 31 spectral bands in the visible-near infrared (0.3–6.6 µm), corresponding to the stellar incoming radiation, and 40 in the thermal infrared (1–50 µm), corresponding to the planet emission wavelengths. This corresponds to a spectral resolution of R=10 for both regimes. We used 16 gauss points for the cumulated distribution function (g) integration, 20 for the N₂-CO₂ atmospheres.

The correlated-k tables used are summarized in Table 3. The correlated-k tables of N₂+CH₄+C₂H₄ and N₂+CH₄+CO₂ have been calculated by combining single gas opacities from Exo-mol (Hill et al. 2013; Mant et al. 2018; Hill et al. 2013; Chubb et al. 2021; Yurchenko et al. 2020) using exo_k. For N₂+CO₂ mixtures, original correlated-k tables have been calculated from line-by-line calculations of the high-resolution absorption spectra of CO₂ broaden by N₂ (Chaverot et al., in prep.). Absorptions have been computed following a Voigt profile with empirical corrections in the far wings (e.g., Tran et al. 2018) based on laboratory measurements (Burch et al. 1969; Perrin & Hartmann 1989; Tran et al. 2011). For all atmospheric compositions, we also included the effect of collision-induced absorptions (CIAs) of CO₂–CO₂ (Gruszka & Borysow 1997; Tran et al. 2024) N₂–N₂ (Karman et al. 2019), CH₄–CH₄ (Karman et al. 2019), CO₂–CH₄ (Turbet et al. 2020b) and CH₄–N₂ (Karman et al. 2019).

2.1.4 Convection and cloud formation

The subgrid-scale dynamical processes (dry turbulent mixing and convection) are taken into account through a convective adjustment scheme that adjusts unstable temperature profile to the dry adiabat (Forget et al. 1999). For the study of the collapse of the CO₂ atmospheres, we activated the CO₂ condensation scheme (Forget et al. 2013). This scheme is performed when a grid cell reaches 100% of saturation. CO₂ when condensing can form CO₂-ice clouds and/or surface frost. Cloud ice particles can sediment and accumulate at the surface, and the surface ice layer can itself re-sublimate in the atmosphere, depending on the thermodynamic conditions. The CO₂ ice albedo at the surface was set to 0.5 (Forget et al. 2013), although this value has a limited impact given that CO₂ ice forms in the nightside, where there is no starlight to reflect. The CO₂ aerosols thus formed were radiatively active, and the number mixing ratio of cloud condensation nuclei fixed to 10⁵ kg^-1, but tested the sensitivity by running simulations at 10² kg^-1 and 10⁸ kg^-1, and also by doing a test in which we turned off the radiative effect of CO₂ ice clouds.

2.1.5 Hazes and aerosols

For some simulations, we implemented and tested the effect of high-altitude aerosols. We first explored the effect of idealized hazes (Ducrot et al. 2025). For this, we added a haze term to the correlated-k table of a pure CO₂ atmosphere, of the shape $$\sigma =f_{\mathrm{h}\mathrm{a}\mathrm{z}\mathrm{e}}\kappa ({\lambda }_0/\lambda {)}^2(M_{\mathrm{C}{\mathrm{O}}_2}/{\mathcal{N}}_a)$$(1)

with κ = 0.5 cm²/g and λ₀ = 1 µm (Ducrot et al. 2025). The haze vertical mixing ratio profile was assumed to be constant. We varied different factors (haze mixing ratio profile, haze factor f_haze, see Equation (1)) in 1D (Appendix C.1).

We also explored more realistic optical properties for different hazes and aerosols (Table 3). In this case the radiative transfer code computes the aerosols’ contribution from their absorption coefficient, single scattering albedo and asymmetry factor. These parameters depend on the radius distribution of the particles. For both exo_k and the Generic PCM, or all aerosols, we prescribed an upper and lower pressure layer, and the aerosols follow the atmospheric scale height between them. The parameters (upper and lower pressure layers, volume mixing ratio, and size distribution of aerosols) are detailed is the dedicated sections (Section 3.3.2, Appendix C.1).

2.2 1D radiative-convective simulations with exo_k

We used the exo_k python library (Leconte 2021) that includes a 1D atmospheric evolution model (Selsis et al. 2023), in order to explore a wide range of possibilities in a limited amount of time. We performed ~80 simulations in which we varied atmospheric composition (different proportions of N₂, CH₄, CO₂, C₂H₄, C₂H₆, NH₃, SO₂) and total atmospheric pressure (from 0.1 to 10 bar). We used the same stellar, planetary and orbital parameters as described above for the Generic PCM (Table 2). We also used the same correlated-k tables as in the Generic PCM (Table 3), also with a spectral resolution of 10, ranging between 0.3 and 50 µm). We only tested dry atmosphere cases.

We parameterized the 3D redistribution for the incoming stellar flux with a multiplication factor f, that spans from 2/3 (no redistribution) to 1/4 (full redistribution): $$ISR=f\times \sigma T_s^4{\left(R_s/a\right)}^2$$

2.3 Post-processing tools

We computed the emission spectra of the modeled 1D atmospheres with exo_k, and for the 3D atmospheres we used Pytmosph3R (Falco et al. 2022), a python library based on exo_k that computes observables (eclipse spectra, phase curves) directly from the Generic PCM 3D outputs. When used in emission or phase curve mode, Pytmosph3R computes the emission for each column of the Generic PCM output, weighted by the angle of the observer. We used a stellar spectrum from the SPHINX model grid (Iyer et al. 2023) interpolated to TRAPPIST-1’s parameters. This spectrum, different from the one used in the GCM, was chosen for his closeness to the measured flux at 15 µm, allowing for a more accurate planet-to-star flux ratio – a criterion that is not sensitive in the GCM simulations, the stellar flux in the thermal infrared having a negligible effect on the GCM results.

We used correlated-k tables from the same opacity data (Table 3) at higher spectral resolution (R ~ 100) for the computation of the observables. To compare with the observations, we integrated the flux in the MIRI observed bands (F1280W and F1500W) with the available JWST MIRI filters documentation. We also computed the planet-to-star flux ratio by dividing the computed planet flux by the stellar flux. For this, we normalized the stellar spectrum with the 15 µm MIRI measurements of the stellar flux. At 12.8 µm, we find that the computed stellar flux is slightly lower than the measured flux (<2%).

2.4 Comparison to data

We compared our synthetic observables to the available JWST emission observations. First, to the eclipse depths of TRAPPIST-1 b at 15 µm (Greene et al. 2023) and 12.8 µm (Ducrot et al. 2025) (Fig. 3); second, to the 15 µm thermal phase curve of b extracted from the observed combined phase curve of TRAPPIST-1 b and c (Gillon et al. 2025). Gillon et al. (2025) performed four independent analyses, each one divided in several subanalyses with different assumptions. We chose to compare our results to analysis #1-MG because the analysis of MG was highlighted in the paper and only the subanalysis #1 allowed for an atmosphere on TRAPPIST-1 b. We note that contrary to Ducrot et al. (2025), our best-fit airless albedo is not 0.2, but ~0.4 (see Section 3.1.4). For the phase curve data, shown in Fig. 4, we chose to compare three parameters extracted from the phase curve: the maximal flux, the minimal flux, and the peak offset to evaluate the match between our simulations and the observation. The extracted phase curve of the planet b is quite dependent on the analysis pipeline – especially in this case where one must extract two signals (planet b and c) from one phase curve. Furthermore, any extracted phase curve is only one of the thousands of outputs of the MCMC analysis. On the other hand, the three variables extracted correspond to the posterior probability distribution functions of the analysis, and are thus much more robust. These quantities are also physically interpretable (amplitude corresponds to heat redistribution, offset to shift of the hottest hemisphere). We note that due to the complexity of the GCM simulations, we cannot perform a bayesian retrieval like with a simpler, 1D model. We simply tried to assess the match of our different simulations with the observations. For this, we computed the χ² as $${\chi }^2=\sum _i\frac{(O_i-M_i{)}^2}{{\sigma }_i^2}$$(2)

where O_i are the data points of the observations, σ_i their 1σ uncertainty, and M_i the models results.

For the eclipse, the two data points correspond to the eclipse depths measured at 12.8 and 15 µm ((Ducrot et al. 2025). For the phase curve, we extracted the dayside flux, nightside flux and peak offset from the analysis #1- MG of Gillon et al. (2025) and compared each of them.

We note that this χ² is a metric aiming at comparing the match of the different simulated atmospheres with the data, relatively to each other. We cannot provide an absolute quantification of the mismatch that could be compared with other papers, because our approach is different from 1D models, and we are not performing a proper fit, but rather evaluating a few realizations that we wanted to explore.

3 Results

We identify four distinct families of atmospheres that can produce a strong dayside emission matching the secondary eclipse MIRI observations (Greene et al. 2023; Ducrot et al. 2025). These four families are summarized in the sketches of Fig. 2 and listed below:

• Thin, residual atmospheres: the thin atmosphere produces very little heat redistribution and greenhouse effect. In this case, the infrared brightness temperature is very close to the surface temperature, where most of the emission comes from. Here the planet behaves very closely to a bare rock planet.

• Non-condensible, thick, transparent atmospheres: the observer sees through the atmosphere directly the surface emission of the planet. In the absence of condensation and radiatively active species, the surface and the atmosphere only interact through sensible heat exchange, which is limited. Thus, the atmosphere can produce very efficient heat redistribution but does not interact sufficiently with the surface to deviate the surface temperature field from the bare-rock case. Such transparent atmospheres are difficult to probe remotely, except perhaps through Rayleigh scattering at visible wavelengths.

• Thick, greenhouse efficient atmospheres: the atmosphere is thick, and produces a very efficient greenhouse effect. The large scale dynamics of such thick atmospheres produces efficient horizontal mixing, resulting in very efficient heat redistribution. As a result, the day-to-night temperature contrast is very low. In such thick atmospheres, the temperature of the lower atmosphere can far exceed the brightness temperature measured by MIRI, due to the powerful greenhouse effect of the atmosphere. For such atmospheres to fit the eclipse data, the molecules present in the atmosphere must have absorption windows that let through the emission from the deeper, warmer layers of the atmosphere.

• Atmospheres with a strong thermal inversion: the atmosphere possesses species (e.g., CH₄ or high-altitude hazes) that absorb efficiently at the wavelengths of the emission of the host star. A significant part of the stellar radiation is absorbed high in the atmosphere, producing strong shortwave heating, causing a thermal inversion. In this case, the CO₂ present in the atmosphere, with a strong absorption band at 15 µm, can now produce a strong signature in emission.

Fig. 3 shows the synthetic eclipse depths for four arbitrarily selected cases from the four families of atmospheric scenarios as a proof of concept. All of them can match the eclipse depth observations at the 2σ confidence level.

The following subsections describe each of these atmospheric scenarios, show why they match the eclipse depth measurements, and how phase curves can be used to discriminate between them.

Fig. 2 Sketches of the different atmospheric scenarios studied in depth in this work, as seen from the North pole.

Fig. 3 Dayside emission vs wavelength for a bare rock and for a representative of each family studied in this paper. The red data points and associated 1σ error bars are from Ducrot et al. (2025).

Fig. 4 Phase curves of the simulations from Fig. 3 and data from Gillon et al. (2025). The points in transparency correspond to the binned results of one of the MCMC outputs of analysis (#1-MG) of the paper. The nontransparent points are the maximum flux, minimum flux, and peak offset with 1σ uncertainties resulting from different analyses of the paper.

3.1 Residual thin and thick transparent atmospheres

3.1.1 Principle

The first two families of atmospheres that can match the measured eclipse depths are residual, thin atmospheres and thick, transparent atmospheres. To explore these scenarios quantitatively, we have decided to focus on CO₂ -N₂ -dominated atmospheres. We chose this composition for three distinct reasons: (1) first, they were initially proposed as good candidate molecules for a putative TRAPPIST-1 b’s atmosphere (Turbet et al. 2020a; Krissansen-Totton & Fortney 2022); (2) second, these are compositions that have been explored in detail in previous theoretical studies of TRAPPIST-1 b (Ih et al. 2023); and (3) third, because CO₂ is the strongest absorber at 15 µm and a condensable species, which is the most conservative case for our reasoning and makes it a good representative of this category. The principle and mechanisms at stake that we describe here can be transposed to other compositions of residual and/or transparent atmospheres (see Section 4).

To do this, we created a grid of 3D GCM simulations of TRAPPIST-1 b for a wide range of CO₂ and N₂ partial pressures. Our results are described in the following subsections.

Fig. 5 N2–CO2 grid of performed simulations. The points in red are not compatible with the observations; the green points match the data. The orange, purple, and brown lines correspond respectively to the 3σ limit of Ih et al. (2023), and limits from this work from eclipse and phase curve observations. The arrows indicate a collapse of the CO2 present in the atmosphere. The dots correspond to stable atmosphere, stars to unstable ones. The gray area indicates the atmospheres that are unstable and will collapse.

3.1.2 A grid of N₂₊CO₂ atmospheres

We performed simulations varying total pressure (0.01–1 bar) and CO₂ volume mixing ratio (VMR) (1 ppm to 100%), using a grid inspired from Ih et al. (2023). For each 3D simulation, we computed the eclipse depth spectrum, and the phase curve in the 15 µm (F1500W) MIRI filter. Compared to Ih et al. (2023), we added two extra points to specifically study the case of thick, “transparent” atmospheres with a low, constant CO₂ partial pressure of 10^-7 bar, for a total pressure of 1 and 10 bar A specific study of these cases is presented in Appendix A.1. Fig. 5 shows the grid of N₂+CO₂ atmosphere simulations we explored, and summarizes the various observational and theoretical constraints discussed in the following subsections.

3.1.3 Constraints from observations

Fig. 6 shows the observables resulting from the simulations of the grid. In order to compare them to the observations, we computed the reduced χ², provided in Table 1.

3.1.4 Eclipse depths

For the eclipse depths, the χ² is high for all the high-P_CO2, and drops for the atmospheres that correspond to the 3σ limit of Ih et al. (2023) (Fig. 5). It corresponds to all the cases that match the 15 µm data point within the 1σ error bar in Fig. 6. We find that the eclipse depth is mostly dependent on the CO₂ partial pressure, and could not distinguish the different thicker “transparent” atmospheres with constant CO₂ partial pressure but variable total pressure.

In this study, the emission for an airless planet of albedo 0.2 is higher than in Ducrot et al. (2025). In their work, the albedo is assumed to be wavelength-independent, with the emissivity defined as (ϵ_λ = 1 - A). As a result, the surface temperature, governed by the equation (1 - A)IS R = ϵσT⁴ (where ISR represents the incoming stellar radiation) becomes independent of the albedo, and the emission is directly proportional to the emissivity, (1-A). In contrast, with exo_k, we impose an albedo cutoff at λ ≥ 5 µm, which results in an emissivity ϵ_λ = 1 - A_λ equal to 1 in thermal wavelengths. λ = 5 µm was chosen to separate the zone of stellar emission and the zone of thermal planetary emission (Fig. 9). This leads to a surface temperature that depends on the albedo. Although the bolometric planet emission remains proportional to (1-A), thus maintaining the radiative balance, the emission at 15 µm is not. Finally, for the 3D simulations, the albedo is assumed to be wavelength-independent, but the emissivity is set to 1, which gives similar results to exo_k, as the albedo at longer wavelengths does not significantly influence the radiative balance, due to the low stellar emission at these wavelengths, and conversely, the emissivity is important only in the thermal wavelengths.

In consequence the surface albedo of 0.1, corresponding to the assumption made in simulations presented in Ducrot et al. (2025), results in a high eclipse depth for the thinnest atmospheres. The complexity and time-consumption of the GCM simulations does not allow us to explore widely the effect of the surface albedo, but we tested it in 1D for the different cases (see Appendix A.3). We found that a higher albedo lowers the continuum part (where we probe the surface), but not the CO₂ band when it is strong, because in this case we probe the emission directly from the atmosphere. We thus consider that the match with the data would be better with a higher albedo, but it would not change our overall conclusions.

Fig. 6 Synthetic eclipse spectra (left panel) and phase curves (right panel) for all the CO2–N2 atmospheres.

3.1.5 Phase curve

From the phase curve in the 15 µm band, as expected, the atmospheres that present a strong CO₂ feature in eclipse give a very low dayside flux. As we are probing in the atmosphere because of the CO₂ absorption, the redistribution is visible in these atmospheres through the non-zero nightside flux, and even a strong offset for the 1 bar 100% CO₂ case. These atmospheres are thus ruled out again, as expected.

The thinner atmospheres, with a higher dayside flux, low nightside flux and zero peak offset, can be considered an acceptable match, corresponding to the analysis in eclipse.

The only case where the phase curve brings a new constraint is for the thicker, more transparent atmosphere (see Appendix A.1). In this case, even with a very low quantity of CO₂, the weak continuum absorption (CIA and continua) is enough for the atmosphere to absorb and re-emit in the nightside. The probed atmosphere is well enough mixed that the nightside flux is high, and thus, not compatible with the observations.

3.1.6 Constraints from atmospheric collapse

TRAPPIST-1 b is so close to its host star that it is very likely to be tidally locked (Turbet et al. 2020a). Therefore, in the absence of atmosphere, its dayside is expected to be very hot and its nightside very cold. An atmosphere can redistribute the heat from the dayside to the nightside, with an efficiency that depends on its thickness, composition and, in general, on the parameters that can impact large-scale dynamics (Wordsworth et al. 2015; Koll 2022).

When reaching its condensation temperature, the gaseous CO₂ is expected to condense and form clouds, and/or frost on the ground. Therefore, if the nightside is cold enough, it will constitute a cold trap, collapsing the CO₂ of the atmosphere onto the surface (Wordsworth et al. 2015; Turbet et al. 2018).

The Generic PCM simulations we performed include CO₂ condensation, allowing us to test this prediction. We thus used the model to explore for which cases (pressure, CO₂ mixing ratio) the collapse happens, and if the collapse is total or if a secondary equilibrium is achieved (Soto et al. 2015).

Details on a pure CO₂ simulation collapse. We find that for a pure CO₂ case, the atmosphere collapses as soon as the initial pressure is 0.1 bar or below. This result is not impacted by the number mixing ratio of cloud condensation nuclei and radiative effect of CO₂ ice clouds. For this end member case, we performed the simulation of the collapse down to a few pascals (Fig. 7, solid black line), demonstrating that there is no secondary equilibrium possible. We computed the observables at different steps of the collapse, that are displayed in Appendix A.2.

Simulations showed that the ice is mainly produced as a result of frost formation and not precipitation. We note that it implies that the collapse strongly depends on the temperature of the surface at the cold traps, which has been shown to be strongly dependent on the strength of turbulence on the nightside (Auclair-Desrotour et al. 2022; Joshi et al. 2020).

Due to global dynamics, the cold traps are located close to the terminator, at mid-latitudes. As the atmosphere collapses, the greenhouse effect and redistribution become less efficient, the surface temperature drops, and the cold traps become more and more pronounced.

We tested the impact of the radiative effect of CO₂ ice aerosols parametrization. Changing the ratio of the cloud condensation nuclei and turning off the radiative effect of the aerosols did not change our results.

Collapse for a N₂-CO₂ atmosphere. For the N₂ atmospheres with 1% to 1 ppm of CO₂, partial pressure of CO₂ is lower so its condensation temperature decreases, and the N₂ participates in the redistribution of heat. Thus, the collapse is prevented for 0.1 bar atmospheres. For thinner atmosphere, when these processes are not sufficient to avoid collapse, CO₂ condenses, but not N₂. We called it a partial collapse (Fig. 5). In this case, the total pressure does not drop, but the CO₂ gas content drops. These atmospheres are thus collapsing towards an atmosphere with fewer CO₂ (for 0.1 bar total pressure, 1ppm is stable) eventually producing a transparent-type atmosphere, in better agreement with the JWST MIRI observations. Therefore, the case of a thin, N₂ atmosphere with a small amount of CO₂, which matches available data, can be considered as a resulting atmosphere of the CO₂-richer atmospheres, and not as implausible as it would look at first sight.

Fig. 7 Total pressure of a pure CO2 atmosphere against time, for the threshold case of a 0.1 bar atmosphere. In a tidally locked case, the atmosphere collapses totally in approximately ten Earth years. If we set a circulation of the substellar point to one Earth year, the 0.1 bar atmosphere becomes stable, and a partially collapsed atmosphere vaporizes all the CO2 back to the atmosphere.

3.1.7 Processes that can stop atmospheric collapse

Non purely synchronous orbit. Like a lot of exoplanets, TRAPPIST-1 b is close enough to its star to be expected to be in synchronous rotation due to tidal interaction with the star. But the TRAPPIST-1 system is so compact that planets gravitationally interact also with each other and this could result in a movement of the substellar point: libration, and a slow circulation of the substellar point. The circulation, if not too slow, can prevent collapse of the atmosphere by re-sublimating the ice caps when exposed to starlight. A 3:2 spin–orbit resonance could also prevent the collapse, but for TRAPPIST-1 b the tidal effects make it improbable (Turbet et al. 2020a; Revol et al. 2024). We tested this hypothesis by simulating a 0.1 bar atmosphere of pure CO₂ (collapsing in the tidally locked case), with a substellar point circulation (i.e., synodic day) of 1 Earth year. A longer duration would have been interesting to investigate but too computationally expensive and time-consuming. First, we started like the tidally locked case, with a warm (500 K), isothermal, dry atmosphere, and in another case, we started from a partially collapsed atmosphere. In both tests, the reached equilibrium state was a dry, non-collapsed atmosphere. The timescale for the collapse of a 0.1 bar atmosphere of pure CO₂ being around 10 Earth years (Fig. 7), a circulation of the substellar point with a rotation period of 1 Earth year is enough to prevent collapse. We can expect that the characteristic time of the substellar drift must be not longer than the one of the collapse in order to prevent it. The complete rotation of the substellar point has already been suggested by several studies. Makarov et al. (2018) suggested possible spin-orbit resonances state in the spin evolution of the planets b, d, and e, before synchronization. Chaotic circulation of the substellar points has been explored by Vinson et al. (2019), Shakespeare & Steffen (2023), and Chen et al. (2023). They attribute the large librations rising from the chaotic motion in their simulations, coupled with a simplistic model of tides (e.g., Constant Time Lag model (Hut 1981)). In addition, Revol et al. (2024) computed the synodic periods for each of the TRAPPIST-1 planets, finding a synodic day of about 70 Earth years for planet b.

In the end, for TRAPPIST-1 b, escaping atmospheric collapse thanks to the circulation of the substellar point is quite unlikely since simulations show that its circulation is not expected to be faster than one complete rotation of the substellar point per ~70 Earth years (Revol et al. 2024). However, the circulation of the substellar point is expected to be more important for the outer planets and should be taken into account for future GCM studies.

Tidal heat flux. On TRAPPIST-1 b, the star-planet tidal interactions may induce an internal flux. The eccentricity from Agol et al. (2021) constrains this flux between 4 × 10^-2 and 500 W m^-2. An internal heat flux can impact CO₂ collapse in two ways: 1) if it is strong enough, it can prevent any condensation at the surface, and 2) it can limit the height of the ice caps because of the basal melting, thus limiting the amount of CO₂ that can be condensed at the surface.

First, we can compute the equilibrium temperature at the nightside induced by an internal heat flux: $$\sigma T_{\mathrm{e}\mathrm{q}}^4=F_{\mathrm{i}\mathrm{n}\mathrm{t}}.$$

If this temperature is superior to the condensation temperature of CO₂, CO₂ condensation cannot occur, and thus no collapse. Thus, we can determine the minimum partial CO₂ pressure, corresponding to the pressure at which the condensation temperature is equal to the computed equilibrium temperature.

Second, the CO₂ ice cap can be limited by basal melting (Turbet et al. 2017). The maximum volume of CO₂ ice is limited to an ice cap of the surface of the nightside and the maximum height fixed by basal melting. Thus, the maximum pressure of CO₂ that can condense corresponds to $$P_{base}=\frac{1}{2}g\rho h_{max}$$

with g the gravity of the planet, and ρ the density of CO₂ ice.

We computed for different orders of magnitudes of the internal flux, the minimum pressure of CO₂ from 1), and the maximum condensed CO₂ from 2). Results are in Table 4.

From this table, we identified three regimes:

• (a)

  For a flux of ∼100 W m^-2 and higher, the condensation pressure is higher than the boundary pressure at which the atmosphere starts collapsing P_collapse. Therefore, there is no collapse possible.

• (b)

  For a flux of ∼1 W m^-2 and lower, the internal flux does not prevent condensation, and the maximum amount of CO₂ trapped correspond to a CO₂ pressure higher than P_collapse, and thus there is no constraint brought by the internal flux.

• (c)

  For intermediate cases around 10 W m^-2, process (2) tells us that the amount of CO₂ condensed is limited to ∼0.1 bar, which is below P_collapse ∼ 0.3 bar (for a pure CO₂ atmosphere). Thus, if there is an initial pressure between 0.1 bar and 0.3 bar, the atmosphere will not collapse completely, and there could be a remaining atmosphere, thinner than the initial limit P_collapse.

3.2 Thick greenhouse-efficient atmosphere

3.2.1 Principle

This family of atmospheres represents high-redistribution, high-greenhouse effect cases. To produce an efficient greenhouse effect, the gases present in the atmosphere must have efficient absorption in the thermal infrared, where the planet emission is strong. This way, the atmosphere should prevent the thermal radiation to directly escape to space. In other words, it means that the atmosphere becomes optically thick very high in the atmosphere, producing a brightness temperature representative of the higher, likely colder atmospheric layers, and not of the hot surface and/or deep atmosphere. The more absorbent the atmosphere is, the higher atmospheric layers we probe in secondary eclipses, and the colder the brightness temperature we see.

In principle, this family of atmosphere should therefore be incompatible with the measured secondary eclipses at 12.8 and 15 µm. However, this family of atmosphere can theoretically generate a very high flux at the wavelengths measured by MIRI, if these wavelengths coincide with atmospheric absorption windows. This also requires that the atmosphere does not absorb too efficiently in the near infrared to let the stellar flux reach the deep atmosphere. Fig. B.1 shows a proof-of-concept example where we set up an atmosphere for TRAPPIST-1 b with artificial opacities (high opacity set in the thermal infrared, but low opacities near 12.8 and 15 µm). This thought experiment, detailed in Appendix B.1, is designed here simply to illustrate the mechanism at play.

3.2.2 Gases with the right spectroscopic properties

We looked for a composition that could produce such an effect. We do not focus on the plausibility of the composition, as there are still few constraints regarding planetary evolution. We allow ourselves the freedom to explore all compositions, with the sole criterion being compatibility with the data. Among a list of possible molecules (CH₄, CO, CO₂, C₂H₂, C₂H₄, C₂H₆, HCN, H₂O, H₂, H₂S, H₂O₂, NH₃, N₂, PH₃, SO₂) we selected the species that had a weaker absorption around 15 µm. C₂H₄ and C₂H₆, i.e., reduced gases and potential photochemistry products of CH₄, are the only molecules with such a feature (Fig. B.2). Thus, to present this mechanism, an atmosphere must likely be reduced, and contain at least one of these gases. Other gases, transparent or with limited greenhouse effect, can be present in addition.

After an exploration in 1D with exo_k with these two gases and N₂ and CH₄, we identified a composition that maximized the dayside emission: N₂ as the background gas, CH₄ (20%) and C₂H₄ (0.2%). This composition has been chosen for the 3D Generic PCM simulations. We tried this composition with a total pressure of 10 bar.

We note that a larger quantity of CH₄ would make the atmosphere absorb too much in the near infrared, where TRAPPIST-1 emits a significant fraction of its flux, and act as an antigreenhouse effect, preventing the stellar flux to go deep in the atmosphere (see Section 3.3).

3.2.3 Constraints from the secondary eclipse

We show in Fig. 8 the results and observables of the simulation of N₂–CH₄–C₂H₄, with P_tot = 10 bar. Due to the strong greenhouse effect, the surface of the planet is quite hot: more than 600 K (compared to the equilibrium temperature of ∼400 K). The thermal profiles of the dayside and nightside are superimposed in the deep atmosphere, because of the highly efficient redistribution. In the upper layers, the redistribution is less efficient, and the thermal profiles of dayside and nightside diverge. The layers from which the main part of the 15 µm band pass emission comes from (shown in shadowed area) are between 1 bar and 0.01 bar, because even with the opacity “window”, the atmosphere becomes optically thick due to the opacity scaling with atmospheric pressure, and one cannot probe all the way down to the surface. Thus, the eclipse depth is around 600 ppm for 15 µm, corresponding to a brightness temperature of only ∼425 K, and 450 ppm/∼415 K for 12.8 µm.

3.2.4 Constraints from the phase curve

Still in Fig. 8, the wind and temperature map at the emission layers shows strong eastwards winds, moving heat towards the nightside of the planet. These winds have also the effect of moving the hottest point eastwards the substellar point.

The modeled phase curve reproduces these features: it is globally flat, because of heat redistribution, and the maximum of the phase curve is shifted slightly ahead of the time of the eclipse, because of the shift of the hottest point. In conclusion, such an atmosphere is compatible with the eclipse depth measurements, but is ruled out by the phase curve.

3.3 Atmospheres with a strong thermal inversion

3.3.1 Principle

In first approximation, Earth’s atmosphere is transparent in the visible, and opaque in the IR, allowing most of the light from the Sun to go through all the atmosphere, warming the surface, which re-emits this energy in the infrared, and the atmosphere itself absorbs this light. This process tends to make the atmosphere warmer close to the surface, and colder as we go up in the atmosphere. On Earth, the ozone layer absorbs in the UV, heating the atmosphere directly from the top. It creates a so-called thermal inversion: the atmosphere goes hotter as we go up, in what is called the stratosphere.

In the case of TRAPPIST-1 b, the star is colder than the sun, and its spectrum is thus shifted towards infrared wavelengths. A significant part of the emission of the star is in the near infrared. Thus, a strong thermal inversion can happen quite easily, if the atmosphere contains some kind of absorber in the stellar emission range. Aerosols, but also molecules like CH₄, can absorb in the near infrared. In this section, we explore thermal inversion created by different kinds of components.

If there is a thermal inversion, a molecule with an absorption band at the observed wavelength will result in observations probing higher in the atmosphere, as usual, but this time, it means that one will probe higher temperatures. The absorption band of the molecule is thus expected to be observed as an emission band (Fig. 9).

Fig. 8 Thick atmosphere case (10 bar, N2 atmosphere with 20% of CH4 and 0.2% of C2H4). Upper panels: synthetic observables and observed data (eclipse depth vs wavelength, left, and phase curve, right). Lower panels: temperature profile (left) and map (right). The temperature map corresponds to a sum of the atmospheric layers, weighted by their contribution to the emission at 15 µm. We added the similarly weighted wind vectors as black arrows. The shadowed area in the right panel corresponds to the layers contributing to 95% of the emission. The difference between the 1D and 3D temperature profiles is discussed in Appendix B.3.

3.3.2 Composition

Thermal inversion with gas. We tested an atmosphere with a high amount of CH₄, to create the inversion, and some CO₂, to be seen in emission. We supposed a redistribution factor of 1/4 (total redistribution) for the 1D simulations. The results of the best case in 1D are shown in Fig. 9. The chosen composition is N₂ as background gas, 40% of CH₄, and 0.4% of CO₂. The 3D study of this atmosphere is shown in Fig. 10. The eclipse depth in 3D is much higher than in 1D, because the redistribution in the probed region is not total, contrary to the assumption made in 1D. We note that such an amount of CH₄ is expected to produce hazes in the upper layers, and that these hazes are expected to contribute to the thermal inversion.

Thermal inversion with simplified hazes. Our simplified model for hazes, from Ducrot et al. (2025), is described in 2.1.5. We find that using the value of f_haze = 7.0 × 10^-4 (Equation (1)) as in Ducrot et al. (2025), we probe a very high-altitude region of the atmosphere, where the heat redistribution is inefficient. Despite the strong winds, at such low pressure the radiative timescale (${\tau }_{\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{C_{p}P}{g\sigma T^3}\sim {10}^2$ s at P=10 Pa, with C_p the heat capacity of the atmospheric gas, g the gravity of the planet and σ the Stefan-Boltzmann constant) is very short compared to the advection timescale (${\tau }_{\mathrm{a}\mathrm{d}\mathrm{v}}=\frac{R_p}{v}\sim {10}^4$s, with R_p the radius of the planet and v the wind horizontal speed). Therefore, the dayside at this height is much hotter than the 1D prediction that assumed a full redistribution. In order to match the data we thus tried two different options: 1. Add a simplified single-scattering albedo of 0.5 by lowering the incoming stellar flux of a factor 2 (Model 1, corresponding to Haze high in Gillon et al. (2025) and 2. Lower the haze factor f_haze to f_haze = 3.0 × 10^-5, while adding a small simplified single-scattering albedo of 0.2 (Model 2, corresponding to Haze low in Gillon et al. (2025). These cases are shown in Fig. 11.

Thermal inversion with more realistic aerosols. We also tested for aerosols for which we had more realistic optical properties: Titan-like tholins, H₂ SO₄, and martian dust. A 1D study of the impact of the different parameters can be found in Appendix C.1. For all aerosols tested, the single-scattering albedo was high in the near infrared, preventing a strong heating of the high atmosphere. Thus, a strong quantity of aerosols was needed to get a strong flux at 12.8 and 15 µm. The most efficient aerosols for the thermal inversion are the largest, and only the dust particles allowed for a strong flux at 15 µm compatible with the eclipse depth. For all the cases tested, the aerosols are very high and in large quantities, making the high atmosphere opaque at all wavelengths. At these low pressures, the CO₂ does not have a strong radiative effect, and is not seen in emission. Instead, the emission from the atmosphere corresponds to a blackbody at the temperature of the high atmosphere. We selected two cases to test in 3D: first, a layer of dust (with an effective radius of 100 µm) between pressures 1000 and 1 Pa, with a constant volume mixing ratio (VMR) of 3.5 × 10^-18, and second, a layer of tholins (effective radius 1 μm, between 10 and 1 Pa, VMR=1. × 10-13). We note that for the dust case, the opacity is so important that there is not enough flux to go through the atmosphere, resulting in an isotherm rather than a thermal inversion (see Appendix C.1).

Fig. 9 1D case of thermal inversion (10 bar, N2+40% CH4+0.4%CO2). Left: opacity of the atmosphere with the zone of emission of the star and the planet. Middle: temperature profile and the layers of emission of the atmosphere at 15 µm. Right: eclipse depth. The methane absorbs partly in the star-emission region, creating the thermal inversion in the temperature profile. The CO2 absorbs in the 15 µm band, and thus the emission layers at 15 µm are high enough in the atmosphere, where the temperature is high. Thus, the eclipse depth is higher.

Fig. 10 Thermal inversion: 1 bar, N2+40% CH4+0.4%CO2 . The panels are the same as in Fig. 8.

Fig. 11 Thermal inversion: CO2+ simplified hazes. The panels are the same as in Fig. 8. Model 1 and Model 2 for the hazes are described in the text. The temperatures and winds of the lower panels correspond to the Model 1 case, which is the model closest to the data.

3.3.3 Constraints from the secondary eclipse

For the inversion with gas and with simplified hazes, one can see that contrary to all other simulations, the CO₂ band is not seen in absorption, but in emission. This is consistent with the eclipse observations, since the 15 µm brightness temperature (478 ± 27 K) is higher than the 12.8 one (424 ± 28 K). For the inversion with realistic aerosols, the CO₂ is not visible as the atmosphere becomes optically thick too high.

For the simplified hazes case, the haze factor and singlescattering albedo were fit to obtain a satisfying match with the observations (Fig. 11). This case, although fine-tuned, gives the best agreement to the data (Table 1).

For both cases from aerosol optical properties (dust and tholins), the aerosol layer in the high atmosphere is so opaque that it mimics a surface, and the emission spectrum is similar to a bare rock (Fig. 12).

Fig. 12 Thermal inversion: CO2+aerosols with realistic optical properties. The panels are the same as in Fig. 8. The parameters for dust and tholins are described in the text. The temperatures and winds of the lower panels correspond to the dust case, which is the model closest to the data. The difference between the 1D and 3D thermal profiles is discussed in Appendix C.2.

3.3.4 Constraints from the phase curve

For all cases, we probe high enough in the atmosphere to have a non-flat phase curve: despite the efficient mixing in the deeper atmosphere, the upper atmosphere presents a hotter dayside because the heat redistribution is less efficient there. For the CH₄ case (Fig. 10) and the simplified haze Model 2 (Fig. 11), though, the probed pressures includes part of the well-mixed atmosphere, resulting in a low amplitude, high minimum temperature in the phase curve, but also a strong offset, because of the eastwards shift of the hottest point due to the super-rotation.

For the simplified haze Model 1 (Fig. 11) and more realistic aerosols, the probed region is so high, that all the redistribution is inefficient there, and the phase curve presents a quasi-zero nightside flux, a very high dayside flux, and no peak offset, thus matching the observations if the upper atmosphere is hot enough. Thus, to be compatible with the observations, the thermal inversion must be strong, and one needs a strong enough absorption at 15 µm to probe very high in the atmosphere. The hazes compatible are therefore abundant very high in the atmosphere and very dense. For the realistic aerosol optical properties (Fig. 12), the probed temperature is slightly too low to match the data. The realistic aerosols have a high single-scattering albedo, reflecting part of the stellar flux in space, so the upper layers temperature is lower than in the simplified cases, and the emission is weaker.

In summary, the emission signature is influenced by the type, abundance, and location of the absorbing species within the atmosphere. Three cases stand out: (1) When the stellar absorber is weak, such as CH₄, the emission originates from deeper atmospheric layers, with CO₂ detected in emission but resulting in a phase curve with low amplitude. (2) For more efficient absorbers, such as those in the simplified haze Model 1, the emission comes primarily from the upper atmosphere, leading to a phase curve with high amplitude mimicking the bare rock signal. (3) Finally, if the stellar absorber is highly efficient and abundant, it can render the entire atmosphere optically thick. Under these conditions, CO₂ will no longer be observed in emission, and the resulting emission signature will resemble that of an airless rock.

4 Discussions

4.1 Limitations in the model

The GCM allows for more realistic simulations, but they are too complex and expensive to pretend to explore the whole parameter space of exoplanet atmospheres. Among the parameters that we arbitrarily fixed and the hypotheses made, we note some that could have an impact on our results:

• The internal heat flux of the planet, set to zero in our simulations, could prevent the atmospheric collapse for thin atmospheres, as already discussed in Section 3.1.6;

• The surface albedo changes the flux at 15 and 12.8 µm, in particular for optically thin atmospheres. We discuss its impact on the eclipse depth in Appendix A.3;

• The hazes we added in our simulations are highly idealized: we prescribed the vertical position and quantity of the hazes in the atmosphere. A more in-depth study of the formation of hazes would be necessary to assess whether the assumptions we made about haze properties are possible;

• Moreover, a CO₂ -CH₄ -rich atmosphere is not expected to be thermochemically stable (Krissansen-Totton et al. 2018; Watanabe & Ozaki 2024). However, for this study we only tried to identify gas mixtures matching the data, without any prior on the expected composition. We recall that this is the first time as many observations are obtained on this category of planet (relatively temperate, Earth-size exoplanet around M-star), and leave the question of plausibility for future work;

• Apart from the N₂–CO₂ case, we did not study the formation and impact of clouds on our simulated atmospheres. At the pressure and temperature range studied, H₂O is the only condensable species. Condensation of water around TRAPPIST-1 planets have been studied already in Turbet et al. (2023).

4.2 Limitations in the data analysis approach

For a transit spectrum or a secondary eclipse, one fits the timedependent data with a well-known light curve shape, and extracts the transit or eclipse depth from it. The extraction of the phase curve is not as straightforward. The shape of the phase curve depends on the properties of the planet and its hypothetical atmosphere. Thus, one must make physical assumptions regarding the planet to parameterize the shape of the phase curve, so the data reduction is not as well decoupled from the physical interpretation as in the case of transits and eclipses.

Moreover, in our case, the observed phase curve from Gillon et al. (2025) contains the contribution of TRAPPIST-1 b and c. Thus, it is delicate to extract the contribution of TRAPPIST-1 b alone in our study. The most relevant approach would be to perform a joint fit of the two planets, but this is out of the scope of this study.

4.3 Comparison with previous studies

Lim et al, 2023. We computed synthetic transit spectra for the atmospheres that match all the data in emission, in order to check that these are compatible with the information available from transmission (see Appendix D). Our signal is typically below 100 ppm which is undetectable from the data of Lim et al. (2023).

Thus, these scenarios can be neither ruled out nor verified from transmission spectroscopy information.

However, it is important to keep this in mind for future observations of TRAPPIST-1 system (de Wit et al. 2024), in particular observation strategies using TRAPPIST-1 b considered airless as a probe to correct for stellar contamination (JWST GO 6456).

Ih et al. (2023). First, we compared our eclipse spectra from Fig. 6 with the ones from Ih et al. (2023) at same pressure and composition and obtained good agreement for atmospheres with little CO₂. For the atmospheres that are optically thick (e.g., pure CO₂ atmospheres), the flux in the CO₂ band is lower in our simulations, probably because of differences between our temperature profiles (the higher atmosphere is colder in our simulations).

Though we studied in depth the condensation of CO₂, we did not perform simulations for the other compositions in Ih et al. (2023). N₂, CO, O₂ and CH₄ should not condense under TRAPPIST-1 b conditions. SO₂ and H₂O, though, are more condensable than CO₂ (Fray & Schmitt 2009; Turbet 2018), and Ih et al. (2023) found that a 1 bar O₂ atmosphere with 100 ppm H₂O and a 0.1 bar O₂ atmosphere with 100 ppm SO₂ would fit the observations. Considering the surface temperature in our 100 ppm CO₂ simulations at these pressures, we find that the minimum surface temperature are below the condensation point of SO₂ and H₂O, and therefore most of the SO₂ and/or H₂O is expected to collapse for a tidally locked case.

Ducrot et al. (2025). In Ducrot et al. (2025), the authors introduced the hypothesis of a thick, hazy atmosphere that allowed seeing CO₂ in emission. They expected an efficient redistribution for such a thick atmosphere. We find that with such an atmosphere, the redistribution is indeed efficient in the deep atmosphere, but in the higher, less dense atmosphere, which is actually the part that is probed at 15 µm, the heat redistribution becomes less efficient because of the very short radiative timescale at low pressure, allowing for a much brighter dayside, and a high amplitude phase curve (Fig. 11), an effect which was not represented by the 1D study.

4.4 Perspectives for the study of TRAPPIST-1 b

We highlighted different atmospheres that could produce observables compatible with our data available. Are there additional emission observations of TRAPPIST-1 b that could distinguish between the remaining cases? The thin, residual atmospheres are by construction very similar to a bare rock, and cannot be distinguished by observation. For the atmospheres with thermal inversion, depending on the aerosols, we saw that there can be a CO₂ emission band or no spectral features in emission. Only the case of a thermal inversion with CO₂ in emission could be investigated by more observations.

A phase curve at 12.8 µm could in principle help distinguish between some thermal inversion scenarios (with a lighter amplitude at 12.8 µm than 15 µm) and the airless case. However, the hazy atmosphere is still not transparent at these wavelengths, thus not resulting in a totally flat phase curve, and simulations of observations show that the signal-to-noise ratio would be too weak to distinguish the different scenarios (Fig. 13).

Secondary eclipse observations far from the CO₂ band, where some hazy models are expected to give much less flux than airless models, could be interesting. We used the JWST Exposure Time Calculator (ETC) to compute the expected precision on the occultation depth of TRAPPIST-1 b at 10 and 18 µm. Using the following formula we can compute the precision on the occultation depth: $$\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\frac{\sqrt{(1/N_{in}+1/N_{out})}}{S/N_{point}}.$$(3)

Here N_in is the number of points during the occultation and N_out is the number of points outside the occultation, during a single visit.

For MIRI filter F1800W with SUB256 (to maximize the number of groups per integration) we obtain a S/N per point equal to 600 and an exposure time equal to 89.56 s for 299 groups per integration. We compute a precision of 371 ppm for one visit. For F1000W we compute a precision of 97 ppm in one visit (S/N of 868.49 and exposure time of 13.18s). Fig. 14 shows the predicted observations for five visits. We compute that we would need at least eight visits to distinguish between a bare-rock scenario and the hazy, CO₂ emission case with an observation at 18 µm, but only two visits at 10 µm. This is to distinguish between two extreme cases, however, a hazy atmosphere emission can stand between these two cases, and even mimic an airless scenario (see Section 3.3.3).

Fig. 13 Phase curves in the 12.8 µm MIRI filter for an airless body and the remaining atmospheric families, along with the corresponding simulated observational points.

Fig. 14 Dayside emission for an airless body and the remaining atmospheric families. The red points are from Ducrot et al. (2025). The purple points are predicted observation points for five visits for an airless body of albedo 0.2.

5 Conclusion

Starting from the secondary eclipse observations, we performed a study of the atmospheres on TRAPPIST-1 b that would produce a high brightness temperature at 12.8 and 15 µm. After exploring the parameter space in 1D, we found several families of atmospheres that do so, and selected representatives of each family to study in 3D. The 3D study brought several results that could not be predicted in 1D:

• Thin condensable atmospheres: we identified the process of atmospheric collapse, through cold trap condensation, that allowed us to rule out some atmospheres that are thus not stable. Cross-checking this information with the observations allowed us to narrow the window of possible thin N₂ -CO₂ atmospheres (Fig. 5). This result can be extended to other compositions (Section 4.3);

• Thick transparent atmospheres: we found that even a very low opacity, such as N₂-N₂ collision-induced absorption, is enough to make the atmosphere radiatively active and emit a flux in the nightside, making such an atmosphere incompatible with the phase curve data (Appendix A.1);

• Thick greenhouse-efficient atmospheres: we highlighted a type of thick highly redistributive atmosphere that gives a dayside emission that is similar to bare rock, but that could be ruled out with the phase curve observation (Fig. 8). Atmospheres with these characteristics are reduced atmospheres (containing CH₄ and C₂H₄ or C₂H₆) in order to provide an atmospheric window for radiation to escape from deeper layers to give the high flux observed by JWST photometry;

• Atmosphere with thermal inversion: for this kind of atmosphere already proposed in Ducrot et al. (2025), we found that the 1D hypothesis of a full redistribution is not accurate since the flux is coming from high in the atmosphere where the heat redistribution becomes inefficient. Furthermore, this effect – if the hazes have the right altitude and properties – can result in a phase curve with a very strong day–night amplitude and no offset, compatible with the observations (Fig. 11). Only atmospheres with a thermal inversion and hazes or aerosols can match the phase curve data.

In summary, we showed that an airless planet and a thin residual atmosphere are not the only scenarios compatible with the secondary eclipse observations. We proposed three families of thick atmospheres (transparent atmospheres with almost no condensing or absorbing species; reduced atmospheres with an opacity window at 15 µm; and atmospheres with a thermal inversion) that also match the eclipse data. The phase curve observation showing a very cold nightside and no peak offset was needed to rule out some cases (thick reduced atmospheres; thick transparent atmospheres; and atmospheres with a thermal inversion but no hazes). Taking into account all the emission observations available, we showed that TRAPPIST-1 b might be a bare rock or possesses a very thin atmosphere poor in absorbing and/or condensing gases; alternatively, it could possess an atmosphere with a thermal inversion and hazes, although this case seems unlikely and fine-tuned.

Though this study was dedicated to TRAPPIST-1 b, as this planet had an unprecedented amount of observed data in emission (eclipse depths at two different wavelengths and phase curves), the methodology of this study is applicable to other rocky targets to be observed with JWST. Moreover, our results (e.g., atmospheric collapse, question of heat redistribution) can qualitatively apply to other close-in rocky exoplanets observed in emission. In particular, we call attention to the degeneracy of a high eclipse depth at 15 µm, making it hazardous to conclude that a planet is airless just from this photometric point.

Acknowledgements

This work has made use of the Infinity Cluster hosted by Institut d’Astrophysique de Paris; we thank Stéphane Rouberol for his work on the cluster. We acknowledge support from the Tremplin 2022 program of the Faculty of Science and Engineering of Sorbonne University. We acknowledge support from the High-Performance Computing (HPC) resources of Centre Informa-tique National de l’Enseignement Supérieur (CINES) under the allocations Nos. A0140110391 and A0160110391 made by Grand Équipement National de Calcul Intensif (GENCI). We thank the Generic PCM team for the teamwork development and improvement of the model. We acknowledge support from BELSPO BRAIN (B2/212/PI/PORTAL).

References

Appendix A CO₂-N₂

A.1 Thick transparent atmosphere

We studied the influence of P_tot for a constant partial pressure of CO₂ in order to test the thick-transparent atmosphere case. Fig. A.1 shows the observables for different tests of atmospheres from the “transparent” family, and the temperature distribution for the 10 bar, N₂+0.1ppm CO₂ atmosphere. One can see that the surface is well decoupled from the atmosphere, since the temperature profile of the atmosphere is almost the same in the nightside and dayside, although the surface presents strong differences between day and night (lower left panel).

For a fully transparent atmosphere, the heat redistribution should not be detectable in the phase curve as the emission comes from the surface. Yet, atmospheres composed of real gases can not be completely transparent, and even a small opacity allows for the atmosphere to absorb and re-emit part of the radiation from the surface, where the heat is redistributed and the nightside warmer. Thus, one can see that the 10 bar atmosphere phase curve has a weaker amplitude, with a nightside flux of almost 400 ppm. Part of the absorption is caused by collision-induced absorption of N₂, and part of it by the remaining CO₂ present in the atmosphere. When we remove these contributions (“transparent” case in the figure), the phase curve is similar to a thin atmosphere.

Fig. A.1 Upper panels: Observables of the (“transparent” atmospheres. The CO2 partial pressure is 1.e-7 bar, except for the (“transparent” case where the observables have been calculated for PCO2 =1.e-15 bar. For the “no CIA” case, the collision-induced absorptions have not been included in the opacity computation. Lower panels: Temperature profile and map for the 10 bar case. The atmosphere is very well mixed, but there is no efficient coupling with the surface, and the surface has a hot dayside and cold nightside.

A.2 Pure CO₂

We show in Fig. A.2 the observables and temperature for the case of a collapsing, pure CO₂ atmosphere, with the boundary of the CO₂ ice cap.

Fig. A.2 Pure CO2 case. Upper panels: Observables (right: phase curve at 15 µm, left: eclipse depth). The dotted lines represent non-stable atmospheres; these are snapshots of the atmosphere during its collapse. Lower panels: Physical quantities. Temperature profiles and temperature map for the 10 Pa case (compatible with the observations, but not stable). In the lower right panel, the red contour shows the CO2 ice caps. The red dot is the substellar point. In the lower left panel, the red line corresponds to the CO2 condensation curve.

A.3 Effect of surface albedo

We studied the effect of the surface albedo for several atmospheres from the CO₂ -N₂ grid in Fig. A.3. If the CO₂ absorption band is strong (a,c), the atmosphere is optically thick at 15 µm and the surface albedo does not matter at this wavelength. For cases with little CO₂ proportion (b,d), there is no saturation effect and the albedo lowers the whole spectrum equally. We can then suppose that for low CO₂ proportion, increasing the albedo can decrease the eclipse depth for 12.8 and 15 µm, although for high CO₂ proportion, increasing the albedo would decrease the eclipse depth at 12.8 µm, but not at 15 µm.

Fig. A.3 Effect of the surface albedo for different cases of N2-CO2 atmospheres (results of 1D simulations).

Appendix B Thick atmosphere with an opacity window

B.1 Proof of concept

Figure B.1 shows a 1D proof-of concept of the family of the thick atmospheres. We designed an unrealistic opacity of the atmosphere meeting our requirements: a high opacity in the range of thermal emission of the planet, a window around the probed wavelengths, and a weaker opacity at the range of emission of the star, to let the stellar flux reach the surface of the planet. We ran the 1D radiative transfer code exo_k and tuned the values of the opacity to reach the eclipse observation points. Thanks to the absorption window, the main part of the emission at 12.8 and 15 µm comes from deep in the atmosphere, where the temperature corresponds to the measured brightness temperature of the observations.

Fig. B.1 1D study of thick atmosphere with greenhouse effect with an opacity window at the observed wavelengths. Left: Idealized opacity vs wavelength. The yellow window corresponds to the zone of emission of the star, the red window to the planet. Middle: Temperature profile. The window corresponds to the layers of 95% of the emission. Right: Corresponding eclipse depth.

B.2 Opacities

We show in Fig. B.2 the opacities of different molecules. We look for molecules that absorb in most of the planet emission range, but less in the MIRI filters at 12.8 and 15 µm. C₂H₄ and C₂H₆ present these features.

Fig. B.2 Opacities of molecules. The dotted lines show the MIRI filters where the molecules must not absorb too much in order to produce the right effect. The planet blackbody emission range is in red.

B.3 3D effect

We provide the thermal profiles in 1D and 3D of different simulations in Fig. B.3 in order to illustrate the 3D effect resulting in a mean thermal profile colder than 1D computations. In 1D, the deep atmosphere is convective, and the thermal profiles follows an adiabat. In addition to the usual exo_k simulations, we performed two 1D PCM simulation, one starting from the exo_k thermal profile (hot start) and the other from the 3D mean profile (cold start), in order to make sure the difference between the 3D Generic PCM simulations and the 1D simulations comes from 3D effects. The thermal profile of the 1D PCM indeed converges towards a profile close to the exo_k profile.

In 3D, the deep atmosphere is well mixed horizontally, but the poles are colder than the equator region. This leads to subsidence at the poles, forcing an adiabat that is colder than in the 1D case. This changes the circulation regime, counter-intuitively making the surface temperature colder at the equator than at the poles, and preventing convection at the equator region.

Fig. B.3 Thermal profiles for different simulations of a 10-bar N2 atmosphere with 20% CH4 and 0.2% C2H4.

Appendix C 1D study of hazes

C.1 Parameter study

We conducted a study in 1D with exo_k. First, we varied different parameters for of the simplified model presented in 2.1.5. Then, we compared the model with different aerosol optical properties (extinction coefficient, single scattering albedo, and asymmetry factor) from laboratory measurements, and finally we varied some parameters for two selected aerosol types. The parameters of each 1D simulations are summarized in Table C.1.

1. Impact of the vertical distribution profile of the hazes (Fig. C.1a). We find that the most important is that the hazes are present at the top of the atmosphere. If the upper layer of the atmosphere is hazeless, the CO₂ present in this layer is seen in absorption. The presence or absence of hazes in the deeper atmosphere is not important.

2. Impact of the factor f_haze (Fig. C.1b). This factor corresponds to the quantity of hazes present in the atmosphere. We find, like Ducrot et al. (2025), that 7.e-4 is optimal, as a lower quantity of hazes results in a thermal inversion deeper in the atmosphere, and the CO₂ band probes in the upper atmosphere. Thus, we confirm that the case we present is fine-tuned.

3. Comparison with real aerosol optical properties (Fig. C.1c). For all the different aerosols that we could test, we do not get a CO₂ emission band. Martian dust results in a thermal inversion allowing for a high eclipse depth.

4. Impact of the single-scattering albedo (Fig. C.1d). With a single scattering albedo fixed to 0.5, independent of the wavelength, we can reproduce some features of our simplified model.

5. Impact of the size of the aerosols (Fig. C.1e). We find, for two different types of aerosols (tholins and dust), that the largest aerosols give inversions upper in the atmosphere, and thus highest eclipse depth.

6. Impact of the volume mixing ratio (Fig. C.1f). We varied the volume mixing ratio (VMR) for the same two types of aerosols. For the dust case, the highest VMR leads to the highest eclipse depth, although there is a saturation at ~1 × 10^-15. For tholins, for VMR ~1 × 10^-12 there is an inversion, but it is not strong enough to give a high eclipse depth.

With these results, we identify that the realistic single scattering albedo of the aerosols, that is very high in the shortest wavelengths (Fig. C.2), is responsible for the absence of CO₂ seen in emission by reflecting a very large part of the stellar emission, prohibiting a strong heating of the upper atmosphere.

Fig. C.1 Temperature profile and eclipse depth for different parameterizations of aerosols. Each panel illustrates the influence of a parameter.

In all our cases, either there is no strong thermal inversion in the upper layers and CO₂ is seen in absorption, or the hazes seem to be optically thick at all wavelength, resulting in an isotherm in the deep atmosphere instead of a thermal inversion, and a thermal emission that reproduces a blackbody radiation with the temperature of the upper atmosphere.

In conclusion, only dust type aerosols are able to reproduce a strong eclipse depth, and no aerosol reproduces the CO₂ emission band that is seen with the simplistic haze model.

Fig. C.2 Single scattering albedo vs wavelength for different aerosol types. The radii of aerosols are 10 × 10 6 for tholins and H2SO4 and 10 × 10 4 for dust.

C.2 Convergence issues

Some of these simulations contain a high amount of aerosols, and only a very small part of the radiation goes through the atmosphere. Therefore, the deep atmosphere is getting almost no flux, and its convergence is extremely long. The time step of the simulations is necessarily short, because of the small radiative timescale at the high atmosphere, whereas in the deep atmosphere the evolution timescale is way longer than the time step. exo_k comes with several acceleration routines, making it possible to converge this type of atmosphere. However, the Generic PCM does not, and our 3D simulation of the dust scenario is not converged in the deep atmosphere. This does not impact our conclusion on observables as the thermal emission comes from the upper parts of the atmosphere, that are converged. Fig. C.3 shows that the 1D exo_k profile is similar, when not converged, to the 3D profile.

Fig. C.3 Thermal 1D profile for the dust case, converged vs not converged. The initial profile is a 370 K isotherm.

Appendix D Transit spectra

We computed transit spectra for representatives of the families that match the eclipse depths and phase curve: N₂-CO₂ atmospheres (here, 1 bar with 1ppm of CO₂), and atmospheres with thermal inversion (here, the hazes from the simplified hazes Model 1). The hazes distinguish from the bare rock and N₂-CO₂ case by a slope in the shortest wavelengths. However, this feature is of the order of dozens of ppm which is for now too dim to be observed. Nevertheless, these features could be of importance for future observations of the system, especially when TRAPPIST-1 b is used as a probe to correct for stellar contamination.

Fig. D.1 Transit spectra computed for scenarios that match all the available emission observations. The transit spectrum has been offset to have a minimum value of 0.

All Tables

All Figures

Fig. 1 Exoplanets under consideration for the Rocky Worlds DDT. Planets above the cosmic shoreline are expected to have lost their atmosphere. The circled exoplanets have already been observed in emission during JWST’s first three cycles. The exoplanets for which the observations concluded there is no atmosphere are in red. We also added the Solar System terrestrial planets. The XUV irradiation is computed from Zahnle & Catling (2017).

In the text

	Fig. 2 Sketches of the different atmospheric scenarios studied in depth in this work, as seen from the North pole.
In the text

	Fig. 3 Dayside emission vs wavelength for a bare rock and for a representative of each family studied in this paper. The red data points and associated 1σ error bars are from Ducrot et al. (2025).
In the text

	Fig. 4 Phase curves of the simulations from Fig. 3 and data from Gillon et al. (2025). The points in transparency correspond to the binned results of one of the MCMC outputs of analysis (#1-MG) of the paper. The nontransparent points are the maximum flux, minimum flux, and peak offset with 1σ uncertainties resulting from different analyses of the paper.
In the text

Fig. 5 N2–CO2 grid of performed simulations. The points in red are not compatible with the observations; the green points match the data. The orange, purple, and brown lines correspond respectively to the 3σ limit of Ih et al. (2023), and limits from this work from eclipse and phase curve observations. The arrows indicate a collapse of the CO2 present in the atmosphere. The dots correspond to stable atmosphere, stars to unstable ones. The gray area indicates the atmospheres that are unstable and will collapse.

In the text

	Fig. 6 Synthetic eclipse spectra (left panel) and phase curves (right panel) for all the CO2–N2 atmospheres.
In the text

	Fig. 7 Total pressure of a pure CO2 atmosphere against time, for the threshold case of a 0.1 bar atmosphere. In a tidally locked case, the atmosphere collapses totally in approximately ten Earth years. If we set a circulation of the substellar point to one Earth year, the 0.1 bar atmosphere becomes stable, and a partially collapsed atmosphere vaporizes all the CO2 back to the atmosphere.
In the text

Fig. 8 Thick atmosphere case (10 bar, N2 atmosphere with 20% of CH4 and 0.2% of C2H4). Upper panels: synthetic observables and observed data (eclipse depth vs wavelength, left, and phase curve, right). Lower panels: temperature profile (left) and map (right). The temperature map corresponds to a sum of the atmospheric layers, weighted by their contribution to the emission at 15 µm. We added the similarly weighted wind vectors as black arrows. The shadowed area in the right panel corresponds to the layers contributing to 95% of the emission. The difference between the 1D and 3D temperature profiles is discussed in Appendix B.3.

In the text

Fig. 9 1D case of thermal inversion (10 bar, N2+40% CH4+0.4%CO2). Left: opacity of the atmosphere with the zone of emission of the star and the planet. Middle: temperature profile and the layers of emission of the atmosphere at 15 µm. Right: eclipse depth. The methane absorbs partly in the star-emission region, creating the thermal inversion in the temperature profile. The CO2 absorbs in the 15 µm band, and thus the emission layers at 15 µm are high enough in the atmosphere, where the temperature is high. Thus, the eclipse depth is higher.

In the text

	Fig. 10 Thermal inversion: 1 bar, N2+40% CH4+0.4%CO2 . The panels are the same as in Fig. 8.
In the text

	Fig. 11 Thermal inversion: CO2+ simplified hazes. The panels are the same as in Fig. 8. Model 1 and Model 2 for the hazes are described in the text. The temperatures and winds of the lower panels correspond to the Model 1 case, which is the model closest to the data.
In the text

	Fig. 12 Thermal inversion: CO2+aerosols with realistic optical properties. The panels are the same as in Fig. 8. The parameters for dust and tholins are described in the text. The temperatures and winds of the lower panels correspond to the dust case, which is the model closest to the data. The difference between the 1D and 3D thermal profiles is discussed in Appendix C.2.
In the text

	Fig. 13 Phase curves in the 12.8 µm MIRI filter for an airless body and the remaining atmospheric families, along with the corresponding simulated observational points.
In the text

	Fig. 14 Dayside emission for an airless body and the remaining atmospheric families. The red points are from Ducrot et al. (2025). The purple points are predicted observation points for five visits for an airless body of albedo 0.2.
In the text

Fig. A.1 Upper panels: Observables of the (“transparent” atmospheres. The CO2 partial pressure is 1.e-7 bar, except for the (“transparent” case where the observables have been calculated for PCO2 =1.e-15 bar. For the “no CIA” case, the collision-induced absorptions have not been included in the opacity computation. Lower panels: Temperature profile and map for the 10 bar case. The atmosphere is very well mixed, but there is no efficient coupling with the surface, and the surface has a hot dayside and cold nightside.

In the text

Fig. A.2 Pure CO2 case. Upper panels: Observables (right: phase curve at 15 µm, left: eclipse depth). The dotted lines represent non-stable atmospheres; these are snapshots of the atmosphere during its collapse. Lower panels: Physical quantities. Temperature profiles and temperature map for the 10 Pa case (compatible with the observations, but not stable). In the lower right panel, the red contour shows the CO2 ice caps. The red dot is the substellar point. In the lower left panel, the red line corresponds to the CO2 condensation curve.

In the text

	Fig. A.3 Effect of the surface albedo for different cases of N2-CO2 atmospheres (results of 1D simulations).
In the text

	Fig. B.1 1D study of thick atmosphere with greenhouse effect with an opacity window at the observed wavelengths. Left: Idealized opacity vs wavelength. The yellow window corresponds to the zone of emission of the star, the red window to the planet. Middle: Temperature profile. The window corresponds to the layers of 95% of the emission. Right: Corresponding eclipse depth.
In the text

	Fig. B.2 Opacities of molecules. The dotted lines show the MIRI filters where the molecules must not absorb too much in order to produce the right effect. The planet blackbody emission range is in red.
In the text

	Fig. B.3 Thermal profiles for different simulations of a 10-bar N2 atmosphere with 20% CH4 and 0.2% C2H4.
In the text

	Fig. C.1 Temperature profile and eclipse depth for different parameterizations of aerosols. Each panel illustrates the influence of a parameter.
In the text

	Fig. C.2 Single scattering albedo vs wavelength for different aerosol types. The radii of aerosols are 10 × 10 6 for tholins and H2SO4 and 10 × 10 4 for dust.
In the text

	Fig. C.3 Thermal 1D profile for the dust case, converged vs not converged. The initial profile is a 370 K isotherm.
In the text

	Fig. D.1 Transit spectra computed for scenarios that match all the available emission observations. The transit spectrum has been offset to have a minimum value of 0.
In the text