© The Authors 2024

1 Introduction

Ground-based high resolution observations (R ≥ 25 000) of exoplanets have recently advanced our ability to constrain the chemical, physical, and dynamical properties of their atmospheres. These observations exploit the presence of thousands of individual planetary absorption and emission lines within a given wavelength range. This sensitivity enables us to unambiguously identify atmospheric species using techniques such as narrowband transmission spectroscopy or the cross-correlation (CC) method (Snellen et al. 2010). In addition to these powerful methods, advanced techniques at high resolution, such as CC-to-logL mapping or direct likelihood evaluation (Brogi & Line 2019; Gibson et al. 2020), have now allowed us to go beyond simply identifying species to constraining the abundance of atmospheric species, vertical temperature-pressure (T-P) profiles, and dynamical properties via direct atmospheric retrievals for both emission and transmission spectra. Multiple recent studies have successfully utilised these techniques to constrain the relative abundances of atomic and molecular species, enabling precise constraints on atmospheric C/O ratios and metallicities at high resolution (Line et al. 2021; Brogi et al. 2023; Boucher et al. 2023; Ramkumar et al. 2023; Smith et al. 2024). These constraints may offer valuable insights into how these complex objects formed and evolved to their current atmospheric composition and orbital architecture.

These techniques are particularly tailored to close-in giant planets, such as the ultra-hot Jupiters (UHJs). Their extreme proximity to their host stars inevitably ensures a relatively large planet-to-star radius ratio, as well as a short orbital period, which is crucial for the feasibility of observations of transiting systems. Also, under such extreme conditions, the planet’s constituent atomic and molecular species undergo thermal dissociation and ionisation (Parmentier et al. 2018) such that heavier refractory species are readily observable as they have not condensed out of the atmosphere. This makes optical high resolution spectroscopy particularly powerful as these refractory elements possess a large number of absorption and emission lines in this wavelength regime. Strong lateral winds and vertical mixing are also thought to prevent the cold-trapping of condensates in the lower atmosphere of the cooler nightside, circulating them to the hotter dayside where they are vaporised (Mikal-Evans et al. 2022). Furthermore, as these species have not yet condensed out of the gaseous phase, the observed atmospheric composition and thus metal enrichment of UHJs may allow us to link these observations directly to theoretical models of planet formation and migration (Öberg et al. 2011; Madhusudhan et al. 2014; Mordasini et al. 2016; Mollière et al. 2022; Lothringer et al. 2021). In recent years, a plethora of refractory species have been observed in the atmospheres of multiple hot Jupiters and UHJs (Hoeijmakers et al. 2018; Nugroho et al. 2020b; Pino et al. 2020).

These close-in, tidally locked objects undergo significant rotation during transit, such that we are sensitive to distinct longitudinal regimes of the atmosphere as a function of orbital phase (Wardenier et al. 2021a). The potential variation in the abundance and dynamical properties of these species as a function of orbital phase has opened a new avenue for characterising these objects, by exploring their intrinsic 3D nature. Their tidally locked orbits also ensure a permanently irradiated day-side and a cooler nightside, with drastically differing chemical and dynamical properties; this is expected to drive asymmetries in the chemical abundance and T-P profile across the limbs of the planet. This effect was first observed in the atmosphere of WASP-76b, in which an apparent ‘kink’ in the CC trail of Fe was observed, hinting at the existence of an asymmetric Fe abundance in the atmosphere across the limbs of the planet, as different regions rotated into and out of view (Ehrenreich et al. 2020). However, Savel et al. (2022) have shown that this apparent discrepancy in the Fe abundance across the limbs of the planet is likely explained by a global effect, such as asymmetric cloud formation or a longitudinal temperature asymmetry. This hypothesis is further supported by the presence of a similar kink in the CC trails of similar neutral atomic species in the atmosphere of WASP-76b (Kesseli et al. 2022; Pelletier et al. 2023). These species form at similar altitudes as Fe, suggesting that a global process is likely influencing their common apparent radial velocity offsets.

Similar asymmetries have also been found in other UHJs. Prinoth et al. (2022) found absorption signals from a range of neutral and ionised metals in the atmosphere of WASP-189b, with line positions that vary significantly between species, suggesting an inhomogeneity in their spatial distributions. Offsets in velocity between species were also detected in the atmosphere of WASP-121b; this is particularly true for ionised calcium (Maguire et al. 2023), which is thought to exist as part of an extended outflowing envelope and as such is present in a different dynamical regime than the aforementioned neutral species. This outflowing envelope of ionised species was further supported by the existence of artificially ‘scaled-up’ Fe II features relative to Fe, suggesting the model assumption of hydrostatic equilibrium had broken down for Fe II. Further analysis of this dataset by Seidel et al. (2023) found a potential asymmetric super-rotating Na jet, which stems from the planet’s hotter dayside post mid-transit. However, as this is only a partial transit of WASP-121b, further analysis is required to rule out its existence as part of a global day-to-night flow.

This varying longitudinal representation of an exoplanetary atmosphere has recently been modelled directly at high resolution with the combination of advanced rotational broadening kernels (Gandhi et al. 2022; Boucher et al. 2023) and distinct transit light curves for the leading morning limb and trailing evening limb. As a planet transits its host star, the morning limb passes the stellar disc before the evening limb. For tidally locked planets, the morning limb is permanently redshifted, whereas the evening limb is permanently blueshifted. This process, however, is further complicated by the existence of complex wind patterns, such as super-rotational equatorial jets or radially outflowing winds, which are often degenerate with planetary rotation. We can use this information to separate the morning and evening contributions of planetary absorption lines in velocity space and separately model each limb, and thus constrain the atmospheric properties (abundances, T-P profiles, dynamics, etc.) of the morning and evening limbs separately.

Gandhi et al. (2023, hereafter G2023) recently utilised this technique to constrain the atmospheric abundance of multiple metals across the limbs of WASP-76b. The abundances were found to be broadly consistent longitudinally, supporting the aforementioned hypothesis presented by Savel et al. (2022).

Pelletier et al. (2023, hereafter P2023) recently utilised MAROON-X (M dwarf Advanced Radial velocity Observer Of Neighboring eXoplanets) observations to retrieve a broadly stellar composition of refractory elements for WASP-76b. However, they find substellar abundances relative to Fe for alkali metals, as well as V and Cr, suggesting that a combination of chemical processes is depleting the abundance of several species. They also detect VO in the atmosphere of WASP-76b with both ESPRESSO (Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations) and MAROON-X observations, which has long been theorised as a thermal inversion agent that drives high-altitude temperature inversions in the upper atmospheres of UHJs but has until now remained elusive.

Similarly to the spatially resolved analysis of exoplane-tary atmospheres at high resolution, first presented by Gandhi et al. (2022), it has been demonstrated that the distinct transit light curves caused by asymmetric limbs enable the extraction of separate transmission spectra for morning and evening limbs directly from low resolution transit time-series measurements (Espinoza & Jones 2021; Grant & Wakeford 2022). These updated modelling approaches will enable the characterisation of the 3D nature of these objects at both low and high resolution directly from transmission spectroscopy observations. Coupled with the enhanced sensitivity offered by the next-generation space telescopes and ultra-stable ground-based echelle spectro-graphs, this will also provide a more detailed understanding of the intricate chemical and physical processes at play in their atmospheres.

