The GAPS Programme at TNG XXXVIII.Five molecules in the atmosphere of the warm giant planet WASP-69b detected at high spectral resolution

Context.Aims. The ﬁeld of exo-atmospheric characterisation is progressing at an extraordinary pace. Atmospheric observations are now available for tens of exoplanets, mainly hot and warm inﬂated gas giants, and new molecular species continue to be detected revealing a richer atmospheric composition than previously expected. Thanks to its warm equilibrium temperature (963 ± 18 K) and low-density (0.219 ± 0.031 g cm − 3 ), the close-in gas giant WASP-69b represents a golden target for atmospheric characterisation. Methods. With the aim of searching for molecules in the atmosphere of WASP-69b and investigating its properties, we performed high-resolution transmission spectroscopy with the GIANO-B near-infrared spectrograph at the Telescopio Nazionale Galileo. Results. We observed three transit events of WASP-69b. During a transit, the planetary lines are Doppler-shifted due to the large change in the planet’s radial velocity, allowing us to separate the planetary signal from the quasi-stationary telluric and stellar spectrum. Conclusions. Considering the three nights together, we report the detection of CH 4 , NH 3 , CO, C 2 H 2 , and H 2 O, at more than 3 . 3 σ level. We did not identify the presence of HCN and CO 2 with conﬁdence level higher than 3 σ . This is the ﬁrst time that ﬁve molecules are simultaneously detected in the atmosphere of a warm giant planet. These results suggest that the atmosphere of WASP-69b is possibly carbon-rich and characterised by the presence of disequilibrium chemistry.


