11institutetext: CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Datun Road A20,
Beijing 100101, China 22institutetext: School of Astronomy and Space Science,
University of Chinese Academy of Sciences, Beijing 100049, China
33institutetext: Key Laboratory of Planetary Sciences, Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210033, PR China
44institutetext: Department of Astronomy, University of Science and Technology
of China, Hefei 230026, China 55institutetext: Physics Department, University
of Warwick, Coventry CV4 7AL, United Kingdom 66institutetext: Centre for
Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, UK
77institutetext: Departamento de Astronomía, Universidad de Chile, Camino El
Observatorio 1515, Las Condes, Santiago, Chile 88institutetext: CAS South
America Center for Astronomy, National Astronomical Observatories, Chinese
Academy of Sciences, Datun Road A20, Beijing 100101, China

# Atmospheric composition of WASP-85Ab with ESPRESSO/VLT observations

Zewen Jiang , 1122 Wei Wang E-mail<EMAIL_ADDRESS>Guo Chen 33 Fei Yan 44
Heather M. Cegla UKRI Future Leaders Fellow5566 Patricio Rojo 77 Yaqing Shi ,
1122 Qinlin Ouyang , 1122 Meng Zhai 88 Yujuan Liu 11 Fei Zhao 11 Yuqin Chen 11

(Received xxx ; accepted xxx)

Transit spectroscopy is the most frequently used technique to reveal the
atmospheric properties of exoplanets, while that at high resolution has the
advantage to resolve the small Doppler shift of spectral lines, and the trace
signal of the exoplanet atmosphere can be separately extracted. We obtain the
transmission spectra of the extrasolar planet WASP-85Ab, a hot Jupiter in a
2.655-day orbit around a G5, $V=11.2$ mag host star, observed by high-
resolution spectrograph ESPRESSO at the Very Large Telescope array for three
transits. We present an analysis of the Rossiter-McLaughlin effect on
WASP-85A, and determine a spin-orbit angle
${\lambda=-16.155^{\circ}}^{+2.916}_{-2.879}$, suggesting that the planet is
in an almost aligned orbit. Combining the transmission spectra of three
nights, we tentatively detected H$\alpha$ and Ca II absorption with
$\gtrapprox 3\sigma$ via direct visual inspection of the transmission spectra
with the Center-to-Limb variation and the Rossiter-McLaughlin effects removed,
which still remain visible after excluding the cores of these strong lines
with a 0.1Å mask. These spectral signals seems likely to origin from the
planetary atmosphere, but we can not fully exclude their stellar origins. Via
the cross-correlation analysis of a set of atoms and molecules, Li I is
marginally detected at $\sim 4\sigma$ level, suggesting that Li might be
present in the atmosphere of WASP-85Ab.

###### Key Words.:

planetary Systems – planets and satellites: individual: WASP-85Ab – planet and
satellites: atmospheres – methods: observational – techniques: spectroscopic

##  1 Introduction

Characterization of exoplanets is one of the fastest-growing and most exciting
fields in astronomy in recent years. 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, high-resolution spectroscopy (HRS) becomes widely applied to
search for atoms and molecules whose transmission spectra have either several
strong individual lines or a few dense forests of spectral lines. HRS 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 H2O, etc. For example, CH4 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\pm 0.065$ $M_{\rm Jup}$, a radius of $1.240\pm 0.030$
$R_{\rm Jup}$ and an equilibrium temperature $T_{\rm eq}$ of 1452 K with an
orbital period of $\sim$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.

