High-resolution transmission spectroscopy of the hot-Saturn HD 149026b

Advances in modern technologies enable the characterisation of exoplanetary atmospheres, most efficiently exploiting the transmission spectroscopy technique. We performed visible (VIS) and near infrared (nIR) high-resolution spectroscopic observations of one transit of HD 149026b, a close-in orbit sub Saturn exoplanet. We first analysed the radial velocity data, refining the value of the projected spin-orbit obliquity. Then we performed transmission spectroscopy, looking for absorption signals from the planetary atmosphere. We find no evidence for H$\alpha$, Na $\scriptstyle\rm I$  D2$-$D1, Mg $\scriptstyle\rm I$ and Li $\scriptstyle\rm I$ in the VIS and metastable helium triplet He $\scriptstyle\rm I$(2³S) in the nIR using a line-by-line approach. The non-detection of HeI is also supported by theoretical simulations. With the use of the cross-correlation technique, we do not detect Ti $\scriptstyle\rm I$, V $\scriptstyle\rm I$, Cr $\scriptstyle\rm I$, Fe $\scriptstyle\rm I$ and VO in the visible, and CH₄, CO₂, H₂O, HCN, NH₃, VO in the nIR. Our non-detection of Ti $\scriptstyle\rm I$ in the planetary atmosphere is in contrast with a previous detection. We performed injection-retrieval tests, finding that our dataset is sensitive to our Ti $\scriptstyle\rm I$  model. The non-detection supports the Ti $\scriptstyle\rm I$  cold-trap theory, which is valid for planets with $T_{\rm eq}<$ 2200 K like HD 149026b. Even if we do not attribute it directly to the planet, we find a possibly significant Ti $\scriptstyle\rm I$  signal highly redshifted ($\simeq$+20 $\mathrm{km\,s}^{-1}$) with respect to the planetary restframe. Redshifted signals are also found in the Fe $\scriptstyle\rm I$  and Cr $\scriptstyle\rm I$  maps. While we can exclude an eccentric orbit to cause it, we investigated the possibility of material accretion falling onto the star, possibly supported by the presence of strong Li $\scriptstyle\rm I$  in the stellar spectrum, without finding conclusive results. The analysis of multiple transits datasets could shed more light on this target.

1 Introduction

With advancements in instruments and facilities, the scientific community’s focus has shifted from the discovery of exoplanets to their characterization, with one of the main goals the understanding of their atmospheric compositions. Presently, transmission spectroscopy is the best technique to characterize the properties of the atmospheres of exoplanets (e.g., Birkby, 2018; Madhusudhan, 2019). During transits, exoplanetary atmospheres absorb specific wavelengths of the radiation emitted by their host stars, providing insights into their elemental composition in the stellar light coming to us. High-resolution spectroscopy (resolving power R $\gtrsim$ 20,000) enables us to glean a wealth of information about exoplanet parameters and atmospheric content. This technique allows us to pinpoint atoms and molecules thanks to their unique firm isolating their distinct high cross-section lines. As an example, we are able to identify Na $\scriptstyle\rm I$  doublet (e.g., Wyttenbach et al., 2015), Ca $\scriptstyle\rm II$ H & K lines (e.g., Yan et al., 2019), H$\alpha$ (e.g., Jensen et al., 2012) and the metastable He $\scriptstyle\rm I$  triplet (e.g., Allart et al., 2018). Alternatively, the use of the cross-correlation function (CCF) (e.g., Snellen et al., 2010; Brogi et al., 2012; Hoeijmakers et al., 2018; Borsa et al., 2022; Prinoth et al., 2023) offers another powerful method to detect species for which strong lines can not be identified. These types of analyses are mainly employed in the characterization of the atmospheres of the so called hot Jupiters (Bell & Cowan, 2018), i.e., gas giants in close-in orbit. Due to their proximity to the host star, these planets experience intense XUV (10-912 Å) and FUV (912-2,585 Å) irradiation, leading to the heating and possibly inflation of their atmospheres (e.g., Owen, 2019; Biassoni et al., 2024).

Among the class of gas giants, a particularly interesting case is that of HD 149026b. The planet was first detected using the radial velocity method by Sato et al. (2005), who measured a radius and mass (in Jupiter’s units) of $R_{p}=0.725\pm 0.05\,R_{J}$, $M_{p}=0.36\pm 0.02\,M_{J}$, respectively. HD 149026b is then more massive but smaller in size than Saturn, despite receiving intense irradiation from its parent star, a fact that should lead to an inflated atmosphere (Sato et al., 2005). The planet’s high density, such as $\rho_{p}=1.18_{-0.30}^{+0.38}$  g cm^-3 (Wolf et al., 2007), $\rho_{p}=1.59_{-0.36}^{+0.38}$  g cm^-3 (Torres et al., 2008), $\rho_{p}=2.1\pm 0.8$  g cm^-3 (Southworth, 2010), suggests a core rich in metals, as proposed by Sato et al. (2005), who estimated a core mass $\approx 67\,M_{\oplus}$, with $M_{\oplus}$ is the Earth’s mass. Subsequent studies echoed these high core estimations. Fortney et al. (2006) claimed a significant presence of heavy elements in both the core and envelope, estimating values $\approx 60-93\,M_{\oplus}$; Burrows et al. (2007) inferred a core mass ranging between $80-110\,M_{\oplus}$; Baraffe et al. (2008) and Carter et al. (2009) estimated a heavy element core mass of $\approx 60-80\,M_{\oplus}$ and $\approx 45-70\,M_{\oplus}$, respectively.

