© The Authors 2024

1 Introduction

Over the last few decades, the field of exoplanet research has grown rapidly, revealing that extrasolar systems are very common and extremely diverse in mass, radius, temperature, and orbital parameters. Thanks to increasingly efficient ground- and space-based surveys, today we are able to explore the exoplanet compositions and atmospheres in an ever-larger sample.

Puffy planets (i.e. planets with large radii and very low densities) constitute some of the most favourable targets for characterisation through the transmission spectroscopy technique (e.g. Sedaghati et al. 2016; Allart et al. 2020; Colón et al. 2020; Czesla et al. 2022). Due to their large radius and low surface gravity, they are expected to present a high atmospheric pressure scale height (H), which is the characteristic length scale of the atmosphere. The amplitude of an absorption signal in transmission spectroscopy, that is, the change in measured transit depth, is proportional to H ($2H{R}_{\text{p}}/R_{\star }^2$, Brown 2001). Transmission spectral signals are typically on the order of 1 H to ~5 H in size. Thus, if the transit depth can be measured to about 1 H in precision with sufficient spectral resolution, detectable spectral features would begin to appear. Since puffy planets are characterised by atmospheres with higher H values, they are expected to present stronger absorption signals (assuming the host star’s brightness and the signal-to-noise ratio are high enough). Some of the puffy planets that have been explored so far, such as KELT-11 b (Pepper et al. 2017), WASP−17 b (Anderson et al. 2010), and WASP−127 b (Lam et al. 2017), even show H ≳ 1900 km (compared to ~10 km on Earth or ~27 km on Jupiter). In these cases, the signal in transmission for 1 H is on the order of 10⁻²−1%, compared to ~2 × 10⁻⁵% on Earth and ~10⁻³% on Jupiter.

HAT-P-67 b (Zhou et al. 2017) is a puffy gas giant planet transiting a rapidly rotating (v sin i = 35.8 ± 1.1 km s⁻¹) F-sub-giant star every ~4.8 days at the orbital distance of 0.065 au. It is one of the largest (R_p ~ 2 R_J) and lowest density (ρ ~ 0.05 g cm⁻³) planets found to date. Due to the rapid rotation of the host star, the radial velocity (RV) technique does not allow for precise determination of the semi-amplitude of the planet K_p nor (consequently) its mass, for which there is only an upper limit (M_p < 0.59 M_J). A lower limit (M_p > 0.056 M_J) has been applied by Zhou et al. (2017) when assuming the planet is not undergoing Roche-lobe overflow. The host star belongs to a binary stellar system; however, the M dwarf companion (HAT-P-67B, Mugrauer 2019), has a projected separation of about 3400 au and is thus not a source of contamination for observations.

The low density and high irradiation of HAT-P-67 b (it receives approximately two times the incident flux of a zero-age main-sequence star) also results in a bloated atmosphere with a large H of ~3500 km (assuming an H₂−He mixture of near-solar composition atmosphere). This makes the planet another good candidate for transmission spectroscopy studies.

Recently, the atmosphere of HAT-P-67 b was explored by Bello-Arufe et al. (2023) through the analysis of one full transit acquired with CARMENES (Quirrenbach et al. 2016). The authors reported the detection of Call and Nai (with a statistical significance of 13.2 σ and 4.6 σ, respectively), which they ascribed to the planetary signal. Moreover, they found strong variability in the Hα line and in the He triplet, suggesting the possible presence of an extended planetary outflow. This extended atmosphere was confirmed by Gully-Santiago et al. (2023), who reported an absorption depth up to 10% in the stellar He I triplet, thanks to a series of observations taken with the Habitable Zone Planet Finder Spectrograph (HPF, Mahadevan et al. 2012). Gully-Santiago et al. (2023) also derived an increase of the stellar radius (2.65 ± 0.12 R_J) in order to match the updated Gaia DR3 distance (8.7% farther than previously estimated by Zhou et al. 2017). However, this update does not entail any major changes to the rest of the parameters derived by Zhou et al. (2017).

HAT-P-67 b is one of the targets of the atmospheric sample of the Global Architecture of Planetary Systems (GAPS; Guilluy et al. 2022). Thanks to the availability of the HARPS-N spectrograph (Cosentino et al. 2012) mounted at the Italian Telescopio Nazionale Galileo (TNG) in La Palma, we have collected and analysed four transits of HAT-P-67 b. The GIARPS observing mode (Claudi et al. 2018) allowed us to gather simultaneous observations both in the visible (VIS) with HARPS-N (0.39-0.68 µm, R ≃ 115 000) and in the near infrared (NIR) with GIANO-B (0.95−2.45 µm, R ≃ 50000), thereby providing us with the capability to investigate the system architecture and the exoplanet atmosphere in a wider wavelength range.

A description of the observations used in this work is presented in Sect. 2. Taking advantage of the multiple Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) light curves and the high-precision HARPS-N RV time series, we redetermined the orbital parameters and analysed the RML effect, shown in Sect. 3. We also tried to explore the activity level of the host star using different indicators in Sect. 4. We present the characterisation of the planetary atmosphere through analysis of the transmission spectrum both in the VIS and NIR in Sect. 5. Finally, we discuss our findings in Sect. 6.

2 Observations

For the analysis of the HAT-P-67 system, we used both high-resolution VIS and NIR spectroscopy with the GIARPS (GIANO−B + HARPS−N) observing mode of the TNG as well as photometry from TESS. A summary of the observations and stellar parameters adopted in this work for analysis is reported in Tables 1 and 2, respectively.

2.1 Spectroscopy

We observed four transits of HAT-P-67 b, three of which in the GIARPS mode, in the framework of the GAPS large program (Covino et al. 2013). The exposure time (t_exp) of the observations of HAT-P-67 b was fixed to 600 s for HARPS-N and 300 s for GIANO-B, yielding runs of 46, 31, 46, and 42 VIS spectra and 80, 50, and 68 NIR spectra during each transit. However, the very long duration of the transit (~7 h) resulted in quite a low number of out-of-transit spectra with respect to the total acquired spectra (only 23 out of 165 in the VIS and 29 of 198 in the NIR). To determine whether a spectrum was (fully) intransit or out-of-transit, we considered half of t_exp in addition to the Barycentric Julian date (BJD) of the observation. Sky spectra were retrieved simultaneously with the science observations thanks to a dedicated fiber, named fiber B, pointing at a fixed position at around 10 arcsec from the target star, ensuring the same atmospheric conditions in both spectra.

