© ESO 2020

1 Introduction

Gas giant exoplanets that orbit sufficiently closely to their host star exhibit elevated equilibrium temperatures that significantly alter their atmospheric structure when compared to cooler hot Jupiters. The atmospheres of these so-called ultra-hot Jupiters (UHJs) are subject to thermal dissociation of molecules and partial thermal of atomic species. At temperatures over 2500− 3000 K, atomic hydrogen becomes the dominant atmospheric constituent (Lothringer et al. 2018; Kitzmann et al. 2018; Parmentier et al. 2018). The dissociation of hydrogen reduces the mean molecular weight, thus increasing the pressure scale height, and electrons liberated by the thermal ionisation of metals (mostly alkalis) combine with atomic hydrogen to form H⁻, which is a strong source of continuum opacity (Arcangeli et al. 2018). Indeed, the dominance of atomic hydrogen and the resulting importance of H⁻ in shaping the radiative properties of the atmosphere may be a defining characteristic of UHJs.

Ultra-hot Jupiters are prime targets for observations that target the transmission spectrum of the atmosphere during transit events. The atmospheric scale height is large due to the high temperature and low mean molecular weight, while the short orbital period allows transits to be observed frequently. From a theoretical point of view, the atmospheres of the hottest UHJs reside in an interesting chemical regime. An absence of complex molecular chemistry precludes the formation of aerosol particles and opaque cloud decks, which have mostly defied chemical characterisation when encountered in the atmospheres of cooler hot Jupiters. Combined with the fast rates of the chemical reactions that are prevalent at high temperatures, the absence of complex molecular chemistry simplifies the theoretical interpretation of transmission and emission spectroscopy of these objects, making them especially amenable for a detailed chemical characterisation.

The UHJ KELT-9 b is currently the hottest planet in this class (Gaudi et al. 2017). With an equilibrium temperature of over 4000 K and a night-side temperature greater than 3000 K (Wong et al. 2020), its atmosphere is expected to be nearly fully dissociated and in chemical equilibrium (Lothringer et al. 2018; Kitzmann et al. 2018; Parmentier et al. 2018), thereby epitomising the defining characteristics of UHJs. The atmospheres of cooler UHJs span a transition: from being dominated by molecular hydrogen at lower temperatures to being mostly dissociated and partially ionised (Lothringer et al. 2018). Due to the large day-to-night side temperature contrast, the atmosphere of any single UHJ may span multiple regimes itself, where the atmosphere on the strongly heated day side may be thermally dissociated whereas the cooler night side is instead governed by condensation processes and molecular chemistry (Parmentier et al. 2018; Ehrenreich et al. 2020).

WASP-121 b orbits a relatively bright (V = 10.5) F6V star in a 1.27 day period orbit (Delrez et al. 2016). The equilibrium temperature was estimated at 2358 ± 52 K (Delrez et al. 2016), prompting Evans et al. (2016) to search for the presence of TiO/VO and an atmospheric inversion layer (Hubeny et al. 2003; Burrows et al. 2007), using the WFC3 instrument aboard the Hubble Space Telescope (HST). Transit observations between 1.12–1.64 μm and optical spectro-photometry indicated clear absorption by water that is typically observed in hot Jupiter transmission spectra (Sing et al. 2016), as well as additional opacity near 1.2 μm tentatively ascribed to a combination of TiO/VO and FeH (Evans et al. 2016). Secondary-eclipse observations using the G141 grism of WFC3 (1.12–1.64 μm) revealed emission by water on the day side hemisphere of the planet, indicating the presence of a temperature inversion (Evans et al. 2017). An additional source of emission near 1.2 μm was tentatively attributed to vanadium-oxide (VO).

To confirm the importance of VO, water and an inversion layer (Evans et al. 2018; Mikal-Evans et al. 2019, 2020) obtained repeated HST observations of the transmission spectrum and the secondary eclipse using the STIS and WFC 3 instruments. The optical transmission spectrum displays rich variation, with multiple features consistent with VO absorption that Evans et al. (2018) could reproduce by assuming an isothermal T-P profile at 1500 K and a metallicity equivalent to 10× to 30× solar. Absorption bands of TiO appeared to be muted in the transmission spectrum, which was explained by Evans et al. (2018) as evidence of condensation of Ti-bearing species, which commences at higher temperatures than condensation of V-bearing species, producing for example, calcium titanates (Lodders 2002) while VO remains in the gas phase. Mikal-Evans et al. (2019) observed the day side emission spectrum with the G102 grism of WFC3 (0.8–1.1 μm), augmenting their earlier observations with the G141 grism. The G102 spectrum does not show the VO bands expected to be present there, and this led Mikal-Evans et al. (2019) to question the interpretation that the 1.2 μm feature is caused by VO emission. The secondary eclipse was observed at 2 μm (Kovács & Kovács 2019) and at optical wavelengths with the TESS instrument. These were analysed together with the preceding Hubble, Spitzer, and ground-based observations to yield tighter constraints on the atmospheric structure, composition and overall system parameters (Bourrier et al. 2020a; Daylan et al. 2019). These studies found that the hottest point on the day side exceeds a temperature of 3000 K, that the atmosphere is inverted on the day side, and a metallicity that is consistent with solar (Bourrier et al. 2020a) or slightly elevated (Daylan et al. 2019). Although the chemical retrievals follow different strategies (equilibriumversus free-chemistry), both indicate that a depletion of TiO relative to VO is needed to explain the observed emission spectrum, supporting the earlier findings by Mikal-Evans et al. (2019). Recently, Mikal-Evans et al. (2020) obtained new secondary-eclipse observations using the G141 grism of WFC3. Although confirming the presence of emission by H₂O, a joint analysis with their previous WFC3 observations did not reproduce the emission feature at 1.2 μm, prompting the authors to entirely discard their previous interpretation of emission caused by VO.

At shorter wavelengths covered by the STIS instrument, Evans et al. (2018) observed a steep increase of the transit radius, which they propose to explain as being caused by the NUV absorption bands of an unknown absorber, and explore the possibility of the SH molecule. SH has been proposed to be a significant by-product of photo-dissociation of H₂ S (Zahnle et al. 2009, 2016). Despite being a highly reactive radical, its abundance could exceed 1 ppm at the millibar level where transit transmission spectroscopy is sensitive. Evans et al. (2018) are able to reproduce the NUV slope with atmospheric models spanning temperatures of 1500–2000 K and SH abundances of 20–100 ppm, but note that an unambiguous identification of the SH molecule is beyond the potency of these low-resolution spectra. The transmission spectrum was observed at UV wavelengths using Swift/UVOS between 200 and 270 nm, yielding a tentative excess in the photometrictransit depth that evidences strongly absorbing metal ions at high altitudes (Salz et al. 2019). Two transits observed using the high-resolution E230M echelle grating of HST/STIS between 228 and 307 nm revealed strong absorption lines by Fe II and Mg II (Sing et al. 2019). These absorption line depths are significantly greater than the transit radius of the planet’s Roche lobe, indicating that these heavy elements may be part of a hydrodynamic outflow.

The planet atmosphere has also been observed at high spectral resolution at optical wavelengths, with the HARPS (Bourrier et al. 2020b; Cabot et al. 2020) and UVES (Gibson et al. 2020; Merritt et al. 2020) spectrographs. These observations yielded confident detections of atomic metals, including Fe I and Na I, as well as H − α absorption – the latter of which is further evidence for the existence of an extended outflowing envelope. Given that Fe I is notably richin absorption lines at blue-optical and NUV wavelengths, these observations suggest that it could be responsible for the heating required to generate the observed temperature inversion in the upper atmosphere (Gibson et al. 2020; Pino et al. 2020), and may also explain the observed slope towards NUV wavelengths (Lothringer et al. 2020).

Merritt et al. (2020) further investigated the UVES spectra published by Gibson et al. (2020) in search for TiO and VO absorption. These authors report non-detections of both molecules, with an upper limit on the TiO abundance of [TiO ] ≲−9.3, consistent with the retrieval of Evans et al. (2018) and the interpretation that TiO is condensed out of the gas phase. Merritt et al. (2020) also establish an upper limit on the VO abundance, but note that existing line-list are likely not accurate enough for application with high-resolution spectroscopy.

High-resolution spectrographs with $R=\frac{\lambda }{\Delta \lambda }=\frac{c}{\Delta v}~{10}^5$ like HARPS and UVES are typically able to measure the centroid velocities of atmospheric signatures with sensitivities close to 1 km s⁻¹, making such observations sensitive to global atmospheric dynamics (Snellen et al. 2010; Brogi et al. 2016; Flowers et al. 2019). Bourrier et al. (2020b), Cabot et al. (2020) and Gibson et al. (2020) independently report blueshifts of − 5.2 ± 0.5, $-3_{-1}^{+3}$ and − 4.4 ± 0.6 km s⁻¹, indicating the presence of a wind that carries these atoms from the day side to the night side at high altitudes.

In this paper we present an analysis of the high-resolution transmission spectrum of WASP-121 b, using the three transit events observed with the HARPS spectrograph first described by Bourrier et al. (2020b). We apply a classical differential analysis targeting the strong sodium (Na I) doublet following the strategy of Wyttenbach et al. (2015) and Seidel et al. (2019) as well as the cross-correlation technique (Snellen et al. 2010) to search for additional atomic and ionised metals, following the analyses of Hoeijmakers et al. (2018a, 2019). Section 2 describes the observations and details the transmission spectroscopy and cross-correlation analyses. Section 3 presents the detected species and discusses their implications for the chemical composition and structure of the atmosphere of WASP-121 b. Section 4 enumerates the conclusions with reference to the most important figures and tables. Appendices B–D provide supporting figures of the cross-correlation procedure, all obtained cross-correlation functions and a detailed description of our bootstrap methods used to assess detection robustness.

Fig. 1 S/N at the centre of order 56 (near the location of the Na I doublet) in all three nights as a function of orbital phase. The vertical lines indicate the start and end of the transit. The second night, in blue, shows markedly lower S/N during transit. For this reason, it is rejected during the focused analysis of the Na doublet.

2 Observations and analysis process

2.1 HARPS observations

We obtained spectra during three transits of the hot gas giant planet WASP-121b around its host star WASP-121 (spectral type F6, V = 10.4), with the HARPS instrument at the ESO 3.6 m telescope in La Silla Observatory, Chile (Mayor et al. 2003). The observations were performed on the 31st Dec. 2017, the 9th Jan. 2018 and 14th Jan. 2018 as part of the HEARTS survey (ESO programme: 100.C-0750; PI: Ehrenreich). A log of the observations is provided in Table 1. Flat field and wavelength calibration frames are obtained during the daily afternoon calibration. All science observations are performed with fibre A on the target and fibre B on the sky.

