© The Authors 2024

1 Introduction

The majority of known exoplanets have orbital and physical characteristics that are different from those of the Solar System planets. This is the case, for example, for hot and warm Neptunes, which are planets with a mass similar to that of Neptune and with equilibrium temperatures of T_eq ≳ 1000 K and T_eq ≲ 1000 K, respectively. The existence of such planets so close to their host stars (orbital periods P_orb ≲ 10 days) provides a unique opportunity to study atmospheric physical and chemical conditions that cannot be studied in the Solar System.

The study of exoplanetary atmospheres makes a crucial contribution to the exoplanet characterisation process. For example, by knowing the chemical composition of the primary atmosphere of an exoplanet, it is possible to constrain its formation and evolution path based on the study of relative elemental abundances (e.g. Öberg et al. 2011; Madhusudhan et al. 2017; Madhusudhan 2019; Banzatti et al. 2020; Bitsch et al. 2022; Pacetti et al. 2022).

By studying exoplanetary atmospheres in the near-infrared (NIR), it is possible to probe deep layers (down to pressures P ≈ 0.1 bar) where molecular species dominate the atmospheric composition and absorb the IR light through thousands of rotational-vibrational transitions. Over the past few years, the large number of transiting exoplanet discoveries and the possibility to use space instruments have made low-resolution spectroscopy (LRS) the most used technique for exoplanetary atmospheric characterisation (Madhusudhan 2019); in particular, for probing the terminator region of planetary atmospheres via the transmission-spectroscopy technique. This method is based on measuring how the planetary effective radius varies with wavelengths during transit because of the absorption of the stellar light by the chemical species present in the atmosphere of the target.

An alternative technique for probing exoplanetary atmospheres using high spectral-resolution instruments is high-resolution spectroscopy (HRS) from ground-based observatories (see Birkby 2018 for a review). For what concerns the NIR studies, different molecular species have been identified in the atmosphere of hot Jupiters with the HRS technique, such as CO (Snellen et al. 2010), H₂O (Birkby et al. 2013), TiO (Nugroho et al. 2017), HCN (Hawker et al. 2018), CH₄ (Guilluy et al. 2019), NH₃ and C₂H₂ (Giacobbe et al. 2021), and OH (Nugroho et al. 2021). Having improved the data analysis approach (e.g. employing the principal component analysis to remove the telluric and stellar contaminations and performing an optimal selection of the spectral orders), it is now possible to simultaneously detect multiple molecular species in the atmospheres of both hot and warm giant planets (e.g. Giacobbe et al. 2021; Guilluy et al. 2022; Carleo et al. 2022).

With respect to LRS, HRS is more sensitive to line position than depth, allowing a self-calibration of data, and has the advantage of combining the signal of thousands of spectrally resolved absorption or emission lines rather than bands from a chemical species. On the other hand, LRS can take advantage of the two space telescopes currently in operation (HST and JWST), which allow observers to avoid telluric contamination. A future combination of these two complementary techniques will improve the amount of information that can be inferred about the physics and chemistry of exoplanetary atmospheres (Brogi & Line 2019).

The ideal targets for atmospheric studies are close-in giant exoplanets given the higher planet–star radius and flux contrast. In the case of Neptune-size planets, the relatively small radius and the smaller atmospheric scale height combined with the possible presence of clouds or hazes reduce the expected amplitude of the atmospheric absorption, making the atmospheric characterisation of these targets more difficult. Indeed, there are only a few works in the literature reporting the detection of molecular species in the atmosphere of warm and hot-Neptunes (i.e. Fraine et al. 2014; Benneke et al. 2019; Bézard et al. 2022; Kreidberg et al. 2022; Brande et al. 2022; A-thano et al. 2023; Mikal-Evans et al. 2023), and most of them report the detection of water vapour obtained with LRS using data from HST/WFC3 (A-thano et al. 2023 also report the presence of titanium oxide on HAT-P-26 b, while Bézard et al. 2022 only report the presence of methane on K2-18 b).

A remarkable target for atmospheric studies is HAT-P-11 b, a warm Neptune-size exoplanet (the first of this class of planets discovered with transit searches) orbiting a K4V-class star (Bakos et al. 2010). Fraine et al. (2014), Tsiaras et al. (2018), Chachan et al. (2019), and Cubillos et al. (2022) detected the presence of water vapour in its atmosphere at low resolution by analysing transmission spectra from HST. With the same data, Welbanks et al. (2019) estimated an abundance of . The analyses presented by Chachan et al. (2019; HST/WFC3+STIS) and by Cubillos et al. (2022; HST/WFC3) also suggest the presence of methane; however, when also considering the Spitzer observations, both these analyses found no evidence for the presence of this molecule due to an offset between the Spitzer and HST transit depths. The detection of atmospheric absorption in the He metastable 1083 nm triplet during transit (Allart et al. 2018; Mansfield et al. 2018) also provided complementary constraints on the size of the planet’s upper atmosphere (extending beyond two planetary radii) and the corresponding mass-loss rate (the planet has only lost up to a few percent of its mass over its history, leaving its bulk composition largely unaffected). Due to the large brightness of its host star (H = 7.131 ± 0.021 mag, Cutri et al. 2003) and the presence of already detected chemical species (i.e. He and H₂O) in its atmosphere, this target provides a great opportunity to characterise the atmosphere of warm Neptunes. In addition, the orbit of HAT-P-11 b is eccentric and Sanchis-Ojeda & Winn (2011) estimated a high obliquity angle between the orbital plane of the planet and the equatorial plane of its host star, indicating a quasi-polar orbit. This planet therefore offers a rare occasion to explore complex planetary evolution paths, such as the one that led to its current orbital configuration.

To study the atmosphere of the planet HAT-P-11 b at high spectral resolution, precise and accurate knowledge of its orbit (in particular the eccentricity e and the planetary argument of periastron ω_p) is mandatory. Multiple analyses of the orbital parameters of HAT-P-11 b can be found in the literature (i.e. Bakos et al. 2010; Southworth 2011; Knutson et al. 2014a; Huber et al. 2017; Allart et al. 2018; Yee et al. 2018). The most recent is that of Yee et al. (2018), which is based on the analysis of the radial-velocity (RV) data of HAT-P-11. In particular, these authors obtained the following values for the eccentricity and the planetary argument of periastron: . The most precise estimation of the orbital parameters was reported by Huber et al. (2017). These latter authors simultaneously modelled the planetary transits and secondary eclipses in the Kepler data, obtaining the following values: . Although these two sets of estimations are compatible with each other at the 2 σ level, mainly due to the large uncertainty on the values by Yee et al. (2018), it is important to note that the two values of e differ by about 0.05, with the analysis of Huber et al. (2017) pointing towards a higher value of the eccentricity. Even such a small inaccuracy on the value of e still produces an RV shift of roughly ≈ 10 km s⁻¹ (larger than the GIANO-B spectral resolution of 6 km s⁻¹) during the transit, with severe impacts on the atmospheric characterisation using the HRS method. To improve the accuracy of the orbital solution, we decided to determine the HAT-P-11 b orbital parameters through independent analysis of both the Kepler photometric data and the RV measurements. Moreover, by doing so, we reviewed the physical parameters of HAT-P-11 b and the architecture of the HAT-P-11 planetary system. All the details of this preliminary analysis are described in the following subsections. The main parameters of the HAT-P-11 planetary system are summarised in Table 1.

In this work, we report a review of the physical and architectural properties of the HAT-P-11 planetary system and the results of the study of four transit events of HAT-P-11 b recorded with GIANO-B, the high-resolution NIR échelle spectrograph mounted at the 3.58 m Telescopio Nazionale Galileo (TNG), in order to probe the presence of eight molecular species in its atmosphere with the transmission HRS as part of the Global Architecture of Planetary Systems (GAPS) Project¹ and in particular as part of the exoplanetary atmospheres characterisation subprogramme, described in Guilluy et al. (2022).

The paper is organised as follows: in Sec. 2 we review the HAT-P-11 planetary system; in Sec. 3 we describe the high-resolution transit observations and data analysis process in detail, and discuss the results of the atmospheric characterisation. Finally, our conclusions and future perspectives are reported in Sec. 4.

2 Revisitation of the HAT-P-11 planetary system

2.1 Kepler light-curve data analysis

We downloaded the Kepler light curves of HAT-P-11 b from the Mikulski Archive for Space Telescopes (MAST²). These consist of short-cadence (60 s) light curves from 14 quarters out of the total 17 observed from 2 May 2009 to 11 May 2013. The short-cadence light curves of each quarter are subdivided into tranches of three except for quarters 0 and 1 with a single and quarter 17 with two light curves. This results in a total of 37 separate light curves to be analysed independently.

2.1.1 Transits and occultations

We trimmed the Kepler light curves around the transits and occultations, whose central times are predicted using the ephemeris of Huber et al. (2017). The width of each time interval is equal to three times the transit duration T₁₄ and is large enough to include the time of the eclipses for both the orbital solutions discussed above. We discard all the transits and occultations that, due to gaps in the Kepler light curve, do not cover the whole time interval of 3 T₁₄.

Many transits of HAT-P-11 b show clear signs of starspot crossings (see, e.g. Sanchis-Ojeda & Winn 2011; Béky et al. 2014; Morris et al. 2017b; Scandariato et al. 2017 for a detailed analysis of the transit anomalies of HAT-P-11 b), which may bias the fit of the light curve. For this reason, we adopted an iterative approach aimed at selecting the transits with minimum evidence of anomalies. Using the orbital period derived by Huber et al. (2017), we phase-folded the planetary transits and we processed the phase-folded data using a running median average. The width of the running window was 15 s, which is less than the cadence of the light curves. This guarantees that the averaged transit profile is negligibly time-smoothed by our approach.

Once the average transit profile was obtained, we computed the Median Absolute Deviation (MAD) of the data with respect to it. We then rejected all the transits with at least one data point located more than 6 MAD above the averaged transit profile. We iteratively repeated this process until no additional transit was rejected. At the end of this process, we noticed that a few transits passed our selection despite showing correlated noise due to either bad data detrending or stellar activity. We therefore refined the selection of the transits in the following way. For each transit, we first computed the residuals with respect to the averaged profile. Then, we associated each transit with the standard deviation of the corresponding residuals. Finally, we rejected the 10% transits with the largest standard deviations. This produced a final list of 64 bona fide transits free from anomalies above the noise level. Of course, the possible presence of non-crossed starspots can influence the stellar flux level and therefore the measured transit depth; however, this effect mainly affects the value of the planetary radius (with an over-estimation of ≈1 %³) rather than the orbital parameters and is therefore negligible for the main scope of this work.

For the occultations, we adopted a similar approach. The only difference is that we do not expect anomalies during the eclipses, and therefore the first iterative selection of “unspotted” light curves is skipped. The final selection is made up of 196 occultations.

In summary, the dataset that we used to fit the orbit of the planet is composed of two subsets of data, centred on the transits and occultations respectively. To save computation time, we fit the phase-folded and rebinned data (1 min). The timestamps of the rebinned light curve were defined so that the time of transit provided by Huber et al. (2017) corresponds to the origin of the time axis.

The transit profile was modelled with the quadratic limb darkening (LD) law provided by Mandel & Agol (2002), with the reparametrisation of the coefficients proposed by Kipping (2013). Similarly, the occultation profile was modelled following Mandel & Agol (2002), but assuming that the planetary dayside is uniformly bright (Singh et al. 2022; Scandariato et al. 2022). Since we have operated a different data rejection with respect to Huber et al. (2017), we re-derived all the orbital parameters except for the orbital period, which we fixed to the best estimate of Huber et al. (2017) to phase-fold the data. The free parameters of the model are the stellar density ρ_⋆, the time of transit T₀, the planet-to-star radius ratio R_p/R_⋆, the impact parameter b, the LD coefficients q₁ and q₂, the e_c = e cos ω_⋆ and e_s = e sin ω_⋆ parameters (where e and ω_⋆ are the orbital eccentricity and the stellar argument of periastron, respectively) and the occultation depth δ_ecl.

In the model, we also included a jitter term and a renormalisation coefficient independently for the transit and eclipse subcurves. The two jitter terms take into account the fact that the two subcurves have different noise properties, being the combination of a different number of transits/eclipses. The renormalisation coefficients fix the preliminary normalisation of the transit and eclipse subcurves computed by the extraction pipeline.