In the first run (N1), there are hints of stellar activity, which we discuss in Sect. 4. The second run (N2) covers only the second half of the transit and is characterised by a strong variability in seeing (1−3 arcsec), but it is the run with the highest number of out-of-transit spectra. During the third run (N3), GIANO-B was offline, so we collected only HARPS-N observations, which present some spectra with a lower signal-to-noise ratio (S/N), probably due to thin clouds and the presence of some calima. The fourth run (N4) does not cover the egress and the post-transit phase and is characterised by a lower S/N compared to the other nights (see Fig. A.1, right panel). We decided to remove the first exposure of N1 since the corresponding RV measurement deviates from the expected pattern, creating a distortion of the fit. For transmission spectrum analysis only, we also discarded the exposure of N3 with the lowest S/N (<20, at the centre of the 53rd order, with the sodium doublet).

2.2 Photometry

The time span covered by our spectroscopic survey almost completely overlaps with the TESS (Ricker et al. 2014) observations of HAT-P-67. This allowed us to determine the ephemeris of HAT-P-67 b close to the epochs of our GIARPS data and thus best suited to the extraction of the transmission signal of HAT-P-67 b. TESS observed the HAT-P-67 system with a 2-min cadence in sector 24 (S24; from 2020 April 16 to 2020 May 12, six transits observed), sector 26 (S26; from 2020 June 9 to 2020 July 4, six transits), sector 51 (S51; from 2022 April 23 to 2022 May 18, three transits), sector 52 (S52; from 2022 May 19 to 2022 June 12, five transits), and sector 53 (S53; from 2022 June 13 to 2022 July 8, four transits). Of the analysed TESS sectors, only S26 covers one of the simultaneously observed transits with GIARPS, which corresponds to the half transit retrieved during the second night.

Using the package lightkurve (Lightkurve Collaboration 2018), we retrieved the Pre-search Data Conditioning Single Aperture Photometry (PDCSAP), which was corrected for instrumental systematics and for contamination from some nearby stars (Smith et al. 2012; Stumpe et al. 2012, 2014). We took into consideration only the photometry with a good quality flag.

3 Orbital parameters

We redetermined the HAT-P-67 b transit parameters. Given the difficulty in the interpretation of the long-term stellar variability (Sect. 4), we trimmed segments of the LCs centred on the transit events and as wide as three times the expected transit duration. Each photometric segment was normalised to the median out-of-transit flux. A total of 22 transits were extracted from the TESS LCs.

We adopted the same Bayesian approach described in Scandariato et al. (2022): we maximised the likelihood of a model that includes the combined transit fit along with a GP to de-trend against long-term stellar or instrumental systematics. The presence of long-term trends is evidenced by the PSD of the data, which monotonically increases with decreasing frequency. We thus used a Matérn 3/2 kernel for the GP (see Fig. 1), as it is characterised by a PSD similar to the data. This approach is less time-consuming than fitting a polynomial trend for each transit. Moreover, it greatly reduces the dimensionality of the model. That is, there are two free parameters of the GP model compared to 22 × (deg + 1) parameters for the polynomial de-trending, where deg is the polynomial degree used for the de-trending.

The transit profile was computed using the quadratic limb darkening (LD) law of Mandel & Agol (2002) with the reparametrisation of the LD coefficients of Kipping (2013). In the model, we also included a re-normalisation factor and a jitter term to fit the white noise not included in the nominal photometric uncertainties. We sampled the parameter space in a Markov chain Monte Carlo (MCMC) framework using the Python emcee package version 3.1.3 (Foreman-Mackey et al. 2013). We used 44 walkers, corresponding to four times the number of free parameters. We ran the chains for 50 000 samplings, enough to ensure convergence following the criterion described in Goodman & Weare (2010). We used flat priors (listed in Table 3) for all the fitting parameters. We ran the code in the HOTCAT computing infrastructure (Bertocco et al. 2020; Taffoni et al. 2020). The results of the Markov Chain Monte Carlo (MCMC) fit are listed in Table 3. The best-fitting model is overplotted on the phase-folded data in Fig. 1.

We then used the RV time series to detect the RML effect (see, e.g., Albrecht et al. 2022) and measure the sky-projected angle between the stellar spin axis and the planet orbital axis. The data were reduced using the HARPS-N Data Reduction Software (DRS) version 3.7 through the Yabi web application, which is hosted at the Italian center for Astronomical Archive (IA2)¹. The RV measurements were obtained using a G2 mask template and a cross-correlation function (CCF) width of 250 km s⁻¹, with a step of 0.25 km s⁻¹. The lists of the RVs are presented in Tables A.1–A.4 together with the stellar activity index log $R_{\text{HK}}^{\prime }$ (see Sect. 4). We did not consider the observations taken on N2, as they do not cover a full transit and were affected by highly variable seeing (see Table 1 and Fig. A.1). The RML effect modelling and the RV fitting were carried out using a code we developed within the MATLAB software environment. A detailed explanation of the working principles of the code can be found in Esposito et al. (2017).

For the fit of the RML effect, most of the relevant parameters were adopted from the literature (M_★, R_★, v sin i; see Table 2) and from our analysis of the TESS LCs (T₀, P_orb, i, R_p; see Table 3). The only parameters that were left free to vary were the barycen-tric RV (υ_sys) and the projected spin-orbit angle (λ). Given the large uncertainty on the planet’s mass and considering that significant RV variations induced by stellar activity are expected on a timescale of a few hours, we also added an RV linear trend (LT) as a free parameter.

We separately fitted the three transit RV time series taken on 2020 May 26, 2021 June 10, and 2023 July 10. For the first (second, third) transit, the best-fit results are as follows: λ = 6.4 ± 7.5 deg (8 ± 11 deg, 10 ± 12 deg); υ_sys = −1968 ± 22 m s⁻¹ (−2225 ± 33 m s⁻¹, −2250 ± 45 m s⁻¹); and LT = 0.0057 ± 0.0015 m s⁻² (0.0021 ± 0.0020 m s⁻², 0.0010 ± 0.0020 m s⁻²). The measure of υ_sys found for N1 is statistically far from the values found for the other nights. This offset is likely due to a higher stellar variability (see Sect. 4). The three RV time series with the best-fit models superimposed are shown in Fig. 2.

In order to improve the precision on Λ, we performed a fit to the Doppler shadow in the same way as in Borsa et al. (2021b).

Indeed, for relatively fast rotators, such as HAT-P-67, the tomography method can better constrain Λ than RVs. The Doppler shadow model was taken from EXOFASTv2 (Eastman 2017) and fitted to the data in a Bayesian framework by employing a differential evolution MCMC (DE-MCMC) technique (Ter Braak 2006; Eastman et al. 2013), running ten DE-MCMC chains of 50 000 steps and discarding the burn-in until convergence was reached. We fixed T₀ and P_orb as in Table 3 as well as the quadratic limb darkening parameters (µ₁ = 0.43 and µ₂ = 0.25, taken from ExoCTK²) and υ_sys = −2000 ms⁻¹. We note that the vsys value cannot be well constrained by the Doppler tomography fit for fast rotators, which is contrary to what happens for RVs, and that changing this value within the differences found by the RVs analysis of the different transits does not affect the results. We left i, a/R_★, Rp/R_★ as free parameters with values and error bars in Table 3 as priors. The υ sin i (which includes macroturbulence) and λ parameters were left free with uniform priors. The medians and the 15.86% and 84.14% quantiles of the posterior distributions were taken as the best values and 1 σ uncertainties. We independently fitted all four transits, finding υ sin i = 39.2 ± 0.5, 36.7 ± 0.9, 38.3 ± 0.7, 39.5 ± 0.8 km s⁻¹ and λ =1.0 ± 0.6, 5.5 ± 1.1, 1.2 ± 0.8, 3.2 ± 0.9 degrees for night one to four, respectively. When fitting the four transits together, we found υ sin i = 39.3 ± 0.4 km s⁻¹ and λ = 2.2 ± 0.4 degrees. As expected, the best values of λ are both more accurate and precise than the ones found by fitting the RVs. We note that the v sin i value is compatible with that found by Zhou et al. (2017) within 3σ.