Table 1: Summary of the WASP-85Ab transit observations | Date | VLT/UT | Number of spectra | Exp. time | Airmass range | Mean S/N | Program ID
---|---|---|---|---|---|---|---
| (UT Time) |  | Total | In-transit | Out-of-transit | (s) |  | (@550 nm) |
Night 1 | 2021-02-13 | UT3 | 47 | 21 | 26 | 400 | $1.1-1.7$ |  $\sim$64 | 0106.D-0853
Night 2 | 2021-02-21 | UT4 | 41 | 22 | 19 | 400 | $1.1-1.6$ |  $\sim$63 | 0106.D-0853
Night 3 | 2021-03-17 | UT2 | 39 | 21 | 18 | 400 | $1.1-1.5$ |  $\sim$61 | 0106.D-0853

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 Section 2. We measure the
orbital architecture of the WASP-85A system analyzed based on the Rossiter-
McLaughlin (RM) effects in Section 3. The high-resolution transmission spectra
are presented in Section 4, followed by our cross-correlation analysis of the
spectra. The discussion and conclusions are presented in Section 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 (H.M Cegla et al.
in prep). While in this study, we will mainly target to reveal the composition
of WASP-85Ab atmosphere.

ESREESSO is a fibre-fed ultra-stable échelle high-resolution spectrograph,
mounted at the 8.2m 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\times 1$ binning and
slow readout mode, which yields a spectral resolving power $R$ of $\sim$140
000 and a wavelength range of 380$\sim$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 to 2 hours out-of-transit baselines. 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 $\Phi$ from
$-0.062$ to $+0.031$. In T2, 41 spectra were acquired including 22 in-transit
and 19 out-of-transit spectra with $\Phi$ from $-0.043$ to $+0.037$. In T3,
the number of exposures collected in total is 39 with $\Phi$ 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 $\sim 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 determine 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 S-index 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 worthy to
note 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 $\sim
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 ratios and weaker planet signals. On the other hand, most of
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 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.

Figure 1: The temporal variations of the airmass, S/N@550 nm, 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. Table 2: Molecfit parameters used for telluric correction Parameter | Value | Description
---|---|---
ftol | $10^{-9}$ |  ${\chi}^{2}$ tolerance
xtol | $10^{-9}$ | Tolerance for the molecfit
|  | fitted variables
fit_cont | 1 | Continuum fitting flag
cont_n | 3 | Degree of polynomial continuum
fit_res_gauss | 1 | Gaussian kernel
res_gauss | 3.5 | Kernel size (pixels)
kernfac | 6.0 | Kernel size measured in units of
|  | the kernel FWHM
list_molec | H2O, O2 | Molecules to be synthetised
Table 3: Physical and orbit parameters of the WASP-85 system

Description | Symbol | Value
---|---|---
Stellar Parameters |  |
$V$ magnitude | $m_{\rm v}$ |  $11.2\pm 0.011$ mag
Effective temp | $T_{\rm eff}$ |  $6112\pm 27$ K
Surface gravity | log $g_{\star}$ |  $4.48\pm 0.11$ cgs
Metallicity | [Fe/H] |  $0.00\pm 0.05$ dex
Stellar mass | $M_{\star}$ | $1.09\pm 0.08\,M_{\odot}$
Stellar radius | $R_{\star}$ | $0.935\pm 0.0023\,R_{\odot}$
Right ascension | R. A. | $11^{\rm h}43^{\rm m}38.01^{\rm s}$
Declination | Dec | $+06^{\circ}33^{\prime}49.4^{\prime\prime}$
Planet Parameter |  |
Planet mass | $M_{\rm p}$ | $1.265\pm 0.065\,M_{\rm Jup}$
Planet radius | $R_{\rm p}$ | $1.24\pm 0.03\,R_{\rm Jup}$
Planet density | $\rho$ |  $0.660\pm 0.020$ g cm-3
Eequilibrium temp | $T_{\rm eq}$ |  $1452\pm 6$ K
Radius ratio | $R_{\rm p}/R_{\star}$ | $0.0187\pm 0.00002$
Orbit Parameters |  |
Epoch$-2450000$ | $T_{\rm c}$ |  $6847.472856\pm 0.000014$ BJD
Semi-amplitude1 | $K_{\star}$ |  $173.3\pm 1.8$ m s-1
Period | $P$ |  $2.6556777\pm 0.00000044$ d
Transit duration | $T_{14}$ |  $0.10817\pm 0.00002$ d
Ingress duration | $T_{12}$ |  $0.013037\pm 0.000002$ d
Semi-major axis | $a$ |  $0.039\pm 0.001$ AU
Inclination | $i$ |  $89.69^{+0.11}_{-0.03}$ deg