Introduction
Ground-based, high-resolution (resolving power R ≥20 000) spectroscopy (HRS) in the near-infrared (NIR) represents an effective approach to investigate exoplanet atmospheres. Indeed, at high spectral resolution, thousands of resolved molecular lines can be disentangled, as they are Doppler shifted by tens of km s −1 due to the planet motion along its orbit from the telluric (Earth's atmosphere) and stellar lines that are (quasi-)stationary in wavelength.
Hot (T eq > 1000 K) and warm (T eq ≤ 1000 K) inflated giant planets are the most suitable targets for atmospheric characterisation (Seager & Deming 2010). Recently, Giacobbe et al. (2021) reported the simultaneous detection of six molecules in the atmosphere of one of the best-studied hot-Jupiters (HJs), HD 209458b. As only two molecules (CO and H 2 O) had been detected at the same time in exoplanetary atmospheres in the past (e.g. Tsiaras et al. 2018), the work presented in Giacobbe et al. (2021) revealed a previously unknown chemical complexity. Thus this discovery raises the question of whether the complexity of HD 209458b's atmosphere is unique or if other Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundacion Galileo Galilei of INAF at the Spanish Observatorio Roque de los Muchachos of the IAC in the frame of the programme Global Architecture of the Planetary Systems (GAPS). exo-atmospheres show such a rich molecular composition. Such molecular diversity enables one to derive atmospheric abundances and thus elemental ratios, such as the carbon-to-oxygen (C/O) ratio, can then be obtained, allowing us to derive important clues on the formation and migration histories of hot giant planets (e.g. Madhusudhan 2012;Turrini et al. 2021b). For example, Giacobbe et al. (2021) retrieved a C/O ratio ≥ 1, meaning that HD 209458b formed far from its present location and subsequently migrated inwards.
With a warm temperature and low density (see Table 1), the close-in transiting planet WASP-69b represents a suitable laboratory for atmospheric characterisation studies. Furthermore, the low equilibrium temperature T eq ≤1000 K places the planet in a different regime compared to that of the HJ HD 209458b (T eq ∼1500 K): the equilibrium chemistry shifts in favour of CH 4 rather than CO and disequilibrium phenomena (i.e. vertical mixing and photochemical processes) are more likely to alter the molecular abundances compared to the atmospheres of hotter planets (e.g. Moses 2014).
The atmosphere of WASP-69b has been studied with both space-based low-resolution (LR) and ground-based highresolution (HR) in-transit (transmission spectroscopy) observations. LR observations performed with WFC3 on board the Hubble Space Telescope (HST) unveiled the presence of H 2 O and possibly aerosols (Tsiaras et al. 2018;Fisher & Heng 2018). Sodium has also been detected through HR optical observations (e.g. Casasayas- Barris et al. 2017). Moreover, HR and LR observations covering the NIR He i (1083.3 nm) revealed that the planet has an extended and escaping atmosphere, where the escape is driven by the intense X-ray and extreme ultraviolet radiation received from the active K5V host star (e.g. Nortmann et al. 2018;Vissapragada et al. 2020).
In this work, we present transmission spectroscopy of WASP-69b using HR observations in the NIR, searching for the absorption of multiple molecular species. We observed the WASP-69 system within the long-term observing programme at the Telescopio Nazionale Galileo (TNG) telescope 'GAPS2: the origin of planetary systems' -awarded to the Italian Global Architecture of Planetary System (GAPS) Collaboration (Poretti et al. 2016). Part of this programme is devoted to the characterisation of exoplanetary atmospheres at HR. Briefly, we are carrying out transmission and emission spectroscopy of twentysix hot Jupiters and four hot Neptunes and/or sub-Neptunes in a wide range of T eq (625-2500 K), and hosted by stars of different spectral types (A to M with K < 10.5 mag), to accomplish the following: (i) identify atomic species to estimate the abundance of refractory elements (Pino et al. 2020) and/or to probe atmospheric expansion or escape, for instance from the Hα ) and/or metastable helium triplet (Guilluy et al. 2019); (ii) detect molecular compounds and estimate the atmospheric C/O ratio and metallicity, which in turn provide constraints on the planet formation and evolution scenarios (Giacobbe et al. 2021); (iii) derive temperature-pressure profiles (Pino et al. 2020;Borsa et al. 2022); and (iv) detect the atmospheric Rossiter-McLaughlin effect (Borsa et al. 2019;Rainer et al. 2021). Fig. 1 places WASP-69b into the broader context of the GAPS atmospheric survey. The symbol size is proportional to the expected atmospheric signal in the K band: where we account for the stellar magnitude in K band (F = 10 − magK 2.5 ), the atmospheric scale height (H), the planetary (R P ), and the stellar radius (R ). The most suitable targets to perform atmospheric characterisation studies are labelled. As this plot shows, WASP-69b represents one of the most suitable candidates to perform atmospheric studies (2R P H/R 2 ∼ 283 ppm 1 , Brown et al. 2001 ) in the temperature range (≤1000 K) we aimed to probe within this study. This paper is organised as follows: in Sect. 2 we present the data we collected and in Sect. 3 we describe the analysis that we performed to extract the planetary transmission spectrum and, therefore, the molecules responsible for absorption from the raw GIANO-B data. Finally, we discuss our results and draw the conclusions in Sect. 4.

Data sample
We observed the WASP-69 system with the TNG telescope in GIARPS configuration by simultaneously measuring HR spectra in the optical (0.39-0.69 µm) and NIR (0.95-2.45 µm) with the HARPS-N (R ≈ 115 000) and GIANO-B (R∼50,000) spectrographs. Our data encompass a total of three primary transits of WASP-69 b, which were obtained at UT 24 July 2019, UT 09 August 2020, and UT 28 October 2021. The GIANO-B observations were taken according to an ABAB nodding pattern   Table 1) are marked with vertical dashed lines. Black stars on the last night (28 October 2021) indicate the discarded AB couples for low S/N. (Claudi et al. 2017) for an optimal subtraction of the thermal background noise and telluric emission lines. The observations consist of spectra taken before, during, and after the planetary transit; Fig. 2 shows the signal-to-noise ratio (S/N) averaged over the entire GIANO-B spectral coverage as a function of the planet's orbital phase and of the airmass for each of the three nights. The observations during the last night (28 October 2021) were affected by several thin clouds (cirri), so we decided to discard the AB couples of observations that exhibit a very low S/N compared to the others if the lowest S/N in the couple is less than 15, see Fig. 2). A log of the observations is reported in Table A.1.  (40) 1,2 * Derived from a, P, and i as 2πa P sin i.