4 Stellar activity

Stellar activity can mimic spurious features in the retrieved transmission spectra (e.g. Oshagh et al. 2014; Apai et al. 2018). To detect potential stellar effects in the analysed nights, we first measured the log $R_{\text{HK}}^{\prime }$ chromospheric activity index (Fig. 3 and Tables A.1–A.4). We extracted it from the HARPS-N spectra through the YABI platform (Hunter et al. 2012) using a B-V colour index of 0.441 mag. We derived average values of −4.651 ± 0.001, −4.675 ± 0.003, −4.743 ± 0.004, and −4.704± 0.007 for N1, N2, N3, and N4, respectively. The obtained values are moderately larger than the solar value at the maximum of activity (−4.75 and −4.905 according to Dumusque et al. 2011; Egeland et al. 2017, respectively). Compared with the other nights, N1 presents slightly higher activity, while N4 is characterised by greater uncertainty, likely due to the highest S/N and airmass (Fig. A.1). However, it is important to note that the chromospheric log $R_{\text{HK}}^{\prime }$ index is thought to be significantly depressed in stars with transiting giant planets. This is because of the absorption by a circumstellar torus produced by planetary evaporation, which is particularly strong in planets with a very low surface gravity, such as HAT-P-67 b. Considering the data and the linear models relating the chromospheric index with the inverse of the surface gravity (Lanza 2014; Fossati et al. 2015), the true value of log $R_{\text{HK}}^{\prime }$ of HAT-P-67 could be larger by at least 0.4−0.5, that is, it can approach −4.2.

The occurrence of phenomena associated with stellar activity during N1 is also reflected in the corresponding CCF linear profile (Fig. 4, left panels), which is much more distorted and time-varying than N3 and N4 (N2 is not considered here since it only covers half a transit). The distortion of the stellar CCF also strongly affects the measured RVs. Indeed, the systemic velocity obtained fitting the RVs on N1 presents an offset of ~200 m s⁻¹ compared to N3 and N4 (see Sect. 3).

The origin of this distortion is most likely some kind of stellar variability. To support this hypothesis, for each night we Doppler-shifted the individual CCFs to the stellar rest frame and computed the average CCF. Then we divided each CCF by the average one and built the two-dimensional stack of the residuals, shown in the right panels of Fig. 4. In addition to the evident Doppler shadow due to the planetary transit (the tilted trace), some vertical patterns can be seen in N1, even outside the transit. Being at rest in the stellar rest frame, these patterns cannot be ascribed to planetary effects. We thus argue that they arise from the evolution of the stellar CCF profile during the transit. By applying the SpotCCF tool (Di Maio et al. 2023), we could model the deformation of the CCF profile. Assuming differential stellar rotation, we found compatibility with a dark spot signature only in the hypothesis that the stellar rotation axis is inclined about 44 degrees from the line of sight. Despite the star having a relatively low temperature (T_eff ~ 6406 K), another possible explanation could be the presence of non-radial stellar pulsations (e.g. Rieutord et al. 2023). Unfortunately, the time series we had available did not allow us to investigate this further.

We searched for evidence of stellar activity in the TESS LCs. To analyse the photometric variability, we clipped out the in-transit photometry and computed the generalised Lomb–Scargle periodogram (Zechmeister & Kürster 2009; Ferraz-Mello 1981) to detect any periodic signal. We found a clear periodicity at ~5.4 days in S24 and S26 (year 2020) with an amplitude on the order of 1 mmag. In S51, S52, and S53 (year 2022) we did not find any clear indication of periodic signals (Fig. 5). As a sanity check, we performed the same analysis using the TESS Simple Aperture Photometry (SAP) and obtained the same results.

Using R_* indicated by Gully-Santiago et al. (2023) (2.65 ± 0.12 R_⊙) and in the scenario of an equator-on star, our estimate of the stellar rotation period led to an equatorial rotation velocity of 25 ± 7 km s⁻¹. Due to its large uncertainty, this value is consistent within 1 σ of the υ sin i reported by Zhou et al. (2017) (30.9 ± 2 km s⁻¹) and within 1 σ from our derived υ sin i. This supports the hypothesis that the periodicity of ~5.4 days found in the periodogram (Fig. 5) is close to the real stellar rotation period. If, conversely, we assume an inclination of ~44 degrees for the stellar rotation axis, we then derive a rotation period of ~3.7 days, which does not correspond to any clear peak in the periodogram of the TESS photometry. We thus postulate that the star is most likely seen in an equator-on configuration.

Fig. 1 Results of the fit of the transit light curves (LCs). Top: phase-folding of the de-trended data (top panel) and the corresponding O–C diagram (bottom panel). For clarity, in each panel, we show the binned data with black dots. Bottom: power spectral analysis of the TESS photometry. The grey dashed line shows the power spectral density (PSD) of the data after removing the best-fit transit model, while the blue solid line represents the combined PSD of the best-fit Gaussian process (GP) (the shoulder at f ≲ 1 day−1) and the white noise in the data (the plateau at f ≳ 10 day−1). The black line is the PSD of the residuals of the best fit, showing that the power excess at low frequencies has been effectively removed by the GP in the model.

Fig. 2 Rossiter-McLaughlin effect analysis. Top: RV time series taken during the transit on 2020 May 26 (N1). The best-fit model is superimposed, and the corresponding residuals are shown in the lower panel. The dashed black lines indicate the points of the first and fourth contact of the transit, while the dotted black lines represent the points of the second and third contact. Middle: same but for the transit on 2021 June 10 (N3). Bottom: same but for the transit on 2023 July 10 (N4).

Fig. 3 Activity indicator log $R_{\text{HK}}^{\prime }$ as a function of the time for each night. The continuous black lines indicate the points of first and fourth contact of the transit, while the dashed black lines represent the points of second and third contact.