Močnik et al. (2016), Brown et al. (2014)1

##  3 Data analysis

###  3.1 Telluric correction

The telluric correction was firstly 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.

Fig. 10 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$\alpha$ spectral line region (Bottom
panel) on the night of 2018 October 31, which illustrate that the telluric
features have been corrected quite well. The specific parameters used in this
work for telluric correction are listed in Table 2. 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 will overlay on 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 $\upsilon sin\textit{i}_{\star}$ and the
angle between the orbital plane and the apparent equatorial plane, namely the
spin-orbit angle $\lambda$.

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 use Markov chain Monte
Carlo (MCMC) algorithm which is available in emcee (Foreman-Mackey et al.,
2013) and employ 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 with 9 parameters, including the planet’s orbital period
$P$, the time of transit midpoint $T_{\rm c}$, the inclination of planetary
orbit plane $i$, the inclination of the stellar rotation axis $i_{\star}$, the
ratio of planetary and stellar radius $R_{\rm p}/R_{\star}$, the semi-major
axis $a$, the linear limb-darkening coefficient $\epsilon$, the stellar
angular rotation velocity $\Omega$, the sky-projected angle $\lambda$ which is
the angle between the stellar rotation axis and the normal of planetary orbit
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. Note that there are four
RV points in T3 deviating from the RV trend which 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_{\star}=K_{\star}\,{\rm sin}(2\pi\phi)+V_{\rm bary}+V_{\rm sys}$ (1)

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

We set $\lambda$, $\epsilon$, $\Omega$, and $i_{\star}$ as free parameters for
our model, with their priors listed in Table 4. The remaining parameters
including $a$, $P$, $i$, $\textit{T}_{c}$ and $R_{p}/R_{\star}$ are fixed to
the specific values as listed in Table 3. 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 4, whereas the retrieved
best model is compared with the data set in Fig. 2, with an rms residual of
$\sim 1.877$ m s-1. The derived projected spin-orbit angle $\lambda$ =
-16.155∘${}^{+2.916}_{-2.879}$, implying the planetary orbit plane is roughly
aligned to the stellar equator. We also derive the projected stellar rotation
velocity $v{\rm sin}\,i_{\star}=2.98^{+0.050}_{-0.043}$ $\rm km\,s^{-1}$, the
linear limb dark coefficient $\epsilon$ =$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. 11 in detail.

Figure 2: The RM effect fit using the GP+PyAstronomy package on the ESPRESSO
RVs taken from the spectrum headers. _Top panel:_ The 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:_ The residuals between the data and the model predictions.
Table 4: Parameters derived from the RV curve fitting.

Description | Symbol | Prior | Posterior |
---|---|---|---|---
Projected spin-orbit angle | $\lambda$ | U(-50, 50) |  $-16.155^{+2.916}_{-2.879}$ deg |
Linear limb dark coefficient | $\epsilon$ | U(0.5, 1) | $0.855^{+0.014}_{-0.014}$ |
Projected stellar rotation velocity | $\Omega$ | U($3\rm e^{-6}$, $1\rm e^{-5}$) |  $0.00000468^{+1.57\rm e^{-7}}_{-8.59\rm e^{-8}}$ rad s-1 |
Inclination of stellar rotation axis | $i_{\star}$ | U(70, 110) |  $91.6947^{+12.7688}_{-14.7692}$ deg |

Notes:U(a,b) represents a uniform distribution with a low and high limit of a
and b, respectively.

