Tighter constraints on the atmosphere of GJ 436 b from combined high-resolution CARMENES and CRIRES+ observations

A. Peláez-Torres; A. Sánchez-López; L. Nortmann; M. López-Puertas; E. González-Álvarez; H. M. Tabernero; C. Jiang; D. Revilla; G. Morello; J. Orell-Miquel; E. Pallé; P. J. Amado; J. A. Caballero; I. Ribas; A. Reiners; A. Quirrenbach; D. Cont; S. Dreizler; A. Fernández-Martín; A. P. Hatzes; Th. Henning; F. Lesjak; D. Montes; A. Schweizer; T. Trifonov; F. Yan

doi:10.1051/0004-6361/202557570

Home

All issues

Volume 705 (January 2026)

A&A, 705 (2026) A256

Full HTML

Open Access

Issue		A&A Volume 705, January 2026


Article Number		A256
Number of page(s)		21
Section		Planets, planetary systems, and small bodies
DOI		https://doi.org/10.1051/0004-6361/202557570
Published online		27 January 2026

A&A, 705, A256 (2026)

Tighter constraints on the atmosphere of GJ 436 b from combined high-resolution CARMENES and CRIRES⁺ observations

A. Peláez-Torres¹^★, A. Sánchez-López¹, L. Nortmann², M. López-Puertas¹, E. González-Álvarez³, H. M. Tabernero⁴^,5, C. Jiang⁶^,7, D. Revilla¹, G. Morello¹^,8, J. Orell-Miquel⁹, E. Pallé⁶^,7, P. J. Amado¹, J. A. Caballero³, I. Ribas⁵^,4, A. Reiners², A. Quirrenbach¹⁰, D. Cont¹¹^,12, S. Dreizler², A. Fernández-Martín¹³, A. P. Hatzes¹⁴, Th. Henning¹⁵, F. Lesjak¹⁶^,2, D. Montes¹⁷, A. Schweizer¹⁸, T. Trifonov¹⁰^,19 and F. Yan²⁰

¹ Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, Genil, 18008 Granada, Spain
² Institut für Astrophysik und Geophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
³ Centro de Astrobiología (CSIC-INTA), ESAC, Camino bajo del castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain
⁴ Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
⁵ Institut de Ciències de l’Espai (CSIC-IEEC), Campus UAB, c/ de Can Magrans s/n, 08193 Bellaterra, Barcelona, Spain
⁶ Instituto de Astrofísica de Canarias (IAC), 38200 La Laguna, Tenerife, Spain
⁷ Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
⁸ INAF - Palermo Astronomical Observatory, Piazza del Parlamento 1, 90134 Palermo, Italy
⁹ Department of Astronomy, University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA
¹⁰ Landessternwarte, Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany
¹¹ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr 1, 81679 München, Germany
¹² Exzellenzcluster Origins, Boltzmannstrasse 2, 85748 Garching, Germany
¹³ Centro Astronómico Hispano en Andalucía (CAHA), Observatorio de CalarAlto, Sierra de los Filabres, 04550 Gérgal, Spain
¹⁴ Thüringer Landessternwarte Tautenburg, Sternwarte 5, 07778 Tautenburg, Germany
¹⁵ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
¹⁶ Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
¹⁷ Departamento de Física de la Tierra y Astrofísica and IPARCOS-UCM (Instituto de Física de Partículas y del Cosmos de la UCM), Facultad de Ciencias Físicas, Universidad Complutense de Madrid, 28040 Madrid, Spain
¹⁸ Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany
¹⁹ Department of Astronomy, Faculty of Physics, Sofia University “St Kliment Ohridski”, 5 James Bourchier Blvd., 1164 Sofia, Bulgaria
²⁰ Department of Astronomy, University of Science and Technology of China, Hefei 230026, PR China

^★ Corresponding author: This email address is being protected from spambots. You need JavaScript enabled to view it.

Received: 6 October 2025
Accepted: 24 November 2025

Abstract

Context. Transmission spectra of Neptune-sized exoplanets are frequently observed to be featureless at low-to-mid resolutions from space; whereas high-altitude clouds can mute spectral features, high atmospheric metallicities can also result in compressed envelopes, where low scale heights may also yield undetectable signatures.

Aims. We aim to study the atmospheric properties of the warm Neptune GJ 436 b by combining a set of five transit events observed with the CARMENES spectrograph with one transit from CRIRES⁺ so as to provide the most constrained results possible at high resolution.

Methods. We removed telluric and stellar signals from the data using SysRem and potential planetary signals were investigated using the cross-correlation technique. Following standard procedures for undetected species, we performed injection recovery tests and Bayesian retrievals to place constraints on the detectability of the main near-infrared absorbers. In addition, we simulated ELT/ANDES observations by computing end-to-end in silico datasets with EXoPLORE.

Results. No molecular signals were detected in the atmosphere of GJ 436 b, which is consistent with previous studies. Combined CARMENES-CRIRES⁺ injection-recovery and Bayesian retrieval analyses show that the atmosphere is likely covered by high-altitude clouds (~1 mbar) at low and intermediate metallicities or, alternatively, is very metal-rich (≳ 900× solar), which would suppress spectral features without invoking clouds. Simulations of ELT/ANDES observations suggest a boost by nearly an order of magnitude to the upper limit in the photon-limited regime, reaching 0.1 mbar at 10-300× solar metallicities.

Conclusions. The joint analysis of all useful transit observations from CARMENES and CRIRES⁺ provides the most stringent constraints to date on the atmospheric properties of GJ 436 b. Complementary CCF-based and retrieval approaches consistently indicate that the atmosphere is either cloudy or highly metal enriched. Any weak near-infrared absorption lines, if present, are likely to be below current detection limits. However, according to our simulations, these features may be revealed with ELT/ANDES even in single-transit observations.

Key words: methods: observational / methods: statistical / techniques: spectroscopic / planets and satellites: atmospheres / planets and satellites: individual: GJ 436 b

© The Authors 2026

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Planets similar to and smaller than Neptune present a great challenge for atmospheric characterization due to their generally low transit signals, when compared to the widely studied hot and ultra-hot Jupiters (Bean et al. 2010; Berta et al. 2012; Kreidberg et al. 2014; Benneke et al. 2019a, 2024; Gandhi et al. 2020; Madhusudhan et al. 2021, 2023; Dash et al. 2024; Grasser et al. 2024). Following the discovery of GJ1214b (8.14M_⊕; Charbonneau et al. 2009), a benchmark sub-Neptune orbiting an M dwarf, a thorough analysis of its atmosphere using the Hubble Space Telescope (HST) revealed a featureless transmission spectrum (Kreidberg et al. 2014). Since then, multiple Neptune- and sub-Neptune-like planets have shown similar properties (Demory et al. 2016; Guilluy et al. 2021; Alam et al. 2022; Kreidberg et al. 2022; Brande et al. 2022). This lack of spectroscopic features is often explained by the presence of high-altitude cloud decks, highly compressed (i.e. high-mean-molecular-weight) atmospheres, or even by these exoplanets being airless bodies (Knutson et al. 2014a,b; Wakeford et al. 2018; Benneke et al. 2019a; Kreidberg et al. 2022). Indeed, rocky planets with significant atmospheres and gas giants in the Solar System generally have widespread cloud coverage (Irwin 2009; Marley et al. 2013).

However, observations with the James Webb Space Telescope (JWST; Gardner et al. 2006) using the NIRCam instrument have revealed the presence of water vapour (H2O), methane (CH₄), sulphur dioxide (SO₂), and carbon monoxide (CO) in the atmosphere of the Neptune-like planet GJ 3470 b, which has an estimated atmospheric metallicity of approximately 100× solar. This contrasts with prior HST and Spitzer observations, which suggested an atmospheric metallicity close to solar (Benneke et al. 2019a; Beatty et al. 2024). Similarly, highly metal-enriched envelopes have been inferred for several sub-Neptunes, including TOI-270 d (50% metals, ~500×solar; Benneke et al. 2024), GJ1214b (~1000×solar; Kempton et al. 2023; Schlawin et al. 2024; Ohno et al. 2025), and K2-18b (–100-300× solar; Charnay et al. 2021). In the case of K2-18b, an initial interpretation of its HST Wide Field Camera 3 (WFC3) transmission spectrum suggested the presence of H₂O, but this signal was later determined to be CH₄ with JWST (Benneke et al. 2019b; Barclay et al. 2021; Madhusudhan et al. 2023), whose absorption bands are difficult to distinguish from those of water vapour at the resolving power of WFC3. Besides CH₄, these recent JWST observations have led to the first confirmed detection of CO₂ in K2-18 b (Benneke et al. 2019b; Barclay et al. 2021; Madhusudhan et al. 2023).

Another interesting target in the Neptune-like range is GJ 436 b (Butler et al. 2004). Comprehensive studies of the GJ 436 planetary system have revealed its unusual properties. Namely, this transiting warm Neptune possesses a steeply polar orbit (Bourrier et al. 2022), as well as a long hydrogen tail, observed in Lyα, which results from its escaping atmosphere (Bourrier et al. 2016). Spectroscopic observations have yielded non-detections of the metastable helium triplet He I λ1083 nm (Nortmann et al. 2018; Guilluy et al. 2024; Masson et al. 2024), which is possibly caused by a low He concentration (Rumenskikh et al. 2023), and prevents further characterization of the escaping gas. Interestingly, a recent study by Revilla et al. (2025) seems to provide evidence of star-planet magnetic interactions, which lead to a constraint on the planet’s magnetic field and magnetospheric radius (Loyd et al. 2023; Vidotto et al. 2023). Contrary to popular belief, the presence of a planetary magnetic field does not necessarily protect a planet from stellar wind erosion but it can even increase the rate of ion escape (Ramstad & Barabash 2021).

The atmosphere of GJ 436 b was also the subject of early investigations with Spitzer transit and eclipse observations, which led to contrasting claims about the presence of CO and CH₄, as well as the influence of stellar variability on those measurements (Beaulieu et al. 2011; Knutson et al. 2011). Based on chemical-equilibrium modelling, Stevenson et al. (2010) predicted that the atmosphere of GJ 436 b could be rich in CO, with a significant depletion of CH₄. This is somewhat unusual as CH₄is expected to be the dominant carbon-bearing species in warm hydrogen-rich atmospheres (Seager 2013). Vertical mixing, combined with a high atmospheric metallicity (>10× solar), could lead to enhanced CO abundances, while a depletion of CH4 may arise via non-equilibrium chemical processes (Madhusudhan & Seager 2011). According to the model of Moses et al. (2013), the claimed composition could be consistent with a high metal-licity in the range of ∼230-2000× solar. High metallicities lead to a low atmospheric scale height (H ∝ m⁻¹, where m is the mean molecular mass), which reduces the amplitude of spectral features. However, reanalyses of the Spitzer data indicate a constant transit depth across both time and wavelength (3-24 μm; Lanotte et al. 2014; Morello et al. 2015). Furthermore, optical and near-infrared transmission spectra obtained with HST did not reveal any significant absorption features (Knutson et al. 2014a; Lothringer et al. 2018). More recently, GJ 436 b’s atmosphere was analysed at low resolution spectroscopy (LRS) with JWST using Near Infrared Camera (NIRCam) and Mid Infrared Instrument (MIRI) by Mukherjee et al. (2025), revealing a tentative detection of CO₂ (2σ). However, the absence of strong transmission features in previous high-resolution spectroscopy (HRS) studies suggests the presence of clouds. These findings have led to a general consensus that GJ 436 b’s atmosphere could contain a high-altitude cloud deck that mutes most of the spectral features across the near- (NIR) and mid-infrared (MIR) wavelength ranges.