Fig. 4 Cross-correlation function for N1 (first row), N3 (second row), and N4 (third row) extracted from the DRS using a G2 mask. Left: cross-correlation function linear profiles as a function of time. The black line indicates the average CCF. Right: two-dimensional maps of the CCFs as a function of orbital phase and RV in the stellar rest frame. The four horizontal black lines show the times of transit contacts. The straight white line shows the expected planetary Keplerian velocity Kp using a planetary mass value of ~0.34 MJ (Zhou et al. 2017).

Fig. 5 Generalised Lomb-Scargle periodogram of the five TESS sectors and the joint LC (bottom box). The vertical dashed line marks the orbital period of HAT-P-67 b reported in Table 3.

5 Atmospheric characterisation

5.1 Extraction of the transmission spectra in the visible range

The standard data reduction in the VIS range is performed using the HARPS-North dedicated DRS, which produces both 2D and 1D spectra. For each night of observation, we analysed the 2D spectra using the SLOPpy³ (Spectral Lines Of Planets with python) pipeline (Sicilia et al. 2022). The SLOPpy pipeline is a standard user-friendly and reliable tool that is optimised for spectral reduction and the extraction of planetary transmission spectra in the VIS obtained from high-resolution observations. First, SLOPpy applies several data reduction steps that are required to correct the input spectra for sky emission, atmospheric dispersion, and the presence of telluric features. To that end, the pipeline subtracts the sky spectrum, which is simultaneously retrieved with the science observations, from the stellar spectrum (see Sect. 2.1). After dividing each observation with a reference spectrum, the pipeline models this ratio with either a low-order polynomial or a spline and finally divides each observation by this model. Among the different approaches implemented in the pipeline, we decided to apply the one that uses the atmospheric transmission code Molecfit (Smette et al. 2015; Kausch et al. 2015).

After applying the aforementioned data reduction steps, each in-transit observation was divided by a master-out spectrum (M_out, that is, the integration of the exposures out-of-transit acquired before the ingress and after the planet’s egress). In this way, the pipeline removed the stellar contribution, and in principle, the residuals should contain the exoplanet atmospheric signal. The M_0ut was built by moving and combining all the out-of-transit spectra to the stellar rest frame. The wavelength shift depends on the barycentric Earth RV (BERV) and the υ_sγs of the star. While the BERV values are provided by the DRS in the header of the FITS files, we derived υ_sγs (−2.234 ± 0.027 km s⁻¹) by taking the weighted average of the two values found from the fit of the RVs of N3 and N4 (we excluded the value found for N1, as it is likely to be contaminated by higher stellar variability; see Sect. 4).

We did not consider the reflex motion of the star since, being a fast rotator, even a wavelength shift on the order of one pixel (corresponding to an RV shift of ~0.8 km s⁻¹ in the case of HARPS-N, i.e. the RV variation due to a planet of ~8 M_J) does not change the shape of the spectrum noticeably due to the large broadening of the spectral lines. Even assuming the planet mass upper limit (M_p = 0.59 M_J, Zhou et al. 2017), the maximum stellar RV change is ~±0.05 km s⁻¹. Thus, the out-of-transit spectra can be co-added without taking into account the stellar reflex motion due to the planet (K_★ = 0 km s⁻¹).

Figure 6 shows M_out in the region of the sodium doublet, obtained by combining all four nights. The presence of interstellar lines is evident, as expected given its large distance from Earth (~320 pc). Any small variation of these lines can mimic a false signature. However, not correcting by the stellar reflex motion and assuming that interstellar lines remain at the same spectral position and are totally stable during the night, they would be automatically removed when dividing each spectrum by M_out in the stellar rest frame (e.g. Casasayas-Barris et al. 2018).

The two main effects altering the transmission spectra, namely the center-to-limb variation (CLV) and the RML effect, were also taken into account (simultaneously). Based on Spectroscopy Made Easy (SME, Piskunov & Valenti 2017) and using a line list from the VALD database (Ryabchikova et al. 2015) and Kurucz ATLAS9 (Kurucz 2005) models, we obtained synthetic stellar spectra at different limb angles (ranging from zero to one). The SLOPpy pipeline divides the observed transmission spectrum by a synthetic transmission spectrum computed including the RML and CLV effects in the stellar models but without the planetary absorption. This step requires the knowledge of some parameters such as the projected spin-orbit angle (λ), the differential rotation rate (a), and the stellar inclination (i_*). For the first parameter, for which Zhou et al. (2017) measured an upper limit of 12°, we considered the value found through the tomography method (2.2 ± 0.4°). For most of the stars, α has not been measured reliably, and in this case, the default choice was to exclude the differential rotation in the model and to assume a rigid-body rotation of the star (α = 0). Regarding i_★, on the other hand, considering what we stated in Sect. 4, we expected the star to be close to equator-on (i_★ ≃ 90°). Nevertheless, when rigid-body rotation is assumed, the i_★ value is not relevant when modelling the CLV and RML effects simultaneously.

Fig. 6 Composite master-out spectrum normalised to unity around the Na I doublet lines. Deep and narrow interstellar features peak out from the wider stellar lines.

5.2 Cross-correlation analysis

We investigated the atmospheric composition of HAT-P-67 b using the cross-correlation technique. We used high-resolution model templates at the temperature of 1903 Κ (Zhou et al. 2017) for the following atoms and molecules: Cal, Cri, Fel, Fe II, H₂O, Κ I, Mg I, Na I, Ti I, TiO (using line lists from Exo-mol and Plez), V I, VO, and Y I. The synthetic models of the studied species were generated using petitRADTRANS (Mollière et al. 2019) assuming equilibrium chemistry, solar abundance from Asplund et al. (2009), and stellar and planetary parameters from Tables 2 and 3. Additionally, in order to simulate the continuum opacity produced by H⁻, we added a cloud layer at P₀ = 1 mbar. In the final step, the high-resolution synthetic spectra were convolved to the HARPS-N resolution using instrumental (Gaussian kernel) broadening instrBroadGaussFast from PyAstronomy⁴ (Czesla et al. 2019).

We cross-correlated the planetary transmission spectra extracted in the VIS range with all the templates listed above and then combined the CCFs to obtain the K_p−υ_sγs maps. We did not detect any robust signal except for Cr I, Na I, Fe I, and Ti I.

We remark that these species are commonly found in the transmission spectra of Hot Jupiters (HJs) (see e.g. Hoeijmakers et al. 2018; Ishizuka et al. 2021; Scandariato et al. 2023). The stacked CCFs showed, night by night, an absorption feature stretching between the first and fourth contacts (t₁ and t₄) and shifting in the velocity space according to the expected planetary Keple-rian motion (Figs. A.2–A.5). The corresponding K_p−υ_sys maps confirmed this evidence by showing a clear absorption feature near the expected Kp and vsys of the system (Fig. 7). Assuming that the median absolute deviation of the K_p−υ_sys maps is a good estimate of the noise, the S/N of the absorption features is then above 6, 21, 9, and 11 for Cr I, Fe I, Na I, and Ti I respectively, when combining all nights.