We adopted a maximum-likelihood Bayesian approach where the data were fitted by running a Monte Carlo sampling of the parameter space. The parameter space was defined by the priors listed in Table 2. We remark that we have used uniform priors for all the parameters. The priors on e_c and e_s have been conveniently set to span a large range around the expected values and include both the orbital solutions of Huber et al. (2017) and Yee et al. (2018), while saving computation time.

For the log-likelihood maximisation, we first searched the parameter space for the global maximum position using the Python package PyDE⁴. Afterwards, we sampled the posterior probability distribution of the model parameters using the emcee package version 3.1.3 (Foreman-Mackey et al. 2013). Given the demand for resources for the model fitting, we ran the code in the HOTCAT computing infrastructure (Bertocco et al. 2020; Taffoni et al. 2020). We let the chains run for 250 000 steps, long enough to ensure formal convergence. The best-fitting model of the light curve, together with the corresponding residuals, is shown separately in Figs. 1 and 2, for the transit and the occultation, respectively. The list of the free parameters and their corresponding priors and best-fitting 1σ confidence interval is given in Table 2. Our estimates are consistent with previous analyses within 2σ.

2.1.2 Albedo and equilibrium temperature

The day-side flux of an exoplanet is a combination of reflected starlight off the planet’s illuminated hemisphere and its thermal irradiation. The former is parameterised by the geometric albedo A_g while the latter is parameterised by the brightness temperature T_d(Δλ), which is a measure of the day-side temperature in a given wavelength interval Δλ (Santerne et al. 2011; Singh et al. 2022). Consequently, the occultation depth can be expressed as the following: (1)

where h is the Planck constant, k_B the Boltzmann constant, c the speed of light, d_sec the distance of the planet from the star during the secondary eclipse and is the stellar Kurucz flux (Castelli & Kurucz 2003) (computed for T_eff = 4750 K, log g = 4.5, and [Fe/H] = 0.2). Both the planetary and the stellar fluxes are integrated over the Kepler passband Δλ with the corresponding response function Ω_λ.

For the derived occultation depth, the relationship between the geometric albedo and the brightness temperature is shown in Fig. 3. For HAT-P-11 b, the thermal contribution to the observed depth in optical passbands is practically negligible given the low temperature of its atmosphere. Therefore, the occultation depth is the result of a highly reflective atmosphere. We derive a geometric albedo of corresponding to the occultation depth of ppm reported in Table 2. We therefore confirm and improve the results obtained by Huber et al. (2017). Following Han et al. (2014), we assume that the geometric albedo and Bond albedo are related via . We use this A_b to estimate the planet’s day-side equilibrium temperature⁵ as a function of time following (Cowan & Agol 2011): (2)

where ε is the heat re-circulation efficiency and d(t) the distance of the planet from the star as a function of the time (to take into account the eccentricity of the orbit).

We considered two extreme scenarios: one with inefficient (ε = 0) and another with complete (ε = 1) heat re-circulation. We report in the plot the equilibrium temperature estimates at the occultation position. As a result of the varying stellar irradiation, the temperature estimates at the periastron are: K and K for ε = 0 and ε = 1, respectively, and at the apoastron are K and K for ε = 0 and ε = 1, respectively. During the transit, the planet-to-star separation distance is d_tr = 13.33 ± 0.26 R_⋆ and therefore, the corresponding day-side temperatures are K and K assuming no changes in planetary albedo throughout the eccentric orbit. The pairs of T_d values we report represent the extremes of the range in which the planet’s equilibrium temperature is expected to be at different planet positions along its orbit. At the transit position, we consider the scenario ε = 1, that is, a uniform temperature distribution throughout the planet, so that we can use T_eq = K as the equilibrium temperature around the day-night terminator to build our atmospheric transmission spectrum models described in Sect. 3.3.

Fig. 1 Transit light curve analysis. Top panel: phase-folded Kepler transit light curve of HAT-P-11 b binned by 1 min, based on 64 bona fide transits (see text). Our best-fitting model of the transit in the Kepler bandpass is overplotted in grey. Bottom panel: residuals of the fit.

Fig. 2 Occultation light curve analysis. Top panel: phase-folded Kepler occultation light curve of HAT-P-11 b binned by 1 min, based on 196 occultations (see text). Our best-fitting model of the occultation in the Kepler bandpass is overplotted with a black solid line. The data have been binned for clarity (black dots with error bars). Bottom panel: residuals of the fit.

Fig. 3 Geometric albedo (Ag) estimated as a function of the day-side temperature for the measured occultation depth . The cyan lines represent the 1σ uncertainty curves. The dash-dotted and dashed lines represent the variation of Ag with varying Td for the 2 heat re-circulation cases (ε) considered, computed at d = dsec = 14.83 ± 0.30 R⋆ (that is, during the occultation). The 2 orange shaded regions correspond to the intersection of the curves, identifying the average HAT-P-11 b day-side temperatures for the 2 extreme scenarios of ε, computed at d = dsec. These 2 values of the equilibrium temperature with the associated uncertainties are reported in the legend.

2.2 Radial-velocity data analysis

We also analysed 180 publicly available radial velocities of HAT-P-11, which were obtained with the HIRES at Keck spectrograph by Bakos et al. (2010) and Yee et al. (2018), after discarding the in-transit measurements to avoid the Rossiter-McLaughlin effect (Rossiter 1924; McLaughlin 1924) and three outliers at the observing epochs 4334.9662, 4957.0433 and 7933.0122 BJD_tdb – 2 450 000, which were identified in the residuals of our RV model (Sect. 2.2.2) through the Chauvenet’s criterion (e.g. Bonomo et al. 2023).

2.2.1 Planet c or magnetic activity cycle?

The HIRES RVs show a clear long-term trend (Fig. 4, top panel) which was attributed by Yee et al. (2018) to a second planet companion, HAT-P-11 c, with P_orb ~ 3410 d (9.3 yr), M_p sin i ≈ 1.6 M_J, and e ≈ 0.6, and, to a lesser extent, to the stellar magnetic activity cycle. Despite the very similar behaviour of the trend in both the S-index and the RVs (Fig. 4), the stellar activity cycle was not deemed sufficient by Yee et al. (2018) to account for the RV variations in the long term for three main reasons (see their Sect. 3 for more details): (i) the large RV semi-amplitude (~30 m s⁻¹) of the long-term signal compared to semi-amplitudes of ≲ 10 m s⁻¹ observed by Lovis et al. (2011) for magnetic activity cycles in ~300 solar-type stars; (ii) the presence of a shift of ~500 days between the minimum of the S-index and that of the RVs (Fig. 4); and (iii) the relatively weak correlation between the S-index and RV measurements with a Pearson’s coefficient of ~0.34.

In our view, these three motivations do not provide strong evidence that the long-term signal is planetary in origin. Indeed, concerning (i), HAT-P-11 is considerably more active than the stars in the sample studied by Lovis et al. (2011), with a log of −4.35 (Morris et al. 2017a) higher than the typical log of −4.8 −5.0 of that stellar sample. Since one of those stars, namely HD 21693, shows a semi-amplitude of ~10 m s⁻¹ for log = −4.89 (see Fig. 16 in Lovis et al. 2011), a semi-amplitude of ~30 m s⁻¹ for the RV variation associated with the activity cycle is certainly possible for an unusually active star such as HAT-P-11 (Morris et al. 2017a).

Regarding (ii), detailed studies of the correlation between RV and S-index measurements by Meunier et al. (2019) (therefore subsequent to Yee et al. 2018) showed that a combination of geometrical effects (stellar inclinations and butterfly diagrams) and variations of magnetic activity level over time may easily produce hysteresis patterns, and hence temporal shifts of a few hundreds of days in the minima of the RV and S-index variations (see Fig. 8 in Meunier et al. 2019). For example, the minimum of the long-term RV variations of the above-mentioned star HD 21693, which are caused by the magnetic activity cycle, also leads the minimum of the log variations by ~500 d, similarly to HAT-P-11 (Lovis et al. 2011; Meunier et al. 2019). Given the high inclination i_⋆ = 100 ± 2 deg of the host star HAT-P-11 and a starspot latitudinal distribution similar to the solar butterfly diagrams, as unveiled from the occultations of starspots by HAT-P-11 b during transits (Morris et al. 2017b), the temporal difference in the minima of the RV and S-index variations could be due to the hysteresis patterns described by Meunier et al. (2019). After all, the fact that the magnetic activity cycle from the S-index was found to have the same periodicity as the hypothetical planet c (Morris et al. 2017a) remains suspicious.

Last but not least, even the absence of a strong correlation between the S-index and RV measurements does not necessarily lean towards the planetary origin of the RV long-term signal; this is because a relatively low correlation, at least partly, ensues from the temporal shift between the S-index and RV variations, while it is much higher in the first ~1000 d of observations (Knutson et al. 2014b). Since we observed HAT-P-11 in GIARPS mode, we extracted the HAT-P-11 RVs and activity indicators from the HARPS-N spectra of the four transits of HAT-P-11 b to look at their behaviour. For that purpose, we used the online v3.7 data reduction software and cross-correlated the HARPS-N spectra with a K5 V synthetic stellar template (Pepe et al. 2002). The variations of the HARPS-N RVs, S-index and full width at half maximum (FWHM) of the cross-correlation function show an almost identical behaviour (Fig. 5): the Pearson’s correlation coefficient is 0.94 between RVs and S-index, and 0.98 between RVs and FWHM. This suggests that the ~9-10 yr long-term RV signal is more likely due to the stellar activity cycle than to the long-period eccentric companion HAT-P-11 c. Moreover, if the RV measurements of the last transit night actually caught the minimum of the activity cycle, given that the HARPS-N RV peak-to-peak variation of ~60 m s⁻¹ is the same as that observed by Yee et al. (2018), temporal shifts between the RV and S-index minima as caused by hysteresis phenomena may not have occurred in the current activity cycle.

Even though our data suggest that the magnetic activity cycle is the most plausible origin of the long-term RV signal, HIPPARCOS-Gaia absolute astrometry still provides hints that a long-period companion may actually exist. The catalogues of astrometric accelerations produced by Brandt (2018, 2021) and Kervella et al. (2019, 2022) indicate the presence of a proper motion anomaly (PMA) at the mean Gaia epoch, whose signal-to-noise ratio (S/N) grows from S/N ≃ 2.1 to S /N ≃ 4.7 and from S/N ≃ 2.8 to S/N ≃ 4.3 between the Gaia DR2 and Gaia EDR3 editions of the former and latter catalogue, respectively. Indeed, Xuan & Wyatt (2020) used the HIPPARCOS-Gaia DR2 PMA values in combination with the Bakos et al. (2010) RVs of HAT-P-11 to constrain the true mass and inclination of the putative companion HAT-P-11 c. Figure 6 shows the HIPPARCOS-Gaia DR2 and DR3 PMA sensitivity curves based on Eq. (15) of Kervella et al. (2019), along with the minimum-mass value of HAT-P-11 c derived by Yee et al. (2018) and the best-fit true mass obtained by Xuan & Wyatt (2020). The HIPPARCOS-Gaia DR3 PMA sensitivity curve indicates that, at the orbital separation of HAT-P-11 c, a companion inducing a statistically significant PMA should have a mass of ~3.5 M_Jup. In the Xuan & Wyatt (2020) analysis, the true mass value of HAT-P-11 c falls below the PMA sensitivity curve, with a companion having true mass equal to the minimum mass from Yee et al. (2018) compatible at ~1.4 σ. This is somewhat surprising as such a companion is not expected to produce a PMA with S/N ≳ 3. As the PMA technique heavily relies on the constraints from the RVs in order to successfully provide inferences on the mass and inclination of a companion, it is therefore possible that, if the long-term modulation in the RVs is actually dominated by the activity cycle, then HAT-P-11 c exists at larger separation and with a different mass than those inferred by Xuan & Wyatt (2020).

Fig. 4 HIRES RV (top panel) and Call S-index (bottom panel) measurements of HAT-P-11. The two-time series show almost identical long-term variations with a shift of ~500 days in the minimum.

Fig. 5 HARPS-N RV (top panel), Call S-index (middle panel), and FWHM (bottom panel) measurements of HAT-P-11 during the four HAT-P-11 b transit nights for atmospheric characterisation. The three time series are highly correlated showing the same long-term trend, with no apparent shifts in the minima of variations. The vertical dashed line indicates the predicted periastron time of the hypothetical planet c, and the grey area shows its 1σ error bar accounting for the uncertainty on the orbital period from Yee et al. (2018). Note: The HARPS-N radial velocities were divided by their median of −63420.4 m s−1, to make a straightforward comparison with the HIRES radial velocities in Fig. 4.