Data analysis
Here we report the steps we performed in order to recover the planetary signal from the GIANO-B raw images.
(i) Basic data reduction. We processed the raw GIANO-B spectra with the GOFIO pipeline (Rainer et al. 2018) to perform flat-fielding, bad pixel correction, background subtraction, optimal extraction of the 1D spectra, and a preliminary wavelength calibration using a uranium-neon lamp. We corrected small temporal variations in the wavelength solution due to the imperfect stability of the spectrograph (Giacobbe et al. 2021), and we refined the initial wavelength solution by employing a HR transmission spectrum of the Earth's atmosphere generated via the ESO Sky Model Calculator 2 . As some spectral regions have too many or too few telluric lines, this calibration refinement was not possible for all the orders (e.g. Brogi et al. 2018;Guilluy et al. 2019;Fossati et al. 2022). In Table A.2, we report the discarded orders for each individual night of observation.
(ii) Telluric and stellar lines' removal. We then removed the telluric contamination and the stellar spectrum by employing the principal component analysis (PCA) and linear regression techniques (Giacobbe et al. 2021). We selected the optimal number of components N opt following the procedure described in Giacobbe et al. (2021). More precisely, N opt varies with the considered order depending on the root mean square (rms) variation. We selected the number for which the first derivative between the last and penultimate component decreases by less than σ white N −1/2 opt , where σ white is the standard deviation of the full matrix assuming pure white noise. Since in two nights out of three (i.e. UT 24 July 2019 and UT 09 August 2020) the relative velocity shift between telluric and planetary lines was close to zero, we masked the telluric components to avoid contamination in our final spectra. Following Brogi et al. (2018), we thus masked the data assigning zero flux to those spectral channels corresponding to telluric lines deeper than 20%. Finally, we visually inspected the outcome for spurious telluric residuals. (iii) Search for molecules through cross-correlation. For each considered molecule, we combined the contribution of thousands of lines in the spectral coverage of GIANO-B by performing the cross-correlation (CC) of the residual spectra with isothermal transmission models generated using the GENE-SIS code (Gandhi & Madhusudhan 2017) adapted for transmission spectroscopy (Pinhas et al. 2018) (see Table 2 for the adopted line list). Prior to CC, the models were convolved to the GIANO-B instrument profile (a Gaussian profile with FWHM 5.4 km s −1 ). Models were calculated between the 100 and 10 −8 bar in pressure, at a T eq = 963(18)K, between 0.9 and 2.6µm, with a constant spacing of 0.01 cm −1 . We included collision-induced absorption from H 2 -H 2 and H 2 -He interactions (Richard et al. 2012). For each night we performed the CC for an optimal selection of GIANO-B orders (Giacobbe et al. 2021) and every phase over a lag vector corresponding to planet radial velocities (RVs) in the range −252 ≤ RV ≤ 252 km s −1 , in steps of 3 km s −1 . This range was chosen in order to cover all possible RVs of the planet. We then co-added the CC functions over the selected orders, nights, and the orbital phase after shifting them in the planet rest frame by assuming a circular orbit (e.g. Brogi et al. 2012Brogi et al. , 2014Brogi et al. , 2018. In doing so, we considered a range of planet RV semi-amplitudes 0 ≤ K p ≤ 200 km s −1 , in steps of 1.5 km s −1 . We then quantified the confidence level of our detec-Article number, page 3 of 10 A&A proofs: manuscript no. 43854corr tion by computing the S/N detection map (e.g. Brogi et al. 2012Brogi et al. , 2014Brogi et al. , 2018. For each investigated molecule, we divided the total CC matrix by its standard deviation (calculated by excluding the CC peak). We also used the Welch t-test to compare the 'intrail' values (CC values that carry signal) of the 2D CC matrix aligned in the planet rest frame to those 'out-of-trail' (CC values away from the WASP-69b's RV). We then converted the corresponding t value into the statistical significance at which these two distributions are not drawn from the same parent distribution. We applied this method for the same range of rest-frame velocity V rest and the planetary semi-amplitude K P adopted before. The final S/N and significance maps for the three nights combined together are shown in the upper and in the middle panels of Fig. 3, respectively. These results show the advantage of a multi-night approach; indeed even if the data quality does not allow one to always detect the signal in each individual night (see Fig. A.1), the detection significance increases when transits are co-added.
The detections in CC have been obtained with isothermal models of one molecule at a time with a volume mixing ratio (VMR) large enough to reveal a given molecule, if present in the atmosphere; it has indeed been demonstrated (e.g. Gandhi et al. 2020) that the CC is only weakly dependent on the absolute VMR. As a good compromise, we found a VMR∼10 −4 .
(iv) Confidence intervals via likelihood. The CC framework as applied in the previous analysis implicitly assumes that the signal amplitude and the S/N are uniform across the considered spectral range. However, this assumption is only valid at a first-order level. To investigate this aspect more thoroughly, we converted the CC values into LikeliHood (LH) mapping values using the approach proposed by Brogi & Line (2019) and used by Giacobbe et al. (2021). In this formalism, the LH is indeed set up by taking into account the model's line depth and the S/N order-by-order and spectrum-by-spectrum. Retrieving a peak in the LH at the correct position of the (V rest , K P ) space is hence an additional confirmation of our detection. We thus computed a log-LH function for each selected order, spectrum, and RV shift of the model across the K P versus V rest space. We used a lag vector corresponding to planet RVs in the range −20 ≤ RV ≤ 20 km s −1 , in steps of 1 km s −1 . To take into account any possible modification of the planetary signal due to the PCA process, on each theoretical model, we performed the exact telluric removal process applied to the observations. The final log-LH for each (V rest , K P ) couple was obtained by adding all the log-LHs over each order, night, and observed spectrum. Subsequently, as in Giacobbe et al. (2021), we introduced the line-intensity scaling factor S (−1.1<log S <1.1, in steps of 0.1) to understand how much the used models have to be scaled to match the observed data. For a model perfectly matching the data, S is exactly 1 (log S =0), that is to say the model and the data have the same average line amplitude compared to the local continuum. Scaling factors smaller than 1 indicate that spectral lines are too strong compared to the observations and, for example, the presence of a grey-cloud model can be invoked to account for that, as in Giacobbe et al. (2021). Furthermore, as in the single species models, the continuum can be approximate, and the introduction of the log S allows us to account for the muting effects of other species. On the contrary, in the CC framework, since the line depth does not enter into the calculation, the log S is not necessary. We then converted the LH values into confidence intervals by using the LH-ratio test between the maximum log-LH and each LH value in the (V rest , K P , and log S ) space. The LH confidence intervals' maps corresponding to the log S , maximising the detection, are shown in the bottom panel of Fig. 3.
We then evaluated the significance of our LH detection by performing a LH-ratio test between the maximum log-LH and the median of the 2D log-LH matrix values for the best-fit scaling factor which exhibited a statistical confidence greater than 1. For each considered molecule, Table 2 reports the results of our analysis in both the CC and the LH frameworks. We consider a molecule as detected if altogether the S/N, the detection significance, and the LH significance are higher than 3σ.

