© The Authors 2023

1 Introduction

The characterization of exoplanets is among the fastest-growing and most exciting fields in modern-day astronomy. The information encoded in the exoplanet atmosphere can provide critical insights into myriad atmospheric processes, as well as the formation and evolutionary history of the planet (Madhusudhan 2019). Thanks to their relatively high equilibrium temperatures and large atmosphere scale heights, hot Jupiters (HJs) have relatively prominent atmosphere signals. This makes HJs the best targets to carry out detailed characterization of the physical and chemical properties of exoplanet atmospheres. Among all the powerful tools for the atmospheric characterization of exoplanets, highresolution spectroscopy (HRS) has become widely applied in the search for atoms and molecules, whose transmission spectra have either several strong individual lines or a few dense forests of spectral lines. It is sensitive to the change of exoplanet radial velocity, since the host star and Earth moves much slower than a planet, so that the atmospheric signals of an exoplanet are separated from the stellar and telluric signals. For the same reason, HRS can be used to characterize the planet’s atmosphere due to its sensitivity to the depth, shape, and position of the planet’s spectral lines. Compared to low-resolution spectroscopy, HRS has the advantage to probe above cloud decks where the cores of the strongest spectral lines are formed (Birkby 2018).

The first ground-based detection of CO was achieved by HRS in the hot Jupiters HD 209458b (Snellen et al. 2010). Plenty of atoms and molecules have been detected in the atmospheres of many HJs including Na, K, Li, Ca, Cr, Mg, He, Ti, Fe, Fe+, and Ca+ as well as H₂O, and so on. For example, CH₄ and Na are detected in the atmospheres of HD 209458b (Giacobbe et al. 2021; Charbonneau et al. 2002), Cr in WASP-189b (Prinoth et al. 2022), K in WASP-121b (Merritt et al. 2021), Sc and Ti+ in HD 189733b (Lira-Barria et al. 2022), Ti and TiO in HD 149026b (Ishizuka et al. 2021). All these species are detected by the HRS method.

The planet WASP-85Ab (Brown et al. 2014) was identified in both the Super-WASP and WASP-South photometry independently and confirmed by spectroscopic follow-up using the SOPHIE spectrograph (Perruchot et al. 2008) mounted on the 1.93-m telescope of the Observatoire de Haute Provence. It is a classic HJ with a mass of 1.265 ± 0.065 M_Jup, a radius of 1.240 ± 0.030 R_Jup and an equilibrium temperature, T_eq, of 1452 K with an orbital period of ∼2.65 d (Močnik et al. 2016). Thus, it is an interesting target for atmospheric characterisation given its moderately high equilibrium temperature and the low density. Up to now, no such study has been performed yet. WASP-85Ab orbits a G5V star in a binary star system. The host star WASP-85A is active with starspots detected by Močnik et al. (2016) from their analysis on the K2 short-cadence data of WASP-85A.

This paper is organized as follows. We present the details of three transit observations of WASP-85Ab using the Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO) in Sect. 2. We measure the orbital architecture of the WASP-85A system analyzed based on the Rossiter-McLaughlin (RM) effects in Sect. 3. The high-resolution transmission spectra are presented in Sect. 4, followed by our cross-correlation analysis of the spectra. The discussion and conclusions are presented in Sect. 5.

2 Observations and data reduction

Three transits of WASP-85A were observed in February and March 2021 with ESPRESSO (Pepe et al. 2021) under the ESO programs 0106.D-0853 (PI: H. M. Cegla). They aimed to detect and characterize stellar differential rotation, and convection-induced variations and to determine spin-orbit alignment, which can help us to better understand the properties of WASP-85A (Cegla et al., in prep.). In the present study, our main aim is to reveal the composition of WASP-85Ab’s atmosphere.

ESREESSO is a fibre-fed ultra-stable échelle high-resolution spectrograph, mounted at the 8.2 m Very Large Telescope at European Southern Observatory in Cerro Paranal, Chile (Pepe et al. 2021). The observations were taken in the high resolution 1-UT mode and data were read out in 2 × 1 binning and slow readout mode, which yields a spectral resolving power R of ∼140 000 and a wavelength range of 380∼788 nm. During the observations, fiber A was pointed to the target star, while fiber B was pointed to the sky for simultaneous monitoring of the sky emission.

The three transits were observed on the following dates: 2021 Feb. 13 (Night 1, the first transit, hereafter T1), Feb. 21 (Night 2, the second transit, T2), and Mar. 17 (Night 3, the third transit, T3), each covering an entire transit plus 1–2 h out-of-transit baselines. The individual exposure time for each transit is 400 s. In T1, a total of 47 spectra were obtained with 21 in transit and 26 out of transit covering the orbital phase Φ from −0.062 to +0.031. In T2, 41 spectra were acquired including 22 in-transit and 19 out-of-transit spectra with Φ from −0.043 to +0.037. In T3, the number of exposures collected in total is 39 with Φ ranging from −0.036 to +0.044. The numbers of the in- and out-of-transit spectra are 21 and 18, respectively. Details of the three observations are summarized in Table 1.