High resolution spectroscopy has emerged as one of the most effective techniques for characterizing exoplanet atmospheres, both in transit and dayside geometries, across optical and infrared wavelengths (e.g. Snellen 2025). Furthermore, HRS is highly sensitive to small Doppler shifts in the spectral lines and to their shape, which enables us to gain information about exo-atmospheric circulation (Snellen et al. 2010; Miller-Ricci Kempton & Rauscher 2012; Flowers et al. 2019). HRS characterization of exoplanet atmospheres has enabled the detection of single species (Charbonneau et al. 2002; Snellen et al. 2008, 2010; Brogi et al. 2012; Birkby et al. 2013; de Kok et al. 2013; Lockwood et al. 2014; Casasayas-Barris et al. 2021; Stangret et al. 2022) and multiple molecules in the same atmosphere (Hawker et al. 2018; Cabot et al. 2019; Kesseli et al. 2020; Giacobbe et al. 2021; Cont et al. 2021; Sánchez-López et al. 2022), and has eventually evolved into the development of Bayesian retrievals (Brogi et al. 2017; Brogi & Line 2019; Gibson et al. 2020, 2022; Cont et al. 2022; Maguire et al. 2024; Blain et al. 2024). These retrieval techniques are capable of constraining relevant parameters (e.g. atmospheric elemental abundances, temperature structure and winds, and the planet’s rotational velocity, among others) using Bayesian statistics. Their application has yielded robust statistical significance in detecting molecular signatures from close-orbiting gas giants observed by transmission and emission spectroscopy (Lesjak et al. 2023; Cont et al. 2024; Nortmann et al. 2025). Furthermore, low-resolution retrievals (Madhusudhan & Seager 2009; Madhusudhan et al. 2014; Evans et al. 2016; Pinhas et al. 2018; Wakeford et al. 2018; Chachan et al. 2019; Chubb et al. 2020) or photometric data (Cont et al. 2024, 2025) can be combined with high-resolution studies to provide tighter constraints on the retrieved parameters (Brogi et al. 2017; Boucher et al. 2023; Smith et al. 2024). In contrast to low-resolution studies, HRS data allow us to distinguish the core and wings of spectral lines, which makes this technique ideal to attempt detections in hazy exo-atmospheres where only line cores may protrude above the continua (e.g. Pino et al. 2018a,b; Sánchez-López et al. 2020).

A recent study by Grasser et al. (2024) using the highresolution spectrograph CRIRES⁺, on European Southern Observatory’s (ESO) Very Large Telescope (Kaeufl et al. 2004; Follert et al. 2014; Dorn et al. 2014, 2023) did not detect transmission signals of H2O, CH4, or CO in the atmosphere of GJ 436 b in the 1431-1837 nm spectral range covered by the H1567 setting. However, they were able to place upper limits on the atmospheric composition by performing injection recovery tests of model absorption signals, which reveals that CRIRES⁺ observations should have detected water vapour if GJ 436 b presented a cloud deck at pressures higher than 10 mbar (lower altitudes) and if, simultaneously, the atmospheric metallicity was in the range of 1-300× solar.

In this paper, we compile all available datasets observed with the Calar Alto high-Resolution search for M dwarfs with Exo-earths with Near-infrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2014, 2018) and all observed with the CRyogenic high-resolution InfraRed Echelle Spectrograph (CRIRES⁺; Kaeufl et al. 2004; Follert et al. 2014; Dorn et al. 2014) for GJ 436 b (six transits in total) to explore the atmospheric composition of this exoplanet. To this end, we used the cross-correlation technique over multiple transit events.

The paper is structured as follows. In Sect. 2, we describe the observations. In Sect. 3, we refine both the stellar and planetary parameters of the GJ 436 system. The methods used for the analysis of the spectral data are detailed in Sect. 4, including those used for identifying the molecular signatures. We present and discuss our results in Sect. 5, including future prospects for atmospheric characterization of cloudy sub-Neptunes with the ArmazoNes high Dispersion Echelle Spectrograph (ANDES) (Marconi et al. 2022; Palle et al. 2025) on the Extremely Large Telescope (ELT), and we outline our main conclusions are given in Sect. 6.

Fig. 1

Spectral coverage of water vapour transmission models for GJ436 b of the CARMENES NIR channel (grey) and the CRIRES⁺ H1567 setting (red). Similar figures for CH₄ and CO can be found in Appendix A.1.

2 Observations

2.1 CARMENES

We analysed five transit datasets of GJ 436 b observed with the CARMENES instrument, installed at the 3.5 m telescope of the Centro Astronómico Hispano en Andalucía (CAHA). CARMENES has two spectral channels: the optical channel (VIS), covering the wavelength range of 0.52-0.96 μm in 55 spectral orders (R≈ 94 600), and the near-infrared channel (NIR), covering from 0.96-1.71 μm in 28 orders (R≈ 80 400). The CARMENES NIR spectral coverage includes that of the CRIRES⁺'s H1567 setting, as shown in Fig. 1. Each CARMENES channel has two fibres: fibre A was used to observe the target star, while fibre B was placed on the sky (88 arcsec to the west) to identify sky emission lines. Here, we focus on the observations conducted using the NIR channel, since this is the spectral region where the ro-vibrational lines from the spectroscopically active molecules are located. All of our datasets were reduced by the CARMENES pipeline, caracal v2.20 (Caballero et al. 2016). In Sect. 3, we revised all stellar and planetary parameters required for this study. The five datasets used for GJ 436 b (see Table 1 and Fig. 2) covered the pre-, peri-, and post-transit phases of GJ 436 b, except for the pre-transit of the fifth night (see Fig. 3). A sixth (03 December 2017) and a seventh (09 April 2018) night were available, but Due to the lack of in-transit spectra from the sixth night, and the high relative humidity and low signal to noise ratio (S/N) of the seventh, we excluded them from our analysis. The S/N across all five used nights ranges from ∼60 to ∼100. However, during certain exposures on the first and last nights, it fell below 60. Consequently, we excluded four exposures from the first night and three from the fifth to avoid low-quality spectra. During the course of all five nights, the star moved from airmass 1.02 to 2.06. We also excluded eight exposures from the second night and four from the third, as they were taken at an airmass greater than 1.75. In addition, our fourth night showed a 30 minutes gap (from φ = –0.040 to φ = –0.027) due to a guiding issue. Additionally, the S/N and relative humidity in the first part of that night exhibited irregularities. To ensure consistency, we also excluded that part of the night from our analysis.

2.2 CRIRES⁺

We analysed an additional transit dataset observed with the CRIRES⁺ spectrograph, installed at the ESO’s Very Large Telescope (VLT) at Cerro Paranal. The observation was carried out on 23 January 2023 with the H1567 wavelength setting and it was previously analysed by Grasser et al. (2024). This setting covers the 1.49 −1.78 μm range over eight spectral orders, at a resolving power of R ≈ 100 000. Further details on the observing conditions of this night are provided in Table 1. Two additional transit events were observed with CRIRES⁺ but, as described by Grasser et al. (2024), their signal-to-noise ratio is rather low due to issues in the adaptive optics system.

The CRIRES⁺ raw data for GJ 436 b were retrieved from the ESO Archive¹ and reduced using the standard CR2RES pipeline routines (version 1.4.4) through the ESO Recipe Execution Tool (EsoRex, version 3.13.8). This pipeline includes dark calibration (not used during nodding subtraction), bad-pixel masking, flat calibration, wavelength calibration using a combination of the Fabry-Pérot etalon and uranium-neon lamp, and 1D spectral extraction. The observation consisted of 7.5 nodding cycles, following the “A_1,2,3 B_1,2,3 - B_4,5,6 - A_4,5,6” nodding pattern. In order to perform nodding subtraction, we regrouped the nodding pairs as A₁B₁, A₂B₂, ... A₆B₆. The spectra were extracted using the optimal extraction algorithm, with a fixed extraction height of 20 pixels (∼8 times the median full width at half maximum of the point spread function along the slit) and no-light rows subtracted. Also, we discarded the bluest order #9 (1.43-1.46 μm) due to incomplete wavelength coverage and strong telluric absorption. The observation was conducted in in adaptive-optics mode mode and that, due to seeing conditions (∼0.65″), the slit was not evenly illuminated with light in the dispersion direction during the observation. Although this might cause a slight spectral shift between A and B positions (as described in the CR2RES user manual²), no significant spectral shift (< 1 kms⁻¹) was detected when cross-correlating the A and B spectra.

3 Estimation of stellar and planetary parameters

3.1 Stellar parameters

We computed the stellar atmospheric parameters, namely T_eff, log g, and [Fe/H], on the CARMENES template spectrum corrected for telluric absorption with SteParSyn³ (Tabernero et al. 2022) using the line list and model grid described by Marfil et al. (2021b). The luminosity was computed from the integration of the spectral energy distribution following Cifuentes et al. (2020) using photometry from B to W4 alongside the Gaia DR3 parallax (Gaia Collaboration 2023). The stellar radius R_* was determined from the Stefan-Boltzmann law using the STEPARSYN T_eff and the stellar mass M_* from the linear mass-radius relation given by Schweitzer et al. (2019). Values for both stellar atmosphere parameters and derived stellar parameters are compiled in Table 2. These values agree within 1σ with previous determinations (Marfil et al. 2021a; Passegger et al. 2022; Cristofari et al. 2022; Antoniadis-Karnavas et al. 2024; Melo et al. 2024, just to cite a few recent examples), but we preferred to use our own homogeneous determination of stellar and planetary parameters, as presented below.

Table 1

Analysed spectroscopic observations of GJ436.

Fig. 2

Evolution of the airmass, relative humidity, and mean S/N per spectrum over the entire CARMENES NIR range. The orange shaded area marks the in-transit times. Excluded spectra are represented with open symbols.

Fig. 3

Schematic of transmission spectroscopy observations of GJ 436 system used in this work. The orbital phase coverage for the different transits is shown in the lower part.

Table 2

Calculated and selected parameters for GJ436.

3.2 Planetary parameters

3.2.1 TESS photometric data

The Transiting Exoplanet Survey Satellite (TESS) observed GJ 436 in 2 min short-cadence integrations in two different sectors, S22 and S49, in March 2020 and March 2022, respectively. All sectors were processed by the Science Processing Operations Center pipeline (SPOC, Claudi et al. 2016) and searched for transiting planet signatures with an adaptive wavelet-based transit detection algorithm (Jenkins et al. 2010). In a first analysis of sector 22, GJ 436 was announced on 15 April 2020 as a TESS object of interest (TOI) via the dedicated MIT TESS data alerts public website⁴ where the planet was identified at an orbital period of 2.64397±0.00003 d. The planet was fitted with limb-darkening transit models by SPOC (Li et al. 2019) and successfully passed all the diagnostic tests performed on them (Twicken et al. 2018). No further transiting planet signatures were detected in the residual light curve by the SPOC. The planet candidate, GJ 436 b, has an estimated radius of 3.76±0.12 R_⊙ and a depth of 6480±56 ppm.

We downloaded the light curve files for the two sectors from the Mikulski Archive for Space Telescopes. The first step was to verify that the simple aperture photometry (SAP) and pre-search data conditioning SAP (PDCSAP) fluxes automatically computed by the pipeline are useful for scientific studies. For that, we checked that no additional bright source is contaminating the aperture photometry. Figure A.2 displays the target pixel files (TPFs) of GJ 436 for the two sectors using the publicly available tpfplotter⁵ code (Aller et al. 2020), which overplots the Gaia Data Release 3 (DR3) catalogue on top of the TESS TPFs. We confirmed that there are no additional Gaia sources within the photometric aperture around GJ 436 automatically selected by the pipeline. Therefore, we considered the extracted TESS light curves to be free of contamination from nearby stars.

3.2.2 HIRES, HARPS, and ESPRESSO spectroscopic data

GJ 436 b is a well-characterized planet that was initially detected with HIRES by Butler et al. (2004). The planet has an orbital period of 2.64 days, a minimum mass of 23 M_⊕, and an eccentricity of 0.15. Gillon et al. (2007) later confirmed it as a transiting planet with a size and mass comparable to Neptune. Trifonov et al. (2018) combined 113 CARMENES precise Doppler measurements with those available in the literature from HIRES and HARPS to derived new orbital parameters for the system. The updated Keplerian parameters of GJ 436 b, based on the modelling of all 638 Doppler measurements, were a planet semiamplitude K_b = 17.38ms⁻¹, period of P_b = 2.644 days, and eccentricity e_b = 0.152. The inclination constraints from the transit (i = 85.80±0.25 deg) yield a planetary mass of m_b = 21.4±0.2M_⊙ and semi-major axis of a_b = 0.028±0.001 au. In our new analysis, we also included the ESPRESSO RV data used by Bourrier et al. (2022) to analyse the Rossiter-McLaughlin (RM) signal. We took advantage of the RM anomaly induced by GJ 436 b on the RV residuals measured during the two transits observed with ESPRESSO (see Fig. 2 of Bourrier et al. 2022) to exclude the RV phases affected by the RM throughout all the archival spectroscopic data collected here.