Since the spectral lines of these atomic species are also present in the stellar spectrum, it is unclear whether these detections are due to either stellar residuals in the transmission spectra or planetary absorption lines. To test the origin of the aforementioned detection, we checked whether the stellar spectrum has been effectively removed by the reduction pipeline. To do this, we computed a cross-correlation mask for HAT-P-67 using T_eff = 6500 K, log g = 4.0, [Fe/H] = 0, and the theoretical line list provided by the VALD database. Then, we removed all the CrI, Fe I, Na I, and Ti I spectral lines from the mask in order to leave only the species that do not lead to any detection. Finally, we cross-correlated the transmission spectra, where in the ideal case there should be no signature of the stellar spectrum, with the modified mask. We found that there is a residual trace in the stacked CCFs (Fig. A.6) along the expected position of the planetary absorption that closely follows the stellar Doppler shadow (see Fig. 4). This trace, translated to the K_p−υ_sys plane, generates an absorption feature close to the expected planetary Kp. The most obvious explanation for this result is that the stellar spectrum has not been completely removed from the in-transit spectra. To mitigate this problem, we tried to redo the extraction by varying the relevant parameters (K_p, υ sin i, λ) within their uncertainties. Unfortunately, the leftover of the stellar spectrum never faded out.

We thus obtained clear evidence that the removal of the stellar spectrum during the extraction of the transmission spectra is unable to completely remove the stellar signature, casting serious doubts on the reliability of the planetary Cr I, Fe I, Na I, and Ti I detections seen in the transmission spectrum. This does not exclude the fact that a genuine planetary signal exists. Unfortunately, however, the fact that the planetary trace is expected to follow the Doppler shadow complicates the unequivocal attribution of the signal.

Fig. 7 Planes of Kp−υsys obtained by combining all nights for the Cr I (top left), Fe I (top right), Na I (bottom left), and Ti I (bottom right) templates. In each plot, the dashed cross marks the expected Kp and vsys of the system.

Fig. 8 Transmission spectra of HAT-P-67 b for each night. The spectra are centred around the Hα line in the planetary rest frame (light grey) and also binned 20 times (in black circles). The red line is the MCMC Gaussian fit performed by SLOPpy, while the vertical blue dashed line indicates the rest frame transition wavelength of the Hα line.

5.3 Transmission spectroscopy of the H-alpha line with HARPS-N

As already mentioned, the atmosphere of HAT-P-67 b has recently been explored by Bello-Arufe et al. (2023) through the analysis of one full transit retrieved with CARMENES. In addition to the detection of Na I, which they ascribe to the planetary signal, they also reported a strong absorption near the H-alpha (Hα), the first spectral line in the Balmer series. We searched for the same signal in our optical transmission spectra $(\tilde{\Re })$ extracted from SLOPpy (see Sect. 5.1) by summing all in-transit observations divided by the M_out in the planet’s reference system. The result for each night is shown in Fig. 8. The analysis revealed the presence of a strong absorption feature in N1, N2, and N3, while a clear emission feature is visible in N4.

For a measure of the absorption signals, we decided to apply an MCMC Gaussian fit in the region of the Hα line. A summary of the best-fit parameters obtained for each night is reported in Table 4. By combining all nights, we found an absorption signal with a contrast (c) of ~2.2% and a full-width half maximum (FWHM) of ~38 km s⁻¹. However, the signal was characterised by a very high redshift (~23 km s⁻¹) with respect to the predicted line position. We point out that we get compatible results even if we set a prior on Kp next to the expected value.

Doppler shifts as strong as ~23 km s⁻¹ have never been claimed in the literature. Moreover, global circulation models of the atmospheres of HJs predict the presence of zonal winds flowing from the planetary dayside to the nightside (e.g. Komacek & Showman 2016; Parmentier & Crossfield 2018; Roman et al. 2021). This would lead to a blueshifted atmospheric signal, which is in contrast with our finding. We thus postulate that the planetary origin of the Hα absorption signal is unlikely. Nonetheless, a more thorough investigation, though beyond the scope of this paper, is needed to give an appropriate interpretation of the absorption feature seen during the transits.

Regarding N4, as explained in Appendix B, we found that the emission feature originates from around mid-transit and seems to be at rest in the stellar reference system. As shown in Fig. 3, N4 presents higher variability in the chromospheric activity index log $R_{\text{HK}}^{\prime }$ right at the centre of the transit until the end of the observations. This leads us to state that the observed emission signal is probably due to some stellar activity event.

5.4 Transmission spectroscopy of the near infrared He I triplet with GIANO-B

Similar to the optical analysis, in order to separate the potential planetary He I signal from the stellar contribution, we performed transmission spectroscopy on the GIANO-B spectra using the approach outlined in Guilluy et al. (2023, 2024).

We employed the GOFIO pipeline (Rainer et al. 2018) to extract the spectra from the raw GIANO-B images and a preliminary wavelength calibration using a U–Ne lamp spectrum as a template in the vacuum wavelength frame. We then refined this initial wavelength solution by employing the same two-step approach described in our previous works (e.g. Guilluy et al. 2020; Giacobbe et al. 2021), which consists of aligning all the spectra to the Earth’s atmospheric rest frame, assuming it to be the frame of the observer (disregarding any ~10 m s⁻¹ differences due to winds). We thus focussed on order #39, which includes the He I triplet.

We employed Molecfit to remove the telluric absorption contamination. Additionally, we masked the emission telluric line at around 1083.43 nm following the methodology described in Guilluy et al. (2023, 2024).

We moved the telluric-corrected spectra to the stellar rest frame using parameters listed in Tables 3 and 2. We normalised each spectrum to the continuum by median division, neglecting the spectral region near the He I triplet and excluding spectra with significantly lower signal-to-noise ratios compared to the other exposures. Similar to Gully-Santiago et al. (2023), we then created a master-before, master-in, master-after spectra by averaging the before-transit (i.e. with an orbital phase smaller than t₁), after-transit (greater than t₄), and in-transit (between t₂ and t₃) frames. Through visual comparison of these master spectra, an absorption feature was readily discernible both in the in-transit and before-transit spectra at the position of the HeI triplet. This is particularly evident when the nights are considered together (last two panels of Fig. 9). This conspicuous excess absorption (~8%) is visually compatible with what is found in Gully-Santiago et al. (2023; see their Fig. 8).

Since not all of the observed nights provided complete transit coverage, as seen in Fig. A.1, we decided to consider all transits collectively. This is particularly crucial for N1, where the only available comparison stellar spectrum was obtained before the transit. However, considering that Gully-Santiago et al. (2023) identified a helium tail preceding the planet by several hours, we were unable to employ these spectra for comparison due to the potential influence of the planetary tail. In the following parts of this section, we present the results we obtained by considering only the two nights with the highest S/N, namely N1 and N3. For completeness, the findings obtained for all the investigated nights are illustrated in the Appendix (see Fig. A.7 and Table A.5).

