* Colpi  _et al._ (1993) M. Colpi, S. L. Shapiro, and S. A. Teukolsky, Astrophys. J. 414, 717 (1993).
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astronomy and Astrophysics 380, 151?167 (2001).
  * Read  _et al._ (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009).
  * Paschalidis and Stergioulas (2017) V. Paschalidis and N. Stergioulas, Living Rev. Rel. 20, 7 (2017), arXiv:1612.03050 [astro-ph.HE] .
  * Bhattacharyya  _et al._ (2010) A. Bhattacharyya, I. N. Mishustin, and W. Greiner, Journal of Physics G: Nuclear and Particle Physics 37, 025201 (2010).
  * Yasutake  _et al._ (2011) N. Yasutake, G. F. Burgio, and H. J. Schulze, Physics of Atomic Nuclei 74, 1502?1507 (2011).
  * Duez  _et al._ (2005) M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev. D72, 024029 (2005), arXiv:astro-ph/0503421 [astro-ph] .
  * Etienne  _et al._ (2010a) Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Phys. Rev. D 82, 084031 (2010a).
  * Etienne  _et al._ (2015) Z. B. Etienne, V. Paschalidis, R. Haas, P. Moesta, and S. L. Shapiro, Class. Quant. Grav. 32, 175009 (2015), arXiv:1501.07276 [astro-ph.HE] .
  * Allen  _et al._ (2001) G. Allen, D. Angulo, I. T. Foster, G. Lanfermann, C. Liu, T. Radke, E. Seidel, and J. Shalf, CoRR cs.DC/0108001 (2001).
  * Schnetter  _et al._ (2004) E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quant. Grav. 21, 1465 (2004), arXiv:gr-qc/0310042 [gr-qc] .
  * Schnetter  _et al._ (2006) E. Schnetter, P. Diener, E. N. Dorband, and M. Tiglio, Class. Quant. Grav. 23, S553 (2006), arXiv:gr-qc/0602104 [gr-qc] .
  * Shibata and Nakamura (1995) M. Shibata and T. Nakamura, Phys. Rev. D52, 5428 (1995).
  * Baumgarte and Shapiro (1999) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, 024007 (1999), arXiv:gr-qc/9810065 [gr-qc] .
  * Bona  _et al._ (1995) C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. Lett. 75, 600 (1995), arXiv:gr-qc/9412071 [gr-qc] .
  * Alcubierre  _et al._ (2001) M. Alcubierre, B. Bruegmann, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D64, 061501 (2001), arXiv:gr-qc/0104020 [gr-qc] .
  * Alcubierre  _et al._ (2003) M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D67, 084023 (2003), arXiv:gr-qc/0206072 [gr-qc] .
  * Etienne  _et al._ (2008) Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro, K. Taniguchi, and T. W. Baumgarte, Phys. Rev. D77, 084002 (2008), arXiv:0712.2460 [astro-ph] .
  * Etienne  _et al._ (2010b) Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Phys.Rev. D82, 084031 (2010b).
  * Etienne  _et al._ (2012) Z. B. Etienne, V. Paschalidis, Y. T. Liu, and S. L. Shapiro, Phys.Rev. D85, 024013 (2012).
  * Espino and Paschalidis (2019) P. Espino and V. Paschalidis, Phys. Rev. D99, 083017 (2019), arXiv:1901.05479 [astro-ph.HE] .
  * Newman and Penrose (1962) E. Newman and R. Penrose, J. Math. Phys. 3, 566 (1962).
  * Penrose (1963) R. Penrose, Phys. Rev. Lett. 10, 66 (1963).
  * Reisswig  _et al._ (2011) C. Reisswig, C. D. Ott, U. Sperhake, and E. Schnetter, Physical Review D 83 (2011), 10.1103/physrevd.83.064008.
  * Bozzola (2021) G. Bozzola, Journal of Open Source Software 6, 3099 (2021).
  * Guercilena  _et al._ (2017) F. Guercilena, D. Radice, and L. Rezzolla, Comput. Astrophys. Cosmol. 4, 3 (2017), arXiv:1612.06251 [gr-qc] .
  * Manukovskii (2010) K. V. Manukovskii, Astronomy Letters 36, 191 (2010).
  * Yudin  _et al._ (2020) A. V. Yudin, T. L. Razinkova, and S. I. Blinnikov, Astronomy Letters 45, 847 (2020), arXiv:2005.04428 [astro-ph.HE] .
  * Lattimer and Schramm (1974) J. Lattimer and D. Schramm, Astrophys. J. Lett. 192, L145 (1974).
  * East  _et al._ (2012) W. E. East, F. Pretorius, and B. C. Stephens, Physical Review D 85 (2012), 10.1103/physrevd.85.124009.
  * East  _et al._ (2015) W. E. East, V. Paschalidis, and F. Pretorius, Astrophys. J. 807, L3 (2015), arXiv:1503.07171 [astro-ph.HE] .
  * East  _et al._ (2016) W. E. East, V. Paschalidis, F. Pretorius, and S. L. Shapiro, prd 93, 024011 (2016), arXiv:1511.01093 [astro-ph.HE] .
  * Glendenning  _et al._ (1997) N. K. Glendenning, S. Pei, and F. Weber, Phys. Rev. Lett. 79, 1603 (1997).
  * Zdunik  _et al._ (2006) J. Zdunik, M. Bejger, P. Haensel, and E. Gourgoulhon, Astron. Astrophys. 450, 747 (2006), arXiv:astro-ph/0509806 .
  * Hessels  _et al._ (2006) J. W. T. Hessels, S. M. Ransom, I. H. Stairs, P. C. C. Freire, V. M. Kaspi, and F. Camilo, in  _American Astronomical Society Meeting Abstracts #207_ , American Astronomical Society Meeting Abstracts, Vol. 207 (2006) p. 209.07.
  * Papitto  _et al._ (2014) A. Papitto, D. F. Torres, N. Rea, and T. M. Tauris, A&A 566, A64 (2014), arXiv:1403.6775 [astro-ph.HE] .
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Moore  _et al._ (2015) C. J. Moore, R. H. Cole, and C. P. L. Berry, Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Paschalidis  _et al._ (2015) V. Paschalidis, W. E. East, F. Pretorius, and S. L. Shapiro, Phys. Rev. D92, 121502 (2015), arXiv:1510.03432 [astro-ph.HE] .
  * Tsokaros  _et al._ (2019) A. Tsokaros, M. Ruiz, L. Sun, S. L. Shapiro, and K. Uryū, Phys. Rev. Lett. 123, 231103 (2019), arXiv:1907.03765 [gr-qc] .
  * Ruiz  _et al._ (2019) M. Ruiz, A. Tsokaros, V. Paschalidis, and S. L. Shapiro, Phys. Rev. D99, 084032 (2019), arXiv:1902.08636 [astro-ph.HE] .
  * Nathanail (2020) A. Nathanail, The Astrophysical Journal 892, 35 (2020).
  * Ruiz  _et al._ (2020) M. Ruiz, A. Tsokaros, and S. L. Shapiro, Phys. Rev. D 101, 064042 (2020), arXiv:2001.09153 [astro-ph.HE] .


11institutetext:  Lund Observatory, Division of Astrophysics, Department of
Physics, Lund University, Box 43, 221 00 Lund, Sweden 22institutetext:
Department of Physics and Trottier Institute for Research on Exoplanets,
Université de Montréal, Montreal, QC, Canada 33institutetext: University of
Bern, Physics Institute, Division of Space Research & Planetary Sciences,
Gesellschaftsstr. 6, 3012, Bern, Switzerland 44institutetext: Space Telescope
Science Institute, Baltimore, MD 21218, USA 55institutetext: Department of
Astronomy & Astrophysics, The University of Chicago, Chicago, IL, USA
66institutetext: European Southern Observatory, Alonso de Córdova 3107,
Vitacura, Región Metropolitana, Chile 77institutetext: Material Sciences and
Applied Mathematics, Malmö University, 205 06, Malmö, Sweden 88institutetext:
ETH Zurich, Institute for Particle Physics and Astrophysics; Wolfgang-Pauli-
Strasse 27, CH-8093 Zürich, Switzerland 99institutetext: Landessternwarte,
Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, 69117,
Heidelberg, Germany 1010institutetext: Observatoire de la Côte d’Azur, CNRS
UMR 7293, BP4229, Laboratoire Lagrange, F-06304 Nice Cedex 4, France

# Time-resolved transmission spectroscopy of the ultra-hot Jupiter WASP-189 b

B. Prinoth 11 H. J. Hoeijmakers 11 S. Pelletier 22 D. Kitzmann 33 B. M. Morris
44 A. Seifahrt 55 D. Kasper 55 Heidi H. Korhonen 66 M. Burheim 1177 J. L. Bean
55 B. Benneke 22 N. W. Borsato 11 M. Brady 55 S. L. Grimm 3388 R. Luque 55 J.
Stürmer 99 B. Thorsbro 101011<EMAIL_ADDRESS>
Ultra-hot Jupiters are tidally locked with their host stars, dividing their
atmospheres into a hot dayside and a colder nightside. As the planet moves
through transit, different regions of the atmosphere rotate into view,
revealing different chemical regimes. High-resolution spectrographs can
observe asymmetries and velocity shifts and offer the possibility for time-
resolved spectroscopy. The ultra-hot Jupiter WASP-189 b has recently been
found to possess a rich transmission spectrum with evidence for atmospheric
dynamics and chemical inhomogeneity. In this study, we search for other atoms
and molecules in the planet’s transmission spectrum and investigate asymmetric
signals. We analysed and combined eight transits of the ultra-hot Jupiter
WASP-189 b collected with the HARPS, HARPS-N, ESPRESSO, and MAROON-X high-
resolution spectrographs. Using the cross-correlation technique, we searched
for neutral and ionised atoms as well as oxides, and we compared the obtained
signals to model predictions. We report significant detections for H, Na, Mg,
Ca, Ca+, Ti, Ti+, TiO, V, Cr, Mn, Fe, Fe+, Ni, Sr, Sr+, and Ba+. Of these, Sr,
Sr+, and Ba+ are detected for the first time in the transmission spectrum of
WASP-189 b. In addition, we robustly confirm the detection of titanium oxide
based on observations with HARPS and HARPS-N using the follow-up observations
performed with MAROON-X and ESPRESSO. By fitting the orbital traces of the
detected species by means of time-resolved spectroscopy using a Bayesian
framework, we inferred posterior distributions for orbital parameters as well
as line shapes. Our results indicate that different species must originate
from different regions of the atmosphere to be able to explain the observed
time dependence of the signals. Throughout the course of the transit, most
signal strengths are expected to increase due to the larger atmospheric scale
height at the hotter trailing terminator. For some species, however, we
instead observed that the signals weaken, either due to the ionisation of
atoms and their ions or the dissociation of molecules on the dayside.

###### Key Words.:

planets and satellites: atmospheres, planets and satellites: gaseous planets,
planets and satellites: individual: WASP-189 b, techniques: spectroscopic

††offprints: Bibiana Prinoth,

##  1 Introduction