3.2.3 CARMENES spectroscopic data

We had additional CARMENES observations of GJ 436 beyond those presented by Trifonov et al. (2018), aimed at refining the planetary mass measurement. We used a total of 176 highresolution spectra from 2016 January 08 to 2024 May 25. Relative RVs were extracted separately for each échelle order using the serval software (Zechmeister et al. 2018). The final VIS RVs per epoch were computed as the weighted RV mean over all échelle orders of the respective spectrograph, with an rms of 11.97ms⁻¹ and a mean error bar of 1.76 ms⁻¹.

3.2.4 Joint photometric and spectroscopic analysis

We derived the planetary parameters of GJ 436 b by performing a joint analysis of photometric and spectroscopic data, combining TESS light curves with RV measurements from HIRES, HARPS, CARMENES VIS, and ESPRESSO. The modelling was carried out using the juliet Python package, which incorporates radvel (Fulton et al. 2018) for Keplerian RV fitting and batman (Kreidberg 2015) for transit modelling. Outliers in the photometric data were filtered using a 3σ clipping procedure.

Instead of directly fitting for the planet-to-star radius ratio (R_p/R_*) and the impact parameter (b), we adopted the r₁ and r₂ parametrization introduced by Espinoza et al. (2019), which allows both parameters to vary between 0 and 1 while fully covering the physically allowed space for R_p/R_* and b. We imposed a prior on the stellar density (ρ_*), rather than on the scaled semimajor axis (a), ensuring a single consistent value of ρ* across the system. To optimize computational efficiency, given the volume of TESS data points and the lack of additional transiting signals, we restricted the photometric dataset to 2.5 hour windows centred on the predicted mid-transit times of GJ 436 b (with transits lasting approximately 1 hour). We defined normal priors for the parameters to the RV and light curve data P, t₀, r₁, r₂, centred on the values previously obtained from independent fits to the radial velocity and light curve data.

Given the generalized Lomb-Scargle periodogram of the RV data our joint model did not include a stellar activity component and omitted the use of a Gaussian process. Moreover, the expected RV semi-amplitude of the planet is sufficiently large that any potential stellar activity signal would have a negligible impact on the results. Based on the evidence reported in the literature, we modelled the system with one Keplerian signal: a 2.64-day transiting planet also detected in the RV data.

The posteriors from the one-planet fit that we adopted as planetary parameters for the GJ 436 system are presented in Table 3. The folded light curves, with the two TESS sectors combined and phased with the orbital period of the transiting planet, are shown in Fig. A.3. The corner plot displaying the posterior distributions of some of the planetary parameters obtained from the joint fit is presented in Fig. A.4. The resulting RV model folded in phase is displayed in Fig. A.5.

The combined analysis of the TESS, HARPS, HIRES, ESPRESSO, and CARMENES data yields an orbital period of 2.6438972±0.0000004days. From the RV amplitude of 17.3 $_{- 0.16}^{+ 0.16}$ $^{+0.16}_{-0.16}$ ms⁻¹ and an orbital inclination of 86.52 $_{- 0.45}^{+ 0.49}$ $^{+0.49}_{-0.45}$ deg, we derived a true planetary mass of 21.02 $_{- 0.72}^{+ 0.69}$ $^{+0.69}_{-0.72}$ M_⊙, and radius of 3.96 $_{- 0.18}^{+ 0.19}$ $^{+0.19}_{-0.18}$ R_⊙ for the transiting planet GJ436b.

An additional analysis including the complete TESS light curves, fitting a GP simultaneously to detrend the light curves while modelling the transit signal, was performed. The obtained values were fully consistent with those reported in Table 3, within 1σ error bars. The planetary parameters previously reported in the literature (e.g. Trifonov et al. 2018; Maxted et al. 2022; Bourrier et al. 2022) are also compatible within the 1σ level.

4 Methods

4.1 Normalization, outliers, and masking

Ground-based observations of exo-atmospheres suffer from an inherent disadvantage that stems from the varying conditions of Earth’s atmosphere (e.g. seeing and precipitable water vapour (PWV) variability). Such variations during the observations induce a fluctuating baseline in the spectra (Figs. 4a and 4b). Thus, to ensure a consistent level across all spectra, we performed for both CARMENES and CRIRES+ an order-by-order normalization by fitting a third-degree polynomial to the pseudocontinuum (see Fig. 4c).

Similar to the approach used by Kesseli & Snellen (2021), we removed contamination from cosmic rays by applying a 3σ-clipping procedure where we flagged values that significantly deviate from the mean of the pixel’s time series. These values were corrected by interpolating over their nearest neighbours in wavelength. Next, we masked the most opaque Earth’s atmospheric transmittance (telluric) windows, where 80% of the normalized flux was absorbed for CARMENES. This value was relaxed to 85% in the case of CRIRES⁺, following Grasser et al. (2024). A safety window of 5 pixels was added to each flagged pixel to further limit potential contamination of strong telluric line wings. The percentage of masked pixels varies significantly for different spectral orders and it is shown, for each night, in Fig. A.6. For the CRIRES⁺ dataset, we removed the first and last 100 pixels of each spectral order.

Table 3

Calculated and selected parameters for GJ 436 b.

4.2 Telluric and stellar correction

Telluric and stellar features within the spectra are orders of magnitude stronger than the Doppler-shifted excess absorption from the planetary atmosphere. In order to disentangle the exo-atmospheric signal, we used SysRem (Tamuz et al. 2005; Mazeh et al. 2007), an iterative principal component analysis algorithm that accounts for unequal uncertainties for each data point. SysRem has been successfully applied in similar studies in the past to reveal signatures in the atmosphere of several exoplanets (de Kok et al. 2013; Birkby et al. 2013, 2017; Sánchez-López et al. 2019; Nugroho et al. 2020; Cont et al. 2022; Nortmann et al. 2025). In addition, recent work by Maguire et al. (2024) suggests that SysRem yields a better performance for our purposes, compared to other common telluric-removal techniques (e.g. Molecfit; Smette et al. 2015; Kausch et al. 2015).

Each spectral order is affected by different levels of contamination, essentially depending on the spectral shape of the telluric transmittance. As a result, the number of required SysRem passes required is order-dependent and a criterion needs to be set to determine it. This is a difficult step that can lead to unintended biases and even to spurious planet-like signals, depending on the employed criterion (see e.g. Cabot et al. 2019; Cheverall et al. 2023). Following the approach used by Herman et al. (2020, 2022); Deibert et al. (2021); Ridden-Harper et al. (2023); Rafi et al. (2024); Parker et al. (2025), we studied the variance of the residual spectral matrices after applying SysRem and halted the algorithm when the standard-deviation difference between two consecutive passes (∆σ) was below a given threshold (1%) and started to plateau. This metric therefore encapsulates the percentage change in the standard deviation of the data for the application of each SysRem iteration. Thus, the ∆σ metric is defined as $Δ σ = \frac{^{(i - 1)} σ -^{(i)} σ}{^{(i - 1)} σ},$ $\Delta \sigma=\frac{{ }^{(i-1)} \sigma-{ }^{(i)} \sigma}{^{(i-1)} \sigma},$ (1)

where ^(i-1)σ and ⁽ⁱ⁾σ are the standard deviation of the residuals in a spectral order before and after the i-th pass of SysRem respectively. In practice, this method helps recognize when ∆σ plateaus for each spectral order, indicating that the major spectral variations have been removed and no further passes are required. This method should preserve signals close to the noise level such as the exo-atmospheric absorption, while also being a model-independent approach, in contrast to injection recovery methods employed in the past (Alonso-Floriano et al. 2019; Sánchez-López et al. 2019, 2022).

In the defined metric, a plateau is typically reached between two to five SysRem passes for CARMENES datasets and five to seven for CRIRES⁺ dataset, depending on the spectral order, as illustrated in Figs. A.7 and A.8. During the first SysRem passes, telluric and stellar features still remain (panel D of Fig. 4), suggesting the need for further cleaning. For intermediate passes, the exo-atmospheric signal,if present, is expected to be buried in the noise (Fig. 4E).

Fig. 4

Steps of the data preparation applied to a representative spectral region of the dataset obtained on the night of February 04, 2017. The vertical axis in panel A represents flux in arbitrary units, while in all other panels it represents time, expressed as phase. Panel A: original spectrum observed at mid-transit, where telluric absorption lines from H₂O can be identified. Panel B: original spectral matrix extracted using CARACAL, where major S/N differences are observed between the spectra (horizontal stripes), and prominent telluric H₂O emission can be identified. Panel C: normalized spectra, where a 3σ clipping has been applied to filter out outliers and where telluric masks have been applied (white, excluded pixels). All spectra are now normalized to a common continuum. Panel D: resulting spectra after one SysRem pass, where major residuals can be observed, especially from emission lines. Panel E: spectra after nine SysRem passes, where most of the telluric contribution has been removed and any exoplanet signal, if present, is buried in the noise.

4.3 Model atmospheres

We modelled the exo-atmospheric transmission spectra using petitRADTRANS⁶ (pRT), a Python-based package for calculating transmission and emission spectra of exoplanets through radiative transfer (Mollière et al. 2019). We assumed H₂-He dominated atmospheres where Rayleigh scattering and collision-induced absorption are included (Figueira et al. 2009; Nettelmann et al. 2010). The same preparation steps performed on the real data were reproduced for the models to ensure that the same distortions are introduced.

We used pRT to model the absorption of the most relevant spectroscopically active species in the covered NIR range. Specifically, we considered H₂O and CO using the line lists presented by Rothman et al. (2010), NH₃ by Sousa-Silva et al. (2014), and CH₄ by Hargreaves et al. (2020). Model transmission spectra were generated on a two-dimensional grid spanning the pressure level of a potential cloud deck (p_c) and the metallic-ity of the atmosphere (Z). This latter step was facilitated by the algorithm easyCHEM⁷ (Lei & Mollière 2025), which provides with the mass fractions of the atmospheric compounds and the mean molecular weight (thus, the atmospheric scale height) as a function of the metallicity. Our p_c grid spanned from 0.1 mbar to 1 bar, while the Z grid covered from 1× solar to 1000× solar metallicity, assuming a solar C/O ratio of 0.55 (Asplund et al. 2009). Grid values were evenly spaced on a logarithmic scale, resulting in a total of 100 templates per molecule. Whereas only the respective species’ opacity is used to compute the transmission spectra in petitRADTRANS, easyCHEM does consider multiple species to assess the atmospheric compositions and scale height. Figure 5 illustrates the variations of the transit depth in the models due to changes in both p_c and Z.

4.4 Cross-correlation technique

Even after removing the telluric and stellar contributions, any potential atmospheric lines from GJ 436 b remain buried below the noise level. By using the cross-correlation technique, it is possible to co-add the contribution from potentially hundreds of these lines in the form of a cross-correlation function (CCF) peak (e.g. Snellen et al. 2010; Brogi et al. 2012; de Kok et al. 2013; Birkby et al. 2013, 2017). We performed the cross-correlation between the residual matrices obtained in the previous steps R_ij and the transmission models (templates) m_j calculated using pRT. In order to explore a wide range of velocities with respect to the Earth, we Doppler-shifted the templates in a range of velocities v (–325 to +325 km s⁻¹) with respect to the Earth’s rest frame by linear interpolation. We used 1.3 km s⁻¹ (CARMENES) and 1 kms⁻¹ (CRIRES⁺) intervals, as set by the average velocity step size between the instruments pixels. CCFs were obtained individually for each spectrum, i, forming a cross-correlation matrix, CCF(v, i), as follows: $C C F (v, i) = \sum_{j} \sum_{λ} \frac{R_{i, j} (λ) m_{j} (λ, v)}{{\hat{σ}}_{i, j}^{2} (λ)},$ $CCF(v,i) = \sum_j \sum_\lambda \frac{R_{i,\,j}(\lambda)\,m_j(\lambda,\,v)}{\hat \sigma_{i,\,j}^2 (\lambda)},$ (2)