We considered as a comparison spectrum, an average spectrum of all the observations acquired after the transit S_aft(λ). We thus derived the individual transmission spectra, T(λ, i), by dividing each spectrum by S_aft(λ). Finally, we linearly interpolated the transmission spectra in the planet’s rest frame. The upper panels of Fig. 10 display the 2D transmission spectroscopy map in the star (panel a) and planet (panel b) rest frame, where an absorption signal is visible. With the goal of monitoring the variation of the He I signal during transit, we additionally performed spectrophotometry of the helium triplet within a passband of 0.075 nm centred at the peak of excess absorption in the planet rest frame (Allart et al. 2019). The computed transit light curve is presented in the panel c of Fig. 10. Finally, we obtained the fully in-transit transmission spectrum in the planet’s rest frame T_mean(λ) by averaging the transmission spectra with an orbital phase between t₂ and t₃.

Following Guilluy et al. (2024), we estimated the contrast c of the excess absorption at the position of the He I triplet by fitting a Gaussian profile to T_mean(λ) with the DE-MCMC method, varying the peak position, the FWHM, the peak value (c) and the offset for the continuum. We also accounted for the presence of correlated noise by employing GP regression within the same DE-MCMC tools, using a covariance matrix described by a squared exponential kernel (Bonomo et al. 2023). Finally, we considered possibly uncorrelated noise by introducing a jitter term σ_j. Panel d of Fig. 10 shows the final T_mean corrected with the GP model and with the best-fit Gaussian model overplotted. We detected a He I contrast c of ${5.56}_{-0.30}^{+0.29}\%$ at 19 σ with a velocity shift of $\upsilon ={10.0}_{-1.2}^{+1.1}\hspace{0.17em}\text{km}\hspace{0.17em}{\text{s}}^{-1}$. The best-fit parameters from the DE-MCMC Gaussian analysis are listed in Table 5, while the posterior distribution and the application of the GP model are reported in the Appendix in Figs. A.8 and A.9, respectively.

We then translated c into an effective planetary radius R_eff ~ 3 R_p (e.g. Chen et al. 2018), and the extension beyond the planet’s Roche lobe radius (2.7 R_p, Gully-Santiago et al. 2023) indicates that the planet is evaporating. We found a restricted Jeans escape parameter Λ~18 (Fossati et al. 2017). The value obtained supports the presence of an extended atmosphere with significant mass loss, as one would expect from seeing metastable helium absorption beyond the Roche lobe. Finally, we derived the quantity ${\delta }_{{\text{R}}_{\text{p}}}/H_{\text{eq}}$ (Nortmann et al. 2018), which represents the number of scale heights probed by the atmosphere in the spectral range under consideration (see Table 5).

We advise to take caution with our He I detection and emphasise that our measured signal could be contaminated by potential pseudo-signals of stellar activity, caused by the planet’s passage during transit over a non-uniform stellar disc (e.g. Salz et al. 2018; Guilluy et al. 2020, 2024). The variability observed in the Hα line (Sect. 5.3), likely attributed to stellar activity, is a warning signal. However, in the He I analysis, we considered all the nights together, primarily due to the absence of post-transit observations in both N1 and N4, which are crucial if a pre-transit tail exists. For this reason, we were unable to analyse possible He I night-to-night variability due to stellar activity.

Fig. 9 Master-after, master-before, and master-in spectra in the star frame. Panels a–c represent the GIANO-B nights considered individually. Panel d shows all three transits together, while in panel e only the two observations with the highest S/N are considered. Vertical lines indicate the position of the He I. In- and before-transit absorption is visible by eye.

Fig. 10 Transmission spectroscopy near the He I triplet. Panels a and b show the 2D transmission spectroscopy maps in the stellar and planet rest frames, respectively. Tilted and vertical lines denote the position of the He I triplet, while horizontal dotted lines mark the transit contact points, t1, t2, t3, and t4. Panel c shows the spectroscopic light curve computed in a 0.075 nm band centred around the peak of absorption. Vertical dashed lines mark the transit contact points. Panel d presents the full in-transit averaged transmission spectrum.

6 Summary and conclusions

HAT-P-67 b is one of the lowest-density exoplanets known to date (ρ ~ 0.05 g cm⁻³), and it is a prime target for transmission spectroscopy. In the framework of the GAPS programme, we observed four transits of HAT-P-67 b using the GIARPS mode of the TNG, aiming at the simultaneous VIS and NIR study of the exoplanet atmosphere.

We first derived a new orbital solution for the HAT-P-67 system using archival TESS photometry and found it to be compatible within 2 σ with Gully-Santiago et al. (2023), who analysed the same TESS sectors. We combined our orbital solution with the RV times series to analyse the RML effect (Rossiter 1924; McLaughlin 1924). From the fit, we estimated the systemic velocity of the system and the projected spin-orbit angle. In order to improve the precision on λ, we applied the Doppler tomographic analysis, finding a value of 2.2±0.4° when combining all four transits. Our result is compatible with Zhou et al. (2017), who measured an upper limit of 12° through the same technique. The derived value indicates an aligned planetary orbit and suggests that the planet has likely migrated to its current orbit via tidal interactions with a protoplanetary gas disc (Lin et al. 1996; Ward 1997). It is unlikely that the outer stellar companion is directly involved in the migration of the planet through, for example, the Kozai mechanism because of its very large separation (3400 au). Nevertheless, some combinations of the original planetary orbit (e.g. ≥1 au) and binary orbit may allow for a phase of high eccentricity and very close periastron followed by tidal circularisation on the current orbit within the age of the system. The tidal timescale for the damping of the obliquity, computed for a planetary mass of 0.32 M_J (Gully-Santiago et al. 2023), is about 0.6 Gyr, while for the mass upper limit of 0.59 M_J, it is about 0.3 Gyr. We adopted a strong tidal interaction assuming a modified tidal quality factor of the star $\left(Q_{\star }^{\prime }\right)$ of 2.5 × 10⁴, as in the case of Kepler-1658, which is another system with a similar sub-giant star (Vissapragada et al. 2022). Therefore, an alternative interpretation could be that the planet mass is close to 0.59 M_J, and the low obliquity has been reached only relatively recently, when the star was in the final phase of its main-sequence evolution and experienced a strong tidal interaction with HAT-P-67 b. This is also in agreement with a stellar rotation period close to the orbital period because tides would act both to damp the obliquity and to synchronise the rotation of the star.

Our first night of observation shows signatures of higher stellar variability than on the other nights analysed. As a matter of fact, compared with the remaining part of our survey, we found an RV offset of ~200 m s⁻¹, a higher log $R_{\text{HK}}^{\prime }$ activity index, and a distorted CCF profile. Furthermore, although the TESS photometry does not cover our first night of observation, we nevertheless extracted a clear periodicity consistent with stellar rotation from the available TESS LCs, spanning a period fairly close to the one observed. By modelling the CCF profiles, we found that the observed stellar variability does not appear to be due to stellar spots, which would have generated a different periodicity than what was instead extracted from TESS photometry. It is more likely that the source of this variability is related to the presence of non-radial pulsations on the stellar surface.