Figure 1: Schematic of the contribution of different regions of the planetary
atmosphere to the observed absorption signal over the course of the transit
from left to right. As the planet rotates over the course of the transit,
different regions become observable and dominate the contributions. Generally,
the leading terminator (3, 4) is colder than the trailing terminator (1,2).
The atmosphere at the trailing terminator is therefore expected to be more
inflated (puffy) and generally more susceptible to transmission spectroscopy.
The terminator regions that dominate the observed absorption signal are
denoted with black circles. The colouring of the atmosphere is taken from Lee
et al. (2022) to reflect the temperature of the atmosphere at a pressure of
$1\text{\,}\mathrm{m}\mathrm{b}\mathrm{a}\mathrm{r}$. The shape of the
atmosphere with an increased radius on the dayside is adopted from Wardenier
et al. (2021). We note that the schematic is not to scale. The changing
viewing angle from the start to the centre of the transit for WASP-189 b is
$\alpha=\arcsin\left(\sqrt{1-b^{2}}\frac{R_{\ast}}{a}\right)\approx
11,33^{\circ}$, where $b$ is the impact parameter ($b=$
$0.433^{+0.014}_{-0.015}$, Deline et al. 2022), $R_{\ast}$ is the stellar
radius, and $a$ is the semi-major axis.

Our understanding of planetary formation and evolution has been transformed by
the discovery and characterisation of exoplanets over the past three decades.
The hottest of these exoplanets are gas giant planets, often termed ultra-hot
Jupiters, and they reside close to their host stars and typically have orbital
periods of less than 10 days (Fortney et al. 2021). As a result, ultra-hot
Jupiters are subject to extreme variations in the temperature and chemistry of
their atmosphere since they are thought to enter tidal confinement soon after
their formation (Showman et al. 2009; Arras & Socrates 2010).

This tidal confinement, or tidal locking, of ultra-hot Jupiters divides their
atmosphere into a permanently irradiated hot dayside and a cooler, permanently
dark nightside. On the hot dayside, the high temperatures trigger processes
such as thermal dissociation of most molecules, partial thermal ionisation of
atomic species, and thermal inversions by metals and their oxides (Lothringer
et al. 2018; Kitzmann et al. 2018; Parmentier et al. 2018; Arcangeli et al.
2018). On the colder nightside, we expect H2 recombination (Bell & Cowan
2018), cloud formation by condensation, rain out, and cold traps of certain
elements to the lower region of the planets, all of which are invisible to the
distant observers (Spiegel et al. 2009; Helling et al. 2021; Komacek et al.
2022; Kesseli et al. 2022; Hoeijmakers et al. 2022; Gandhi et al. 2023;
Pelletier et al. 2023). The close proximity of ultra-hot Jupiters to their
host stars and their inflated atmospheres make them favourable targets for the
study of atmospheric chemistry and dynamics using high-resolution transmission
and emission spectroscopy. Studies of ultra-hot Jupiters at optical
wavelengths include, amongst others, the studies of KELT-9 b (Hoeijmakers et
al. 2018, 2019; Cauley et al. 2019; Pino et al. 2020, 2022; Borsato et al.
2023), KELT-20 b (Casasayas-Barris et al. 2018, 2019; Hoeijmakers et al.
2020a; Stangret et al. 2020), WASP-76 b (Ehrenreich et al. 2020; Tabernero et
al. 2021; Seidel et al. 2021; Kesseli et al. 2022; Pelletier et al. 2023),
WASP-121 b (Hoeijmakers et al. 2020b, 2022; Borsa et al. 2021), WASP-189 b
(Yan et al. 2020; Prinoth et al. 2022), WASP-33 b (Cauley et al. 2021; Yan et
al. 2021), and WASP-19 b (Sedaghati et al. 2021). In addition, the
implementation of General Circulation Models using Monte Carlo simulations has
significantly contributed to shaping our understanding of the atmospheric
dynamics of ultra-hot Jupiters (Beltz et al. 2021; Wardenier et al. 2021;
Beltz et al. 2022; Lee et al. 2022; Savel et al. 2022; Beltz et al. 2023;
Savel et al. 2023).

During a transit event, a tidally locked planet that is sufficiently close to
its host star will rotate significantly during transit in such a way that
different regions of the atmosphere become observable as the planet moves in
front of the star on ingress, passes through the centre during mid-transit,
and then leaves the line of sight on egress (Wardenier et al. 2021). In the
case of some ultra-hot Jupiters, the rotation angle during transit is on the
order of 20 to 40 degrees (e.g. WASP-189 b: $\approx$ 23∘; WASP-76 b:
$\approx$ 28∘; WASP-121 b: $\approx$ 30∘; KELT-9 b: $\approx$ 36∘), causing
different regions of the planetary atmosphere to contribute to the observed
absorption spectrum. At the very beginning of the transit event, it is the hot
side of the leading (morning) terminator that dominates the observed
absorption, as the line of sight passes through more of the hotter regions at
the leading terminator than at the trailing (evening) terminator (see Fig. 1).
In the absence of any atmospheric dynamics, apart from co-rotation, this would
result in a redshifted signal as the planet rotates counterclockwise.
Similarly, towards the end of the transit, the absorption arises predominantly
from the hot side at the trailing terminator, blueshifted by the planet’s
rotation. These effects lead to deviations in the apparent orbital velocity of
the planet (Ehrenreich et al. 2020; Kesseli & Snellen 2021; Wardenier et al.
2021; Kesseli et al. 2022; Prinoth et al. 2022; Pelletier et al. 2023).
Deviations in the observed systemic velocity of such planets suggest the
presence of global dayside-to-nightside winds, while significant line
broadening suggests the presence of a vertical outflow (e.g. Seidel et al.
2021; Kesseli et al. 2022; Prinoth et al. 2022). Variations in the detected
systemic and orbital velocities of different chemical species suggest that
their absorption originates in different layers in the atmosphere (Prinoth et
al. 2022; Maguire et al. 2023). These variations also suggest the presence of
condensation and rain out on the nightside (Ehrenreich et al. 2020; Kesseli et
al. 2022).

Detections of absorbing species in the atmospheres of exoplanets are nowadays
routinely made using the cross-correlation technique (Snellen et al. 2010). A
key observable is obtained by stacking the time series of cross-correlation
functions in the co-moving frame of the planet. In single exposures, the
signal of the atmosphere may be too faint to be detected, whereas stacking in
the rest frame of the planet increases the signal-to-noise ratio. While it is
possible to optimise the exposure time to ensure a relatively high signal-to-
noise ratio while avoiding smearing for a single transit event itself (Boldt-
Christmas et al. subm.), a high signal-to-noise ratio can also be achieved by
combining the data from multiple transit events. If the signal-to-noise ratio
is sufficiently high, stacking in the rest frame of the planet is not
necessary to see the atmospheric signature. In this case, the absorption
signal can be studied using time-resolved spectroscopy to reveal empirical
evidence for the presence of dynamics and chemical regimes probed at different
times during the transit. While stacking in the rest frame of the planet
effectively averages over the contributions from the different terminator
regions, time-resolved spectroscopy allows the asymmetry in line shape and
strength over the course of the transit to be studied and linked to the
underlying chemistry and dynamics.

Due to the strong thermal gradient caused by tidal locking, different regions
of the atmosphere may exhibit distinctly different chemical regimes. The
hotspot may be significantly offset from the sub-stellar point (Showman &
Polvani 2011), making the trailing terminator hotter than the leading
terminator, as it has a larger atmospheric scale height and deeper absorption
lines when seen in transmission. In this region, the atmosphere may be hot
enough to trigger the ionisation of atoms at certain altitudes. As an example,
the ionisation of Ba to Ba+ requires less energy than the ionisation of
lighter elements (e.g. Ca and Mg to Ca+ and Mg+, respectively) within the same
group. This is due to the fact that the effective nuclear charge experienced
by the outer s-electron is constant within the group, whereas the orbital
radius increases with the weight of the elements. Similarly, less energy is
needed to ionise neutral atoms than singly-ionised atoms. The screening on the
outer electron in an ionised atom is less efficient compared to the neutral
atom, resulting in a larger effective nuclear charge and thus higher
ionisation energy (Froese-Fischer et al. 1997). Similarly, species that are
present on the nightside of the trailing terminator may condense out of the
gas phase on the cooler nightside. Depending on the veracity of atmospheric
dynamics, they may either get remixed on the dayside or be cold trapped (e.g.
Merritt et al. 2021; Pelletier et al. 2023). Provided sufficient signal-to-
noise, such effects may be observable using time-resolved spectroscopy. In
this work, we present such observations of WASP-189 b, an ultra-hot Jupiter.

The ultra-hot Jupiter WASP-189 b is a gas giant planet orbiting the bright
A-type star HD 133112 (V = 6.6 mag; Anderson et al. 2018). It is one of the
most intensely irradiated planets known to date, with an equilibrium
temperature of $2641\pm 34\text{\,}\mathrm{K}$ (Anderson et al. 2018), and it
resides on a 2.7 day long orbit that is most likely tidally circularised and
synchronised. Recent observations by CHEOPS yielded constraining measurements
of the planet radius, ephemeris, and a dayside temperature of
$3400\text{\,}\mathrm{K}$ (Lendl et al. 2020). Follow-up and refinement of the
system parameters have been given by Deline et al. (2022). The atmosphere of
WASP-189 b is inverted (Yan et al. 2020, 2022), based on the detection of
emission by atomic iron and CO, which may in part be due to the presence of
titanium oxide in its atmosphere (Prinoth et al. 2022). Observations with
HARPS (Cosentino et al. 2012) and HARPS-N (Mayor et al. 2003) on the ESO 3.6m
telescope on La Silla and the 3.58m telescope on the Telescopio Nazionale
Galilei on La Palma, respectively, have revealed Fe, Fe+, and Ti (Stangret et
al. 2022; Prinoth et al. 2022), as well as Ti+, V, Mg, Mn, Cr, and TiO
(Prinoth et al. 2022). Tentative detections include H and Ca+ (Stangret et al.
2022) as well as Cr+, Sc+, Na, Ni, and Ca (Prinoth et al. 2022). These
detections already showed evidence of dynamical effects (such as winds)
through deviations in the apparent systemic and orbital velocities. Generally,
the detected species appear blueshifted with respect to the expected systemic
velocity due to prevailing day-to-night winds.

In this work, we present the combined analysis of the high-resolution
transmission spectrum of WASP-189 b using eight transit events observed with
the HARPS, HARPS-N, ESPRESSO, and MAROON-X spectrographs. We apply the cross-
correlation technique (Snellen et al. 2010) in order to search for atomic and
ionised metals as well as molecules, and we fit the traces in the cross-
correlation functions to determine spectral and orbital parameters using time-
resolved spectroscopy. Section 2 describes the observations and data reduction
procedures. Section 3 describes the telluric and preparatory corrections and
further details the analysis process for the cross-correlation technique and
trace fitting of this dataset, which presents the basis of our findings. In
Section 4, we present the results of the analysis and discuss their
implications. Conclusions are presented in Section 5.

##  2 Observations and data reduction

Figure 2: Signal-to-noise ratio versus orbital phase of each of the eight
transit time series indicated in Table 1. The vertical dashed lines indicate
the start and the end of the transit. The signal-to-noise ratio was calculated
by averaging over all orders.

We observed one transit of the ultra-hot Jupiter WASP-189 b using the ESPRESSO
echelle spectrograph at the ESO Very Large Telescope in Chile (Pepe et al.
2021) and two transits with the MAROON-X high-resolution optical spectrograph
at the 8.1 m Gemini-North observatory in Hawaii (Seifahrt et al. 2018, 2020).
The transits were observed on 2021-06-04 with ESPRESSO (programme ID 107.22QF,
PI: Prinoth) and on 2022-04-03 and 2022-06-02 with MAROON-X (programme ID
GN-2022A-FT-208, PI: Pelletier).