Fig. 6 Sensitivity of the PMA technique to companions of given mass and orbital separation orbiting HAT-P-11. The black long-dashed and solid curves correspond to the combinations of mass and orbital radius explaining the PMA values at the mean Gaia DR2 and DR3 epochs, respectively. The shaded light blue region corresponds to the 1σ uncertainty domain of the DR3 PMA, while the shaded magenta region encompasses the 1 σ uncertainty of the DR2 PMA. The red diamond indicates the separation and minimum mass of the HAT-P-11 c companion proposed by Yee et al. (2018), while the red hexagon corresponds to the true mass value determined by the Xuan & Wyatt (2020) analysis.

2.2.2 Radial-velocity modelling and improved mass determination for HAT-P-11 b

In the lack of strong evidence that the RV long-term trend is caused by the planet c with the orbital parameters given in Yee et al. (2018) for the reasons explained above, we modelled the HIRES RVs with a Keplerian orbit for the transiting planet HAT-P-11 b only, which has six free parameters: T₀, P, e_c, e_s, the RV semi-amplitude, K_⋆, and the RV zero point, γ.

To account for non-stationary stellar variations produced by magnetic activity phenomena, we used Gaussian-process (GP) regression (e.g. Haywood et al. 2014; Haywood 2015; Grunblatt et al. 2015) with three different kernels, namely the squared-exponential (SE) kernel: (3)

the quasi-periodic (QP) kernel (4)

and a third (QP-SE) kernel that is the sum of the QP and SE kernels in such a way as to simultaneously model the stellar activity variations on both short-term rotation timescales (Eq. (4)) and long-term activity cycle timescales (Eq. (3)), namely (5)

where t and t′ are the epochs at two different RV observations, h is the semi-amplitude of the correlated noise, λ is the correlation decay timescale, P_rot is the period of the quasi-periodic variations, w is the inverse complexity harmonic parameter, is the formal uncertainty of the RV point at time t, and σ_jit is the uncorrelated jitter term, which would absorb any extra white noise not modelled by either the Keplerian or the GP.

We employed Bayesian differential evolution Markov chain Monte Carlo (DE-MCMC; Ter Braak 2006; Eastman et al. 2013; Bonomo et al. 2015) techniques to derive the posterior distributions of the model parameters, by using the same prescriptions for the number and convergence of the DE-MCMC chains given in Eastman et al. (2013) and Ford (2006). We imposed Gaussian priors on T₀ and P from the Kepler transit ephemeris and uniform priors on K_⋆ and γ as well as on the GP hyper-parameters h, λ, P_rot, and w, with the boundaries specified in Table 3 for each GP kernel. As for the eccentricity and stellar argument of peri-astron, we ran two analyses per kernel: in the first one, we used uninformative priors on e and ω_⋆, while in the second one we adopted Gaussian priors from the results of the modelling of the optical secondary eclipse (Sec. 2.1). We took the medians and the 15.87%-84.14% quantiles of the posterior distributions as the values and 1 σ uncertainties of the fitted and derived parameters.

In the first analysis with uniform priors on e and ω, we determined e = 0.277 ± 0.025 and ω_⋆ = 26.8 ± 8.6 with the SE kernel, e = 0.288 ± 0.022 and ω_⋆ = 27.9 ± 7.0 with the QP kernel, and e = 0.290 ± 0.021 and ω_e = 29.2 ± 6.7 with the QP-SE kernel. These are consistent with the e and ω values derived in Sec. 2.1 at ≲ 1.5 σ and ≲ 2.3 σ, respectively. In the second analysis, we found e = 0.2608 ± 0.0086 and ω_⋆ = 13.5 ± 2.7 deg with the SE kernel, e = 0.2638 ± 0.0090 and ω_⋆ = 14.2 ± 2.6 deg with the QP kernel, and e = 0.2654 ± 0.0091 and ω_⋆ = 14.7 ± 2.6 deg with the QP-SE kernel. These e and ω determinations are closer to the values from the secondary eclipse (~0.6 σ and ~0.8 σ, respectively) as expected from the use of the Gaussian priors on them. The radial-velocity semi-amplitude, K_⋆, does not vary from the first to the second analysis for a given kernel, but was found to be slightly higher for the SE kernel, that is, K_⋆ = 11.20 ± 0.50 m s⁻¹, to be compared to K_⋆ = 10.75 ± 0.41 m s⁻¹ and K_⋆ = 10.78 ± 0.42 m s⁻¹ for the QP and QP-SE kernels, respectively (see Table 3).

By using the Bayesian information criterion (BIC) as a proxy for the Bayesian evidence, we found that the model with the QP-SE kernel is the most favoured, while that with the SE kernel is highly disfavoured. We therefore adopted the orbital parameters of the former (QP-SE) model (Table 3), which has also a more physical rationale because the rotation and activity cycle signals were modelled with two different (QP and SE) kernels. On the other hand, the QP kernel had to adapt to fit the activity cycle long-term variation in addition to the rotational signal, with its hyper-parameters h and λ taking intermediate values between h_rot and h_cycle, and and λ_cycle, in the third kernel (Eq. (5)). We note that the GP models with both the QP and QP-SE kernels properly retrieved a stellar rotation period of P_rot ~ 32-33 days, close to P_rot ~ 29-30 days as estimated from the Kepler photometry, despite the large uniform prior adopted (see Table 3). We show the best-fit GP+Keplerian models as a function of time in the left panel of Fig. 7, and the Keplerian orbit due to HAT-P-11 b as a function of the orbital phase, after the removal of the GP activity model, in the right panel of the same figure.

For the rest of our analysis, we decided to use the orbital solution that we obtained from the analysis of transits and occultations due to the higher precision/accuracy in the determination of e and ω. We combined the stellar parameters, the transit parameters from the Kepler light curve (Sect. 2.1.1), and the RV parameters to derive a mass of M_p = 0.0787 ± 0.0048 M_J (M_p = 25.0 ± 1.5 M_⊕) and a mean density of ρ_p = 1.172 ± 0.085 g cm⁻³, for HAT-P-11 b. Finally, by knowing both the mass of the star and the planet, we computed the value of the planetary RV semi-amplitude K_p, which we used for the atmospheric characterisation in Sec. 3. All the derived parameters of the HAT-P-11 planetary system are reported in Table 1.

3 Atmospheric characterisation of HAT-P-11 b at high spectral resolution

3.1 Observations and data reduction

Four transits of HAT-P-11 b were simultaneously observed with the GIANO-B (wavelength range: 950-2450 nm, spectral resolving power R ≈ 50000) and the HARPS-N (wavelength range: 383-693 nm, spectral resolving power R ≈ 115 000) high-resolution spectrographs in the GIARPS at TNG configuration (Claudi et al. 2017) during the following nights: 7 July 2019; 18 June 2020; 19 September 2020; 13 June 2023. We only used the NIR (GIANO-B) observations for the present work. A total of 240 spectra were collected during the four observing nights (60 during the first one, 60 during the second one, 58 during the third one, and 62 during the fourth one), each with an exposure time of 200 s. The observations were performed with the nodding acquisition mode ABAB, where target and sky spectra were taken in pairs while alternating between two nodding positions along the slit (A and B) separated by 5″, allowing an optimal subtraction of the detector noise and background. All the observations were scheduled in order to obtain spectra before, during, and after the transit with airmass between 1 and 2. The measured mean signal-to-noise ratio (S/N) per spectrum, averaged across the entire spectral range and dataset, is between 48 and 59. In Table F.1 a schematic log of the observations is reported.

GIANO-B spectra cover the Y, J, H, K spectral bands in 50 spectral orders. The raw spectra were dark-subtracted and extracted using the GOFIO pipeline Python 3 version (Rainer et al. 2018). Although GOFIO also performs a preliminary wavelength calibration using U-Ne lamp spectra as a template, the mechanical instability of the instrument causes the wavelength solution to change during the observations. Since the U-Ne lamp spectrum is only acquired at the end of the observations to avoid persistence on the camera, the wavelength solution of the spectra determined by GOFIO is not sufficiently accurate and is expected to shift and jitter between consecutive exposures. In order to correct this shift, the spectra have been aligned to a common reference frame via cross-correlation with a time-averaged observed spectrum of the target used as a template. Thanks to this correction, we achieved a residual scatter in the measured peak position of the cross-correlation function (i.e. a residual shift of the spectra) well below 0.3 km s⁻¹ (approximately 1/10th of a pixel) for most of the spectral orders.

As these observations were performed from the ground, the spectra are contaminated by the presence of telluric lines (i.e. absorption lines due to the chemical species present in the Earth’s atmosphere). However, the telluric spectrum provides a good wavelength-calibration source, since the lines’ position does not change with time and the lines’ wavelength is well known. Refined wavelength calibration is made by matching a set of telluric lines in the time-averaged observed spectrum with a high-resolution model of the Earth transmission spectrum generated via the ESO Sky Model Calculator (Noll et al. 2012), and solving for the pixel-wavelength relation with a fourth-order polynomial fit. The spectral orders that showed either heavily saturated telluric lines or a high residual drift (> 0.4 pixels) have been excluded from the rest of the analysis. In particular, the excluded orders are: 8, 9, 10, 23, 24, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49 (in the GIANO-B spectra, order 0 is the reddest and order 49 is the bluest).

Fig. 7 Radial-velocity data analysis. Left panels: HIRES radial velocities of HAT-P-11 (black circles) showing a clear long-term trend likely due to the stellar magnetic activity cycle. The best-fit models including both the Keplerian signal of HAT-P-11 b and Gaussian-process regression with squared-exponential (top; Eq. (3)), quasi-periodic (middle; Eq. (4)), and quasi-periodic+squared-exponential (bottom; Eq. (5)) kernels are indicated with red solid lines. Right panels: HAT-P-11 RVs (black circles) phase-folded with the transit ephemeris of HAT-P-11 b, after removing the long-term trends modelled with Gaussian processes and the same kernels as in the corresponding left panels. The Keplerian eccentric model is displayed with a red solid line. We note the smaller and smaller scatter in the residuals as we move from top (SE kernel) to bottom (QP+SE kernel).

3.2 Telluric and stellar spectra removal

At this stage of the analysis, the planet spectrum is overshadowed by the stellar and telluric spectra. However, the planet’s orbital velocity has a non-zero radial component during transit (υ_p,⊥ ≈ 10 km s⁻¹); consequently, while telluric and stellar lines are stationary or quasi-stationary (the stellar barycentric radial velocity changes by few m s⁻¹ during transit events) in wavelength, the planet spectrum experiences a detectable change in Doppler shift during the ~140 min of transit. This property can be used to disentangle the planetary signal from the stationary components that we have to remove.

In this work, in order to remove telluric and stellar spectra, a principal component analysis (PCA) was conducted, after having masked the deepest absorption lines. The idea behind the PCA technique, which has been successfully applied in the past by several authors (e.g. de Kok et al. 2013; Damiano et al. 2019; Piskorz et al. 2017; Giacobbe et al. 2021; Guilluy et al. 2022; Carleo et al. 2022), is to identify common trends in the spectra as a function of time (in this case represented by telluric and stellar lines in different spectral channels) and remove them. Before computing the PCA, we performed an optimal selection of the spectral orders, for each molecule and for each night, to discard the orders that do not contain enough signal (molecular lines) and/or are strongly contaminated by telluric and stellar lines, following the method explained by Giacobbe et al. (2021). The PCA was performed on a total of K_i data matrices per night (where K_i is the number of selected spectral orders for the ith night). For each (M-rows; N-columns) data matrix (M = 58-62 images, N = 2048 pixels), the covariance matrix of the data was computed. Then, the covariance matrix was diagonalised by computing its eigenvectors (the Principal Components) ordered by their contribution to the variance of data (represented by the value of their associated eigenvalues). Following the procedure described in detail by Giacobbe et al. (2021), after having selected the appropriate number of principal components (between 9 and 23, depending on the quality of the night and the spectral order) that are supposed to mainly describe the telluric and stellar contaminations, these signals were reconstructed and removed from the data. The stages of the GIANO-B data reduction process are described with an example in Appendix A, while the details of the spectral orders’ selection and PCA procedures can be found in Appendix D.

3.3 Planet signal extraction via cross-correlation