A combination of both a high core mass and high stellar metallicity (${[Fe/H]}=0.36\pm 0.05$, Sato et al., 2005) appears to be crucial prerequisites for a large atmospheric metallicity of the planet (Zhang et al., 2018). Several atmospheric models have been proposed to unravel the composition of HD 149026b. Fortney et al. (2006) observed a hot stratosphere attributed to the absorption of stellar flux by TiO and VO, while Stevenson et al. (2012), analysing Spitzer secondary eclipse observations, adopted a chemical equilibrium model and suggested the presence of large amounts of CO and CO₂, with moderate heat redistribution, high metallicity, and absence of thermal inversion on the dayside. Zhang et al. (2018) analysed Spitzer phase curves observations at 3.6 and 4.5 $\mu$m and inferred a high albedo, possibly explained by the presence of reflective cloud layers in the planet’s upper atmosphere. Recently analyses based on transmission spectroscopy conducted with the High Dispersion Spectrograph (Noguchi et al., 2002) on the Subaru telescope were performed by Ishizuka et al. (2021). They reported a tentative $\simeq 4.4\sigma$ detection of Ti $\scriptstyle\rm I$, and a marginal $\simeq 2.8\sigma$ detection of Fe $\scriptstyle\rm I$, alongside non detections of Sc $\scriptstyle\rm I$, V $\scriptstyle\rm I$, Cr $\scriptstyle\rm I$, Mn $\scriptstyle\rm I$, Co $\scriptstyle\rm I$  and TiO. The presence of Ti $\scriptstyle\rm I$  without TiO suggests a supersolar C/O ratio (Ishizuka et al., 2021). Recently, Bean et al. (2023) analysed the HD 149026b dayside emission obtained with the James Webb Space Telescope (JWST). Similarly to prior studies, they concluded that the planet exhibits a highly super-solar metallicity ($[M/H]=2.09^{+0.35}_{-0.32}$) and a high carbon-to-oxygen ratio ($[C/O]=0.84\pm 0.03$). These results further support the detection of CO₂ and H₂O in HD 149026b. The same dataset was analysed by Gagnebin et al. (2024) using a 1D radiative-convective-thermochemical equilibrium models from which new values of planet atmospheric metallicity are estimated to be $\simeq 1.15$ if VO isn’t included in the models, and $\simeq 1.30$ if it is. These values are $\simeq$ 10 times smaller than the results of Bean et al. (2023).

Spinelli et al. (2023) updated the HD 149026b parameters by exploiting more precise GAIA DR2 (Gaia Collaboration et al., 2018) distances, inferring a mass and radius of $M_{p}=0.28\pm 0.03\,M_{J}$ and $R_{p}=0.74\pm 0.02\,R_{J}$, respectively. These values yield a planetary density of $\rho_{p}=0.86\pm 0.09$   g cm^-3, consistent with the lowest estimates of previous works, e.g., $\rho_{p}=0.750^{+0.089}_{-0.077}$ g cm^-3 as in Bonomo et al. (2017), and $\rho_{p}=0.85^{+0.10}_{-0.09}$ g cm^-3 as in Carter et al. (2009). Given the results by Spinelli et al. (2023), we decided to shed more light on the atmospheric content of HD 149026b with high-resolution transmission spectroscopy. To this aim, we collected data from two transits at the Telescopio Nazionale Galileo (TNG) equipped with high-resolution visible and near-infrared spectrographs, i.e., the High Accuracy Radial velocity Planet Searcher North (HARPS-N) (R $\simeq$ 115,000) and GIANO-B (R $\simeq$ 50,000), respectively.

In this work, we describe our observations, data analysis and scientific results. The paper is organized as follows: details regarding the observations can be found in § 2, while § 3 and § 4 cover the modelling of Rossiter-McLaughlin effect, center-to-limb variations and data analyses. Discussion of the results of the data analysis and concluding remarks are detailed in the final § 6.

Figure 1: HARPS-N and GIANO-B SNR for the second night evaluated at 5,500 and 16,300 Å, respectively. The red dots represent the SNR of the spectra used for transmission spectroscopy analysis, while the blue ones are discarded. The two vertical orange lines denote the T₁ and T₄ transit contact points.

2 Observations

We observed one transit of HD 149026b on June 14th, 2022 with the HARPS-N and GIANO-B high-resolution spectrographs, mounted at the Telescopio Nazionale Galileo (Program A45TAC_30, PI Borsa). We exploited the GIARPS configuration (Claudi et al., 2016), which allows us to observe simultaneously with the two spectrographs. A second scheduled transit on 3rd April 2022 was lost due to bad weather. The wavelength range of the acquired spectra span from visible (VIS, 3,900 - 6,900 Å, HARPS-N) to near-infrared (nIR, 9,400 - 24,200 Å, GIANO-B). With GIANO-B, we collected 78 spectra with exposures of 200 sec per nodding position (ABAB pattern). Meanwhile with HARPS-N we observed a total of 64 spectra with an exposure time of 300 sec. All the spectra were reduced with the standard DRS pipelines v3.7 and v1.6.1 for HARPS-N (Cosentino et al., 2012) and GIANO-B (Rainer et al. (2018); Harutyunyan et al. (2018)), respectively. A summary of the observations is shown in Table 1.

During the observations, we lost about half an hour at the transit ingress because of a problem with the instrument guiding system. A strong loss of flux was experienced also during the transit, which was solved by performing a repointing procedure. Figure 1 shows the signal-to-noise ratio (SNR) obtained with HARPS-N and GIANO-B at 5,500 and 16,300 Å, respectively. For our analysis, we discarded all the spectra with low SNR with respect to the adjacent ones (see Fig. 1), ultimately performing our analysis on a total of 61 (75) HARPS-N (GIANO-B) spectra.

Table 1: log of the observations performed with the two instruments^(a).

Instrument	# spectra (Out/In)	¡SNR¿
HARPS-N	61 (30/31)	52.0
GIANO-B	75 (37/38)	53.6

(a) Parenthesis specify the spectra collected in and out of transit. The average SNR is taken at 5,500 Å and 16,300 Å for HARPS-N and GIANO-B respectively.

Table 2: physical parameters for the HD 149026 system.