The ESPRESSO observations were performed in high-resolution mode with a
spectral resolution of $R=\frac{\lambda}{\Delta\lambda}=$140000$$ covering the
wavelength range between $$380$-$788\text{\,}\mathrm{n}\mathrm{m}$$. It is
known that the readout speed of $500$ kpx s-1 associated with the unbinned
(1x1) readout mode introduces correlated noise at the bias level that is much
larger than with the binned (2x1) readout mode. Therefore, we chose the binned
readout mode, which is associated with a readout speed of $100$ kpx s-1. Due
to the initial poor weather conditions, the telescope dome was not opened
until a few minutes before the start of the transit, causing the ingress to be
partially missed. The observations of the in-transit spectra were obtained
with UT-2. Due to a technical problem with UT-2 towards the end of the
transit, the observations were continued with UT-1, causing the egress of the
planetary transit to be lost. The out-of-transit observations were obtained
with UT-1. As described in the ESPRESSO user manual, ESPRESSO spectra are
influenced by two different interference patterns induced in the Coudé train
optics.111https://www.eso.org/sci/facilities/paranal/instruments/espresso/ESPRESSO_User_Manual_P108_v2.png
It is known that approximately sinusoidal wiggles with periods of 30 and 1 Å
and amplitudes of $\sim 1\%$ and $\sim 0.1\%$, respectively, become visible
when spectra are recorded at different telescope positions and then divided by
each other (Allart et al. 2020). The high-frequency pattern (1 Å) is below the
photon-noise level, while the low-frequency pattern was corrected as discussed
in Section 3.2. Science observations were performed with fibre A on the target
and fibre B on the sky.

Table 1: Overview of observations.

Date | Instrument | PI, Prog. # |  $\\#$ Spectraa |  texp [s] b | Wavelength |  Resolution c | Mirror
---|---|---|---|---|---|---|---
|  |  |  |  |  coverage (wlc) [nm] c |  |  size [m] c
2018-03-26 | HARPS | Hellier, 0100.C-0847 | 39 (15/24) | 600 | 383-693 | 120,000 | 3.6
2019-04-14 | HARPS | Hoeijmakers, 0103.C-0472 | 126 (69/57) | 200 | 383-693 | 120,000 | 3.6
2019-04-25 | HARPS | Hoeijmakers, 0103.C-0472 | 107 (52/55) | 200 | 383-693 | 120,000 | 3.6
2019-05-06 | HARPS-N | Casasayas-Barris, CAT19A_97 | 112 (69/43) | 200 | 383-690 | 115,000 | 3.6
2019-05-14 | HARPS | Hoeijmakers, 0103.C-0472 | 122 (68/54) | 200 | 383-693 | 120,000 | 3.6
2021-06-04 | ESPRESSO | Prinoth, 107.22QF | 55 (39/16) | 300 | 380-788 | 140,000 | 8.2
2022-04-03 | MAROON-X | Pelletier, GN-2022A-FT-208 | 57 (38/19) | b: 200 | b: 490-720 | 80,000 | 8.1
|  |  |  | r: 160 | r: 640-920 |  |
2022-06-02 | MAROON-X | Pelletier, GN-2022A-FT-208 | 57 (40/17) | b: 200 | b: 490-720 | 80,000 | 8.1
|  |  |  | r: 160 | r: 640-920 |  |

Note: (a) In parentheses, in-transit and out-of-transit spectra, respectively.
(b) Exposure times differ for the blue and red arm (channel) of MAROON-X. (c)
Adapted from Kitzmann et al. (2023) and references therein.

The observations performed with ESPRESSO were reduced using the ESPRESSO data
reduction workflow (pipeline version 2.3.3) provided by ESO and executed via
EsoReflex (version
3.13.5).222https://www.eso.org/sci/software/pipelines/espresso/espresso-pipe-
recipes.html The pipeline provides both one- and two-dimensional products with
and without sky subtraction. The one-dimensional products contain a stitched,
blaze-corrected spectrum, corrected for Earth’s velocity around the barycentre
of the Solar system (BERV) with both vacuum and air wavelength solutions.
Within the scope of this paper, we generally denote the wavelength solution in
air, unless mentioned otherwise. The two-dimensional products are saved per
exposure. Each file contains all orders, including the wavelength solution in
air.

The observations with MAROON-X cover the wavelength range between $490$ and
$900\text{\,}\mathrm{n}\mathrm{m}$ with two arms (red and blue) at a spectral
resolution of $\frac{\lambda}{\Delta\lambda}\approx$85000$$. The blue and red
arms cover the wavelength ranges between $490$ and
$720\text{\,}\mathrm{n}\mathrm{m}$ and $640$ and
$920\text{\,}\mathrm{n}\mathrm{m}$ respectively. The observations cover the
full transit, including baseline observations before and after the transit.

The MAROON-X observations were reduced using the standard pipeline (Seifahrt
et al. 2020), which extracts the detector images into one-dimensional
wavelength-calibrated spectra order by order for each of the exposures in the
time series. For each of the channels (arms), the pipeline outputs a high-
resolution spectral time series in the form of $N_{\rm exp}\times N_{\rm
order}\times N_{\rm pixel}$ for each transit sequence. Respectively, $N_{\rm
exp}$, $N_{\rm order}$, and $N_{\rm pixel}$ denote the number of exposures,
orders, and pixels for each detector. The pipeline provides the wavelength
solution in vacuum, but we changed it to air for consistency between the
observations.

In addition to the newly obtained observations with ESPRESSO and MAROON-X, we
used the transit observations of WASP-189 b with HARPS-N (Cosentino et al.
2012) and HARPS (Mayor et al. 2003) previously published in Anderson et al.
(2018), Stangret et al. (2022), and Prinoth et al. (2022). We used the reduced
data as described in Prinoth et al. (2022). A log of the time series
observations used in this work is provided in Table 1, and an overview of the
coverage of the observations is provided in Fig. 2.

##  3 Methods

We follow the methodology as in Prinoth et al. (2022) and discuss deviations
and necessary clarifications here in detail. We further discuss the attempts
to correct the problems with the interference patterns in the ESPRESSO data.

###  3.1 Telluric correction

Figure 3: Extracted one-dimensional spectra and wavelength inclusion regions
for telluric correction with molecfit for HARPS/ HARPS-N (top), ESPRESSO
(middle), and MAROON-X (bottom). We note here that due to the simple de-
blazing of the MAROON-X spectra, the continuum is not flat but instead has a
slope. Molecfit takes care of this by fitting a higher-order polynomial for
the continuum. Using these wavelength inclusion regions, we modelled the
telluric profile for each exposure in each time series independently. MAROON-X
observes two arms that overlap. The blue and red arms cover the wavelength
ranges between $490$ and $720\text{\,}\mathrm{n}\mathrm{m}$ and $640$ and
$920\text{\,}\mathrm{n}\mathrm{m}$, respectively. We treated the two arms
separately for our analysis due to different exposure times.

Earth’s atmospheric telluric lines were removed using molecfit (Smette et al.
2015; Kausch et al. 2015). It was applied to the one-dimensional spectra to
create a model for the telluric transmission spectrum of the entire wavelength
range. In the case of ESPRESSO, the reduction pipeline provides stitched and
blaze-corrected one-dimensional spectra, and all header information in these
files is constructed to be compatible with molecfit. For the MAROON-X
observations, we manually stitched and de-blazed the spectrum in every order
to create one-dimensional spectra for all the exposures. In addition, some of
the keywords are different from ESO standards and need to be specified in the
parameter file to run molecfit.333An example of the parameter file with
comments and a script that puts the data in the right format to be read by
tayph and later by molecfit in the framework of tayph can be found on GitHub
at https://github.com/bprinoth/MESOX.git.

The telluric model computed with molecfit was run separately for each exposure
in each time series. Regions with strong H2O and O2 absorption bands were used
to fit the model by selecting wavelength regions with telluric lines
surrounded by a flat continuum with no stellar lines. The inclusion regions
are shown in Fig. 3. The resulting telluric models were interpolated to the
same wavelength grid as the individual spectral orders and then divided out to
remove telluric absorption lines. Regions where the correction left visible
residuals were later masked out manually. We note here that due to the simple
de-blazing of the MAROON-X spectra, the continuum is not flat but has a slope
instead. This was taken care of by molecfit by fitting a higher-order
polynomial for the continuum. As shown in Fig. 3, there are not many telluric
absorption bands in the blue arm of MAROON-X; notably only a strong O2 band
around $630\text{\,}\mathrm{nm}$ ($\gamma$ band) and some weaker H2O bands
towards the end of the wavelength coverage of the arm can be seen.

###  3.2 Preparatory corrections

The individual spectra were Doppler-shifted to the rest frame of the host star
by correcting for the Earth’s velocity around the barycentre of the Solar
system as well as the radial velocity of the star caused by the orbiting
planet. This yielded a stellar spectrum with a constant velocity shift
consistent with the systemic velocity of $-24.452\pm 0.012$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. Following Hoeijmakers et al. (2020a),
we applied an order-by-order sigma clipping algorithm with a running median
absolute deviation computed over sub-bands of the time series with a width of
40 pixels and rejected $5\sigma$-outliers. We further manually flagged
spectral columns with visible systematic noise, where the telluric correction
of deep lines presented residuals, most notably around known O2 and H2O bands
at redder wavelengths. For the red arm of the MAROON-X observations, outlier
rejection as well as manual flagging affected, at most, 25.56% of the pixels
due to strong telluric bands being entirely masked out. For the other
observations, less than 2% of the spectral pixels were affected. For the
colour correction, we equalised the shape of the broadband continuum of each
exposure in the time series, order by order, by fitting a polynomial of degree
three.

To account for the induced pattern in the ESPRESSO observations that became
clearly visible when dividing the in-transit spectra by the out-of-transit
spectra (see Fig. 10), we introduced additional terms in the form of

$p_{0}\sin\left(\frac{2\pi x}{p_{1}}\right)+p_{2}\cos\left(\frac{2\pi
x}{p_{3}}\right),$ (1)

where $p_{i}$ are the free parameters and $x$ the pixels per order. Dividing
out the obtained fit removes any colour-dependent variations (polynomial) as
well as the low-frequency pattern (sinusoidal). The sine and cosine terms in
Eq. (1) are only used on the ESPRESSO time series. The improvement by the
inclusion of the sinusoidal fit was later shown to be insignificant, as the
high-pass filter at the end of the cross-correlation workflow removes most
residual variations (see Section 3.4). A comparison of the cross-correlation
results for Fe with and without the sinusoidal component is shown in Fig. 11.
We decided not to include this sinusoidal correction term when conducting our
analysis.

###  3.3 Cross-correlation templates