As the faint planetary signal is very dispersed by the highresolution spectrograph, the single spectral lines are below the noise value (S/N_line ≲ 1) of the residual signal obtained after the telluric removal. However, there are thousands of planetary lines observed at the same time in the large wavelength range of GIANO-B whose signal can be co-added, resulting in a boost in terms of S/N for the planet signal proportional to the square root of the number of lines observed, (Birkby 2018). The information contained in such a large number of lines can be combined by cross-correlating the residual data with template transmission spectra of the planet’s atmosphere.

In particular, different models of the atmospheric transmission spectrum of HAT-P-11 b have been simulated with PetitRADTRANS (Mollière et al. 2019). All the simulated models assumed an isothermal atmosphere at the equilibrium temperature of the exoplanet at the transit epoch, assuming a heating re-distribution efficiency of ε = 1 (T_eq = 699 K), as described in Sect. 2.1.2, and are computed between 10 bar and 10⁻⁸ bar in pressure. The models assume constant-with-altitude abundance (volume mixing ratio) profiles. The different models, one for each molecule that we wanted to test the presence of, assumed fixed values of the volume mixing ratios (VMR) of molecular hydrogen = 0.855 and helium VMR_He = 0.145 and assumed a VMR_molecule = 10⁻³ for the molecule to test (single-species models). Although these values do not match any specific chemical scenarios, this was the simpler framework we could adopt to probe the presence of a particular molecule. We investigated the presence of eight molecular species, very common in the atmosphere of hot giant planets: water vapour (H₂O), methane (CH₄), ammonia (NH₃), acetylene (C₂H₂), hydrogen cyanide (HCN), carbon monoxide (CO), carbon dioxide (CO₂), and hydrogen sulfide (H₂S).

To probe the presence of a particular molecule, after having performed the PCA, we computed the cross-correlation function (CCF) between the data and the template spectrum associated with the molecule. The CCF was evaluated shifting the model in wavelength on a fixed grid of RV lags from −270 km s⁻¹ to +270 km s⁻¹, in steps of 0.1 km s⁻¹. The numeric computation was performed using the C_CORRELATE PXY(L) IDL function⁶, with null lag (L = 0), since the RV lags were applied to the wavelength array associated to the model and then the model was interpolated (via spline interpolation) on the same wavelength array of data before computing the CCF (i.e. with lag = 0). For every night and exposure, the CCFs calculated for each selected spectral order are co-added to obtain a single CCF for each exposure of each night.

Thanks to the high-resolution spectroscopy technique, it is possible to measure the Doppler shift of the spectral lines due to the planet’s orbital motion. In this way, in order to be sure that a particular spectral feature is produced by a molecule in the atmosphere of an exoplanet, the signal should have a Doppler shift that ‘follows’ the planetary movement and therefore we should observe that the peak of the CCF moves in wavelength as time passes according to the planetary motion-induced Doppler shift. We assumed that the measured Doppler shift of the planetary spectral lines is made of three-velocity components: (6)

where V_bary is the velocity induced by Earth’s motion around the barycentre of the Solar System (barycentric velocity), V_sys is the centre of mass velocity of the star-planet system with respect to the Earth (systemic velocity), and V_p is the planet RV. The time-dependent contribution V_p can be expressed as a function of two of the planetary orbital parameters (i.e. the eccentricity e and the argument of periastron ω_p), the RV semi-amplitude K_p and the true orbital anomaly ν(t): (7)

The radial velocity semi-amplitude K_p can be expressed as a function of the orbital eccentricity e, orbital inclination i, orbital period P_orb, mass of the planet M_p and mass of the star M_⋆: (8)

For our analysis, it is convenient to re-express K_p isolating the term containing the eccentricity and grouping all the others in the constant : (9)

In Fig. 8, we report, as an example, a plot of the CCF values as a function of the planetary orbital phase computed between the models of H₂O and CO and the data of the second observing night (18 June 2020). As it can be seen, the CCF peak trail is not visible by eye, since also at this stage of the analysis the planetary signal is too faint. This kind of plot helps in checking the telluric (stellar) spectrum-removal procedure since any strong residual signal due to a non-optimal telluric (stellar) subtraction would produce a spurious vertical trail of CCF peaks at a radial-velocity RV = 0 km s⁻¹ (RV = V_sys, neglecting the few m s⁻¹ motion induced by the planet onto the star) when data are cross-correlated with the H₂O (CO) templates. A benefit of the eccentric orbit of HAT-P-11 b is that the planetary radial velocity during the transit always remains strictly negative (~ 30 km s⁻¹ at the transit midpoint). This, combined with the high V_sys, shifts the planetary signal at ~100 km s⁻¹ far from the signal of the telluric lines (that is at RV = 0 km s⁻¹), further reducing spurious contaminations due to the Earth’s atmosphere.

If in the atmosphere of an exoplanet a molecule is present, its spectral signal has a null Doppler shift measured in the exoplanet rest frame (V_rest = 0 km s⁻¹), in the absence of a Doppler shift induced by atmospheric dynamics. It follows that, after having subtracted the barycentric, systemic and planetary RV from the CCF trails (that is, after having moved to the exoplanet rest frame), all the CCF peaks should align at V_rest = 0 km s⁻¹.

By subtracting different orbital solutions (in this work we explored a range of K_p values) from the CCF peaks trail, a different alignment of the CCF peaks is obtained. By co-adding the CCF values in phase (only considering the in-transit phases) into a single CCF for each trial K_p, the planetary signal as a function of the rest-frame velocity V_rest and K_p is maximised and it is possible to build the so-called K_p – V_rest map, in which, in the case of the detection of a molecule, a strong peak of the signal at the expected K_p and V_rest = 0 km s⁻¹ is observed. We took advantage of the high sampling of the CCF (larger than the GIANO-B pixel scale of 2.7 km s⁻¹) for a precise shift of the CCF trails into the planetary rest frame. However in order to avoid the use of correlated data points in our analysis, we binned the CCF values in radial velocity using a bin width of 2.7 km s⁻¹. We took the median of the 27 CCF values in each radial velocity bin as the value of the CCF associated with each bin, before co-adding the CCF values in phase for each trial K_p.

In this work, we explored a range of K_p values by varying the value of between 0 km s⁻¹ and 200 km s⁻¹ in steps of 2.7 km s⁻¹, having fixed the eccentricity and the other orbital parameters to the values reported in Table 1. This means that the corresponding explored range of K_p is [0; 207]km s⁻¹ (this is the K_p range reported in the K_p – V_rest maps in the next section and it is computed from using Eq. (9) with the value of e from Table 1). From Table 1, the expected value of K_p is = 123.4 ± 9.9 km s⁻¹ and, consequently, the expected value of is = 119.3 ± 9.5 km s⁻¹. However, the exploration of a large parameter space offers a strong diagnostic on all sources of noise and allows us to verify that no other spurious signal produces a significant detection near the planet’s rest frame position.

Fig. 8 Examples of CCF values as a function of planetary orbital phase computed with data from the second observing night (18 June 2020) and model containing only H2O (top panel) or CO (bottom panel) lines. The horizontal dashed lines represent the transit ingress and egress while the dash-dotted line represents the expected CCF peak trail due to the planetary motion as measured in the observer rest frame. The expected CCF peak trail in transit is not represented for clarity. As it can be seen, due to the faintness of the signal, the CCF peak trail is not visible by eye. This kind of plots serves as a visual check of any remaining telluric and stellar residuals, in this case showing no residuals and signifying that these are adequately corrected by the PCA.

Fig. 9 Signal-to-noise ratio Kp – Vrest maps for the probed chemical species: H2O, CH4, NH3, C2H2, HCN, CO, CO2, and H2S. Each Kp – Vrest map shows the S/N of the cross-correlation of the GIANO-B spectra (4 transits combined) with isothermal atmospheric models, as a function of the planet’s RV semi-amplitude (Kp) and the planet’s rest-frame velocity (Vrest). The S/N is computed by dividing the peak value of the cross-correlation function at each Kp by the standard deviation of the noise far from the peak, as described in the text. Negative S/N values correspond to anti-correlation. The vertical and horizontal white dashed lines correspond to Vrest = 0 km s−1 and the expected Kp value (), respectively.

3.4 Results

We first co-added the data of all the 4 observing nights and built the signal-to-noise ratio (S/N) K_p – V_rest maps for the different probed molecules, obtained by cross-correlating data with models and dividing the result by the standard deviation of the noise far from the peak (|V_rest| > 25 km s⁻¹), to search for signals following the expected planetary RV (potential detections). Then we computed the significance of each potential detection by performing a Welch t-test (Welch 1947) on two samples of CCF values: the former far (|V_rest| > 25 km s⁻¹) from the planet’s restframe velocity (“out-of-trail”) and the latter near to it (|V_rest| < 3 km s⁻¹, “in-trail”). The test rejects the null hypothesis (H₀) that the two samples have the same mean (and therefore that the CCF signal produced by the planet is only a statistical fluctuation of the background signal value) at a certain significance level that we adopted as the significance (σ) of our detections. Our significance calculations are based on the hypothesis of uncorrelated noise, which has been shown to be a valid approximation in previous works (Brogi et al. 2018; Guilluy et al. 2019). In order to build the significance K_p – V_rest maps, the Welch t-test was performed on CCF ‘in-trail’ and ‘out-of-trail’ distributions centred at the different V_rest explored, for each trial K_p.

In Fig. 9, we report the S/N K_p – V_rest maps for the different tested molecules, obtained by cross-correlating data with the different models and following the procedure explained in Sec. 3.3. As it can be seen, we obtain a signal around the expected planetary position in the K_p – V_rest maps with S/N > 3 for 4 molecular species: H₂O (S/N = 5.1), CH₄ (S/N = 4.8), NH₃ (S/N = 5.3), and CO₂ (S/N = 3.0).

In Fig. 10, we report the significance K_p – V_rest maps for the different tested molecules, obtained by performing the Welch t-test on in-trail and out-of-trail CCF distributions. As it can be seen from the maps, the peak of the signal for the four potentially detected species has a significance of 3.4 σ (H₂O), 2.6 σ (CH₄), 5.0 σ (NH₃), and 3.2 σ (CO₂). For the other four probed species, we measure no significant chemical signature at the planetary RV, and therefore we consider them as non-detections and focus our attention and following analyses on the signals of H₂O, CH₄, NH₃, and CO₂ only. The distributions of CCF values in-trail and out-of-trail used for computing the significance of the detections via the Welch t-test are reported in Fig. B.1 for all the probed species.

It is interesting to notice the orientation of the signals in the K_p – V_rest maps, which is typically vertical for atmospheric transmission studies (e.g. Giacobbe et al. 2021) and sloped for atmospheric emission studies (e.g. Line et al. 2021). This is due to the eccentric orbit of HAT-P-11 b that makes the planetary RV to be not symmetrical around 0 km s⁻¹ during the transit (it is always strictly negative), as it happens for planets on circular orbit observed in emission (i.e. at orbital phases different from 0 or 0.5).

Finally, for the four selected species, we computed how much the Doppler signature of the signals is in accordance with the expected planetary RV. To do this, we built contour plots of the detection significance, defining the 1 σ and 2 σ significance intervals as the regions of the K_p – V_rest maps where the significance drops by 1 σ and 2 σ, respectively, with respect to the significance peak, and looked at where the expected position of the atmospheric signal in the K_p – V_rest maps was with respect to those intervals. In Fig. 11, we report the contour plots of the detection significance for the four chemical species. As it can be seen, the signals we measure have a Doppler signature compatible with the planetary one at < 1 σ, for all four molecules.

Before claiming any detection, we checked the reliability of our results by performing a further test, described in the following. For each of the four chemical species, we combined the matrices of the CCF as a function of the orbital phase (see, e.g. the one in Fig. 8) of the four observing nights in a single CCF matrix with the rows sorted in crescent orbital phase. Then, we randomly shuffled the CCF order in phase (this corresponds to shuffling the sequence of observed spectra in time, including the out-of-transit ones) 250 times. For each shuffle we built both the S/N and the significance K_p – V_rest maps in a restricted interval of K_p = [89.4; 156.5] km s⁻¹ (corresponding to ± 34 km s⁻¹) and V_rest = [−10; 10] km s⁻¹, in order to test the presence of spurious signals that have not a planetary origin but can produce significant features around the expected signal position in the maps due to some peculiar time-correlated noise. We repeated this test for the four selected chemical species and for each species we studied the distributions of the 43 750 values of S/N and 43 750 values of significance obtained in the chosen K_p – V_rest interval combining the 250 permutations.