Parameter	Symbol	Value [Unit]
Stellar Parameters
Age (a)		$2.60\pm 0.20\,\,$[Gyr]
Effective Temperature (g)	$T_{eff}$	$6075\,\,$[K]
Spectral Class (c)		G0
Stellar Mass (h)	$M_{\star}$	$1.09\pm 0.13\,\,$[M⊙]
Stellar Radius (h)	$R_{\star}$	$1.47\,\,$[R⊙]
Projected Rotation Speed (f)	$v\sin{i}$	$6.10_{-0.48}^{+0.47}\,\,$[km s-1]
Surface Gravity (b)	$\log{g}$	$4.37\pm 0.04\,\,$[cm s-2]
Metallicity (a)	$[Fe/H]$	$0.36\pm 0.08\,\,$ [dex]
Planetary Parameters
Planet Mass (h)	$M_{p}$	$0.28\pm 0.03\,\,$[MJ]
Planet Radius (h)	$R_{p}$	$0.74\pm 0.02\,\,$[RJ]
Equilibrium Temperature (d)	$T_{eq}$	$1634^{+90}_{-23}\,\,$[K]
Planet Density (h)	$\rho_{p}$	$0.86\pm 0.09\,\,$[g cm-3]
Orbital Parameters
Epoch (f)	$T_{0}$	$2454597.702816162_{-0.004317}^{+0.003742}\,\,$[BJD]
Period (a)	$P$	$2.8758916\pm 0.0000014\,\,$[days]
Transit Duration (e)	$T_{14}$	$3.23\pm 0.15\,\,$[h]
Systemic Velocity (f)	$\gamma$	$-18.0319_{-0.0008}^{+0.0008}\,\,$[km s-1]
Radial Velocity Amplitude (f)	$K_{\star}$	$38.3\pm 3.1\,\,$[m s-1]
Semi-major Axis (h)	$a$	$0.041\pm 0.002\,\,$ [AU]
Inclination (f)	$i$	$84.1130^{+0.5259}_{-0.4986}\,\,$[deg]
Eccentricy (a)	$e$	$<0.013\,\,$
Projected Obliquity (f)	$\lambda$	$2.5040_{-5.3317}^{+5.3640}\,\,$[deg]
(a) Bonomo et al. (2017); (b) Stassun et al. (2017); (c) Ment et al. (2018); (d) Torres et al. (2008)
(e) Albrecht et al. (2012); (f) This work; (g) Gaia Collaboration et al. (2018); (h) Spinelli et al. (2023)

3 Rossiter-McLaughlin effect

During a transit, the partial occultation of the stellar disk by the planet induces a distortion in the observed stellar spectrum, leading to apparent variations of the derived stellar radial velocities (RVs). This specific phenomenon is commonly referred to as the Rossiter-McLaughlin (RML) effect (Rossiter, 1924; McLaughlin, 1924). The amplitude and shape of the RML effect are mainly dependent on the projected stellar rotation velocity ($v\sin{i}$), projected spin-orbit inclination, impact parameter and $R_{p}/R_{\star}$ ratio, where $R_{\star}$ is the stellar radius. An analysis of the RML effect enables the extraction of these planetary and stellar parameters which are later used in the transmission spectroscopy (Sect. 4.3). In order to model the RML effect, we employed the CaRM code, a semi-automated tool (Cristo et al., 2022) that employs an MCMC algorithm. Details of our analysis of the RML effect are presented in Appendix A. Globally, the values of projected rotational speed, systemic velocity, radial velocity amplitude and inclination we obtain are in agreement with literature results. In particular, we refine the projected spin-orbit angle $\lambda=2.5\pm 5.3$ deg, in agreement with the literature value ($\lambda=12\pm 7$ deg, Albrecht et al., 2012) but with smaller error-bars.

Figure 2: Result from the RML fitting procedure performed with CaRM. The red line is the best fit model, the grey shadows are a random sample of the posterior distribution. The yellow regions show the ingress and egress of the transit, while the orange region shows where the planet is fully transiting the stellar disk.

4 Transmission spectroscopy

The main goal of our observations is to look for the atmosphere of HD 149026b using high-resolution transmission spectroscopy. We started our analysis from the merged 1D spectra (s1d) provided by the HARPS-N DRS. Similarly, we analysed the 71^th GIANO-B spectral order searching for the presence of the metastable helium triplet He $\scriptstyle\rm I$(2³S) at $\simeq 10,830$ Å. We chose to work with the echelle spectral orders (ms1d) in the nIR band instead of the merged s1d spectra to improve the normalisation of the GIANO-B spectra.

4.1 Telluric correction and wavelength calibration

We performed correction for telluric H₂O and O₂ lines on all HARPS-N spectra using Molecfit v. 4.3.1 (Smette et al., 2015; Kausch et al., 2015), following the guidelines specified in Allart et al. (2017) but slightly modifying the wavelength ranges of the correction. The left panel of Fig. 3 provides an example of telluric correction obtained using Molecfit, near the H$\alpha$ line.

We used Molecfit v. 4.3.1 also to correct the GIANO-B spectra, but following a different recipe. The GIANO-B wavelength calibration may be slightly inconsistent between the different echelle orders due to the characteristics of the U-Ne calibration lamp (the number of useful emission lines varies greatly between the orders). While Molecfit accounts for some wavelength shift in the observed spectrum, trying to correct the whole merged GIANO-B spectra resulted in p-Cygni-like residuals due to the slight misalignment of the echelle orders. For this reason we worked with the ms1d GIANO-B spectra, where the echelle orders are still separated, and we ran Molecfit on each order independently. Another complication arose from the fact that the nIR wavelength range is heavily affected by the telluric absorption: due to the combination of the wide wavelength range of GIANO-B and the large telluric contribution, the computational time needed to model and remove the telluric lines from all the 50 orders of a single GIANO-B spectrum is quite high. In order to speed up the process, we computed the model only once every ten spectra, assuming that the atmospheric conditions did not vary too much in this time. To account for any possible variation due to the position of the star on the slit (A and B position of the nodding observational strategy), we worked separately on the A and B spectra, computing one model for every ten A spectra and one for every ten B spectra. While we globally removed H₂O, O₂, CO₂, N₂O, CH₄ and CO from the whole spectrum, working on the separate orders allowed us to model only the atmospheric molecules present in each echelle order. Right panel of Fig. 3 provides an example of telluric correction near the He $\scriptstyle\rm I$(2³S) lines.

Figure 3: Example of Molecfit telluric correction. The left figure refers to the HARPS-N telluric correction in proximity of the H$\alpha$ line, while right figure to the GIANO-B telluric correction in proximity of the He $\scriptstyle\rm I$(2³S). Both spectra are taken from the first exposure.

4.2 Transmission spectrum extraction