WASP-85A has a K0-dwarf stellar companion namely WASP-85B with a V-mag difference of 0.7 mag and an angular distance of −1.5″ (Brown et al. 2014), larger than the ESPRESSO fiber size of 1″. Assuming a 2D Gaussian PSF with seeing sizes of 1″ and 0.85″, we determined that WASP-85B’s flux contamination to 85A is less than 0.44 and 0.14%. This suggests that at least for the T2 and T3 data, which were mostly taken with seeing better than 0.85″, should be mostly clean from WASB-85B’s contamination.

The temporal variations of the airmass, S/N ratio, seeing and S-index during the three observing nights are shown in Fig. 1. The stellar activity level of WASP-85A as represented by the Sindex is the lowest in T3 and the highest in T2, suggesting that the data taken in T3 should be the least affected by stellar activity while those in T2 are affected most. It is also noteworthy that the median and variance of the seeing size in T1 are larger than those in T2 and T3, with quite a number exceeding 1″. Given that the fiber size of ESPRESSO is 1″ and the distance between the WASP-85A and B is ∼1.5″ (Brown et al. 2014), at least some of the spectra taken in T1 may bear both slit-loss and light contamination from the stellar companion, which will lead to lower S/N and weaker planet signals. On the other hand, most of the T2 & T3 data points should have only minor slit loss and should not be affected by the companion’s contamination. Combining the above discussion, and the fact that planet signal is the strongest or only observable in T3, we use the T3 data as the main data set in this work.

For data analysis, we used the one-dimensional (1D) sky-corrected spectra processed by the ESPRESSO reduction pipeline (version 2.3.4, downloaded from the ESO advanced science archive). The raw spectral images were corrected for bias, dark, flat, and bad pixels before being extracted into 1D spectra.

Fig. 1 Temporal variations of the airmass, S/N@550nm, seeing and S-index of the three transit observations from top to bottom, respectively. The green dashed lines indicate the first and fourth contacts of the transits.

3 Data analysis

3.1 Telluric correction

The telluric correction was first conducted by subtracting the fiber B on-sky spectra from the fiber A on-target spectra for telluric emission features, which had been performed by the ESPRESSO reduction pipeline. We then corrected the telluric absorption imprinted on the obtained spectra using the ESO software Molecfit version 1.5.7 (Smette et al. 2015), following Allart et al. (2017).

Molecfit is based on synthetic modeling of the Earth’s atmospheric transmission with a line-by-line radiative transfer model, which has been widely used recently in high-resolution transmission spectroscopic studies. Allart et al. (2017) removed telluric features in the spectrum using Molecfit in order to search for water vapor of HD 189733b; Hoeijmakers et al. (2019) detected absorption of Na I, Cr II, Sc II, and Y II in the atmosphere of an ultra-hot Jupiter, KELT-9b, after removing telluric contamination; Tabernero et al. (2021) modeled the telluric transmission spectrum using Molecfit to find the absorption of Li I, Na I, Mg I, Ca II, K I, and Fe I.

Figure A.1 presents comparisons of a sample original spectrum and its corresponding spectrum with telluric correction applied around the Na D1 & D2-line doublet region (top panel) and Hα spectral line region (bottom panel) on the night of 2018 October 31, illustrating that the telluric features have been corrected quite well. The specific parameters used in this work for telluric correction are listed in Table 2. We note that the obtained ESPRESSO spectra are given in the solar system barycentric rest frame, while Molecfit requires terrestrial rest frame spectra as input. Therefore, we shifted the observed spectra to the terrestrial rest frame considering the barycentric Earth radial velocity (BERV) for a better correction of telluric absorption.

3.2 Rossiter-McLaughlin analysis

When a planet transits its rotating host star, the portion of the stellar disk blocked by the planet varies with time, resulting in anomalous radial-velocity variation that lay over the Doppler reflex motion. This effect was first pointed out by Rossiter (1924) and McLaughlin (1924), and is known as the RM effect (Triaud 2018). Queloz et al. (2000) reported the first discovery of the RM effect on an exoplanet, by analyzing the anomaly in the radial velocity curve after removing the influence of Doppler reflex motion. From this anomaly, they determined υ sin i_⋆ as well as the angle between the orbital plane and the apparent equatorial plane, namely the spin-orbit angle, λ.

Two main methods are widely used to measure the radial velocities (RVs), one is based on least-squares fitting by carrying out the template matching of the observed spectra (Butler et al. 1996); the other utilizes the cross-correlations function (CCF) of observation spectra with the binary mask (Baranne et al. 1996; Pepe et al. 2002). In this work, the RVs and their uncertainties determined using the CCF method by the ESPRESSO pipeline are adopted.

To analyze and model the RM effect based on the RVs, we used a Markov chain Monte Carlo (MCMC) algorithm, available in emcee (Foreman-Mackey et al. 2013), and we employed the RM effect model presented in Ohta et al. (2005) implemented in PyAstronomy (Czesla et al. 2019) as modelSuite.RmcL, which performs the analytical fit for the RM effect of a circular orbit. The Rmcl model incorporates nine parameters, including the planet’s orbital period, P, the time of transit midpoint, T_c, the inclination of planetary orbital plane, i, the inclination of the stellar rotation axis, i_⋆, the ratio of planetary and stellar radius, R_p/R_⋆, the semi-major axis, a, the linear limb-darkening coefficient, ∈, the stellar angular rotation velocity, Ω, and the sky-projected angle, λ (angle between the stellar rotation axis and the normal of planetary orbital plane). Before modeling the RM effect, we estimated and removed the corrected noise in the observed RVs with Gaussian process (GP, Foreman-Mackey et al. 2017) and corrected for the Doppler reflex motion using Kepler’s Laws. The rest RVs are used as input for the RM effect model. We note that there are four RV points in T3 deviating from the RV trend that are excluded from the RM analysis.