The results of this test are reported in Table 4. As it can be seen, for all the 4 molecular species, 95% of the test yields signals with S/N < 3 and significance ≤2.5 σ in the selected K_p – V_rest interval, while 99.73% of the test yields S/N ≤ 4.6 and significance ≤4.2 σ, in the same interval.

In Figs. C.1 and C.2, we report the distributions of S/N and significance values obtained in this test, respectively, for the four selected chemical species.

After having performed this test, for each of the four molecular species we computed the probability (p-values) of randomly drawing from the distributions reported in Figs. C.1 and C.2, the S/N and significance values that we measure. To compute the p-values for the S/N (significance) value, we summed the occurrences of S/N (significance) greater than the S/N (significance) measured and divided them by 43 750. We report the p-values in Table 5. As it can be seen, for all the four molecular species we obtain a p-value ≤2.5% for the S/N and ≤3.8% for the significance, so there is less than 5% of probability that these signals are due time-correlated noise.

In Table 5, we summarise the significance level of the signal of the four molecular species of interest obtained with different statistical methods. From these results, we conclude that we have a statistically robust detection of NH₃, which is the most significant signal (5 σ) and the one with the lowest p – values (< 0.04%). We consider the H₂O as a detection, even though it is less statistically robust than the NH₃ one, since it has already been detected at low-resolution multiple times in the atmosphere of HAT-P-11 b (see the literature cited in Sec. 1) and therefore we can be more confident about the reliability of the signal we measure. Since the t-test significance for the signal of CH₄ is <3 σ, the significance K_p – V_rest map presents a second comparable peak at V_rest ≈ 20 km s⁻¹, and the associated p-value is the highest among the four chemical species, even if the signal has a S/N = 4.8 and a Doppler signature that is compatible with the planetary one at <1 σ, the detection of this molecular species remains tentative as it is not sufficiently robust from a statistical point of view. Even though the signal from CO₂ has both the S/N and significance ≥ 3 σ and the RV trail is compatible with the planetary one at <1 σ, it has the lowest S/N and the highest S/N p-value among the four selected species. For these reasons, we consider the detection of CO₂ as tentative too. As it can be seen, for what concerns these two tentative detections, even though we obtain signals at the expected planetary RV, they are not sufficiently statistically robust and we suggest conducting further studies to unambiguously assess the presence of CH₄ and CO₂ in the atmosphere of HAT-P-11 b.

To further assess the robustness of our analysis for what it concerns the impact of the number of principal components removed by the PCA (changing with the spectral order and night) on the final results, we also repeated the whole analysis removing a fixed number of principal components with the PCA for all the selected spectral orders per-molecule, an for all the observing nights. We performed this test twice, the first time we removed 9 principal components (the minimum number of principal components removed in our work among the different spectral orders and nights), the second time we removed 23 components (the maximum number of principal components removed in our work among the different spectral orders and nights). We obtain that none of these 2 extreme conditions changes our interpretation of which chemical species we detect (i.e. H₂O and NH₃), and which chemical species we tentatively detect and need further investigations (i.e. CH₄ and CO₂), even if their S/N and significance slightly change (at less than 1 σ level), as expected.

Finally, to assess the robustness of the final results (in particular the detection of H₂O and NH₃) in relation to the impact of the spectral orders’ selection procedure on the CCF analysis, we performed an additional test. In this test, we repeated the analysis refining the orders’ selection procedure to assess the presence of possible spurious signals near the expected planetary radial velocity that could lead to possible false positive detections. The description of this test and the results are reported in Appendix E. We do not observe relevant changes in the results, further confirming the conclusion we reached with the main analysis.

In summary, in this work, we report the detection of two molecular species (i.e. H₂O and NH₃) in the atmosphere of HAT-P-11 b and the tentative detection of two others (i.e. CH₄ and CO₂), whose presence has to be assessed by further studies.

Fig. 10 Significance Kp – Vrest maps for the tested chemical species: H2O, CH4, NH3, C2H2, HCN, CO, CO2, and H2S. Each Kp – Vrest map shows the significance of the cross-correlation signal of the GIANO-B spectra (4 transits combined) with isothermal atmospheric models, as a function of the planet’s RV semi-amplitude (Kp) and the planet’s rest-frame velocity (Vrest). The significance is computed with a Welch t-test on two samples of cross-correlation values, that is, far from and near to the planet’s RV respectively (see the text for more details). The vertical and horizontal white dashed lines correspond to Vrest = 0 km s−1 and the expected Kp value (), respectively.

Fig. 11 Contour plots of the detection significance of the four selected chemical species. The solid (dashed) lines represent the 1 σ (2 σ) interval around the peak value of the significance (marked with a red cross). These intervals are computed as described in the text. The point of the Kp – Vrest map in which the 2 black dashed lines (the horizontal one corresponds to the expected Kp value, while the vertical one corresponds to Vrest = 0 km s−1) cross each other represents the expected detection significance peak position in the case in which the detected signal has a planetary origin. As it can be the significance peak position is compatible with the planetary origin hypothesis at better than 1 σ for all the 4 species.

3.5 Discussion

The significance of the detections reported in this work is lower than that of the species detected in other hot and warm Jupiters’ atmospheres, which even reached the 10σ level (e.g. Giacobbe et al. 2021). This is mainly due to the lower atmospheric signal level of warm Neptune-size exoplanets (the atmospheric signal level is , where H_s is the atmospheric scale height, and therefore it is ~55 ppm⁷ for HAT-P-11 b) and underlines the difficulty in probing the atmospheric features of this class of exoplanets.

3.5.1 Detected species in the HAT-P-11 b planetary context

Our detection of water vapour in the atmosphere of HAT-P-11 b is in accordance with the results of Fraine et al. (2014), Tsiaras et al. (2018), Chachan et al. (2019), and Cubillos et al. (2022). The signal from the methane, which we measure at the expected planetary radial velocity, supports the results obtained by Chachan et al. (2019) and Cubillos et al. (2022), who suggested the presence of CH₄ from the analysis of HST data however, as told in the previous section, the significance of the signal is very low and therefore we cannot confirm the presence of this molecular species in the atmosphere of HAT-P-11 b and further studies are required.

The hot and warm Neptunes’ atmospheric chemistry and dynamics are different from the ones that act in the hot and warm Jupiters’ atmospheres, mainly because of the smaller radius and mass and the possible higher metallicity of the former (Moses et al. 2013). The atmospheric composition of individual exoplanets also depends on their effective temperature, formation history, atmospheric evolution, orbital parameters, and irradiation environment, so it is difficult to predict their exact atmospheric properties. However, for what concerns the chemical composition, some general trends with temperature, metallicity and C/O ratio can be found, as shown in the works of Moses et al. (2013), who particularly focused on hot Neptunes, Madhusudhan (2012) and Fonte et al. (2023).

Even though with this work we are not able to constrain HAT-P-11 b atmospheric physical and chemical properties (such as the elemental abundances), we can still make some general considerations about our detections. In particular, under thermochemical equilibrium, the presence of hydrocarbons, like HCN and C₂H₂, is particularly favoured in carbon-rich environments (C/O ≳ 1) at high temperatures (T ≳ 1000 K). This is in line with our non-detections of these two chemical species, given the relatively low-temperature atmosphere (T_eq ~ 700 K) of the target we analysed. Chachan et al. (2019) measured an atmospheric C/O ratio close to unity . A C/O value close to unity together with a low-temperature environment favours the formation of ammonia, which is the molecule that we detected with the highest significance and the second nitrogen-bearing species we probed. In addition, at temperatures lower than T ≲ 1300 K, the formation of H₂O and CH₄ is favoured, in accordance with our detection of H₂O and tentative detection of CH₄, considering the HAT-P-11 b equilibrium temperature.

Another interesting result of our work is the tentative detection of CO₂ in the atmosphere of HAT-P-11 b. Carbon dioxide is an important indicator of the metal enrichment of the atmosphere of exoplanets and therefore it can give important information about the formation processes of the primary atmospheres of gas giants. Indeed several models for warm gas giant atmospheres predict that CO₂ is one of the molecules most strongly enhanced with increasing atmospheric metallicity, becoming detectable for metallicities greater than ≈10 times that of the Sun (Lodders & Fegley 2002; Moses et al. 2013). Until now, only JWST Transiting Exoplanet Community Early Release Science Team et al. (2023) have firmly detected the presence of CO₂ in the atmosphere of an exoplanet (the warm Saturn WASP-39 b) at low spectral resolution, while at high spectral resolution Carleo et al. (2022) tentatively detected it in the atmosphere of the warm Jupiter WASP-80 b. Therefore, our tentative detection of CO₂ represents the first hint of the presence of this molecule in the atmosphere of a warm Neptune and, if confirmed, could point towards a super-solar metallicity for the atmosphere of HAT-P-11 b. However, as underlined before, further studies are needed and only a statistical comparison between different atmospheric chemical-physical models (for example through atmospheric retrievals), could robustly determine the atmospheric chemical-physical characteristics of this target.

For the case of HAT-P-11 b, the estimated metallicity is Z < 4.6 Z_⊙ at 2 σ level (Chachan et al. 2019) with Z_⊙ the solar metallicity, that is, Z < 2.3 Z_⋆ referred to the host-star metallicity Z_⋆ = 2.0 ± 0.2 Z_⊙ (Bakos et al. 2010), while Welbanks et al. (2019) estimated a substellar H₂O/H metallicity. At this level of metallicity, CO₂ should be scarce (see Moses et al. 2013, Fig. 5), but if we consider the 3 σ confidence level upper limit on metallicity (Z < 86 Z_⊙ = 43 Z_⋆) by Chachan et al. (2019), the CO₂ abundance could increase up to the point of being detectable with our kind of analysis.

It is worth noting that the cross-correlation analysis is more sensitive to molecular lines’ position in wavelength rather than their depths with respect to the continuum. This means that a more significant detection for a particular chemical species does not necessarily imply that species is more abundant with respect to the other ones because it could arise for example from a denser forest of lines and consequently from a stronger correlation peak. In addition, our non-detections do not necessarily imply the absence of such chemical species, for which further investigations are needed. Finally, for what concerns the non-detections, the K_p – V_rest significance maps show spurious signals that are not related to the planetary signal. They are typically aliases generated by the autocorrelation function of the template that can also be seen in the maps of the detected species far from the planetary RV.

3.5.2 Atmospheric chemical modelling

In order to better characterise the atmospheric properties of HAT-P-11 b, we tried to find the chemical scenarios that are most compatible with the detection of H₂O and NH₃, and possibly CH₄ and CO₂ (while also considering a non-detection of the other species) in the atmosphere of HAT-P-11 b, by exploring with a grid of theoretical models the possible radiative state and elemental compositions. To this end, we employed the Pyrat Bay (Cubillos & Blecic 2021) modelling framework to compute transmission spectra from an atmosphere in radiative and thermochemical equilibrium. The code iterates over a two-stream radiative-transfer calculation until the atmosphere converges towards a stable radiative equilibrium solution at each layer (Heng et al. 2014; Malik et al. 2017). The end products of this approach are the temperature and composition profiles of an atmosphere in radiative equilibrium, from which we can produce transmission-spectrum models to determine which species have observable spectral features.

The inputs of the models are the known system parameters, a stellar spectrum (Castelli & Kurucz 2003), and the atmospheric elemental composition. We explored a range of plausible scenarios by varying the elemental composition over a range of metallicities (from [M/H] = −1.0 to 2.0), C/O ratios (from 0.1 to 1.5), N/O ratios from (0.14 to 0.85), and over a range of heat radiation regimes by varying the β_irr = (1 – A_b)/f parameter (from β_irr = 0.5 to 1.0), where A_b is the Bond albedo and f is the day-night heat redistribution efficiency.

The atmospheric model spans a pressure range from 100 to 10⁻⁹ bar, and a wavelength grid ranging from 0.3 to 30 μm sampled at a resolving power of R = 15000, sufficient to contain the bulk of the stellar and planetary fluxes. The chemical network includes 45 neutral and ionic species that are the main carriers of H, He, C, N, O, Na, Si, S, K, Ti, V, and Fe. The opacity sources include line-list data for CO, CO₂, and CH₄ from HITEMP (Rothman et al. 2010; Li et al. 2015; Hargreaves et al. 2020), and H₂O, HCN, NH₃, and C₂H₂ from ExoMol (Polyansky et al. 2018; Chubb et al. 2020; Yurchenko et al. 2011; Harris et al. 2006, 2008; Coles et al. 2019). We preprocessed these large data sets with the REPACK package (Cubillos 2017) to extract the dominant line transitions. Additionally, we included Na and K opacities (Burrows et al. 2000); H, H₂, and He Rayleigh opacities (Kurucz 1970); H₂–H₂ and H₂–He collision-induced absorption (Borysow et al. 1988, 1989, 2001; Borysow & Frommhold 1989; Borysow 2002; Jørgensen et al. 2000), and H⁻ free-free and bound-free opacity (John 1988). Once we obtained a grid of radiative-equilibrium atmospheric models, we post-produced transmission spectra at the GIANO-B spectral resolution.