Table 5: Parameters derived from the GP analysis and the estimated noise of the three transit observations. Parameter | T1 | T2 | T3
---|---|---|---
$K_{\star}$ (m s-1) | $170^{+24}_{-15}$ | $157^{+4}_{-0}$ | $165^{+20}_{-12}$
$V_{\rm sys}$ (km s-1) | $13.531^{+0.004}_{-0.005}$ | $13.530^{+0.001}_{-0.001}$ | $13.519^{+0.024}_{-0.005}$
GP noise med. | 1.000105 | 0.999998 | 1.000064
GP noise std. | 0.000367 | 0.000087 | 0.000115

##  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$\alpha$, H$\beta$, 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 Section 3, the observed spectra have been corrected for
telluric contamination. Then the spectra are normalized to their corresponding
continuum level with a two-degree spline curve fit with the regions around
strong absorption lines excluded. We also apply a sigma-clipping rejection
algorithm on the normalized spectra and replace 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 shift the spectra to the stellar
rest frame, taking into account BERV, $V_{\rm sys}$ and the stellar reflex
motion induced by the planet $V_{\rm reflex}$. BERV is obtained directly from
the files header information of the ESPRESSO spectra, while $V_{\rm sys}$ is
derived from the RM analysis as described in Section 3.2 and $V_{\rm reflex}$
is calculated using the orbital parameters in Table 3. Then the aligned out-
of-transit spectra are averaged with their mean signal-to-noise ratio (SNR) as
their weights to construct a master out-of-transit spectrum, which is 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 planet rest frame using the planet radial velocity $V_{\rm p}$ determined
from $K_{\rm p}$ and the orbit phase of a given time with the following
formula:

$\Re(\lambda)=\sum_{in}{\frac{F_{in}(\lambda)}{F_{out}}}\bigg{|}_{\textnormal{{Planet\
RV\ Shift}}}\ -1$ (2)

Finally, the transmission spectrum is the SNR-weighted sum of in-transit
residuals.

Note that there is a wiggle pattern in the derived transmission spectrum, with
an amplitude of $\sim$0.075% and a period of $\sim$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. Fig. 3 shows an example of such correction.

Figure 3: An 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.

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 use the
Spectroscopy Made Easy tool (SME,Valenti & Piskunov (1996)) to compute the
theoretical stellar spectra at 21 different limb-darkening angles ($\mu$)
using the MARCS and VALD3 line lists (Ryabchikova et al., 2015). We employ
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_{\star}\times
0.01\,R_{\star}$, each element owning specific parameters including $v\,{\rm
sin}\,i_{\star}$, $\mu$, and $\theta$, 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 $\mu$ spectrum.

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

Figure 4: The 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 The strong atomic lines

We then performed a visual inspection on the phase-resolved transmission
spectra around each known strong feature, which will allow us to confirm or
deny the presence of the atoms of interest Tabernero et al. (2021). An
integrated SNR-weighted transmission spectrum will be obtained and fitted with
a Gaussian function to determine the center, depth, and width of each line.
Fig. 5 depicts the 2D phase-resolved transmission spectra and the integrated
spectrum around H$\alpha$, showing a clear detection of H absorption from the
WASP-85Ab’s atmosphere. Fig. 12 shows the same plots for other atomic lines,
including the Na I, Mg I, Li I, H$\alpha$, H$\beta$, 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_{\lambda}$, which can be derived with $\sqrt{1+h/\delta}$ $R_{p}$, where
$h$ is the line depth corresponding to different atoms, and $\delta$ is the
transit depth of the planet.

As shown in Table 6 and Fig. 5, H$\alpha$ 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 $\sim$0.625 $\pm$ 0.06 %, a full
width at half maximum (FWHM) of $\sim 23.806\pm 2.50$ km s-1, a $R_{\lambda}$
of $1.028\pm 0.003R_{\rm p}$. Clear absorption is detected for the Ca II
doublet in the integrated transmission spectra (cf. the two bottom panels of
Fig. 12). FWHM and $R_{\lambda}$ values are similar with those of H$\alpha$.
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
$\sim-4.805$ and 5.974 km s-1, 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 15, while the 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$\sigma$ detection for
Li I and no clear detection for the other species including Na I, Mg I, and
H$\beta$. For those with poor Gaussian fits, i.e., Mg I and H$\beta$, upper
limits of the absorption depths are estimated and listed in Table 6. As shown
in Fig. 14 and Fig. 15, 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
Section 2.