Results and discussion
We characterise the atmosphere of WASP-69b showing a multispecies molecular detection in the planet transmission spectrum.
Our results indicate the presence of CH 4 , NH 3 , CO, C 2 H 2 , and H 2 O, at more than a 3.3σ level ( Table 2). We did not detect the presence of HCN or CO 2 with a confidence level higher than 3σ (see Table 2). The retrieved K p values for all molecules are compatible within 1σ with the nominal K p value in the S/N and t-test maps (see Table 2 and Fig. 3, top and middle panels). We find that the same thing occurs for NH 3 , CO, and C 2 H 2 when switching to the LH framework, while the peaks of the CH 4 and H 2 O signals are slightly shifted, though consistent with the expected K p at the 2σ level (see Fig. 3, bottom panel). Since the LH function is more sensitive to the variance of the spectra, these small shifts may be due to the contribution of the other molecules, which are not modelled simultaneously when considering one molecule at a time.
We also checked whether our CO detection might be spurious, as stellar CO lines can still be an important contaminant of the planetary transmission spectrum (Brogi et al. 2016). To resolve the possible ambiguity, we evaluated the expected K P at which a spurious signal of stellar origin should have been maximised as K Pstar ∼ vsini sin(2πφ 1 ) ∼29.2±5.3 km/s (by using the values in Table 1). While the above equation provides an estimate of the expected position of stellar residuals in the (V rest ,K P ) map, the actual position will depend on the interplay between the stellar and the planetary signal, and it might also be affected by any partial removal of stellar residuals via PCA. Thus, one should not expect stellar residuals to arise exactly at the K Pstar computed above. In fact, there are cases in the literature where stellar CO causes the planet signal to be detected at higher K P than the orbit of the planet allows (e.g. Chiavassa & Brogi 2019, Fig.7), which seems counter-intuitive at first.
With the LH framework applied, we could also infer to which extent the model is a plausible representation of the data. However, we measured low log S values for several molecules, as single-species models have a higher variance (deeper lines) than an actual spectrum, because they do not contain the additional opacity from other species. We remark that the assumption of a fixed VMR composition has an impact on the line strengths S , and thus it can lead to ambiguities in the scaling factors' interpretation. Consequently, the comparison between the S values for the different molecules is not straightforward, meaning that the relative abundances between the different species cannot be estimated easily. The scaling factors obtained here ( Table 2) thus must be considered representative and taken with caution. The large number of individual species' detections (and nondetections) on WASP-69b provides important constraints on the physical processes that govern the planet's atmosphere. Using the Pyrat Bay (Cubillos & Blecic 2021) modelling framework (see Appendix Sect. B), we explored physically-plausible atmospheric scenarios in radiative and thermochemical equilibrium that could be consistent with the observed molecular detections. Given the striking difference in composition that is expected for scenarios with elemental ratios of C/O<1 or C/O>1 (see, e.g. Tsuji 1973;Madhusudhan 2012;Moses et al. 2013), we adopted the C/O ratio and the atmospheric metallicity as the main proxies to explore different plausible compositions for WASP-69b. We found that very high-metallicity regimes ( 10× solar) appear disfavoured by the observations because, for C/O<1 scenarios, the non-detected CO 2 molecule would present detectable features since its abundance increases faster with metallicity than for the other species. For C/O>1 scenarios at high metallicities, the observed H 2 O features would not be detected since they would be obscured by the CH 4 absorption, due to the much lower abundance of H 2 O relative to CH 4 . When comparing solar metallicity models with C/O ratios lower and greater than one (Fig. 4), the observations seem to be more consistent with the latter. First, while transmission models for both regimes show clear spectral features of H 2 O, CH 4 , and CO, the stronger detection of CH 4 than H 2 O may favour the C/O>1 scenario (in which CH 4 is predominant). Then, the NH 3 molecule presents a wide-range yet weak features across the spectrum in both C/O regimes, thus, it does not particularly favour any of the scenarios. However, the C 2 H 2 is hard to explain in equilibrium conditions since its estimated abundance does not yield detectable features in any of the above scenarios. A way to resolve this conundrum is by including disequilibrium chemistry (see e.g. Moses 2014, and references therein), which would significantly enhance the abundances of NH 3 and C 2 H 2 by dragging material from the deep interior to the upper atmospheric layers (vertical quenching), by photochemical production, or both. Certainly, the much greater abundance of C 2 H 2 for the C/O>1 case makes the detection of this molecule more plausible for this scenario than for C/O<1 (a difference of ∼5 orders of magnitude for the cases shown in Fig. 4), as long as disequilibrium processes help to bridge the gap between the WASP-69b detections and the equilibrium theoretical predictions.
To give a more quantitative interpretation of the physical and chemical conditions of WASP-69b's atmosphere, we computed and tested three different scenarios that could be a possible representation of WASP-69b's atmosphere: -(i) solar metallicity and solar C/O (C/O=0.55); -(ii) solar metallicity and carbon rich (C/O=2.0); and -(iii) solar metallicity, carbon rich (C/O=2.0), and quenched C 2 H 2 and NH 3 to VMRs of 10 −6 − 10 −7 .
Article number, page 5 of 10 A&A proofs: manuscript no. 43854corr  Notes: From left to right: the investigated molecule, the adopted line list, the CC result with theoretical models both in the S/N and the significance framework, the LH findings, and the status of the detection (with and , indicating a detection and a non-detection, respectively). For each framework, the planet orbital semi-amplitude (Kp 0 )as well as the velocity in the planet rest frame (V rest0 ) of the CC and LH peak are reported. The S/N, the significance (σ), and the scaling factor (log S ) maximising the detection are also present. We used a grid of log S values, thus they are shown without an error bar. As opposed to the single-molecule analyses from Section 3, here we generated more realistic 'full spectrum' models combining opacities from all species. To quantify which scenario describes the observations better, we applied the Wilks theorem (Wilks 1938) by computing the LH-ratio test relative to the LH-maxima (following the procedure of Giacobbe et al. 2021). Our analysis (see Table A.3) indicates that scenario (iii) is marginally (∼ 1σ) preferred compared to scenario (ii), and it is significantly favoured (∼ 3.5σ) over scenario (i). The atmosphere of WASP-69 thus seems to be better described by a scenario where quenching of NH 3 and CH 4 is taking place. However, we need to stress two aspects: (1) this result does not imply that (iii) is the bestfitting scenario, but it is the preferred one within the tested models; and (2) we assumed a fixed solar metallicity, which might not be the most representative of the planet's atmosphere. Dealing with (1) and (2) implies a retrieval analysis that goes beyond the scope of this paper.
The emerging picture provides some preliminary indications on WASP-69b's formation history once put in the context of the host star. WASP-69's metallicity is super-solar (Anderson et al. 2014), with recent works suggesting metallicity values 2-3 times higher than that of the Sun (Magrini et al. 2022). The atmospheric scenario favoured by our modelling efforts points to a solar metallicity for WASP-69b's atmosphere, meaning that the atmospheric metallicity is sub-stellar. The combination of substellar planetary metallicity and planetary C/O>1 is the signpost of giant planets whose content of heavy elements is dominated by the accretion of gas with a limited contribution of solids (Turrini et al. 2021a,b). Furthermore, this combination suggests that the chemical structure of WASP-69b's native circumstellar disc was inherited from the parent molecular cloud (Pacetti et al. 2022). We also note that cloud formation processes offer an alternative interpretation of these results. In an atmosphere with prominent cloud formation, condensation locks and transports heavy elements such as carbon or oxygen across an atmosphere, which can significantly modify the local gas-phase carbon-to-oxygen ratios. For example, for a solar-like composition, oxygen depletion due to condensation can be strong enough to raise C/O ratios up to ≈ 0.7 (Bilger et al. 2013). For carbon-rich atmospheres, on the other hand, carbon depletion can make an atmosphere locally more oxygen rich, approaching values of C/O ≈ 1.0 at the regions with the strongest carbon depletion (Helling et al. 2017).
While the picture described above will need to be validated by future studies, the simultaneous presence of several chemical species constitutes a very important step forward, as so far only a few molecules have been firmly detected simultaneously in the atmosphere of warm giant planets (e.g. Tsiaras et al. 2018;Fisher & Heng 2018). Our findings reinforce the notion that a rich atmospheric chemistry is not the sole dominion of hotter planets such as HD 209458b (Giacobbe et al. 2021). The frontier in the characterisation of exoplanetary atmospheres is expanding further, with even better prospects in the future when stronger constraints on molecular abundances, and hence on metallicity and elemental ratios (e.g. Brogi et al. 2017;Brogi & Line 2019;Guilluy et al. 2022), will be placed based on systematic joint analyses of HR ground-based observations gathered with GIANO-B, CARMENES, Spirou, NIRPS, and CRIRES+, for instance, and LR space-borne spectra such as those obtained with HST, JWST, and Ariel. Discarded orders in the basic data reduction step of the analysis, i.e. alignment in the telluric rest frame and refinement of the wavelength solution.