where σ̂(λ) are the propagated uncertainties of the data and the summations loop over wavelength (λ) and spectral order j, so as to co-add the potential information contained in our full spectral coverage. The resulting total CCF(v, i) in the Earth’s rest frame is illustrated in Fig. 6, where the trace of GJ 436 b would be expected to appear during transit (i.e. between the horizontal dash-dotted red lines) along the planetary RVs with respect to the Earth, v_p. To properly describe the particular orbit of GJ 436 b, we modified the classic definition of v_p to account for an eccentric orbit as follows: $v_{p} (ϕ) = K_{p} \cdot [\cos (ν (ϕ) + ω_{p}) + e \cdot \cos (ω_{p})] - v_{bary} + v_{sys},$ $v_{\mathrm{p}}(\phi)=K_{\mathrm{p}} \cdot [\cos(\nu(\phi) + \omega_{\rm p})+e \cdot \cos(\omega_{\rm p})]-v_{\text {bary }}+v_{\text {sys}},$ (3)

where φ is the orbital phase, ν(φ) is the true anomaly, ω_p the argument of periastron, e is the eccentricity, v_bary the barycentric velocity due to Earth’s motion around the Solar System’s barycentre, v_sys the systemic velocity of the star-planet system, and K_p the radial velocity semi-amplitude of the planet, defined as $K_{p} = \frac{2 π a}{P_{o r b} \sqrt{1 - e^{2}}} \sin (i),$ $K_{\mathrm{p}}=\frac{2 \pi a}{P_{\rm orb} \sqrt{1-e^2}} \sin(i),\end{equation}$ (4)

where a is the orbital parameter (the semi-major axis), P_orb the orbital period, and i the inclination.

The eccentric υ_p has an average shift of about 20 kms⁻¹ and a slight change in slope compared to the circular velocity. This change highlights the importance of accurate orbital characterization, as failing to account for the correct orbital configuration could result in incorrect interpretations of potential signals (Basilicata et al. 2024; Grasser et al. 2024). It is at this point that we visually inspected the cross correlation matrices in the Earth’s rest-frame for all spectral orders and for each night separately. We discarded from the analysis those orders where the telluric contamination was not properly removed by SysRem (last column of Table 1). That is, when a trail of high CCF values along the vertical at 0 km/s is observed (i.e. where correlation with uncorrected telluric-H2O lines is expected to appear). We used the same list of H₂ O-based discarded orders for all CCF analysis, regardless of the studied molecule, as a conservative approach to minimize telluric noise and opaque wavelength windows as much as possible. As it is the usual approach, we also excluded spectral orders where the target molecule does not have significant spectral features (e.g. Parker et al. 2025).

Next, we co-added all in-transit CCFs along the time axis (i.e. sum over i in Eq. (2)) to maximize any potential signal. To do this, we first Doppler-shift CCF(υ, i) to the exoplanet’s restframe by using Eq. (3), assuming K_p as unknown and using a range of values from –300 to 300 km s⁻¹. Thereby, an absorption with a planetary origin should only appear around the expected K_p of GJ 436 b (i.e. the true exoplanet rest frame) in the form of a CCF peak.

Fig. 5

Transit depth of H₂O for GJ 436 b generated with pRT as a function of cloud deck pressure level (top panel, fixed metallicity of 10× solar) and atmospheric metallicity (bottom panel, clear atmosphere) as multiples of solar. For clarity, an offset of 0.002 in transit depth was applied to the transmission models in both panels. Clouds at higher altitudes mute more H₂O lines.

Fig. 6

Cross-correlation analysis of potential H₂O signals in the transmission spectrum of GJ 436 b observed with CARMENES in the near-infrared on the night of February 02, 2017. We show the crosscorrelation matrix, CCF (v,i), in the Earth’s rest frame as a function of the velocity Doppler shifts applied to the template (horizontal axis) and the planet’s orbital phase (vertical axis). The cross-correlation was derived using a 10× solar metallicity template with a cloud deck at 10 mbar, co-adding all useful spectral orders. Dashed-dotted black horizontal lines mark the start and end of the transit, while the dotted red line indicates the expected exoplanet velocities with respect to Earth, accounting for eccentricity.

4.5 Assessment of potential signals

4.5.1 S/N for CCF studies

Following common signal significance evaluation approaches, we computed a S/N map as a function of the exoplanet restframe velocity and K_p (top panel of Fig. 7; Birkby et al. 2017; Brogi et al. 2018; Alonso-Floriano et al. 2019; Sánchez-López et al. 2019; Nugroho et al. 2021; Cont et al. 2022, 2024). For each K_p, we divided each CCF point by the CCF’s standard deviation, excluding a ±20 km/s window around it. The division ensured that we avoided inflating the standard deviation by including potential CCF-peak wings. Potential signals are hence evaluated in a common K_p-v_rest map, as shown in the top panel of Fig. 7. In this case, the map does not reveal any significant signals (e.g. at S/N> 4) at the expected K_p (bottom panel of Fig. 7). All the observed patterns could be well explained by noise fluctuations or weak telluric residuals. This then leads to a non-detection of water vapour in GJ 436 b when using this particular CARMENES dataset observed on February 02, 2017. The same methods are applied subsequently to all transit datasets (Fig. A.9).

To enhance the recovery of a potential planetary signal, we combined the information from each observed night by creating a merged CCF map. In this process, we sorted the frames from all datasets along the time axis according to their corresponding orbital phases following the approach of Cont et al. (2022, 2025). The respective v_bary of each spectrum was considered at the time of Doppler-shifting the CCF matrix from the Earth’s rest-frame to the planet’s. This procedure hence yielded a single CCF map, which incorporates all the data from the different nights and is used to compute the K_p-v_rest and S/N maps in the same way as for a single-night dataset.

Fig. 7

Significance analysis for potential water vapour cross-correlation signals on the night of February 02, 2017. Top panel: signal-to-noise ratio map as a function of radial velocity in the exoplanet’s rest frame (v_rest, horizontal axis) and the projected orbital velocity semi-amplitude, K_p (vertical axis). Horizontal red lines denote the expected K_p, and vertical red lines mark the expected zero velocity, assuming eccentric orbits. Bottom panel: one-dimensional cross-correlation function at the expected K_p (102.3 km s⁻¹) of GJ436 b for the eccentric case.

4.5.2 Bayesian retrieval framework

To investigate the multidimensional parameter space of the transmission models, we employed the nested sampling algorithm MultiNest (Feroz et al. 2013), along with its Python wrapper PyMultiNest (Buchner 2016), to sample a three-dimensional space defined by the logarithm of the atmospheric metallicity in units of the solar value, log(Z/Z_⊙), the logarithm of the pressure level of the cloud deck, log(p_c), and the uncertainties-scaling parameter β. We adopted uniform priors on all free parameters, as summarized in Table 4, with a broad coverage of plausible atmospheric conditions. We set per-night retrievals with 300 live points each and a non-constant efficiency mode (e.g. Blain et al. 2024). These parameters and prior ranges served as coordinates for the algorithm, which assesses the goodness-of-fit of the model to the data by evaluating a log-likelihood function log(L). Here, we followed the framework of Gibson et al. (2020, 2022) and Maguire et al. (2024), with $\log L (F ∣ θ) = - \frac{1}{2} \sum_{N, λ, t} {(\frac{P (F (λ, t)) - P (M_{θ} (λ, t))}{β \hat{σ} (λ, t)})}^{2} - N \log (β),$ $\log \mathcal{L}(F \mid \theta) = - \frac{1}{2} \sum_{N,\,\lambda,\,t} \left(\frac{\mathcal{P}(F(\lambda,\,t)) - \mathcal{P}(M_\theta\,(\lambda,\,t))}{\beta\,\hat \sigma(\lambda,\,t)} \right)^2 - N \log(\beta),$ (5)

where θ = [log(Z/Z_⊙), log₁₀(p_c),β] are the parameter vectors, P represents a mathematical operator including all operations performed to the data in order to enable exo-atmospheric detections (i.e. preparations steps such as normalization and telluric correction), F represents our data spectral matrix, M_⊙ represents the model template computed from the parameters θ, β is a scaling parameter for the uncertainties, and where the sum loops over all useful nights N, spectral points and orders (λ), and intransit spectra (t). The same preparation steps performed on the real data were applied to the model so as to ensure the same distortions were introduced in both sets, hence increasing their likeness. Consequently, this method determines the best-fitting model to the data using Bayesian statistics (Brogi & Line 2019; Gibson et al. 2020, 2022; Cont et al. 2022; Blain et al. 2024).

Importantly, when running the retrieval, we slightly changed the normalization step described by Sect. 4.1 to follow exactly the same steps described in Gibson et al. (2022), which we summarize below. This change was to ensure that no statistical bias is introduced, while performing a significantly faster normalization step in all log(L) evaluations. In practice, this entails starting from the outlier-corrected CARMENES data, performing a normalization dividing each order by the median spectrum over time, and applying SysRem order by order following the criterion discussed above. For the determined number of SysRem passes for each order, we stored the resulting SysRem column vectors U and precomputed the SysRem -filtering projector P that was used to prepare models before comparing them to the data (Gibson et al. 2020, 2022; Cont et al. 2025). Essentially, a filtered model is obtained by dividing it by its median (M_norm), and then computing $M_{f i l t} = P (M_{n o r m}) = U (Λ U)^{†} (Λ M_{n o r m}),$ $\mathcal{M_{\rm filt}} = \mathcal{P}(\mathcal{M_{\rm norm}}) = U(\Lambda U)^{\dagger}(\Lambda \mathcal{M_{\rm norm}}),$ (6)

where Λ is the inverse of the mean σ over wavelength, and where the matrix of SysRem column vectors has an additional column of ones to account for potential offsets introduced by media division at the time of normalizing the data. Since the projector P = U(ΛU)^†Λ is precomputed, each retrieval iteration is greatly accelerated. Essentially, each likelihood evaluation entails a call to easyCHEM to provide the atmospheric mass fractions and mean molecular weight, a call to pRT to obtain the transmission spectrum with a given cloud contribution and including the opacities from H₂, He, H₂O, CH₄, NH₃, CO, and CO₂, followed by a preparation of the templates as per Eq. (6), and the subsequent log-likelihood computation from Eq. (5).

Each night’s dataset was retrieved independently to produce weighted posterior samples and per-night multivariate kernel density estimates (KDEs; non-parametric estimators of the posterior constructed from the weighted samples). For each night we computed the marginalized one- and two-dimensional posteriors as usual, where the posterior is formally defined as $\log P (θ ∣ F) = \log L (F ∣ θ) + \log π (θ) - \log Z,$ $\log P(\theta \mid F) = \log \mathcal{L}(F \mid \theta) + \log \pi(\theta) - \log \mathcal{Z},$ (7)

with L(F|θ) the likelihood of the data F given parameters θ, prior π(θ), and Z the Bayesian evidence. To obtain the combined posteriors for all nights, we built a candidate set of parameter vectors θ by drawing from the per-night KDEs, from the priors, and from the pooled per-night samples. Subsequently, we used the KDEs to evaluate the per-night posterior densities for each θ, denoted p_i(θ), and combined them in log-space according to $\log w (θ) = \sum_{i = 1}^{N_{n i g h t s}} \log p_{i} (θ) - (N_{n i g h t s} - 1) \log π (θ),$ $\log w(\theta) = \sum_{i=1}^{N_{\rm nights}} \log p_i(\theta) - (N_{\rm nights}-1)\,\log \pi(\theta),$ (8)

which corrects for the prior being included in each per-night posterior. The weights w(θ) were normalized and used to resample the candidate pool, producing combined samples from which medians and 68% credible intervals (16th, 50th, and 84th percentiles) and the final corner plot were derived. To ensure broad coverage of parameter space, the candidate pool combined KDE draws, prior draws, and pooled per-night samples. This procedure approximates the joint posterior for the combined dataset while preserving per-night posterior structure and avoiding prior double-counting.