Assuming a circular planetary orbit, the RV curve with the baseline trend removed can be expressed with the following formula: $$V_{★}=K_{★}\sin \left(2\pi \varphi \right)+V_{\text{bary}}+V_{\text{sys}}$$(1)

where V_⋆ is the RV of the host star, K_⋆ is the RV semi-amplitude, ϕ is the orbital phase of the planet, V_bary is the barycentric velocity because of the Earth’s revolution, V_sys is the radial velocity of targeted exoplanet system. Among them, V_⋆, ϕ and V_bary are functions of time, while the values of K_⋆ and V_sys for each night are derived from our GP analysis as listed in Table 3. Here, we use the K_⋆ value shown in Table 4 and fixed it in the next steps.

We set λ, ∈, Ω, and i_⋆ as free parameters for our model, with their priors listed in Table 5. The remaining parameters including a, P, i, T_c, and R_p/R_⋆ are fixed to the specific values listed in Table 4. In order to explore a wide parameter space, we set the walkers to 200 with 50 000 steps. The first 10 000 steps were discarded as burn-in in order to probe the parameter space and get settled into the maximum density in the MCMC chain. The thus obtained parameters from the MCMC fitting are shown in Table 5, whereas the retrieved best model is compared with the data set in Fig. 2, with an rms residual of −1.877 m s⁻¹. The derived projected spin-orbit angle $\lambda =-16.155{°}_{-2.879}^{+2.916}$, implying the planetary orbital plane is roughly aligned to the stellar equator. We also derive the projected stellar rotation velocity $\upsilon \sin i_{\star }={2.98}_{-0.043}^{+0.050}$ km s⁻¹, the linear limb dark coefficient $ϵ={0.855}_{-0.014}^{+0.014}$. These values are used as input to model the center-to-limb variation (CLV) effect and RM effect in order to eliminate the influence caused by the occultation of the host star by the planet. The derived posterior distributions for the three transits are given in Fig. A.2 in detail.

Fig. 2 RM effect fit using the GP+PyAstronomy package on the ESPRESSO RVs taken from the spectrum headers. Top panel: RV curve with Kepler motion and red noise removed, showing the RM effect of WASP-85Ab of the three transits represented by the red, green, and blue circles with error bars. The best-fit model using the PyAstronomy package is shown in the solid line. Bottom panel: residuals between the data and the model predictions.

4 Transmission spectrum analysis

The atmosphere composition of HJs plays an important role in understanding atmospheric dynamics and planetary evolution. Here we mainly search for the most abundant atoms including Na I, Mg I, Hα, Hβ Ca II H, Ca II K, and Li I. We adopt the method outlined in Wyttenbach et al. (2015) and Casasayas-Barris et al. (2019) to extract transmission spectra of individual lines. As described in Sect. 3, the observed spectra were corrected for telluric contamination. Then the spectra were normalized to their corresponding continuum level with a two-degree spline curve fit with the regions around strong absorption lines excluded. We also applied a sigma-clipping rejection algorithm on the normalized spectra and replaced the cosmic ray hits with the mean value of all the other spectra at each wavelength (Allart et al. 2017).

In order to align the stellar lines, we first shifted the spectra to the stellar rest frame, taking into account BERV, V_sys, and the stellar reflex motion induced by the planet, V_reflex. Specifically, BERV is obtained directly from the files header information of the ESPRESSO spectra, while V_sys is derived from the RM analysis (as described in Sect. 3.2) and V_reflex is calculated using the orbital parameters (in Table 4). Then the aligned out-of-transit spectra are averaged with their mean signal-to-noise ratio (S/N) as their weights to construct a master out-of-transit spectrum, which is then used to divide each individual spectrum in order to remove stellar light contribution so that the absorption signal from the planetary atmosphere may be exposed (Stangret et al. 2021). Next, we shifted each divided in-transit spectrum to the planetary rest frame using the planet radial velocity, Vp, determined from Kp and the orbit phase of a given time with the following formula: $$\Re \left(\lambda \right)={\displaystyle \sum _{\text{in}}\frac{F_{\text{in}}\left(\lambda \right)}{F_{\text{out}}}|{}_{\text{Planet\hspace{0.17em}RV\hspace{0.17em}Shift}}-1,}$$(2)

Finally, the transmission spectrum is the S/N-weighted sum of in-transit residuals.

We note that there is a wiggle pattern in the derived transmission spectrum, with an amplitude of ∼0.075% and a period of ∼30 Å, which might potentially affect the final spectral signature analysis and need to be corrected. These wiggles were first reported in Allart et al. (2020) and are possibly induced by an interference pattern caused by the Coudé train optics as pointed out by Tabernero et al. (2021). To correct the wiggles, we remove a sinusoidal trend by fitting a set of cubic splines to the obtained transmission spectra. Figure 3 shows an example of such correction.