We used the HARPS Data reduction pipeline (DRS, version 3.8), of which the primary science products are the individually extracted Echelle orders from the 2D Echellogram, as well as blaze-corrected, stitched and resampled one-dimensional spectra. The wavelength solutions provided by the DRS are in air, in the rest frame of the observatory. The signal-to-noise ratio (S/N) as recorded by the DRS at the centre of order #56 (which covers the sodium doublet) for each of the three nights are plotted in Fig. 1. Ground-based transmission spectroscopy at optical wavelengths of exoplanet atmospheres requires the correction of telluric contamination, caused by H₂O and O₂ absorption bands. Telluric contamination was corrected using the molecfit package (Smette et al. 2015; Kausch et al. 2015), following previous work by, for example, Allart et al. (2017), Hoeijmakers et al. (2019) and Seidel et al. (2019) for an application to the Na-D doublet.

We applied molecfit to the one-dimensional spectra created by the pipeline to create a model of the telluric transmission spectrum over the entire wavelength range of HARPS for each spectrum in the time-series. These telluric models were then interpolated onto the wavelength solutions of each of the respective spectral orders and divided out, which yielded individually corrected Echelle orders. This correction is valid despite the fact that the Echelle orders are not blaze-corrected, because the telluric transmission spectrum and the blaze correction are both multiplicative operations. We visually inspected the spectra and found the correction to be effective down to the noise level at most wavelengths¹.

After telluric correction, the spectra were Doppler-shifted to place the host star in a constant rest-frame. To this end we performed two velocity corrections simultaneously: the Earth’s velocity around the barycentre of the solar system, and the radial velocity of the star induced by the gravitational effect of the orbiting planet, leaving the stellar spectra at a constant velocity shift set by the systemic velocity of ~ 38 km s⁻¹. This combined velocity shift constitutes the only re-interpolation of the extracted spectra during this cross-correlation analysis. Missing edge values were masked (see Sect. 2.3.1). We proceeded to perform two independent analyses of the transmission spectrum of WASP-121: A focused analysis of the narrow waveband around the sodium D-lines (Sects. 2.2 and 3.1), and a cross-correlation analysis targeting metals and molecules with lines spread out over the waveband of HARPS (Sects. 2.3 and 3.2).

2.2 Narrow-band transmission spectroscopy

The in-transit spectra were divided by the mean of the out-of-transit spectra which constitutes a master out-of-transit spectrum. This implicitly removed the blaze, the stellar continuum flux and stellar absorption lines, and yields a time-series of transmission spectra in the rest-frame of the star. These were shifted to the planetary rest-frame using the known orbital parameters (Bourrier et al. 2020b) and averaged in time. The transmission spectra for the individual nights were then combined to create a transmission spectrum averaged over three transits. A more detailed description of this approach can be found in Wyttenbach et al. (2015). The normalisation processes, where all spectra are brought to the same flux level for a correct calculation of the in- and out-of-transit master, followed that of Seidel et al. (2019).

2.3 Cross-correlation analysis

Due to the high equilibrium temperature of this planet the transmission spectrum features absorption lines of atomic metals other than sodium, as already noted by Bourrier et al. (2020a), Sing et al. (2019), Gibson et al. (2020), and Cabot et al. (2020). We use the cross-correlation method (Snellen et al. 2010) to search for an ensemble of neutral atoms and ions, following the strategy applied in previous work (Hoeijmakers et al. 2019).

To implicitly retain information about the absolute flux recorded by the spectrograph, we have modified the analysis used by Hoeijmakers et al. (2018a, 2019) to be applied directly to the extracted spectral orders, as opposed to the blaze-corrected, stitched and resampled one-dimensional spectra that are also created by the data-reduction pipeline. Under the assumption that the noise is dominated by photon noise, retaining the absolute flux implicitly weighs spectral pixels according to their variance when performing the cross-correlation (see Sect. 2.3.3)².

First, the spectral orders (each with a width of 4096 pixels) of the time-series were sorted in time and grouped into two-dimensional matrices of 4096 × N values, where N is the number of exposures obtained per night, equal to 63, 45 and 50 respectively. Because the wavelength-solution computed by the HARPS pipeline is constant during the time-series, this re-ordering did not require resampling of the extracted spectra.

2.3.1 Masking and colour correction

We performed a two-step process to mask outlying flux values and spectral regions that are otherwise affected by artefacts. We first performed order-by-order sigma-clipping by computing a running standard deviation over a sub-region measuring 40 × N values and flagging 5σ outliers from the mean as missing data. Secondly, we visually inspected all orders and select spectral channels (i.e. 1 × N columns) where detector artefacts, regions with excessive noise (e.g. at the edges of the bluest orders), stellar residuals (e.g. in the cores of deep absorption lines) or imperfectly removed tellurics are apparent. These columns were also flagged.Both masking steps affected a total of 1.1, 2.7 and 1.6% of all the values in each of the three nights respectively.

Unless an entire spectral column is missing from the time-series, missing values would affect the evaluation of the cross-correlation function, effectively introducing a time-dependence of non-astrophysical origin³. We therefore distinguish between missing columns and isolated missing values. Missing columns were ignored, as these spectral channels do not contribute to the cross-correlation function at any time in the numerator of the cross-correlation function. Conversely, isolated missing values were interpolated over in the spectral direction. Columns of which over 20% of values are missing were ignored from the analysis entirely. After masking, the spectral orders were duplicated and a model spectrum of the planet was injected into these data to enable model comparison (see Sect. 2.3.2). After model injection, these contaminated spectra were treated in the same way as the spectra in which no model was injected, following earlier practice (Hoeijmakers et al. 2015, 2018a, 2019).

Because the broad-band continuum of the HARPS spectrograph varies slightly over the course of the observations, we performed a colour-correction of the spectral orders by normalising the mean flux of each spectral order over wavelength to the mean flux ofeach spectrum over time. This enforces fixed ratios between the mean flux levels between the orders, consistent with the combined effect of the typical response function of the spectrograph, the transmission of the Earth’s atmosphere and the intrinsic spectrum of the star. Because this normalisation step removes information about the average flux of the exposures as a function of time, we weighted the cross-correlation functions of the time series accordingly at the end of the analysis (see Sect. 2.3.3).

Fig. 2 Abundance profiles of selected species at a temperature of 2000 K and 20× solar metallicity, as computed by FastChem. Solid lines correspond to atomic species, dashed lines to molecules, and dashed-dotted lines to atomic hydrogen and H−. Transmission spectroscopy is principally sensitive to pressures below the milli-bar level (e.g. Kitzmann et al. 2018).

2.3.2 Model spectra and cross-correlation templates

To create model spectra for purpose of model injection and as cross-correlation templates, we followed the same procedure as Hoeijmakers et al. (2019): the atmosphere of the planet was assumed to be isothermal and in chemical and hydrostatic equilibrium, with elemental abundances corresponding to a fixed metallicity value. The chemical abundance profiles were computed using FastChem (Stock et al. 2018), and the radiative transfer was performed as described in Gaidos et al. (2017), using opacity functions derived from line-lists provided by VALD and Exomol for atoms and molecules respectively (Ryabchikova et al. 2015; Tennyson et al. 2016). We included opacity from all atomic neutrals and single ions with atomic numbers up to 40, due to diminishing abundance for higher atomic numbers. Although FastChem produces abundance profiles for hundreds of molecules, we only considered molecular opacity caused by TiO (McKemmish et al. 2019), VO (McKemmish et al. 2016), SH (Gorman et al. 2019), H₂O and H₂. This choice is partly motivated by diminishing abundances of larger molecules, but also by a sparsity in reliable line-list data for most molecules.

In this way, we produced a single spectrum for each absorbing species (including continuum absorption) assuming a temperature of 2000 K, to be used as cross-correlation templates. For the purpose of model comparison, we additionally produced four model spectra assuming a temperature of 1500, 2000 or 3000 K, a metallicity of 20× solar (corresponding to the median value found by Evans et al. 2018), and including the opacity of all species, but artificially disabling the effect of TiO opacity to account for the possible effects of TiO condensation, as listed in Table 2. The abundance profiles for the model at 2000 K are shown in Fig. 2, and models at 1500 K and 2000 K are shown in Fig. 3.

The four spectra that are used for model injection were convolved with a broadening kernel to simulate the effects of tidally locked rotation of the planet (v_eq = 7.0 km s⁻¹), the resolution of the spectrograph and the change of the radial velocity of the planet between the start and end of each exposure, following the approach by Brogi et al. (2016). They were shifted to the average radial velocity of the planet during the exposure (assuming v_orb = 221.1 km s⁻¹ and v_sys = 38.043 km s⁻¹), and multiplied with a normalised model of the transit light-curve⁴ before injection into each spectrum via multiplication.

The model spectra that are used as cross-correlation templates were continuum-subtracted by fitting a low-order polynomial. Any residual values that are smaller than 1 × 10⁻⁴ times the amplitude of the deepest spectral line in the waveband were clipped to zero, ensuring a constant continuum. Given that the deepest spectral lines in the optical have typical depths of ~ 1 ×10⁻³ times the flux of the star, this procedure affects only spectral lines at the ~ 1 ×10⁻⁷ level. Then,the template was convolved with a Gaussian kernel with a FWHM width of 2.72 km s⁻¹, matching the resolving power of the spectrograph. This ensured that the template is never undersampled when interpolated onto the wavelength grid of the observed spectra⁵. All templates developed for this analysis are shown in Fig. 4.

The above templates do not contain lines by ionised metals because thermal ionisation is insignificant at 2000 K. To test for the presenceof absorption lines by ions, we used templates computed for the transmission spectrum of KELT-9 b (Hoeijmakers et al. 2019)at a temperature of 4000 K and solar metallicity. These templates are publicly available via the CDS⁶.

2.3.3 Cross-correlation

The cross-correlation operation (Eq. (1)) was applied to the extracted orders without correcting the blaze function nor the stellar absorption lines. Therefore, the correlation coefficients c(v, t) constitute flux measurements as the number of photon counts registered by the detector (averaged over the time series due to the colour-correction described above). The functional form of the cross-correlation coefficient c(v, t) as an averaging operator is: $$c(v,t)={\sum }_{i=0}^{N_x}{x}_i(t)T_i(v),$$(1)

where x_i(t) are each of the spectral points in all the echelle orders of the spectrum obtained at time t, T_i (v) are the corresponding values of the template Doppler shifted to a radial velocity, v. T takes on non-zero values inside spectral lines of interest, and is normalised such that ${\sum }_{i=0}^N{T}_i(v)=1$. The sum includes the spectral points of all echelle orders, eliminating the need for explicit order-by-order weights as applied in earlier works (Hoeijmakers et al. 2018a, 2019), as well as computationally expensive filtering steps applied to the full high-resolution spectrum. The transit-depth associated with the atmosphere of planet is obtained from c(v, t) by dividing all the cross-correlation functions of the time-series by their out-of-transit average, reminiscent of studies that use the cross-correlation functions produced by the pipeline directly (Bourrier et al. 2020a; Ehrenreich et al. 2020).