This combination implicitly assumes that different nights are statistically independent (so the joint likelihood factorizes), have the same parameter definition, prior ranges, and functional form for the model and preparation pipeline (even if individual models or SysRem projectors differ for different spectral orders and nights), and that the per-night KDEs adequately approximate the night-wise posteriors. Due to this, it was possible for us to combine the retrieval information for log(Z/Z_⊙) and log₁₀(p_c), but this approach was not used for the β parameter, since it explicitly reflects per-night noise properties. Thus, we opted for presenting the retrieval of β separately for each night in Sect. 5.5.

Table 4

Prior ranges and posterior medians (68% credible intervals).

5 Results and discussion

We studied the presence of H₂O, CH₄, CO, and NH₃ in the atmosphere of GJ 436 b by applying the cross-correlation technique to five useful CARMENES datasets and one from CRIRES⁺. In the following, we present the results of direct CCF, CCF-based injection recovery of artificial signals and Bayesian retrievals.

Fig. 8

Detectability maps obtained from injection recovery tests, showing upper limits on the abundances of the main near-infrared absorbers in the atmosphere of G 436 b. Specifically, we investigated H₂O (left panel), CH₄ (middle panel), and CO (right panel) across a range of metallicities and cloud deck altitudes. Dashed lines indicate the S/N = 5 contours (‘SN5’, S/N = 3 for CO denoted as ‘SN3’) obtained from CARMENES (white, five transits combined), CRIRES⁺ (red, one transit), and combining the results from both instruments (green). The dotted pink lines mark the S/N = 5 (S/N = 3 for CO) contours from Grasser et al. (2024), referred to as ‘G24’, using the same CRIRES⁺ dataset. For CO, an additional contour with a potential offset corrected is also shown as ‘G24 offset corrected’ (see discussion in Sect. 5.4). The dash-dotted brown line corresponds to the 99.7% upper limit from Fig. 3 in Knutson et al. (2014a), referred to as ‘K14’. The latter values were extracted from their respective figures. The dashed purple line represents the 68% upper limit placed by Mukherjee et al. (2025), referred to as ‘M25’. The area below each curve represents the discarded atmospheric parameter values.

5.1 Non-detection of molecular signals with CCF

After performing the CCF analysis with an atmospheric model with p_c = 10 mbar and Z = 10× solar, we did not detect any signal of H₂O, CH₄, or CO in GJ 436 b’s atmosphere with either CARMENES or CRIRES⁺ (Fig. A.10). We additionally searched for the expected weak spectral features of NH3 , but no signals were obtained either (not shown). From this result and previous knowledge of the existence of a H₂ - He atmosphere in this planet, we can think of three possible scenarios: (i) GJ 436 b does not have significant amounts of H2O, CH4 , or CO in its atmosphere (which seems unlikely), (ii) our assumed model is significantly different from the real conditions in the atmosphere, a problem that is tackled by using Bayesian retrievals in Sect. 5.5, or (iii) the atmosphere is veiled by high-altitude clouds that suppress spectral features to the extent that even high-resolution spectroscopy cannot discern the line cores. We assess this point by discussing the results from different methodologies.

5.2 Upper limits from injection recovery tests

In the spectral range studied, H2O and CH4 exhibit prominent absorption features in the form of line forests, whereas CO shows only a relatively weak band around 1.5 μm with far fewer lines. As we did not detect the presence of any of these molecules in our data we proceeded to place limits on their detectability. To this end, we performed signal-injection tests using a grid of H2O, CH4 , and CO templates with varying pc from 0.1 mbar to 1 bar, and Z ranging from 1× solar to 1000× solar. The model signals were injected before the normalization step at the expected velocities of the exoplanet with respect to the Earth (Eq. (3)). We then applied our preparation and CCF methods for each injected case, and used the same model as a template for cross-correlation in that grid point. The resulting S/N of the recovered signals are shown in Fig. 8. We note that, in this process, the non-injected CCF is subtracted from the injected one to remove any potential noise structures at the injected position, hence ensuring the significance represented better the recovery of the injected atmosphere.

Focusing on the H₂O detectability map (Fig. 8, left), which is the most restrictive case of the three molecules, we combined the five CARMENES transits and determined that the most plausible scenario for GJ 436 b’s atmosphere is the presence of clouds at pressures lower (higher altitudes) than the 10 mbar level if the metallicity is low (1-2× solar), below the 1 mbar level for intermediate metallicities (10-100× solar), or cloud-free above a metallicity of ∼600× solar. For any other case, our detection thresholds indicate that, if real, we should have observed hints of H₂O absorption in our analysis. Interestingly, we observed CARMENES-derived upper limits (after combining transits) to be significantly tighter than the particular single-transit results from CRIRES⁺. We note that our CRIRES⁺ upper limits for water vapour agree rather well with the previous analysis of this dataset by Grasser et al. (2024). The difference between the results of both instruments for H₂O and the other species is further explored in Sect. 5.4.

When studying CH4 (Fig. 8, middle), we find CARMENES constraints to be significantly stronger than those derived with CRIRES⁺, and both CH₄ curves are generally less stringent than our H₂O upper limits. This does not come as a surprise, since methane’s intrinsically weaker lines result in a larger fraction of the parameter space being undetectable. However, we find significantly more restrictive constraints than reported in Grasser et al. (2024) when studying the same CRIRES⁺ dataset, a discrepancy that is addressed in Sect. 5.4. From CARMENES results, CH₄upper limits discard cloudless or low-altitude cloud (p_c > 6 mbar) scenarios for any metallicity < 60× solar.

Regarding CO (Fig. 8, right), we performed our analysis only for CARMENES and in the useful spectral orders where this molecule shows absorption lines (CARMENES orders 39 and 38), which are also covered by CRIRES⁺ orders 4 and 5. Since the CO injection-recovery tests were already performed for CRIRES⁺ by Grasser et al. (2024), and due to these H-band CO lines being relatively weak, we opted for not repeating this study using CRIRES⁺. Indeed, from Fig. 8, a detection threshold of S/N = 5 is not met for any combination of p_c-Z, and we could only show the S/N = 3 contours as a common reference. At first glance, we observed an apparent difference in sensitivity, with CRIRES⁺ constraints being slightly more stringent than CARMENES’. However, when slicing Fig. 8 (right panel) at different Z and p_c values (i.e. reproducing Fig. 6 from Grasser et al. 2024), we observed a consistent offset of about 0.4 between the baseline S/N levels of our CO signals and theirs (Fig. A.11). We were able to reproduce a nearly identical offset with the CARMENES data when using all spectral orders overlapping with the full CRIRES⁺ H1567 coverage. That is, when including also the spectral orders with no CO lines (which were originally discarded in our analysis). The inclusion of all CRIRES⁺ H1567 orders in Grasser et al. (2024) has been confirmed with the authors (Grasser; priv. comm.). After considering this offset (see “G24 offset corrected curve”), we found that the CO detection limits are nearly identical in both studies (see Sect. 5.4), but only allowed us to exclude a region with metallicities between 10× solar and 100× solar, with cloud decks ranging from 1000 to 100 mbar. In other words, CO studies did not improve the H₂O constraints either.

Overall, our CCF injection recovery tests are dominated by the most-restrictive water vapour case. We find that our results are consistent with the high-altitude cloud deck and/or high metallicity scenarios inferred by Figueira et al. (2009), Moses et al. (2013), Knutson et al. (2014a), Grasser et al. (2024), and Mukherjee et al. (2025) for GJ 436 b’s atmosphere.

5.3 CARMENES and CRIRES⁺ combined

The tightest possible constraints on the detectability of GJ 436 b’s atmosphere using cross-correlation are obtained by combining CRIRES⁺ and CARMENES observations. For this, we restricted our analysis to water vapour, as it is the case that yields the most stringent upper limits (Fig. 8, green).

Interestingly, the combined constraints with both instruments are very similar to those placed with CARMENES only. This might indicate that the signal-detectability in GJ 436 b’s atmosphere has a physical limit imposed by pc. That is, if the planet’s atmosphere has clouds at altitudes higher than the 1 mbar-level (which is a sensible assumption given ours and previous results), a few additional transit datasets with current instrumentation might not be help us disentangle molecular lines that could be almost completely muted. Similarly, the detectability of H₂O becomes extremely challenging or unfeasible both at the solar and 1000× solar metallicity boundaries, due to either its low abundance or highly compressed atmospheres. Hence, we would not expect these upper limits to significantly change until the arrival of the next-generation of ground-based telescopes and high-resolution spectrographs.

5.4 A study of instrumental capabilities

The different constraints between our CARMENES and CRIRES⁺ analyses are a combination of multiple factors, involving different telescope apertures (3.5 m for CAHA/CARMENES versus 8.2 m for the VLT/CRIRES⁺), instrumental properties (e.g. resolving power or total wavelength coverage), and observing conditions (e.g. different challenges involving the telluric correction). The first of these can be compensated, in principle, by combining several CARMENES transits. Despite the challenges of preparing multiple datasets, this allowed us to build a similar S/N to that of a single CRIRES⁺ observation at the VLT. That is, a direct comparison between single-transit datasets from CARMENES and CRIRES⁺ would provide significantly worse results for the former, given that the VLT does collect over five times more photons per exposure.

After combining all CARMENES observations, the broader spectral coverage of its NIR channel, compared to a single CRIRES⁺ setting, does indeed contribute substantially to obtaining more restrictive results. In order to illustrate this, we repeated the injection recovery tests shown in Fig. 8 for H₂O and CH₄ , but including only the CARMENES orders that approximately overlap with those of the CRIRES⁺ H1567 setting (see Figs. 1 and A.1), obtaining the injection-recovery maps shown in Fig. A.12. This approach, while still hampered by differences in other instrumental properties and observing conditions, allowed us to evaluate the impact of H₂O and CH₄ absorption bands in the 0.96-1.45 μm range. For H₂O, the results from Fig. A.12 show that CARMENES-derived upper limits worsened and became weaker than those of CRIRES⁺ at the low- and high-metallicity cases, as expected from discarding many useful spectral points. The more restrictive upper limit from CRIRES⁺ in these cases is likely due to a combination of its higher resolution (R=100 000) being able to resolve weaker H₂O cores at low Z, and more compressed atmospheres at high Z. However, we observe nearly identical constraints for intermediate metallicities (from 10 to 100× solar) and ∼3 mbar cloud decks, for which the water vapour content is significant and both instruments disentangle the signal with a similar performance (albeit combining multiple CARMENES observations). In the case of CH4 , we found a similar effect: CARMENES upper limits became more relaxed when reducing the spectral points included in the analysis, bringing the new S/N = 5 curve closer to CRIRES⁺ results. In this case, CARMENES-restricted constraints were still slightly more stringent. Since both instruments cover the strongest methane lines in the H-band between 1.6 and 1.7 μm with almost overlapping wavelength coverage, the origin of the difference is challenging to determine. Given the similar (though slightly worse) results obtained for H2O, a potential explanation might be a more difficult telluric removal in five different datasets for CARMENES, with respect to only one for CRIRES⁺ This issue might then be relaxed when studying CH4 , revealing a slightly better performance for CARMENES. However, providing definite evidence for this hypothesis is beyond the scope of this work.

Interestingly, we found the order-restricted CARMENES results to be now more similar to our CRIRES⁺-derived upper limits for CH4 , but both remained notably more restrictive than the curve presented in Grasser et al. (2024). Since we used the same CRIRES⁺ dataset and given the very consistent results for H₂O (between our analyses and the upper limits in Grasser et al. 2024), the origin of the difference is somewhat uncertain. Lastly, for CO, the CARMENES spectral coverage does not introduce additional spectral lines than those covered by CRIRES⁺ (see Fig. A.1) and thus, the results of Fig. 8 already depict comparable results for both instruments (after considering the potential CO offset in Grasser et al. 2024).

5.5 Bayesian retrievals