We extracted the transmission spectrum following the guidelines outlined in Wyttenbach et al. (2015), performing the procedure in a short wavelength interval bracketing each of the atmospheric lines investigated (see Table 3). As a first step, all the spectra were normalized and shifted to the stellar reference frame. Then, we generated an average reference stellar spectrum (Master-Out, $M_{\rm out}$), by performing a weighted mean of all the out-of-transit spectra. Each spectrum was then divided by $M_{\rm out}$ and the residual spectra $R_{i}$ derived accordingly. These $R_{i}$ were further normalized to eliminate any trend possibly given by imperfect continuum atmospheric dispersion corrections. In order to maximize any potential planetary signal (if present), all $R_{i}$ were shifted into the planet reference frame, considering a circular orbit with parameters form Table 2. Then we computed a weighted mean of all the in-transit residual spectra to obtain the final average transmission spectrum. These transmitted spectra have shorten wavelength range with respect to the initial intervals chosen for each element analysed; this is due to the Doppler shifts and subsequent interpolation procedures adopted during the transmission spectroscopy analysis that bring divergence at the edges of the wavelength ranges (Table 3). After each interpolation the wavelength bin width was maintained constant at the value of 0.01 Å.

Table 3: interval used for the transmission spectroscopy analysis.

Element	Initial $\lambda$ range [Å]	Final $\lambda$ range [Å]
MgI	5,125-5,195	5,140-5,190
NaI D2-D1	5,850-5,930	5,885-5,905
H$\alpha$	6,525-6,595	6,550-6,580
LiI	6,670-6,720	6,700-6,715
HeI(23S)	10,800-10,880	10,825-10,840

Figure 4: CLV+RM models (red) compared with transmission spectra (light grey) around the H$\alpha$ and Na $\scriptstyle\rm I$  D doublet wavelength regions. Green lines show the position of the theoretical values of H$\alpha$, Na $\scriptstyle\rm I$  D2 and D1 lines. The CLV+RML model contamination is well within the data noise.

4.3 Single lines analysis

We extracted the transmission spectrum in the wavelength ranges of Mg $\scriptstyle\rm I$  triplet (5,167.32 - 5,172.68 - 5,183.60 Å), Na $\scriptstyle\rm I$  D2-D1 (5,889.95 - 5,895.92 Å), H$\alpha$ (6,562.81 Å), Li $\scriptstyle\rm I$  (6,707.76 Å) in the optical band, and of the He $\scriptstyle\rm I$(2³S) triplet (10,829.09 - 10,830.25 - 10,830.34 Å) in the infrared band.

During a transit, one must consider that the host star is not an homogeneously bright disc, rather its surface brightness changes with the the distance from disc centre (centre-to-limb variations, CLV). Moreover, the star rotates. These stellar properties affect the emitted spectrum occulted during transit through CLV and RML effect, respectively, that, in turn can affect the transmission spectrum in various ways, e.g., by modifying the shape of line profiles and/or by causing false atmospheric detections (e.g., Yan et al., 2017; Borsa & Zannoni, 2018; Casasayas-Barris et al., 2020). Given our target’s characteristics of low projected rotational velocity and small $R_{p}/R_{\star}$ ratio, we expect minimal-to-low contamination from RML and CLV. In order to verify if this is indeed the case, we modelled CLV and RML as outlined in Yan et al. (2017) on few lines (H$\alpha$, Na $\scriptstyle\rm I$  doublet and Mg $\scriptstyle\rm I$  triplet), using the same methodology employed in Borsa et al. (2021). The star was modelled as a disc, and mapped on fine grid with 0.01 $R_{\star}$ resolution. At each point of the grid we calculated the angle between the normal to the stellar surface and the line of sight, $\mu=\cos{\theta}$, and the projected rotational velocity (rescaling the $v\sin{i}$ and using rigid body rotation). A spectrum was then assigned to each point of the grid by a quadratic interpolation on $\mu$, and Doppler-shifting according to the stellar rotation the model spectra created using the tool Spectroscopy Made Easy (SME, Piskunov & Valenti, 2017), with the line list from the VALD database (Ryabchikova et al., 2015) and ATLAS9 (Kurucz, 1993; Castelli & Kurucz, 2003) stellar atmospheric models. The model spectra were created with null rotational velocity for 21 different $\mu$ values at the resolving power of HARPS-N. Using the orbital information from Table 2, we then simulated the transit of the planet, calculating the stellar spectrum for different orbital phases as the average spectrum of the non-occulted modelled sections. In the last step, we divided each spectrum for a master stellar spectrum calculated out of transit, obtaining the information of the relevance of CLV+RML effects at each in-transit orbital phase. We then moved everything in the planetary rest frame, calculating the simulated CLV+RML effects on the transmission spectrum.

Figure 5: H$\alpha$ and Na $\scriptstyle\rm I$ D2-D1 transmission spectra. The grey points are our transmission spectra, the red line the best Gaussian fit while the black dots are a resample of the transmission spectra with bin size of 0.3 Å.

In Fig. 4 we show the comparison between the transmission spectra of H$\alpha$ and Na $\scriptstyle\rm I$  D doublet and the simulated CLV+RML contamination. The latter is well within the noise of the data, confirming expectations. Since CLV+RML are negligible in the transmission spectrum, we ignore them in our subsequent analyses.

We performed a Gaussian fit around all the lines of interest, using a linear regression method (see footnote in Appendix A). In the visible band, our analysis led to non-detection for all the considered lines, i.e., H$\alpha$, Na $\scriptstyle\rm I$  D2-D1, Mg $\scriptstyle\rm I$  triplet and Li $\scriptstyle\rm I$ . Examples of the results are shown in Fig. 5. We estimated upper limits of such lines by calculating the standard deviation within $\pm 0.5$ Å from their reference wavelengths (see Table 4).

Indeed, the depth of stellar lines strongly affects our ability to put strong limits on the atmospheric heights. We used $M_{\rm out}$ to measure the stellar flux intensity within the core of all these lines relative to the continuum. As shown in Table 4, most of the transmission signals we were looking for are affected by a stellar flux which is quite low. As a result, the SNR of the transmission spectrum of the planet is constrained by such low values of the flux, potentially making any possible planetary signal lost in the noise.

We also surveyed the nIR band in search of the metastable helium triplet He $\scriptstyle\rm I$(2³S) at $\simeq 10,830$ Å. This state of Helium predominantly appears in planets orbiting K-type stars (Oklopčić, 2019), and its abundance, and relative absorption signal, depends upon various factor, e.g., the hardness ratio ($L_{\rm XUV}/L_{\rm FUV}$) of the stellar spectrum, the planetary orbital distance, the planetary Hill’s radius, and the He-to-H abundance ratio in the planet’s outer atmosphere (e.g., Biassoni et al., 2024).