The CLV and RM effect affect stellar lines profile and may produce additional time-correlated signals in the transmission spectra. These two effects are both induced by the joint effect of planet occultation and stellar rotation. We adopt the method as presented in Yan & Henning (2018), Chen et al. (2020), and Casasayas-Barris et al. (2019) to model the stellar spectra at different transit positions, which considers both the RM and CLV effects. We used the Spectroscopy Made Easy tool (SME, Valenti & Piskunov 1996) to compute the theoretical stellar spectra at 21 different limb-darkening angles (µ) using the MARCS and VALD3 line lists (Ryabchikova et al. 2015). We employed the solar abundance and local thermodynamical equilibrium (LTE) for the calculation of stellar spectra.

The stellar disc is then divided into elements of size 0.01 R_⋆ × 0.01 R_⋆, each element owning specific parameters, including υ sin i_⋆, µ, and θ, which is the angle between the normal to each element and the line of sight. The position of the planet relative to the stellar disk is calculated assuming a uniform velocity during the transit. Then, the synthetic spectrum during transit is calculated by integrating all the surface elements, which are not obscured by the planet (Yan & Henning 2018), considering the obscured elements, the proper RV shift, and the corresponding interpolated µ spectrum.

The next step is to simulate the CLV and RM effects and correct them. We constructed the master out-of-transit spectrum of the synthetic spectrum and divided each synthetic spectrum by it and then we shifted each in-transit residual of the synthetic spectrum to the planetary rest frame and calculated the sum of all in-transit residuals which only contains the CLV and RM effects. These simulated spectra were re-scaled to match the observed spectra and are subtracted to eliminate the influence of CLV and RM effects. Figure 4 shows the CLV and RM effects in the residual spectra around the Na D1 and D2 lines.

Fig. 3 Example showing the wiggle pattern in one transmission spectrum taken in T2 in blue, the fitted model in red, and the wiggle-corrected spectrum in green.

Fig. 4 Simulated CLV and RM effects overplotted in red on the transmission spectrum. The transmission spectrum is shown in gray, and the spectrum binned by 0.1 Å is shown in black.

4.1 Strong atomic lines

We went on to perform a visual inspection on the phase-resolved transmission spectra around each known strong feature, which allows us to confirm or deny the presence of the atoms of interest Tabernero et al. (2021). An integrated S/N-weighted transmission spectrum was obtained and fitted with a Gaussian function to determine the center, depth, and width of each line. Figure 5 depicts the 2D phase-resolved transmission spectra and the integrated spectrum around Hα, showing a clear detection of H absorption from the WASP-85Ab’s atmosphere. Figure A.3 shows the same plots for other atomic lines, including the Na I, Mg I, Li I, Hα, Hβ, Ca II H, and K lines. The Gaussian fit results are listed in Table 6. In addition, Table 6 also lists the observed Doppler shift of the line center (V_wind), which traces the planetary winds towards the observer from the morning hemisphere to the evening hemisphere and effective wavelength-dependent planet radius, R_λ, which can be derived with $\sqrt{1+h/\delta }{R}_{\text{p}}$, where h is the line depth corresponding to different atoms and δ is the transit depth of the planet.

As shown in Table 6 and Fig. 5, Hα absorption is marginally visible in the 2D spectrum (the left panel) and is clearly detected in the integrated spectrum (the right panel), with a depth of ∼0.625 ± 0.06%, a full width at half maximum (FWHM) of ∼23.806 ± 2.50 km s⁻¹, a R_λ of 1.028 ± 0.003 R_p. Clear absorption is detected for the Ca II doublet in the integrated transmission spectra (cf. two bottom panels of Fig. A.3). Values for the FWHM and R_λ are similar to those of Hα. However, it is not visible in the H line 2D spectrum and is only marginally noticeable for the K line. The derived V_wind for the H & K lines are approximately −4.805 and 5.974 km s⁻¹, respectively, while the obvious inconsistency may arise from the T1 and T2 spectra that suffer large uncertainties due to stellar activity and/or companion’s flux contamination. As shown in the top two rows of Fig A.6, while the H & K lines in T3 are both close to the expected line centers, they are obviously deviated from the line centers (either blue-shifted or red-shifted). There seems to be 3σ detection for Li I and no clear detection for the other species including Na I, Mg I, and Hβ. For those with poor Gaussian fits, namely, Mg I and Hβ, the upper limits of the absorption depths are estimated and listed in Table 6. As shown in Figs. A.5 and A.6, the T3 spectrum provides the most significant signal for all the species, while the T1 spectrum provides the least information. This is reasonable according to the discussion on the observing conditions in Sect. 2.

We note that the Na I D2 absorption line in T1 is abnormal, bearing huge error bars (as shown in Fig. A.5). Through a visual inspection of every 2D raw images, we find that the Na I D2 line in nearly half of the T1 spectra is contaminated by two bright spikes, while the corresponding 1D spectra without sky subtraction show strong emission at 589 nm (examples are given in Fig. A.7). These spikes appear simultaneously in both the target and sky fibers and some of them are saturated, leading to close-to-zero counts and thus large uncertainties after sky subtraction. Given the fact that these bright spikes only present around Na I D2 at 589 nm, it is quite likely that they arise from projected laser signal leaked from the VLT Laser Guide star.