In general, our grid of HAT-P-11 b models indicates that the strong H₂O absorption bands dominate the transmission spectra at most wavelengths. The scenarios that are more consistent with the detected species required super-solar metallicity, since this enhances the CO₂/H₂O abundance ratio, which in turn makes CO₂ detectable. Similarly, the lower-heat models improve the detection of both NH₃ and CH₄, since both of these species are more abundant than at higher temperatures. However, NH₃ is never the dominant nitrogen bearer at these pressure and temperature conditions. Consequently, only for the highest N/O ratios tested (N/O=0.85) the NH₃ absorption lines impact the transmission spectra (Fig. 12), this is well above the solar N/O ratio of 0.14 (Asplund et al. 2021). Alternatively, one can invoke disequilibrium-chemistry processes to enhance the abundance of NH₃ with respect to the other species (Moses 2014). Specifically, vertical quenching can transport NH₃ from deeper layers where it is significantly more abundant (Fig. 12, bottom-left panel). Combining these constraints with future observations of HAT-P-11 b (e.g. JWST GO proposals 2950 and 4150) has the potential to place stronger constraints on the planet’s atmospheric composition and further infer which physical processes shape its atmospheric properties.

The right-hand panels in Fig. 12 also highlight the complementary potential of combining low- and high-resolution observations for atmospheric characterisation. For a molecule to be detectable, its absorption lines must be at least at the same level as the dominant species in the atmosphere; in this case, they overlap with the H₂O lines, as shown by the coloured curves on the right panels of Fig. 12. While this can be the case for several species at high resolution, instead, at low resolution (black curves) many of these individual line features are washed out and blend into the absorption of the dominant species. Therefore, with sufficient S/N, high-resolution observations can pick up molecular features that may go undetected at low resolution.

Fig. 12 Most favourable atmospheric chemical models. Top-left panel: HAT-P-11 b atmospheric temperature-pressure profile corresponding to the most favourable chemical-equilibrium model that shows spectral features from all the detected species, including the tentatively detected ones (Model 1: [M/H]=2.0, C/O=0.9, N/O=0.85, βirr=0.5; solid line) and to a grid model modified including NH3 vertical quenching at the 10 bar level (Model 2: [M/H]=1.7, C/O=0.59, N/O=0.14, βirr=0.5; dashed line). Bottom-left panel: volume mixing ratios of the dominant chemical species (see legend) as a function of the atmospheric pressure corresponding to Model 1 (solid lines) and to Model 2 (dashed lines). Right panels: theoretical transmission spectra of HAT-P-11 b for atmospheres in thermochemical equilibrium (top-right panel) and with NH3 quenching (bottom-right panel). The coloured curves show the contribution to the synthetic spectra of the four species detected by GIANO-B, including the tentatively detected ones (see legend). The black curve shows a lower-resolution (R = 500) spectrum combining the absorption from all atmospheric species in the model.

3.5.3 How the possible presence of planet c could have influenced the formation and the atmospheric composition of planet b

The study of the atmospheric chemical composition of close-in sub-Neptunes and Neptunes, such as HAT-P-11 b, can give important constraints about giant-planet formation in the pebble-accretion scenario. In particular, in this scenario, the formation of these kinds of planets in the inner part of the protoplanetary disc depends on the flux of water ice-rich pebbles drifting inward from the outer regions of the protoplanetary disc (Bitsch et al. 2022). At the water ice line, the ice on the pebbles evaporates enriching the inner part of the disc with water vapour. Substructures and cavities in the protoplanetary disc, for example, due to the formation of a Jupiter-size exoplanet, create pressure bumps that might block the inward migration of pebbles (Banzatti et al. 2020). In this way the presence of a rapidly growing external Jupiter-size companion and its position with respect to the water ice line could influence the chemical composition of the inner protoplanetary disc and therefore of the close-in Neptune: if the giant outer companion forms outside the water ice line it could block the water-rich pebbles before they cross the water ice line and therefore it prevents them from enriching the inner disc with water; if the giant outer companion forms inside the water ice line, the inward migrating pebbles might be blocked after they cross the water ice line and their ice evaporates enriching the inner disc with water. In the first case the close-in Neptune would form in a “dry” environment, in the second case it would form in a “wet” environment accreting a water-rich (up to a few per cent) atmosphere (Bitsch et al. 2022). A third scenario is also possible: the close-in Neptune forms beyond the water ice line by the accretion of a large amount of water ice before migrating inwards. In this case, the Neptunian planet would contain up to 50% of its total mass in water ice (’very wet’ Neptune) and the formation time and location of the giant companion would not matter. Also carbon-bearing species can be blocked by growing giant planets, reducing the metal enrichment of the inner disc and therefore the atmosphere of the close-in Neptune. Of course, it is important to underline that protoplanetary discs evolve over time by cooling down, shifting the ice lines inwards over time and therefore both the time and the position of the forming planets inside the protoplanetary disc play an important role in this scenario (e.g. Eistrup et al. 2016, 2022; Eistrup & Henning 2022).

In addition, a late accretion of planetesimals or planet-planet scattering could further influence the atmospheric composition of the Neptunian planet. For example, the high eccentricity of the orbit of HAT-P-11 b could be explained by the interaction with a possible external companion (if any) through planet-planet scattering migration. During the migration, HAT-P-11 b might have been enriched in heavy elements mainly through the accretion of planetesimals present in the protoplanetary disc instead of pebbles (blocked by the external companion). In the orbital region populated by HAT-P-11 b, planetesimals are more effective in delivering O than N to the accreting planet (Turrini et al. 2021). As their accretion acts to lower the high N/O ratio of giant planets formed by pebble accretion, this scenario favours the second atmospheric model we considered in Sect. 3.5.2. As it can be seen, the presence of an external companion could have played a crucial role in the HAT-P-11 b formation, influencing its atmospheric composition, as also highlighted by Chatziastros et al. (2024). Of course, the investigation of the presence of an outer giant companion and its characterisation, as well as the retrieval of the elemental abundance and the elemental ratios (e.g. C/O, N/O) of the chemical species present in the atmosphere of HAT-P-11 b, are crucial to probe the formation history of the planetary system and to enlarge our knowledge about the planetary formation process.

3.5.4 Planetary radial-velocity semi-amplitude shift investigation

As can be seen from the K_p – V_rest maps and from the contour plots, even though the signals of all the molecules for which we obtain a detection (including the tentative ones) cross the expected K_p and V_rest values, they span a large range of K_p and V_rest values, of the order of ≈ 100 km s⁻¹ and ≈ 30 km s⁻¹, respectively. The possible presence of winds and the planetary rotation, combined with other complex atmospheric dynamics effects involving global circulation⁸, introduces distortions of molecular line profiles during the transit event that produce shifts of the detection signal peak in K_p and V_rest (e.g. Wardenier et al. 2021; Kesseli et al. 2022; Pelletier et al. 2023). In particular, for the case of HAT-P-11 b, Allart et al. (2018) found evidence for a high-altitude wind flowing from the day-side to the night-side of the planet at a velocity of υ_wind ≈ 3 km s⁻¹. For what concerns the planetary rotation speed, due to its eccentric orbit, the most probable rotation configuration for HAT-P-11 b is the pseudo-synchronous one (Hut 1981), in which the planet rotation is synchronous at periastron, implying a rotation period P_rot that is shorter than the orbital one by about 29% (for an eccentricity of e = 0.2577). In the case of HAT-P-11 b, this rotation period would correspond to an equatorial rotation velocity of υ_eq ≈ 1 km s⁻¹. As it can be seen, these kinds of effects have typical velocity scales of a few km s⁻¹, which is about one order of magnitude smaller than our precision. This implies that we are not sensitive to them. With future more precise measurements, it will be possible to investigate the complex atmospheric dynamics of this warm Neptune.

Even though our results do not allow us to probe the presence of atmospheric dynamics effects, we observe (see Fig. 11) a systematic shift of the signal peak at K_p values smaller than the expected one (with the exception of the signal of NH₃ that is at the expected K_p). In order to quantify this shift, we made a weighted average of the K_p values corresponding to the peak of significance for each of the detected molecules (including the tentatively detected one), assuming as the extremes of the 1 σ uncertainty intervals the K_p values where the significance drops by 1 σ from the peak (after having marginalised over the V_rest values). The averaged value of the measured K_p is: K_pmeas = 107 ± 29 km s⁻¹, which is compatible at 0.54 σ with the expected value reported in Table 1.

We note that the large uncertainty in K_p, due to the relatively small change in planet RV during transit (a few tens of km s⁻¹), makes the systematically observed downward shift in K_p of the order of ≈ 10 km s⁻¹ not significant on a statistical basis. However, we decided to investigate the possible presence of a systematic effect by studying how the uncertainties on e and ω_p affect the retrieved value of K_p and therefore the position of the peak of significance in the K_p – V_rest maps. These two parameters are the two main sources of uncertainties in determining the planetary orbital configuration and the atmospheric Doppler signature (Montalto et al. 2011).

We generated N = 500 000 couples of values of orbital eccentricity and planetary argument of periastron (e; ω_p), extracted randomly from two Gaussian distribution (one for e and one for ω_p) with mean equal to the values of these two parameters that we used in this work (ω = 0.2577; ω_p = 192.0 deg) and standard deviation equal to the 1 σ uncertainties on these values. For each (e; ω_p) couple we computed the planetary RV curve during the transit combining Eq. (7) and Eq. (9): (10)

fixing the value to the expected one ( = 119.3 km s⁻¹). In this way, we simulate different possible orbital configurations according to our limited knowledge of the (e; ω_p) parameters. Subsequently, we fitted the value for each of the N RV curves keeping fixed the (e; ω_p) values to the ones we used in this work.

Following this procedure, we obtained a distribution of retrieved values that we transformed into a distribution of K_p values (using the e value reported in Table 1) in order to compare them with the K_p value obtained from the significance maps. This distribution is reported in Fig. 13. With this procedure we simulated the effect of the choice of (e; ω_p) used in this work on the K_p value estimation from the peak of the CCF signal, taking into account our limited knowledge of the orbital solution.

The RV semi-amplitude value is strictly related to the slope of the planetary RV curve. In our data analysis, the slope of the planetary RV curve is not directly fitted but, since it influences the alignment of the CCF trails, it changes the position of the peak of significance in the K_p – V_rest maps and therefore the value of K_p we can estimate from the maps.

The first result worth noting is that the K_p distribution is asymmetric with a long tail extending towards K_p > 150 km s⁻¹ values. In terms of the properties of the distribution, its median value coincides with the expected value and by computing the 1 σ interval of this distribution, the estimated value of K_p can be obtained: . This is consistent with the expected value but, as it can be seen, the distribution has not a peak at the expected K_p value, instead, the modal value is: K_Pmode = 114 km s⁻¹. This means that the most probable measured K_p value is ≈ 10 km s⁻¹ smaller than the expected one, which is consistent with our measured K_pmeas value. From this investigation, two important considerations can be made: 1) even if the values of eccentricity and argument of periastron are very precise (≈1%), the K_p value that can be extracted from our kind of data analysis cannot be known at better than ≈25%, due to the relatively small change in planet RV during transit; 2) this investigation suggests the presence of a systematic effect that we do not fully understand that makes more probable to observe a shift of the signal of the order of ≈ 10 km s⁻¹ towards K_p values lower than the expected one, at least in the case of the analysed target, related to the uncertainties about the orbital solution adopted. This has to be taken into consideration for future atmospheric dynamics studies, since in that case precise but also accurate measurements are needed.

Fig. 13 Distribution of Kp values fitted on randomly generated planetary RV curves with properties described in the text. The vertical black solid line represents the median value of the distribution, while the 2 vertical black dashed lines represent the 2 quantile values delimiting the 68% (1 σ) interval of Kp around the median. The vertical dash-dotted blue line represents the expected Kp value: = 123.4 km s−1. The median and the expected values are coincident and therefore they are not visibly distinguishable.

4 Conclusion