We applied the same technique used for the visible band to extract the transmission spectrum close to the He $\scriptstyle\rm I$(2³S) lines. The result shows substantial noise and a sinusoidal pattern (Fig. 6), which is a known artefact in GIANO-B transmission spectra (e.g. Guilluy et al., 2023). In order to improve the normalization, we fitted a sinusoidal function of the form $A\,\cdot\,\sin(x)\,+\,B\,\cdot\,\cos(bx)\,+\,mx\,+\,C$, for which we divided the transmission spectrum (grey curve, Fig. 6). During the analysis, we took into account the two telluric OH doublets (10,829.46-10,829.15 Å and 10,831.38-10,831.29 Å, wavelengths in air, Oliva et al. (2013)), whose emission can notably affect this wavelength region (e.g., Guilluy et al., 2023). We shifted them in the stellar reference system after correcting for the systemic velocity and the barycentric velocity of the Earth. In this reference system these lines fall at 10,829.84-10829.54 Å and 10,831.77-10,831.67 Å respectively. To remove them, we masked completely the regions between [10,832.38 - 10,832.89] Å and [10,834.53 - 10,834.81] Å.

Similarly to the analysis in the visible band, the helium profile was fitted with two Gaussian profiles like the Na $\scriptstyle\rm I$  doublet. The doublet lines at $10,830.25-10,830.34$ Å were considered blended (Kirk et al., 2022). The prior on the center of the doublet, $\mu=10,830.29$ Å was set between $\pm 0.2$ Å to prevent significant deviation from the fit due to the high noise level in the spectrum. The fit hinted at an absorption pattern near the He $\scriptstyle\rm I$(2³S) doublet (figure 7), but the profile lacked a clear definition, and its deviation from the continuum was only at $\simeq 2.3\,\sigma$.

The absence of a significant He $\scriptstyle\rm I$(2³S) signal in the planetary atmosphere was further investigated theoretically. Using the 1D photo-ionization hydrodynamic code ATES (Caldiroli et al., 2021) and the Transmission Probability Module (TPM) (Biassoni et al., 2024), we modelled the outflow of HD 149026b, estimating its mass loss rate alongside the number density of He $\scriptstyle\rm I$(2³S) and the transmission spectrum. To achieve this simulation, we employed the HD 149026 spectrum provided by Behr et al. (2023) rescaled at the planetary orbital distance assuming a circular orbit. In the simulation we used planetary and stellar parameters from Table 2. To match our observations, we convolved the simulated absorption profile with the GIANO-B instrumental resolution (R = 50,000, assumed as Gaussian), considering also the planet’s rotation under the assumptions of tidal locking and a radius equal to $R_{\rm eff}=R_{p}\sqrt{1+h/\delta}$, where $h$ represents the theoretical absorption depth at the He $\scriptstyle\rm I$(2³S) doublet core derived from TPM, and $\delta=(R_{p}/R_{\star})^{2}$ denotes the transit depth. The simulation produces a mass-loss rate of $\approx 5\times 10^{10}$ g s^-1, resulting in an absorption depth of He $\scriptstyle\rm I$(2³S) of $\simeq 0.7\%$ (blue curve in Figure 7), of the same order of the error-bars in the data.

Figure 6: He $\scriptstyle\rm I$(2³S) transmission spectrum. Left panel: in grey $F_{\rm in}/F_{\rm out}$ obtained with the transmission spectroscopy analysis. A sinusoidal shape is evident. Red line is the sinusoidal fit obtained using the relation described in the text. Right panel: the grey transmission spectrum divided by the sinusoidal red fit. Black dots represent the binned spectra with equally spaced bins of 0.3 Å, while vertical dashed green lines show the position of the theoretical He $\scriptstyle\rm I$(2³S) lines in vacuum.

Figure 7: He $\scriptstyle\rm I$(2³S) normalized transmission spectrum in comparison with the double Gaussian fit and TPM simulation.

Table 4: HD 149026b transmission spectroscopy upper limits.

Element	Theoretical line in Air	Flux in the core of the line	Upper limit (1$\sigma$)	Instrument
	[Å]	[%]	[%]
Mg $\scriptstyle\rm I$	5,167.32	16.01	0.86
	5,172.68	12.71	0.76	HARPS-N
	5,183.60	11.71	0.65
Na $\scriptstyle\rm I$ D2	5,889.95	10.03	0.45	HARPS-N
Na $\scriptstyle\rm I$ D1	5,895.92	12.15	0.40	HARPS-N
H$\alpha$	6,562.81	16.70	0.55	HARPS-N
Li $\scriptstyle\rm I$	6,707.76	89.44	0.29	HARPS-N
He $\scriptstyle\rm I$(23S)	10,829.09	100	0.23
	10,830.25	96.23	0.34	GIANO-B
	10,830.34

4.4 Cross-correlation analysis

Figure 8: Flux variation for all the exposures for the 59^th HARPS-N spectral order (6,538.7 - 6,612.3 Å). The strong absorption line at pixel $\simeq$ 1,200 is the stellar H$\alpha$ line. Top panel (before the normalisation) shows the 59^th spectral order along all the exposures. Bottom panel shows the result after the normalisation.

Line-by-line transmission spectroscopy is very effective in detecting absorption lines with large cross-sections in exoplanetary atmospheres. However, it often fails in detecting the numerous, intrinsically faint lines expected from atoms such as, e.g., $\rm Fe$, $\rm Ti$, $\rm V$, or from molecules. In order to increase SNR, the cross correlation function (CCF) is indeed a more powerful tool (e.g., Snellen et al., 2010; Brogi et al., 2012; Hoeijmakers, H. J. et al., 2019; Borsa et al., 2021).

In our analysis, we employed the CCF technique to investigate the presence of Ti $\scriptstyle\rm I$, V $\scriptstyle\rm I$, Cr $\scriptstyle\rm I$, Fe $\scriptstyle\rm I$  and VO in the visible band, and CH₄, CO₂, H₂O, HCN, NH₃, and VO in the nIR. For this purpose, we constructed the templates employing petitRADTRANS (Mollière, P. et al., 2019), a python package that allows to calculate the planetary transmission and emission spectra. We constructed models by using the parameters of Table 2, with an isothermal atmospheric profile with T = 2000 K, a metallicity [Fe/H] = 0.36 (Ishizuka et al., 2021) and a continuum pressure of 10 mbar, which is consistent with the white light transit depth. The models, generated as expected variation in planetary radius with wavelength, were then translated in ($R_{p}/R_{\star}$)² and convolved with the HARPS-N resolving power. We normalized the templates and convolved them also to match the planet’s rotational velocity, assuming tidal locking.