Using the cross-correlation technique, we searched for the presence of atomic or molecular species in the optical transmission spectrum of HAT-P-67 b. Despite the expected high absorption signal, we did not get any robust detection with the exception of the Cr I, Fe I, Na I, and Ti I spectral lines, finding a formal statistical significance of ~6 σ, 21 σ, 9 σ, and 11 σ, respectively, in the combined signals. However, when we tested the robustness of the detections, we could not confirm their planetary origin. By carrying out a series of tests, we realised that we were not able to completely remove the stellar contribution from the in-transit spectra, leading to an unreliable transmission signal. After checking that the same analysis technique works well on other targets with clear detections (e.g. KELT-20 b), we assumed that the models used for the data reduction, especially the removal of the contribution of the RML and CLV effects, poorly match the real spectrum of HAT-P-67. The use of 3D models could improve the characterisation of the exoplan-etary atmosphere. Evidence for this was recently reported by Canocchi et al. (2023), who applied 3D non-LTE synthetic spectra to estimate the stellar RM+CLV effects in transmission spectra of solar-like planet hosts, showing that 1D models seem to overestimate the CLV signature. Even earlier, Chiavassa & Brogi (2019) showed that removing the stellar spectrum using 3D radiative hydrodynamical simulations leads to a significant improvement in planet detectability, both in solar-type and K-dwarf stars.

Due to the strong irradiation from the host star combined with internal heating within the planet and its very low escape velocity υ_ese ~ 25 km s⁻¹, Zhou et al. 2017), HAT-P-67 b’s atmosphere has been undergoing an intense radius inflation. This should have led to the formation of an extended hydrogen atmosphere. Gas giants with masses below Jupiter and temperatures above 1800 K, such as HAT-P-67 b (T_eq ~ 1900 K), are so inflated and puffed up that they are all on unstable evolutionary paths that will eventually lead to Roche-lobe overflow and the evaporation and loss of the planet’s atmosphere (Batygin et al. 2011). In the VIS wavelength range, the Hα line is a powerful probe for the escaping atmosphere (Yan & Henning 2018; Borsa et al. 2021a; Czesla et al. 2022). On HAT-P-67 b, our results for the analysis of the first three nights are compatible with what was found by Bello-Arufe et al. (2023), that is, a highly red-shifted (~23 km s⁻¹) strong absorption signal (~2.5%). However, this is an unexpected finding from the current global circulation models of HJs. We thus argue that such an absorption signal can hardly be ascribed to a planetary origin. Only for the last night of observations did we find a clear emission signal in the transmission spectrum. The correlation with the chromospheric index log $R_{\text{HK}}^{\prime }$ leads us to believe that the origin of the emission feature is most likely stellar activity.

Gully-Santiago et al. (2023) found no detectable variability in other hydrogen lines but confirmed the presence of a large and variable He I tail preceding the planet. Thanks to the GIARPS observing mode of the TNG, we were able to extract the transmission spectrum in the region of the NIR He I triplet. Using the average of all the spectra acquired after the transit as the master-out spectrum, we found a clear absorption signal. Applying a Gaussian fit to the full in-transit mean transmission spectrum, we estimated a contrast of the excess absorption of ${5.56}_{-0.30}^{+0.29}\%$ (19.0 σ). This value corresponds to an effective planetary radius of ~3 R_p (~6 R_J), indicating that the planet’s atmosphere is evaporating. Assuming a scale height (H) of ~3500 km and considering 1 σ uncertainties, we found that the He I atmosphere probes a number of scale heights varying between ~13 and ~84 H.

In conclusion, despite HAT-P-67 b being one of the most favourable targets for transmission spectroscopy studies, there are factors, such as stellar variability and observational constraints (i.e. the long transit duration), that make analysis challenging and therefore less accurate. Furthermore, the fact that the orbital trace of the planet is almost identical to that of the Doppler shadow complicates the unequivocal attribution of the found signal to either a planet or a star.

Acknowledgements

The authors acknowledge financial contribution from PRIN INAF 2019 and from the European Union–Next Generation EU RRF M4C2 1.1 PRIN MUR 2022 project 2022CERJ49 (ESPLORA). We thank the anonymous referee for their thoughtful comments which helped to improve the quality of this work. This paper includes data collected by the TESS mission, which are publicly available from the Mikulsky Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by the NASA’s Science Mission directorate. D.S. thanks L. Pino for discussions and insights. L.M. acknowledges support from the MIUR-PRIN project no. 2022J4H55R.

Appendix A Additional figures and tables

Fig. A.1 Variation of airmass (left) and S/N (right) for each night of HAT-P-67 b. The S/N was extracted from the FITS header of the HARPS-N spectra on the 53rd order, which contains the sodium feature. The continuous and the dashed vertical black lines represent the four points of contact of the transit.

Fig. A.2 Cross-correlation signatures of Cr I. Left column–Phase stack of the cross-correlation between the transmission spectra and the Cr I mask. The nights N1, N2, N3, and N4 are shown from top to bottom, respectively. In each plot, the orange dashed lines mark the orbital phases corresponding to the second and third contacts. Similarly, the red dashed lines mark the first and fourth contacts. The white solid line traces the expected planetary signal. Right column–Kp−υsys planes corresponding to the CCFs shown in the left column. In each plot, the dashed cross marks the expected Kp and υsys of the system.

Fig. A.3 Same as Figure A.2 but for the Fe I mask.

Fig. A.4 Same as Figure A.2 but for the Na I mask.

Fig. A.5 Same as Figure A.2 but for the Ti I mask.

Fig. A.6 Same as Figure A.2 but for the VALD mask devoid of the Cr I, Fe I, Na I, and Ti I spectral lines.

Fig. A.7 Same as Fig. 10 but considering all the available GIANO-B nights.

Fig. A.8 Posterior distribution of the investigated parameters in the DE-MCMC analysis of the He I triplet. The excess of absorption c [%], offset [%], peak position, and FWHM correspond to the parameters we used in the Gaussian fit, while the jitter term σj, the semi-amplitude of the correlated noise h, and the correlation length λ were used to parametrize the SE kernel within the GP.

Fig. A.9 Gaussian process correction. Right panel: Transmission spectrum centred on the He I triplet (in the planet rest frame) with the GP regression model overplotted and with the 1 σ uncertainty intervals (in blue) and the Gaussian+GP model (in red). Left panel: Final transmission spectrum after removing the GP model. The error intervals for the Gaussian fit were computed by displaying 1000 Gaussian fits within the 1 σ uncertainties of the derived parameters spanning the 16%-84% quantiles. Vertical black dotted lines indicate the position of the He I triplet.