Like in previous work, we subtracted an empirical model of the Doppler shadow that is caused by the obscuration of part of the rotating stellar disc during the transit of the planet. We used the double- Gaussian model introduced by Bourrier et al. (2018) for the M dwarf GJ 436, which consists of the sum of two Gaussians with a positive central core and negative side-lobes (see Fig. 5). The resulting model multiplied by a scaling factor is subtracted from the cross-correlation function obtained for each template, minimising the sum of the squared residuals while ignoring cross-correlation values at velocities that coincide with the expected radial velocity of the planet. We also applied a high-pass filter with a width of 70 km s⁻¹ that removes residual broad-band variations that were not already removed during the colour-correction stage (Hoeijmakers et al. 2018a, 2019).

Finally, the cross-correlation functions were shifted to the expected rest-frame velocity of the planet (based on its assumed orbital velocity and ephemeris), weighted according to the mean flux in their corresponding exposures (as already mentioned in Sect. 2.3.1), masked according to whether they are expected to be in-transit (1.0) or out-of-transit (0.0), and co-added. This yields a time-averaged, one-dimensional cross-correlation function in the rest-frame of the planet, for each of the three nights. These were finally averaged again without application of weights. This entire analysis is illustrated in Fig. B.1.

Fig. 3 Three of the four models of the transmission spectrum of WASP-121 b assuming chemical equilibrium, a metallicity of 20× solar, a temperature of 1500 (blue) or 2000 K, with (yellow) and without (green) TiO opacity, used for model injection and comparison. Inset panels: models with (right) and without (left) the contribution of TiO opacity near 517 and 622 nm, both atthe high native resolution of these models (light colour), as well as broadened to include the effects of the instrumental resolution and rigid body rotation assumed for the planet (dark colour). VO absorption bands are only evident when TiO opacity is removed (the bands that remain visible in the green and blue models near, for example, 550 and 575nm), which would otherwise be masked by much stronger TiO bands. The fourth model (at 3000 K) is not plotted here, but is shown separately in Fig. 12.

3 Results and discussion

3.1 Transmission spectra