For the CCF analysis, we started performing similar steps as done for the line-by-line analysis on the s1d spectra, but using the e2ds spectra (e.g., Stangret et al., 2020). This allows us to perform a better normalization when working on the whole spectral range of the spectrograph. We employed Molecfit to correct the 1D s1d HARPS-N spectra and applied the retrieved telluric profile to correct the e2ds spectral orders (e.g., Hoeijmakers et al., 2020). To enhance the quality of the CCF analysis, we excluded the first five orders and the last one from all exposures due to their low SNR, focusing on the wavelength range 4,000-6,840 Å. The e2ds spectral orders were then reformatted into matrices, each of them contain all the exposures of the same order (Fig. 8, top panel). As we are working with ground-based instruments, we must account for the flux variation over time. Therefore, we normalized each exposure by dividing it by the mean value across the pixels (Fig. 8, lower panel).

We then started our cross-correlation function calculation,

$$CCF(v,t)=\frac{\sum_{k}{x_{k}(t)\cdot T_{k}(v)}}{\sum_{k}T_{k}(v)},$$ (1)

where $v$ is the velocity at which the template is shifted and $k$ refers to the pixels in our one-dimensional spectrum $x$ that contains all the stacked normalised and telluric corrected spectral orders, while $T$ represents the template. The CCF is done in the range -200 $\mathrm{km\,s}^{-1}$¡ v ¡ 200 $\mathrm{km\,s}^{-1}$, in steps of 1 $\mathrm{km\,s}^{-1}$. The resulting CCF was shifted to the stellar reference frame, interpolated within the RV range -150 $\mathrm{km\,s}^{-1}$¡ v ¡ 150 $\mathrm{km\,s}^{-1}$in steps of 1 $\mathrm{km\,s}^{-1}$and normalized by the CCF $M_{\rm out}$. For a more accurate normalization, we further divided this CCF by the median values across the exposures. To eliminate low-frequency fluctuations and smoothing the CCF, we conducted a Fast Fourier Transform cutting any velocity beyond 100 $\mathrm{km\,s}^{-1}$for each exposure. The resulting final CCFs for the inspected elements are given in Fig. 10 and Fig. 11. The blurred region at phase $\simeq\,-0.02$ is caused by the lack of data, since we discarded them to their low SNR.

To detect any possible planetary atmospheric signal, we generated the $K_{p}-V_{\rm sys}$ maps by shifting all the CCFs in different planetary reference systems. We considered circular orbits with various $K_{p}$ velocities in the range 0-300 $\mathrm{km\,s}^{-1}$, with a step of 1 $\mathrm{km\,s}^{-1}$, and averaged all the in-transit residual CCFs. $K_{p}-V_{\rm sys}$ maps are then divided by the standard deviation of the overall map, after excluding a region of $\pm 100$ $\mathrm{km\,s}^{-1}$for both RV and $K_{p}$ values centered at the expected $K_{p}$ value (156.19 $\mathrm{km\,s}^{-1}$) and 0 $\mathrm{km\,s}^{-1}$  RV. This normalization yields $K_{p}-V_{\rm sys}$ maps in terms of SNR as shown in Fig. 10 and Fig. 11. The blurred regions in the CCFs at the beginning of the transit are due to the lack of data for these orbital phases (see Fig. 1 and Sect. 2).

From the analysis of our data, we do not find evidence for the presence of any of the searched species (SNR $>4$) in the atmosphere of the planet in proximity of the expected planetary signal. The non-detection of Ti $\scriptstyle\rm I$, in particular, is in contrast with the previous detection of this element by Ishizuka et al. (2021). In order to verify if with our data we are sensitive to the expected signals, we performed the injection of our planetary atmospheric models into the data prior to any analysis, using a planetary orbital velocity opposite to the theoretical one (e.g., Pelletier et al., 2021). The e2ds spectra were then analysed following the same prescription described at the beginning of Sect. 4.4. CCF and $K_{p}-V_{\rm sys}$ maps containing the injected models are shown in Fig. 10. In this way, we retrieved the injected signal of Ti $\scriptstyle\rm I$  and Fe $\scriptstyle\rm I$  at the expected $K_{p}$ and RV value at $\sim 4.1\,\sigma$ ($K_{p}=-174$ $\mathrm{km\,s}^{-1}$, RV = 0 $\mathrm{km\,s}^{-1}$) and $\sim 3.3\,\sigma$ ($K_{p}=-152$ $\mathrm{km\,s}^{-1}$, RV = $-1$ $\mathrm{km\,s}^{-1}$) respectively. For V $\scriptstyle\rm I$  and Cr $\scriptstyle\rm I$  the data are not sensitive to our planetary model to give us a clear detection.

Since the results of CCF analysis may strongly depend on the line lists used to create the templates (e.g., Gandhi et al., 2020), we conducted a check on our Ti $\scriptstyle\rm I$  non-detection by creating templates with different line lists. We used the linelists from Kurucz (Kurucz, 1993), VALD (Piskunov et al., 1995), and NIST downloaded from the DACE database and converted in petitRADTRANS format to create new models, and the isothermal templates at 2000 K, 2500 K, and 3000 K from Kitzmann et al. (2023). We then re-performed the CCFs by using all these templates, obtaining very similar results and non-detections each time.

We performed CCFs also in the infrared channel on the ms1d GIANO-B spectral orders, searching for the presence of molecules such as CH₄, CO₂, H₂O, HCN, NH₃, and VO. Also in this channel we performed telluric correction using Molecfit, with a dedicated procedure as described in Sect. 4. Many orders are not usable due to their low SNR and the intrinsic transmissivity of the Earth’s atmosphere. Therefore, we masked the flux of the following GIANO-B diffraction orders (81, 80, 79, 69, 68, 67, 66, 57, 56, 55, 54, 53, 52, 51, 43, 42, 41, 40, 39, 38, 37, 32) before the CCF analysis for each exposure. Figure 9 shows an example of two different orders where the Molecfit telluric correction was applied successfully and not, respectively.