Appendix B Cross-correlation analysis of the Hα line

The transmission spectrum of the Hα line revealed the presence of a clear emission feature during N4 (see Fig. 8). To try to understand the origin of this feature, we replicated the same CCF analysis discussed in Sect. 5.2 by creating a binary mask containing only the Hα line. The result for each night is shown in Fig. B.1. As expected, the cross-correlation analysis revealed the presence of a strong absorption feature in N1, N2, and N3 well aligned in the stellar rest frame and with a S/N above five in all three nights. The peaks of these signals are not in correspondence with the expected K_p of the system. From the cross-correlation analysis of N4, we observed the emission feature starting from around mid-transit. Moreover, even on this night, the signal does not follow the track expected by the planet but is well-aligned in the stellar reference system.

Fig. B.1 Same as Fig. A.2 but for a binary mask containing only the Hα line.

References

All Tables

All Figures

Fig. 1 Results of the fit of the transit light curves (LCs). Top: phase-folding of the de-trended data (top panel) and the corresponding O–C diagram (bottom panel). For clarity, in each panel, we show the binned data with black dots. Bottom: power spectral analysis of the TESS photometry. The grey dashed line shows the power spectral density (PSD) of the data after removing the best-fit transit model, while the blue solid line represents the combined PSD of the best-fit Gaussian process (GP) (the shoulder at f ≲ 1 day−1) and the white noise in the data (the plateau at f ≳ 10 day−1). The black line is the PSD of the residuals of the best fit, showing that the power excess at low frequencies has been effectively removed by the GP in the model.

In the text

Fig. 2 Rossiter-McLaughlin effect analysis. Top: RV time series taken during the transit on 2020 May 26 (N1). The best-fit model is superimposed, and the corresponding residuals are shown in the lower panel. The dashed black lines indicate the points of the first and fourth contact of the transit, while the dotted black lines represent the points of the second and third contact. Middle: same but for the transit on 2021 June 10 (N3). Bottom: same but for the transit on 2023 July 10 (N4).

In the text

	Fig. 3 Activity indicator log $R_{\text{HK}}^{\prime }$ as a function of the time for each night. The continuous black lines indicate the points of first and fourth contact of the transit, while the dashed black lines represent the points of second and third contact.
In the text

Fig. 4 Cross-correlation function for N1 (first row), N3 (second row), and N4 (third row) extracted from the DRS using a G2 mask. Left: cross-correlation function linear profiles as a function of time. The black line indicates the average CCF. Right: two-dimensional maps of the CCFs as a function of orbital phase and RV in the stellar rest frame. The four horizontal black lines show the times of transit contacts. The straight white line shows the expected planetary Keplerian velocity Kp using a planetary mass value of ~0.34 MJ (Zhou et al. 2017).

In the text

	Fig. 5 Generalised Lomb-Scargle periodogram of the five TESS sectors and the joint LC (bottom box). The vertical dashed line marks the orbital period of HAT-P-67 b reported in Table 3.
In the text

	Fig. 6 Composite master-out spectrum normalised to unity around the Na I doublet lines. Deep and narrow interstellar features peak out from the wider stellar lines.
In the text

	Fig. 7 Planes of Kp−υsys obtained by combining all nights for the Cr I (top left), Fe I (top right), Na I (bottom left), and Ti I (bottom right) templates. In each plot, the dashed cross marks the expected Kp and vsys of the system.
In the text

	Fig. 8 Transmission spectra of HAT-P-67 b for each night. The spectra are centred around the Hα line in the planetary rest frame (light grey) and also binned 20 times (in black circles). The red line is the MCMC Gaussian fit performed by SLOPpy, while the vertical blue dashed line indicates the rest frame transition wavelength of the Hα line.
In the text

	Fig. 9 Master-after, master-before, and master-in spectra in the star frame. Panels a–c represent the GIANO-B nights considered individually. Panel d shows all three transits together, while in panel e only the two observations with the highest S/N are considered. Vertical lines indicate the position of the He I. In- and before-transit absorption is visible by eye.
In the text

Fig. 10 Transmission spectroscopy near the He I triplet. Panels a and b show the 2D transmission spectroscopy maps in the stellar and planet rest frames, respectively. Tilted and vertical lines denote the position of the He I triplet, while horizontal dotted lines mark the transit contact points, t1, t2, t3, and t4. Panel c shows the spectroscopic light curve computed in a 0.075 nm band centred around the peak of absorption. Vertical dashed lines mark the transit contact points. Panel d presents the full in-transit averaged transmission spectrum.

In the text

	Fig. A.1 Variation of airmass (left) and S/N (right) for each night of HAT-P-67 b. The S/N was extracted from the FITS header of the HARPS-N spectra on the 53rd order, which contains the sodium feature. The continuous and the dashed vertical black lines represent the four points of contact of the transit.
In the text

Fig. A.2 Cross-correlation signatures of Cr I. Left column–Phase stack of the cross-correlation between the transmission spectra and the Cr I mask. The nights N1, N2, N3, and N4 are shown from top to bottom, respectively. In each plot, the orange dashed lines mark the orbital phases corresponding to the second and third contacts. Similarly, the red dashed lines mark the first and fourth contacts. The white solid line traces the expected planetary signal. Right column–Kp−υsys planes corresponding to the CCFs shown in the left column. In each plot, the dashed cross marks the expected Kp and υsys of the system.

In the text

	Fig. A.3 Same as Figure A.2 but for the Fe I mask.
In the text

	Fig. A.4 Same as Figure A.2 but for the Na I mask.
In the text

	Fig. A.5 Same as Figure A.2 but for the Ti I mask.
In the text

	Fig. A.6 Same as Figure A.2 but for the VALD mask devoid of the Cr I, Fe I, Na I, and Ti I spectral lines.
In the text

	Fig. A.7 Same as Fig. 10 but considering all the available GIANO-B nights.
In the text

	Fig. A.8 Posterior distribution of the investigated parameters in the DE-MCMC analysis of the He I triplet. The excess of absorption c [%], offset [%], peak position, and FWHM correspond to the parameters we used in the Gaussian fit, while the jitter term σj, the semi-amplitude of the correlated noise h, and the correlation length λ were used to parametrize the SE kernel within the GP.
In the text

Fig. A.9 Gaussian process correction. Right panel: Transmission spectrum centred on the He I triplet (in the planet rest frame) with the GP regression model overplotted and with the 1 σ uncertainty intervals (in blue) and the Gaussian+GP model (in red). Left panel: Final transmission spectrum after removing the GP model. The error intervals for the Gaussian fit were computed by displaying 1000 Gaussian fits within the 1 σ uncertainties of the derived parameters spanning the 16%-84% quantiles. Vertical black dotted lines indicate the position of the He I triplet.

In the text

	Fig. B.1 Same as Fig. A.2 but for a binary mask containing only the Hα line.
In the text