1D atmospheric study of the temperate sub-Neptune K2-18b

The atmospheric composition of exoplanets with masses between 2 and 10 M$_\oplus$ is poorly understood. In that regard, the sub-Neptune K2-18b, which is subject to Earth-like stellar irradiation, offers a valuable opportunity for the characterisation of such atmospheres. Previous analyses of its transmission spectrum from the Kepler, Hubble (HST), and Spitzer space telescopes data using both retrieval algorithms and forward-modelling suggest the presence of H$_2$O and an H$_2$--He atmosphere, but have not detected other gases, such as CH$_4$. We present simulations of the atmosphere of K2-18 b using Exo-REM, our self-consistent 1D radiative-equilibrium model, using a large grid of atmospheric parameters to infer constraints on its chemical composition. We show that our simulations favour atmospheric metallicities between 40 and 500 times solar and indicate, in some cases, the formation of H$_2$O-ice clouds, but not liquid H$_2$O clouds. We also confirm the findings of our previous study, which showed that CH$_4$ absorption features nominally dominate the transmission spectrum in the HST spectral range. We compare our results with results from retrieval algorithms and find that the H$_2$O-dominated spectrum interpretation is either due to the omission of CH$_4$ absorptions or a strong overfitting of the data. Finally, we investigated different scenarios that would allow for a CH$_4$-depleted atmosphere. We were able to fit the data to those scenarios, finding, however, that it is very unlikely for K2-18b to have a high internal temperature. A low C/O ratio ($\approx$ 0.01--0.1) allows for H$_2$O to dominate the transmission spectrum and can fit the data but so far, this set-up lacks a physical explanation. Simulations with a C/O ratio $<$ 0.01 are not able to fit the data satisfactorily.