We performed a Bayesian retrieval of our six transit datasets obtained with CARMENES and CRIRES⁺, exploring a wider parameter space than with direct CCF and using a different metric to verify our previous results (see also Table 4). With Fig. 9, we confirm the non-detection of species in the atmosphere of GJ 436 b with Bayesian retrievals, as best-fitting models present primarily a low-pressure (high-altitude) cloud deck (pressures < 10 mbar) that mute any potential feature regardless of the atmospheric metallicity. For high-metallicity scenarios, at log₁₀(Z/Z_⊙) ≳ 500× solar, the atmospheric scale height becomes too small, and the transit signal is undetectable in our data regardless of the presence of clouds. For low and intermediate metallicities, we observe a general correlation between log₁₀(Z/Z_⊙) and log₁₀(p_c): the higher the former is, the higher the altitude of the cloud deck needs to be to mute stronger H₂O lines. In other words, this means our retrieval only finds good fits to the combined datasets when no prominent spectral lines are present in the templates, reflecting that no exo-atmospheric lines are detected in the data.

The retrieval thus ratified the constrains put by CCF-based injection recovery tests. In Fig. 9 we overplot the CCF injection-recovery curve corresponding to a detection threshold of S/N = 5. The region lying beneath that curve thus defines combinations of metallicity and cloud-top pressure for which our CCF pipeline would have produced a detectable signal, while the posterior region returned by the Bayesian retrieval encloses those parameter combinations that are consistent with the observed spectra (i.e. effectively no-signal scenarios). These two metrics are therefore complementary and mutually reinforcing: the retrieval-allowed area occupies parameter space that lies predominantly outside the CCF-detectable domain, meaning that the models consistent with the data would not have been picked up by our CCF search at S/N = 5.

Regarding the uncertainties-scaling parameter, we find it to range between 0.75 and 0.80 for CARMENES data (see Table 4), which is in line with previous retrievals of this instrument using this or similar frameworks (Yan et al. 2022; Guo et al. 2024; Blain et al. 2024; Cont et al. 2022, 2025). This means that we find only a slight overestimation of uncertainties by CARACAL, or slight overcorrection of spectral features by our SysRem -based preparation pipeline. However, a retrieval on the CRIRES⁺ dataset reveals a β_CRIRES = 1.432 ± 0.001, which indicates a larger discrepancy between the pipeline uncertainties we derived and the scatter in the prepared CRIRES⁺ data. This suggests either an underestimation of noise by the standard error propagation in our pipeline or additional residual systematics not being captured by the model (e.g. telluric residuals). Indeed, when preparing the data using the methods of Gibson et al. (2022) prior to the retrieval, we observed some telluric residuals. However, the retrieved β_CRIRES did not decrease significantly when increasing the number of SysRem passes. Hence, it is likely that the issue might be a combination of both some telluric residuals and slight underestimation of the uncertainties by our raw-data reduction pipeline.

Fig. 9

Posterior distributions from atmospheric retrievals of six transit datasets of GJ 436 b, five from CARMENES, and one from CRIRES⁺. Panel A: combination of all useful transit datasets from both instruments. This yields the marginalized one-dimensional posteriors shown in the diagonal panels for log₁₀(Z/Z_⊙) and log₁₀(p_c), with vertical lines at the 16th, 50th, and 84th percentiles. The off-diagonal panel shows the two-dimensional posterior density estimated from the combined samples, with contours indicating enclosed-probability regions (we display 68.3%, 95.5%, and 99.7% levels). For comparison, the CCF injection-recovery curve (S/N = 5, yellow) is overplotted in the log₁₀(p_c) vs. log₁₀(Z/Z_⊙) panel. Panel B: marginalized one-dimensional posteriors for the β parameter separately for CARMENES (five transits) and CRIRES⁺ (one transit).

5.6 Detectability of H₂O with ANDES

The absence of molecular signatures in our CARMENES and CRIRES⁺ analyses motivated an assessment of what nextgeneration facilities might reveal for a planet like GJ 436 b. To this end, we generated a suite of simulated ELT/ANDES observations using EXoPLORE, an end-to-end high-resolution spectroscopic simulator, based on the methods described by Blain et al. (2024)⁸. This framework makes use of easyCHEM and petitRADTRANS for the chemical modelling of the atmosphere and radiative transfer computations, respectively. The tool produces synthetic, high-resolution time-series spectra that incorporate the planned YJH-bands wavelength coverage and S/N estimates derived from the ANDES exposure time calculator (ETC, integrations of 30 s, see Sanna et al. 2024), a resolving power of R = 100000, and a parametrised treatment of telluric evolution with airmass. The current implementation does not yet include additional sources of time-correlated noise nor the impact of rapid PWV fluctuations. Hence, the detectability values reported here should be interpreted as optimistic upper limits on the achievable performance.

For consistency with our observational analysis, we explored a two-dimensional grid in atmospheric metallicity and cloud-top pressure. For each grid point (log Z/Z_⊙, p_c), EXoPLORE generated a complete synthetic transit dataset that includes fundamental spectral contributions, namely:

We evaluated the telluric absorption using transmission curves from ESO’s Skycalc Tool (Noll et al. 2012; Jones et al. 2013), at the airmass of each synthetic exposure.
The modelled exoplanet contribution containing H₂, He, H₂O, CO₂, CO, CH₄, NH₃, and HCN was injected after applying the appropriate orbital Doppler shift (Eq. (3)) and was scaled by the instantaneous transit depth computed with BATMAN (Kreidberg 2015), adopting a uniform limb-darkening profile and the orbital geometry of GJ 436 b (Table 3).
The stellar spectrum was taken from the PHOENIX library (Husser et al. 2013), selecting a template consistent with the properties of GJ 436, T_eff = 3600 K, log g = 4.5, [Fe/H] = +0.5, and [a/M] = +0.2, closely following (within the database choices) the values reported in Mukherjee et al. (2025) and our own values (Table 2).
Spectral noise was randomly generated for each wavelength and frame, using Python’s functions np.random.default_rng(seed=12345) and rng.normal(), where the standard deviation was computed for each bin as σ(λ) = 1/(S/N)_ETC(λ).

The resulting in-silico time series of spectra were then processed through the same cross-correlation pipeline used for the real observations (Sect. 4.4) and the recovered signal at the expected K_P-υ_wind was recorded as the predicted detectability for that atmospheric scenario. Figure 10 presents the resulting detectability map for water vapour assuming a single transit of GJ 436 b, also comparing the in-silico predictions with the empirical results presented in this work.

The simulations show that the large collecting area of the ELT would provide a substantial boost in sensitivity relative to current facilities. The predicted significance is highest, as expected, in the region already discarded by the data analysed here, where faint cloud opacity (cloud tops at a few millibars or deeper) yield H₂O detections at S/N ≳ 10 from a single transit. The major improvement we observed came hence from the sensitivity to weak signals originating above clouds, as we would expect to detect H₂O in atmospheres for p_c > 0.1 mbar and 10 < Z < 100× solar at a S/N ∼ 5, which is almost an order-of-magnitude improvement over current instrumentation. However, the map reflects that atmospheres with extreme metal enrichment will remain challenging, since the reduced scale height suppresses the H₂O signal.

Therefore, our simulations indicate that ELT/ANDES will substantially enhance the recovery of molecular features in low mean-molecular-weight atmospheres, where the transmission signal is still in the photon-limited regime. In such cases, the increased collecting area boosts the visibility of line cores emerging above high-altitude clouds. For atmospheres with very high metallicities though, the compressed scale height still leads to weak line core absorption (see, also, Fig. 5). Thus, even the ELT would require multiple transits to secure a detection. This is also consistent with both the ELT simulation (one transit) and the analysed data combination (five CARMENES transits plus one CRIRES⁺ dataset) yielding similar significances at intermediate and high-metallicities (i.e. at the signal-limited regime). Nevertheless, the ability to probe high-metallicity atmospheres by combining multiple ELT datasets represents a new avenue for probing these exoplanets, one that is not reachable with current facilities.

Fig. 10

Simulated H₂O detectability map for GJ 436 b produced with EXoPLORE assuming a single transit observed with ANDES. The white and red contours indicate recovered S/N levels of 10 and 5 with in-silico data, respectively. To illustrate the improvement over current constraints, we show the S/N = 5 contour (orange curve) derived from combining the CARMENES and CRIRES⁺ data (Fig. 8). Although based on simulations, the map illustrates the significant improvement in sensitivity that ANDES and the ELT are expected to provide for probing cloudy sub-Neptune atmospheres.

6 Conclusions

We present the analysis of six transit datasets of the warm Neptune GJ 436 b observed with the high-resolution spectrographs CARMENES (five events) and CRIRES⁺ (one event) at NIR wavelengths. After applying the CCF technique, performing injection recovery tests, and running Bayesian retrievals, we were able to establish the most stringent constraints on the detectability of H₂O, CH₄, and CO in the atmosphere of GJ 436 b from ground-based observations to date. The main conclusions derived from our work are:

Our CCF studies of the CARMENES and CRIRES⁺ datasets did not recover any significant signals for the studied species in this exoplanet’s atmosphere, which is in line with findings by Knutson et al. (2014a), Grasser et al. (2024), and Mukherjee et al. (2025). Nevertheless, molecular species such as H₂O, CH₄, NH₃, and CO are expected in GJ 436 b’s atmosphere, which is thought to be metal rich (Moses et al. 2013; Knutson et al. 2014a);
We performed CCF-based injection recovery studies so as to better assess the detectability of atmospheric compounds with our combination of instruments, nights, and techniques. We ruled out atmospheric scenarios for GJ 436 b with low-altitude (p_c < 1 mbar) or no cloud deck coverage and medium to low metallicities (Z < 600× solar), as such atmospheres would have been detectable in our datasets. These values agree well with the upper limits derived by models (Moses et al. 2013) and observations at high and low spectral resolution (Knutson et al. 2014a; Grasser et al. 2024; Mukherjee et al. 2025). The most likely explanation for the lack of observable spectral features in the atmosphere of GJ 436 b is the presence of a cloud deck at higher altitudes than < 1 mbar pressure level, and/or a very high atmospheric metallicity, beyond 600× solar;
A joint analysis of five CARMENES transits of GJ436b provides tighter constraints for H₂O and CH₄ than a single CRIRES⁺ observation, essentially due to the broader coverage of the former instrument. For CO, however, both datasets yield similar (high) constraints since they cover the same isolated CO band. Restricting the CARMENES orders to match CRIRES⁺ H1567 coverage resulted in CARMENES upper limits being similar, albeit slightly less restrictive, to those obtained with CRIRES⁺. This occurs especially at the lowest and highest metallicities explored (i.e. weaker H₂O lines), potentially due to the higher resolving power;
The tightest CCF-based upper limits to date in the atmosphere of GJ 436 b are obtained by combining all available CARMENES and CRIRES⁺ datasets and performing injection recovery tests of H2O, the strongest absorber in the NIR range covered by both CARMENES and CRIRES⁺. This analysis determines that the atmosphere presents a cloud deck at altitudes higher than the 1 mbar level (regardless of the atmospheric metallicity), higher than the 10 mbar level for a solar metallicity, or presenting extremely high metallicities (≳ 900× solar), regardless of cloud presence;
Bayesian retrievals of the six transit datasets confirm that no prominent atmospheric features are detectable in GJ 436 b. The combined posteriors show that only low-pressure (highaltitude) cloud decks or very high-metallicity atmospheres are consistent with the data. The regions allowed by the retrieval fall outside the CCF-detectable domain, which means that the data effectively rule out models that would produce detectable spectral signals. Therefore, these two metrics provide a mutually consistent picture of an essentially featureless transmission spectrum;
Our ELT/ANDES simulations indicate that the leap to 40 m-class telescopes will significantly expand the atmospheric parameter space accessible for sub-Neptune characterization, even in the presence of clouds. For GJ 436 b, molecular features arising above cloud decks at pressures of 0.1 −1 mbar could be detectable in a single transit across a broad metallicity range (10 < Z < 300× solar), which is unattainable with present-day instruments. At the same time, atmospheres with high metal enrichment (Z > 400× solar) may remain intrinsically difficult targets due to their compressed scale heights, and meaningful detections in that regime are likely to require the combination of multiple ELT transits. Taken together, these results highlight both the new opportunities and the persistent physical limits that will shape atmospheric studies in the ELT era;
For completeness, we derived precise and homogeneous stellar and planetary parameters from TESS, CARMENES, HIRES, HARPS, and ESPRESSO data, which are consistent with previous determinations but add extra robustness to our work.