Note that the Na I D2 absorption line in T1 is abnormal, bearing huge error
bars as shown in Fig. 14. Through 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. 16. 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 perform
Empirical Monte Carlo (EMC) simulations (Redfield et al., 2008) to estimate
the effect of systematics, which is widely used previously (Wyttenbach et al.,
2015; Casasayas-Barris et al., 2019; Seidel et al., 2019; Allart et al.,
2020). The idea of this method is to artificially create new data set by
randomizing original data, to verify whether the investigated signals still
exist. We explore three scenarios covering in-in, in-out, out-out for each
strong line, with 10000 iterations for the $\sim$1.5 Å passband. The results
are shown in Fig. 17. The in-out distributions of Li I, H${\alpha}$, Ca II H
and Ca II K exhibit excess absorption, while the absorption depths in the in-
in and out-out scenarios center at zero. The EMC simulations suggest that the
signals of Li I, H${\alpha}$, Ca II H & K are likely to be created by the
transits and may origin from the planet. In addition, we also mask the line
cores with a width of 0.1 Å, where the SNRs are low and the RM+CLV residuals
may still exist. These results are shown in Fig. 13 and the fitted result is
listed in Table 6.

Figure 5: The phase-resolved 2D transmission spectra based on the combination
of three nights around the H$\alpha$ 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 planet
rest frame (PRF) assuming $K_{\rm p}$= 159.76 km s-1, and the vertical black-
dashed lines indicate the position of expected signal from the exoplanet
atmosphere in the planet rest frame. The right panel shows the combined
integrated transmission spectrum of H$\alpha$. The grey line represents the
transmission spectra in the planet rest frame which 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$\alpha$ line is shown in red. The
dashed green vertical line represents the static position of H$\alpha$ at
vacuum wavelength.  Table 6: The derived parameters of the atomic lines from
the 3-night combined transmission spectrum1

Line | $\lambda$ | h | Vwind | FWHM | Rλ |
---|---|---|---|---|---|---
| [nm] | [%] | [km s-1] | [km s-1] | [$R_{\rm p}$] |
Ca II K | 393.478 | 2.283 $\pm$ 0.389 | 8.536 $\pm$ 2.528 | 28.584$\pm$ 6.144 | 1.512$\pm$ 0.145 |
Ca II K (masked) | - | 1.769 $\pm$ 0.372 | 8.952 $\pm$ 3.359 | 30.705$\pm$ 8.180 | 1.413$\pm$ 0.148 |
Ca II H | 396.959 | 2.410 $\pm$ 0.476 | $-4.201\pm 2.378$ | 23.646 $\pm$ 5.713 | 1.535 $\pm$ 0.175 |
Ca II H (masked) | - | 1.729 $\pm$ 0.497 | $-3.495\pm 3.105$ | $21.294\pm 7.458$ | 1.405$\pm$ 0.199 |
Mg I | 457.238 |  $\leq 0.017$ $\pm$ 0.035 | - | - |  $\leq 1.005$ $\pm$ 0.01 |
H$\beta$ | 486.271 |  $\leq 0.050$ $\pm$ 0.010 | - | - |  $\leq 1.014$ $\pm$ 0.006 |
Na I D22 | 589.158 |  $\leq 0.069$ $\pm$ 0.073 | - | - |  $\leq 1.019$ $\pm$ 0.040 |
Na I D1 | 589.756 |  $\leq 0.067$ $\pm$ 0.077 | - | - |  $\leq 1.018$ $\pm$ 0.043 |
H$\alpha$ | 656.461 | 0.341 $\pm$ 0.120 | $-1.786\pm 3.098$ | 18.890 $\pm$ 7.109 | 1.091 $\pm$ 0.060 |
H$\alpha$ (masked) | - | 0.350 $\pm$ 0.14 | $-3.025\pm 2.431$ | $12.899\pm 5.659$ | $1.096\pm 0.072$ |
Li I | 670.961 | 0.089$\pm$ 0.033 | 4.223$\pm$ 6.167 | 33.951$\pm$ 14.425 | 1.025$\pm$ 0.018 |
Li I (masked) | - | $0.080\pm 0.031$ | $5.554\pm 6.514$ | $38.307\pm 15.243$ | $1.025\pm 0.017$ |

  * 1