For atoms and ions, we used the cross-correlation templates provided by
Kitzmann et al. (2023), and for TiO, we produced a separate cross-correlation
template based on the Toto line list (McKemmish et al. 2019, states updated on
2021-08-25444 https://www.exomol.com/data/molecules/TiO/48Ti-16O/Toto/). For
the cross-correlation templates, the atmosphere was assumed to be in
isothermal and hydrostatic equilibrium, with temperatures of
$2500\text{\,}\mathrm{K}$ for TiO, $3000\text{\,}\mathrm{K}$ for atoms, and
$4000\text{\,}\mathrm{K}$ for ions. The temperature of
$4000\text{\,}\mathrm{K}$ is motivated by the fact that most atoms are not
significantly ionised below this temperature. For atoms, we used templates at
a temperature just above the equilibrium temperature of
$2641\text{\,}\mathrm{K}$ (Anderson et al. 2018). The templates and models
were broadened to the full width at half maximum of $2.7$, $2.14$, and
$3.53\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, matching the approximate
line-spread functions of HARPS/ HARPS-N, ESPRESSO, and MAROON-X, respectively.
A two-dimensional cross-correlation map was computed for each individual
species template as described in Prinoth et al. (2022). We used all the
available templates in Kitzmann et al. (2023) with sufficient absorption lines
at the considered temperatures, resulting in 90 cross-correlation templates
available for this study, as indicated in Fig. 12.

###  3.4 Cleaning steps

During transit the planet partially obscures areas of the rotating star disc,
introducing residual spectral lines when performing differential transmission
spectroscopy (Doppler shadow). In most cases, the Doppler shadow feature
occurs at a different apparent radial velocity than the planet signature such
that it can be removed without significantly affecting the planetary trace in
the cross-correlation map (for the case of the overlapping orbital trace and
Doppler shadow of MASCARA-1 b, see Casasayas-Barris et al. 2022). Residing on
a nearly polar orbit (Anderson et al. 2018), WASP-189 b introduces a Doppler
shadow feature that is nearly constant in velocity space, appearing vertically
in the cross-correlation map. It only overlaps with the planetary trace during
egress, at the end of the transit. The overlapping region was excluded when
co-adding the spectral time series to ensure that any residuals of the
correction of the Doppler shadow do not affect the results. Because of the
strong Fe absorption in the host star, we constructed empirical models for the
Doppler shadow features of Fe and Fe+ for atoms and ions respectively
(previously done in Prinoth et al. 2022) by fitting a function of the
following form to each row in the two-dimensional cross-correlation map (i.e.
each exposure $i$):

$\displaystyle\begin{split}f(v,i)&=A_{1}A_{\rm
poly}\exp\left(-\frac{(v-v_{\ast,i})^{2}}{2(W_{1}W_{\rm poly})^{2}}\right)\\\
&+A_{2}\exp\left(-\frac{(v-v_{\ast,i}-v_{\rm
offset})^{2}}{2W_{2}^{2}}\right)+C,\end{split}$ (2)

where $A_{i}$ are the constant amplitudes of the two Gaussians, $W_{i}$ are
the constant widths of the two Gaussians, and $C$ is a constant offset.
Furthermore, $A_{\rm poly}$ and $W_{\rm poly}$ describe two polynomials of
degree four to account for the centre-to-limb variations and the gravity
darkening of the host star visible in the Doppler shadow features. These
polynomials effectively allowed us to vary both the width and the amplitudes
of the first Gaussian component such that the bright emission-like feature at
the edges can also be removed successfully. The second Gaussian component
corrects the wider absorption-like feature and additionally offers the
possibility to be offset by a velocity shift of $v_{\rm
offset}\sim$46.55\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$, which
corresponds to half the projected rotational velocity of the star (see Table
2).

Table 2: Summary of stellar and planetary parameters of the WASP-189 system
adopted in this study.

Planetary parameters
---
Planet radius ($R_{\rm p}$) [$R_{\rm Jup}$] |  $1.600^{+0.017}_{-0.016}$ |  [3]
Planet mass ($M_{\rm p}$) [$M_{\rm Jup}$] |  $1.99^{+0.16}_{-0.14}$ |  [1]
Eq. temperature ($T_{\rm eq}$) [$\mathrm{K}$] |  $2641\pm 31$ |  [2]
Density ($\rho$) [$\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$] |  $0.6220^{+0.0769}_{-0.0365}$ |  [1]
Surface gravity ($\log g_{\rm p}$) |  $18.8^{+2.1}_{-1.8}$ |  [1]
Stellar parameters
Star radius ($R_{\ast}$) [$R_{\odot}$] |  $2.36\pm 0.030$ |  [1]
Star mass ($M_{\ast}$) [$M_{\odot}$] |  $2.030\pm 0.066$ |  [1]
Proj. rot. velocity ($v\sin{I_{\ast}}$) [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] |  $93.1\pm 1.7$ |  [1]
Systemic velocity ($v_{\rm sys}$) [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] |  $-24.452\pm 0.012$ |  [2]
Orbital and transit parameters
Transit centre time ($T_{0}$) [${\rm BJD}_{\rm TT}$ \- 2450000] |  $8926.5416960^{+0.0000650}_{-0.0000640}$ |  [1]
Orbital semi-major axis ($a$) [au] |  $0.05053\pm 0.00098$ |  [1]
Scaled semi-major axis ($a/R_{\ast}$) |  $4.600^{+0.031}_{-0.025}$ |  [1]
Orbital inclination ($i$) [∘] |  $84.03\pm 0.14$ |  [1]
Projected orbital obliquity ($\lambda$) [∘] |  $86.4^{+2.9}_{-4.4}$ |  [1]
Eclipse duration ($T_{14}$) [h] |  $4.3336^{+0.0054}_{-0.0058}$ |  [1]
Radius ratio ($R_{p}/R_{\ast}$) |  $0.07045^{+0.00013}_{-0.00015}$ |  [1]
RV semi-amplitude ($K$) [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] |  $0.182\pm 0.013$ |  [1]
Period ($P$) [d] |  $2.7240330\pm 0.0000042$ |  [2]
Eccentricity |  0.0 |  [1]
Derived parameters†
Orbital velocity $v_{\rm orb}$ [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] |  $201.8\pm 3.9$ |
Proj. orbital velocity ($K_{\rm p}=v_{\rm orb}\sin i$) [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] |  $200.7\pm 4.9$ |
Approx. scale height ($H$) [$\mathrm{km}$] |  $579^{+65}_{-56}$ |
Transit depth of $H$ ($\Delta F/F$) [$\times 10^{-5}$] |  $4.94^{+0.57}_{-0.49}$ |

Note: † The projected orbital velocity is given by $K_{\rm p}=v_{\rm orb}\sin
i=\frac{2\pi a\sin i}{P}$. The approximate scale height is given by
$H=\frac{RT_{\rm eq}}{\mu g_{\rm p}}$, where $R$ is the gas constant and $\mu$
is the mean molecular weight of the atmosphere. In the case of WASP-189 b, we
calculated the mean molecular weight using FastChem (Stock et al. 2018, 2022)
and found $\mu=$2.016\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$$, which
is in agreement with a mostly hydrogen and helium dominated atmosphere. The
transit depth of a scale height is calculated using $\frac{\Delta
F}{F}=\frac{A}{\pi R_{\ast}^{2}}=\frac{2HR_{\rm p}+H^{2}}{R_{\ast}^{2}}$. All
uncertainties were calculated assuming Gaussian error propagation. References:
[1] Lendl et al. (2020), [2] Anderson et al. (2018), [3] Deline et al. (2022).

The local radial velocity $v_{\ast}$ in Eq. (2) is given by the parameters of
the system, such as the semi-major axis $a$, the stellar radius $R_{\ast}$,
the phase of the transiting planet on its orbit $\phi$, the projected
rotational velocity of the star $v\sin I_{\ast}$, the projected orbital
obliquity $\lambda$, and the inclination of the system $i$ as stated in Table
2. It takes the following form (Cegla et al. 2016)

$\mathbf{v_{\ast}(\phi)}=\frac{a}{R_{\ast}}\left(\sin
2\pi\mathbf{\phi}\cos\lambda+\cos 2\pi\mathbf{\phi}\cos
i\sin\lambda\right)v\sin I_{\ast}.$ (3)

We scaled the empirical models for Fe and Fe+ to fit the Doppler shadow of the
individual species and divided out of the time-resolved two-dimensional cross-
correlation map. We expected the projected rotational broadening of the
tidally locked planetary atmosphere to be

$\sigma_{\rm rot}=v_{\rm rot}\sin i=\frac{2\pi R_{p}}{P}\sin
i=3.04\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ (4)

using the parameters in Table 2. While fitting the models, the planetary
signature was protected by masking out the radial velocity range of the planet
at each orbital phase, allowing for a width of
$16\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$of the mask to also ensure
that broadened lines and potential winds of the order of 10
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ of the planetary atmosphere do not
contribute to the fit.

After the removal of the Doppler shadow, we performed a Gaussian high-pass
filter with a width of $100\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$
that removes broadband structures present in the spectral direction, which
also removes the pseudo-sinusoidal variations introduced by the reflections in
the Coudé train visible in the ESPRESSO data. All the correlated structures
and systematic noise constant in radial velocity that remained present after
applying the Gaussian high-pass filter are not related to the planetary
atmospheres and could thus be removed through detrending of the spectral
channels individually following (Prinoth et al. 2022; Borsato et al. 2023).
These structures are mainly caused by aliases between strong lines of the
cross-correlation templates and the stellar absorption lines consecutively
obscured by the transiting planet. They appear at constant velocities, thus
vertically in the cross-correlation maps, and are again caused by the nearly
polar alignment of the planetary orbit.

###  3.5 Model injection

In order to predict the expected signal strength, we injected models into the
data before performing the cross-correlation analysis. We followed the same
procedure as in Prinoth et al. (2022) and used the same models at a
temperature of $2500\text{\,}\mathrm{K}$, close to the planetary equilibrium
temperature of $2641\text{\,}\mathrm{K}$ (Anderson et al. 2018), and
$3000\text{\,}\mathrm{K}$, close to the dayside temperature (Lendl et al.
2020), assuming an isothermal atmosphere in chemical and hydrostatic
equilibrium of solar metallicity. The models include the continuum and
accurate line depths and profiles. Prior to cross correlation, we injected
model spectra into the observed spectra. When subtracting the true cross-
correlation map (without injection) from the one with the injected signal, the
Doppler shadow, aliases, and correlated noise were cancelled. This effectively
left a residual signature signifying the predicted line depth and shape of the
models.

Figure 4: Cross-correlation results in the $K_{p}-V_{\rm sys}$ space for the
species detected in WASP-189 b’s atmosphere during eight transit events with
HARPS, HARPS-N, ESPRESSO, and MAROON-X combined. Each panel shows the
significance of the detection for a given species, and the solid lines
indicate the expected location of the signal if we assume a symmetric planet
with a static atmosphere. Deviations from this location can be indicative of
chemical asymmetries and dynamical effects. The lower-right panel shows the
cross-correlation result in $K_{p}-V_{\rm sys}$ space for a template
containing all species at a temperature of $3000\text{\,}\mathrm{K}$. As
previously seen in Prinoth et al. (2022), the atmospheric signal is located at
a systemic and projected orbital velocity smaller than the true values,
indicating a global dayside-to-nightside flow.  Figure 5: One-dimensional
cross-correlation results extracted at the projected best-fit orbital velocity
compared to the injected signal at $K_{\rm p}=$ $200.7\pm 4.9$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The model predictions show the
expected signal strength using models at $2500\text{\,}\mathrm{K}$ (in teal)
and $3000\text{\,}\mathrm{K}$ (in magenta) in chemical and hydrostatic
equilibrium of solar metallicity, previously used in Prinoth et al. (2022).
The shaded regions indicate the standard deviation away from the peak. Neutral
atoms are generally well described by the model at $2500\text{\,}\mathrm{K}$.
The dashed vertical line shows the expected systemic velocity of $-24.452\pm
0.012$ $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The model prediction of TiO
has been scaled by a factor of six to match the predicted signal strength in
Prinoth et al. (2022) for the purpose of comparing the datasets from the
instruments in Fig. 6.

###  3.6 Velocity-velocity maps

We then converted the resulting cleaned two-dimensional cross-correlation maps
into velocity-velocity maps ($K_{p}-V_{\rm sys}$ diagrams) by shifting them
towards the expected rest frame of the planet ($K_{\rm p}=$ $200.7\pm 4.9$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, $v_{\rm sys}=$ $-24.452\pm
0.012$$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$), assuming values of projected
orbital velocity between $0$ and
$400\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ in steps of
$1\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. In each step, the exposures
were weighted according to the mean flux of the corresponding spectra. The
$K_{p}-V_{\rm sys}$ diagrams of the different observations were then combined
by stacking and weighting based on the injected signals, where predicted (see
Section 3.5), or based on the expected flux in the covered wavelength regions
predicted by the template, if the species was not predicted by model
injection. In the case of weighing by injection, this accounts for observing
conditions, wavelength coverage, and signal-to-noise properties of the
spectrographs. If weighted by the flux from the templates, we would only
account for the different wavelength coverage.

###  3.7 Time-resolved spectroscopy

The $K_{p}-V_{\rm sys}$ diagrams represent an average over the entire limb of
the planet, though we know the signal will vary at each location on the limb,
so it is challenging to identify the origin and nature of the contribution of
the different regions of the atmosphere to the observed absorption signature
(Wardenier et al. 2021). With the high signal-to-noise from the combined eight
individual transit time series, it is possible to fit the planetary trace in
the two-dimensional cross-correlation functions by means of time-resolved
spectroscopy. As these cross-correlation functions are provided in the rest
frame of the star, the planetary atmosphere traces a nearly linear part of the
sinusoid of the radial velocity curve given as

$\bm{v_{\rm planet}}(\bm{\phi})=v_{\rm
orb}\sin\left(2\pi\bm{\phi}\right)\sin{i}+v_{\rm sys},$ (5)

where $v_{\rm orb}$ is the orbital velocity of the planet, $\bm{\phi}$ is the
orbital phase, $i$ is the inclination of the system, and $v_{\rm sys}$ is the
systemic velocity. We modelled the two-dimensional absorption feature in
transit as

$\displaystyle\begin{split}{\rm\mathbf{CCF}}&=\mathbf{A}(\bm{\phi},F_{\rm
start},F_{\rm end})\cdot\exp{\left(-\frac{(\bm{v}-\bm{v_{\rm
planet,1}})^{2}}{2\sigma^{2}}\right)}\cdot\bm{f}\\\
&+\mathbf{A}(\bm{\phi},F_{\rm start},F_{\rm
end})\cdot\exp{\left(-\frac{(\bm{v}-\bm{v_{\rm
planet,2}})^{2}}{2\sigma^{2}}\right)}\cdot\bm{f_{\rm inv}}+C,\end{split}$ (6)

with $\bm{v_{\rm planet,i}}$ as defined in Eq. (5) for the two halves of the
transit independently. We allowed the orbital velocity ${v_{\rm orb}}$ to vary
between the two halves of the transit but required the systemic velocity
$v_{\rm sys}$ to be the same for both halves, and thus it was fit at the
centre of the transit. The terms $f$ and $f_{\rm inv}$ approximate two
Heaviside step functions defined as

$\displaystyle\bm{f}$ $\displaystyle=\frac{1}{2}\tanh{k\bm{\phi}}+\frac{1}{2}$
(7) $\displaystyle\bm{f_{\rm inv}}$
$\displaystyle=-\frac{1}{2}\tanh{k\bm{\phi}}+\frac{1}{2}$ (8)

which divide the transit into the first ($\bm{f}$) and second half
($\bm{f_{\rm inv}}$). We chose the $\tanh$ prescription of the Heaviside step
function because our framework requires the model to be fully auto-
differentiable to work with NumPyro and Jax (Bradbury et al. 2018; Bingham et
al. 2018; Phan et al. 2019). To account for the steep cut at the centre of the
transit, we chose $k=1000$. The amplitude $\mathbf{A}(\bm{\phi},F_{\rm
start},F_{\rm end})$ is given as a function of the phases of the transit, and
the cross-correlation signal strength are indicated with $F_{\rm start}$ and
$F_{\rm end}$ at the start and the end of the transit, respectively. This is
given as

$\mathbf{A}(\bm{\phi},F_{\rm start},F_{\rm end})=\frac{F_{\rm start}\phi_{\rm
end}-F_{\rm end}\phi_{\rm start}}{\phi_{\rm end}-\phi_{\rm
start}}-\frac{F_{\rm end}-F_{\rm start}}{\phi_{\rm end}-\phi_{\rm
start}}\bm{\phi}.$ (9)

The ratio $\frac{F_{\rm start}}{F_{\rm end}}$ provides a first-order estimate
for the increase ($>1$) or decrease ($<1$) in signal strength throughout the
transit. The model in Eq. (6) has the following free parameters and priors:

  * •

The signal strength at the start of the transit $F_{\rm
start}\sim\mathcal{U}(0,1000)$ in parts per ten thousand, where $\phi_{\rm
start}=-0.034=0.966$.

  * •

The signal strength at the end of the transit $F_{\rm
end}\sim\mathcal{U}(0,1000)$ in parts per ten thousand, where $\phi_{\rm
end}=0.034$

  * •

A constant offset $C\sim\mathcal{U}(-3,3)$ in ppm.

  * •

The projected orbital velocities $K_{\textrm{p, i}}\sim\mathcal{U}(50,350)$ in
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, roughly centred on the expected
projected orbital velocity of $200.7\pm 4.9$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$.

  * •

The systemic velocity $v_{\textrm{sys}}\sim\mathcal{U}(-35,-15)$ in
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, roughly centred on the known systemic
velocity of $-24.452\pm 0.012$$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$.

  * •

The Gaussian line width $\sigma_{\rm w}\sim\mathcal{U}(1,14)$ in
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The expected rotational broadening
due to the planetary co-rotation is
$\approx$$3\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, as given by
$v_{\rm rot}=\frac{2\pi R_{p}\sin i}{P}$ using the parameters stated in Table
2.

We sampled from these prior distributions and evaluated the likelihood in a
Bayesian framework using a No-U-Turn Sampler (see Betancourt 2017, for a
review). Our implementation of this framework is based on NumPyro and Jax
(Bradbury et al. 2018; Bingham et al. 2018; Phan et al. 2019), as previously
done for emission in Hoeijmakers et al. (2022). We chose 200 warm-up samples
and 600 samples over 20 chains running in parallel. When fitting the trace, we
ignored the contributions from the region of the Doppler shadow entirely to
avoid biasing the results by potential over- or under-correction of the
stellar signal.

Apart from a first-order estimate of the change in signal strength through the
ratio $\frac{F_{\rm start}}{F_{\rm end}}$, the free parameters describe the
shape and location of the trace, $v_{\textrm{sys}}$ describes the velocity of
the atmosphere at mid-transit, and a deviation from the expected systemic
velocity is indicative of the location of the species in the atmosphere and
dynamics. We allowed different orbital velocities in the two halves of the
transit (see e.g. Ehrenreich et al. 2020; Kesseli & Snellen 2021; Kesseli et
al. 2022, for the shape of Fe feature of WASP-76 b). Lastly, $C$ describes the
continuum of the two-dimensional cross-correlation function.

##  4 Results and discussion

###  4.1 Cross-correlation analysis

Figure 4 shows the 17 detected chemical species in the transmission spectrum
of WASP-189 b after combining eight transit time series obtained with HARPS,
HARPS-N, ESPRESSO, and MAROON-X. We detected H, Na, Mg, Ca, Ca+, Ti, Ti+, TiO,
V, Cr, Mn, Fe, Fe+, Ni, Sr, Sr+, and Ba+. The detection significances
including uncertainties are listed in Table 3. While Prinoth et al. (2022)
tentatively detected Cr+ and Sc+ at a velocity consistent with the projected
orbital velocity of the planet, we were not able to confirm this detection
when including three more transit observations. Apart from Cr+ and Sc+, we
confirm all detections and tentative detections in Prinoth et al. (2022) and
Stangret et al. (2022). Additionally, we detected Sr, Sr+, and Ba+, which have
not been detected in the atmosphere of this planet previously. We report
detections for both Sr and Sr+, but we note that the signal shapes are
irregular (more broadened) in comparison to the rest of the detections. This
requires further investigation.

###  4.2 Model comparison

Figure 5 shows the expected significance when injecting the models introduced
in Section 3.5. The absorption lines of most neutral atoms are broadly
consistent with a model at a temperature of $2500\text{\,}\mathrm{K}$ in local
thermodynamic equilibrium (LTE) and hydrostatic and chemical equilibrium. This
agreement was already observed in Prinoth et al. (2022). As previously
observed in WASP-121 b (Hoeijmakers et al. 2020a), the signal strength for Ni
is higher than predicted by any of the two models. The strong absorptions of
Fe+ and Ca+ are inconsistent with this class of models, where Fe+ is not
predicted to be observable, while Ca+ is predicted to be weaker, even at a
model at $3000\text{\,}\mathrm{K}$. Departures from model predictions,
especially in the case of ions, are likely reflective of the unknown
temperature-pressure profiles at high altitudes, the possible breakdown of
hydrostatic equilibrium for the deepest lines, and non-LTE effects (see also
Prinoth et al. 2022; Ehrenreich et al. 2020; Kesseli et al. 2022; Savel et al.
2022). The broad shape of the Ca+ feature is consistent with absorption at
high altitudes as part of an extended outflow (Yan et al. 2019; Seidel et al.
2021; Tabernero et al. 2021; Maguire et al. 2023; Pelletier et al. 2023). We
found that Ba+ is consistent with the model prediction at
$2500\text{\,}\mathrm{K}$ and Sr+ at $3000\text{\,}\mathrm{K}$, while Ti+
requires a model at a higher temperature to match the observed signal. In
general, the predictions for ionised species show that this class of models
does not reliably predict the absorption depths observed in ultra-hot
Jupiters.

###  4.3 Robust confirmation of TiO

Figure 6: Cross-correlation results in $K_{p}-V_{\rm sys}$ space for TiO in
WASP-189 b’s atmosphere during eight transit events with HARPS & HARPS-N,
ESPRESSO, and MAROON-X. From left to right: MAROON-X alone, MAROON-X
restricted to ESPRESSO wavelength coverage (wlc), and ESPRESSO alone (see
Table 1). The dashed lines indicate the expected location of the signal if we
assume a symmetric planet with a static atmosphere. The model predictions show
the expected signal strength using models at $2500\text{\,}\mathrm{K}$ (in
teal) and $3000\text{\,}\mathrm{K}$ (in magenta) in chemical and hydrostatic
equilibrium of solar metallicity, previously used in Prinoth et al. (2022),
and have been scaled by a factor of six for clarity to show that the model
discrepancy is equivalent in the data from the three different instruments
despite their different wavelength coverage. The predicted signals are a
measure of the sensitivity of the dataset and show that we expect MAROON-X to
have a stronger signal. We note that discrepancies with Prinoth et al. (2022)
stem from a different weighting algorithm based on the predicted signal
strength (through model injection) applied here. This weighting algorithm
favours the signal of better nights, increasing the signal strength overall.
The results show that the signal is recovered in the new datasets alone as
well but with a lower relative flux compared with expectations from the model
injection. Figure 7: Change in signal strength between the start and the end
of the transit as a function of the observed deviation from the expected
systemic velocity. The error bars denote 1-$\sigma$ uncertainties calculated
via Gaussian error propagation. The grey vertical region indicates the true
systemic velocity range ($-24.452\pm
0.012$$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$). We note that Sr and Sr+ are
not included in this figure due to large uncertainties. We also note that Fe+
is an outlier (see Fig. 9). The hotter atmosphere rotating into view (see Fig.
1) causes most signal strengths to increase over the course of transit.
However, further ionisation (Ca+, Ni, Ba+, and H) or dissociation (TiO) on the
dayside can also cause some signals to decrease over time.

The combination of the two MAROON-X and the ESPRESSO observations robustly
confirms the detection for TiO in the transmission spectrum of WASP-189 b.
Combining the data of all eight transit time series, we detected TiO at 6.65
$\sigma$. The observed significance is in agreement with the model at
$2500\text{\,}\mathrm{K}$ if scaled by a factor of six (needed to bring the
model prediction in agreement with Prinoth et al. 2022). Despite the
additional noise caused by interference patterns in the ESPRESSO data in the
form of internal reflections in the optical path, it is possible to detect and
thus confirm TiO (see Fig. 6). The detection of TiO in the MAROON-X and
ESPRESSO data combined is mostly driven by the MAROON-X data.

Figure 8: Ionisation fraction and energies for the atoms detected in the
atmosphere of WASP-189 b. The ionisation curves (top panel) were computed
using FastChem, assuming chemical equilibrium at a pressure of 1 mbar and
solar metallicity. Group 1 elements (hydrogen + alkali metals) are plotted in
blue, and group 2 elements (alkaline earth metals) are plotted in orange.
Transition metals are plotted in teal. The dashed, dash-dotted, and dotted
lines indicate that the detected signal strength of the neutral atom, ionised
atom, and oxide weakens over the course of the transit, respectively. Solid
lines indicate that the signal strengths of the detections either strengthen
over the course of the transit or remain constant, and thus behave as
expected. The ionisation curves move towards lower temperatures with
decreasing pressure (increasing altitude). Additionally, the fractions
decrease at high temperatures for alkaline earth metals due to a second
ionisation. The ionisation energies (bottom panel) are from the NIST database
(Kramida et al. 2022) which predict that Na, Ba, Ca, Ba+, and Sr are all
easily ionised. Regarding H, it may significantly ionise at very high
altitudes. Because of the constant effective nuclear charge within groups 1
and 2, the increasing orbital radius causes the ionisation energy to decrease
for higher atomic numbers $Z$ (within the group). Within the same period
(row), the ionisation energy increases due to the stronger attractive force of
the nucleus. Both trends (within period and group) are generally not true for
transition metals due to electron shells not necessarily being filled (and
ionised) consecutively (Froese-Fischer et al. 1997).

###  4.4 Time-resolved spectroscopy

Combining the eight transit time series enabled the investigation of the time
dependence of signals through time-resolved spectroscopy, as described in
Section 3.7. Figures 14–29 show the results of the orbital trace fitting, and
the best-fit parameters are summarised in Table 3.

We observed changes in the signal strength of nine of the detected species, as
shown in Fig. 7. The atmosphere was expected to be colder and less inflated at
the leading terminator (east side), whereas both temperature and scale height
increased when at the trailing terminator (west side), see Fig. 1. The
observations towards the end of the transit are therefore dominated by the
light that has passed through the trailing terminator, so an increase in
signal strength is consistent with the increased scale height there. With the
exception of H, Ba+, Ca+, Ni, and TiO, all detected species either showed an
increase in signal strength over the course of the transit or appeared to be
constant.

A decrease in signal strength may be caused by a loss of the absorbing species
towards the hot trailing terminator. We propose that this can be explained
through either ionisation or dissociation on the dayside. Figure 8 shows the
ionisation temperatures and energies for the detected atomic species for the
first and second ionisation states. The ionisation curves were computed using
FastChem (Stock et al. 2018, 2022), assuming chemical equilibrium. The
ionisation energies were taken from the NIST database (Kramida et al. 2022).
The temperature on the dayside is likely high enough to thermally ionise
species with relatively lower ionisation energies. Particularly notable are
Ba+ and Ni, which show a decrease in signal strength as hotter regions rotate
into view, where they are more susceptible to ionisation (secondary ionisation
in the case of Ba+). We infer that the absorption signal of Ba+ predominantly
originates from the hot region of the leading terminator (region 4 in Fig. 1),
which rotates out of view over the course of the transit. The same argument
may also hold for Ni, as it is expected to ionise on the dayside through to
Ni+, and for H and Ca+, as they are easily ionised to H+ and Ca^2+,
respectively, at higher temperatures. As Ni+ predominantly absorbs at
wavelengths below the coverage of ESPRESSO (bluest wavelength range of the
spectrographs used), a search for it is currently limited by the blue cut-off
of the spectrograph. The signal strength of TiO decreases over the course of
the transit, while Ti+ increases and Ti appears to be constant. We hypothesise
that TiO dissociates on the dayside (Cont et al. 2021), while Ti partly
ionises to Ti+, leaving the signal strength of Ti roughly constant. At the
same time, molecular dissociation at high temperatures may also act to
replenish the neutral reservoir. For example, Ca is one of the most refractory
elements (Lodders 2003), though it also has a relatively low ionisation
temperature. This means that the signal strength of neutral calcium on the
opposite limbs results from a balance between depletion by ionisation and
creation from molecular dissociation. More advanced chemical models including
condensates may be needed to self-consistently explain these time-dependent
signals.

Azevedo Silva et al. (2022) recently reported the detection of Ba+ in the
transmission spectrum of WASP-121 b ($T_{\rm eq}=$2358\pm
52.\text{\,}\mathrm{K}$$, Delrez et al. 2016)) and WASP-76 b ($T_{\rm
eq}=$2160\pm 40.\text{\,}\mathrm{K}$$, West et al. 2016), and Borsato et al.
(2023) detected it in the transmission spectrum of KELT-9 b ($T_{\rm
eq}=$3921\pm 180.\text{\,}\mathrm{K}$$, Borsa et al. 2019). We fit the
atmospheric trace of Ba+ obtained by Borsato et al. (2023) to show that this
signal also weakens over the course of the transit (see Fig. 30).