In this work, we revisited the HAT-P-11 planetary system. Using Kepler and HIRES at Keck archival data, we refined the orbital and physical parameters of HAT-P-11 b. We further showed that the long-term RV signal with a semi-amplitude of 30 m s⁻¹ and periodicity of ~9-10 yr is more likely due to the stellar activity cycle than to the presence of the planet HAT-P-11 c (Yee et al. 2018). Nonetheless, the HIPPARCOS–Gaia difference in propermotion anomaly suggests that an outer-bound companion might still exist, though with S/N ≲ 5. The continuation of RV monitoring and, even more importantly, a combined analysis of RV and future DR4 Gaia astrometric data, will help to characterise this possible companion.

This review of the HAT-P-11 planetary system allowed us to perform a consistent analysis of the atmosphere of planet b at high spectral resolution. Moreover, HAT-P-11 b represents a remarkable target as it can provide a better understanding of the atmospheres of warm Neptunes, which remain to be explored in detail.

In particular, we probed the presence of eight chemical species in the atmosphere of HAT-P-11 b by cross-correlating template atmospheric models with data taken with the NIR GIANO-B high-resolution spectrograph at the 3. 58 m TNG telescope during four planetary transits. We detect the presence of two molecular species, H₂O and NH₃, with a S/N of 5.1 and 5.3, and a significance of 3.4 σ and 5.0 σ, respectively. The signals from these molecules have a Doppler signature compatible with the planetary one at <1 σ. We also tentatively detect the presence of CH₄, whose signal has a S/N of 4.8 but a significance level of 2.6 σ, and of CO₂, with a S/N of 3.0 and a significance level of 3.2 σ. Further studies are necessary to confirm the presence of CH₄ and CO₂ in the atmosphere of HAT-P-11 b. These results constitute the first simultaneous observation of multiple molecular species in the atmosphere of a warm Neptune-type planet.

Our results are in accordance with the previous H₂O detections made by Fraine et al. (2014), Tsiaras et al. (2018), Chachan et al. (2019), and Cubillos et al. (2022) and with the suggestion of the presence of methane by Chachan et al. (2019) and Cubillos et al. (2022). We hereby enlarge the number of chemical species known to be present in the atmosphere of this exoplanet, which allows us to start exploring plausible physical conditions for the atmosphere of HAT-P-11 b. Our models suggest two scenarios that are more in accordance with the observations: the first model describes an atmosphere in chemical equilibrium with supersolar metallicity and enhanced C/O and N/O ratios relative to solar values; the second model describes an atmosphere with disequilibrium chemistry (i.e. NH₃ vertical quenching), lower metallicity, and C/O and N/O ratios close to solar values.

We note that the significance peak of the detected species (including the tentatively detected ones) is shifted by ≈ 10 km s⁻¹ towards K_p values smaller than the expected one ( = 123.4 km s⁻¹). Indeed the average value of the planetary RV semi-amplitude estimated from our significance maps is K_Pmeas = 107 ± 29 km s⁻¹, which is compatible at 0.54 σ with the expected one due to the large uncertainty. We show how a small error (≈ 1%) on the eccentricity and argument of periastron parameters translates into a large uncertainty (≈25%) on the retrieved K_p parameter due to the small change in RV during the transit, and can introduce the observed systematic K_p shift.

The next step in our analysis will be to statistically retrieve the chemical physical properties of the atmosphere of HAT-P-11 b (e.g. the elemental abundances and the temperature-pressure profile) in order to better characterise its atmosphere and to constrain its formation path. Finally, a future combination of high-and low-resolution observations will improve the characterisation of the atmosphere of this target and our knowledge about warm Neptunes.

Acknowledgements

This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. The authors acknowledge the financial contribution from PRIN INAF 2019, and from the European Union – Next Generation EU RRF M4C2 1.1 PRIN MUR 2022 project 2022CERJ49 (ESPLORA) and project 2022J4H55R. A.S. acknowledges support from the Italian Space Agency (ASI) under contract 2018-24-HH.0 “The Italian participation to the Gaia Data Processing and Analysis Consortium (DPAC)” in collaboration with the Italian National Institute of Astrophysics. The following internet-based resources were used in research for this paper: the ESO Digitized Sky Survey; the NASA Astrophysics Data System; the SIMBAD database operated at CDS, Strasbourg, France; and the arχiv scientific paper preprint service operated by the Cornell University. A.M. acknowledges partial support from the ASI-INAF agreement no. 2018-16-HH.0 (THE StellaR PAth project). P.E.C. is funded by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship, program J4595-N.

Appendix A Stages of the GIANO-B data reduction process

Fig. A.1 Example of the stages of the GIANO-B data reduction process applied to real GIANO-B spectra over a short wavelength interval (x-axis). The planetary orbital phase is reported on the y-axis. From top to bottom panel: (a) the extracted spectra; (b) residual spectra after having normalised each spectrum (each row) by its median value to correct baseline flux differences between the spectra that are due, for example, to variable transparency of the atmosphere, imperfect telescope pointing, or instability of the stellar point spread function; (c) residuals after each spectral channel (each column) had its mean subtracted; (d) residuals after having masked the strongest or saturated telluric lines and having removed the telluric and stellar spectra with the PCA; as it can be seen also other spurious effects such as the alternating pattern visible at certain wavelengths likely related to the instrument A-B nodding are effectively removed with PCA; (e) residuals after having applied a high-pass filter to each row of the residual matrix (this step mitigates any residual correlation between different spectral channels).

Figure A.1 shows an example of the stages of our data reduction process applied to real GIANO-B spectra, over a short wavelength interval, from the extracted spectra to the residuals after the PCA application.

Appendix B Distribution of the cross-correlation function values

Figure B.1 shows the distribution of CCF values far (|V_rest| > 25 km s⁻¹) from the planet’s rest-frame velocity (“out-of-trail”) and near to it (|V_rest| < 3 km s⁻¹, “in-trail”) for the 3 detected molecular species. The Welch t-test used to estimate the significance of the detections was performed on these distributions.

Fig. B.1 Distribution of CCF values far (|Vrest| > 25 km s−1) from the planet’s rest-frame velocity (“out-of-trail”, OOT, dashed lines) and near to it (IVrest| < 3 km s−1, “in-trail”, INT, solid lines) for the probed chemical species.

Appendix C Test on the reliability of the detections

In order to check the reliability of the detections, we combined the matrices of the CCF as a function of the orbital phase of the four observing nights in a single CCF matrix with the rows sorted in the crescent orbital phase. We then randomly shuffled the CCF order in phase 250 times. For each shuffle, we built both the S/N and t-test significance K_p – V_rest maps in an interval of K_p = [89.4; 156.5] km s⁻¹ ( ± 34km s⁻¹) and V_rest = [−10; 10] km s⁻¹, in order to test the presence of spurious signals due to correlated noise that can produce significant features around the expected signal position in the maps. We repeated this test for all four detected chemical species for a total of 250 combinations, 43 750 values of S/N, and 43 750 values of significance for each species. In Fig. C.1 and Fig. C.2 we report the distributions of S/N and significance values obtained in the selected K_p – V_rest interval, respectively, for the 4 detected species.

Fig. C.1 Distribution of S/N obtained in the selected Kp – Vrest interval by shuffling 250 times the time order of the observed spectra, for each of the four selected chemical species. The vertical dashed (dash-dotted) lines represent the borders of the intervals enclosing the 95% (99.73%) of the S/N values, corresponding to the 2.5% − 97.5% (0.135% − 99.865%) quantiles of the distributions. These intervals are reported in Table 4.

Fig. C.2 Distribution of significance (from the t-test) obtained in the selected Kp – Vrest interval by shuffling 250 times the time order of the observed spectra, for each of the four selected chemical species. The vertical dashed (dash-dotted) lines represent the borders of the intervals enclosing the 95% (99.73%) of the significance values, corresponding to the 95% (99.73%) quantile of the distributions. These intervals are reported in Table 4. Since the 50% of the significance values (~20 000 values) are <0.2 σ for all the 4 chemical species, the plots of the distributions are limited on the y-axis in the interval [0; 650] for clarity.

Appendix D Orders’ selection and PCA procedure

In this section we briefly describe the orders’ selection procedure and the PCA procedure we applied to data.

We followed Giacobbe et al. (2021) for the spectral orders’ selection procedure, per molecule and per night. In particular, we injected single-species model planetary spectra at planetary radial velocity into the observations just before the PCA procedure. We then searched for these artificial signals and measured their significance order by order, by using the same procedure as for the molecule detection. The injected models were computed as those used for cross-correlation, except that they were amplified to be detectable at > 3 σ potentially in each order. Therefore an order was selected when the most significant signal was recovered within ±3 km s⁻¹ from V_rest = 0 km s⁻¹ and ±30 km s⁻¹ from the expected K_p, with significance greater than 3 σ. We underline that the injected signal is not used to optimise the number of PCA components via its recovery maximisation, but it is used only with the aim of selecting the orders where molecular absorption for a single species is likely to occur and/or the contamination by telluric residuals and other systematic effects (e.g. the wavelength-dependent efficiency of GIANO-B) is minimised.

Once we selected the spectral orders, we applied the PCA procedure. This procedure was applied to data of each spectral order selected for a particular molecular species, for each night. The aim of the PCA procedure is to find common spectral trends in the different wavelength channels (i.e. telluric and stellar lines) as a function of time and remove them. It consists of computing the principal components of the data matrix of each spectral order. In particular, each spectral order is an (M-rows; N-columns) data matrix (M = 58 − 62 spectra, N = 2048 pixels) and the principal components correspond to the eigenvectors of its covariance matrix. For each eigenvector, the associated eigenvalue denotes the contribution of that component to the data variance.

First, we normalised each spectrum to its median value, to correct baseline flux variations. Then, we masked the spectral channels with stronger or saturated telluric lines. To do this, we subtracted the time-averaged spectrum to all the spectra (i.e. we subtracted to each spectral channel its mean value), computed the standard deviation of data in each spectral channel and its median value (σ_m), and then we masked the spectral channels that had a standard deviation σ_λ > 1.5 · σ_m· After subtracting the time-averaged spectrum from all the spectra and masking highly contaminated spectral channels, we reduced data to a null-mean variable by subtracting to each spectrum its mean value computed on the different spectral channels. Then, the principal components were computed using the PCOMP IDL routine⁹. The output eigenvectors were ordered in decreasing contribution to the variance and we selected the minimum number of components that together contributed up to 70% to the data-variance. Once the number of principal components was selected, the matrix that should mainly describe the telluric and stellar contaminations was built via a linear combination of the principal components and removed from the original data matrix, obtaining an ideal residual telluric- and stellar-free matrix. Finally, we applied a high-pass filter to each row of the residual matrix to remove any possible residual correlation between different spectral channels. The goodness of the telluric (or stellar) lines removal was evaluated by a visual inspection of the residual matrix of each spectral order and of the cross-correlation values as a function of planetary orbital phase computed between data and the H₂O (or CO) model (see, e.g. Fig. 8).

In Table D.1 we report the selected spectral orders and the number of principal components removed by PCA for each selected spectral order, for the different molecular species and nights.

Appendix E Test on the reliability of the orders’ selection procedure

The spectral orders’ selection procedure described in the previous section is important to select the orders where molecular absorption for a single species is likely to occur considering different systematic effects (e.g. night-by-night variations of the observing conditions, telluric contamination, efficiency of wavelength calibration). To probe the possible presence of spurious positive signals near the expected planetary radial velocity that could lead to possible false positive molecular signals through the orders’ selection procedure we adopted, we performed the following test.

For each of the two detected molecules (i.e. H₂O and NH₃) and for each night, we started from the spectral orders’ lists selected with the procedure described in Appendix D (and reported in Table D.1) and refined the selection in the following way:

• step 1) for each of the selected spectral order, we applied the PCA and computed the signal-to-noise ratio (S/N) K_p – V_rest maps by cross-correlating the spectra of that order with the single-species model for the tested molecule;

• step 2) for each of the selected spectral order, we injected the single-species model for the tested molecule into the data at the planetary radial velocity and then applied the PCA and computed the signal-to-noise ratio (S/N) K_p – V_rest maps by cross-correlating them with the model itself;

• step 3) finally, for each spectral order, we computed the change (ΔS/N) in the S/N of the CCF peak at the planetary radial velocity obtained in the two previous steps (i.e. data alone and data+injected model) and selected the orders that showed a ΔS/N ≥ 2 σ.