In Sect. 2 we present our ESPRESSO observations and data reduction, the removal of spectral distortions in the form of centre-to-limb variations (CLVs) and the Rossiter-McLaughlin (RM) effect, and the removal of stellar and telluric contamination with SYSREM. In Sect. 3 we detail our forward model atmosphere and our rotational broadening kernel, and demonstrate our analysis methods. These methods include the CC technique, as well as outlining our full atmospheric retrieval framework. In Sect. 4 we discuss our findings before summarising our results and conclusions in Sect. 5.

2 ESPRESSO observations

We observed one full primary transit of WASP-76b with the ESPRESSO spectrograph (Pepe et al. 2021) at the Very Large Telescope (VLT) on 2019 October 18, hereafter referred to as T3. In addition to this observation, we included two further archival primary transits, observed on 2018 September 2 and 2018 October 30, respectively. These archival transits of WASP-76b were first analysed by Ehrenreich et al. (2020), and have been analysed extensively in the subsequent years (Tabernero et al. 2021; Seidel et al. 2021; Gandhi et al. 2022, Gandhi et al. 2023; Kesseli et al. 2022; Pelletier et al. 2023). They are hereafter referred to as T1 and T2, respectively. Further details regarding these observations are given in Table 1. The data were reduced using the ESPRESSO Data Reduction Software (DRS) v3.1.0, and the non-blaze-corrected S2D spectra were extracted.

2.1 Cleaning and blaze correction

The S2D spectra were then cleaned and blaze corrected following Merritt et al. (2020) to which we refer the reader. Similarly to Maguire et al. (2023), the residual array in our blaze correction was smoothed with a narrower filter than in our previous works (a median filter with a width of 101 pixels and a Gaussian filter with a standard deviation of 200 pixels) in order to minimise the effects caused by the ‘wiggle’ pattern seen in ESPRESSO spectra (Allart et al. 2020). To ensure these wiggles – alongside the relatively aggressive preprocessing steps and stellar and telluric removal procedures – did not inadvertently distort our underlying exoplanetary signal, we conducted extensive injection and recovery tests (see Sect. 3.3 for further details).

2.2 Removal of CLVs and the RM effect

In order to optimise the information content of these datasets, it was imperative that we preprocessed the data appropriately, removing any distortions caused by stellar and telluric contamination, as well as secondary effects. With transit observations taken with an ultra-stable high resolution spectrograph such as ESPRESSO, we are also sensitive to distortions in stellar line shapes in transit caused by CLVs and the RM effect. As has become common practice in the reduction of high resolution transmission spectroscopy observations (Yan et al. 2017, 2019; Casasayas-Barris et al. 2017, 2018; Turner et al. 2020; Nugroho et al. 2020b; Maguire et al. 2023), we modelled and removed these distortions of the stellar spectral lines, caused by both CLVs and the RM effect, from the data, assuming stellar parameters outlined in Table B.1. It is worth noting that the efficacy of this procedure relies heavily on a precise estimate of the projected spin-orbit angle, λ, which has a relatively large uncertainty. A detailed description of this procedure can be found in Maguire et al. (2023).

2.3 Removal of stellar and telluric contamination

In addition to the in-transit distortions caused by CLVs and the RM effect, the contamination caused by stellar and telluric lines was also removed. Firstly, we shifted each spectrum into the stellar rest frame via the stellar radial velocity, computed via the orbital parameters outlined in Table B.1. This ensured that the stellar and telluric lines are (quasi-)stationary in time. These (quasi-)stationary features were then fit with the iterative algorithm SYSREM, first applied to high resolution spectroscopy by Birkby et al. (2013) and has now become a standard practice in the field (e.g. Nugroho et al. 2017, 2020a,b; Birkby et al. 2017; Gibson et al. 2019, 2020, 2022; Merritt et al. 2020; Serindag et al. 2021; Maguire et al. 2023; Ramkumar et al. 2023). Before applying SYSREM, each order was divided by its median spectrum, after which the SYSREM model is then subtracted from the data. This procedure is akin to fitting the data directly with SYSREM and dividing through by the resultant model; however, subtraction allows for a faster model filtering procedure whilst retaining the accuracy required for our analysis (Gibson et al. 2022). This model filtering procedure then in turn replicates any distortions caused by SYSREM on our atmospheric model (see Fig. 1).

In order to determine the ‘optimum’ number of SYSREM iterations for each dataset, an atmospheric model, including all species outlined in Sect. 3.1, was injected with a negative K_p value, where K_p is the planet’s radial velocity semi-amplitude, such that it was well separated from the real exoplanet signal. We then performed a likelihood mapping analysis, and the number of iterations that resulted in the optimum detection significance was chosen. In order to compute the detection significance, we followed Gibson et al. (2020), and computed a log-likelihood distribution (see Sect. 3.4) as a function of the model scale factor, α, conditioned on the optimumβ K_p, and Δv_sys values, where β and Δv_sys are the noise scale factor and the systemic velocity offset, respectively. The expected value, or mean, of this distribution was computed, and divided by its standard deviation, resulting in a detection significance. This, in essence, computes the number of standard deviations the maximum signal deviates from an alpha of 0 (i.e. a non-detection). The optimum number of SYSREM iterations for each dataset was 16, 18, and 22 for T1, T2, and T3, respectively. It is worth noting that this procedure is highly sensitive to the injected model, as a model with many absorption lines in orders dominated by stellar and telluric lines will, in theory, require a higher number of iterations to optimally extract the exoplanet signal. Furthermore, separating each order and subsequently determining the optimum number of SYSREM iterations per order would also minimise the inevitable distortion caused by SYSREM on the exoplanet’s signal (Nugroho et al. 2017).

3 Analysis

3.1 Forward model

For each of the transits, we performed traditional CC analysis, as well as a full atmospheric retrieval, with 1D atmospheric models computed using IRRADIATOR (Gibson et al. 2022), including cross-sections from a number of atomic and molecular species (Fe, V, Cr, Mg, Na, and VO). The cross-sections for Fe, V, Cr, and Mg were obtained via the Kurucz line list (Kurucz & Bell 1995). Similarly, the Na cross-sections were obtained via the VALD (Vienna Atomic Line Database) line list (Piskunov et al. 1995), whereas the VO cross-sections were obtained via the Exo-Mol line list (VOMYT; McKemmish et al. 2016). We defined 70 atmospheric layers, equidistant in log space, ranging from 10² bar to 10⁻¹² bar, across which we computed a T-P profile using the parametric model from Guillot (2010) with input parameters T_irr, T_int, ĸ_IR, and γ, corresponding to the irradiation temperature, internal temperature, mean infrared opacity, and the ratio of visible-to-infrared opacity, respectively. After solving for hydrostatic equilibrium, with variable gravity, the opacity of each predefined layer of our atmosphere was computed by weighing each species’ cross-section on their respective volume mixing ratios (VMRs), χ_species, assuming a constant VMR with altitude. In addition, we included a H₂ VMR, χ_ray, along with H₂ opacities from Dalgarno & Williams (1962), to account for Rayleigh scattering, as well as a further continuum opacity source caused by an opaque cloud deck at a pressure level P_cloud, below which the transmission spectrum model is truncated. In practice, we fitted for the logarithms of strictly positive parameters, with the exception of temperature. Therefore, in total, our 1D transmission spectrum model has 6 + N_species input parameters (T_irr, T_int, log₁₀(ĸ_IR), log₁₀(γ), log₁₀(P_cloud), log₁₀(χ_ray), log₁₀(χ_species)).