Introduction
While fairly common among the thousands of exoplanets discovered to date 1 (see Figure 1), the atmospheres of super-Earths and sub-Neptunes -with masses between 2 and 10 M ⊕ -are poorly understood. By essentially using the mass and radius of the planets, it has been well-established that planets with low masses (< 2 M ⊕ ) must primarily be made up of iron and silicates, and generally with a thin atmosphere. On the other hand, planets with high masses (> 10 M ⊕ ) must retain a thick atmosphere composed mainly of H 2 and He, representing a significant portion of the planet mass (Chen & Kipping 2016). Within the transition between the two populations, however, studies with models struggle to give a clear answer (Valencia et al. 2013;Fulton et al. 2017;Zeng et al. 2019;Otegi et al. 2020), exhibiting a degeneracy between massive rocky planets, ocean planets, and small gaseous planets. Thus, spectral observations are crucial for characterising the atmosphere of these objects and to better constrain their internal composition.
The recently discovered, transiting exoplanet K2-18b (see Table 1) offers not only the opportunity to retrieve spectroscopic data on the atmosphere of a sub-Neptune, but also to study such atmospheres under nearly Earth-like conditions. Indeed, the stellar irradiance received by K2-18b (E p, e = 1368 +114 −107 W·m −2 ) is very close to that of the Earth. Nine transits were acquired us-1 exoplanet.eu ing the Wide Field Camera 3 on the Hubble Space Telescope (HST/WFC3), as well as data from Kepler K2 and Spitzer IRAC channels 1 and 2 (Benneke et al. 2019a). These data have already been analysed by several teams (Benneke et al. 2019a;Tsiaras et al. 2019;Madhusudhan et al. 2020;Scheucher et al. 2020;Bézard et al. 2020). Also, Scheucher et al. (2020) and Bézard et al. (2020) used self-consistent models (respectively, '1D-TERRA' and 'Exo-REM') to analyse the data , while the other used free retrieval algorithms. The first five investigations concluded that there is a presence of H 2 O as well as a significant amount of H 2 -He. They also derived upper limits for the abundance of CH 4 , which in all cases was lower than ≈ 3.5% at a 99% confidence level. In contrast, Bézard et al. (2020) found that the HST/WFC3 spectrum is dominated by CH 4 absorption and found abundances of CH 4 and H 2 O, respectively, between 3% and 10% and 5% and 11% at a 1-σ confidence level. They assumed an H 2 -dominated atmosphere and varied the atmospheric metallicity, but their simulations did not include H 2 O clouds and they did not simulate non-solar C/O ratios.
Given the irradiance of the planet and the presence of H 2 O in the atmosphere, the question of the existence of liquid H 2 O is naturally posed. Benneke et al. (2019a) found that a cloud layer was needed to reproduce the data. They then used a selfconsistent model and found that liquid H 2 O could condense at the right pressure to explain this cloud layer. Contrary to Benneke et al. (2019a), Madhusudhan et al. (2020) did not find com-  (1.005 E ⊕, e ) 2 Notes. ( * ) We actually use a 10 5 -Pa radius of 16400 km in our simulations. See text. The radius from (Benneke et al. 2019a) corresponds to a pressure of 10 3 Pa. pelling evidence for clouds or hazes in the atmosphere. Using an interior model, Madhusudhan et al. (2020) indicated that if the planet had a small rocky core and a thin H 2 /He atmosphere, an ocean of liquid H 2 O could exist. However, Scheucher et al. (2020) ruled out this possibility, arguing that an H 2 O ocean would partially evaporate in the atmosphere, giving a spectrum that would be incompatible with the data.
In the present work, we use Exo-REM, our self-consistent one-dimensional (1D) atmospheric model adapted for transiting exoplanets to study K2-18b atmospheric composition assuming a thick H 2 -He atmosphere. We first present the extension of Exo-REM (Baudino et al. 2015(Baudino et al. , 2017Charnay et al. 2018) to irradiated planets and detail the calculations of the absorption crosssections. We then expand on the work of Bézard et al. (2020) by including clouds and studying their formation, and by investigating the effects of irradiation, the eddy diffusion coefficient, internal temperature, and metallicity on the atmosphere properties.

Model
Exo-REM is a 1D radiative-equilibrium model first developed for the simulation of young gas giants far from their star and brown dwarfs (Baudino et al. 2015(Baudino et al. , 2017Charnay et al. 2018). Fluxes are calculated using the two-stream approximation assuming hemispheric closure. The radiative-convective equilibrium is solved assuming that the net flux (radiative + convective) is conservative. The conservation of flux over the pressure grid is solved iteratively using a constrained linear inversion method. We take into account Rayleigh scattering from H 2 , He, and H 2 O, as well as absorption and scattering by clouds -calculated from the extinction coefficient, single scattering albedo, and asymmetry factor interpolated from pre-computed tables for a set of wavelengths and particle radii (see Charnay et al. 2018).
We made several upgrades to this model in order to extend its simulation capabilities to high-metallicity (up to ≈ 1000 times the solar metallicity) and irradiated planets, namely sub-Neptunes and moderately hot Jupiters (equilibrium temperature  ). The scatter plot shows only the 454 planets for which the mass and the radius are known within ±15%. The mass-radius curves for an Earth-like interior (orange) and a pure water ice planet (blue) are taken from Seager et al. (2007). lower than 2000 K). These upgrades are detailed in the following sections. The Exo-REM source code and the k-coefficients we used (see Section 2.1.3) are available online 2 .

Collision-induced absorptions
The collision-induced absorptions (CIA) included in our simulations and their references are given in Table 2. We added the contributions of the H 2 O-H 2 O CIA and, since spectroscopic data for the H 2 O-H 2 CIA were not available, of the H 2 O-air CIA. The reasons for this addition are detailed in Section A.

Absorption cross-sections
We calculated the absorption cross-sections of CH 4 , CO, CO 2 , FeH, H 2 O, H 2 S, HCN, K, Na, NH 3 , PH 3 , TiO, and VO at 25 pressure levels equally spaced in the log-space between 0.1 and 10 7 Pa. At each of these pressure levels, we calculated the crosssections at temperatures 100,150,200,250,300,400,500,600,800,1000,1200,1500,2000,2500, and 3000 K -except for NH 3 , the reasons for which are detailed below. We calculated line absorption up to a given distance (∆ν) from the line center, using the same procedure as described in Baudino et al. (2015). We used a sub-Lorentzian line profile with a χ factor, based on Burch et al. (1969) and Hartmann et al. (2002), for, respectively, CO 2 and then all the other species. The cross-sections were calculated over the wavenumber range displayed in Table 3, at a resolution that is similar to the line half width.
When line lists of individual isotopes were available, we merged them by multiplying line intensities in order to reproduce the isotopic ratio found for Jupiter. Otherwise, we used the default isotopic ratio given by the database. We used H 2 and He pressure-broadened halfwidths (γ) and temperature exponents (n) whenever they were available. When both were available, we used a mixture of 90% H 2 and 10% He, which roughly corresponds to the standard Solar System He/H ratio (Lodders 2019). More details are given in Table 3 and Table 4. The specificities of some species are listed below.
CH 4 : We included the contribution from 12 CH 4 , 13 CH 4 , and CH 3 D, all taken from the TheoReTS database (Rey et al. 2017), which is more accurate than the ExoMol database. We used a 12 C/ 13 C ratio of 89 (Niemann et al. 1998) and a D/H ratio of 2×10 −5 (Lellouch et al. 2001). We recall that the TheoReTS line list of CH 3 D stops at 6500 cm −1 . The line list provided by the TheoReTS database separates the 'strong lines', which have an intensity of > 10 −26 cm·molecule −1 , from the very-high-spectraldensity 'weak lines', which are regrouped in 'super lines'. This separation allows for much faster cross-section calculations, but at the price of losing information on the quantum numbers of individual lines. Hence, we had to use the same γ and n for all lines.
FeH, TiO and VO: H 2 or He γ and n were not available for these molecules and we used the parameters from CO.
K and Na : Given that these species have very high intensity lines (up to ≈ 10 −13 cm·molecule −1 ), we needed to extend our ∆ν up to 9000 cm −1 in order to correctly account for far wing absorption. Following Burrows et al. (2000), we used a Voigt profile up to a detuning frequency of 20 × (T/500) 0.6 cm −1 for K and 30 × (T/500) 0.6 cm −1 for Na. Beyond that detuning frequency, we used the profile described by Baudino et al. (2015).
NH 3 : Cross-sections were calculated up to 1500 K because the line list we used lacks completeness above this temperature. We took a 14 N/ 15 N ratio of 500 (Furi & Marty 2015).
PH 3 : The line list we used is not complete above 1000 K. However, Sousa-Silva et al. (2014) provide the PH 3 partition function up to 3000 K. This allows us to use the percentual loss of completeness to estimate the proportion of missing opacity, as suggested in the cited work.

k -coefficients
We used these high-resolution absorption coefficients to calculate k-coefficients according to the method described in Baudino et al. (2015), using 16 Gauss-Legendre quadrature points (eight for values of the cumulative distribution of the absorption coefficients between 0 and 0.95, and eight between 0.95 and 1). Four sets of k-coefficients were calculated, at resolutions of 0.5, 20, 200, and 2000 cm −1 , although only the 20 cm −1 step one was used in the main part of this study.

Radiative-convective equilibrium model: stellar irradiance
We added the planetary averaged stellar irradiance E ↓, e,ν reaching the planet to Exo-REM. We used stellar spectra from the BT-Settl model (Allard et al. 2012). We chose a spectrum modelled at an effective temperature T S = 3500 K, log 10 (g[cm · s −2 ]) = 5, a null metallicity and no alpha enhancement. We neglected the spectral dependency on surface gravity and metallicity. Then we interpolated the BT-Settl spectrum on a spectral grid with a 0.1 cm −1 wavenumber step. This interpolated spectrum was then convolved at the resolution of Exo-REM flux calculations (see Section 3) to obtain the modelled radiosity J S, e,ν . Then we obtained E ↓, e,ν from:   (11) References.   (VMR) profiles of our absorbing species (see Section 2.1.2), are susceptible to the formation of a major cloud layer, or substantially impact the molar mass of the atmosphere (according to Lodders 2010).
Most of the equilibrium abundances are derived from the equation of conservation of each species, using the standard Gibbs free energy of formation ∆G • f listed in Chase (1998). For PH 3 (g), we used the ∆G • f from Lodders (1999), while for MnS(cr,l) and ZnS(cr,l), we used the values from Robie & Hemingway (1995), and for CaTiO 3 (cr) we used values from Wood-field et al. (1999). We used a grid of temperatures between 200 and 4000 K, with a 100 K step, to map our ∆G • f . When the temperature range of the references used were smaller than our target temperature range, we used a linear extrapolation to fill up our grid. The ∆G • f at a given temperature are linearly interpolated from this grid. We calculated the saturation pressure of the following species directly: H 2 O (Wagner & Pruss 1993;Wagner et al. 1994;Lin et al. 2004;Fray & Schmitt 2009), NH 3 (Fray Lodders (2002), but the impact on our simulations was expected to be limited. We simulated the formation of Fe-Ni alloys by treating the Fe independently while keeping the Ni(g)/Fe(g) ratio constant. We considered the formation of H 3 PO 4 instead of P 4 O 6 , based on the results of Wang et al. (2016).
We included non-equilibrium processes in the CH 4 -CO, CO-CO 2 , N 2 -NH 3 -HCN chemical systems following the approach of Zahnle & Marley (2014). The chemical quenching level is determined by equating a reaction timescale to the mixing time H 2 /K zz , where H is the atmospheric scale height and K zz the eddy mixing coefficient. Below the quenching level, the abundance profiles of the relevant species are governed by thermochemical equilibrium while, above this level, the mixing ratios of the quenched species are held constant. We used the chemical reaction timescales given by Eqs. 12-14 for CH 4 -CO, Eq. 44 for CO-CO 2 , Eq. 32 for NH 3 -N 2 , and Eq. 40 for HCN-NH 3 -N 2 in Zahnle & Marley (2014). In the case of PH 3 , we assume that its conversion to H 3 PO 4 is inhibited, as observed in the giant planets of our solar system. We compared our chemical model in the case of K2-18b (taking identical temperature profile, elemental abundances, and eddy diffusion coefficient) against the model described in Venot et al. (2020) 4 . For CO, HCN, and NH 3 , we found VMR differences lower than 20% with our model in the HST sensitivity region, and lower than 2% in the case of CH 4 and H 2 O. The difference is larger for CO 2 , with our VMR being 25% higher than the value found in O. Venot's model. We consider these differences to be satisfactory.

Methodology
The star and planetary parameters used in our simulations are displayed in Table 1. A summary of our model grid is displayed in Table 6. We took a planetary radius at 10 5 Pa of R p = 16400 km, which is slightly different from the value used by Benneke et al. (2019a) (who define R p as the radius at 10 3 Pa), because it is closer to the radius we find when fitting the observed data. To compare our model spectrum with the data, we applied an offset on R p in the calculation of the transit depth, such that the χ 2 is minimised.
The internal temperature resulting from the residual heat of formation of the planet was calculated following (Rogers & Seager 2010): where L int is obtained from (Rogers & Seager 2010): where a 1 = −12.46±0.05, a M p = 1.74±0.03, a R p = −0.94±0.09, a t p = −1.04 ± 0.04, t p is the age of the planet, and the astronomical constants are defined by the International Astronomical Union (IAU, Mamajek et al. 2015). According to a gyrochronological model from Guinan & Engle (2019), the age of K2-18 can be estimated at 2.4 ± 0.6 Gyr. Assuming that K2-18b ended its formation a few Myr after the formation of its star (so that t * ≈ t p ), like what happened in our solar system, we obtained from Equation 2 T p, int ≈ 83 +8 −6 K. Rounding down this value, we chose 80 K as our nominal internal temperature.
The net fluxes are calculated from 40 to 30000 cm −1 , with a step of 20 cm −1 . The atmospheric grid consists of 81 levels equally spaced in the log-space between 0.1 and 10 7 Pa. We imposed a correlation length of 0.5 pressure scale height to the solution temperature profile in order to avoid non physical oscillations (see Baudino et al. 2015). We considered cloud radiative effects and scattering only from H 2 O clouds. The NH 3 and NH 4 SH clouds would form in atmospheres that are colder than what is expected for K2-18b, while NH 4 Cl clouds are too thin to have a significant impact, and other clouds are condensing too deeply into the atmosphere (see Section 4.2). Cloud vertical mass mixing ratios are calculated assuming equilibrium between the downward flux of falling particles and the upward flux of gas and cloud particles due to advection and turbulent mixing (see Charnay et al. 2018).
When we change the irradiance of the planet by a factor of k, we do so by multiplying a p by 1/ √ k (see Equation 1). We proceed this way to make model comparison easier. Indeed, a p impacts only the irradiance, while the two other parameters, T * , eff and R * , affects respectively the stellar spectrum and both the stellar spectrum and the transmission spectrum.
We estimated K zz using the values derived by Exo-REM (see Charnay et al. 2018). We found that typical K zz values for K2-18b range from 10 6 to 10 9 cm 2 ·s −1 , with the highest values found in the convective layers. Given the uncertainties on the estimation of K zz , we enlarged this range to 10 5 -10 10 cm 2 ·s −1 . We also found that a quenching of our species often occurs just below the uppermost convective layer (see Section 4.3), hence, we set our nominal K zz value at 10 6 cm 2 ·s −1 .
We used a fractional area covered by clouds f c = 0.15 for the calculation of the temperature profile, and f c = 1.0 for the calculation of the transmission spectra. To ensure numerical stability, the cloud mean particle radius, r, was fixed at a constant value (50 µm) that roughly corresponds to the maximum one predicted by the Exo-REM self-consistent cloud model at 300 (Z/H) and nominal irradiation (see Charnay et al. 2018, and Section 4.2), unless stated otherwise. We call this our nominal model. We explore the effects of changing some of these parameters in the following sections.

Metallicity and irradiance
Here, we assume that K2-18b has retained a relatively thick H 2 -He atmosphere. This atmosphere can be enriched in heavy elements -compared to the initial proto-stellar nebula -during the formation of the planet via collisions with planetesimals (Fortney et al. 2013). There could also be a contribution from the erosion of the primordial core (Iaroslavitz & Podolak 2007).
In Figure 2, we show our best fit for the dataset of Benneke et al. (2019a) at nominal irradiation, along with the HST data reduced by Tsiaras et al. (2019). The contributions of the different opacity sources to this spectrum are represented in Figure 3. The temperature profile and VMR of this spectrum are represented respectively in Figure 4 and Figure 5. Because the stellar spectrum we used is slightly more intense than a black body between ≈ 1 and 2.5 µm, the temperatures we obtain in the upper atmosphere are slightly higher (≈ 5 K at 1 kPa) than what would be  0.5-1.5 1 K zz (cm 2 ·s −1 ) 10 5 -10 10 10 6 r (µm) 20-600 50 T p, int (K) 45-200 80 obtained with a black body at the effective temperature of the star. The prevalence of CH 4 absorptions over H 2 O absorptions is discussed in Section 5.1.
In Figure 6, we display the χ 2 of our nominal models for the datasets from Benneke et al. (2019a) and Tsiaras et al. (2019), including their respective reduction of the same HST raw data, as well as in both cases the K2 and Spitzer data from Benneke et al. (2019a). We performed simulations of K2-18b for 1, 3,10,30,50,75,100,125,150,175,200,300,400,500, and 1000 times (Z/H) and for irradiations between 0.5 and 1.5 time the nominal irradiation, with a step of 0.1. We consider a model as statistically accurate if it can reproduce the data within the 1σ confidence level (68%). Since there are 20 data points, we have 20 to 1 degrees of freedom, so the 1σ confidence level corresponds to a χ 2 of 21.36. While the dataset from Tsiaras et al. (2019) allows, according to our interpolation, for metallicities ≥ 65 (Z/H) within K2-18b irradiation 3σ uncertainties, the dataset from Benneke et al. (2019a) is much more restrictive, allowing only a metallicity between 100 and 200 (Z/H) at 0.9 times the nominal irradiation, or a metallicity between 150 and 200 (Z/H) at or above the nominal irradiation. These values are within the most common range of metallicity predicted by Fortney et al. (2013) for planets with radius in the 2-4 R ⊕ range (between 100 and 400+ (Z/H) ). At nominal irradiation, our best fit against Benneke et al. (2019a) and Tsiaras et al. (2019) data is located respectively at 175 (Z/H) (χ 2 = 21.07) and 150 (Z/H) (χ 2 = 17.24). In both cases, there is no H 2 O cloud formation, and we derive a Bond albedo of respectively 0.017 and 0.018.
The relationship between temperature, metallicity, and goodness of fit presented in Figure 6 can be explained as follows. The amplitudes of the features in the transmission spectra are correlated with the abundance of absorbers in the atmosphere as well Article number, page 6 of 22 D. Blain et al.: 1D atmospheric study of the temperate sub-Neptune K2-18b as to its scale height, which can be written as where R is the gas constant, and T , µ and, g are respectively the temperature, the molar mass and the gravity in an atmospheric layer. As the metallicity decreases, two concurrent effects are occurring. The diminution of the VMR of absorbers in the atmosphere will of course decrease the spectrum amplitude. However, decreasing the VMR of heavy species will also leads to a decrease of µ, and thus to an increase of the amplitudes. This latter effect is dominant in our simulations for metallicities ≥ 10 (Z/H) , as illustrated in Figure 7. At lower metallicities, heavy elements contribute to less than 10% of µ, so the effect of less intense gas absorptions becomes dominant, flattening the transmission spectrum. If we decrease irradiation, the atmosphere gets colder and H 2 O clouds start to form and to thicken, flattening the spectrum at transparent wavelengths where part of the transmitted stellar light reaches the condensation level. As stellar irradiation decreases, the cloud forms lower and lower in the atmosphere, eventually ending below the region probed by transit spectroscopy where it can no longer directly affect the spectrum.
The condensation also removes more gaseous H 2 O, decreasing µ and H 2 O absorption. On the other hand, increasing the irradiation leads to an increase of the temperatures, and, hence, of the scale height and of the absorption amplitudes, as shown in Figure 8.

Clouds
In Figure 9, we show that H 2 O cloud can form at nominal irradiation and metallicities above 200 (Z/H) . At 300 (Z/H) , the cloud forms at 1 kPa. However, this cloud cannot form even at 1000 (Z/H) if the irradiation is more than 10% higher than the nominal irradiation. As a general tendency, the H 2 O cloud forms lower in the atmosphere as the irradiation decreases and the metallicity increases. This is because as the irradiation decreases, the temperature decreases, so the temperature profile crosses the H 2 O condensation profile at higher pressures. Also, as the metallicity increases, the partial pressure of H 2 O increases, shifting its condensation profile towards higher temperatures.
The effect of H 2 O clouds at nominal irradiation on the temperature profile is insignificant, but its effect on the transmission spectrum is quite important, as shown in Figure 10. The cloud shields the transit spectrum from the star light passing below its altitude, with a particularly strong effect in the visible range. The cloud also has a slight effect on the Bond albedo, as shown in Figure 11. Under 200 (Z/H) , increasing the metallicity decreases the Bond albedo, because there is more and more gas absorption. At 300 (Z/H) , adding clouds increases the Bond albedo from 0.017 to 0.018. This effect then increases with metallicity, because the cloud thickens and, thus, it is able to reflect more light. We note that above 300 (Z/H) , the Bond albedo increases with metallicity even without clouds. This is due to the removal of gaseous H 2 O from the atmosphere.
We were unable to make Exo-REM simulations converge at nominal irradiation when using its self-consistent cloud mode (Charnay et al. 2018). In this mode, the mean particle radius is determined from a comparison of the timescales of the microphysical processes governing the formation and growth of cloud particles. This is likely due to the temperature profile barely crossing the H 2 O saturation profile, so that a slight variation of temperature in the solution profile can have a major effect on cloud formation. Nevertheless, we were able to determine that the H 2 O cloud particles, according to Exo-REM, should be between 20 and 50 µm at 300 (Z/H) , for a calculated K zz of ≈ 9 × 10 7 cm 2 ·s −1 at the layer of condensation. The result is essentially the same if we fix the sedimentation parameter ( f sed ) between 1 and 5, which is typical of the clouds in the solar system (see Charnay et al. 2018). If we fix K zz at 10 10 cm 2 ·s −1 and f sed at 5, we obtain a maximum mean cloud particle radius of 600 µm.
In Figure 12, we show the effect of the mean cloud particle radius on the transmission spectrum. As the mean radius increases, the effect of the cloud is less and less visible. At 600 µm, it is almost indistinguishable from the spectrum without clouds.
We found that in our simulations, liquid H 2 O clouds can form at irradiations lower than 80% of the nominal irradiation, which is slightly lower than the 3σ lower uncertainty of K2-18b irradiation. However, according to Figure 6, this case is not favoured by the data, and is above the 1σ confidence level against Benneke et al. (2019a) data. It is not surprising that our results for the likelihood of liquid H 2 O differ from Benneke et al. (2019a), even though they also used a self-consistent model in their demonstration. Indeed, they assumed an albedo of 0.3, which is much higher than what we found. They probably also included a much lower amount of CH 4 , which significantly reduces the stellar heating.
We note from Figure 4 that aside from H 2 O, clouds of NH 4 Cl, KCl, ZnS and Na 2 S are condensing within our pressure grid. Clouds below KCl (included) form too low in the atmosphere to significantly impact the temperature profile above the convective layers, so they can safely be ignored. This is not the case for NH 4 Cl. According to our simulations, NH 4 Cl removes ≈ 30% of the total amount of Cl in the atmosphere, with the remaining 70% being removed mainly via the condensation of KCl, RbCl, and CsCl. The latter two are not included in Exo-REM, the condensation of RbCl and CsCl removing less than 0.15% of the Cl, assuming a standard composition. While the NH 4 Cl cloud is, to our knowledge, rarely mentioned in the exoplanet literature, it offers an explanation to the lack of Cl-bearing species in the upper atmosphere of the giant planets of our solar system (Fouchet et al. 2004;Teanby et al. 2006). Nevertheless, we found that with our nominal model, NH 4 Cl clouds with r = 5 µm -roughly the minimum value found by the Exo-REM self-consistent cloud model on the layer with the most particles -have essentially no impact on the temperature profile and a marginal impact on the visual part of the transmission spectrum 5 , lowering the χ 2 against Benneke et al. (2019a) data by ≈ 0.2. Therefore, they can probably be neglected for K2-18b, but they might play a more important role in the transmission spectra of slightly hotter planets.

Eddy diffusion coefficient and internal temperature
In Figure 13, we represented the temperature profiles and VMRs obtained with our upper and lower uncertainty boundaries on the internal temperature of K2-18b (assuming t p between 1 and 13 Gyr), as well as our selected nominal value (80 K). It appears that the internal temperature has only a minor effect on the temperature profile, except in the deepest layers of the atmosphere.  (2017). This is because the received flux (the stellar irradiance) is two to three orders of magnitude higher than the internal flux. Hence, above the convective layers (here at ≈ 20 kPa), the difference of temperature between our two extreme simulations is less than 40 K (8 K on average). However, below the convective layers (≈ 200 kPa), the atmosphere is opaque, thus, the stellar irradiance can no longer heat the atmosphere and the internal heating become prominent. Consequently, at the bottom of our pressure grid, the difference in temperature between our two extreme models reaches ≈ 500 K. This affects the species abundances, in the upper atmosphere as well if the quench level of a species is at a pressure level below the region of stellar heating.
In Figure 14, we show the effect of the eddy diffusion coefficient on the VMR of the main C-, N-, and O-bearing species at 175 (Z/H) . The VMR at thermochemical equilibrium are also displayed. The main effect of increasing K zz is to shift down the quench levels of the species.
In Figure 15, we display the VMR of CH 4 , CO, and H 2 O for K zz = 10 6 , 10 8 , and 10 10 cm 2 ·s −1 as a function of the age of the planet and the corresponding internal temperature (see Equation 2). As the internal temperature decreases (or as the planet gets older), the impact of the K zz on these VMR gets smaller. If the K zz is low ( 10 6 cm 2 ·s −1 ), the chemical composition of the atmosphere essentially does not change. On the other hand, under vigorous mixing, the chemical composition of the upper atmosphere can be significantly altered over time. With K zz 10 10 cm 2 ·s −1 , if the planet is young ( 1 Gyr or T p, int 100 K), CO could become the dominant O-bearing species and remain more abundant than CH 4 for several Gyr. In parallel, K zz could also decrease over time, as the heat to be dissipated by convection decreases. Moreover, if the quench level switches from a convective layer to a radiative one, K zz could drop by several orders of magnitude, rapidly changing the chemistry of the upper atmosphere if this occurs when the planet is still young.
We compared our simulations within our range of K zz , at T p, int of 45, 80, and 115 K and at metallicities of 50,75,100,125,150,175,200, and 300 (Z/H) . The best-fit parameters against Benneke et al. (2019a) data for each tested internal temperature can be found in Table 7. We found no significant differences at constant metallicity in terms of goodness of fit for the tested T p, int , the χ 2 value varying by less than 0.7 across all K zz values. In each case, our best fit is located at 175 (Z/H) . The slight difference in goodness of fit at T p, int = 45 K compared to the other tested internal temperatures is due to H 2 O cloud condensation occurring in the former case and not in the latter. χ 2 values above 200 and below 150 (Z/H) are systematically above the 1σ confidence level. We note that even at 175 (Z/H) , T p, int = 115 K and K zz = 10 10 cm 2 ·s −1 , we retrieve a CH 4 VMR above 1.5% in the upper atmosphere. Hence, internal heating from a relatively young (t p ≈ 1 Gyr) K2-18b and Wavelength (m)

Considering whether H 2 O or CH 4 is the dominant absorber in the HST transit spectra
Our nominal model, presented in Section 4, shows that CH 4 should be the main C-bearing species in the atmosphere, and, despite H 2 O being more abundant than CH 4 , the latter should be the main contributor to the absorption features of the transmission spectrum. In particular, CH 4 and H 2 O have an overlapping band at 1.4 µm in the HST-WFC3 spectral range. This is outlined in Bézard et al. (2020), who showed that at this wavelength and for sub-Neptunes with effective temperatures lower than ≈ 600 K and high atmospheric metallicities, CH 4 should be a stronger absorber than H 2 O in transit spectroscopy due to the numerous  To investigate this discrepancy, Bézard et al. (2020) compared their best fit Exo-REM spectrum (in which CH 4 absorption dominates over H 2 O) and an H 2 O-only spectrum (in which all other absorbers were removed) to the HST data reduced by Tsiaras et al. (2019). The χ 2 associated with the two models are essentially the same (0.91 and 0.93). In addition, applying the retrieval algorithm TauREx 3 (Al-Refaie et al. 2019) to these data, Bézard et al. (2020) found solutions favouring a CH 4 -rich atmosphere whose absorption at 1.4 µm is dominated by CH 4 . Thus, the apparent disagreement with the analysis of Tsiaras et al. (2019) simply arises from the fact that these authors did not consider an atmosphere with significant amounts of CH 4 in the three scenarios they investigated.
Regarding the Benneke et al. (2019a) dataset, we used the same technique as Bézard et al. (2020): we took the H 2 O contri-bution of our best-fit spectrum, displayed in Figure 3, and compared it with the data. In the same way, we made another spectrum with only the contributions of the absorbers included in the retrieval analysis of Benneke et al. (2019a) except CH 4 (i.e. H 2 O, CO, CO 2 , and NH 3 , HCN contribution is negligible and was not included). We will refer to this latter spectrum as 'no CH 4 '. We stress that none of these spectra are directly the results of self-consistent simulations. They are displayed in Figure 17. We found that our H 2 O-only spectrum does provide a superior fit compared with our best fit Exo-REM model. The fit is improved in the 1.15-1.20 µm spectral range, around 1.6 µm, and for the Spitzer data points. With our 'no CH 4 ' spectrum, the fit is also better than our best-fit model, but less so than our H 2 Oonly spectrum due to the absorptions of NH 3 in the HST spectral range and the CO and CO 2 absorptions in the Spitzer spectral range. We note that our VMRs of NH 3 , N 2 , CO, and CO 2 of respectively 0.04%, 1.02%, 2.73%, and 0.96% (see Figure 5), are all lower than the '2σ (97.5%)' upper limits given by Benneke et al. (2019a) for these species (respectively 13.5%, 10.9%, 7.45% and 2.4%). The spectrum with only H 2 O absorption is strongly overfitting the Benneke et al. (2019a) data, with a χ 2 of 13.01. On the other hand, our best-fit model, with a χ 2 of 21.07, is still within the 1-σ confidence interval (χ 2 < 21.36). The 'no CH 4 ' spectrum is characterised by a slight overfitting, providing a better χ 2 (18.83) compared to our best fit, but the removal of CH 4 absorption is artificial.
We also applied TauREx 3 to the Benneke et al. (2019a) data. This algorithm uses the nested sampling code Multinest (Feroz et al. 2009) to explore the parameter space and find the best fit corresponding to a given spectrum. In this retrieval analysis, we used 500 live points, an evidence tolerance of 0.5, and the crosssections provided by the TauREx website 7 at a resolution power of 15000 and between 0.3 and 15 µm. We simulated the atmosphere of K2-18b using an isothermal temperature profile with 100 atmospheric layers between 10 −3 and 10 6 Pa. We took into account molecular absorptions from CH 4 , CO, CO 2 , H 2 O, and NH 3 , CIA from H 2 -H 2 and H 2 -He, the contribution to µ of N 2 , Rayleigh scattering and spectrally gray clouds. The planet radius was allowed to vary by ±10% of its value determined by Benneke et al. (2019a), while the cloud top was allowed to be between 10 6 and 1 Pa. We chose to set bounds on the temperature and the species VMR, respectively, between 200 and 400 K, and between 10 −10 and 0.   butions are shown in Figure 18. We also tested models removing all molecular absorptions except H 2 O or CH 4 , removing clouds, or including all absorptions except that of H 2 O. For comparison, we also computed a 'no active gas' model including only the cloud contribution, that is, a flat spectrum. Finally, we calculated the χ 2 and corresponding σ confidence of each model using the same technique as in Section 4 and the 'binned' spectrum output from TauREx 3. Our results are summarised in Table 8. A small discussion on the results is available in Appendix B.
We found that solutions with low amounts of CH 4 are unambiguously favoured: the 'all absorbers with clouds' model we used is favoured over the 'no H 2 O with clouds' model at 273:1 (3.79σ, see Benneke & Seager 2013). The difference in logevidence between our CH 4 -only models and our 'H 2 O-included' models (i.e. all models including at least the molecular absorption of H 2 O) is, at worst of 4.16, indicating that the latter are favoured at ≥ 64:1 (3.35σ). Moreover, the 2σ (95.4%) upper limit of CH 4 is 0.009% with our 'all absorbers with clouds' model, even lower than what was inferred by Benneke et al. (2019a) (0.248%). This confirms that an H 2 O-dominated spectrum is a far better fit to the Benneke et al. (2019a) dataset than a CH 4 -dominated spectrum, but not that the latter must be rejected. Indeed, looking at the χ 2 , we find that all of the spectra obtained from our TauREx 3 models that do not include H 2 O are slightly below the 1σ confidence level, with χ 2 values and CH 4 VMR close to our Exo-REM best fit. In contrast, all the H 2 O-included TauREx 3 spectra strongly overfit the spectrum, with χ 2 ≤ 10.50, which is much less than the number of data points, and likely much less than the number of effective degrees of freedom. We interpret this as an H 2 O-only spectrum leading to a significant overfit of the dataset.
The overfit that we identified could be explained by an overestimation of the HST error bars. However, the Benneke et al. (2019a) and Tsiaras et al. (2019) HST spectra that we analysed come from independent reduction methods (Tsiaras et al. 2018) and both show similar transit depth uncertainties. Moreover, the custom HST pipeline used by Tsiaras et al. (2019) is described in Tsiaras et al. (2018) to exhibit 'nearly photon-noise limited' performance, meaning that the spectral uncertainties should be close to the theoretical minimum. It can also be seen, for example, in T p, int = 45 K T p, int = 80 K T p, int = 115 K 10 4 10 3 10 2 10 1 10 0 Volume Mixing Ratio 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Pressure (Pa)  are not compatible at the 1-σ confidence level, which points to additional systematic errors. These are probably linked to uncertainties in some parameters used in the data reduction(orbital parameters, limb-darkening coefficients, or calibration, as suggested in Tsiaras et al. 2018). For all of these reasons, we regard this possibility as unlikely.
To summarise, of the two datasets available to us, only the one from Benneke et al. (2019a) strongly favours a CH 4 -depleted atmosphere. This favouritism is likely due to free retrieval algo-rithms searching for the solution that best fits the data. But if such a solution is strongly overfitting, as is the case for an H 2 Oonly spectrum against the dataset from Benneke et al. (2019a), any other solution that is considered to be more physically acceptable may also appear as statistically unlikely. However, if a strong depletion of CH 4 in the atmosphere of K2-18b is not a straightforward scenario from a chemical standpoint and ends up overfitting the data, it is indisputably a statistically valid one 10 2 10 1 10 0 Volume Mixing Ratio  and, thus, it must be considered. In the next section, we try to find a model that could explain such a depletion.

CH 4 -depleted scenarios for K2-18b
While our standard models, presented in Section 4, are statistically able to reproduce the observed spectrum of K2-18b, here we test some scenarios allowing for CH 4 VMR to be as small as retrieved by all the other teams who analysed the data so far. Indeed, our nominal model gives a CH 4 VMR of ≈ 5 +1 −2 % (for metallicities between 65 and 500 (Z/H) ), while Benneke et al. (2019a), Madhusudhan et al. (2020) and Scheucher et al. (2020) give, respectively, a 2σ upper limit of 0.248%, a 99% upper limit of 3.47%, and an upper limit of 460 ppm. Benneke et al. (2019b) propose three explanations to this depletion: (i) a high internal temperature, either from residual heat of formation or tidal heat-ing, (ii) a low C/O ratio resulting from planetary formation process, and (iii) a catalytic destruction of CH 4 by photolysis. The latter cannot be simulated using Exo-REM, so we will focus on the first two possibilities.

High internal temperature
Following the trend presented in Section 4.3, it is possible to reach a CH 4 VMR that is compatible with what was found by Benneke et al. (2019a) by considering T p, int as high as 200 K. Using this internal temperature, we simulated atmospheres at 30,50,75,100,125,150,175,200,300,400, and 500 (Z/H) and within our selected range of K zz . The goodness of fit for our models compared to the datasets from Benneke et al. (2019a) and Tsiaras et al. (2019) is displayed in Figure 19. We note that H 2 O never dominates the transmission spectrum in the 1.355-1.415 µm window. Our best fit against Tsiaras et al. (2019) data (χ 2 = 16.78) is located at 150 (Z/H) and 10 7 cm 2 ·s −1 , while our best fit against Benneke et al. (2019a) data (χ 2 = 20.70) is located at 150 (Z/H) and 10 10 cm 2 ·s −1 . There is no H 2 O cloud condensation in the latter scenario. The goodness of fit is again slightly better than for our nominal model. The VMR are displayed in Figure 20. In this configuration, we obtain VMR for CH 4 , CO, CO 2 , NH 3 , and H 2 O in the upper atmosphere at levels of 0.283%, 6.589%, 0.280%, 0.006%, and 1.657%, respectively, which are all below or close to the 2σ upper limit found by Benneke et al. (2019a).
Reaching such a high internal temperature, however, would require K2-18b to be a very young planet ( 1 Gyr old according to Equation 2) or a planet experiencing strong tidal heating, or both. According to the internal temperature model we used and the age of the system determined by Guinan & Engle (2019), it seems unlikely that the residual heat of formation could be sufficient to explain a high internal temperature. To assess the possibility of strong tidal heating on K2-18b, we can make a few considerations. The K2-18 system is composed of at least another planet: K2-18 c (Cloutier et al. 2019), orbiting closer to its star than K2-18b. Both planets may have a moderately high eccentricity that is rapidly evolving due to the secular interactions between the two objects (Gomes & Ferraz-Mello 2020). In that configuration, a lot of orbital energy can indeed be dissipated, but only in the innermost planet 8 , that is, K2-18 c. Thus, there should be no tidal heating on K2-18b.

Low C/O
To test the possibility of a low atmospheric C/O ratio on K2-18b, we use the results from the simulations presented in Section 4.4. The VMR of the most abundant absorbers in our best fit against Benneke et al. (2019a) are represented in Figure 21. The amount of CH 4 in the upper atmosphere we obtain in this best-fit simulation is 1.23%, incompatible with the 2σ upper limit found by Benneke et al. (2019a) (0.248%), but compatible with the 99% upper limit found by Madhusudhan et al. (2020) (3.47%). However, simulations at metallicities between 50 and 150 (Z/H) at 5 (C/H) and at 20 (Z/H) and 1 (C/H) are statistically able to reproduce the data while allowing the CH 4 VMR to be close or below 0.248%, with a minimum at 0.101%. There is no simulation within the 1σ confidence level for C/H 1 (C/H) (CH 4 VMR 0.05% at 50 (Z/H) ). Accordingly, C/O must be 0.01 (0.02 (C/O) ) in order to obtain a satisfactory fit of the data. We note that C/H must be 5 (C/H) (or C/O 0.1 (C/O) ) in order for Heavily C-depleted atmospheres (C/O 0.01) do not seem to agree with the data. Lowering the abundance of C in the atmosphere naturally decreases the contribution of CH 4 to the transmission spectrum, and improves the fit of Benneke et al. (2019a) data, as discussed in Section 5.2. A side effect, however, is that because CH 4 has a strong greenhouse effect, the temperatures decrease significantly, favouring the condensation of H 2 O, dampening the absorption features of the latter, and worsening the fit around 1.4 µm. This side effect explains why our results are quite different from the results of other teams who used free retrieval algorithms that do not include any physical constraint on the model atmosphere. Moderately-C-depleted atmospheres (C/O ≈ 0.10), on the contrary, are clearly favoured, but in that case CH 4 still dominates the absorption features of the transmission spectrum. Moreover, the precise phenomenon that could lead to such a depletion has, to our knowledge, yet to be discovered.

Observational perspectives
For K2-18b, acquiring higher resolution spectra would permit to discriminate between the three scenarios we identified. We used our 0.5 cm −1 resolution k-coefficient set to derive high-resolution spectra of K2-18b and show some spectral ranges in Figure 22 that could prove interesting in that regard. We note that NH 3 dominates the 1.45-1.56 µm spectral range for any of the following three scenarios.
Firstly, in the nominal scenario, CH 4 should dominate over most of the infrared (IR) transmission spectrum, with particularly strong lines around 1.6 µm (end of the H band), or in the 3.0-4.0 µm spectral range (L band). CO would dominate the spectrum in the 4.6-5.0 µm spectral range (M band) but not in the 2.3-2.5 µm spectral range (end of the K band). Isolating H 2 O lines would be challenging in the IR, but might be possible around 1.16 µm and 1.37 µm (beginning and end of the J band). However, it should dominate around 823 and 924 nm.
Secondly, in the low C/O scenario, contrary to the nominal scenario, H 2 O should overall dominate the IR spectrum over CH 4 . CH 4 may still dominate the spectrum at the end of the H band (1.6-1.8 µm), and in the K and L bands, depending on the intensity of the C-depletion. However, in contrast to the nominal scenario, it should have a minor contribution in the J band. CO might still dominate the M band, again depending on the extent of the of C-depletion.
Thirdly, in the scenario involving high internal temperature or incomplete chemistry, H 2 O should again dominate most of the IR spectrum. CO should dominate in the M band and possibly even in the K band, in contrast to the low C/O scenario. The CH 4 spectral contribution should be similar to the low C/O scenario.
The James Webb Telescope (JWST) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) mission, with their broader spectral range and precision higher than the data we studied here, will probably be very helpful to discriminate between these scenarios. In complement, we could acquire higher resolution spectra from ground-based instruments -such as the CRyogenic InfraRed Echelle Spectro-graph+ (CRIRES+) or the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) mounted on the Very Large Telescope (VLT) -or by using the incoming generation of large telescopes (Extremely Large Telescope, Thirty Meter Telescope, etc.) in order to detect the absorption of individual lines and, thus, to unambiguously detect species.
As an example, we calculated the signal to noise ratio (S/N) we could obtain for a detection of CH 4 in our nominal case using CRIRES+ at a resolution power of 100 000 in the 1.5-1.8 µm spectral range. This spectral range has the advantage of presenting a lot of deep CH 4 lines (50-100 ppm) in a band for which CRIRES+ has a high sensitivity 9 . We calculated a synthetic transmission spectrum at a resolving power of 100 000 us-  Figure 3, with an offset of the 10 5 -Pa level of +98 km. Cyan: H 2 O, CO, CO 2 , and NH 3 total contribution of the transit spectrum displayed in Figure 3, with an offset of the 10 5 -Pa level of +57 km. Black: K2, HST and Spitzer from Benneke et al. (2019a). The dotted red line is our best-fit spectrum with all absorptions (see Section 4). The χ 2 of these spectra against the data is indicated in parentheses.
restrial absorptions and the photon noise of the star, using the following: where λ is the discretised wavelength, λ min and λ max are respectively the minimum and maximum wavelengths of the considered spectral range, δ CH 4 and δ tot are, respectively, the transit depths that consider only the absorption from CH 4 and the absorption from all absorbers; the bar indicates a smoothing by a boxcar of width 16 resolution elements (also used by Brogi et al. 2016, for their CRIRES observation of HD 189733), N is the number of calculated points per resolution element within the spectral range, and S/N * is the expected CRIRES+ S/N on K2-18 for one hour of observation at R = 100 000. From the H limiting magnitude given in eso.org, we derive a S/N of 2500 for K2-18 at R = 50 000 per spectral dispersion element for a 1-hr integration. At R = 100 000, assuming no loss of the star signal, the S/N should thus be 2500/ √ (2) = 1800. Hence, we can obtain a total CH 4 signal of 1500 ppm and, thus, from Equation 5, we have an optimistic value for S/N CH 4 of ≈ 2.7 for one hour of observation with CRIRES+ at R = 100 000. We note that the transit duration of K2-18b is approximately three hours. A full transit would give S/N CH 4 ≈ 4.7, which is sufficient for detecting CH 4 . However, this estimate is probably optimistic as we did not take into account the telluric absorption and the photon noise from the star.
Observing CO from the ground would, on the other hand, prove to be extremely difficult. Even at R = 100 000, the CO lines overlap with those of the star. The Doppler shift between the planet and the star during the transit is in the range ± 0.7 km·s −1 , which is too low to use a cross-correlation method as described by, for example, Snellen et al. (2010).