These findings confirm GJ 436 b among the class of subNeptunes with intrinsically challenging transmission spectra. Looking ahead, our ELT/ANDES simulations demonstrate that the forthcoming generation of 40 m-class facilities will open up a substantial new discovery space and potentially enable routine access to molecular features that remain beyond the reach of current instruments. As such, GJ 436 b will become a valuable benchmark for testing the capabilities of next-generation observatories and for deepening our understanding of the atmospheric diversity of warm Neptunes.

Acknowledgements

We thank the anonymous referee for their very helpful comments on this manuscript. We also thank Natalie Grasser for kindly sharing her results with us for comparison purposes. CARMENES is an instrument at the Centro Astronómico Hispano en Andalucía (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). Funding for CARMENES has been provided by the Max-Planck-Institut für Astronomie (MPIA), Consejo Superior de Investigaciones Científicas (CSIC), European Regional Development Fund (ERDF), Ministerio de Ciencia, Innovación y Universidades (MICIU), Deutsche Forschungsgemeinschaft (DFG) and the members of the CARMENES Consortium (https://carmenes.caha.es). We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the MICIU and the ERDF “A way of making Europe” through projects PID2022-141216NB-I00, PID2022-137241NB-C4[1:4], PID2021-125627OB-C3[1:2], RYC2022-037854-I, and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Andalucía (CEX2021-001131-S) and Institut de Ciències de l’Espai (CEX2020-001058-M). This work was also funded by the Centres de Recerca de Catalunya (CERCA), a Fraunhofer-Schwarzschild fellowship at the Ludwig-Maximilians-Universität München and DFG under Germany’s Excellence Strategy (EXC 2094 390783311), and the European Union through European Research Council (ERC) Consolidator grant EvapourA-TOR (No. 101170037). Views and opinions expressed are those of the author only and do not necessarily reflect those of the European Union or the ERC. Neither the European Union nor the granting authority can be held responsible for them.

References

Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5 [NASA ADS] [CrossRef] [Google Scholar]
Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Alonso-Floriano, F. J., Sánchez-López, A., Snellen, I. A. G., et al. 2019, A&A, 621, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Antoniadis-Karnavas, A., Sousa, S. G., Delgado-Mena, E., Santos, N. C., & Andreasen, D. T. 2024, A&A, 690, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [CrossRef] [Google Scholar]
Barclay, T., Kostov, V. B., Colón, K. D., et al. 2021, AJ, 162, 300 [NASA ADS] [CrossRef] [Google Scholar]
Basilicata, M., Giacobbe, P., Bonomo, A. S., et al. 2024, A&A, 686, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bean, J. L., Miller-Ricci Kempton, E., & Homeier, D. 2010, Nature, 468, 669 [NASA ADS] [CrossRef] [Google Scholar]
Beatty, T. G., Welbanks, L., Schlawin, E., et al. 2024, ApJ, 970, L10 [Google Scholar]
Beaulieu, J. P., Tinetti, G., Kipping, D. M., et al. 2011, ApJ, 731, 16 [Google Scholar]
Benneke, B., Knutson, H. A., Lothringer, J., et al. 2019a, Nat. Astron., 3, 813 [Google Scholar]
Benneke, B., Wong, I., Piaulet, C., et al. 2019b, ApJ, 887, L14 [Google Scholar]
Benneke, B., Roy, P.-A., Coulombe, L.-P., et al. 2024, arXiv e-prints [arXiv:2403.03325] [Google Scholar]
Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, ApJ, 747, 35 [NASA ADS] [CrossRef] [Google Scholar]
Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35 [Google Scholar]
Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017, AJ, 153, 138 [NASA ADS] [CrossRef] [Google Scholar]
Blain, D., Sánchez-López, A., & Mollière, P. 2024, AJ, 167, 179 [NASA ADS] [CrossRef] [Google Scholar]
Boucher, A., Lafreniére, D., Pelletier, S., et al. 2023, MNRAS, 522, 5062 [NASA ADS] [CrossRef] [Google Scholar]
Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bourrier, V., Zapatero Osorio, M. R., Allart, R., et al. 2022, A&A, 663, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Brande, J., Crossfield, I. J. M., Kreidberg, L., et al. 2022, AJ, 164, 197 [CrossRef] [Google Scholar]
Brogi, M., & Line, M. R. 2019, AJ, 157, 114 [Google Scholar]
Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502 [Google Scholar]
Brogi, M., Line, M., Bean, J., Désert, J. M., & Schwarz, H. 2017, ApJ, 839, L2 [Google Scholar]
Brogi, M., Giacobbe, P., Guilluy, G., et al. 2018, A&A, 615, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Buchner, J. 2016, PyMultiNest: Python interface for MultiNest, Astrophysics Source Code Library [record ascl:1606.005] [Google Scholar]
Butler, R. P., Vogt, S. S., Marcy, G. W., et al. 2004, ApJ, 617, 580 [NASA ADS] [CrossRef] [Google Scholar]
Caballero, J. A., Guàrdia, J., López del Fresno, M., et al. 2016, SPIE Conf. Ser., 9910, 99100E [Google Scholar]
Cabot, S. H. C., Madhusudhan, N., Hawker, G. A., & Gandhi, S. 2019, MNRAS, 482, 4422 [NASA ADS] [CrossRef] [Google Scholar]
Casasayas-Barris, N., Orell-Miquel, J., Stangret, M., et al. 2021, A&A, 654, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Chachan, Y., Knutson, H. A., Gao, P., et al. 2019, AJ, 158, 244 [Google Scholar]
Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377 [Google Scholar]
Charbonneau, D., Berta, Z. K., Irwin, J., et al. 2009, Nature, 462, 891 [NASA ADS] [CrossRef] [Google Scholar]
Charnay, B., Blain, D., Bézard, B., et al. 2021, A&A, 646, A171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cheverall, C. J., Madhusudhan, N., & Holmberg, M. 2023, MNRAS, 522, 661 [NASA ADS] [CrossRef] [Google Scholar]
Chubb, K. L., Min, M., Kawashima, Y., Helling, C., & Waldmann, I. 2020, A&A, 639, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cifuentes, C., Caballero, J. A., Cortés-Contreras, M., et al. 2020, A&A, 642, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Claudi, R., Benatti, S., Carleo, I., et al. 2016, Proc. SPIE, 9908, 99081A [Google Scholar]
Cont, D., Yan, F., Reiners, A., et al. 2021, A&A, 651, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cont, D., Yan, F., Reiners, A., et al. 2022, A&A, 668, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cont, D., Nortmann, L., Yan, F., et al. 2024, A&A, 688, A206 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cont, D., Nortmann, L., Lesjak, F., et al. 2025, A&A, 698, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cristofari, P. I., Donati, J. F., Masseron, T., et al. 2022, MNRAS, 516, 3802 [NASA ADS] [CrossRef] [Google Scholar]
Dash, S., Brogi, M., Gandhi, S., et al. 2024, MNRAS, 530, 3100 [Google Scholar]
de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Deibert, E. K., de Mooij, E. J. W., Jayawardhana, R., et al. 2021, AJ, 161, 209 [NASA ADS] [CrossRef] [Google Scholar]
Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016, MNRAS, 455, 2018 [NASA ADS] [CrossRef] [Google Scholar]
Dorn, R. J., Anglada-Escude, G., Baade, D., et al. 2014, The Messenger, 156, 7 [NASA ADS] [Google Scholar]
Dorn, R. J., Bristow, P., Smoker, J. V., et al. 2023, A&A, 671, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [Google Scholar]
Evans, T. M., Sing, D. K., Wakeford, H. R., et al. 2016, ApJ, 822, L4 [NASA ADS] [CrossRef] [Google Scholar]
Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2013, arXiv e-prints [arXiv:1306.2144] [Google Scholar]
Figueira, P., Pont, F., Mordasini, C., et al. 2009, A&A, 493, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Flowers, E., Brogi, M., Rauscher, E., Kempton, E. M.-R., & Chiavassa, A. 2019, AJ, 157, 209 [Google Scholar]
Follert, R., Dorn, R. J., Oliva, E., et al. 2014, SPIE Conf. Ser., 9147, 914719 [Google Scholar]
Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gandhi, S., Brogi, M., & Webb, R. K. 2020, MNRAS, 498, 194 [Google Scholar]
Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485 [Google Scholar]
Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205 [NASA ADS] [CrossRef] [Google Scholar]
Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215 [Google Scholar]
Gibson, N. P., Nugroho, S. K., Lothringer, J., Maguire, C., & Sing, D. K. 2022, MNRAS, 512, 4618 [NASA ADS] [CrossRef] [Google Scholar]
Gillon, M., Pont, F., Demory, B. O., et al. 2007, A&A, 472, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Grasser, N., Snellen, I. A. G., Landman, R., Picos, D. G., & Gandhi, S. 2024, A&A, 688, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guilluy, G., Gressier, A., Wright, S., et al. 2021, AJ, 161, 19 [Google Scholar]
Guilluy, G., D’Arpa, M. C., Bonomo, A. S., et al. 2024, A&A, 686, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guo, B., Yan, F., Nortmann, L., et al. 2024, A&A, 687, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hargreaves, R. J., Gordon, I. E., Rey, M., et al. 2020, ApJS, 247, 55 [NASA ADS] [CrossRef] [Google Scholar]
Hawker, G. A., Madhusudhan, N., Cabot, S. H. C., & Gandhi, S. 2018, ApJ, 863, L11 [NASA ADS] [CrossRef] [Google Scholar]
Herman, M. K., de Mooij, E. J. W., Jayawardhana, R., & Brogi, M. 2020, AJ, 160, 93 [Google Scholar]
Herman, M. K., de Mooij, E. J. W., Nugroho, S. K., Gibson, N. P., & Jayawardhana, R. 2022, AJ, 163, 248 [NASA ADS] [CrossRef] [Google Scholar]
Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, Astron. Astrophys., 553, A6 [Google Scholar]
Irwin, P. G. J. 2009, Giant Planets of Our Solar System (Springer) [Google Scholar]
Jenkins, J. M., Chandrasekaran, H., McCauliff, S. D., et al. 2010, SPIE Conf. Ser., 7740, 77400D [Google Scholar]
Jones, A., Noll, S., Kausch, W., Szyszka, C., & Kimeswenger, S. 2013, A&A, 560, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kaeufl, H.-U., Ballester, P., Biereichel, P., et al. 2004, SPIE Conf. Ser., 5492, 1218 [NASA ADS] [Google Scholar]
Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kempton, E. M. R., Zhang, M., Bean, J. L., et al. 2023, Nature, 620, 67 [NASA ADS] [CrossRef] [Google Scholar]
Kesseli, A. Y., & Snellen, I. A. G. 2021, ApJ, 908, L17 [Google Scholar]
Kesseli, A. Y., Snellen, I. A. G., Alonso-Floriano, F. J., Mollière, P., & Serindag, D. B. 2020, AJ, 160, 228 [NASA ADS] [CrossRef] [Google Scholar]
Knutson, H. A., Benneke, B., Deming, D., & Homeier, D. 2014a, Nature, 505, 66 [Google Scholar]
Knutson, H. A., Madhusudhan, N., Cowan, N. B., et al. 2011, ApJ, 735, 27 [NASA ADS] [CrossRef] [Google Scholar]
Knutson, H. A., Dragomir, D., Kreidberg, L., et al. 2014b, ApJ, 794, 155 [Google Scholar]
Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69 [Google Scholar]
Kreidberg, L., Mollière, P., Crossfield, I. J. M., et al. 2022, AJ, 164, 124 [CrossRef] [Google Scholar]
Lanotte, A. A., Gillon, M., Demory, B. O., et al. 2014, A&A, 572, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lei, E., & Mollière, P. 2025, J. Open Source Softw., 10, 7712 [Google Scholar]
Lesjak, F., Nortmann, L., Yan, F., et al. 2023, A&A, 678, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506 [NASA ADS] [CrossRef] [Google Scholar]
Lockwood, A. C., Johnson, J. A., Bender, C. F., et al. 2014, ApJ, 783, L29 [Google Scholar]
Lothringer, J. D., Benneke, B., Crossfield, I. J. M., et al. 2018, AJ, 155, 66 [NASA ADS] [CrossRef] [Google Scholar]
Loyd, R. O. P., Schneider, P. C., Jackman, J. A. G., et al. 2023, AJ, 165, 146 [NASA ADS] [CrossRef] [Google Scholar]
Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24 [NASA ADS] [CrossRef] [Google Scholar]
Madhusudhan, N., & Seager, S. 2011, ApJ, 729, 41 [Google Scholar]
Madhusudhan, N., Crouzet, N., McCullough, P. R., Deming, D., & Hedges, C. 2014, ApJ, 791, L9 [Google Scholar]
Madhusudhan, N., Piette, A. A. A., & Constantinou, S. 2021, ApJ, 918, 1 [NASA ADS] [CrossRef] [Google Scholar]
Madhusudhan, N., Sarkar, S., Constantinou, S., et al. 2023, ApJ, 956, L13 [NASA ADS] [CrossRef] [Google Scholar]
Maguire, C., Sedaghati, E., Gibson, N. P., Smette, A., & Pino, L. 2024, A&A, 692, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Marconi, A., Abreu, M., Adibekyan, V., et al. 2022, in Ground-based and Airborne Instrumentation for Astronomy IX, 12184, SPIE, 720 [Google Scholar]
Marfil, E., Tabernero, H. M., Montes, D., et al. 2021a, A&A, 656, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Marfil, E., Tabernero, H. M., Montes, D., et al. 2021b, VizieR Online Data Catalog: J/A+A/656/A162 [Google Scholar]
Marley, M. S., Ackerman, A. S., Cuzzi, J. N., & Kitzmann, D. 2013, in Comparative Climatology of Terrestrial Planets, eds. S. J. Mackwell, A. A. Simon-Miller, J. W. Harder, & M. A. Bullock, 367 [Google Scholar]
Masson, A., Vinatier, S., Bézard, B., et al. 2024, A&A, 688, A179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77 [NASA ADS] [CrossRef] [Google Scholar]
Mazeh, T., Tamuz, O., & Zucker, S. 2007, in Astronomical Society of the Pacific Conference Series, 366, Transiting Extrapolar Planets Workshop, eds. C. Afonso, D. Weldrake, & T. Henning, 119 [Google Scholar]
Melo, E., Souto, D., Cunha, K., et al. 2024, ApJ, 973, 90 [Google Scholar]
Miller-Ricci Kempton, E., & Rauscher, E. 2012, ApJ, 751, 117 [CrossRef] [Google Scholar]
Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67 [Google Scholar]
Morello, G., Waldmann, I. P., Tinetti, G., et al. 2015, ApJ, 802, 117 [NASA ADS] [CrossRef] [Google Scholar]
Moses, J. I., Line, M. R., Visscher, C., et al. 2013, ApJ, 777, 34 [Google Scholar]
Mukherjee, S., Schlawin, E., Bell, T. J., et al. 2025, ApJ, 982, L39 [Google Scholar]
Nettelmann, N., Kramm, U., Redmer, R., & Neuhäuser, R. 2010, A&A, 523, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Noll, S., Kausch, W., Barden, M., et al. 2012, A&A, 543, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388 [Google Scholar]
Nortmann, L., Lesjak, F., Yan, F., et al. 2025, A&A, 693, A213 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, ApJ, 898, L31 [Google Scholar]
Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJ, 910, L9 [CrossRef] [Google Scholar]
Ohno, K., Schlawin, E., Bell, T. J., et al. 2025, ApJ, 979, L7 [Google Scholar]
Palle, E., Biazzo, K., Bolmont, E., et al. 2025, Exp. Astron., 59, 29 [Google Scholar]
Parker, L. T., Mendonça, J. M., Diamond-Lowe, H., et al. 2025, MNRAS, 538, 3263 [Google Scholar]
Passegger, V. M., Bello-García, A., Ordieres-Meré, J., et al. 2022, A&A, 658, A194 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pinhas, A., Rackham, B. V., Madhusudhan, N., & Apai, D. 2018, MNRAS, 480, 5314 [Google Scholar]
Pino, L., Ehrenreich, D., Allart, R., et al. 2018a, A&A, 619, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pino, L., Ehrenreich, D., Wyttenbach, A., et al. 2018b, A&A, 612, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, SPIE Conf. Ser., 9147, 91471F [Google Scholar]
Quirrenbach, A., Amado, P. J., Ribas, I., et al. 2018, SPIE Conf. Ser., 10702, 107020W [Google Scholar]
Rafi, S. A., Nugroho, S. K., Tamura, M., Nortmann, L., & Sánchez-López, A. 2024, AJ, 168, 106 [NASA ADS] [CrossRef] [Google Scholar]
Ramstad, R., & Barabash, S. 2021, Space Sci. Rev., 217, 36 [Google Scholar]
Revilla, D., Amado, P. J., Luque, R., et al. 2025, Science, accepted [Google Scholar]
Ridden-Harper, A., Nugroho, S. K., Flagg, L., et al. 2023, AJ, 165, 170 [NASA ADS] [CrossRef] [Google Scholar]
Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
Rumenskikh, M. S., Khodachenko, M. L., Shaikhislamov, I. F., et al. 2023, MNRAS, 526, 4120 [NASA ADS] [Google Scholar]
Sánchez-López, A., Alonso-Floriano, F. J., López-Puertas, M., et al. 2019, A&A, 630, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sánchez-López, A., López-Puertas, M., Snellen, I. A. G., et al. 2020, A&A, 643, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sánchez-López, A., Landman, R., Molliere, P., et al. 2022, in Bulletin of the American Astronomical Society, 54, 108.03 [Google Scholar]
Sanna, N., Martins, B. C., Martins, A. d. M., et al. 2024, in Ground-based and Airborne Instrumentation for Astronomy X, 13096, SPIE, 1299 [Google Scholar]
Schlawin, E., Ohno, K., Bell, T. J., et al. 2024, ApJ, 974, L33 [Google Scholar]
Schweitzer, A., Passegger, V. M., Cifuentes, C., et al. 2019, A&A, 625, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Seager, S. 2013, Science, 340, 577 [NASA ADS] [CrossRef] [Google Scholar]
Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Smith, P. C. B., Line, M. R., Bean, J. L., et al. 2024, AJ, 167, 110 [NASA ADS] [CrossRef] [Google Scholar]
Snellen, I. A. G. 2025, ARA&A, 63, 83 [Google Scholar]
Snellen, I. A. G., Albrecht, S., de Mooij, E. J. W., & Le Poole, R. S. 2008, A&A, 487, 357 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049 [Google Scholar]
Sousa-Silva, C., Hesketh, N., Yurchenko, S. N., Hill, C., & Tennyson, J. 2014, J. Quant. Spec. Radiat. Transf., 142, 66 [Google Scholar]
Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2022, A&A, 662, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Stevenson, K. B., Harrington, J., Nymeyer, S., et al. 2010, Nature, 464, 1161 [NASA ADS] [CrossRef] [Google Scholar]
Tabernero, H. M., Marfil, E., Montes, D., & González Hernández, J. I. 2022, A&A, 657, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466 [Google Scholar]
Trifonov, T., Kürster, M., Zechmeister, M., et al. 2018, A&A, 609, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502 [Google Scholar]
Vidotto, A. A., Bourrier, V., Fares, R., et al. 2023, A&A, 678, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wakeford, H. R., Sing, D. K., Deming, D., et al. 2018, AJ, 155, 29 [Google Scholar]
Yan, F., Reiners, A., Pallé, E., et al. 2022, A&A, 659, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

¹

https://archive.eso.org/cms/eso-data.html

²

https://ftp.eso.org/pub/dfs/pipelines/instruments/cr2res/cr2re-pipeline-manual-1.6.10.pdf

³

https://github.com/hmtabernero/SteParSyn/

⁴

https://tess.mit.edu/toi-releases/

⁵

https://github.com/jlillo/tpfplotter

⁶

https://petitradtrans.readthedocs.io

⁷

https://easychem.readthedocs.io

⁸

A separated publication and an open-access framework for EXoPLORE are in development.

Appendix A Additional plots

Fig. A.1

Spectral coverage comparison of CH₄ (top) and CO (bottom) transmission spectrum models for GJ 436 b, as observed by the CARMENES NIR channel (grey) and the CRIRES⁺ H1567 setting used by Grasser et al. (2024, red).

Fig. A.2

TPF of GJ 436 (cross symbol) in TESS Sectors 22 (left) and 49 (right). The electron counts are colour-coded. The TESS optimal photometric aperture per sector used to obtain the SAP fluxes is marked with red squares. The Gaia DR3 sources with G-band magnitudes down to 8 mag fainter than GJ 436 are labeled with numbers (GJ 436 corresponds to number 1) and their scaled brightness based on Gaia magnitudes is shown by red circles of different sizes (see figure inset). The pixel scale is 21 arcsec pixel⁻¹.

Fig. A.3

TESS light curves folded in phase with the orbital periods of GJ436 b. Binned data are plotted as white circles. The best transit model is shown with a black line.

Fig. A.4

Posterior distributions of the principal fitted planetary parameters of the GJ 436 system as obtained from the joint transit and RV fit, assuming eccentric orbit. The (black) vertical dashed lines indicate the 16% and 84% quantiles that were used to define the optimal values and their associated 1σ uncertainty. The blue lines stands for the median values (50% quantile) of each fitted parameter.

Fig. A.5

GJ436 HARPS, HIRES, ESPRESSO, and CARMENES RVs of GJ 436 (coloured dots) and the best model (black line) from the joint photometric and spectroscopic fit folded in phase with the orbital period of GJ 436 b, assuming eccentricity in the model.

Fig. A.6

Percentage of masked pixels per spectral order for each observation night of the CARMENES (left) and CRIRES⁺ (right) datasets. Higher percentages correspond to spectral orders containing the strongest telluric absorption bands, primarily due to water vapour.

Fig. A.7

Illustration of the selection criterion used for order-wise SysRem optimization for the CARMENES datasets. Here, the values of ∆σ (vertical axis) are plotted as a function of the SysRem pass (horizontal axis) across all spectral orders for observations from all five nights. Red stars mark the SysRem pass at which a plateau is reached, indicating the point at which the algorithm is halted. Spectral order number is indicated in the upper right corner of each cell. For clarity, an offset was added to the per-night curves, as indicated in the legend.

Fig. A.8

Same as Fig. A.7 but for the CRIRES⁺ dataset.

Fig. A.9

Same as Fig. 7 but for each individual CARMENES transit analysed and for the molecules studied. First row corresponds to water vapour, second row to methane, and third to carbon monoxide.

Fig. A.10

Cross-correlation maps in S/N units for potential atmospheric signals as a function of v_rest and K_p for the primary expected absorbing species in the atmosphere of GJ 436 b. For CARMENES (five combined transits; top row) and CRIRES⁺ (one transit; bottom row) we explored, water vapour (left), methane (middle), and carbon monoxide (right, only for CARMENES). Horizontal and vertical red lines indicate the expected K_p and v_rest eccentric orbital models. These example maps were derived using a 10× solar metallicity template with a cloud deck at 10 mbar. No molecular signals were detected in our cross-correlation analyses.

Fig. A.11

S/N slices from Fig. 8 showing the S/N (vertical axis) as a function of the injected-model metallicity (horizontal axis). In the top panel, we show this result for CH₄, and CO is shown in the bottom panel. Different colours correspond to the different fixed p_c. Dashed lines correspond to S/N slices at the same fixed p_c from Fig. 6 of Grasser et al. (2024). Due to the unavailability of the original data, these values were extracted using WebPlotDigitizer (Version 5.2).

Fig. A.12

Same as Fig. 8 but for H₂O (left) and CH₄ (right) but limiting the CARMENES spectral range to closely match that of CRIRES⁺.

All Tables