The line center wavelength is in vacuum, the depth of absorption line $h$, the
Doppler shift of line center $V_{\rm wind}$, line width (FWHM) and the
effective planetary radius $R_{\lambda}$

  * 2

The Na I D2 measurement is obtained on the transmission spectra generated from
only the T2 and T3 data.

Figure 6: The CCF results derived by cross-correlating the template spectra
produced by petitRADTRANS with the observed transmission spectra. _First
panels_ : The 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_ : the $K_{p}$-$\Delta V_{\rm
sys}$ maps in the range of $-50\sim 350$ km s-1. The signal is expected to
appear around the intersection of two white dotted lines. _Fourth panels_ :
the SNR plots at the expected $K_{p}$ position.

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

###  4.2 Cross-correlation function analysis

Figure 8: The 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 addition to the direct measurements of strong lines, we employ 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 exoplanets, 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, Y I and
the molecules TiO and VO, and the spectrum template is computed by
petitRADTRANS code (Mollière et al., 2019). As input parameters, we assume the
WASP-85Ab radius is 1.24$R_{\rm J}$ (Močnik et al., 2016), the planet mass is
1.265$M_{\rm J}$, and thus log $g_{\rm p}=3.329$ cgs. We apply an isothermal
temperature of 1500 K and assume 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 $\pm 100$ km s-1 with a step of 0.5
km s-1, and calculate CCF for each residual spectrum and model spectrum. The
cross-correlation coefficients are calculated by the following formula:

$\centering
c(v,t)=\frac{\sum_{i}^{N}x_{i}(t)T_{i}(v)}{\sum_{i}^{N}T_{i}(v)},\@add@centering$
(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
2-dimensional 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_{\rm p}$ and system velocity $V_{\rm 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 which
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
Section 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. 8.

The obtained 2D CCF map with CLV+RM effects uncorrected 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. 6, where the inclined white dashes line indicates the expected track of
planet signal. The third panels show the $K_{\rm p}-\Delta V_{\rm sys}$ map,
in which the planet signal, if exists, should appear in the intersection of
white dashed lines, i.e, with a velocity close to $V_{\rm sys}$ and an RV
semi-amplitude $K_{\rm p}$ = 159.76 $\pm$ 4.09 km s-1. The fourth panels show
the SNR of the corresponding species at the expected $K_{\rm p}$. From Fig. 6,
only Li I is clearly detected with a maximum SNR of $\sim 4.5$, where the peak
position corresponds to $K_{\rm p}$ of $\sim 134.8$ $\pm$ 88.43 km s-1 and
$\Delta V_{\rm sys}\sim$ -4.6 $\pm$ 11.29 km s-1. Some remaining features in
the 2D maps may be associated with stellar activity (Stangret et al., 2021)
and propagated noise. The CCF maps and SNR plots for Mg I, Ti I, V I, Y I and
TiO and VO are shown in Fig. 7, with no detection reported.

##  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,
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^{\circ}}^{+2.916}_{-2.879}$, 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.

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

Among them, direct inspection revealed a $\sim 10\sigma$ absorption signature
at H$\alpha$ but it is not possible to confirm its origin as either stellar,
planetary, or else. In the meanwhile, the Ca II H & K lines are also
tentatively detected, with quite similar FWHM and consistent $R_{\lambda}$.
The potential spectral signals detected from H I and Ca II may arise from the
planetary atmosphere but we can not exclude yet the possibility of being
induced by stellar activity. Another species that may be detected is Li I,
which has 4.5 $\sigma$ significance at the estimated $K_{\rm p}$ velocity in
the $K_{\rm p}-\Delta V_{\rm sys}$ map. We note there is an offset of $21\pm
4.3$ km s-1 between the estimated $K_{\rm p}$ and the retrieved $K_{\rm p}$.
This discrepancy may be due to the large uncertainty of $K_{\rm p}$, given
that the Li CCF signal is quite extended in the $K_{\rm 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
be caused by the low SNRs 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 signals left. As shown in Fig. 13, the absorption features
of H${\alpha}$, Ca II H&K and Li I are still visible, although slightly
shallower. Therefore, the tentative detection of these three elements are
unlikely to be caused by the RM+CLV effects or the low SNRs 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 SNRs, 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). Borsa et al. (2021) and
Tabernero et al. (2021) present 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 $\sim$55 $M_{\rm Jup}$ do not deplete this element during their lifetime
(Chabrier & Baraffe, 2000; Baraffe et al., 2015). Such detection, if
confirmed, is an important step, which can help understand 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 estimate that the 5$\sigma$ upper limits of the excess absorption of the
non-detected atoms and molecules, which are listed 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. Fig. 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.

Table 7: Number of lines used in the CCF calculations and upper limits derived for the investigated species. Species | Number of lines | Detection limits
---|---|---
Mg I | 6 | 527.3 ppm
Li I | 8 | 323.9 ppm
Fe I | 34 | 341.8 ppm
Cr I | 45 | 264.1 ppm
K I | 53 | 364.2 ppm
Y I | 99 | 283.8 ppm
V I | 154 | 301.6 ppm
Ti I | 159 | 377.1 ppm
TiO | 7832 | 27.6 ppm
VO | 8280 | 22.3 ppm

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_{\rm
eq}$ and small $R_{\rm p}$ and thus small Transmission Spectroscopy Metric
(Kempton et al., 2018), and 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 three
nights barely 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 SNR 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 eliminate the
influence of the variation of stellar line profile during transit.

We also note that the three nights data set does not comply with each other
well. It seems the last night has 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.

###### 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. MZ, YQS, QLOY are supported
by the Chinese Academy of Sciences (CAS), through a grant to the CAS South
America Center for Astronomy (CASSACA) in Santiago, Chile. HMC acknowledges
support from a UKRI Future Leaders Fellowship (MR/S035214/1).

## References

  * Allart et al. (2017) Allart, R., Lovis, C., Pino, L., et al. 2017, A&A, 606, A144
  * Allart et al. (2020) Allart, R., Pino, L., Lovis, C., et al. 2020, A&A, 644, A155
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Birkby (2018) Birkby, J. L. 2018, Exoplanet Atmospheres at High Spectral Resolution
  * Borsa et al. (2021) Borsa, F., Allart, R., Casasayas-Barris, N., et al. 2021, A&A, 645, A24
  * Bouvier (2008) Bouvier, J. 2008, A&A, 489, L53
  * Brown et al. (2014) Brown, D. J. A., Anderson, D. R., Doyle, A. P., et al. 2014, arXiv e-prints, arXiv:1412.7761
  * Butler et al. (1996) Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, PASP, 108, 500
  * Casasayas-Barris et al. (2022) Casasayas-Barris, N., Borsa, F., Palle, E., et al. 2022, A&A, 664, A121
  * Casasayas-Barris et al. (2019) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2019, A&A, 628, A9
  * Chabrier & Baraffe (2000) Chabrier, G. & Baraffe, I. 2000, ARA&A, 38, 337
  * Charbonneau et al. (2002) Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377
  * Chen et al. (2020) Chen, G., Casasayas-Barris, N., Pallé, E., et al. 2020, A&A, 635, A171
  * Chen et al. (2018) Chen, G., Pallé, E., Welbanks, L., et al. 2018, A&A, 616, A145
  * Czesla et al. (2019) Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages, Astrophysics Source Code Library, record ascl:1906.010
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, The Astronomical Journal, 154, 220
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Giacobbe et al. (2021) Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205
  * Hoeijmakers et al. (2019) Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165
  * Ishizuka et al. (2021) Ishizuka, M., Kawahara, H., Nugroho, S. K., et al. 2021, AJ, 161, 153
  * Kempton et al. (2018) Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401
  * Lira-Barria et al. (2022) Lira-Barria, A., Rojo, P. M., & Mendez, R. A. 2022, A&A, 657, A36
  * Madhusudhan (2019) Madhusudhan, N. 2019, ARA&A, 57, 617
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
  * Merritt et al. (2021) Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2021, MNRAS, 506, 3853
  * Močnik et al. (2016) Močnik, T., Clark, B. J. M., Anderson, D. R., Hellier, C., & Brown, D. J. A. 2016, The Astronomical Journal, 151, 150
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67
  * Ohta et al. (2005) Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96
  * Pepe et al. (2002) Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632
  * Perruchot et al. (2008) Perruchot, S., Kohler, D., Bouchy, F., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70140J
  * Prinoth et al. (2022) Prinoth, B., Hoeijmakers, H. J., Kitzmann, D., et al. 2022, Nature Astronomy, 6, 449
  * Queloz et al. (2000) Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, 359, L13
  * Redfield et al. (2008) Redfield, S., Endl, M., Cochran, W. D., & Koesterke, L. 2008, ApJ, 673, L87
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005
  * Seidel et al. (2019) Seidel, J. V., Ehrenreich, D., Wyttenbach, A., et al. 2019, A&A, 623, A166
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  * Stangret et al. (2021) Stangret, M., Pallé, E., Casasayas-Barris, N., et al. 2021, A&A, 654, A73
  * Tabernero et al. (2021) Tabernero, H. M., Zapatero Osorio, M. R., Allart, R., et al. 2021, A&A, 646, A158
  * Triaud (2018) Triaud, A. H. M. J. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 2
  * Valenti & Piskunov (1996) Valenti, J. A. & Piskunov, N. 1996, A&AS, 118, 595
  * Wyttenbach et al. (2015) Wyttenbach, A., Ehrenreich, D., Lovis, C., Udry, S., & Pepe, F. 2015, A&A, 577, A62
  * Yan & Henning (2018) Yan, F. & Henning, T. 2018, Nature Astronomy, 2, 714

##  Appendix A Additional figures

Figure 10: The effect of 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: The telluric correction of spectral line around Na D1 $\&$ D2 lines.
Bottom panel: The telluric correction of spectral line around H$\alpha$ lines.
Figure 11: Posterior distribution of the best-fit RM model using PyAstronomy
and GP.

Figure 12: Same as Fig. 5,but for Na I, Mg I, H${\beta}$, Li I, Ca II H and Ca
II K.

Figure 13: The masked phase-resolved 2D transmission spectra based on the
combination of three nights for H$\alpha$, 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 planet rest frame (PRF)
assuming $K_{\rm p}$= 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.

Figure 14: Transmission spectra of Na I, Mg I, H${\alpha}$, H${\beta}$ for
each night. Left (T1), middle (T2), right (T3). The grey line represents the
raw transmission spectra in the planet rest frame 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.

Figure 15: Same as Fig 14, but for Ca II H,Ca II K and Li I. Figure 16: This
left panel shows that Na I D2 absorption line in the raw data image of T1 is
contaminated by two bright spikes which leads to the abnormal flux around Na I
D2 line. The right panel plots the abnormal counts of no sky-subtracted flux
colored in blue and sky-subtracted flux shown in red.

Figure 17: The results of EMC simulation for Na I, Mg I, Li I, H${\alpha}$,
H${\beta}$, Ca II H and Ca II K combined with 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; blue is the out-out
distribution. The in-out distribution of Li I, H${\alpha}$, Ca II H and Ca II
K exhibits excess absorption.