Fig. 1 2D Na model convolved with the broadening kernel outlined in Sect. 3.2. Upper left: rotational broadening kernel outlined in Eq. (1). The solid lines red and blue curves show the final morning and evening halves, respectively, after convolution with an instrumental profile, shown in grey. Upper right: asymmetric Na model, with varying broadening and atmospheric parameters between the morning and evening limbs. The combined model is overplotted with a dashed black line. Middle: 1D model projected to a 2D orbital phase vs wavelength array for a single ESPRESSO order. This array is then weighted by the limb-darkening model described in Gandhi et al. (2022). This ensures solely morning or evening contributions during ingress or egress, respectively. Lower: 2D orbital phase vs wavelength array after SYSREM model filtering, highlighting the distortions of the planetary lines after preprocessing. The dotted black curve shows the planetary velocity to which the model is projected.

3.2 Broadening kernel

Recent studies have modelled the varying longitudinal representation of an exoplanetary atmosphere directly with the use of realistic rotational broadening kernels (Gandhi et al. 2022; Boucher et al. 2023). At high resolution, we are not only sensitive to the positions of individual absorption lines, but also the subtle variations in their line shapes caused by winds and rotation (Seidel et al. 2020; Hoeijmakers et al. 2020). Furthermore, with transmission spectroscopy, we are only sensitive to a thin annulus of the planet’s atmosphere, namely the terminator. For a rotating annulus of inner radius r₁ and outer radius r₂ = r₁ + (d • r₁), the broadening profile can be modelled by a kernel of the form