In order to further confirm the existence of these signals, we performed empirical Monte Carlo (EMC) simulations (Redfield et al. 2008) to estimate the effect of systematics, a technique that has been widely used in previous studies (Wyttenbach et al. 2015; Casasayas-Barris et al. 2019; Seidel et al. 2019; Allart et al. 2020). The concept of this method is to artificially create new data set by randomizing original data to verify whether the investigated signals still exist. We explored three scenarios covering in-in, in-out, and out-out for each strong line, with 10000 iterations for the ∼1.5 Å passband. The results are shown in Fig. A.8. The in-out distributions of Li I, Hα, Ca II H, and Ca II K exhibit excess absorption, while the absorption depths in the in-in and out-out scenarios are centered at zero. The EMC simulations suggest that the signals of Li I, Hα, Ca II H, and K are likely to have been created by the transits and may originate from the planet itself. In addition, we also masked the line cores with a width of 0.1 Å, where the S/N values are low and the RM+CLV residuals may still exist. These results are shown in Fig. A.4 and the fitted result is listed in Table 6.

Fig. 5 Phase-resolved 2D transmission spectra based on the combination of three nights around the Hα line without and with RM+CLV correction applied in the left panel. The horizontal black-dashed lines indicate the beginning and end of the transit. The inclined black-dashed line presents the expected trace of signal from the exoplanet atmosphere. The in-transit part of this line is not plotted so as not to interfere with the visibility of the trace signal. The middle panel shows the 2D transmission spectra in the planetary rest frame (PRF) assuming Kp = 159.76 km s−1, and the vertical black-dashed lines indicate the position of expected signal from the exoplanet atmosphere in the PRF. The right panel shows the combined integrated transmission spectrum of Hα. The grey line represents the transmission spectra in the PRF that has been corrected for the CLV+RM effects, the black dotted line is the binned version with a bin size 0.1 Å, while the best Gaussian fit of the Hα line is shown in red. The dashed green vertical line represents the static position of Hα at vacuum wavelength.

4.2 Cross-correlation function analysis

In addition to the direct measurements of strong lines, we employed the CCF analysis to explore the presence of species or enhance their detection. The CCF method is to extract multiple spectral signals from a certain species by cross-correlating the template spectrum with only this species included with the observed transmission spectra. This technique is widely used for the high-resolution transmission spectroscopic study of exo-planets, which can fully utilize all the spectral signals from multiple lines of a certain species to maximize the detectability of that species. In this work, the species to be explored include Na I, K I, Ca I, Cr I, Fe I, Li I, Mg I, Ti I, V I, and Y I as well as the molecules TiO and VO, and the spectrum template is computed by petitRADTRANS code (Mollière et al. 2019). As input parameters, we assumed the WASP-85Ab radius is 1.24 R_J (Močnik et al. 2016), the planetary mass is 1.265 M_J, and, thus, log ɡ_p = 3.329 cgs. We applied an isothermal temperature of 1500 K and assumed a solar abundance for the determination of the volume mixing ratio (VMR) of various species. The template spectrum generated by petitRADTRANS is convolved to match the resolution of ESPRESSO. We set the radial velocity at the range of ±100 km s⁻¹, with a step of 0.5 km s⁻¹, and calculate CCF for each residual spectrum and model spectrum. The cross-correlation coefficients are calculated by the following formula: $$c\left(\upsilon ,t\right)=\frac{{\text{Σ}}_i^N{x}_i\left(t\right)T_i\left(\upsilon \right)}{{\text{Σ}}_i^N{T}_i\left(\upsilon \right)},$$(3)

where T_i(v) is the template shifted to a radial velocity of v, x_i(t) is the transmission spectrum at time t, c(v, t) is a 2D matrix, which is the function of t and v. If the investigated species does exist, the signal will appear at the position of the estimated orbital velocity, K_p, and system velocity, V_sys (Snellen et al. 2010; Hoeijmakers et al. 2019).

The deformation of the stellar line profile caused by CLV and RM effects can also mask the planetary absorption signal in the cross-correlation maps that contains RM+CLV and planetary atmosphere signal together. It is necessary to correct these effects by modeling stellar lines at each phase. As described in Sect. 3, we used SME and the VALD3 line list to compute theoretical stellar spectra that take into account the CLV and RM effects. The simulated stellar spectra are then cross-correlated with template spectra with a certain species, and the resulting 2D CCF map can be used as a proxy for the influence of the RM+CLV effects and should be subtracted from the CCF map of the observed transmission spectra. The CCF map of the CLV+RM model CCF is shown in Fig. 6.