With this kind of selection, we selected the orders that not only allowed to retrieve the injected signal but also showed a significant increase of the S/N of the CCF peak when the model itself was injected into the data, further minimising the possible presence of spurious systematics in the final dataset that could lead, for example, to false positive signals. In Fig. E.1 we report the ΔS/N for the different selected orders and nights. As it can be for most of the selected spectral orders the increase of S/N is > 2 σ.

Given the 2 σ threshold we discarded the following orders from the original lists:

• for H₂O: order 3 (night 1), order 15 (night 3);

• for NH₃: order 35 (night 1), order 21 (night 2), order 5 (night 4).

We then repeated the telluric removal procedure and the CCF analysis using the updated selected orders’ lists and built the K_p – V_rest maps (Fig. E.2), obtained by combining the signal of the four nights. As it can be seen, we confirm the detection of both the molecules obtained in the main analysis: for H₂O we obtain a CCF peak at the expected planetary radial velocity with a S/N= 5.1 and a t-test significance of 3.5 σ; for NH₃ we obtain a CCF peak at the expected planetary radial velocity with a S/N= 4.8 and a t-test significance of 4.0 σ. As it can be seen, removing the orders that showed a ΔS/N < 2 σ, the result of this test does not lead to significative differences (≤1 σ) with the conclusion reached in our work, giving more robustness to our analysis.

Appendix F Log of the HAT-P-11 b GIANO-B observations

In Table F.1 we report the log of the GIANO-B observations of the four transits of HAT-P-11 b.

Fig. E.1 Change in S/N (ΔS/N) due to the injection of the model of H2O (left) and NH3 (right) into the data, as a function of the selected spectral orders in the different nights. The black horizontal dashed line represents the 2 σ threshold.

Fig. E.2 Signal-to-noise ratio (left panels) and Welch t-test significance (right panels) Kp – Vrest maps for the 2 tested chemical species (top for H2O, bottom for NH3) obtained by combining data from the four observing nights excluding spectral orders with ΔS/N < 2 σ. The point where the 2 white dashed lines cross each other is the expected position of a signal with planetary origin.

References

All Tables

All Figures

	Fig. 1 Transit light curve analysis. Top panel: phase-folded Kepler transit light curve of HAT-P-11 b binned by 1 min, based on 64 bona fide transits (see text). Our best-fitting model of the transit in the Kepler bandpass is overplotted in grey. Bottom panel: residuals of the fit.
In the text

	Fig. 2 Occultation light curve analysis. Top panel: phase-folded Kepler occultation light curve of HAT-P-11 b binned by 1 min, based on 196 occultations (see text). Our best-fitting model of the occultation in the Kepler bandpass is overplotted with a black solid line. The data have been binned for clarity (black dots with error bars). Bottom panel: residuals of the fit.
In the text

Fig. 3 Geometric albedo (Ag) estimated as a function of the day-side temperature for the measured occultation depth . The cyan lines represent the 1σ uncertainty curves. The dash-dotted and dashed lines represent the variation of Ag with varying Td for the 2 heat re-circulation cases (ε) considered, computed at d = dsec = 14.83 ± 0.30 R⋆ (that is, during the occultation). The 2 orange shaded regions correspond to the intersection of the curves, identifying the average HAT-P-11 b day-side temperatures for the 2 extreme scenarios of ε, computed at d = dsec. These 2 values of the equilibrium temperature with the associated uncertainties are reported in the legend.

In the text

	Fig. 4 HIRES RV (top panel) and Call S-index (bottom panel) measurements of HAT-P-11. The two-time series show almost identical long-term variations with a shift of ~500 days in the minimum.
In the text

Fig. 5 HARPS-N RV (top panel), Call S-index (middle panel), and FWHM (bottom panel) measurements of HAT-P-11 during the four HAT-P-11 b transit nights for atmospheric characterisation. The three time series are highly correlated showing the same long-term trend, with no apparent shifts in the minima of variations. The vertical dashed line indicates the predicted periastron time of the hypothetical planet c, and the grey area shows its 1σ error bar accounting for the uncertainty on the orbital period from Yee et al. (2018). Note: The HARPS-N radial velocities were divided by their median of −63420.4 m s−1, to make a straightforward comparison with the HIRES radial velocities in Fig. 4.

In the text

Fig. 6 Sensitivity of the PMA technique to companions of given mass and orbital separation orbiting HAT-P-11. The black long-dashed and solid curves correspond to the combinations of mass and orbital radius explaining the PMA values at the mean Gaia DR2 and DR3 epochs, respectively. The shaded light blue region corresponds to the 1σ uncertainty domain of the DR3 PMA, while the shaded magenta region encompasses the 1 σ uncertainty of the DR2 PMA. The red diamond indicates the separation and minimum mass of the HAT-P-11 c companion proposed by Yee et al. (2018), while the red hexagon corresponds to the true mass value determined by the Xuan & Wyatt (2020) analysis.

In the text

Fig. 7 Radial-velocity data analysis. Left panels: HIRES radial velocities of HAT-P-11 (black circles) showing a clear long-term trend likely due to the stellar magnetic activity cycle. The best-fit models including both the Keplerian signal of HAT-P-11 b and Gaussian-process regression with squared-exponential (top; Eq. (3)), quasi-periodic (middle; Eq. (4)), and quasi-periodic+squared-exponential (bottom; Eq. (5)) kernels are indicated with red solid lines. Right panels: HAT-P-11 RVs (black circles) phase-folded with the transit ephemeris of HAT-P-11 b, after removing the long-term trends modelled with Gaussian processes and the same kernels as in the corresponding left panels. The Keplerian eccentric model is displayed with a red solid line. We note the smaller and smaller scatter in the residuals as we move from top (SE kernel) to bottom (QP+SE kernel).

In the text

Fig. 8 Examples of CCF values as a function of planetary orbital phase computed with data from the second observing night (18 June 2020) and model containing only H2O (top panel) or CO (bottom panel) lines. The horizontal dashed lines represent the transit ingress and egress while the dash-dotted line represents the expected CCF peak trail due to the planetary motion as measured in the observer rest frame. The expected CCF peak trail in transit is not represented for clarity. As it can be seen, due to the faintness of the signal, the CCF peak trail is not visible by eye. This kind of plots serves as a visual check of any remaining telluric and stellar residuals, in this case showing no residuals and signifying that these are adequately corrected by the PCA.

In the text

Fig. 9 Signal-to-noise ratio Kp – Vrest maps for the probed chemical species: H2O, CH4, NH3, C2H2, HCN, CO, CO2, and H2S. Each Kp – Vrest map shows the S/N of the cross-correlation of the GIANO-B spectra (4 transits combined) with isothermal atmospheric models, as a function of the planet’s RV semi-amplitude (Kp) and the planet’s rest-frame velocity (Vrest). The S/N is computed by dividing the peak value of the cross-correlation function at each Kp by the standard deviation of the noise far from the peak, as described in the text. Negative S/N values correspond to anti-correlation. The vertical and horizontal white dashed lines correspond to Vrest = 0 km s−1 and the expected Kp value (), respectively.

In the text

Fig. 10 Significance Kp – Vrest maps for the tested chemical species: H2O, CH4, NH3, C2H2, HCN, CO, CO2, and H2S. Each Kp – Vrest map shows the significance of the cross-correlation signal of the GIANO-B spectra (4 transits combined) with isothermal atmospheric models, as a function of the planet’s RV semi-amplitude (Kp) and the planet’s rest-frame velocity (Vrest). The significance is computed with a Welch t-test on two samples of cross-correlation values, that is, far from and near to the planet’s RV respectively (see the text for more details). The vertical and horizontal white dashed lines correspond to Vrest = 0 km s−1 and the expected Kp value (), respectively.

In the text

Fig. 11 Contour plots of the detection significance of the four selected chemical species. The solid (dashed) lines represent the 1 σ (2 σ) interval around the peak value of the significance (marked with a red cross). These intervals are computed as described in the text. The point of the Kp – Vrest map in which the 2 black dashed lines (the horizontal one corresponds to the expected Kp value, while the vertical one corresponds to Vrest = 0 km s−1) cross each other represents the expected detection significance peak position in the case in which the detected signal has a planetary origin. As it can be the significance peak position is compatible with the planetary origin hypothesis at better than 1 σ for all the 4 species.

In the text

Fig. 12 Most favourable atmospheric chemical models. Top-left panel: HAT-P-11 b atmospheric temperature-pressure profile corresponding to the most favourable chemical-equilibrium model that shows spectral features from all the detected species, including the tentatively detected ones (Model 1: [M/H]=2.0, C/O=0.9, N/O=0.85, βirr=0.5; solid line) and to a grid model modified including NH3 vertical quenching at the 10 bar level (Model 2: [M/H]=1.7, C/O=0.59, N/O=0.14, βirr=0.5; dashed line). Bottom-left panel: volume mixing ratios of the dominant chemical species (see legend) as a function of the atmospheric pressure corresponding to Model 1 (solid lines) and to Model 2 (dashed lines). Right panels: theoretical transmission spectra of HAT-P-11 b for atmospheres in thermochemical equilibrium (top-right panel) and with NH3 quenching (bottom-right panel). The coloured curves show the contribution to the synthetic spectra of the four species detected by GIANO-B, including the tentatively detected ones (see legend). The black curve shows a lower-resolution (R = 500) spectrum combining the absorption from all atmospheric species in the model.

In the text

Fig. 13 Distribution of Kp values fitted on randomly generated planetary RV curves with properties described in the text. The vertical black solid line represents the median value of the distribution, while the 2 vertical black dashed lines represent the 2 quantile values delimiting the 68% (1 σ) interval of Kp around the median. The vertical dash-dotted blue line represents the expected Kp value: = 123.4 km s−1. The median and the expected values are coincident and therefore they are not visibly distinguishable.

In the text

Fig. A.1 Example of the stages of the GIANO-B data reduction process applied to real GIANO-B spectra over a short wavelength interval (x-axis). The planetary orbital phase is reported on the y-axis. From top to bottom panel: (a) the extracted spectra; (b) residual spectra after having normalised each spectrum (each row) by its median value to correct baseline flux differences between the spectra that are due, for example, to variable transparency of the atmosphere, imperfect telescope pointing, or instability of the stellar point spread function; (c) residuals after each spectral channel (each column) had its mean subtracted; (d) residuals after having masked the strongest or saturated telluric lines and having removed the telluric and stellar spectra with the PCA; as it can be seen also other spurious effects such as the alternating pattern visible at certain wavelengths likely related to the instrument A-B nodding are effectively removed with PCA; (e) residuals after having applied a high-pass filter to each row of the residual matrix (this step mitigates any residual correlation between different spectral channels).

In the text

	Fig. B.1 Distribution of CCF values far (|Vrest| > 25 km s−1) from the planet’s rest-frame velocity (“out-of-trail”, OOT, dashed lines) and near to it (IVrest| < 3 km s−1, “in-trail”, INT, solid lines) for the probed chemical species.
In the text

Fig. C.1 Distribution of S/N obtained in the selected Kp – Vrest interval by shuffling 250 times the time order of the observed spectra, for each of the four selected chemical species. The vertical dashed (dash-dotted) lines represent the borders of the intervals enclosing the 95% (99.73%) of the S/N values, corresponding to the 2.5% − 97.5% (0.135% − 99.865%) quantiles of the distributions. These intervals are reported in Table 4.

In the text

Fig. C.2 Distribution of significance (from the t-test) obtained in the selected Kp – Vrest interval by shuffling 250 times the time order of the observed spectra, for each of the four selected chemical species. The vertical dashed (dash-dotted) lines represent the borders of the intervals enclosing the 95% (99.73%) of the significance values, corresponding to the 95% (99.73%) quantile of the distributions. These intervals are reported in Table 4. Since the 50% of the significance values (~20 000 values) are <0.2 σ for all the 4 chemical species, the plots of the distributions are limited on the y-axis in the interval [0; 650] for clarity.

In the text

	Fig. E.1 Change in S/N (ΔS/N) due to the injection of the model of H2O (left) and NH3 (right) into the data, as a function of the selected spectral orders in the different nights. The black horizontal dashed line represents the 2 σ threshold.
In the text

	Fig. E.2 Signal-to-noise ratio (left panels) and Welch t-test significance (right panels) Kp – Vrest maps for the 2 tested chemical species (top for H2O, bottom for NH3) obtained by combining data from the four observing nights excluding spectral orders with ΔS/N < 2 σ. The point where the 2 white dashed lines cross each other is the expected position of a signal with planetary origin.
In the text