The CCF and $K_{p}-V_{\rm sys}$ maps for the selected models where then retrieved with the same normalisation procedure employed in the VIS band with the HARPS-N e2ds spectral orders. Also in the nIR band, our analysis did not show any detection of exoplanetary absorption signal (Fig. 12).

Figure 9: Comparison between two different GIANO-B ms1d spectral orders, before (navy) and after (red) Molecfit correction. Left panel refers to the 34^th spectral order, where the correction is successfull. Right panel refers to the 41^th spectral order, where correction does not work properly as in this wavelength region telluric lines are saturated. We didn’t consider the 41^th spectral order and similar ones in our analyses and discarded them.

Figure 10: CCF (left) and $K_{p}-V_{\rm sys}$ (center) maps for Ti $\scriptstyle\rm I$  and Fe $\scriptstyle\rm I$. The horizontal white lines in CCF maps correspond to the T₁ and T₄ contact points, while the cyan dash dotted and gold dotted lines mark the expected planetary $K_{p}$ and -$K_{p}$ of the injected signals, respectively. The white slanted line corresponds to the Doppler shadow curve. Cyan and gold arrows in $K_{p}-V_{\rm sys}$ maps aim at $K_{p}$ and -$K_{p}$, respectively. Orange dash dotted and white dotted pointers aim respectively at the minimum value of the maps between RV [-40, +40] $\mathrm{km\,s}^{-1}$,  $K_{p}$ [100,200] $\mathrm{km\,s}^{-1}$and between RV [-40, +40] $\mathrm{km\,s}^{-1}$,  $K_{p}$ [-200,-100] $\mathrm{km\,s}^{-1}$. The red and navy lines in the last column are the 1D projection of the $K_{p}-V_{\rm sys}$ maps evaluated at the minima found around the theoretical and injected $K_{p}$.

Figure 11: CCF results for V $\scriptstyle\rm I$  and Cr $\scriptstyle\rm I$  as in Fig. 10

Figure 12: CCF results for HCN and CO₂ as in Fig. 10. At the end of the transit (phases $\simeq$ 0.02) the noise in the CCF maps is due to the SNR drops for these exposures (see Fig. 1).

5 Lithium in the stellar spectrum

Figure 13: Master-out stellar spectrum in proximity of the lithium doublet at 6,707.76 Å.

When investigating planetary atmospheric signals, we also explored the lithium line at 6,707.76 Å (Sect. 4.3). While we did not find it in the planetary atmosphere, we noticed that it is clearly present in the master-out stellar spectrum (Fig. 13), which is not often the case for stars of this age because of its rapid depletion (e.g., Herbig, 1965). Different age estimates for this system in the literature are in the range 1.9-2.6 Gyr, all estimated based on evolutionary tracks (Sato et al., 2005; Torres et al., 2008; Carter et al., 2009; Bonomo et al., 2017; Ment et al., 2018), with just one single exception (1.2 Gyr, Southworth, 2010). However there is also the possibility of Li enrichment by, e.g., the engulfment or accretion of a close-orbiting substellar companion (e.g., Soares-Furtado et al., 2021). Interestingly, Li et al. (2008) already suggested that the observed high metallicity of HD 149026 may be confined to its surface layer as a consequence of pollution by the accretion of a gas giant or a of a population of smaller-mass rocky planets.

Since lithium is often a good age estimator (e.g., Herbig, 1965; Jeffries et al., 2023), we investigated the possibility that its presence in the stellar spectrum could be at odds with the estimated age of the parent star. We estimated the equivalent width of the lithium 6,707.76 Å line to be $36\pm 2$ mÅ, after deblending it with the close Fe $\scriptstyle\rm I$  line (e.g., Jeffries et al., 2023). We then inserted this value, together with the stellar temperature of Table 2, in the publicly available EAGLES code (Jeffries et al., 2023), which can give an independent stellar age estimate based on these two values. Unfortunately the results are totally unconstrained, likely because HD 149026 is too hot to use lithium as indicator of stellar age. One other option to test the engulfment hypothesis would be to look at the [6]Li/[7]Li isotopic ratio (e.g., Biazzo et al., 2022; Cuntz et al., 2000), but unfortunately our master-out stellar spectrum does not have sufficient SNR.

6 Discussion and conclusion

In this work, we analysed VIS and nIR high-resolution spectroscopy observations of one transit of the exoplanet HD 149026b. After analyzing the Rossiter-McLaughlin effect, which allowed us to refine the projected spin-orbit angle of the planet, we performed transmission spectroscopy. We searched for atomic and molecular species in the atmosphere of the planet by using both line-by-line and CCF techniques. We could not detect any of the searched species, finding only upper limits (Table 4). Indeed, the low flux in the core of the stellar lines, like it is common for G-type stars, does not help the detection of small atmospheric signals.

It is worth noticing that the $K_{p}-V_{\rm sys}$ maps for Ti $\scriptstyle\rm I$, Fe $\scriptstyle\rm I$, and to a lesser extent Cr $\scriptstyle\rm I$, may hint to a statistically significant (SNR=4.7, 4.2 and 3.3 respectively) signal at $\simeq+(10-27)$ $\mathrm{km\,s}^{-1}$relative to the expected position in radial velocity, consistent with the planetary $K_{p}$ value (Fig. 10 and 11). In the literature there has been some debate regarding a possible eccentricity of the HD 149026b orbit, which could indeed bring such a velocity shift for a signal from the planetary atmosphere. Wang & Ford (2011) found an eccentricity $e\simeq 0.19\pm 0.07$. Harrington et al. (2007) detected an earlier occurrence of the secondary eclipse using Spitzer Space Telescope data, approximately at -3 minutes the expected timing. Similarly, analyzing data obtained three years after, Knutson et al. (2009) predicted a secondary eclipse occurring approximately -21 minutes earlier than expected, possibly suggesting an eccentric orbit ($e\cos\omega\simeq-0.0079$ where $\omega$ is the argument of pericenter). However, the most recent analyses of both RVs and light curves of HD 149026b seem to be consistent with no eccentricity (Bonomo et al., 2017; Bean et al., 2023). To further check if even a small eccentricity could cause this shift in any atmospheric signal, we verified that the values of $e$ and $\omega$ parameters from Ment et al. (2018), when varying them within errorbars, can not be the source of a shift this large. By calculating the radial velocity (RV) curves with these parameters, we find curves that have a maximum RV of +6 $\mathrm{km\,s}^{-1}$at mid-transit (phase = 0). When varying the parameters of $e$ and $\omega$ within 3 $\sigma$, the most extreme curve shows an RV of +15.5 $\mathrm{km\,s}^{-1}$at phase = 0, which is still different from the observed signal at +20 $\mathrm{km\,s}^{-1}$. We thus tend to exclude orbital eccentricity as the possible cause of these signals. An alternative, fascinating origin for this highly redshifted ($\simeq+20$ $\mathrm{km\,s}^{-1}$) signal is from material falling on the star. Li et al. (2008) already suggested that this star may be highly metallic as a consequence of pollution by the accretion of planets. We investigated this hypothesis by studying the stellar lithium, but we found inconclusive results. Since we have only one transit observation and cannot confirm it independently, we thus tend to attribute this signal to a spurious fluctuation. We however suggest that this should be further investigated with multiple high-resolution transit observations.