The obtained 2D CCF map with CLV+RM effects uncor-rected and corrected for Li I, Na I, Ca I, K I, Cr I, and Fe I are shown in the first and second panels of Fig. 7, where the inclined white dashes line indicates the expected track of planet signal. The third panels show the K_p − ΔV_sys map, in which the planet signal, if it exists, should appear in the intersection of white dashed lines; namely, with a velocity close to V_sys and an RV semi-amplitude K_p = 159.76 ± 4.09 km s⁻¹. The fourth panels show the S/N of the corresponding species at the expected K_p. From Fig. 7, only Li I is clearly detected with a maximum S/N of ∼4.5, where the peak position corresponds to K_p of ∼134.8 ± 88.43 km s⁻¹ and ΔV_sys ∼ −4.6 ± 11.29 km s⁻¹. Some remaining features in the 2D maps may be associated with stellar activity (Stangret et al. 2021) and propagated noise. The CCF maps and S/N plots for Mg I, Ti I, V I, Y I and TiO, and VO are shown in Fig. 8, with no detection reported.

Fig. 6 CCF map of model spectra generated by SME for Na I on 2021 Feb. 21, with only CLV and RM effects included, which is used to eliminate the influence of the variation of stellar line profile during transit.

5 Discussion and conclusion

We observed three transits of the hot Jupiter WASP-85Ab using the ultra-stable high-resolution spectrograph ESPRESSO. A total of 127 spectra were obtained, with 64 spectra were taken during transits, while the rest 63 spectra were taken out of transits. Telluric contamination was corrected for each spectrum, from which the master out-of-transit spectra were created and the transmission spectra were generated for the three transits.

For the observed RV of the three transits, we performed a joint RM analysis, which can provide a good estimate of the origin of planets’ angular momentum (Ohta et al. 2005). We retrieved the projected obliquity value of $\lambda =-16.155{°}_{-2.879}^{+2.916}$, suggesting that the hot Jupiter WASP-85Ab’s orbit is almost aligned with its host star. The spectral signals induced by the CLV and RM effects were modeled and removed from the obtained transmission spectra and CCF map. The residual transmission spectra were used for the exploration of atomic and molecular lines that originate from WASP-85Ab’s atmosphere, via direct inspection for strong lines or the CCF method for the species with multiple lines. The species we explored include H I,Li I,Na I, Ca II, K I, Mg I, Fe I, Cr I, TiO, and VO.

Among them, a direct inspection revealed a ∼10σ absorption signature at Hα, but it is not possible to confirm its origin as either stellar, planetary, or otherwise. In the meanwhile, the Ca II H & K lines are also tentatively detected, with quite similar FWHM and consistent R_λ. The potential spectral signals detected from H I and Ca II may arise from the planetary atmosphere, but we cannot yet exclude the possibility of being induced by stellar activity. Another species that may be detected is Li I, which has a 4.5σ significance at the estimated K_p velocity in the K_p − ΔV_sys map. We note there is an offset of 21 ± 4.3 km s⁻¹ between the estimated K_p and the retrieved K_p. This discrepancy may be due to the large uncertainty of K_p, given that the Li CCF signal is quite extended in the K_p direction.

We note that there are still some structures visible in the 2D CCF maps, which may be the RM+CLV effect residuals that have not been removed completely or are caused by the low S/Ns in the stellar cores. We thus mask out the line cores with a width of 0.1Å for the tentatively detected species to check whether there are any signals left. As shown in Fig. A.4, the absorption features of Ha, Ca II H&K, and Li I are still visible, although they are slightly shallower. Therefore, the tentative detection of these three elements are unlikely to be caused by the RM+CLV effects or the low S/Ns of the line cores.

The detection should not be affected by stellar activity, because there is no lithium line in the stellar spectrum. For the same reason, the continuum around the lithium lines possesses very high S/Ns, resulting in a relatively strong lithium signal. Lithium was first reported in the atmosphere of WASP-127b by Chen et al. (2018) with low-R spectroscopy but was not confirmed with high resolution in Allart et al. (2020). Then Borsa et al. (2021) and Tabernero et al. (2021) presented the first detection of Li in the atmosphere of WASP-121b and WASP-76b respectively at high resolution. The detection of Li in an exoplanet atmosphere is not unexpected, as substellar objects with masses below ∼55 M_Jup do not deplete this element during their lifetime (Chabrier & Baraffe 2000; Baraffe et al. 2015). Such a detection, if confirmed, is an important step, which can help improve the understanding of planet formation history and lithium depletion in planet-hosting stars (e.g. Bouvier 2008; Chen et al. 2018).

Following the formalism of Wyttenbach et al. (2015) and Allart et al. (2017), we estimated the 5σ upper limits of the excess absorption of the non-detected atoms and molecules, and listed them in Table 7. Only lines with amplitudes higher than 10 ppm are included in the template spectrum, and the corresponding number of lines employed are listed in Table 7 as well. Figure 9 is an example showing how an upper limit is obtained by comparing the theoretical CCF using template spectrum with the measured CCF with the observed transmission spectra.

It is somehow expected that the spectral features seen in the transmission spectra are very few and in general weak, given the relatively low T_eq and small R_p and, thus, the small transmission spectroscopy metric (Kempton et al. 2018); in addition, we assume that the star is active with recurring starspot reported which may induce relatively large noise in the derived spectra and 2D maps, although our measured S-index activity indicator of WASP-85A for the three nights scarcely varies. in addition, the lack and weakness of atomic and molecular features can be due to several technical factors. The first reason is the low S/N of the observed stellar spectrum near the cores of the strong lines, where the photon counts are close to zero. Another factor is the overlap of CLV+RM effect and planet radial-velocity in cross-correlation residuals maps, making it difficult to disentangle the signal of an exoplanet atmosphere, as noted for example by Casasayas-Barris et al. (2022). More precise measurements of spin-orbit angle should in principle be able to better constrain the planet-occulating position, which will help to eliminate the influence of the variation of stellar line profile during transit.