Figure 7 shows the change in signal strength between the start and the end of
the transit as a function of the observed deviation from the true systemic
velocity of $-24.452\pm 0.012$$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ for each
individual species. With the exception of Fe+, all of the detected species are
blueshifted with respect to the true systemic velocity of $-24.452\pm 0.012$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, indicative of dayside-to-nightside
winds and flows at high altitudes. The orbital fit for Fe+ suggests that the
measured systemic velocity at the centre of the transit is redshifted with
respect to the true systemic velocity. This is likely not astrophysical. Our
fit works well for species that do have a roughly constant orbital velocity
($v_{\rm orb,1}\sim v_{\rm orb,2}$) or where the change between the two
velocities is at the centre of the transit. Figure 9 shows that for Fe+, the
change in orbital velocity should happen earlier. The trace of the Fe+ feature
is thus inconsistent with the best-fit model, as the model does not allow for
change at other times than mid-transit.

###  4.5 Future prospects

The strongest signals are exhibited by Ca+ and Ba+, inviting further analysis
using narrow-band spectroscopy (e.g. Seidel et al. 2023). Notably, narrow-band
spectroscopy has been successfully used in investigations of the Na doublet
(e.g. Casasayas-Barris et al. 2018; Seidel et al. 2019, 2021, 2023) and will
be further discussed in Prinoth et al. (in prep.). As for Na, Ca+, and Ba+,
they are expected to behave similarly, with only a few very strong lines in
their transmission spectra.