Summary and conclusions
We analysed the transmission spectrum of K2-18b between 0.43 and 5.02 µm using a combination of K2, HST, and Spitzer data retrieved from Benneke et al. (2019a) and Tsiaras et al. (2019). We studied the effect of irradiation, metallicity, clouds, internal temperature, eddy diffusion coefficient, and C/O ratio using our self-consistent model Exo-REM and assuming an atmosphere primarily composed of H 2 and He. We also analysed the Benneke et al. (2019a) dataset with the retrieval algorithm TauREx 3, and provided a S/N estimation for the detection of CH 4 using CRIRES+.
We found that the data are compatible with a highly metalenriched atmosphere, between 65 and 500 (Z/H) against Tsiaras et al. (2019) data -or between 100 and 200 (Z/H) against Benneke et al. (2019a) data -when assuming a solar C/O ratio, and 40 (Z/H) when assuming a sub-solar C/O ratio. According to our results, the atmosphere of K2-18b appears quite similar to that of Neptune or Uranus and seems to follow the C/H-mass relationship of the giant planets of our solar system, as shown in Figure 23. Most of the extrasolar planets seem to have an O/H ratio below that expected from the C/H-mass relationship for solarsystem giant planets. Part of the explanation could reside in the fact that for most of these planets, the amount of CO and is unknown, leading to an underestimation of the O/H ratio (as pointed out by e.g. Wakeford et al. 2017;Welbanks et al. 2019). A summary of our results and results of other teams is displayed in Figure 24, and the spectra of our best fit for each of the explored scenarios are displayed in Figure 25. We also show that thick to no H 2 O-ice clouds are allowed by the data. Liquid H 2 O clouds are possible on planets similar to K2-18b but receiving at most 80% of its irradiation. With or without clouds, we found that the Bond albedo of the planet should be around 0.02. The other studied parameters are not well constrained by the data. We note that in most of the cases we studied, CH 4 absorption should dominate or be on par with that of H 2 O in the HST spectral window, as first outlined by Bézard et al. (2020).
Combining the work of Bézard et al. (2020) and our own retrievals using TauREx 3, we show that the discrepancy between our self-consistent results and the results from retrieval algorithms can be explained by either the discarding of CH 4 absorptions, as in the case of the Tsiaras et al. (2019)  10 7 10 5 10 3 10 1 Volume Mixing Ratio 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Pressure (Pa)  Accordingly, in addition to our nominal scenario (i.e. Z/H = 65-500 (Z/H) and a solar C/O), scenarios with a CH 4 -depleted atmosphere could satisfactorily fit the observed spectrum. These could be obtained with a high internal temperature ( 200 K) or a low C/O ratio ( 0.01), with H 2 O being the dominant absorber in HST/WFC3 spectral band if C/O 0.1. However, the high internal temperature scenario seems very unlikely, and the Cdepleted scenario requires the existence of an unknown process. Therefore, it seems that a CH 4 -depleted atmosphere, a scenario favoured by all other teams who have analysed the same data so far, is more difficult to defend than a nominal scenario with a standard abundance of CH 4 . Moreover, a spectrum dominated by H 2 O absorptions seems even more unlikely given the relatively small range of self-consistent solutions allowing for this scenario and the fact that it is favoured by retrieval algorithms only because of overfitting. Another possibility is that our thermochem-10 5 10 4 10 3 10 2 10 1 10 0 Volume Mixing Ratio 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Pressure (Pa)  ical model does not accurately describe the chemical conversion between CO and CH 4 in the deep atmosphere, as suggested by Benneke et al. (2019b) in their atmospheric analysis of GJ3470b (noting that GJ3470b is a much warmer planet with an effective temperature of ≈ 800 K). It should be possible to discriminate between these scenarios by acquiring data in different spectral ranges and/or at higher resolution using the incoming JWST and ARIEL space telescopes. In complement, ground-based instruments such as CRIRES+ could be used to detect individual lines of species. We estimate that with our nominal scenario, we could detect CH 4 with CRIRES+ in the H band with an 'optimistic' S/N of 4.7 after three hours of transit observation.
This study highlights the critical need for more precise spectroscopic data on sub-Neptunes. This will allow us to better characterise these intriguing bodies, which have no analogue in our Solar System. 9.30 9.32 9.34 9.36 9.38 9.40 Wavelength (