Figure 6 shows the transmission spectrum of WASP-121 b in Echelle order (#56) of the sodium doublet in the planetary rest frame. The relative depths and detection levels are calculated by fitting a Gaussian to both peaks taking into account the propagated error of the individual spectral samples. The relative depths of the Na I D-lines are hence obtained directly from the differential transmission spectrum without resorting to pass-bands in the blue and red arms of the sodium doublet as was done by Wyttenbach et al. (2015, 2017), and Seidel et al. (2019). These are summarised in Table 3.

3.1.1 Assessment of data quality

Each transitin this data-set is possibly contaminated by spurious occurrences such as stellar spots, instrumental effects or adverse observational conditions, leading to false-positive detections that do not stem from the exoplanet transit itself. To rule out a false-positive detection due to these systematic errors and estimate their likelihood, we performed a bootstrap (or empirical Monte-Carlo, EMC) analysis. We follow the approach in Redfield et al. (2008), where sub-samples of the dataset are randomly selected to be fed to the analysis, creating randomised instances of the transmission spectra. Redfield et al. (2008) explores three scenarios: taking sub-samples only from the in-transit spectra (in-in), only from the out-of-transit spectra (out-out), or from both the in-transit spectra and out-of-transit spectra respectively (in-out). If the detected Na I signature is caused by absorption in the planet atmosphere, we expect to find no signal when dividing the in-in and out-out spectra with each other. More detailed examples of this technique can be found in Wyttenbach et al. (2015) and Seidel et al. (2019).

Apart from assessing the reliability of the detection, the computed distributions serve the secondary purpose of estimating an upper boundary on the false-alarm probability of the signal, accounting for varying observation conditions and instrumental effects (Redfield et al. 2008). The standard deviation of the out-out distribution, which is unaffected by the planetary atmosphere,is used as the error on the measured absorption depth in Redfield et al. (2008). We weigh the standard deviations by the square root of the ratio of out-of-transit spectra to the overall number of spectra taken during each night, thus accounting for the biased sample selection in the out-out scenario (Wyttenbach et al. 2015). This has been explored by Astudillo-Defru & Rojo (2013), who suggest the standard deviation of the in-out scenario as an alternative. To eliminate any influence from the planetary signal, we follow Wyttenbach et al. (2015) in this analysis.

The results of the bootstrap analysis for each of the three nights can be found in Fig. 7 for the spectral order of the sodium doublet as an example. We performed the bootstrap 10 000 times for each scenario in each night, often enough for the standard deviation to not change significantly when the iterations are increased. The in-in scenario is shown in red, the out-out scenario in blue and the in-out scenario in black. As predicted, in all three nights the in-in and out-out scenarios are centred at 0. For both night one and night three, the sodium detection can clearly be seen at the same detection level as in the transmission spectra, with the black histogram offset with respect to 0. However, night two shows no detection of sodium in transmission in the EMC. The S/N in night two for the order containing the sodium doublet (see Fig. 1) shows a significantly lower S/N for the in-transit spectra (S∕N ~ 20) compared to the out-of-transit spectra (S∕N ~ 45), see Table 1 and Fig. 1. This indicates that the feature in the transmission spectrum of night two is likely a false-positive due to systematic noise and not dominated by Gaussian noise. We therefore decided to exclude night two from the transmission spectrum analysis.

Additionally we found residuals of stellar activity in the obtained spectra for all three nights, despite our correction for the stellar spectral lines (see Sect. 2.2), similar to what has been observed by Bourrier et al. (2020b). To correct this contamination, we masked a window of 1.5 km s⁻¹ around the centre of the stellar sodium line in the stellar rest-frame, containing the entire stellar line core, before shifting into the planetary rest frame.

Fig. 4 Cross-correlation templates introduced in this work prior to continuum subtraction and broadening, derived from model spectra of the atmosphere of WASP-121 b. The chemical abundances of the atoms assume chemical equilibrium with elemental abundances equivalent to a metallicity of 20× solar and a temperature of 2000 K. The templates include all atoms with atomic number up to 40 for which line absorption is significant in the HARPS wavelength range (387.4–690.9 nm) under these conditions, as well as the H2O, SH, TiO and VO molecules. SH is included as a byproduct of sulphur photo-chemistry, with a disequilibrium abundance profile artifically set to match the profile presented in Evans et al. (2018). The profile is parametrised as having a VMR of 10−4 for pressureshigher than 10−5 bar, and zero at lower pressures. Note that for clarity, the scale of the y-axis of the molecular absorption spectra is a factor of 2 smaller than the atoms. Also note that the absorption lines of SH are expected to be significantly weaker than those of the atomic metals. Templates of ionic species (Ti II, Cr II and Fe II) are adopted from the analysis of Hoeijmakers et al. (2019), calculated for the atmosphere of KELT-9 b with a temperature of 4000 K and solar metallicity.

Fig. 5 Subtraction of the Doppler shadow and emergence of the signature of the exoplanet atmosphere, following the radial-velocity curve indicated by the white line. Top panel: raw cross-correlation matrix of the first night of observationswith the Fe I template. During the transit, the Doppler shadow emerges as the near-vertical structure around 40 km s−1, because the planet is on a strongly misaligned orbit (Delrez et al. 2016; Bourrier et al. 2020b). Middle panel: best-fit model of the Doppler Shadow described as a two-component Gaussian. Bottom panel: residuals after subtracting the best-fit model from the raw cross-correlation matrix, revealing the signature of the planet atmosphere. Residuals remain in the stellar line core before and after the transit event, but these do not affect the signature of the planet atmosphere, as the latteris constructed by considering in-transit spectra only. The signature of the planet atmosphere appears as the dark slantedfeature indicated by the white line. Note that although the sign of the absorption is negative in this figure, the sign flipped further in the analysis to denote absorption, notably in Figs. 8 and 9.

Fig. 6 Transmission spectrum of WASP-121 b at the location of the Na I D-doublet. The spectrum was constructed from the observations obtained during night one and night three, as night two was ruled out as a spurious signal. Top panel: transmissionspectra at the native resolution of HARPS in grey, binned by 20 for better visibility in black. The transmission spectrum corresponding to the best-fit analytical model (see Sects. 3.1.2 and Appendix A) is shown in red, along with the residual of the data after subtraction in the bottom panel. Blue, green and yellow lines depict two forward-models based on a radially outflowing wind (Seidel et al. 2020), an optically thin torus (Oza et al. 2019; Gebek & Oza 2020) and a hydrodynamically escaping envelope (Wyttenbach et al. 2020). The latter two reproduce the significant difference in depth between the D1 and D2 lines.

Fig. 7 Distributions of the bootstrap analysis for the 12 Åcentral passband. As expected the “in-in” (red) and “out-out” (blue) distributions are centred around zero (no planetary detection) and the randomised “in-out” distribution shows a detection (black) for both night one and night three. Night two shows wide variability in all three scenarios, with the “in-out” distribution also centred at 0.

3.1.2 Possible asymmetry in the sodium doublet

Both sodium lines appear to be significantly broadened compared to what is expected from hydrostatic models, with best-fit Gaussian FWHM values of 66 ± 17 and 47.9 ± 6.6 km s⁻¹. Significant broadening of either or both of the Na I lines has been observed in other hot Jupiters, with FWHM widths in excess of 20 km s⁻¹ observed forHD 189458 b, WASP-49 b, WASP-76 b and WASP-52 b (Wyttenbach et al. 2015, 2017; Seidel et al. 2019; Chen et al. 2020). Seidel et al. (2020) introduce a physical model (MERC) that aims to explain this broadening by invoking radially outflowing winds with speeds near the escape velocity of the planet, as could be expected in the presence of strong atmospheric escape. A forward-model with a vertical wind speed of 30 km s⁻¹, an isothermal, hydrostatic⁷ profile of 3000 K and a uniform Na I abundance of logX_Na = − 2.6 is superimposed on the observed transmission spectrum in Fig. 6. This model qualitatively reproduces the strength and width of the stronger D2 line, but we note a difference between the strengths of the D1 and D2 lines (2.7 × 10⁻³ and 5.6 × 10⁻³ respectively, with a ratio of 2.10 ± 0.62), the former being significantly over-predicted by the MERC model. This difference has also been observed in a recent, independent analysis by Cabot et al. (2020).

Given the limited signal-to-noise of the detected sodium doublet, and a possibility for systematic errors (one of the three nights of observation used in this analysis has been discarded), it is possible that the observed difference between the D-lines includes uncorrected systematic effects, and we propose that future observations of the sodium doublet in WASP-121 b would be valuable to confirm this difference. If the D2 line is indeed significantly deeper than the D1 line, we propose the following astrophysically motivated hypothesis.

3.1.3 An optically thin envelope

Assuming a temperature of 2000 K a mean particle weight of 2.3 u and an atmosphere that is in hydrostratic equilibrium, the depths of the D1 and D2 lines can be expressed in terms of the atmospheric scale height H, corresponding to (14.1 ± 3.8)H and (29.6 ± 3.4)H, respectively, with a difference of (15.5 ± 5.1)H. Isothermal hydrostatic theory predicts that the difference between the depths of the two sodium lines in units of H is proportional to logarithm of the oscillator strengths (Brown 2001; Lecavelier Des Etangs et al. 2008a; Benneke & Seager 2012; Heng et al. 2015; Heng & Kitzmann 2017): $$R({\lambda }_{\text{D}2})-R({\lambda }_{\text{D}1})=H\ln \frac{{\kappa }_{\text{D}2}}{{\kappa }_{\text{D}1}}=H\ln \frac{f_{\text{D2}}}{f_{\text{D1}}}=0.69H.$$(2)

Therefore, an altitude difference of 15 scale heights would not be expected to be observed, and line ratios significantly greater than unity are indeed not typically observed in other hot Jupiters with the exceptions of WASP-69 b (Casasayas-Barris et al. 2017) and possibly HD 189733 b (Wyttenbach et al. 2015).

Under the assumption of an isothermal hydrostatic structure (Brown 2001), the number density profile of Na I particles is modelled as an exponential function of altitude. Depending on the opacity κ, the optical depth crosses a τ ~ 1 surface at some altitude, which is approximately equal to the transit radius at wavelength λ (Lecavelier Des Etangs et al. 2008a; Heng & Kitzmann 2017).

Recently, Gebek & Oza (2020) showed that the Na-D line ratio in the transmission spectrum approaches a value of 2 when the absorbing Na I gas is optically thin (τ ≪ 1) across both absorption line cores: In the optically thin limit, the line depth ratio is proportional to the ratio of oscillator strengths (Draine 2011), meaning ~2 for the Na I D-doublet. In this scenario, the absorption lines can be used to measure the column density of absorbing sodium atoms (Draine 2011).

In Appendix A we present an analytical model of the transmission spectrum of a transiting, optically thin Na I cloud associated with WASP-121 b assuming Gaussian line shapes. We demonstrate the efficacy of this scenario by deriving the total mass of Na I needed to produce the observed absorption lines in the transmission spectrum of WASP-121 b as a function of the radii of the planet (R_p) and the star (R_*) and the observed line depths (W) and widths (Δv): $$M_{\text{Na}}=\pi R_{*}^2\left(1-\frac{R_p^2}{R_{*}^2}\right)\frac{m_{e}c}{\sqrt{\pi }{e}^2}\frac{\sqrt{2}\Delta v}{f{\lambda }_0}{m}_{\text{Na}}W({\lambda }_0).$$(3)

Through Eq. (3), each line provides a best-fit value for M_Na given its Gaussian fit parameters (see Table 3) which yield an error-weighted average of the total mass column of Na I of M_Na = (3.6 ± 0.6) × 10¹⁰ g and δv = 21.37 ± 2.61 (FWHM = 50.32 ± 6.15). This result is independent of the size or shape of the absorbing cloud, as long as it is smaller than the extent of the stellar disc and as long as it is indeed optically thin for a cloud of this size (a smaller cloud needs a larger optical depth to produce a given line depth and vice versa). When assuming that the cloud takes the (projected) shape of a ring around the planet with an outer radius of 2R_p and an inner radius of R_p (see Fig. A.1, A_c = 1.48 × 10²¹ cm²) and filling in Eq. (A.4), we obtain optical depths in the D-line cores of τ_D1 = (5.7 ± 1.6) × 10⁻² and τ_D2 = 0.120 ± 0.016, which means that the optically thin approximation is justified for a cloud that is contained within the Roche radius of the planet (R_L~ 2R_p, Sing et al. 2019).

The presence of Na I at high altitudes may be consistent with the observation that WASP-121 b possess an extended outflowing envelope (Cabot et al. 2020) that contains metals (Sing et al. 2019). As shown above, if the number density of Na I is low but extended enough, such an outflow may be largely optically thin throughout the D-line cores. We fitted the observed sodium lines in WASP-121b by using the PAWN model developed by Wyttenbach et al. (2020). This 1D model retrieves structural parameters of a dynamically escaping outflow from high-resolution transmission spectra. We tested an isothermal hydrodynamic atmospheric structure with a transonic Parker wind solution. The atmosphere is considered in chemical equilibrium and the number densities of the electronic energy levels of sodium are kept free to allow for departures from thermodynamic equilibrium (non-LTE). The best solution is found for a temperature of $T={3700}_{-300}^{+600}$ K, and a mass loss rate of log₁₀Ṁ = 9.0 ± 2.8 g s⁻¹, and is shown in Fig. 6. This model exhibits significant asymmetry in the line-depths, indicating that the optical depth is small at most altitudes. If the Na I doublet is indeed formed in a hydrodynamic outflow, these observations provide further evidence that atmospheric escape could play a non-negligible role in shaping the evolution of WASP-121b.

Scaling the estimated photo-ionisation rate published by Huebner & Mukherjee (2015), Oza et al. (2019) calculate the typical lifetimes of neutral sodium at low pressures in the upper atmospheres of a sample of hot Jupiters and find that these are on timescales of minutes, implying that an extended optically thin cloud of neutral sodium is not stable against photo-ionisation. Gebek & Oza (2020) present a model of Na I envelopes around hot Jupiters, based on the empirical density profile of the Na I torus that surrounds Jupiter. Sputtering and charge exchange processes feed this envelope with fast (~ 10−100) km s⁻¹ neutral sodium atoms, replenishing atoms that are lost due to photo-ionisation while simultaneously providing significant non-thermal line broadening (Wilson et al. 2002, and references therein). Figure 6 includes a forward-model of the resulting transmission spectrum assuming a total Na I mass of 7.2 × 10¹⁰ g and a mean Na I velocity of $<v>=22$ km s⁻¹. Like the optically thin toy-model, this model is also able to reproduce line depth ratios significantly greater than unity.

Regardless of the physical origin of the observed Na I D-line ratio in the present observations, verification of significant Na I line broadening and depth-ratios different from unity as observed in the transmission spectra of some hot Jupiters may provide important clues to the structure and origin of their envelopes, and open the doors to a better understanding of these environments via ab initio modelling (e.g. Steinrueck et al. 2019; Showman et al. 2019) or data driven exploration (e.g. Waldmann et al. 2015; Fisher et al. 2020; Oza et al. 2019; Seidel et al. 2020; Gebek & Oza 2020).

Fig. 8 Summary of co-added cross-correlation functions in velocity–velocity space (top panels) and in the planet rest-frame (bottom panels), for various species. Dashed lines indicate the velocities at which the signature of the planet is expected to occur, given the known systemic velocity and orbital parameters of the planet. Red and yellow lines are cross-correlation functions obtained by injecting two different models into the data prior to cross-correlation. The cross-correlation functions of all surveyed species are provided in Fig. D.1. Note that although the sign of the absorption is positive in this figure, the sign was flipped earlier in the analysis to denote the flux received from the system, notably in Figs. 5 and 6. Also note that the scaling of the vertical axis of the cross-correlation function of Mg I is multiplied by a factor of three, because the line of Mg I is significantly deeper than the other atoms.

Fig. 9 Similar to Fig. 8, but for molecular species. No molecules are detected using the cross-correlation technique. The y-axis values of the 1D cross-correlation functions of H2O and SH are increased by factors of 2 and 10 respectively, to accommodate for their larger scatter. The injected model of SH (green) includes only continuum opacity and line opacity by SH at a temperature of 2000 K.

3.2 Cross-correlation analysis

Figure 8 summarises the cross-correlation functions of the most important atomic species surveyed in this analysis, showing detections of Na I, Mg I, Ca I, V I, Cr I and Fe I, as well as a detection of Ni I and a non-detection of Ti I. Similarly, Figs. 9 and 10 show the cross-correlation functions of TiO, VO, H₂O, SH, Ti II, CrII and Fe II, none of which are detected. We fit the detected species with Gaussian profiles to determine the line positions, depths and widths. Because the cross-correlation is performed in steps of 1 km s⁻¹, and the templates were broadened to a FWHM of 2.7 km s⁻¹ (see Sect. 2.3.2) neighbouring points in the cross-correlation function are correlated with each other up to a range of ~ 3−4 steps. Whenfitting Gaussian profiles, we therefore fitted only to each fourth sample of the cross-correlation function, following the approach by Collier Cameron et al. (2010). These fitting results are shown in Table 4.

We performed two types of bootstrap analysis to confirm the strength and confidence level of each. The first bootstrap method is similar to the one applied to the lines of the sodium doublet in the transmission spectrum, comparing the distributions obtained by randomly selecting subsets of the in and out-of-transit spectra (Sect. 3.1.1). The second bootstrap method uses the wide velocity range over which the cross-correlation is performed to construct a distribution of the systematic noise at Doppler velocities away from the main planetary and stellar signals (i.e. near v_sys. This is achieved by randomly Doppler-shifting each of the cross-correlation functions of the 2D time-series (bottom panel of Fig. 5) of each species before time-averaging the series, and measuring the distribution of Gaussian fit parameters that occur at random. Both bootstrap analyses are detailed further in Appendix C.

In Sect. 2.3.2 we introduce four models that are injected into the data to enable model comparison. The resulting cross-correlation strengths associated with these is overlaid onto the data in Figs. 8–10, and the ratio of the injected cross-correlation strength with the observed strength is provided in Table 4. Given this comparison, we can draw the following conclusions:

• 1.

  As is the case in other ultra-hot Jupiters, Fe I, Mg I and Cr I are among the strongest absorbing species, with Mg I absorption (dominated by the resonant triplet near 517 nm) being optically thick up to altitudes significantly greater than the other species. The effect of line absorption by Fe I and other metals in the high-resolution transmission of WASP-121 b has already been observed by three independent studies by Gibson et al. (2020), Cabot et al. (2020) and Bourrier et al. (2020a). Our analysis confirms the presence of metal absorption claimed by these authors.

• 2.

  Ions are not detected despite the strong detections of Fe II and Mg II presented by Sing et al. (2019) in the NUV. This is because these NUV lines are resonant lines, which have much higher oscillator strengths than any lines of these species in the optical. Moreover, due to their smaller oscillator strengths, the absorption lines probed in the optical become optically thick at higher pressures (i.e. lower altitudes) than those in the NUV, where the degree of thermal ionisation is higher (see Fig. 2). Interestingly, Fe II has recently been detected in the high-resolution transmission spectrum of MASCARA-2 b (Stangret et al. 2020; Hoeijmakers et al. 2020), which has an equilibrium temperature that is approximately 100 K cooler than WASP-121 b. A detection of Fe II was recently reported in Ben-Yami et al. (2020), who use the same data as the present study. This Fe II feature appears to extend over a wide range of values for the assumed orbital velocity at a relatively low significance, and has an average line-depth of 0.23%, significantly greater than the noise level apparent in our cross-correlation function (see Fig. 10). We therefore conclude that this is likely a false-positive detection. Although discrepancies in the results of cross-correlation based studies are sensitive to the choice of template, we note that the Fe II template used in this study has been applied successfully to detect the Fe II feature in MASCARA-2 b, which has a temperature that is similar to WASP-121 b (Hoeijmakers et al. 2020).

• 3.

  The absorption strengths of all atoms are under-predicted by factors between 1.5 and 8. We hypothesise that this is due to the atmosphere being extended beyond what is expected assuming hydrostatic equilibrium. The presence of an extended atmosphere has already been evidenced by observations of ionised metals in the near UV (Sing et al. 2019) and strong H-α absorption in the optical (Cabot et al. 2020). Significant underprediction by hydrostatic models is common in earlier analyses of other ultra-hot Jupiters (e.g. Hoeijmakers et al. 2018a, 2019). We note that most of the detections presented inthis work would not have been possible without the presence of an extended atmosphere, and that the planning of spectroscopic observations of WASP-121 b or similar planets should take this enhancement of atomic spectral lines into account. We therefore advise against the use of model-dependent metrics to predict which atomic species might be detectable using cross-correlation spectroscopy. For example, the metric introduced by Ben-Yami et al. (2020) predicts that Ca I and Mg I should be difficult to observe and these species were consequently excluded from their analysis, even though these are confidently detected in the present study using the same data.

• 4.

  Ti I is not detected, and indeed is not expected to be a stronger absorber than V I according to our model. This is discussed further in the next section.

• 5.

  This analysis is not expected to be sensitive to VO, H₂O nor SH, even though these have significant optical absorption bands. This is primarily due to masking by other lines, which is notably relevant for SH at blue optical wavelengths and TiO/VO, that compete with each other across most of the optical. H₂O is expected to be masked by TiO/VO at all wavelengths. Secondly, lines within molecular bands are closely separated, and form blends when taking into account broadening caused by the finite spectral resolution and the rotation of the planet, which dominates at v_eq = 7.0 km s⁻¹. Finally, although the four molecules in question have rich absorption bands with a relatively large integrated opacity, the opacity of individual lines is small compared to the opacity of single atomic lines. This makes the depths of molecular bands weaker than the depths of atomic lines. All three of these effects are visually evident from Figs. 3 and 4. The consequences for the SH molecule are discussed further in Sect. 3.4.

• 6.

  TiO could have been detected at an expected abundance of ~ 2 × 10⁻⁶ (see Fig. 9), assuming a reliable line list is used. Merritt et al. (2020) presents a deep search for TiO and VO in similar transit observations of WASP-121 b using the UVES spectrograph, and find an upper limit on the TiO abundance of 10^−9.26, far below the TiO abundance expected from our model. This limit is significantly deeper than the sensitivity of these data, primarily due to the fact that UVES has additional coverage of red wavelengths where TiO has strong absorption bands. Merritt et al. (2020) also report an upper limit of 10^−7.88 on the abundance of VO, but note that the most accurate line list presently available for VO (by the Exomol group Tennyson & Yurchenko 2012) is not accurate enough for application of the cross-correlation technique.

Fig. 10 Similar to Fig. 8, but for ions. No ions are detected using the cross-correlation technique.

3.3 Ti-V chemistry and condensation

This analysis provides direct evidence for the existence of neutral, atomic vanadium in the atmosphere of WASP-121 b. A detection of V I contrasts with a non-detection of Ti I, although Ti I is approximately 10× more abundant than V I in the Sun and hence could be expected to be a stronger absorber, as is the case in stars (Asplund et al. 2009). This section aims to elucidate this dichotomy.

Figure 2 shows the abundance profiles resulting from our equilibrium chemical model at 2000 K, 20× solar metallicity. Under these conditions, Ti and V are mostly locked up in their respective oxides, with TiO being an order of magnitude more abundant than VO. However, the transmission spectrum observed by Evans et al. (2018) shows the presence of VO bands, while TiO bands appear to be absent. Evans et al. (2018) suggest that this indicates that TiO is condensed out of the atmosphere, due to it having a higher condensation temperature than VO (Lodders 2002). In this scenario, Ti I would be depleted along with TiO, as chemical equilibrium drives the remaining atomic Ti into its oxide form in response to the decreasing TiO abundance until all Ti bearing species are condensed out of the gas phase.

Conversely, a detection of V I implies that vanadium is not depleted, and therefore a significant amount of vanadium should exist in the form of gaseous VO. Indeed, our equilibrium-chemistry calculation indicates that between the mbar and nano-bar levels where transmission spectroscopy is sensitive, the VMR of VO is expected to be greater than ~ 5 × 10⁻⁷ (see Fig. 2). The current findings therefore support the interpretation by Evans et al. (2018) that the banded structure in the optical transmission spectrum is caused by VO absorption. This conclusion contrasts with the interpretation in Ben-Yami et al. (2020). These authors proposed that a detection of V I is evidence of dissociation of VO that can explain the non-detection of VO by Merritt et al. (2020), who infer an upper limit on the VO VMR of 10^−7.9, which would amount to a depletion compared to a 20× solar metallicity gas at 2000 K in chemical equilibrium. Instead, we argue that if the atmosphere is indeed close to chemical equilibrium, this detection of V I is evidence for the presence of VO, not its depletion, and we agree with Merritt et al. (2020) that a non-detection of VO using high resolution optical spectroscopy cannot be confidently interpreted while line-list inaccuracies persist.

At face value, our non-detection of Ti I is consistent with the interpretation of TiO depletion put forth by Evans et al. (2018), and this position is adopted by Ben-Yami et al. (2020) as well. However, we observe that the injected model spectra (see Fig. 3) predict that Ti I is a weaker absorber than V I, even in chemical equilibrium, and is hence not expected to be detectable despite the fact that Ti is approximately 10× more abundant than V in a gas with solar elemental abundance ratios (Asplund et al. 2009).

To explain this, we investigate the abundances of all Ti and V-bearing species presently included in FastChem, which are shown in Fig. 11. Our model indicates that for a 20× solar metallicity gas at 2000 K in chemical equilibrium, Ti I is more efficiently converted into molecules than V I, to the point that the abundances of Ti I and V I are nearly equal (around 5 × 10⁻⁹, see the red curves in Fig. 11). Although FastChem includes more Ti-bearing species than V-bearing species, over most of the atmosphere these additional Ti-bearing molecules do not represent a significant fraction of the Ti I reservoir, and can therefore not explain the observed disappearance of the abundance difference between Ti I and V I. The exception is TiS, which is more abundant than Ti I below altitudes at a pressure of ~ 1 mbar. Presumably, the VS molecule could have a similarly significant abundance, but it is not included in FastChem, which could lead FastChem to overestimate the V I abundance. However, we note that the abundance profile of Ti I does not seem to react to the steeply decreasing abundance of TiS with altitude. This leads us to conclude that the claim that Ti I and V I have equal abundance is robust against the fact that VS is missing from FastChem.

It is apparent from Fig. 11 that the dioxide molecules TiO₂ and VO₂ represent important fractions of the total Ti and V reservoirs. Especially VO₂ is expected to be approximately twice as abundant as VO. As noted previously, we observe that TiS is an important constituent as well, but that TiH is not, with a diminishing abundance that does not exceed 10⁻¹⁰. Although we presently have no knowledge of the absorption cross-sections of gaseous TiO₂ and VO₂ at high temperatures, our model indicates that these molecules can be important sources of opacity if they are optically active to a similar degree as their mono-oxide counterparts. TiO₂ is known to exhibit electronic transitions at optical wavelengths at cryogenic temperatures (Garkusha et al. 2008). TiS bands have been observed in the spectra of Mira-type AGB stars between 800 nm and 1.35 μm (Jonsson et al. 1992), and VS is measured to be active at near-infrared wavelengths (Bauschlicher & Langhoff 1986). We note that a significant body of experimental observations of the TiS molecules exists in the literature, (e.g. Jonsson & Launila 1993; Ran et al. 1999; Cheung et al. 2000; Pulliam et al. 2010), which may make TiS amenable for detection in future cross-correlation studies.

We expect that spectral retrievals performed on HST/STIS and WFC3 observations of transmission and emission spectroscopy of hot Jupiters (e.g. Evans et al. 2018, a.o.) may yield biased measurements of the chemical composition if significant opacity sources are not taken into account. We therefore emphasise the importance of the investigation of the absorption cross-sections of more complex metal-bearing molecules.

Finally, we note that although chromium and nickel condensation chemistry has not yet been extensively explored in the exoplanet literature, the present observations suggest condensation of these species commences at lower temperatures than titanium. This supports the notion that Ti is one of the first species to condense out of the atmosphere, with a condensation temperatureroughly 200 K higher than Cr and Ni (Lodders 2002, 2003).

Fig. 11 Abundance profiles of all Ti and V-bearing species included in FastChem, at a temperature of 2000 K and a metallicy of 20× solar. All solid lines correspond to Ti-bearing species, while dashed lines correspond to V-bearing species. This figure supports theconclusions that Ti and V are expected to have equal abundance; that both mono-oxides as well as dioxides are important constituents of the atmosphere, as well as metal sulphites; and that Ti II and V II are not important in chemical equilibrium.

3.4 The NUV slope and absorption by SH

Evans et al. (2018) reported the presence of a strong increase of the transit radius towards bluer wavelengths. Such slopes are commonly observed in the transmission spectra of hot Jupiters (e.g. Sing et al. 2016), and are primarily associated with Rayleigh scattering of H₂Or aerosols (Lecavelier Des Etangs et al. 2008a,b; Pont et al. 2008). When caused by H₂ scattering, the slope can be used to measure the atmospheric temperature because the underlying opacity function is known. This scenario is ruled out by Evans et al. (2018), because the strength of the slope in the transmission spectrum of WASP-121 b would require a temperature of 6980 ± 3660 if caused by pure H₂ scattering. This temperature is much higher than the measured day side temperature (Bourrier et al. 2020a), and thermal dissociation would preclude the existence of significant amounts of H₂ (Kitzmann et al. 2018). Evans et al. (2018) therefore hypothesise that this NUV absorption is caused by additional molecular absorption bands, and propose the SH molecule as an explanation because it has strong electronic absorption bands shortwards of ~450 nm and can exist at significant abundances (Zahnle et al. 2009).

We included SH opacity in the model spectra used in the this work, assuming an ad hoc (non-equilibrium) abundance of 10⁻⁴ for pressures higher than 10⁻⁵ bar, approximately matching the profile assumed by Evans et al. (2018). The SH line-list recently published by the exomol group (Gorman et al. 2019) is expected to be calculated at sufficient accuracy for the application of the cross-correlation technique (Yurchenko & Tennyson, priv. comm.). However regardless of line list accuracy, our models indicate that even at an abundance of 10⁻⁴, the absorption bands by SH are significantly weaker than typical lines by atomic metals, causing the SH lines to be masked as already noted by Lothringer et al. (2020). The effect of masking is evident from both Figs. 4 and 3, which show that the increase in transit radius in the SH bands is small compared to those of atomic metal lines, which increase in strength and number towards shorter wavelengths. The result of this masking is also evident from the model comparison after cross-correlation, shown in Fig. 9. Despite the presence of a plenitude of absorption lines in the HARPS bandpass between 373 and ~ 450 nm, the application of the cross-correlation analysis to the model-injected spectra does not result in an expected detection unless the model is free from absorbers other than SH.

Strong absorption by metals with large NUV cross-sections has already been observed as part of previous studies (Sing et al. 2019; Gibson et al. 2020; Cabot et al. 2020), as well as the present work. These detections increase likelihood that the observed NUV opacity can be explained by atomic metal lines alone. Figure 12 shows the transmission spectrum of WASP-121 b as observed by Evans et al. (2018) using the STIS instrument, together with a model that assumes equilibrium chemistry, 20× solar metallicity and an isothermal profile of 3000 K, binned to the resolution of the STIS spectrum. This model shows that if the temperature of the atmosphere of WASP-121 b is elevated (e.g. in the case of a significant thermal inversion), atomic metals alone can account for the observed NUV slope without breaching chemical or hydrostatic equilibrium (Lothringer et al. 2020). Moreover, we note that the anomalously strong metal lines detected by Sing et al. (2019) and Gibson et al. (2020) as well as this study, indicate that the atmosphere is significantly inflated, further strengthening the above hypothesis.

Fig. 12 STIS/G430L transmission spectrum of WASP-121 b as published by Evans et al. (2018) indicated with open circles. A model spectrum assuming chemical equilibrium, a temperature of 3000 K, 20× solar metallicity and SH computed as described in Sect. 2.3.2 and Table 2 is overplotted, both at the native resolution (light-blue) and binned to the resolution of STIS (dark-blue line and circles). The wavelength bin at 395 nm coincides with the Al I doublet, at which location the transit radius is strongly over-predicted compared to the STIS observations. This models predicts that the cores of the Al I doublet become optically thick at transit radii of ~ 5 × 10−3 above the local continuum, which corresponds to an altitude of ~0.3Rp above the planet. At such low pressures, photo-ionisation is expected to be important (e.g. Moses et al. 2011; Kitzmann et al. 2018). At the same time, Al I has one of the lowest ionisation thresholds of the elements, so we predict that equilibrium models may strongly over-predict the strength of the Al I doublet at this wavelength.

3.5 Radial velocity offset and winds

The cross-correlation with the Fe I template can be used as a probe to measure the radial velocity of the star. We select the out-of-transit spectra taken in each of the three nights, and fit the centroid position of the average stellar Fe I line within 15 km s⁻¹ of the core – approximately equal to the projected rotational velocity. In this way, we obtain measurements of the systemic radial velocity of 38.041 ± 0.0044, 38.0393 ± 0.0028 and 38.0487 ± 0.0034 km s⁻¹ for each of the three nights respectively, or 38.0430 ± 0.0021 km s⁻¹ on average.

The absorption lines of the planet appear at systemic velocities of 34 to 37 km s⁻¹, indicative of significant blue shifts of the absorbing species in the planet atmosphere. The signature of Fe I appears most shifted, with an effective shift of 38.0430 ± 0.0021 − 34.37 ± 0.48 = 3.67 ± 0.48 km s⁻¹, consistent with values found in the recent literature, of 4.3 ± 0.56 km s⁻¹ and 5.2 ± 0.5 km s⁻¹ by Gibson et al. (2020) and Bourrier et al. (2020b) respectively. The observation that different atoms exhibit significantly different blue-shifts may be evidence that the atmosphere is dynamically heterogeneous, or that species are distributed heterogeneously around the limb in a fashion similar to what has been observed on WASP-76 b (Ehrenreich et al. 2020).

4 Conclusion

This paper presents an analysis of the optical transmission spectrum of WASP-121 b using spectroscopic observations obtained with the HARPS spectrograph during three transit events. WASP-121 b is an ultra-hot Jupiter in the orbit of a bright (V = 10.4) F6V star with a period of 1.27 days. The planet has been targeted by a series of observational campaigns using the WFC3 and STIS spectrographs aboard the Hubble Space Telescope (Evans et al. 2016, 2017, 2018; Mikal-Evans et al. 2019; Sing et al. 2019), which revealed a rich transmission spectrum with signatures of absorption bands of VO, while TiO appears depleted because of condensation. The NUV spectrum shows a strong slope, which was hypothesised to be caused by the SH radical, which may exist as a photo-chemical by-product. In addition, deep ionic UV lines of Fe II and Mg II as well as H-α are evidence the existence of an extended, possibly escaping envelope (Sing et al. 2019; Cabot et al. 2020).

In this work we search the transmission spectrum for line absorption of the Na I doublet and perform a survey of additional atomic species with rich absorption spectra, as well as TiO, VO, H₂O and SH. We compare the resulting cross-correlation functions with forward models of the transmission spectrum of WASP-121 b assuming chemical equilibrium (computed with FastChem, Stock et al. 2018) at temperatures of 1500, 2000 and 3000 K, and elemental abundance ratios set to 20× the solar metallicity; apart from SH which is a non-equilibrium species. The findings of this work are summarised as follows:

• The Na I doublet is detected in the transmission spectrum at a confidence of 8.6σ, with an average absorption depth of (5.59 ± 0.65) × 10⁻³ times the flux of the host star (see Sect. 3.1.2, and Fig. 6).

• The Na I lines are significantly broadened compared to what is expected from hydrostatic models and rigid-body rotation of the planet when assuming tidal locking, and the line depths show a difference corresponding to 15.5 ± 5.1 times the atmospheric scale height (assuming a temperature of 2000 K and a mean particle weight of 2.3 u). This could indicate that Na I forms an optically thin envelope around the planet (Oza et al. 2019). An analytical description of a homogeneous optically thin cloud with an integrated sodium mass of (3.6 ± 0.6) × 10¹⁰ g reproduces the transmission spectrum well, as does a model by Wyttenbach et al. (2020) of a hydrodynamic outflow, and a model by Gebek & Oza (2020) that assumes an optically thin toroidal Na I envelope (see Sect. 3.1.2, Fig. 6 and Appendix A).

• Via application of the cross-correlation technique we detect neutral magnesium, calcium, vanadium, chromium, iron and nickel (see Sect. 3.2). The line-strengths of all detected species are under-predicted by model spectra that assume chemical equilibrium, which strengthens existing evidence that the atmosphere is not in hydrostatic equilibrium. In addition, the lines of all detected species are blue shifted by 1 to 4 km s⁻¹ compared to the planet rest-frame, which is evidence of atmospheric dynamics (see Table 4), and consistent with previous findings (Bourrier et al. 2020b; Cabot et al. 2020; Gibson et al. 2020).

• The presence of V I implies the presence of VO, which is expected to exist in chemical equilibrium at much higher abundance than V I. These observations therefore support the interpretation by Evans et al. (2018) that molecular bands in the transmission spectrum are caused by VO (see Sect. 3.3 and Fig. 11), instead of being indicative of VO depletion (Ben-Yami et al. 2020). We emphasise that careful model comparison is crucial to interpret the various detections and non-detections of atomic species and their chemical interactions with molecules (see below).

• Ions, including Fe II are searched for using hot (4000 K) templates earlier applied to observations of KELT-9 b (Hoeijmakers et al. 2018a, 2019), but are not detected. Ions are not expected to be important atmospheric constituents in chemical equilibrium (see Figs. 10 and 11). Although strong Fe II lines have recently been observed at high altitudes ~2R_p using HST/STIS in the UV by Sing et al. (2019), these lines are strong resonant lines with higher intrinsic opacities than lines in the optical. Having weaker oscillator strengths, these lines are formed at higher pressures than the resonantlines in the UV, where the temperature and degree of ionisation may be significantly lower. We interpret a recentdetection of Fe II in the optical transmission spectrum of this planet (Ben-Yami et al. 2020) as a likely false-positive, though further observations and comparative analyses may be needed to confirm this.

• Our models indicate that chemical equilibrium drives the abundance of Ti I and V I to near-equal values despite the fact that Ti I is ~ 10× more abundant in the Sun. This makes V I a stronger absorber than Ti I over this wavelength range (see Sect. 3.3 and Fig. 11) and implies that a non-detection of Ti I alone in this dataset is not sufficient to conclude that Ti is depleted from the atmosphere as suggested in Ben-Yami et al. (2020).

• Generally, molecules are difficult to detect due to masking, line broadening and the fact that bands intrinsically have shallower depths (i.e. effective transit radii) than individual atomic lines. No signatures of TiO, VO, H₂O or SH are detected (see Sect. 9).

• Moreover, atoms would not or barely be observable if the atmosphere at the terminator was in hydrostatic equilibrium at 2000 K (see Fig. 8), unless photo-dissociation and vertical mixing are capable of significantly increasing the abundances atomic species at higher altitudes (see Fig. 8). If the temperature in the upper atmosphere is significantly higher than 2000 K, this could act to increase the scale height in hydrostatic equilibrium (e.g. Wyttenbach et al. 2017). However, this would increase the ionisation fraction, in turn reducing opacity by neutral species.

• Without condensation, TiO is expected to reach an abundance of 2 × 10⁻⁶ in chemical equilibrium, which we expect to be barely detectable with the current data (see Fig. 9). Our non-detection of TiO is consistent with a deeper upper limit recently established by Merritt et al. (2020), and both observations are consistent with the hypothesis by Evans et al. (2018) that TiO is subject to condensation and rain-out.

• Titanium and vanadium dioxides and sulphides are expected to be significant atmospheric constituents at 2000 K (see Sect. 3.3 and Fig. 11) and may cause significant absorption at optical or near-infra-red wavelengths (Bauschlicher & Langhoff 1986; Jonsson et al. 1992; Garkusha et al. 2008). Currently, accurate line-lists do not exist for these molecules, impeding the usage of the cross-correlation technique to detect them in future studies.

• Given that atomic metal lines are inflated beyond what is predicted by models in hydrostatic equilibrium at 2000 K, it is plausible that absorption by atomic metals alone is sufficient to explain the strong NUV slope detected with by Evans et al. (2018) with HST/STIS (see Sect. 3.4 and Fig. 12).

Acknowledgements

This work has been carried out in the framework of the PlanetS National Centre of Competence in Research (NCCR) supported by the Swiss National Science Foundation (SNSF). It has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (projects Four Aces, EXOKLEIN and Exo-Atmos with grant agreement numbers 724427, 771620 and 679633, respectively), and STFC Project No. ST/R000476/1. A.W. acknowledges the financial support of the SNSF by grant number P400P2_186765. N.A.-D. acknowledges the support of FONDECYT project 3180063. The analysis presented in this work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna; Ian Crossfields’ Astro-Python Code library and Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). We thank our anonymous referee for his/her thorough reading and constructivecommentary, which helped us to significantly improve this manuscript and refine our understanding.

Appendix A Sodium lines in the optically thin limit

Following the argument expanded by Oza et al. (2019) and Gebek & Oza (2020), we start with the flux emanating from a source with spectrum F₀ (λ) that is blocked by a cloud of with optical depth τ(λ) located in the line of sight to the source (Draine 2011): $$F(\lambda )=F_0(\lambda ){\text{e}}^{-\tau (\lambda )}.$$(A.1)

This is valid if the source is totally covered by the foreground cloud. When considering a cloud with projected area A_c and optical depth τ in front of a larger star with projected area A_* (see Fig. A.1), this equation modifies to: $$F(\lambda )=F_0(\lambda )\frac{A_c}{A_{*}}{\text{e}}^{-\tau (\lambda )}+F_0(\lambda )\left(1-\frac{A_c}{A_{*}}\right).$$(A.2)

Rewriting, the line-depth with respect to the continuum is: $$F_0-F=F_0\frac{A_c}{A_{*}}\left(1-{\text{e}}^{-\tau }\right),$$(A.3)

or in units of relative transit depth, which corresponds to the equivalent width as defined by (Heng et al. 2015): $$\frac{F_0-F}{F_0}\equiv W=\frac{A_c}{A_{*}}\left(1-{\text{e}}^{-\tau }\right)=\frac{A_c}{A_{*}}\tau ,$$(A.4)

in the optically thin limit, meaning τ ≪ 1. The optical depth is either expressed in terms of a cross section multiplied by a column density (with units of cm⁻²), $$\tau =\sigma N,$$(A.5)

or an opacity multiplied by a mass column density (with units of g cm⁻²), $$\tau =\kappa \tilde{M}.$$(A.6)

For a line transition with oscillator strength f, we have: $$\sigma =\frac{\pi e^{2}f}{m_{e}c}\Phi , \text{and} \kappa =\frac{\sigma }{m_{\text{Na}}}.$$(A.7)

Consider a Gaussian line profile Φ, which is broadened due to some non-thermal velocity Δv. We assume that Δv is a constant fitting parameter: $$\Phi =\frac{c}{{\nu }_0}\frac{1}{\Delta v\sqrt{2\pi }} {\text{e}}^{-{\left(\frac{c({\nu }_0-\nu )}{\sqrt{2}\Delta v{\nu }_0}\right)}^2}.$$(A.8)

The coefficient of Φ is set by the requirement that ∫ Φ dν = 1. Rewriting, we obtain the optical depth at the line-centre (Draine 2011): $${\tau }_0=\frac{\sqrt{\pi }{e}^2}{m_{e}c}\frac{Nf{\lambda }_0}{\sqrt{2}\Delta v}.$$(A.9)

Or, in terms of mass column density and the transit depth: $$W({\lambda }_0)=\frac{A_c}{A_{*}}{\tau }_0=\frac{A_c}{A_{*}}\frac{\sqrt{\pi }{e}^2}{m_{e}c}\frac{ \tilde{M}}{m_{\text{Na}}}\frac{f{\lambda }_0}{\sqrt{2}\Delta v}.$$(A.10)

This equation contains four unknowns: A_c, W, Δv and $\tilde{M}$. Already having fit the line-depth and width of the Gaussian line-profile from the observed transmission spectrum (Sect. 3.1.2), we solve for $\tilde{M}$: $$\tilde{M}=\frac{A_{*}}{A_c}\frac{m_{e}c}{\sqrt{\pi }{e}^2}\frac{\sqrt{2}\Delta v}{f{\lambda }_0}{m}_{\text{Na}}W({\lambda }_0).$$(A.11)

Finally, in order to obtain the total mass of sodium in the cloud, we assume that all sodium is in the ground state (essentially ignoring the effect of a non-LTE level population, see Fisher et al. 2020) and multiply the mass column density with the projected area A_c of the cloud, which subsequently cancels: $$M_{\text{Na}}=A_c \tilde{M}=A_{*}\frac{m_{e}c}{\sqrt{\pi }{e}^2}\frac{\sqrt{2}\Delta v}{f{\lambda }_0}{m}_{\text{Na}}W({\lambda }_0).$$(A.12)

Finally, filling in the projected area of the star taking into account the presence of the optically thick planet with radius $R_p^2$ (see Fig. A.1): $$A_{*}=\pi R_{*}^2\left(1-\frac{R_p^2}{R_{*}^2}\right),$$(A.13)

we obtain $$M_{\text{Na}}=\pi R_{*}^2\left(1-\frac{R_p^2}{R_{*}^2}\right)\frac{m_{e}c}{\sqrt{\pi }{e}^2}\frac{\sqrt{2}\Delta v}{f{\lambda }_0}{m}_{\text{Na}}W({\lambda }_0).$$(A.14)

In the optically thin limit, the mass of absorbing sodium atoms is independent of the projected area of the cloud as long as the cloud is homogeneous and optically thin. A dependence on $\pi R_{*}^2$ remains because the same line-depth measured at a larger star requires the presence of more absorbing atoms, and vice versa.

Fig. A.1 Schematic of the geometry of an optically thin cloud with projected area Ac surrounding an optically thick planet with radius Rp, both transiting a star with radius R*, with an effective projected area A*. Equation (A.2) consists of the sum of the covered and uncovered components, Fcovered and Funcovered respectively.

Appendix B Cross-correlation analysis procedure

Fig. B.1 Pre-processing steps applied to obtain cross-correlation functions in the rest-frame of the planet, as described in Sect. 2.3.3. Panel A: an example spectral order of night one as extracted by the HARPS DRS, plotted as a time-series. Panel B: after normalisation by the time-average flux and colour correction. Panel C: cross-correlation template of Fe I after continuum-subtraction, with which the spectra are cross-correlated. Panel D: average flux in the spectral order, that is divided out of the spectra in panel A to obtain the spectra in panel B. Panel E: cross-correlation function of Fe I of the entire time series. Panel F: after removal of the time-average of the out-of-transit cross-correlation functions. Panel G: after removal of the Doppler shadow and application of the high-pass filter. Panel H: after shifting the cross-correlations to the rest-frame of the planet. Panel I: after weighing the cross-correlation functions by the average flux of the time-series (panel D) and selecting only in-transit observations. Panel J: one dimensional cross-correlation function in the rest-frame of the planet, after averaging the cross-correlation functions (panel I) in time.

Appendix C Cross-correlation bootstrap analyses

The purpose of the cross-correlation technique is to average hundreds or thousands of absorption lines to reduce the photon-noise and increase detection sensitivity (Snellen et al. 2010). Provided that sufficient spectra and spectral lines available, this technique is known to be capable to reach 1σ noise levels of ~10⁻⁶ (Hoeijmakers et al. 2018b). There are various effects that may introduce systematic noise structures at this level. Some of these do not correlate with the cross-correlation template (e.g. the presence of uncorrected cosmic-ray hits, uncorrected telluric lines when not correlating with H₂O or O₂ templates), resulting in pseudo-stochastic variations in the cross-correlation function that can occur at random times or at random radial velocities. Other effects may fully or partially correlate with the cross-correlation template, for example, residuals from stellar absorption lines (when correlating with templates of species that are present in both the planetary and the stellar spectrum) or aliasing of the cross-correlation template (i.e. when the template has few, or regularly-spaced lines).

These effects cause structures in the cross-correlation function that are not normally distributed, invalidating the usage of summary statistics such as Gaussian confidence intervals to a certain degree. The presence of such noise sources may also cause signatures that mimic planetary absorption lines, leading to the possibility of false-positive detections. We choose to employ two bootstrap methods to investigate the extent at which systematic noise in the cross-correlation function is capable of creating spurious variations that could be mistaken for planetary absorption lines. Both methods are applied to the cross-correlation functions of each species independently, and both methods assume that spurious signals will not be mistaken for planetary absorption lines if they do not resemble Gaussian line profiles with a certain width.

C.1 Bootstrap method A

The first method is similar to the bootstrap (EMC) method described in Sect. 3.1.1 and Redfield et al. (2008), in that it constructs distributions using random permutations of in and out-of-transit cross-correlation functions. It aims to assess whether the detected absorption signal is uniquely and evenly present in the in-transit CCFs only, so to reject cases where the observed signal exists in (some of) the out-of-transit CCFs, or is only present in a small fraction of the in-transit CCFs.

We start with the CCFs in which the time-averaged out-of-transit CCF and the best-fit model for the Doppler-shadow are removed. These CCFs of all three nights are split into two groups of in-transit and out-of-transit CCFs. We then obtain “master” in and out-of-transit CCFs by taking the time-average (mean CCF) of both groups. The latter of these should be close to zero, because the time-average out-of-transit CCF had already been removed for each night in the rest-frame of the star. This time-average was a constant, so it cannot have introduced or removed signals in/from only some of the CCFs, as is being tested in this analysis.

Finally, we randomly select and average 50% of the in and out-of-transit CCFs, divide these by the master in and out-of-transit CCFs to obtain one-dimensional realisations of in-in, in-out and out-out residuals. To each of these, we fitted a Gaussian at the expected location of the planet absorption line within a window of 30–44 km s⁻¹, corresponding to the observed range of velocities for our detected signals ± 4 km s⁻¹, and a width fixed to the measured value (Table 4). This procedure is repeated 20 000 times and the resulting distributions for each of the detected species are shown in Fig. C.1. These figures confirm that the detected absorption lines indeed uniquely result from the in-transit spectra, and are not likely to be reproducible via the false-positive scenarios we set up.

C.2 Bootstrap method B

The second method aims to assess the distribution of stochastically emerging signals due to variations in the CCF caused by systematic noise that is not correlated with the cross-correlation template. These variations produce structures in the CCF that can occur at any time or radial velocity, independent of the rest-frame velocities of the planet and the star.

Like for method A, we start with the CCFs in which the time-averaged out-of-transit CCF and the best-fit model for the Doppler-shadow are removed, and selecting all CCF values at radial velocities away from the planetary and stellar rest-frames, at times during the planet transit. We select only CCFs obtained during the transit because the statistics of the in-and-out of transit spectra are different. There are more out-of-transit spectra than in-transit spectra, and these tend to be observed at higher airmass. CCF values at velocities occupied by the planet at any instance during the observations (i.e. ~ ± 150 km s⁻¹) are ignored, leaving only CCF values that are not directly related to the presence of the signatures of the planet or the star.

Each CCF in the time-series is then shifted to a random radial velocity drawn from a uniform distribution, after which all CCFs are averaged. This process mimicks the co-addition of in-transit CCFs in the rest-frame of the planet, but only evaluated at velocities away from the planet rest-frame and applying an incoherent set of velocity shifts. This method ensures that any systematic correlations that exist purely in the spectral direction of the various spectra of the time-series are preserved.

The time-averaged, one-dimensional CCF is then cropped into 15 smaller bins, within each of which a Gaussian profile is fit at a random location and with a fixed width of either 10 or 20 km s⁻¹. For each realisation of the 1D CCF, this therefore yields 15 randomised measurements, and the experiment is repeated 5000 times to yield 75 000 realisations of random Gaussian fits to the CCF of each species. The resulting distributions are plotted in Fig. D.1. These figures show that the detected line strengths are significantly stronger than expected from random fluctuations occuring in the CCFs at radial velocities far away from the rest-frame of the system.

Fig. C.1 Distributions of randomly generated instances of the 1D CCF at the rest-frame velocity of the planet. The distributions are constructed by considering random subsets of in- or out-of-transit CCFs, for each of the detected species reported in Sect. 3.2. The detected line amplitude is expected to be zero if out-of-transit CCFs are considered (out-out in orange and if in-transit CCFs are normalised with the master in-transit spectrum (in-in, blue). Only when the in-transit CCFs are normalised by the out-of-transit CCFs (in-out, green) is the measured line amplitude expected to deviate from zero; significantly so for the six species reported.

Appendix D Cross-correlation functions and bootstrap results

Fig. D.1 All cross-correlation functions in Kp-Vsys space (left column), and co-added in the rest-frame of the planet (middle column). The blue shaded area indicates the 1-σ error expected from photon-noise. Right column: distributions of Gaussian fits to random realisations of the 1D-cross-correlation function generated by stacking the time-series of cross-correlation functions with random shifts applied to them. Gaussian functions are fitted with fixed widths of 10 km s−1 (blue) and 20 km s−1 (orange), to match the typical width retrieved from real detections. These distributions serve toillustrate the detection significance and sensitivities for each of the tested species, given the real noise properties of the cross-correlation functions. The blue distributions (10 km s−1) are wider than the orange distributions (20 km s−1) because spurious fluctuations that are wider are less likely to occur at random. Red lines indicate the line-strengths of the detectedspecies.

Fig. D.1 continued.

Fig. D.1 continued.

Fig. D.1 continued.

Fig. D.1 continued.

References

All Tables

All Figures

	Fig. 1 S/N at the centre of order 56 (near the location of the Na I doublet) in all three nights as a function of orbital phase. The vertical lines indicate the start and end of the transit. The second night, in blue, shows markedly lower S/N during transit. For this reason, it is rejected during the focused analysis of the Na doublet.
In the text

	Fig. 2 Abundance profiles of selected species at a temperature of 2000 K and 20× solar metallicity, as computed by FastChem. Solid lines correspond to atomic species, dashed lines to molecules, and dashed-dotted lines to atomic hydrogen and H−. Transmission spectroscopy is principally sensitive to pressures below the milli-bar level (e.g. Kitzmann et al. 2018).
In the text

Fig. 3 Three of the four models of the transmission spectrum of WASP-121 b assuming chemical equilibrium, a metallicity of 20× solar, a temperature of 1500 (blue) or 2000 K, with (yellow) and without (green) TiO opacity, used for model injection and comparison. Inset panels: models with (right) and without (left) the contribution of TiO opacity near 517 and 622 nm, both atthe high native resolution of these models (light colour), as well as broadened to include the effects of the instrumental resolution and rigid body rotation assumed for the planet (dark colour). VO absorption bands are only evident when TiO opacity is removed (the bands that remain visible in the green and blue models near, for example, 550 and 575nm), which would otherwise be masked by much stronger TiO bands. The fourth model (at 3000 K) is not plotted here, but is shown separately in Fig. 12.

In the text

Fig. 4 Cross-correlation templates introduced in this work prior to continuum subtraction and broadening, derived from model spectra of the atmosphere of WASP-121 b. The chemical abundances of the atoms assume chemical equilibrium with elemental abundances equivalent to a metallicity of 20× solar and a temperature of 2000 K. The templates include all atoms with atomic number up to 40 for which line absorption is significant in the HARPS wavelength range (387.4–690.9 nm) under these conditions, as well as the H2O, SH, TiO and VO molecules. SH is included as a byproduct of sulphur photo-chemistry, with a disequilibrium abundance profile artifically set to match the profile presented in Evans et al. (2018). The profile is parametrised as having a VMR of 10−4 for pressureshigher than 10−5 bar, and zero at lower pressures. Note that for clarity, the scale of the y-axis of the molecular absorption spectra is a factor of 2 smaller than the atoms. Also note that the absorption lines of SH are expected to be significantly weaker than those of the atomic metals. Templates of ionic species (Ti II, Cr II and Fe II) are adopted from the analysis of Hoeijmakers et al. (2019), calculated for the atmosphere of KELT-9 b with a temperature of 4000 K and solar metallicity.

In the text

Fig. 5 Subtraction of the Doppler shadow and emergence of the signature of the exoplanet atmosphere, following the radial-velocity curve indicated by the white line. Top panel: raw cross-correlation matrix of the first night of observationswith the Fe I template. During the transit, the Doppler shadow emerges as the near-vertical structure around 40 km s−1, because the planet is on a strongly misaligned orbit (Delrez et al. 2016; Bourrier et al. 2020b). Middle panel: best-fit model of the Doppler Shadow described as a two-component Gaussian. Bottom panel: residuals after subtracting the best-fit model from the raw cross-correlation matrix, revealing the signature of the planet atmosphere. Residuals remain in the stellar line core before and after the transit event, but these do not affect the signature of the planet atmosphere, as the latteris constructed by considering in-transit spectra only. The signature of the planet atmosphere appears as the dark slantedfeature indicated by the white line. Note that although the sign of the absorption is negative in this figure, the sign flipped further in the analysis to denote absorption, notably in Figs. 8 and 9.

In the text

Fig. 6 Transmission spectrum of WASP-121 b at the location of the Na I D-doublet. The spectrum was constructed from the observations obtained during night one and night three, as night two was ruled out as a spurious signal. Top panel: transmissionspectra at the native resolution of HARPS in grey, binned by 20 for better visibility in black. The transmission spectrum corresponding to the best-fit analytical model (see Sects. 3.1.2 and Appendix A) is shown in red, along with the residual of the data after subtraction in the bottom panel. Blue, green and yellow lines depict two forward-models based on a radially outflowing wind (Seidel et al. 2020), an optically thin torus (Oza et al. 2019; Gebek & Oza 2020) and a hydrodynamically escaping envelope (Wyttenbach et al. 2020). The latter two reproduce the significant difference in depth between the D1 and D2 lines.

In the text

Fig. 7 Distributions of the bootstrap analysis for the 12 Åcentral passband. As expected the “in-in” (red) and “out-out” (blue) distributions are centred around zero (no planetary detection) and the randomised “in-out” distribution shows a detection (black) for both night one and night three. Night two shows wide variability in all three scenarios, with the “in-out” distribution also centred at 0.

In the text

Fig. 8 Summary of co-added cross-correlation functions in velocity–velocity space (top panels) and in the planet rest-frame (bottom panels), for various species. Dashed lines indicate the velocities at which the signature of the planet is expected to occur, given the known systemic velocity and orbital parameters of the planet. Red and yellow lines are cross-correlation functions obtained by injecting two different models into the data prior to cross-correlation. The cross-correlation functions of all surveyed species are provided in Fig. D.1. Note that although the sign of the absorption is positive in this figure, the sign was flipped earlier in the analysis to denote the flux received from the system, notably in Figs. 5 and 6. Also note that the scaling of the vertical axis of the cross-correlation function of Mg I is multiplied by a factor of three, because the line of Mg I is significantly deeper than the other atoms.

In the text

	Fig. 9 Similar to Fig. 8, but for molecular species. No molecules are detected using the cross-correlation technique. The y-axis values of the 1D cross-correlation functions of H2O and SH are increased by factors of 2 and 10 respectively, to accommodate for their larger scatter. The injected model of SH (green) includes only continuum opacity and line opacity by SH at a temperature of 2000 K.
In the text

	Fig. 10 Similar to Fig. 8, but for ions. No ions are detected using the cross-correlation technique.
In the text

Fig. 11 Abundance profiles of all Ti and V-bearing species included in FastChem, at a temperature of 2000 K and a metallicy of 20× solar. All solid lines correspond to Ti-bearing species, while dashed lines correspond to V-bearing species. This figure supports theconclusions that Ti and V are expected to have equal abundance; that both mono-oxides as well as dioxides are important constituents of the atmosphere, as well as metal sulphites; and that Ti II and V II are not important in chemical equilibrium.

In the text

Fig. 12 STIS/G430L transmission spectrum of WASP-121 b as published by Evans et al. (2018) indicated with open circles. A model spectrum assuming chemical equilibrium, a temperature of 3000 K, 20× solar metallicity and SH computed as described in Sect. 2.3.2 and Table 2 is overplotted, both at the native resolution (light-blue) and binned to the resolution of STIS (dark-blue line and circles). The wavelength bin at 395 nm coincides with the Al I doublet, at which location the transit radius is strongly over-predicted compared to the STIS observations. This models predicts that the cores of the Al I doublet become optically thick at transit radii of ~ 5 × 10−3 above the local continuum, which corresponds to an altitude of ~0.3Rp above the planet. At such low pressures, photo-ionisation is expected to be important (e.g. Moses et al. 2011; Kitzmann et al. 2018). At the same time, Al I has one of the lowest ionisation thresholds of the elements, so we predict that equilibrium models may strongly over-predict the strength of the Al I doublet at this wavelength.

In the text

	Fig. A.1 Schematic of the geometry of an optically thin cloud with projected area Ac surrounding an optically thick planet with radius Rp, both transiting a star with radius R*, with an effective projected area A*. Equation (A.2) consists of the sum of the covered and uncovered components, Fcovered and Funcovered respectively.
In the text

Fig. B.1 Pre-processing steps applied to obtain cross-correlation functions in the rest-frame of the planet, as described in Sect. 2.3.3. Panel A: an example spectral order of night one as extracted by the HARPS DRS, plotted as a time-series. Panel B: after normalisation by the time-average flux and colour correction. Panel C: cross-correlation template of Fe I after continuum-subtraction, with which the spectra are cross-correlated. Panel D: average flux in the spectral order, that is divided out of the spectra in panel A to obtain the spectra in panel B. Panel E: cross-correlation function of Fe I of the entire time series. Panel F: after removal of the time-average of the out-of-transit cross-correlation functions. Panel G: after removal of the Doppler shadow and application of the high-pass filter. Panel H: after shifting the cross-correlations to the rest-frame of the planet. Panel I: after weighing the cross-correlation functions by the average flux of the time-series (panel D) and selecting only in-transit observations. Panel J: one dimensional cross-correlation function in the rest-frame of the planet, after averaging the cross-correlation functions (panel I) in time.

In the text

Fig. C.1 Distributions of randomly generated instances of the 1D CCF at the rest-frame velocity of the planet. The distributions are constructed by considering random subsets of in- or out-of-transit CCFs, for each of the detected species reported in Sect. 3.2. The detected line amplitude is expected to be zero if out-of-transit CCFs are considered (out-out in orange and if in-transit CCFs are normalised with the master in-transit spectrum (in-in, blue). Only when the in-transit CCFs are normalised by the out-of-transit CCFs (in-out, green) is the measured line amplitude expected to deviate from zero; significantly so for the six species reported.

In the text

Fig. D.1 All cross-correlation functions in Kp-Vsys space (left column), and co-added in the rest-frame of the planet (middle column). The blue shaded area indicates the 1-σ error expected from photon-noise. Right column: distributions of Gaussian fits to random realisations of the 1D-cross-correlation function generated by stacking the time-series of cross-correlation functions with random shifts applied to them. Gaussian functions are fitted with fixed widths of 10 km s−1 (blue) and 20 km s−1 (orange), to match the typical width retrieved from real detections. These distributions serve toillustrate the detection significance and sensitivities for each of the tested species, given the real noise properties of the cross-correlation functions. The blue distributions (10 km s−1) are wider than the orange distributions (20 km s−1) because spurious fluctuations that are wider are less likely to occur at random. Red lines indicate the line-strengths of the detectedspecies.

In the text

	Fig. D.1 continued.
In the text

	Fig. D.1 continued.
In the text

	Fig. D.1 continued.
In the text

	Fig. D.1 continued.
In the text