##  5 Conclusions

In this study, we searched for 90 different species, including atoms, ions,
and various oxides, in the transmission spectrum of WASP-189 b using the
cross-correlation technique on eight transit time series with HARPS, HARPS-N,
ESPRESSO, and MAROON-X. We detected H, Na, Mg, Ca, Ca+, Ti, Ti+, TiO, V, Cr,
Mn, Fe, Fe+, Ni, Sr, Sr+, and Ba+ at high significance. We further confirmed
the detection for TiO and added novel detections for Sr, Sr+, and Ba+. We note
that Sr and Sr+ have a deviating line shape that requires further
investigation.

For each of the detected species, we fit the orbital trace by means of a
Bayesian framework to infer posterior distributions for orbital parameters as
well as the line depth and shape as a function of time. For all species except
Fe+, the fit resulted in a blueshifted signal with respect to the true
systemic velocity. We determined that the apparent redshift for Fe+ is due to
the change in orbital velocity that is offset from the centre of the transit.
While the signal strength of the majority of the species increases over the
course of the transit, as expected when probing the hotter trailing terminator
(increased scale height), the signal strength of Ni, Ba+, Ca+, H, and TiO
decreases instead. Throughout the course of the transit, as hotter regions of
the atmosphere rotate into view, the signal of certain species weakens. We
posit that this is likely caused by ionisation through incoming stellar
radiation on the dayside that causes the atmosphere to heat up significantly
(Ni, Ba+, Ca+, H). For TiO, the weakened signal is likely similarly caused by
dissociation on the hot dayside.

Time-resolved spectroscopy offers a new possibility to study the atmospheres
of hot gas giant exoplanets. By combining several transit observations, the
signal-to-noise ratio becomes sufficiently high to study the absorption of the
atmosphere resolved in the orbital phase, and thus, it is no longer necessary
to average the observations in the rest frame of the planet in order to detect
absorbing species. Instead of averaging over the contributions of the two
terminators, time-resolved spectroscopy enables studying the change in signal
strength and dynamics over the course of the transit as different regions of
the tidally locked atmosphere rotate into view. It is especially useful for
short-orbit planets with plenty of available observations, preferably those
with different spectrographs so that a wide wavelength range can be covered
(e.g. WASP-76 b, WASP-121 b, KELT-9 b, and KELT-20 b.)

Ultimately, these observations need to be interpreted using models that
accurately describe dissociation and ionisation processes (whether thermal or
driven by the radiation field) as well as nightside condensation. Explaining
this three-dimensional complexity of ultra-hot Jupiters will likely require
models that simulate global circulation as well as radiative transfer in three
dimensions.

Figure 9: Two-dimensional cross-correlation function of Fe+ for all eight
transit time series combined. The cross-correlation function is shifted to the
projected orbital velocity found for the second part of the transit, $v_{\rm
orb,2}$ to illustrate the earlier change of orbital velocity. We purposely
limited the colour bar in order to enhance the absorption feature that helps
make the misfit visible. The dot indicates the centre of the transit. The blue
shaded region shows the best-fit result. Table 3: Overview of detections.