Ishizuka et al. (2021) reported a $\simeq 4.4\,\sigma$ detection of Ti $\scriptstyle\rm I$, together a tentative, $\simeq 3\,\sigma$ detection of Fe $\scriptstyle\rm I$. They assumed a circular orbit. Contrary to Ishizuka et al. (2021), our analysis did not detect either Ti $\scriptstyle\rm I$ nor Fe $\scriptstyle\rm I$ in the data in the planetary restframe. By performing injection-retrieval tests, we showed that our dataset is sensitive to the model of Ti $\scriptstyle\rm I$  injected (SNR=4.1), and partially to the Fe $\scriptstyle\rm I$  one (SNR=3.3). Our analysis of the HARPS-N dataset thus excludes the presence of Ti $\scriptstyle\rm I$  in the atmosphere of the planet, in contrast with Ishizuka et al. (2021).

The actual absence of Ti $\scriptstyle\rm I$ in the planetary atmosphere could be in line with the possible presence of the so-called “titanium cold trap”. The cold trap occurs when, as in the case of HD 149026b, the equilibrium temperature is below the threshold for observing titanium in the upper atmosphere. According to Hoeijmakers et al. (2024), in planets with temperatures below $\simeq$ 2200 K, Ti $\scriptstyle\rm I$ condenses on the night side. This condensation could result in its absence throughout the planetary atmosphere, rendering it undetectable in the upper atmospheric layers where transmission spectroscopy is most effective to probe. The cold-trap process can cause titanium to remain trapped on the planet’s night side due to inefficient advection. Alternatively, if titanium is circulated back into the day side, it may reside in high-pressure, low altitude layers where it is hardly detectable. In conclusion, the low equilibrium temperature and absence of Ti $\scriptstyle\rm I$ in the planet’s atmosphere support the hypothesis of a titanium cold trap. Alternatively, the planet could possess a very dense atmosphere with extensive cloud coverage, making the atmosphere optically thick altogether at optical and near-infrared wavelengths.

Appendix A Rossiter-McLaughlin effect fitting

Table 5: set of priors used for the RM analysis.

Uniform Parameters	min	max
Systemic velocity(a) $\gamma\,\,$[$\mathrm{km\,s}^{-1}$]	-18.5	-17.5
Gaussian Parameters	$\mu$	$\sigma$
Radial velocity(b) $K_{\star}\,\,$[$\mathrm{km\,s}^{-1}$]	0.0379	0.005
Inclination(b) $i\,\,$[deg]	84.50	0.6
Stellar rotational velocity(b) $v\sin{i}\,\,$[$\mathrm{km\,s}^{-1}$]	6.00	0.50
$R_{p}/R_{\star}$ (c)	0.050565	0.0050565
Projected obliquity(d) $\lambda\,\,$[deg]	12	7
$\delta T_{0}\,\,[Phase]$ (a)	0	0.1

(a) This work; (b) Bonomo et al. (2017); (c) Spinelli et al. (2023); (d) Albrecht et al. (2012).

The parameters are distributed with a uniform and a normal distribution. The minimum and maximum values correspond to the boundaries of the uniform distribution, while $\mu$ and $\sigma$ correspond to the center and the standard deviation of the normal distribution.

The RVs were extracted by the HARPS-N DRS using a G2 mask with a CCF width of 20 $\mathrm{km\,s}^{-1}$. We used only the RVs extracted from spectra with SNR $>$ 40 (identified as red dots in Fig. 1). CaRM employs an MCMC algorithm through emcee, a python module developed by Foreman-Mackey et al. (2013). Table 5 provides all the priors used in our analysis. The guess parameter for the systemic radial velocity used, $\gamma=-18.035\,\rm{km\,s^{-1}}$, was taken from an approximate RVs data estimation. The other guess parameters we used are: the ratio between the semi-major axis and the stellar radius $a/R_{\star}=5.98$ (Spinelli et al. (2023)), the width of non-rotating star based on the HARPS-N resolving power (R $\simeq 115,000$) $\beta_{0}=2.6\,$ km s^-1, the width of the best Gaussian fit to the stellar CCFs $\sigma_{0}=4.37\,$ km s^-1 determined from a mean out-of-transit (master-out) CCF, the macro-turbolence amplitude $\zeta_{t}=4.53\,$ km s^-1 calculated from the Doyle et al. (2014) calibration valid for our stellar parameters, and the logarithmic jitter amplitude $\log\sigma_{w}=-11.5$.

Within CaRM, we specifically employed the model ARoME, that requires a set of limb-darkening coefficients. CaRM provides them using LDTk (Parviainen & Aigrain 2015), a python package that automates and models the stellar limb-darkening profiles and coefficients employing spectra from the PHOENIX database (Husser et al. 2013). In our analysis a quadratic law was employed, yielding the following coefficients: $ldc1=0.593$, $ldc2=0.123$. Following Cristo et al. (2022), we used CaRM with 50 chains, 1,500 steps of burn-in and 3,000 for the production. The RVs and their best fit are shown in Fig. 2. The fit matches well the data, except for the initial values which show a trend that is likely caused by the telescope guiding issue detailed in Sect. 2. The results obtained from our analysis are reported in Table 2, with the posterior distributions shown in Fig. 14.

Figure 14: Posterior distributions of the radial velocities analysis obtained with CaRM.