$$G(\text{Δ}v)=\{\begin{array}{ll}\frac{a}{d}\hfill & ,\text{ for }{r}_1\le \left|x\right|<r_2\hfill \\ \frac{a-\sqrt{1-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}}{d}\hfill & ,\text{ for }\left|x\right|<r_1\hfill \end{array}$$(1)

where $a=\sqrt{{(1+d)}^2-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}$, d is the thickness of the annulus relative to r₁, Δv is the velocity shift of the line at a position x, and Δv_rot is the shift due to the equatorial velocity (see Appendix A for the full derivation).

This broadening kernel allows us to effectively separate the morning (leading) half of the terminator from the evening (trailing) half of the terminator in velocity, via subtle variations in line shape. In practice, we modelled each half of the terminator (limb) separately with its own 1D atmospheric model and convolved each transmission spectrum with one half of the above broadening kernel. The effect of this convolution is shown in Fig. 1. This procedure allows us to retrieve separate atmospheric parameters, as well as separate rotational velocities Δv_rot and thickness d, for each atmospheric limb. This results in four parameters for our kernel, two for each limb: d_m, d_e, Δv_rot,m, and Δv_rot,e. The subscripts m and e represent the morning and evening limbs, respectively.

In practice, we fixed d_m and d_e to a constant value equivalent to 5H_s following Boucher et al. (2023), where H_s is the atmospheric scale height. We also ensured each half of the kernel normalises to 0.5, such that each 1D model contributed to half the total absorption. This also ensured that any asymmetry in line or transit depth retrieved from the data is caused solely by physical effects, such as an increase in scale height due to a temperature asymmetry between the morning and evening limbs. This asymmetry, or lack thereof, can then be interpreted via varying physical parameters retrieved for the limbs’ individual atmospheric models.

After convolution with our broadening kernels, the independent atmospheric models were summed and linearly interpolated to the 2D wavelength grid of our data (order × wavelength). For each order, the 2D forward model was then Doppler shifted to a planetary velocity given by:

$$v_{\text{p}}=K_{\text{p}}\sin (2\pi \cdot \varphi )+\text{Δ}{v}_{\text{sys}}+v_{\text{bary}}$$(2)

where ϕ and υ_bary are the orbital phase and barycentric velocity, respectively. The wavelength solutions have already been corrected for the barycentric velocity by the ESPRESSO DRS. This resulted in a 3D shifted forward model with the same dimensions as our data (time/phase × order × wavelength). The middle panel of Fig. 1 shows a single order of our model after this projection. The forward model was then filtered using the SYSREM basis vectors via the procedure outlined in Gibson et al. (2022). The result of this model filtering is shown in the bottom panel of Fig. 1.

3.3 Injection and recovery tests

Before conducting a full retrieval using our atmospheric model, we had to conduct injection and recovery tests. These tests ensure that the preprocessing steps applied in Sects. 2.1–2.3 do not inadvertently distort the underlying exoplanetary signal, biasing our retrieved constraints and subsequent interpretations. Furthermore, as we aimed to constrain any asymmetric signals stemming from the morning or evening limbs, we had to ensure that our retrieval framework is sensitive to subtle variations in line shapes, which are already deep beneath the noise of our observations. Thus, we are typically unable to measure individual line shapes directly in order to observe these variations. For each of our datasets, we injected an asymmetrically broadened atmospheric model into a subset of ESPRESSO orders¹, prior to the cleaning and blaze correction steps outlined in Sect. 2.1. This injected model included cross-sections from Fe, Na, V, and VO. We then performed our retrieval analysis on the injected data, similar to that outlined in Sect. 3.4, and found the retrieved model parameters, T-P profiles, and broadening kernels to be in agreement with those injected (see Figs. C.1–C.3).

3.4 Atmospheric retrieval

To perform a full atmospheric retrieval analysis on each of our datasets, we had to compute a forward model, m_i, for a given set of model parameters, θ, outlined in Table 2. Assuming uniform priors, the log-posterior is computed by adding the log-prior to the log-likelihood:

$$\ln ℒ=-N\ln \beta -\frac{1}{2}{\chi }^2$$(3)

given

$${\chi }^2={\displaystyle \sum _{i=1}^N\frac{{\left(f_i-\alpha m_i(\theta )\right)}^2}{{\left(\beta {\sigma }_i\right)}^2}}$$(4)

where f_i is the data, ${\sigma }_i^2$ is the data variance, α is the model scale factor, and β is the noise scale factor (see Gibson et al. 2020 for a detailed derivation).

This was then folded into a Markov chain Monte Carlo (MCMC) framework in order to retrieve posterior distributions for each of the model parameters, for each of the datasets. We use a differential-evolution Markov chain (DEMC) (Ter Braak 2006; Eastman et al. 2013), running an MCMC chain with 128 walkers, with a burn-in length of 300 and a chain length of 600, unless stated otherwise. This results in 115 200 samples of the posterior, of which 38 400 are discarded. The chains’ convergence is assessed via the Gelman-Rubin statistic, after splitting each chain into four independent sub-chains. We also split the chains into groups of independent walkers, and overplot the 1D and 2D marginal distributions, in order to visualise the convergence of our MCMC. The results of each of our retrievals are shown in Figs. D.3-D.5.

We compared the retrieved absolute abundances across each of the transits for the morning (Fig. D.1) and evening (Fig. D.2) limbs. These results were then compared with stellar values, for species with measured stellar abundances (Tabernero et al. 2021), as well as with solar abundances, where measured (Asplund et al. 2009).

Our morning and evening abundance constraints are also compared with the constraints presented by G2023. Due to the relatively large rotation angle of close-in planets such as WASP-76b during transit, G2023 opt to retrieve varying regions of the morning and evening limbs before and after mid-transit (see Gandhi et al. 2022 for furtherdetails). We compared oursimulta-neous retrieval of both limbs with the weighted mean of regions A and B (evening) and C and D (morning) from G2023. It is also worth noting that the constraints obtained by G2023 are from a joint retrieval of T1 and T2, whereas we opted to fit both datasets individually.

Similarly to previous studies, we also compared the relative abundances of our atmospheric species, and find these values to be better constrained (see Figs. 2 and 3). This is due to degeneracies in model transmission spectra, as well as the normalisation of high resolution spectra during preprocessing; as a result we are no longer sensitive to the true continuum of the planet (Birkby 2018).

Fig. 2 Relative abundance comparison across each of the three transits, for the morning limb. The stellar values from Tabernero et al. (2021) are shown as a vertical grey line, with their 1σ contours shaded where measured. The solar values from Asplund et al. (2009) are shown as a vertical orange line, with their 1σ contours shaded. The weighted mean of the retrieved relative abundances in regions C and D (morning) from G2023 are also shown, where measured. Otherwise, the relative abundances measured for the whole terminator from P2023 are shown.

Fig. 3 Relative abundance comparison across each of the three transits, for the evening limb. Shaded regions are as in Fig. 2. The weighted mean of the retrieved relative abundances in regions A and B (evening) from G2023 are also shown, where measured.

Fig. 4 Kp – Δvsys maps of V for T1, T2, and T3. Upper: the V atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals but are un-broadened. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values. Lower: comparison of the CCF of Fe and V, at their optimum Kp value. As the models were un-broadened, this distribution can be used as a proxy for the average line shape of each species in our data, highlighting any apparent asymmetries in the species’ longitudinal distributions.

3.5 Cross-correlation

We confirmed the presence of the species included in our atmospheric model individually via traditional CC analysis. The atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals. The results are presented in Appendix E. We confirm the novel detection of VO in the atmosphere of WASP-76b by P2023 in datasets T2 and T3, independently detecting VO in T3 (see Fig. E.1). As this analysis is the first analysis of the T3 dataset, we also aimed to investigate the apparent kink in the CC trail first shown by Ehrenreich et al. (2020) for T1 and T2 (see Fig. E.3). This apparent kink has previously been independently confirmed via reanalysis of these transits (Kesseli et al. 2022), as well as additional analysis from subsequent observations using HARPS (Kesseli & Snellen 2021) and MAROON-X (P2023).

Traditional CC analysis also allows us to investigate any apparent asymmetry in the abundance of individual species between the morning and evening limbs. We first cross-correlated our data with an un-broadened model, Doppler shifted over a range of Δv_sys values:

$$\mathrm{CCF}\left(\text{Δ}{v}_{\text{sys}}\right)={\displaystyle \sum _{i=1}^N\frac{f_i{m}_i\left(\text{Δ}{v}_{\text{sys }}\right)}{{\sigma }_i^2}}.$$(5)

The cross-correlation function (CCF) for each orbital phase is then shifted to a new planetary velocity given by Eq. (2) for a range of K_p values, and integrated as a function of time or phase. We then plotted the CCF at the optimum K_p value, an example of which is shown in Fig. 4. As the model is un-broadened (i.e. it is not convolved with the broadening kernel, but still includes its natural line widths from thermal broadening) the distribution of the CCF can be used as a proxy for the average line shape in our data, after rotational and instrumental broadening. In other words, the CC is a convolution of the intrinsic line shape in the data with the un-broadened model. By comparing the shape of the CCF of Fe, which was previously detected in both limbs (Gandhi et al. 2022, 2023), with that of other species, we can empirically search for any apparent asymmetry in their CCF as a function of Δυ_sys. This asymmetric CCF could potentially be caused by an asymmetric distribution in their abundance between the blueshifted evening limb and redshifted morning limb (Wardenier et al. 2023). CC trails that deviate from the planet’s Keplerian velocity curve, such as that first seen in Ehrenreich et al. (2020), can also often result in a ‘diamond’ feature in the K_p – Δυ_sys maps, as a warped CC trail summed in time will result in multiple peaks in the map localised at the expected K_p and Δυ_sys values (Nugroho et al. 2020b; Wardenier et al. 2021b; Maguire et al. 2023). Figure 4 shows the results of this test for Fe and V, in which we see a narrower CCF for V relative to Fe, for both T1 and T3, caused by a predominantly blueshifted V line shape.

4 Discussion

4.1 Abundance constraints

Our abundance constraints are directly compared with stellar values (Tabernero et al. 2021) where measured (i.e. Fe, Cr, Mg, and Na), as well with solar values (Asplund et al. 2009) where measured. The individual constraints obtained for the morning and evening limbs are also compared with a similar spatially resolved analysis of WASP-76b, conducted by G2023, for specific species (i.e. Fe, V, Mg, and Cr). For species such as Na and VO, where previous spatially resolved constraints are not available, we directly compared our results with those obtained by P2023. It is important to note that, unlike our approach, P2023 modelled the terminator as a whole. The relative abundance constraints for the morning and evening limbs are shown in Figs. 2 and 3, respectively. Whereas, the absolute abundance constraints for the morning and evening limbs are shown in Figs. D.1 and D.2, respectively.

Fe

Across all three transits, we find consistent Fe abundance constraints in both the morning and evening limbs and find the abundances to be consistent across the limbs of the atmosphere for each of the individual transits. This suggests a stable uniform longitudinal distribution of Fe in the atmosphere of WASP-76b, which is consistent with previous spatially resolved constraints by Gandhi et al. (2022, 2023). This uniform distribution of Fe indicates that a global process, such as an uneven cloud deck or temperature asymmetry, is likely responsible for the asymmetric CC trail previously found (Ehrenreich et al. 2020; Savel et al. 2022). This is further supported by the asymmetric CC trails found for a multitude of species with varying condensation temperatures (Kesseli et al. 2022; Pelletier et al. 2023).

Furthermore, both our morning and evening Fe abundance constraints are consistent with stellar values, which suggests a uniform atmospheric composition similar to that of its parent protoplanetary disc. This is consistent with the findings of G2023, in which Fe/H values of a range of UHJs – corrected for the stellar metallicity of each system – were found to be consistent with their host star values. P2023 also found a similar trend for the abundances of refractory metals in the atmosphere of WASP-76b.

V and VO

V is consistent with solar values and G2023 for T2. Similarly, VO is also consistent with constraints obtained by P2023 for T2 in the morning limb. Both V and VO are unconstrained in the morning limb for both T1 and T3, with upper limits consistent with T2.

In the evening limb, all three transits are consistent, aligning with solar values as well as with both G2023 and P2023. Where constrained, we find V and VO abundances to be consistent across the limbs of the planet.

There is an interesting discrepancy between the morning and evening constraints for both V and VO. For T2, we constrained both species across both limbs. However, for T1 and T3, we only constrained the species in the evening limb. This led to the V/Fe ratio in the morning limb for both T1 and T3 (see Fig. 2) being subsolar. Gandhi et al. (2022) found that, overall, the combined signal is dominated by the more blueshifted evening limb, which is in agreement with previous Global Circulation Model (GCM) predictions (Wardenier et al. 2021b). However, signals from irradiated terminator regions (regions A & D of G2023) typically dominate the combined signal in their respective phase ranges due to larger temperatures and, thus, scale heights. G2023 were also unable to constrain V in the dayside morning region (D). These dominant signals result in tighter constraints for these regions, making it surprising that V is no longer constrained. All three of our observations are of similar quality, with comparable S/Ns when adjusted for the varying exposure times. We monitored the integrated water vapour between observations and found no significant trend that could impact the reliability of our data. This is further supported by the similar constraints shown between each of our injection tests (see Figs. C.1–C.3). We also monitored the variation in airmass between observations and found an increase in airmass for T1 and T3 during ingress, which could explain the dampened morning V/VO signal. However, such an effect would likely reduce the quality of our constraints in the morning limb uniformly across all species, which is not the case, and thus a decrease in data quality for T1 and T3 relative to T2, which might explain worse V/VO constraints in the morning limb, seems unlikely. It is possible that both the V and VO signals are dampened in the cooler morning limb for T1 and T3, due to a combination of recombination and condensation at lower temperatures. However, as our upper limits are consistent with those of T2, as well as G2023, we cannot rule out their presence in the morning limb entirely.

The search for the dampened V signal in the morning limb for both T1 and T3 relative to T2 can be approached more empirically. In Fig. 4, we show a slice through the optimal K_p value for both Fe and V. The figure shows CCFs of un-broadened models that peak at similar Δυ_sys values. Fe has been previously detected in both limbs of the planet by Gandhi et al. (2022), as well as in this work. Thus, in order to identify any asymmetries in the line profile of V, we compared the CCFs of Fe and V. For T2, the CCFs of both Fe and V exhibit a similar shape, characterised by a relatively broad profile. This broad shape is due to contributions from both the blueshifted (evening) and redshifted (morning) components of the line profile. In contrast, for both T1 and T3, the CCF of V is notably narrower relative to Fe and is primarily dominated by the blueshifted (evening) component, perhaps indicating a more localised contribution to the line profile.

Although the origin of this asymmetry in the V and VO distributions across the morning and evening limbs for T1 and T3 is unclear, if present, its discrepancy with the well-constrained V and VO abundances measured for both limbs for T2 indicates an atmosphere with clear variable chemistry, in which a combination of recombination and condensation is altering the V and VO abundances at the cooler morning limb over the course of months to years.

Although the VO line list used, VOMYT (McKemmish et al. 2016), has resulted in detections of VO in the atmosphere of WASP-76b (e.g. P2023 and our Fig. E.1), it has previously led to non-detections of VO in the atmosphere of WASP-121b at high resolution (Hoeijmakers et al. 2020; Merritt et al. 2020).

These non-detections were reported as inconclusive, and this line list was shown to be insufficient in recovering injected VO signals from UVES (Ultraviolet and Visual Echelle Spectrograph) transmission spectroscopy data of WASP-121b (de Regt et al. 2022). Recently, Bowesman et al. (2024) presented an updated MARVEL-ised (Measured Active Rotational-Vibrational Energy Levels; Császár et al. 2007) line list for VO (Hy VO), which includes hyperfine splitting effects and should facilitate more reliable high resolution detections and abundance constraints of VO.

We conducted a CC analysis with the updated line list, with cross-sections computed at T = 2800 K and P = 0.1 mbar, shown in Fig. E.2. In order to determine the optimum parameters for this CC analysis, we ran an MCMC similar to that outlined in Sect. 3.4, using a simple 1D Gaussian broadening kernel, and an isothermal T-P profile. We have now detected VO in T1 via CC. However, our previous detection of VO in T2 has now shifted to a significantly lower Δυ_sys and K_p. Our T3 detection significance has increased with the update line list. It is difficult to assess the accuracy of these new detections/non-detections without performing a full atmospheric retrieval, with a varying T-P profile and 2D broadening kernel. As the updated cross-sections are not readily available across our full T-P grid, an atmospheric retrieval akin to that presented above is impracticable at this time. A quantitative assessment of the updated line list is beyond the scope of this work; however, it is required to ensure its reliability for exoplanetary studies at high resolution going forwards. Despite this added uncertainty in terms of a VO detection, our results still show a potentially variable V abundance, with which there is no uncertainty regarding line list accuracy and, alongside this VO analysis, provides an interesting insight into the chemistry of V-bearing species.

Cr

The Cr/Fe constraints obtained are also consistent between transits, as well as with G2023, for both the morning and evening limbs. Similarly to V and VO, however, we find Cr/Fe to be depleted relative to stellar for T1 and T3 in the morning limb, potentially due to partial recombination and condensation of Cr to other compounds, such as CrH, which has previously been detected in the atmosphere of WASP-31b (Braam et al. 2021; Flagg et al. 2023). P2023 also constrained the abundance of CrH in their 1D retrieval of WASP-76b, but were unable to detect it via CC. Both P2023 and G2023 also found substellar Cr/Fe ratios, within 1σ, consistent with our findings.

We compared relative abundances as these constraints are typically more accurate and reliable at high resolution, and in general for transmission spectra. This is due to the well-known degeneracy between the absolute abundance, χ_species, and the opaque cloud deck pressure, P_cloud and/or the reference pressure. Due to this degeneracy, a model with a higher χ_species and a lower P_cloud could, in principle, produce identical features as a model with a lower χ_species and a higher P_cloud (Benneke & Seager 2012). Furthermore, the normalisation of high resolution spectra during preprocessing means we are only sensitive to line strengths above the planetary continuum, making relative line strengths, and consequently, relative abundance constraints, more precise (Birkby 2018; Gibson et al. 2020).

Mg

Similarly to Fe, our Mg abundance constraints are consistent with both stellar values and G2023 across all three transits. They are also consistent across the limbs of the atmosphere for each of the individual transits. Interestingly, the absolute abundance of Mg is consistent with the upper prior limit (see Table 2) for both the morning (T1 and T2) and evening (T2 and T3) limbs. Similar behaviour has been seen previously in high resolution analysis of WASP-121b with UVES and ESPRESSO (Gibson et al. 2022; Maguire et al. 2023). As Mg II has previously been detected in the upper atmosphere of WASP-121b (Sing et al. 2019), it is thought that Mg may be part of an extended, non-hydrostatic envelope, in which case our model assumption of hydrostatic equilibrium breaks down. To account for this, our model would artificially increase the abundance of Mg, increasing its line strengths. It is difficult to say if this is also the case for WASP-76b, as Mg II has not yet been detected in its atmosphere. Furthermore, Mg has relatively few lines in the ESPRESSO wavelength range, which could also explain its poorer abundance constraints.

Na

We find a substellar Na abundance relative to Fe in the evening limb across all three transits. This is in agreement with P2023, who find substellar abundances for alkali metals (Li, Na, and K), suggesting their lower ionisation energies cause a greater degree of ionisation at the temperatures and pressures probed, driving a decrease in the abundance of neutral alkaline atoms. Similar behaviour was found in the atmosphere of WASP-121b, where the abundance of both Na and Ca relative to Fe were typically lower than their stellar counterparts, alongside the presence of high-altitude Ca II indicates a greater degree of ionisation (Maguire et al. 2023). It is also worth noting that G2023 were unable to constrain Na in both the morning and evening limbs of WASP-76b via an independent joint retrieval of T1 and T2.

The morning Na abundances relative to Fe were also found to be substellar for both T1 and T2, consistent with the evening limb and continuing the trend mentioned above. However, for T3, we find a drastically higher Na/Fe ratio driven by an anomalously high Na abundance. It is possible that Na is also part of an extended outer envelope, similar to Mg, becoming more readily ionised in the upper atmosphere, which is accounted for on the morning limb by an artificial increase in abundance. Another possible explanation is the existence of a high-altitude cloud deck on the morning limb, or similar continuum opacity source such as H⁻ opacity, which is not accounted for in our model. This drives the Na abundance, creating a pseudo-continuum. In order to ensure this artificially high Na abundance does not affect the abundances of the other species in our T3 retrieval, the retrieval was re-run without the presence of Na. We find it has a minimal impact on our retrieved abundances but alters our morning T-P profile constraints significantly. This is to be expected, as a higher Na abundance with a smaller uncertainty will result in a lower temperature with tighter constraints.

4.2 Temperature constraints

For T1, our morning and evening temperatures at 0.1 mbar, T_0.1mbar, vary at the 1σ level (see Table 3 and Fig. D.3). These morning and evening temperature constraints are consistent with the spatially resolved temperature constraints of Gandhi et al. (2022). Above the high-altitude morning cloud deck, we find a morning T-P profile that is cooler than that of the evening limb, which is inverted above ~0.1 mbar.

For T2, we find a cooler evening T-P profile that varies from that of the morning limb at the 2σ level at 0.1 mbar (see Table 3 and Fig. D.4). Our evening temperature at 0 1 mbar is consistent with that of Gandhi et al. (2022). For T2, the $T_{0.1\text{mbar}}^{\text{Morn}}$ is significantly higher than that of Gandhi et al. (2022), and at high pressures, it reaches unphysically high values for WASP-76b. This discrepancy between the morning limb with both previous transits presented in this study, as well as with previous studies, is also likely explained by the presence of V and VO in the morning limb, whereas it is unconstrained for T1 and T3, as VO is known to be a thermal inversion agent that drives high-altitude temperature inversions in the upper atmospheres of UHJs. It is worth noting that transmission spectroscopy studies are typically not very sensitive to the T-P structure of the atmosphere, and therefore we must be cautious not to over-interpret our constraints.

For T3 (see Table 3 and Fig. D.5), we find a hotter evening T-P profile that is inverted above ~0.1 mbar, similar to T1. The evening temperature at 0.1 mbar is consistent with that of Gandhi et al. (2022), and varies from the morning limb within 1σ. The $T_{0.1\text{mbar}}^{\text{Morn}}$ is significantly lower than that of Gandhi et al. (2022). This adverse decrease in temperature could be due to our model neglecting the presence of H⁻ opacity, whose presence is known to cool the atmosphere in this pressure regime (Parmentier et al. 2018). As mentioned in Sect. 4.1, this temperature discrepancy is also driven by the presence of an anomalously high Na abundance on the morning limb of T3. When we excluded Na from our retrievals, we found a much broader morning temperature constraint of $T_{0.1\text{mbar}}^{\text{Morn}}={1950}_{-380}^{+990}\text{\hspace{0.17em}K}$, which is consistent with that found in Gandhi et al. (2022). The evening T-P profile does not change significantly.

4.3 Atmospheric dynamics

For T1, we retrieve K_p = 185.48 ± 1.22 km s⁻¹ and Δυ_sys = −7.70 ± 0.82 km s⁻¹, for T2: K_p = 175.48 ± 1.33 km s⁻¹ and Δυ_sys = −8.77 ± 0.44 km s⁻¹, and for T3: K_p = 181.33 ± 2.14 km s⁻¹ and Δυ_sys = −7.88 ± 0.66 km s⁻¹. Although we constrained different values for each of the transits, we computed the planetary orbital velocity, given Eq. (2), using 10 000 random samples from our posterior distributions of K_p and Δυ_sys, and find υ_p consistent within 1σ in transit. There is also a degeneracy between the global systemic velocity offset, Δυ_sys, and the blueshifted component of the rotational velocity, Δυ_rot,e, with the latter ideally accounting for subtle variations in the individual line shapes caused by additional wind patterns and rotation. When accounting for a broadening of an atmospheric model with R = 200 000, we find broadening parameters of υ_rot,m = 1.24 ± 1.16 km s⁻¹ and υ_rot,e = 2.97 ± 1.61 km s⁻¹ for T1. We also determine υ_rot,m and υ_rot,e for T2 of 2.50 ± 1.45 km s⁻¹ and 0.99 ± 0.72 km s⁻¹, respectively. For T3, we find υ_rot,m and υ_rot,e as 5.77 ± 2.62 km s⁻¹ and 1.33 ± 0.94 km s⁻¹. It is important to note that these distributions are non-Gaussian, and therefore, the constraints may not be accurately represented by the median value and standard deviation alone. The combination of these broadening parameters with a constant blueshifted offset of ~8 km s⁻¹ is in agreement with previous high resolution studies of WASP-76b (Ehrenreich et al. 2020; Pelletier et al. 2023; Gandhi et al. 2022, 2023).

The treatment described above also assumes a single K_p, Δυ_sys, Δυ_rot,m, and Δυ_rot,e value governs the dynamics of each of our atmospheric species. In reality, as seen in previous high resolution studies (Prinoth et al. 2022; Wardenier et al. 2023), velocity offsets between species can occur. In reality, each species is likely to exist in different dynamical regimes, particularly when considering differing atmospheric limbs and/or an extended upper atmosphere. A more robust modelling approach would involve the individual treatment of each species’ dynamics (i.e. varying K_p, Δυ_sys, Δυ_rot,m, and Δυ_rot,e for each species), which we leave as an avenue for future research.

Recent studies using low resolution transmission spec-troscopy with JWST have also allowed us to investigate the 3D nature and varying chemistry across the limbs of these objects. For example, Rustamkulov et al. (2023) found a wavelength-dependent central transit time for WASP-39b, suggesting an asymmetric distribution of atmospheric species across the terminator. Furthermore, it has been shown that by probing the varying transit light curves of each atmospheric limb, particularly during ingress and egress, one can retrieve separate transmission spectra for each half of the terminator (Espinoza & Jones 2021; Grant & Wakeford 2022). These techniques at low resolution, coupled with analysis such as that provided in this study, highlight the enormous promise of utilising low and high resolution transmission spectroscopy observations to investigate the 3D nature of close-in exoplanets, a regime that until now was accessible only with phase curve measurements.

5 Conclusions

We have presented multiple high resolution transit observations of the atmosphere of the UHJ WASP-76b, monitored across several months and years, with the ESPRESSO spectrograph at the VLT. We have modelled the morning (leading) and evening (trailing) limbs of the atmosphere separately with the use of realistic rotational broadening kernels. This allowed us to retrieve separate atmospheric abundances, T-P profiles, and dynamical parameters for each atmospheric limb via our atmospheric retrieval framework.

• We have confirmed the presence of VO, recently discovered in the atmosphere of WASP-76b by Pelletier et al. (2023), via CC in both T2 and T3 of our observations, with a detection significance of 4.0 and 7.2σ, respectively.

• We have found a uniform longitudinal distribution of Fe and Mg across the limbs of the atmosphere, for each of the transits, which is consistent with previous works as well as stellar abundances.

• We have constrained Cr/Fe and Na/Fe ratios that are broadly consistent between the limbs of the atmosphere for each of the transits, with the exception of Na/Fe in the morning limb for T3. These abundance ratios were also found to be substellar in both the morning and evening limbs, hinting at a possible depletion of Cr and Na due to recombination and ionisation, respectively, as suggested in previous works.

• The V/Fe and VO/Fe ratios were also found to be consistent between the limbs of the atmosphere for each of the transits, with V/Fe also consistent with solar in the evening limb. V and VO were, however, unconstrained in the morning limb for T1 and T3, suggesting a possible depletion due to recombination and condensation. This is further supported by an apparent asymmetric CCF for V relative to Fe for both T1 and T3.

The consistent results obtained across multiple observations spanning months to years, as well as with previous studies that used different modelling and retrieval frameworks and/or instruments, affirms the efficacy of high resolution atmospheric retrieval frameworks. Spatially resolved temporal atmospheric surveys such as this one, coupled with low resolution efforts in the era of JWST, will allow us to place further constraints on the time-varying dynamics and chemistry of these worlds, allowing us to infer how they evolved to their current climate.

Acknowledgements

This work relied on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere. The author(s) gratefully acknowledge support from Science Foundation Ireland and the Royal Society in the form of a Research Fellows Enhancement Award No. RGF\EA\180315. S.K.N. is supported by JSPS KAKENHI grant No. 22K14092. We are grateful to the developers of the NumPy, SciPy, Matplotlib, iPython, corner, and Astropy packages, which were used extensively in this work (Harris et al. 2020; Virtanen et al. 2020; Hunter 2007; Pérez & Granger 2007; Foreman-Mackey 2016; Astropy Collaboration 2022).

Appendix A Derivation of the broadening kernel

For a uniformly irradiated disc projected by a rigid rotating sphere, the broadening profile G(Δυ) is given by

$$G(\text{Δ}v)=\frac{I_{\lambda }(\text{Δ}v)}{I_{\lambda }(0)}$$(A.1)

where I_λ(0) is the line intensity at the line centre and I_λ(Δυ) is the line intensity after a Doppler shift caused by the line-of-sight velocity, Δυ, at each χ coordinate along the disc. For a rigid rotating sphere, where there is no latitudinal change in the rotation rate, all elements on the disc surface along the same χ coordinate undergo the same Doppler shift, such that x = Δυ/Δυ_rot, where Δυ_rot is the Doppler shift due to the equatorial velocity (r = r₁). Eq. A.1 is given by the ratio of the length of the chord at constant x, to the length of the chord at x = 0. This expression can be adapted for a uniformly irradiated annulus projected by a rigid rotating spherical shell (see Fig. A.l), such as a planetary atmosphere. Fig A.l shows a rotating annulus with inner radius r₁ = R_p. We set r₁ = R_p = 1, such that the outer radius r₂ = 1 + d, where d denotes the thickness of the atmosphere relative to the planetary radius R_p. Thus, the length of the chord at x = 0 is 2d (see Boucher et al. (2023) for a similar expression with alternative notation, in which the absolute height of the atmosphere z = d· R_p). The length of the chord at constant x varies for different regions of the annulus.

Outer region (r₁ ≤ |x| < r₂)

When (r₁ ≤ |x| < r₂), the outer disc (r = r₂) is not occulted by the planetary disc, meaning the lines simply undergo rotational broadening due to a rigid rotating disc ranging from −r₂ → r₂. Therefore the length of the chord at constant x is simply

$$\begin{array}{l}2y_2=2\sqrt{r_2^2-x_2^2}\hfill \\ =2\sqrt{{(1+d)}^2-x_2^2}\hfill \end{array}$$

As all elements on the disc surface along the same x coordinate undergo the same Doppler shift, this expression becomes

$$2y_2=2\sqrt{{(1+d)}^2-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}$$

Such that Eq. A.1 is

$$\begin{array}{l}G(\text{Δ}v)=\frac{2y_2}{2d}\hfill \\ =\frac{\sqrt{{(1+d)}^2-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}}{d}\hfill \end{array}$$(A.2)

Inner region (|x| < r₁)

However, when |x| < r₁, the outer disc (r = r₂) is occulted by the inner disc (r = r₁), such that the lines undergo rotational broadening due to a rotating spherical shell of inner radius r₁ and outer radius r₂. The length of the chord at constant x is

$$2\left(y_2-y_1\right)=2\left(\sqrt{{(1+d)}^2-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}-\sqrt{1-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}\right)$$

and so

$$\begin{array}{l}G(\text{Δ}v)=\frac{2\left(y_2-y_1\right)}{2d}\hfill \\ =\frac{\sqrt{{(1+d)}^2-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}-\sqrt{1-{\left(\text{Δ}v/\text{Δ}{v}_{\text{rot}}\right)}^2}}{d}\hfill \end{array}$$(A.3)

Fig. A.1 Schematic diagram of a rotating annulus of inner radius r1 and outer radius r2 = r1 + d · r1. A portion of the inner and outer regions described in the text is highlighted. Chords at constant x, of length 2y2 and 2(y2 − y1), are also shown.

Appendix B Stellar parameters

Appendix C Injection and recovery tests

Fig. C.1 Summary of our injection and recovery tests for T1 outlined in Sect. 3.3, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red and grey posterior distributions represent independent sub-chains, of the same MCMC chain, both converging to similar distributions, with the injected values highlighted by horizontal and vertical dashed green lines. The various global, morning, and evening parameters are highlighted. The T-P profiles on the right were computed from 10,000 random samples of the MCMC; the solid curve shows the median profile, and the shaded regions show the 1σ, 2σ, and 3σ contours, for both the morning and evening limbs. The dashed green curve shows the injected T-P profile. This corner diagram was generated using CORNER (Foreman-Mackey 2016).

Fig. C.2 Same as Fig. C.1, but for T2.

Fig. C.3 Same as Fig. C.1, but for T3.

Appendix D Atmospheric retrievals

Fig. D.1 Absolute abundance comparison across each of the three transits, for the morning limb. The stellar values from Tabernero et al. (2021) are shown as a vertical grey line, with their 1σ contours shaded where measured. The solar values from Asplund et al. (2009) are shown as a vertical orange line, with their 1σ contours shaded. The weighted mean of the retrieved abundances in regions C and D (morning) from G2023 are also shown, where measured. Similarly, the retrieved abundances for the whole terminator from P2023 are also shown for Na and VO.

Fig. D.2 Absolute abundance comparison across each of the three transits, for the evening limb. The shaded regions are as in Fig. D.1. The weighted mean of the retrieved abundances in regions A and B (evening) from G2023 are also shown, where measured.

Fig. D.3 Summary of our atmospheric retrieval results for T1 outlined in Sect. 3.4, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red and grey posterior distributions represent independent sub-chains, of the same MCMC chain, both converging to similar distributions. The various global, morning, and evening parameters are highlighted. The T-P profiles on the right were computed from 10,000 random samples of the MCMC; the solid curve shows the median profile, and the shaded regions show the 1σ, 2σ, and 3σ contours, for both the morning and evening limbs. The horizontal dotted lines are the median retrieved cloud deck pressures, log10(Pcloud), below which the model spectrum is truncated. Also plotted are the mean contribution functions for the ESPRESSO wavelength range. This corner diagram was generated using CORNER (Foreman-Mackey 2016).

Fig. D.4 Same as Fig. D.3, but for T2.

Fig. D.5 Same as Fig. D.3, but for T3.

Appendix E Cross-correlation

Fig. E.1 Kp – Δvsys maps of VO for T1, T2, and T3, confirming the novel detection of VO by P2023. The VO atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals and broadened with the respective optimum broadening parameters, using the broadening kernel given by Eq. 1. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values.

Fig. E.2 Kp – Δvsys maps of VO, for T1, T2, and T3, using both the VOMYT (McKemmish et al. 2016) and HyVO (Bowesman et al. 2024) line lists. The VO cross-sections were calculated at T = 2800 K and P = 0.1 mbar. The transmission spectra were computed using an isothermal T-P profile and convolved with a 1D Gaussian kernel. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values.

Fig. E.3 CC trails of various species for T3, shifted to the optimum Kp values, as well as to the Kp value from Ehrenreich et al. (2020), in order to visualise any apparent kink. The atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals but were un-broadened so that any deviation from the planet’s Keplerian velocity curve are visible. For clarity, we only plot species in which there was a clearly visible CC trail.

References

All Tables

All Figures

Fig. 1 2D Na model convolved with the broadening kernel outlined in Sect. 3.2. Upper left: rotational broadening kernel outlined in Eq. (1). The solid lines red and blue curves show the final morning and evening halves, respectively, after convolution with an instrumental profile, shown in grey. Upper right: asymmetric Na model, with varying broadening and atmospheric parameters between the morning and evening limbs. The combined model is overplotted with a dashed black line. Middle: 1D model projected to a 2D orbital phase vs wavelength array for a single ESPRESSO order. This array is then weighted by the limb-darkening model described in Gandhi et al. (2022). This ensures solely morning or evening contributions during ingress or egress, respectively. Lower: 2D orbital phase vs wavelength array after SYSREM model filtering, highlighting the distortions of the planetary lines after preprocessing. The dotted black curve shows the planetary velocity to which the model is projected.

In the text

Fig. 2 Relative abundance comparison across each of the three transits, for the morning limb. The stellar values from Tabernero et al. (2021) are shown as a vertical grey line, with their 1σ contours shaded where measured. The solar values from Asplund et al. (2009) are shown as a vertical orange line, with their 1σ contours shaded. The weighted mean of the retrieved relative abundances in regions C and D (morning) from G2023 are also shown, where measured. Otherwise, the relative abundances measured for the whole terminator from P2023 are shown.

In the text

	Fig. 3 Relative abundance comparison across each of the three transits, for the evening limb. Shaded regions are as in Fig. 2. The weighted mean of the retrieved relative abundances in regions A and B (evening) from G2023 are also shown, where measured.
In the text

Fig. 4 Kp – Δvsys maps of V for T1, T2, and T3. Upper: the V atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals but are un-broadened. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values. Lower: comparison of the CCF of Fe and V, at their optimum Kp value. As the models were un-broadened, this distribution can be used as a proxy for the average line shape of each species in our data, highlighting any apparent asymmetries in the species’ longitudinal distributions.

In the text

	Fig. A.1 Schematic diagram of a rotating annulus of inner radius r1 and outer radius r2 = r1 + d · r1. A portion of the inner and outer regions described in the text is highlighted. Chords at constant x, of length 2y2 and 2(y2 − y1), are also shown.
In the text

Fig. C.1 Summary of our injection and recovery tests for T1 outlined in Sect. 3.3, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red and grey posterior distributions represent independent sub-chains, of the same MCMC chain, both converging to similar distributions, with the injected values highlighted by horizontal and vertical dashed green lines. The various global, morning, and evening parameters are highlighted. The T-P profiles on the right were computed from 10,000 random samples of the MCMC; the solid curve shows the median profile, and the shaded regions show the 1σ, 2σ, and 3σ contours, for both the morning and evening limbs. The dashed green curve shows the injected T-P profile. This corner diagram was generated using CORNER (Foreman-Mackey 2016).

In the text

	Fig. C.2 Same as Fig. C.1, but for T2.
In the text

	Fig. C.3 Same as Fig. C.1, but for T3.
In the text

Fig. D.1 Absolute abundance comparison across each of the three transits, for the morning limb. The stellar values from Tabernero et al. (2021) are shown as a vertical grey line, with their 1σ contours shaded where measured. The solar values from Asplund et al. (2009) are shown as a vertical orange line, with their 1σ contours shaded. The weighted mean of the retrieved abundances in regions C and D (morning) from G2023 are also shown, where measured. Similarly, the retrieved abundances for the whole terminator from P2023 are also shown for Na and VO.

In the text

	Fig. D.2 Absolute abundance comparison across each of the three transits, for the evening limb. The shaded regions are as in Fig. D.1. The weighted mean of the retrieved abundances in regions A and B (evening) from G2023 are also shown, where measured.
In the text

Fig. D.3 Summary of our atmospheric retrieval results for T1 outlined in Sect. 3.4, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red and grey posterior distributions represent independent sub-chains, of the same MCMC chain, both converging to similar distributions. The various global, morning, and evening parameters are highlighted. The T-P profiles on the right were computed from 10,000 random samples of the MCMC; the solid curve shows the median profile, and the shaded regions show the 1σ, 2σ, and 3σ contours, for both the morning and evening limbs. The horizontal dotted lines are the median retrieved cloud deck pressures, log10(Pcloud), below which the model spectrum is truncated. Also plotted are the mean contribution functions for the ESPRESSO wavelength range. This corner diagram was generated using CORNER (Foreman-Mackey 2016).

In the text

	Fig. D.4 Same as Fig. D.3, but for T2.
In the text

	Fig. D.5 Same as Fig. D.3, but for T3.
In the text

Fig. E.1 Kp – Δvsys maps of VO for T1, T2, and T3, confirming the novel detection of VO by P2023. The VO atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals and broadened with the respective optimum broadening parameters, using the broadening kernel given by Eq. 1. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values.

In the text

Fig. E.2 Kp – Δvsys maps of VO, for T1, T2, and T3, using both the VOMYT (McKemmish et al. 2016) and HyVO (Bowesman et al. 2024) line lists. The VO cross-sections were calculated at T = 2800 K and P = 0.1 mbar. The transmission spectra were computed using an isothermal T-P profile and convolved with a 1D Gaussian kernel. The dotted red lines show the expected planetary velocity values, whereas the dashed white lines show the maximum Kp and Δvsys values.

In the text

Fig. E.3 CC trails of various species for T3, shifted to the optimum Kp values, as well as to the Kp value from Ehrenreich et al. (2020), in order to visualise any apparent kink. The atmospheric models were calculated using the optimum parameters from our various atmospheric retrievals but were un-broadened so that any deviation from the planet’s Keplerian velocity curve are visible. For clarity, we only plot species in which there was a clearly visible CC trail.

In the text