|  $F_{\rm start}$ [$\times 10^{-4}$] |  $F_{\rm end}$ [$\times 10^{-4}$] |  $A$ [$\times 10^{-4}$] |  $\sigma_{\rm lw}$ [$\times 10^{-4}$] |  $v_{\rm sys}$ [km s-1] |  $v_{\rm orb,1}$ [km s-1] |  $v_{\rm orb,2}$ [km s-1] |  $C$ [$\times 10^{-6}$] | $\sigma$
---|---|---|---|---|---|---|---|---|---
H |  2604.71 $\pm$ 98.53 |  835.42 $\pm$ 106.81 |  1720.06 $\pm$ 72.66 |  11.19 $\pm$ 0.37 |  -25.29 $\pm$ 0.70 |  197 $\pm$ 12.96 |  157.49 $\pm$ 6.33 |  0.14 $\pm$ 1.24 | 23.67
Na |  169.76 $\pm$ 27.27 |  195.33 $\pm$ 26.16 |  182.54 $\pm$ 18.89 |  6.68 $\pm$ 0.62 |  -27.69 $\pm$ 0.64 |  145.07 $\pm$ 7.76 |  145.07 $\pm$ 7.76 |  0.16 $\pm$ 0.33 | 9.66
Mg |  257.48 $\pm$ 40.41 |  647.45 $\pm$ 47.43 |  452.47 $\pm$ 31.15 |  10.86 $\pm$ 0.64 |  -28.85 $\pm$ 0.79 |  109.18 $\pm$ 7.79 |  203.58 $\pm$ 13.05 |  0.22 $\pm$ 0.53 | 14.52
Ca |  116.96 $\pm$ 17.99 |  146.09 $\pm$ 18.25 |  131.52 $\pm$ 12.82 |  7.43 $\pm$ 0.69 |  -29.55 $\pm$ 1.3 |  148.59 $\pm$ 13.57 |  178.49 $\pm$ 15.63 |  0.02 $\pm$ 0.22 | 10.26
Ca+ |  785.0 $\pm$ 22.61 |  518.62 $\pm$ 25.39 |  651.81 $\pm$ 17.0 |  10.58 $\pm$ 0.23 |  -30.42 $\pm$ 0.42 |  182.27 $\pm$ 4.04 |  184.11 $\pm$ 7.85 |  -0.34 $\pm$ 0.33 | 38.34
Ti |  64.92 $\pm$ 8.88 |  70.53 $\pm$ 11.66 |  67.73 $\pm$ 7.33 |  6.63 $\pm$ 0.57 |  -27.48 $\pm$ 0.93 |  194.52 $\pm$ 12.37 |  215.96 $\pm$ 15.41 |  -0.02 $\pm$ 0.1 | 9.24
Ti+ |  182.29 $\pm$ 34.38 |  465.17 $\pm$ 46.38 |  323.73 $\pm$ 28.87 |  8.43 $\pm$ 0.64 |  -27.7 $\pm$ 1.16 |  146.65 $\pm$ 15.4 |  225.35 $\pm$ 16.45 |  0.05 $\pm$ 0.43 | 11.21
TiO |  4.17 $\pm$ 0.74 |  2.45 $\pm$ 0.67 |  3.31 $\pm$ 0.5 |  4.53 $\pm$ 0.64 |  -31.21 $\pm$ 1.38 |  145.75 $\pm$ 12.61 |  208.8 $\pm$ 16.44 |  -0.0 $\pm$ 0.01 | 6.65
V |  70.8 $\pm$ 11.49 |  69.66 $\pm$ 13.05 |  70.23 $\pm$ 8.69 |  9.21 $\pm$ 0.82 |  -28.74 $\pm$ 1.04 |  207.5 $\pm$ 10.17 |  207.5 $\pm$ 10.17 |  -0.06 $\pm$ 0.14 | 8.08
Cr |  96.51 $\pm$ 12.84 |  97.45 $\pm$ 13.22 |  96.98 $\pm$ 9.22 |  7.52 $\pm$ 0.7 |  -29.05 $\pm$ 1.12 |  165.65 $\pm$ 12.05 |  177.93 $\pm$ 14.58 |  0.05 $\pm$ 0.15 | 10.52
Mn |  180.66 $\pm$ 32.36 |  229.53 $\pm$ 46.52 |  205.1 $\pm$ 28.33 |  8.27 $\pm$ 1.02 |  -31.5 $\pm$ 1.45 |  159.33 $\pm$ 17.61 |  251.63 $\pm$ 28.15 |  0.11 $\pm$ 0.41 | 7.24
Fe |  122.19 $\pm$ 11.02 |  247.64 $\pm$ 12.0 |  184.91 $\pm$ 8.14 |  7.62 $\pm$ 0.29 |  -28.76 $\pm$ 0.52 |  171.65 $\pm$ 6.92 |  195.83 $\pm$ 6.67 |  0.06 $\pm$ 0.13 | 22.7
Fe+ |  454.94 $\pm$ 35.28 |  739.59 $\pm$ 47.25 |  597.26 $\pm$ 29.49 |  8.26 $\pm$ 0.38 |  -23.15 $\pm$ 0.55 |  218.0 $\pm$ 7.09 |  177.55 $\pm$ 7.73 |  0.14 $\pm$ 0.44 | 20.26
Ni |  136.63 $\pm$ 22.26 |  56.0 $\pm$ 28.64 |  96.32 $\pm$ 18.14 |  5.56 $\pm$ 0.65 |  -27.13 $\pm$ 1.61 |  180.04 $\pm$ 13.59 |  222.06 $\pm$ 49.0 |  -0.13 $\pm$ 0.23 | 5.31
Sr |  39.52 $\pm$ 31.54 |  640.32 $\pm$ 55.83 |  339.92 $\pm$ 32.06 |  12.49 $\pm$ 0.8 |  -24.35 $\pm$ 1.04 |  233.02 $\pm$ 11.85 |  233.02 $\pm$ 11.85 |  0.18 $\pm$ 0.61 | 10.6
Sr+ |  172.74 $\pm$ 86.49 |  1432.13 $\pm$ 136.44 |  802.44 $\pm$ 80.77 |  11.32 $\pm$ 0.82 |  -24.12 $\pm$ 0.94 |  229.06 $\pm$ 10.13 |  229.06 $\pm$ 10.13 |  0.27 $\pm$ 1.24 | 9.93
Ba+ |  853.22 $\pm$ 89.18 |  392.28 $\pm$ 103.35 |  622.75 $\pm$ 68.25 |  7.16 $\pm$ 0.53 |  -25.42 $\pm$ 1.14 |  175.01 $\pm$ 9.89 |  182.69 $\pm$ 18.58 |  -0.12 $\pm$ 0.97 | 9.12

Note: The average signal strength $A$ is the mean of the signal strength at
the start and end of the transit, $F_{\rm start}$ and $F_{\rm end}$. We note
that for Na, V, Sr, and Sr+, the orbital velocity was forced to be the same
for both regions of the transit because the posterior distribution reached the
boundaries of the feasible priors. The detection significance $\sigma$ was
calculated from the average signal strength $A$ and its uncertainty.

###### Acknowledgements.

B.P., H.J.H., and N.W.B acknowledge partial financial support from The Fund of
the Walter Gyllenberg Foundation. S.P. and B.B. acknowledge financial support
from the Natural Sciences and Engineering Research Council (NSERC) of Canada
and the Fond de Recherche Québécois-Nature et Technologie (FRQNT; Québec).
M.Br. acknowledges support from the National Science Foundation Graduate
Research Fellowship under Grant No. DGE 1746045. R.L. acknowledges funding
from the University of La Laguna through the Margarita Salas Fellowship from
the Spanish Ministry of Universities ref. UNI/551/2021-May 26, and under the
EU Next Generation funds. B.T. acknowledges the financial support from the
Wenner-Gren Foundation (WGF2022-0041). M.Bu. acknowledges support from the
Swedish Research Council VR through grant no. 2016-04185, and the Crafoord
Foundation. We thank E.K.H Lee for providing us with the temperature map of
WASP-189 b at a pressure of 1 mbar published in Lee et al. (2022). We thank
Per Jönsson for his support in understanding the screening effect. We thank
Romain Allart for discussions about the interference patterns in the ESPRESSO
data. This study makes use of astropy (Astropy Collaboration et al. 2013,
2018, 2022), NumPyro (Bingham et al. 2018),  matplotlib-label-lines (Cadiou
2022), corner.py (Foreman-Mackey 2016), ArviZ (Kumar et al. 2019) and JAX
(Bradbury et al. 2018). This work is based in part on observations collected
at the European Southern Observatory under ESO programme 107.22QF. This
research has made use of the services of the ESO Science Archive Facility.
This work is based in part on data obtained at the international Gemini
Observatory, a programme of NSF’s NOIRLab. The international Gemini
Observatory at NOIRLab is managed by the Association of Universities for
Research in Astronomy (AURA) under a cooperative agreement with the National
Science Foundation on behalf of the Gemini partnership: the National Science
Foundation (United States), the National Research Council (Canada), Agencia
Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia,
Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia,
Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science
Institute (Republic of Korea). The MAROON-X team acknowledges funding from the
David and Lucile Packard Foundation, the Heising-Simons Foundation, the Gordon
and Betty Moore Foundation, the Gemini Observatory, the NSF (award number
2108465), and NASA (grant number 80NSSC22K0117). We thank the staff of the
Gemini Observatory for their assistance with the commissioning and operation
of the instrument. This work was enabled by observations made from the Gemini
North telescope, located within the Maunakea Science Reserve and adjacent to
the summit of Maunakea. We are grateful for the privilege of observing the
Universe from a place that is unique in both its astronomical quality and its
cultural significance.