We also note that the three-night data set do not agree with each other very well. it seems the last night exhibits the most significant spectral signals, which should be related to the combined effect considering the lowest stellar activity indicated by S-index and lower impacts of WASP-85B showed by seeing value. More high-precision observations are necessary in order to further constrain and understand the atmosphere of WASP-85Ab.

Fig. 7 CCF results derived by cross-correlating the template spectra produced by petitRADTRANS with the observed transmission spectra. First panels: 2D CCF maps of Li I, Na I, Ca I, K I, Cr I, and Fe I with CLV+RM effects uncorrected. The white dotted lines mark the beginning and ending position of the transit, and the inclined white lines indicate the expected trace of signal from the planet. Second panels: same as the first panels but with CLV+RM effects corrected. Third panels: Kp-ΔVsys maps in the range of −50 ∼ 350 km s−1. The signal is expected to appear around the intersection of two white dotted lines. Fourth panels: S/N plots at the expected Kp position.

Fig. 8 Same as Fig 7: but for Mg I, Ti I, V I, Y I TiO, and VO.

Fig. 9 Comparison between the measured CCF and theoretical CCF calculated by the template spectrum and observed data.

Acknowledgements

We thank the anonymous reviewer for the constructive comments, thank Lauren Doyle for useful conversations and thank Jens Hoeijmakers for many useful suggestions. This research is the National Natural Science Foundation of China grants No. 11988101, 42075123, 42005098, 62127901, supported by the National Key R&D Program of China No. 2019YFA0405102, the Strategic Priority Research Program of Chinese Academy of Sciences, Grant No. XDA15072113, the China Manned Space Project with NO. CMS-CSST-2021-B12. M.Z., Y.Q.S., Q.L.O.Y. are supported by the Chinese Academy of Sciences (CAS), through a grant to the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. H.M.C. acknowledges support from a UKRI Future Leaders Fellowship (MR/S035214/1).

Appendix A Additional figures

Fig. A.1 Effect of a telluric correction using Molecfit software in the observed spectrum of WASP-85Ab on the night of 2021 Feb 13 (T1) for single exposure. The observed spectrum is shown in blue line and the spectrum after telluric correction is shown in orange line. The dashed green vertical line represent the static position of atom spectral line at vacuum wavelength. Top panel: Telluric correction of spectral line around Na D1 & D2 lines. Bottom panel: Telluric correction of spectral line around the Hα lines.

Fig. A.2 Posterior distribution of the best-fit RM model using PyAstronomy and GP.

Fig. A.3 Same as Fig. 5,but for Na I, Mg I, Hβ, Li I, Ca II H, and Ca II K.

Fig. A.4 Masked phase-resolved 2D transmission spectra based on the combination of three nights for Hα, Li I, Ca II H, and Ca II K lines with RM+CLV correction applied in the left panel. The horizontal black-dashed lines indicate the beginning and end of the transit. The inclined black-dashed line presents the expected trace of signal from the exoplanet atmosphere. The middle panels show the 2D transmission spectra in the planetary rest frame (PRF) assuming Kp = 159.76 km s−1. In the right panels, the integrated transmission spectra in PRF with the CLV+RM effects corrected are shown in grey (original) and black (binned) and the best Gaussian fit are shown in red. The dashed green vertical line represents the static position of each line at vacuum wavelength.

Fig. A.5 Transmission spectra of Na I, Mg I, Hα, and Hβ for each night. Left (T1), middle (T2), and right (T3). The grey line represents the raw transmission spectra in the PRF which has been corrected for CLV+RM effects, The black dotted line represents the binned spectra of 0.1 Å and the Gaussian fitting is shown in red line. The dashed green vertical line represent the static position of atom individual line at vacuum wavelength.

Fig. A.6 Same as Fig A.5, but for Ca II H,Ca II K, and Li I.

Fig. A.7 Na I D2 absorption line shown in the raw data image of T1 appears to be contaminated by two bright spikes which leads to the abnormal flux around Na I D2 line (left). Plot of the abnormal counts of no sky-subtracted flux colored in blue and sky-subtracted flux shown in red (right).

Fig. A.8 Results of EMC simulation for Na I, Mg I, Li I, Hα, Hβ, Ca II H, and Ca II K combined over three nights. These distributions for three EMC scenarios are shown in different colors: red is the in-out distribution; green is the in-in distribution; and blue is the out-out distribution. The in-out distribution of Li I, Hα, Ca II H and Ca II K exhibits excess absorption.

References

All Tables

All Figures

	Fig. 1 Temporal variations of the airmass, S/N@550nm, seeing and S-index of the three transit observations from top to bottom, respectively. The green dashed lines indicate the first and fourth contacts of the transits.
In the text

Fig. 2 RM effect fit using the GP+PyAstronomy package on the ESPRESSO RVs taken from the spectrum headers. Top panel: RV curve with Kepler motion and red noise removed, showing the RM effect of WASP-85Ab of the three transits represented by the red, green, and blue circles with error bars. The best-fit model using the PyAstronomy package is shown in the solid line. Bottom panel: residuals between the data and the model predictions.