## References

  * Allart et al. (2020) Allart, R., Pino, L., Lovis, C., et al. 2020, A&A, 644, A155
  * Anderson et al. (2018) Anderson, D. R., Temple, L. Y., Nielsen, L. D., et al. 2018, arXiv:1809.04897
  * Arcangeli et al. (2018) Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30
  * Arras & Socrates (2010) Arras, P. & Socrates, A. 2010, ApJ, 714, 1
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Azevedo Silva et al. (2022) Azevedo Silva, T., Demangeon, O. D. S., Santos, N. C., et al. 2022, A&A, 666, L10
  * Bell & Cowan (2018) Bell, T. J. & Cowan, N. B. 2018, ApJL, 857, L20
  * Beltz et al. (2022) Beltz, H., Rauscher, E., Kempton, E. M.-R., et al. 2022, AJ, 164, 140
  * Beltz et al. (2023) Beltz, H., Rauscher, E., Kempton, E. M.-R., Malsky, I., & Savel, A. B. 2023, AJ, 165, 257
  * Beltz et al. (2021) Beltz, H., Rauscher, E., Roman, M. T., & Guilliat, A. 2021, AJ, 163, 35
  * Betancourt (2017) Betancourt, M. 2017, The Convergence of Markov chain Monte Carlo Methods: From the Metropolis method to Hamiltonian Monte Carlo
  * Bingham et al. (2018) Bingham, E., Chen, J. P., Jankowiak, M., et al. 2018, Pyro: Deep Universal Probabilistic Programming, arXiv:1810.09538
  * Borsa et al. (2021) Borsa, F., Allart, R., Casasayas-Barris, N., et al. 2021, A&A, 645, A24
  * Borsa et al. (2019) Borsa, F., Rainer, M., Bonomo, A. S., et al. 2019, A&A, 631, A34
  * Borsato et al. (2023) Borsato, N. W., Hoeijmakers, H. J., Prinoth, B., et al. 2023, A&A, 673, A158
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs
  * Cadiou (2022) Cadiou, C. 2022, Matplotlib label lines, language: eng
  * 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., Palle, E., Yan, F., et al. 2019, A&A, 628, A9
  * Casasayas-Barris et al. (2018) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2018, A&A, 616, A151
  * Cauley et al. (2019) Cauley, P. W., Shkolnik, E. L., Ilyin, I., et al. 2019, AJ, 157, 69
  * Cauley et al. (2021) Cauley, P. W., Wang, J., Shkolnik, E. L., et al. 2021, AJ, 161, 152
  * Cegla et al. (2016) Cegla, H. M., Lovis, C., Bourrier, V., et al. 2016, A&A, 588, A127
  * Cont et al. (2021) Cont, D., Yan, F., Reiners, A., et al. 2021, arXiv:2105.10230 [astro-ph]
  * Cosentino et al. (2012) Cosentino, R., Lovis, C., Pepe, F., et al. 2012, 8446, 84461V
  * Deline et al. (2022) Deline, A., Hooton, M. J., Lendl, M., et al. 2022, A&A, 659, A74
  * Delrez et al. (2016) Delrez, L., Santerne, A., Almenara, J. M., et al. 2016, MNRAS, 458, 4025
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, Journal of Open Source Software, 1, 24
  * Fortney et al. (2021) Fortney, J. J., Dawson, R. I., & Komacek, T. D. 2021, Journal of Geophysical Research (Planets), 126
  * Froese-Fischer et al. (1997) Froese-Fischer, C., Brage, T., & Johnsson, P. 1997, Computational atomic structure: an MCHF approach (CRC press)
  * Gandhi et al. (2023) Gandhi, S., Kesseli, A., Zhang, Y., et al. 2023, AJ, 165, 242
  * Helling et al. (2021) Helling, C., Lewis, D., Samra, D., et al. 2021, A&A, 649, A44
  * Hoeijmakers et al. (2020a) Hoeijmakers, H. J., Cabot, S. H. C., Zhao, L., et al. 2020a, A&A, 641, A120
  * Hoeijmakers et al. (2018) Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018, Nature, 560, 453
  * Hoeijmakers et al. (2019) Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165
  * Hoeijmakers et al. (2022) Hoeijmakers, H. J., Kitzmann, D., Morris, B. M., et al. 2022
  * Hoeijmakers et al. (2020b) Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020b, A&A, 641, A123
  * Kausch et al. (2015) Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78
  * Kesseli & Snellen (2021) Kesseli, A. Y. & Snellen, I. 2021, ApJL, 908, L17
  * Kesseli et al. (2022) Kesseli, A. Y., Snellen, I. A. G., Casasayas-Barris, N., Mollière, P., & Sánchez-López, A. 2022, AJ, 163, 107
  * Kitzmann et al. (2018) Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183
  * Kitzmann et al. (2023) Kitzmann, D., Hoeijmakers, H. J., Grimm, S. L., et al. 2023, A&A, 669, A113
  * Komacek et al. (2022) Komacek, T. D., Tan, X., Gao, P., & Lee, E. K. H. 2022, ApJ, 934, 79
  * Kramida et al. (2022) Kramida, A., Yu. Ralchenko, Reader, J., & and NIST ASD Team. 2022, NIST Atomic Spectra Database (ver. 5.10), [Online]. Available: https://physics.nist.gov/asd
  * Kumar et al. (2019) Kumar, R., Carroll, C., Hartikainen, A., & Martin, O. 2019, Journal of Open Source Software, 4, 1143
  * Lee et al. (2022) Lee, E. K. H., Prinoth, B., Kitzmann, D., et al. 2022, MNRAS, 517, 240
  * Lendl et al. (2020) Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94
  * Lodders (2003) Lodders, K. 2003, The Astrophysical Journal, 591, 1220
  * Lothringer et al. (2018) Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27
  * Maguire et al. (2023) Maguire, C., Gibson, N. P., Nugroho, S. K., et al. 2023, MNRAS, 519, 1030
  * Mayor et al. (2003) Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20
  * McKemmish et al. (2019) McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836
  * Merritt et al. (2021) Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2021, MNRAS, 506, 3853
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, 17
  * Pelletier et al. (2023) Pelletier, S., Benneke, B., Ali-Dib, M., et al. 2023, Nature
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96
  * Phan et al. (2019) Phan, D., Pradhan, N., & Jankowiak, M. 2019
  * Pino et al. (2022) Pino, L., Brogi, M., Désert, J. M., et al. 2022, A&A, 668, A176
  * Pino et al. (2020) Pino, L., Désert, J.-M., Brogi, M., et al. 2020, ApJL, 894, L27
  * Prinoth et al. (2022) Prinoth, B., Hoeijmakers, H. J., Kitzmann, D., et al. 2022, Nat Astron, 6, 449
  * Savel et al. (2022) Savel, A. B., Kempton, E. M. R., Malik, M., et al. 2022, ApJ, 926, 85
  * Savel et al. (2023) Savel, A. B., Kempton, E. M. R., Rauscher, E., et al. 2023, The Astrophysical Journal, 944, 99
  * Sedaghati et al. (2021) Sedaghati, E., MacDonald, R. J., Casasayas-Barris, N., et al. 2021, MNRAS, 505, 435
  * Seidel et al. (2023) Seidel, J. V., Borsa, F., Pino, L., et al. 2023, A&A, 673, A125
  * Seidel et al. (2021) Seidel, J. V., Ehrenreich, D., Allart, R., et al. 2021, A&A, 653, A73
  * Seidel et al. (2019) Seidel, J. V., Ehrenreich, D., Wyttenbach, A., et al. 2019, Astronomy and Astrophysics, 623, A166
  * Seifahrt et al. (2020) Seifahrt, A., Bean, J. L., Stürmer, J., et al. 2020, 11447, 114471F
  * Seifahrt et al. (2018) Seifahrt, A., Stürmer, J., Bean, J. L., & Schwab, C. 2018, 10702
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
  * Showman & Polvani (2011) Showman, A. P. & Polvani, L. M. 2011, ApJ, 738, 71
  * 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
  * Spiegel et al. (2009) Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487
  * Stangret et al. (2022) Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2022, A&A, 662, A101
  * Stangret et al. (2020) Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2020, A&A, 638, A26
  * Stock et al. (2022) Stock, J. W., Kitzmann, D., & Patzer, A. B. C. 2022, MNRAS, 517, 4070
  * Stock et al. (2018) Stock, J. W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, MNRAS, 479, 865
  * Tabernero et al. (2021) Tabernero, H. M., Osorio, M. R. Z., Allart, R., et al. 2021, A&A, 646, A158
  * Wardenier et al. (2021) Wardenier, J. P., Parmentier, V., & Lee, E. K. H. 2021, MNRAS, 510, 620
  * West et al. (2016) West, R. G., Hellier, C., Almenara, J. M., et al. 2016, A&A, 585, A126
  * Yan et al. (2019) Yan, F., Casasayas-Barris, N., Molaverdikhani, K., et al. 2019, A&A, 632, A69
  * Yan et al. (2022) Yan, F., Pallé, E., Reiners, A., et al. 2022, A&A, 661, L6
  * Yan et al. (2020) Yan, F., Pallé, E., Reiners, A., et al. 2020, A&A, 640, L5
  * Yan et al. (2021) Yan, F., Wyttenbach, A., Casasayas-Barris, N., et al. 2021, A&A, 645, A22

##  Appendix A Interference pattern

Figure 10: Normalised order showing the interference pattern caused by
internal reflections in the Coudé train of ESPRESSO. The spectra have been
corrected for telluric contamination, outliers have been masked out and colour
correction has been applied. The order is then divided by the out-of-transit
baseline. It shows that UT-1 (out of transit) is likely less affected than
UT-2. Figure 11: Comparison between cross-correlation results for the Fe
detection in the ESPRESSO transit for the sinusoidal correction in Eq. 1.
Left: Cross-correlation result for Fe applying the sinusoidal correction.
Right: Cross-correlation result for Fe without applying the sinusoidal
correction. The peak significance of the detection including the sinusoidal
fit is marginally larger, but insignificantly so.

##  Appendix B Templates

Figure 12: Periodic table of the elements indicating the considered cross-
correlation templates in this work. Pink frames indicate the availability of
sufficient absorption lines at $3000\text{\,}\mathrm{K}$ (neutral atoms) and
$4000\text{\,}\mathrm{K}$ (ions) in the wavelength regions of HARPS & HARPS-N,
ESPRESSO, and MAROON-X. Blue frames indicate the availability of sufficient
absorption lines at $3000\text{\,}\mathrm{K}$ (neutral atoms), but not for the
ionised atom at higher temperatures. Templates for TiO, VO, and AlO are also
available at $2500\text{\,}\mathrm{K}$. The cross-correlation analysis is
performed for all of the available templates as indicated in this table and
for the oxides mentioned above (90 species in total).

##  Appendix C Trace-fit posteriors and best-fit models

Figure 13: Posterior distributions of the model parameters of Eq. 6 for H. The
top-right panel shows the two-dimensional cross-correlation function of all
eight transit times series combined binned to a common phase grid (top) and
the model with parameters equal to the median of their posterior distributions
(bottom).

Figure 14: Same as Fig. 13, but for Na. Note that for Na, the projected
orbital velocity was forced to be equal for both regions of the transit as the
posterior distribution reached the bounds of feasible priors. To avoid
correlations, we opted for fitting only one orbital velocity instead of two.

Figure 15: Same as Fig. 13, but for Mg. Note that we allow two orbital
velocities to be fit despite the weak signal in the first half of the transit.

Figure 16: Same as Fig. 13, but for Ca.

Figure 17: Same as Fig. 13, but for Ca+.

Figure 18: Same as Fig. 13, but for Ti.

Figure 19: Same as Fig. 13, but for Ti+.

Figure 20: Same as Fig. 13, but for TiO.

Figure 21: Same as Fig. 13, but for V. Note that we allow two orbital
velocities to be fit despite the weak signal in the first half of the transit.

Figure 22: Same as Fig. 13, but for Cr.

Figure 23: Same as Fig. 13, but for Mn.

Figure 24: Same as Fig. 13, but for Fe.

Figure 25: Same as Fig. 13, but for Fe+.

Figure 26: Same as Fig. 13, but for Ni.

Figure 27: Same as Fig. 13, but for Sr. Note that the orbital velocity was
forced to be the same in both halves of the transit.

Figure 28: Same as Fig. 13, but for Sr+. Note that the orbital velocity was
forced to be the same in both halves of the transit.

Figure 29: Same as Fig. 13, but for Ba+.

Figure 30: Same as Fig. 13, but for Ba+ in the atmosphere of KELT-9 b. The
signal strength at the end of the transit is weaker than at the start, which
agrees with the prediction that Ba+ should ionise to Ba^2+ on the hot dayside.
Note that we only fit for one orbital velocity.

##  Appendix D Non-detections

Figure 31: Same as Fig. 4, but for non-detections. Figure 32: Same as Fig. 4,
but for non-detections. Figure 33: Same as Fig. 4, but for non-detections.
Figure 34: Same as Fig. 4, but for non-detections.