In the text

	Fig. 3 Example showing the wiggle pattern in one transmission spectrum taken in T2 in blue, the fitted model in red, and the wiggle-corrected spectrum in green.
In the text

	Fig. 4 Simulated CLV and RM effects overplotted in red on the transmission spectrum. The transmission spectrum is shown in gray, and the spectrum binned by 0.1 Å is shown in black.
In the text

Fig. 5 Phase-resolved 2D transmission spectra based on the combination of three nights around the Hα line without and with RM+CLV correction applied in the left panel. The horizontal black-dashed lines indicate the beginning and end of the transit. The inclined black-dashed line presents the expected trace of signal from the exoplanet atmosphere. The in-transit part of this line is not plotted so as not to interfere with the visibility of the trace signal. The middle panel shows the 2D transmission spectra in the planetary rest frame (PRF) assuming Kp = 159.76 km s−1, and the vertical black-dashed lines indicate the position of expected signal from the exoplanet atmosphere in the PRF. The right panel shows the combined integrated transmission spectrum of Hα. The grey line represents the transmission spectra in the PRF that has been corrected for the CLV+RM effects, the black dotted line is the binned version with a bin size 0.1 Å, while the best Gaussian fit of the Hα line is shown in red. The dashed green vertical line represents the static position of Hα at vacuum wavelength.

In the text

	Fig. 6 CCF map of model spectra generated by SME for Na I on 2021 Feb. 21, with only CLV and RM effects included, which is used to eliminate the influence of the variation of stellar line profile during transit.
In the text

Fig. 7 CCF results derived by cross-correlating the template spectra produced by petitRADTRANS with the observed transmission spectra. First panels: 2D CCF maps of Li I, Na I, Ca I, K I, Cr I, and Fe I with CLV+RM effects uncorrected. The white dotted lines mark the beginning and ending position of the transit, and the inclined white lines indicate the expected trace of signal from the planet. Second panels: same as the first panels but with CLV+RM effects corrected. Third panels: Kp-ΔVsys maps in the range of −50 ∼ 350 km s−1. The signal is expected to appear around the intersection of two white dotted lines. Fourth panels: S/N plots at the expected Kp position.

In the text

	Fig. 8 Same as Fig 7: but for Mg I, Ti I, V I, Y I TiO, and VO.
In the text

	Fig. 9 Comparison between the measured CCF and theoretical CCF calculated by the template spectrum and observed data.
In the text

Fig. A.1 Effect of a telluric correction using Molecfit software in the observed spectrum of WASP-85Ab on the night of 2021 Feb 13 (T1) for single exposure. The observed spectrum is shown in blue line and the spectrum after telluric correction is shown in orange line. The dashed green vertical line represent the static position of atom spectral line at vacuum wavelength. Top panel: Telluric correction of spectral line around Na D1 & D2 lines. Bottom panel: Telluric correction of spectral line around the Hα lines.

In the text

	Fig. A.2 Posterior distribution of the best-fit RM model using PyAstronomy and GP.
In the text

	Fig. A.3 Same as Fig. 5,but for Na I, Mg I, Hβ, Li I, Ca II H, and Ca II K.
In the text

Fig. A.4 Masked phase-resolved 2D transmission spectra based on the combination of three nights for Hα, Li I, Ca II H, and Ca II K lines with RM+CLV correction applied in the left panel. The horizontal black-dashed lines indicate the beginning and end of the transit. The inclined black-dashed line presents the expected trace of signal from the exoplanet atmosphere. The middle panels show the 2D transmission spectra in the planetary rest frame (PRF) assuming Kp = 159.76 km s−1. In the right panels, the integrated transmission spectra in PRF with the CLV+RM effects corrected are shown in grey (original) and black (binned) and the best Gaussian fit are shown in red. The dashed green vertical line represents the static position of each line at vacuum wavelength.

In the text

Fig. A.5 Transmission spectra of Na I, Mg I, Hα, and Hβ for each night. Left (T1), middle (T2), and right (T3). The grey line represents the raw transmission spectra in the PRF which has been corrected for CLV+RM effects, The black dotted line represents the binned spectra of 0.1 Å and the Gaussian fitting is shown in red line. The dashed green vertical line represent the static position of atom individual line at vacuum wavelength.

In the text

	Fig. A.6 Same as Fig A.5, but for Ca II H,Ca II K, and Li I.
In the text

	Fig. A.7 Na I D2 absorption line shown in the raw data image of T1 appears to be contaminated by two bright spikes which leads to the abnormal flux around Na I D2 line (left). Plot of the abnormal counts of no sky-subtracted flux colored in blue and sky-subtracted flux shown in red (right).
In the text

	Fig. A.8 Results of EMC simulation for Na I, Mg I, Li I, Hα, Hβ, Ca II H, and Ca II K combined over three nights. These distributions for three EMC scenarios are shown in different colors: red is the in-out distribution; green is the in-in distribution; and blue is the out-out distribution. The in-out distribution of Li I, Hα, Ca II H and Ca II K exhibits excess absorption.
In the text