Spatially resolving the terminator: variation of Fe, temperature, and winds in WASP-76 b across planetary limbs and orbital phase

ABSTRACT

Exoplanet atmospheres are inherently three-dimensional systems in which thermal/chemical variation and winds can strongly influence spectra. Recently, the ultra-hot Jupiter WASP-76 b has shown evidence for condensation and asymmetric Fe absorption with time. However, it is currently unclear whether these asymmetries are driven by chemical or thermal differences between the two limbs, as precise constraints on variation in these have remained elusive due to the challenges of modelling these dynamics in a Bayesian framework. To address this, we develop a new model, HyDRA-2D, capable of simultaneously retrieving morning and evening terminators with day-night winds. We explore variations in Fe, temperature profile, winds, and opacity deck with limb and orbital phase using VLT/ESPRESSO observations of WASP-76 b. We find Fe is more prominent on the evening for the last quarter of the transit, with |$\log (X_\mathrm{Fe}) = {-4.03}^{+0.28}_{-0.31}$|⁠, but the morning shows a lower abundance with a wider uncertainty, |$\log (X_\mathrm{Fe}) = {-4.59}^{+0.85}_{-1.0}$|⁠, driven by degeneracy with the opacity deck and the stronger evening signal. We constrain 0.1-mbar temperatures ranging from |$2950^{+111}_{-156}$| to |$2615^{+266}_{-275}$| K, with a trend of higher temperatures for the more irradiated atmospheric regions. We also constrain a day-night wind speed of |$9.8^{+1.2}_{-1.1}$| km s⁻¹ for the last quarter, higher than |$5.9^{+1.5}_{-1.1}$| km s⁻¹ for the first, in line with general circulation models. We find our new spatially and phase-resolved treatment is statistically favoured by 4.9σ over traditional 1D-retrievals, and thus demonstrate the power of such modelling for robust constraints with current and future facilities.

radiative transfer, methods: numerical, techniques: spectroscopic, planets and satellites: atmospheres, planets and satellites: composition

1 INTRODUCTION

Ultra-hot Jupiters (UHJs), with temperatures in excess of 2000 K, are ideal laboratories for exoplanet study, given their hot extended atmospheres and large signal-to-noise. The temperature for these exoplanets is high enough that refractory species remain gaseous and therefore detectable in the photosphere. These conditions are therefore ideal for studying study a range of physical phenomena at such extreme temperatures, such as the thermal dissociation, rain-out, and ionization of atomic species. In recent years, a number of UHJs have been characterized, with numerous species detected in their atmosphere in both primary and secondary eclipse using ground based high-resolution spectroscopy (HRS) at optical wavelengths (e.g. Hoeijmakers et al. 2019; Seidel et al. 2019; Bourrier et al. 2020; Cabot et al. 2020; Pino et al. 2020; Merritt et al. 2021). These observations are ideal to detect trace species in the atmosphere due to their high sensitivity and resolution of R ≳ 25 000, allowing for clear absorption and emission lines to be extracted from the spectra (see e.g. review by Birkby 2018).

WASP-76 b has had one of the most well-characterized atmospheres UHJs to date with a mass of 0.894M_J, radius of 1.854R_J, and an equilibrium temperature of ∼2160 K with full redistribution (West et al. 2016). Recently, Ehrenreich et al. (2020) showed variation in the strength and position of the Fe signal in WASP-76 b with primary eclipse observations using the ESPRESSO spectrograph on the VLT, and proposed that this variation could be a result of condensation of Fe on the night side of the atmosphere. These observations have also been used to detect numerous atomic and ionic species in the atmosphere (Tabernero et al. 2021; Kesseli et al. 2022) and to resolve more of the asymmetric condensation features for some species (Kesseli et al. 2022). The asymmetry and relative shift of ∼10 km s⁻¹ in the position of the Fe signal during the transit was additionally confirmed by observations with HARPS (Kesseli & Snellen 2021), which were also previously used to detect and study Na in the atmosphere (Seidel et al. 2019; Zák et al. 2019). Both of these observations have also been combined to constrain day-night and vertical winds in the atmosphere by studying the Na doublet (Seidel et al. 2021).

A range of three-dimensional (3D) general circulation models (GCMs) have been used to interpret the asymmetric Fe signal in the atmosphere of WASP-76 b. Such complex models simulate a wide range of physical and chemical processes and provide an unparalleled levels of spatial and dynamical detail for hot Jupiters (e.g. Showman et al. 2009; Dobbs-Dixon & Agol 2013; Rauscher & Menou 2013; Mayne et al. 2014), as well as UHJs (e.g. Tan & Komacek 2019). Importantly, such 3D models can also be used to predict spatial variability in the abundance of the chemical species and temperature structure of the atmosphere and hence in each limb of the terminator, as well as velocity shifts of the spectra due to winds. These models have shown that the Fe asymmetry in WASP-76 b could be a result of cloud formation on the cooler morning side (Savel et al. 2022), or the thermal asymmetry and wind profiles within the atmosphere (Wardenier et al. 2021). Hence, it is currently unclear whether Fe is inherently weaker on the leading limb due to its lower abundance or whether clouds/temperature can affect the signal significantly to alter its strength.

In contrast to GCMs, retrieval frameworks of exoplanetary atmospheres determine the best-fitting parameters for exoplanet atmospheres from observations, often exploring many millions of models and encompassing a wide range in temperature, clouds, and abundances of relevant species. Such models generally assume spatially homogeneous atmospheric conditions throughout the transit. Retrievals have been ubiquitous for interpretation of low-resolution observations in the last decade, most notably for HST spectra in both emission and transmission geometries (e.g. Madhusudhan & Seager 2009; Madhusudhan et al. 2014; Line et al. 2016; Barstow et al. 2017; Benneke et al. 2019; Mikal-Evans et al. 2019; Zhang et al. 2020). However, their use in analysing high-resolution spectra has been relatively recent (e.g. Brogi et al. 2017; Gandhi et al. 2019; Line et al. 2021; Pelletier et al. 2021; Gibson et al. 2022), and has been made possible thanks to the development of cross-correlation to likelihood mappings (Brogi & Line 2019; Gibson et al. 2020), which allow for Bayesian analyses to be performed on the data. However, high-resolution retrievals remain a computational challenge due to the many millions of models required at high resolutions of R ≳ 25 000 and the various analysis techniques used to extract the planetary signal from the noise.

Winds in exoplanetary atmospheres have also recently been detected by a range of ground based high-resolution observations, and result in an observable Doppler shift and/or broadening of the overall spectrum arising from the range in velocities along the line of sight. Snellen et al. (2010) detected the winds and orbital motion of an exoplanet for the first time through high-resolution ground based observations in the infrared. Louden & Wheatley (2015) have constrained the eastward winds in the hot Jupiter HD189733 b through forward modelling of the Na absorption lines, showing that including such dynamics is key for accurate characterization. Furthermore, GCMs have also been shown to provide better fits to high-resolution observations than 1D models that do not incorporate thermal variation and dynamical processes (Beltz et al. 2021). Including winds into atmospheric models is also of particular importance, given that there are also a growing number of high-resolution studies that have constrained these on the terminator (Brogi et al. 2016; Flowers et al. 2019; Seidel et al. 2021). In addition, recent work has also shown that cloud formation on the night side of UHJs is strongly dependent on the atmospheric dynamics (Komacek et al. 2022) and these clouds can affect phase-dependent emission spectra.

Ultimately, what is needed is a model that is capable of interpreting the variability in the observations of a UHJ, such as WASP-76 b, whilst simultaneously incorporating dynamical processes in a statistical framework. We must include as much of the information content of a GCM whilst keeping the model computationally tractable with the minimum number of additional parameters to encompass the relevant physics. Such a retrieval model must therefore be capable of constraining separate Fe abundances, thermal profiles and opacity decks for the morning (leading) and evening (trailing) limbs, and include winds. The winds will act to not only transfer thermal energy and material between the day and night sides of the atmosphere, but also shift the spectrum for each limb according to the velocity profile.

Here we develop HyDRA-2D, a spatially resolved retrieval model that incorporates separate constraints on morning and evening limbs and a wind travelling between the day and night sides of the atmosphere. Such a model represents the first step into multi-dimensional retrievals on high-resolution data. HyDRA-2D allows for a full exploration of the parameter space using Bayesian algorithms whilst capturing the dynamical effects that can strongly influence the observed spectra of exoplanets. We separate the atmosphere into a morning and evening limb, as these regions of the terminator exhibit different velocities due to the planet’s rotation, as well as different thermochemical states. This rotation therefore results in separate spectral contributions for each side.

We use HyDRA-2D to retrieve the two nights of ESPRESSO observations of WASP-76 b (Ehrenreich et al. 2020) and determine the variability in Fe, temperature, opacity, and wind speed across each limb of the planet. Whilst there are additional transit observations of WASP-76 b observed with HARPS (Seidel et al. 2019; Kesseli & Snellen 2021), we restrict ourselves to the VLT observations due to the higher signal-to-noise ratio and spectral resolution available with ESPRESSO, which is vital for discerning the different regions of the planet that are separated by small velocity shifts. We perform two separate retrievals for the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| phase ranges of the observations, with ϕ = 0 corresponding to the mid-transit time. The chosen phase ranges allow us to also explore the variation of atmospheric properties with phase, given that the region of the planetary atmosphere probed varies by ∼30° between ingress and egress (Wardenier, Parmentier & Lee 2022). We exclude the phases between ϕ = −0.02 and +0.02 from our analysis owing to the Doppler shadow of the star crossing the planetary signal during mid-transit. In addition to retrievals with HyDRA-2D, we also perform a set of spatially homogeneous 1D retrievals across the combined phase range as well as for each of the ranges separately to compare the differences in the constrained parameters with model sophistication.

In the next section, we discuss our HyDRA-2D retrieval setup, followed by the results and discussion section where we discuss our constraints for each side of the atmosphere and comparisons with GCMs, and finally end with the concluding remarks.

2 METHODS

In this section, we discuss the setup for HyDRA-2D. We retrieve two nights of observations of the primary eclipse of WASP-76 b with the ESPRESSO spectrograph on the VLT (Ehrenreich et al. 2020). The observations span the optical, where Fe has prominent opacity. As these observations show strong variability in both the strength and position of the Fe signal, we must allow our retrieval to be flexible enough to incorporate such dynamics and variability in a computationally efficient manner. This need for efficiency leads us to a new parametrized approach where 3D effects are encoded as extra parameters in the description of the physical properties of the atmosphere (e.g. temperature profiles, abundances and opacity decks) and the dynamics of the atmosphere (e.g. velocity shifts and broadening of spectral lines). We first discuss the line-by-line opacity calculations for Fe used to compute the high-resolution model spectra. We then discuss how we separate the terminator into two distinct regions and the retrieval parameters for each side. Following this, we discuss how we incorporate the day-night wind and limb darkening, and then describe our analysis of the data and the model sequence.

2.1 Fe opacity

We include Fe opacity through computation of the cross-section on a grid of the relevant pressures and temperatures using the Kurucz line list (Kurucz & Bell 1995). This line list provides the log (gf) value for each transition, where g is the degeneracy and f refers to the oscillator strength. For our work, we first compute the line strength in cm⁻¹/(atom cm⁻²) for each temperature T on our grid, through

$$\begin{eqnarray} S(T) = \frac{gf {\pi} e^2}{m_{\rm e}c^2} \frac{\mathrm{exp}(-c_2E_\mathrm{lower}/T)}{Q(T)} (1-\mathrm{exp}(-c_2\nu _0/T)), \end{eqnarray}$$

(1)

where ν₀ = E_upper − E_lower, with the upper and lower state energies of the transition E_upper and E_lower given in cm⁻¹. The charge of an electron, e, the speed of light c, and the mass of the electron m_e are all given in cgs units, and c₂ = 1.438 7769 cm K. In addition, the partition function Q is defined as

$$\begin{eqnarray} Q(T) = \sum _j g_j \mathrm{exp}(-c_2E_j/T) \end{eqnarray}$$

(2)

for a state j with degeneracy g_j.

We apply line broadening for each of the transitions, which arises from the radiative damping constant (natural broadening) and van der Waals broadening. These result in a Lorentzian line profile at a given frequency ν (in cm⁻¹),

$$\begin{eqnarray} f_{\rm L}(\nu -\nu _0) = \frac{1}{{\pi} }\frac{\gamma _{\rm L}}{(\nu -\nu _0)^2+\gamma _{\rm L}^2}. \end{eqnarray}$$

(3)

Here, γ_L is defined as (Sharp & Burrows 2007)

$$\begin{eqnarray} \gamma _{\rm L} = \frac{1}{4{\pi} c}\bigg (\gamma _{\rm W} N \bigg (\frac{T}{10\,000}\bigg)^{0.3} + \gamma _{\rm N}\bigg), \end{eqnarray}$$

(4)

where γ_W and γ_N are the van der Waals and natural broadening coefficients, respectively, and N is the number density of gas, which, in our case, consists of atomic hydrogen. For our purposes, we ignore the Stark broadening of the line, given that the atmosphere is not expected to be highly ionized.

In addition to Lorentzian broadening, thermal broadening arising from the temperature of the gas will also contribute to the overall profile. This contribution results in a Gaussian line profile (e.g. Hill, Yurchenko & Tennyson 2013),

$$\begin{eqnarray} f_{\rm G}(\nu -\nu _0) &=& \frac{1}{\gamma _{\rm G}\, \root \of {{\pi} }}\,\mathrm{exp}\left(-\frac{(\nu -\nu _0)^2}{\gamma _{\rm G}^2}\right),\\\gamma _{\rm G} &=& \root \of {\frac{2k_{\rm b}T}{m}}\frac{\nu _0}{c}, \end{eqnarray}$$

(5),(6)

where k_b is the Boltzmann constant and m is the mass of an Fe atom. The overall line shape is thus a convolution between the Lorentzian and Gaussian profile, given by the Voigt profile (e.g. Gandhi & Madhusudhan 2017),

$$\begin{eqnarray} f_V(\nu -\nu _0,\gamma _{\rm L},\gamma _{\rm G}) = \int _{-\infty }^{\infty }f_{\rm G}(\nu ^{\prime }-\nu _0)f_{\rm L}(\nu -\nu ^{\prime }){\rm d}\nu ^{\prime }. \end{eqnarray}$$

(7)

Thus, the cross-section as a function of frequency for each transition line is then computed as

$$\begin{eqnarray} \sigma _\mathrm{Fe} = S(T)f_V(\nu - \nu _0,\gamma _{\rm L},\gamma _{\rm G}). \end{eqnarray}$$

(8)

The sum of the contributions from each line in the line list gives the total atomic cross-section of Fe. For each pressure and temperature on our grid, we compute the cross-section at 0.01 cm⁻¹ spacing, corresponding to a resolution of R = 2.5 × 10⁶ at 0.4 |$\mu$|m. Further details on line broadening and computation of the cross-section can be found in Gandhi et al. (2020b).

We compute the overall opacity arising from Fe through the volume mixing ratio, X_Fe, which we leave as a free parameter as listed in Table 1. The extinction coefficient χ (in m⁻¹) is given by

$$\begin{eqnarray} \chi (P,T,\lambda) = X_\mathrm{Fe} N \sigma _\mathrm{Fe}(P,T,\lambda), \end{eqnarray}$$

(9)

where σ_Fe is the total cross-section of Fe and N is the total number density of gas. The contribution to the optical depth dτ arising from an element of the atmosphere of length dz is then

$$\begin{eqnarray} \mathrm{d}\tau = \chi (P,T, \lambda) \mathrm{d}z. \end{eqnarray}$$

(10)

Note that we include two free parameters for the Fe abundance, one for the morning and one for the evening side, as discussed below.

Table 1.

Parameters and uniform prior ranges for HyDRA-2D retrieval.

Parameter	Prior range
	Evening side	Morning side
log (XFe)	−15 → −2	−15 → −2
T0.1 mbar / K	1500 → 3800	1500 → 3800
|$\alpha _1\, /\, \mathrm{K}^{-{1}/{2}}$|	0 → 1	0 → 1
|$\alpha _2\, /\, \mathrm{K}^{-{1}/{2}}$|	0 → 1	0 → 1
log (P1/bar)	−7 → 2	−7 → 2
log (P2/bar)	−7 → 2	−7 → 2
log (P3/bar)	−2 → 2	−2 → 2
log (αhaze)	−4 → 6	−4 → 6
γhaze	−20 → −1	−20 → −1
log (Pcl/bar)	−7 → 2	−7 → 2
ΔKp / kms−1	−20 → 20
ΔVsys / km s−1	−20 → 20
δVwind / km s−1	1 → 20

Notes. We retrieve separate Fe abundances, temperature profile, and opacity/haze decks for each of the morning and evening sides of the terminator. The ΔK_p, ΔV_sys, and δ_wind parameters are the same for both sides (see Fig. 1).

2.2 Separating the morning and evening limbs

For our work, we assume separate physical properties for the leading (or morning side) and trailing limbs (or evening side) of the terminator (see Fig. 1). This is because we are able to separate the spectrally contributing regions of the terminator due to their difference in velocity space, given the high resolution (R = 140 000) and signal-to-noise ratio. We therefore retrieve separate Fe volume mixing ratios for the morning and evening sides of the atmosphere. In addition, we also retrieve separate pressure-temperature profiles and opacity decks/hazes for each half of the terminator.

Figure 1.

Schematic of the transit of WASP-76 b and the retrieval parameters for HyDRA-2D. For each of the evening and morning limbs, we retrieve the gaseous iron abundance with X_Fe, the temperature profile is set through the six parameters T_mbar, α₁, α₂, P₁, P₂ and P₃ using the parametrization in Madhusudhan & Seager (2009), and the opacity and haze are parametrized with α_haze, γ_haze, and P_cl. We also retrieve ΔK_p, ΔV_sys, and δV_wind to constrain the wind profile of the atmosphere (see Section 2.3).

2.2.1 Temperature profile

We parametrize the temperature profile for each side using the method described in Madhusudhan & Seager (2009). This requires breaking the atmosphere into 3 distinct regions, with the top two regions capable of having a varying temperature profile with pressure, and the final deepest layer assuming an isothermal temperature expected to be common at high pressures on hot and UHJs (e.g. Fortney et al. 2008). We thus obtain the following parametrization for T as a function of pressure P,

$$\begin{eqnarray} T(P) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}T_0 + \Big (\frac{\ln {P/P_0}}{\alpha _1}\Big)^2 & \rm{if P\le P_1},\\ T_2 + \Big (\frac{\ln {P/P_2}}{\alpha _2}\Big)^2 & \rm{if P_1 \lt P \le P_3},\\ T_3 & \rm{if P\gt P_3}, \end{array}\right. \end{eqnarray}$$

(11)

and where we enforce P₁ ≤ P₃. We additionally fix P₀ = 10⁻⁷ bar as the top of the model atmosphere. At lower pressures, the assumptions of hydrostatic equilibrium and local thermodynamic equilibrium begin to break down. We allow the value of P₂ to vary between P₀ ≤ P₂ ≤ P₃. Importantly, such a constraint allows for a thermal inversion in layer 2 in addition to isothermal and non-inverted temperature profiles. This is done to allow the retrieval to determine the best-fitting P − T profile without the requirement that the temperature be monotonically decreasing with altitude. By enforcing continuity of the temperature between the three layers, we can reduce this set of equations to just six free parameters, a temperature value at any given pressure, α₁, α₂, P₁, P₂, and P₃. For our retrieval, we retrieve the temperature at 0.1 mbar, T_0.1 mbar, which is representative of a typical photospheric pressure in transmission geometry for such high-resolution observations. We have also verified our results by using the temperature at the top of the model atmosphere as our free parameter, as well as with the 0.1-bar temperature as a free parameter, to ensure this did not impact our constraints. The prior ranges for each parameter is shown in Table 1. We use different sets of parameters for the evening side and morning sides, which ensures that the temperature constraints for one side are independently constrained of the other, and hence overall we have 12 free parameters in the retrieval. Further details on the parametrization of the P − T profile can be found in Madhusudhan & Seager (2009).

2.2.2 Grey opacity deck and haze

In addition to the opacity resulting from the presence of gaseous Fe, other sources of opacity may play a key role in determining the spectrum. These opacities arise from of a number of factors, such as other molecular and atomic species in the atmosphere, and non-Rayleigh scattering from small particles/hazes in the atmosphere. These result in two distinct sources of extinction in the atmosphere, namely that of a wavelength-independent opacity deck and a wavelength-dependent super-Rayleigh haze. For our prescription, we adopt a one sector approach with homogeneous opacity and haze parameters (e.g. MacDonald & Madhusudhan 2017; Welbanks & Madhusudhan 2021) for each of the morning and evening limbs. We add this source of opacity into the overall extinction coefficient of the atmosphere, χ, as

$$\begin{eqnarray} \chi (P,\lambda) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\alpha _\mathrm{haze} N\sigma _0\big (\frac{\lambda }{\lambda _0}\big)^{\gamma _\mathrm{haze}} & \rm{if P\le P_\mathrm{cl}},\\ \infty & \rm{otherwise}, \end{array}\right. \end{eqnarray}$$

(12)

where N is the total number density of gas, σ₀ = 5.31 × 10⁻³¹ m² and λ₀ = 0.35 |$\mu$|m. Hence, the free parameters are α_haze, γ_haze, and P_cl, whose prior ranges are given in Table 1. This therefore results in six free parameters in our retrieval, given that we assume the two sides of the atmosphere are independent. The overall opacity of each side of the atmosphere is thus the sum of the contributions arising from Fe line absorption and the opacity deck/haze. Further details on the opacity deck/haze prescription can be found in Welbanks & Madhusudhan (2021).

2.2.3 Computing the unbroadened spectrum

We compute the unbroadened spectrum using numerical radiative transfer for each side of the atmosphere (e.g. Welbanks & Madhusudhan 2021). This involves computation of the transit depth Δ_e and Δ_m for the evening and morning sides, and is given by

$$\begin{eqnarray} \Delta _{\rm e}(\lambda) &=& \frac{1}{2}\, \frac{1}{R_\mathrm{s}^2}\Bigg (R_\mathrm{p}^2 + 2 \int _{R_\mathrm{p}}^{\infty }b_{\rm e} (1-{\rm e}^{-\tau _{\rm e}(\lambda)}) \, \mathrm{d}b_{\rm e} \\ &&-\, 2\int _{0}^{R_\mathrm{p}}b_{\rm e} {\rm e}^{-\tau _{\rm e}(\lambda)} \, \mathrm{d}b_{\rm e}\Bigg), \end{eqnarray}$$

(13)

$$\begin{eqnarray} \Delta _{\rm m}(\lambda) &=& \frac{1}{2}\, \frac{1}{R_\mathrm{s}^2}\Bigg (R_\mathrm{p}^2 + 2 \int _{R_\mathrm{p}}^{\infty }b_{\rm m} (1-{\rm e}^{-\tau _{\rm m}(\lambda)}) \, \mathrm{d}b_{\rm m} \\&&-\, 2\int _{0}^{R_\mathrm{p}}b_{\rm m} {\rm e}^{-\tau _{\rm m}(\lambda)} \, \mathrm{d}b_{\rm m}\Bigg), \end{eqnarray}$$

(14)

where R_p and R_s refer to the radius of the planet and star respectively, τ refers to the optical depth, and the impact parameter is given by b. The subscripts e and m refer to the evening and morning sides, and the factor of 1/2 accounts for the fact that each side is only integrated for half of the limb. The physical properties of each side such as temperature and Fe abundance are variable, and thus the overall optical depth for each side is different. We set our reference radius |$R_\mathrm{p} = 1.854 \, R_\mathrm{J}$| at 10⁻² bar reference pressure, but the chosen pressure level has no effect on the retrieval given the high-pass filtering in the data analysis removes any offset between different chosen pressures (see Section 2.5). The spectrum for each side is then broadened and combined as discussed below.

2.3 Incorporating winds

Once we have computed the transit spectra for each side of the atmosphere, we must broaden them according to the rotation of the planet and the winds present in the atmosphere. The procedure for this broadening is shown in Fig. 2. This setup was motivated through recent work on general circulation models by Wardenier et al. (2021), which showed that the broadened primary eclipse spectrum is significantly altered and therefore requires careful consideration of both aspects for accurate spectral modelling.

Figure 2.

Schematic of the broadening prescription used for HyDRA-2D. We first compute the unbroadened evening and morning transit depths using the parameters for each side (see Table 1). We then convolve this spectrum with the rotational broadening kernel for each side. The evening side is rotating towards the observer and the morning side is rotating away, and thus are blueshifted and redshifted, respectively. We then shift and broaden both sides by a day-night wind with ΔV_sys, and δV_wind, which are free parameters in the retrieval. We then combine the spectra for each limb after incorporating the stellar limb darkening and model reprocessing, as discussed in Sections 2.4 and 2.5.

2.3.1 Planetary rotation

We first begin by applying the planetary rotation to each side. The evening side, or trailing limb, is rotating towards the observer and the morning side, or leading limb, is rotating away. Therefore, the evening and morning sides will be blueshifted and redshifted, respectively, according to the rotation of the planet. However, the polar regions of the terminator will be rotating more slowly than the equator and hence there is also a contribution arising at lower velocities. For a rotation velocity, V_rot, the full rotation kernel is given by

$$\begin{eqnarray} R(x) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{{\pi} \sqrt{1 - x^2/V_\mathrm{rot}^2}} & \rm{if |x|\lt V_\mathrm{rot}}\\0 & \rm{otherwise}. \end{array}\right. \end{eqnarray}$$

(15)

We only apply half of the kernel for each half of the terminator, as shown in Fig. 2. This convolves our unbroadened spectrum according to the respective kernel for each side, which also results in a net shift for each side. We keep the rigid-body rotation kernel fixed for each model in our retrieval as the rotational velocity of WASP-76 b is well known, given that the planet is expected to be tidally locked (e.g. Showman & Guillot 2002).

An example spectrum highlighting the effect of the broadening is shown in Fig. 3 for a small region near ∼0.4 |$\mu$|m. The rotation results in a split of the overall line profiles in the spectrum, with two distinguishable peaks in the rotationally broadened spectrum arising from contributions from each half of the terminator. In addition, the spectrum is more spread in wavelength and therefore the sharp peak arising from Fe line opacity in the unbroadened spectrum has reduced in extent. Such rotation is therefore vital to account for as it significantly alters the line profile, and may otherwise lead to biased estimates of abundance.

Figure 3.

Example spectra in the ∼0.4-|$\mu$|m range showing the effect of rotational and wind broadening on our model. The black line shows the unbroadened spectrum, which we set to be the same for both the evening and morning sides for this illustration. We apply the rotational broadening as shown in purple, followed by a day-night wind with ΔV_sys = −10 and δV_wind = 8 km s⁻¹ as shown in red. We also show an unbroadened spectrum shifted by −10 km s⁻¹ in blue.

2.3.2 Day-night wind

In addition to the rotation, a day-night wind also contributes to the broadening of the overall spectrum. General circulation models have shown that the highest velocity winds are at the lowest pressures, with deeper regions having almost no wind present (e.g. Miller-Ricci Kempton & Rauscher 2012; Wardenier et al. 2021). Therefore incident rays of light from the host star that pass through the atmosphere first encounter a strong wind at the top of the atmosphere, then pass through to higher pressures with a lower wind speed, and finally exit out of the atmosphere through the lower pressures with higher wind speeds again. This path thus results in a net blueshift as well as a spread in the spectrum due to the range in velocities encountered by the rays of light passing through the atmosphere. In addition, there will also be a spread caused by the wind speed even if it is constant with depth, given that the projected velocity into the observer’s reference frame varies as the path travels through the atmosphere. There may also be an additional contribution arising from variation of the winds with latitude.

To account for these physical mechanisms and processes, we parametrize the wind with a velocity shift from the known systemic velocity, ΔV_sys, and a spread in its value given by the full width at half-maximum (FWHM) δV_wind (see Fig. 2). The ΔV_sys parameter is often included in high-resolution analyses (e.g. Birkby et al. 2013; Pino et al. 2020; Giacobbe et al. 2021), given that the systemic velocity often has some measurement uncertainty that we may be sensitive to. However, in our case, the V_sys parameter is very precisely known at −1.17 ± 0.02 km s⁻¹ (Ehrenreich et al. 2020) and thus we can assume that any deviation ΔV_sys can be attributed to the winds in the atmosphere. We also include one additional parameter, the deviation from the planet’s known orbital velocity, ΔK_p, as this often has a degeneracy with the ΔV_sys value. We therefore have three additional parameters in HyDRA-2D to account for the wind and its uncertainty, ΔK_p, ΔV_sys, and δV_wind as shown in Table 1.

Fig. 3 shows the effect of applying the wind broadening. We now see a net blueshift in the spectrum as a result of the wind ΔV_sys, and an additional spread in the line profile as a result of δV_wind. Note that the spectrum is significantly more spread than the unbroadened spectrum shifted by 10 km s⁻¹, and therefore the wind and rotation have a strong effect on our retrieved results. The resulting spectral-line profile for this broadened model with rotation and wind is similar to the type of line distortions and shifts predicted by GCMs (e.g. Showman et al. 2013; Rauscher & Kempton 2014). We are therefore confident that our parametrization can address the diversity of potential line profiles.

2.4 Stellar limb darkening

For each spectrally broadened half of the terminator, we must also apply limb darkening, given that the brightness of WASP-76’s stellar surface varies significantly in the optical. As the retrievals are over the first and last quarter of the transit, the brightness and therefore the relative contributions of the morning and evening side will thus vary with the orbital phase.

Fig. 4 shows the setup to compute the limb darkening for each side of the atmosphere. We first compute the value of x for a given phase of the observations ϕ, normalized by the stellar radius R_s, given by

$$\begin{eqnarray} x = \frac{a \sin (2{\pi} \phi)}{R_\mathrm{s}}, \end{eqnarray}$$

(16)

where a refers to the semi-major axis of the orbit. We must then add and subtract R_p/R_s for the morning and evening sides, respectively. For simplicity we do not integrate the contribution of the limb darkening over the half-annulus of each side, but take its value at the red and blue points in Fig. 4 for the morning and evening. This is justified, given that there is no significant variability in the limb over the half-annuli and given that the modal contribution to the spectrum comes from these points. For each side, we thus obtain

$$\begin{eqnarray} x_\mathrm{e} = \frac{a\sin (2{\pi} \phi)}{R_\mathrm{s}} - \frac{R_\mathrm{p}}{R_\mathrm{s}}, \end{eqnarray}$$

(17)

$$\begin{eqnarray} x_\mathrm{m} = \frac{a\sin (2{\pi} \phi)}{R_\mathrm{s}} + \frac{R_\mathrm{p}}{R_\mathrm{s}}, \end{eqnarray}$$

(18)

with the morning side given by x_m and the evening side given by x_e. We then compute the radial distance μ of the points from the centre of the star,

$$\begin{eqnarray} \mu _\mathrm{e} = \sqrt{1 - x_\mathrm{e}^2 - b^2}, \end{eqnarray}$$

(19)

$$\begin{eqnarray} \mu _\mathrm{m} = \sqrt{1 - x_\mathrm{m}^2 - b^2}. \end{eqnarray}$$

(20)

The impact parameter b is assumed to be 0.027 (Ehrenreich et al. 2020). We then compute the overall limb-darkening function, given the quadratic limb-darkening coefficients u₁ and u₂:

$$\begin{eqnarray} I_\mathrm{e}(\phi) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 - u_1(1-\mu _\mathrm{e}) - u_2(1-\mu _\mathrm{e})^2 & \rm{if x_\mathrm{e}^2 \le 1 - b^2},\\ 0 & \rm{elsewhere}. \end{array}\right. \end{eqnarray}$$

(21)

Similarly for the morning side, we obtain

$$\begin{eqnarray} I_\mathrm{m}(\phi) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 - u_1(1-\mu _\mathrm{m}) - u_2(1-\mu _\mathrm{m})^2 & \rm{if x_\mathrm{m}^2 \le 1 - b^2},\\ 0 & \rm{elsewhere}. \end{array}\right. \end{eqnarray}$$

(22)

For our work, we have taken u₁ = 0.393 and u₂ = 0.219 (Ehrenreich et al. 2020).

Figure 4.

Schematic showing the setup for the limb-darkening calculations. The outer circle shows the star and the inner circle showing the planet with radius R_p, with its path during the transit given by the dotted line. The impact parameter of the transit is given by b. For the evening and morning side limb-darkening factors, we calculate the brightness at the blue and red points, respectively (see Section 2.4).

We obtain the overall spectrum by simply scaling the model for each side by the value of I and summing (see equation 25). Hence, with this limb-darkening model, we are able to account for the variation in the contributions for each side across the transit; the evening side is the dominant contributor during the latter half of the transit and vice versa for the morning in the first half. In addition, the model ensures that at ingress and egress we obtain a contribution from just one side in the overall spectrum.

2.5 Data analysis and model reprocessing

As the data are the same that were used in Kesseli et al. (2022), the data are cleaned and processed in a similar way. We summarize the steps here (see Fig. 5), but for more details on the observing strategy and initial data reduction, see Ehrenreich et al. (2020), and for more details on the spectral cleaning steps, see Kesseli et al. (2022).

Figure 5.

Example showing the data analysis steps that were performed. The top panel shows a small slice in wavelength of the pipeline-reduced ESPRESSO spectra taken during the transit on 2018 October 30. Each of the 69 rows is a different spectrum that has been shifted to the star’s rest frame. Vertical dark lines are stellar absorption lines. The colour bar shows the pixel values of the spectra. In the middle panel, the spectra have been corrected for changes in signal-to-noise ratio and any difference in the blaze function that occur throughout the night, which removes the horizontal banding seen in the top panel. The bottom panel shows the residuals that are left after the stellar spectrum was removed by dividing by the average out-of-transit spectrum. The exoplanet’s signal is hidden in the noise, and is extracted using cross-correlation.

We downloaded the pipeline reduced ESPRESSO spectra from the Data and Analysis Center for Exoplanets (DACE) data base.¹ Before any other steps, we correct the spectra for telluric absorption lines due to H₂O and O₂ using Molecfit (Smette et al. 2015). Molecfit performs especially well for shallow absorption lines, but is known to return poor results for the line cores of the strongest absorption lines (e.g. Allart et al. 2017). For this reason, we mask all pixels where the transmission through the atmosphere is less than 40 per cent. We then correct for velocity shifts due to the reflex motion of the star (K_*), the barycentric velocity (v_bar), and the system velocity (v_sys). Next, we perform a 5σ clipping to remove any spurious flux spikes due to cosmic rays (top panel of Fig. 5). To remove any broad-band noise present in the spectra while preserving the overall shape, we follow Merritt et al. (2020) and place the spectra on a ‘common’ blaze using a high-pass filter with a width of 200 pixels (∼100 km s⁻¹; middle panel of Fig. 5). Finally, we interpolate each spectrum on to the same wavelength grid and divide each spectrum by the average out-of-transit stellar spectrum to remove the contribution from host star (bottom panel of Fig. 5).

To process our model, we must reproduce the model data cube, beginning first by Doppler shifting the spectrum at each phase by its velocity, V_p. At a given phase ϕ(t), the shift is given by

$$\begin{eqnarray} V_\mathrm{p}(\phi) = \Delta V_\mathrm{sys} + (K_\mathrm{p} + \Delta K_\mathrm{p})\sin (2{\pi} \phi (t)), \end{eqnarray}$$

(23)

where K_p = 196 km s⁻¹ is the planet’s rotation velocity. The ΔV_sys parameter is the day-night wind velocity and represents an overall shift in the planet’s spectrum. We also include ΔK_p to account for uncertainties in the measured value of K_p and because there is often a degeneracy between the K_p and V_sys parameters as a result of the limited phase range that we are exploring during the transit. The new wavelengths of the shifted spectrum are

$$\begin{eqnarray} \lambda (\phi) = \lambda _0 \sqrt{ \frac{1+V_\mathrm{p}(\phi)/c}{1-V_\mathrm{p}(\phi)/c}}, \end{eqnarray}$$

(24)

where λ₀ is the wavelength of the unshifted spectrum. Therefore, the final combined spectrum for the terminator is given by

$$\begin{eqnarray} \Delta _\mathrm{shift}(\phi) = \Delta _{\rm e}^{^{\prime }} (\lambda (\phi)) I_{\rm e}(\phi) + \Delta _{\rm m}^{^{\prime }} (\lambda (\phi)) I_{\rm m}(\phi), \end{eqnarray}$$

(25)

where I_e and I_m refer to the stellar limb-darkening factors, and |$\Delta _{\rm e}^{^{\prime }}$| and |$\Delta _{\rm m}^{^{\prime }}$| refer to the evening and morning side transit depths that have been rotationally and wind broadened and shifted to the new wavelengths λ(ϕ) for each phase. One important final step we take with our model is to apply a high-pass filter, which removes any slowly varying trends in the model and reproduces the data sequence. The high-pass filter also acts to remove the continuum from the spectra, which means that we are only sensitive to the spectral-line features that are generated above the continuum.

Once we have a processed spectrum for the atmosphere, we must compare this to the observations of WASP-76 b and determine the likelihood. We use the cross-correlation to log-likelihood (CC-to-logL) map of Brogi & Line (2019) to compare our model to the data. We do this by computing

$$\begin{eqnarray} \log L = - \frac{N}{2} \log (s_f^2 + s_g^2 - 2R), \end{eqnarray}$$

(26)

where

$$\begin{eqnarray} s_f^2 = \frac{1}{N}\sum _{n=0}^N f^2(n), \end{eqnarray}$$

(27)

$$\begin{eqnarray} s_g^2 = \frac{1}{N}\sum _{n=0}^N g^2(n-s), \end{eqnarray}$$

(28)

$$\begin{eqnarray} R = \frac{1}{N} \sum _{n=0}^N f(n)g(n-s). \end{eqnarray}$$

(29)

Here, f and g are the data and model, respectively, for a model with N spectral points, and s represents the wavelength offset. This form of the likelihood accounts for the overall line shape as well as position; models that poorly match the observations or are scaled incorrectly will be penalized by the R term and thus less favoured. We have two nights of observations for WASP-76 b, and so we sum the likelihood for each night to obtain the overall value. The method in Brogi & Line (2019) also prescribes reprocessing each tested model, i.e. reproducing the effects of the data analysis on the model spectrum. This reprocessing is necessary in the infrared because algorithms used to remove telluric lines are known to alter the exoplanet spectrum as well. However, our analysis in the optical is based on the in-transit/out-of-transit approach and does not distort the planetary signal, and therefore the depth and shape of planet’s spectral lines are preserved relative to the stellar spectrum. We can therefore compare each model to the observations directly, without significant reprocessing of the model. Further details on the likelihood method and derivation can be found in Brogi & Line (2019).

2.6 Retrieval setup

We retrieve the high-resolution observations of WASP-76 b, with two transits obtained on 2018 September 2 and October 30 with the ESPRESSO spectrograph on the VLT (Ehrenreich et al. 2020). These observations are some of the highest signal-to-noise atmospheric observations with high-resolution spectra, obtained at a spectral resolution of R = 140 000 in the optical. We use the 0.4–0.8 |$\mu$|m wavelength range for our retrievals. We set up HyDRA-2D with 23 free parameters, which are given in Table 1. This includes 10 parameters for each of the morning and evening sides, representing the Fe abundance, temperature profile (6 parameters), and opacity and haze deck (3 parameters), and 3 additional parameters that are consistent for both sides, ΔK_p, ΔV_sys, and δV_wind. Each spectral model is generated at a spectral resolution of R = 500 000, assuming a planet radius of |$1.854 \, R_\mathrm{J}$|⁠, planet mass of |$0.894M_\mathrm{J}$|⁠, and stellar radius |$1.756\,{\rm R}_\mathrm{\odot }$|⁠. We use the Nested Sampling algorithm MultiNest (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Buchner et al. 2014) for parameter estimation with HyDRA-2D, adopting 1000 live points.

We perform two separate retrievals, in the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| phase ranges, leaving out the phases in between due to the Doppler shadow of the star. The transit of the planet occurs within −0.0432 ≤ ϕ ≤ 0.0432 and the whole planet is in the transit within the phases of −0.035 ≤ ϕ ≤ 0.035. We choose to use phases up to ±0.04 for our retrieval as at least one side of the planet is in the transit within this range. We have performed tests varying the upper phase range but we did not see any significant difference in the retrieved parameters given that the signal is weaker near ingress and egress due to the dimmer stellar limbs. The two separate retrievals over the phase ranges probe differing regions of the atmosphere due to the planet’s rotation angle of ∼30° during the transit (see Fig. 6). This allows us to explore variations within the morning and evening side as the transit progresses and compare any differences.

Figure 6.

Schematic showing the different regions of the morning (leading) and evening (trailing) limb probed during the transit. During ingress, the stellar light passes through the day side of the morning and the night side of the evening, which then changes as the planet rotates during its transit. During egress, the stellar rays pass through more of the say side of the evening and the night side of the morning.

To compare our results with HyDRA-2D, we also perform a spatially homogeneous retrieval across the combined |$\phi = -0.04 \ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| ranges, which we refer to as the combined-phase 1D-retrieval. We also perform two additional spatially homogeneous retrievals for the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| ranges separately (referred to as the separated-phase 1D-retrievals). Such spatially homogeneous retrievals are well established and have been used to analyse HST as well as ground-based high-resolution observations (e.g. Brogi & Line 2019; Gandhi et al. 2019; Line et al. 2021). We assume uniform Fe abundance across the whole terminator for these 1D-retrievals, with a single pressure-temperature profile for the atmosphere and uniform opacity deck and haze. Hence for each retrieval we have 13 free parameters, including ΔK_p, ΔV_sys, and δV_wind. When applying the spectral broadening to these models we take the full rotational kernel for both halves of the atmosphere. Further details on the retrieval setup for the 1D retrievals can be found in Gandhi & Madhusudhan (2018) and Gandhi et al. (2019).

3 RESULTS AND DISCUSSION

In this section we discuss the results of the HyDRA-2D retrievals on the WASP-76 b observations. We show the constraints for Fe in Fig. 7 and the constraints on the temperature profile and opacity deck in Fig. 8. These indicate that the median values and their errors vary significantly between the evening and morning, and between the retrievals across the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| phase ranges. Table 2 also shows that HyDRA-2D offers a better statistical fit the observations for the |$\phi = -0.04\ {\rm to} \ -0.02$| range over traditional spatially homogeneous approaches. This therefore highlights the need for such modelling capabilities in order to effectively constrain atmospheric parameters and resolve spatial differences between the two sides. On the other hand, we find that for the last quarter of the transit the evening side dominates the signal, with little contribution from the morning side, and hence the 1D retrieval shows almost equal evidence. We discuss the constraints for each of the various parameters in detail below.

Figure 7.

Retrieved Fe constraints from HyDRA-2D for each side of the planet and for each of the retrievals conducted over the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| ranges. The orange and red histograms give the constraints for the evening, and the purple and blue show the constraints for the morning. The blue bar indicates the median and ±1σ error bar on the retrieved abundances. We also show the retrieved constraints from the separated-phase 1D retrievals that assume a homogeneous terminator with the dashed black histogram, with the grey bar indicating its 1σ confidence interval. We additionally plot the solar value of Fe with the red dashed vertical line, with log (X_Fe) = −4.52 (Asplund et al. 2009).

Figure 8.

Constraints on the Fe, temperature profile, and opacity deck for WASP-76 b with HyDRA-2D. At the top of the figure, we show a schematic of the region of the atmosphere probed by the retrievals, with the corresponding constraints for each region highlighted below. For the Fe histograms the blue error bar indicates the 1σ confidence interval and the red dashed line shows the solar value of Fe (see Fig. 7). The bottom panels show the 1σ and 2σ confidence interval for the P – T profile in the darker and lighter shades, respectively, with the median value given by the solid curve. Similarly, we highlight the constraints on the opacity deck and its 1σ and 2σ confidence interval in the darker and lighter grey. Note that this is the uncertainty in the pressure of the top of the opacity deck, and the observed spectrum is not sensitive to higher pressures as the opacity deck is optically thick.

Table 2.

Retrieved Fe abundances with HyDRA-2D for each side of the atmosphere and with the 1D retrievals for each phase range (see Fig. 7).

	|${\phi = -0.04 \ {\rm to} \ -0.02}$|		|$\mathbf {\phi = 0.02\!-\!0.04}$|
	Evening	Morning	Evening	Morning
Combined-phase 1D	|$-4.33^{+0.26}_{-0.27}$|
Separated-phase 1D	|$-4.89^{+0.43}_{-0.45}$|		|$-4.20^{+0.26}_{-0.31}$|
Preference over combined-phase 1D	3.8σ
HyDRA-2D	|$-4.97^{+0.63}_{-0.56}$|	|$-3.83^{+0.51}_{-0.59}$|	|$-4.03^{+0.28}_{-0.31}$|	|$-4.59^{+0.85}_{-1.0}$|
Preference over combined-phase 1D	4.9σ
Preference over separated-phase 1D	3.6σ		–

Notes. We also show the model preference for each retrieval in σ. We multiply the overall evidences of the separated-phase retrievals to allow for comparisons with the combined-phase retrievals. We see no preference for HyDRA-2D over a separated-phase 1D retrieval for the |$\phi = 0.02\!-\!0.04$| range.

3.1 Comparison to spatially homogeneous models

We compare our constraints on Fe with more traditional 1D modelling approaches, where homogeneous terminator properties are assumed. We find that our spatially homogeneous combined-phase 1D retrieval obtains strong constraints on the Fe abundance, with |$\log (X_\mathrm{Fe}) = -4.33^{+0.26}_{-0.27}$| (see Table 2). When we perform separate retrievals for each of the phase ranges, we find that the abundance constraints now show some differences. The separate-phase 1D retrieval retrieves |$\log (X_\mathrm{Fe}) = -4.89^{+0.43}_{-0.45}$| for the |$\phi = -0.04\ {\rm to} \ -0.02$| range. On the other hand, the signal from the last quarter of the transit closely matches the Fe abundance constraint from the combined-phase 1D retrieval. Therefore, the combined-phase retrieval constrains the most significant contribution to the Fe signal, which comes from the latter quarter, consistent with previous studies (Ehrenreich et al. 2020; Kesseli & Snellen 2021). Given the differences in the constrained abundance, it is therefore unsurprising that the separated-phase 1D retrievals are statistically preferred by 3.8σ. This indicates that the overall atmospheric structure is not as well modelled by assuming homogeneous physical properties across the transit.

The HyDRA-2D constraints go further by exploring the variation in Fe over the terminator for each of the phase ranges. The evening sides of the spatially resolved retrievals show remarkably similar constraints to the separated-phase 1D retrievals (see Fig. 7 and Table 2). For both phase ranges the evening side is thus the dominant side of the planet’s Fe signal, which remains true for the |$\phi = -0.04\ {\rm to} \ -0.02$| range even though the evening side’s abundance is constrained to be lower than for the morning. This is because the spectral contribution resulting from a lower abundance is offset by the much deeper opacity deck for the evening side (see Section 3.4 and Fig. 8), thus maintaining a stronger contribution to the overall signal. The statistical preference for HyDRA-2D and the separated-phase 1D retrievals over the combined-phase retrieval is driven primarily by differences in the temperature and the wind speeds (see Sections 3.3 and 3.5), with a weaker effect caused by the chemical inhomogeneities of Fe in the atmosphere. This is consistent with previous work by Beltz et al. (2021) comparing high-resolution observations in emission with general circulation models.

We also calculate terminator-averaged abundances from HyDRA-2D and compare these to our separated-phase 1D retrievals for each of the phase ranges we study, as shown in Fig. 9. We do this to investigate the potential biases that may arise by assuming spatial homogeneity. We calculate the average abundance of the terminator from HyDRA-2D through two methods. First, we use the arithmetic mean of Fe across the two terminators. Our second method employs a more sophisticated temperature-weighted geometric mean, which has been shown to provide a better prescription for averaged terminator abundances for low-resolution observations (Welbanks & Madhusudhan 2022). This calculation does implicitly assume that the cross-section between the two sides does not vary significantly, as we do not weight each side by the overall opacity for each half of the terminator. However, we can justify this given that we do not expect significant variation in the opacity over the range of the temperature constraints.

Figure 9.

Averaged posterior distribution of HyDRA-2D over the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| ranges. We calculate an arithmetic mean (shown in red) and a temperature-weighted geometric mean based on the prescription in Welbanks & Madhusudhan (2022) (shown in blue). We also show the results from the separated-phase 1D retrievals with the dashed black line. The red dashed line indicates the solar value of Fe.

The averaged values are consistent with each other for both phase ranges given the relatively little difference in abundance between the two sides. Our 1D retrieval is in fact most closely matching the evening side abundance from HyDRA-2D (see Fig. 7). This is likely due to the fact that with the high-resolution R = 140 000 observations of ESPRESSO we obtain separate signals from each side of the atmosphere due to their difference in velocity. Hence the signals from each side of the atmosphere are not completely blended together, as would be the case for low-resolution observations, and thus for our high-resolution observations the spatially unresolved retrievals converge to the stronger signal. The last quarter of the transit (⁠|$\phi = 0.02\!-\!0.04$|⁠) shows similar constraints for both methods of averaging the abundance as the constraint from the morning is relatively weak and thus dominated by the evening side.

3.2 Variation in Fe

The Fe is most strongly constrained from HyDRA-2D for the evening side of the atmosphere for the |$\phi = 0.02\!-\!0.04$| range (see Fig. 7). We see a tight well-constrained abundance of |$\log (X_\mathrm{Fe}) = -4.03^{+0.28}_{-0.31}$|⁠, slightly above the solar value of log (X_Fe) = −4.52 (Asplund et al. 2009), but consistent with the metallicity of the host star, [Fe/H] = 0.366. However, our highest abundance estimate for Fe in fact comes from the morning of the |$\phi = -0.04\ {\rm to} \ -0.02$| range, with |$\log (X_\mathrm{Fe}) = -3.83^{+0.51}_{-0.59}$|⁠. Although, it should be noted that this estimate has a higher uncertainty. Both of these regions of the atmosphere with the highest abundance estimates probe more of the day side of the atmosphere, i.e. that is being irradiated due to the geometry of the transit (see Fig. 8). Hence, it is unsurprising that we obtain the tightest constraints for these regions, as the day side is expected to have the strongest Fe signal given the hotter temperatures, as discussed in Section 3.3.

Our morning side constraint for the |$\phi = 0.02\!-\!0.04$| range has the lowest abundance and the weakest constraint, at |$\log (X_\mathrm{Fe}) = -4.59^{+0.85}_{-1.0}$|⁠, as the evening side dominates the overall signal. The stellar limb darkening also acts against the signal from the morning side for the latter half of the transit, as this side occults more of the dimmer limbs of the star. Given the wider uncertainty for the morning side, it is unsurprising that the 1D retrieval shows a similar evidence and constrains an Fe abundance that closely matches the evening. Further observations of the transit of WASP-76 b may constrain the Fe signal from the morning more definitively, particularly with high-resolution (R∼140 000) observations that are able to resolve the leading and trailing limbs.

We find that the abundance constraints for the sides probing more of the night side, i.e. |$\phi = -0.04\ {\rm to} \ -0.02$| evening side and |$\phi = 0.02\!-\!0.04$| morning side, have lower median values by ∼0.8 dex, indicating that some Fe rain-out may be occurring across the terminator, as suggested by Ehrenreich et al. (2020) for the morning side. This rain-out is also consistent with recent work on cloud formation on UHJs by Gao & Powell (2021) that showed that the Fe abundance could potentially decrease by ∼1 dex due to it condensing. However, we note that the chemical abundance differences between the various regions of the atmosphere are generally quite small and broadly consistent with each other, and differences in the temperature and wind speeds have a higher impact on the overall spectrum.

3.2.1 Mass–metallicity relation

Fig. 10 shows our most precise abundance estimate for Fe, that of the evening side for the last quarter of the transit, along with previous measurements of various species in Solar System planets and other exoplanets. As chemical models of the atmospheres of such hot objects have shown (e.g. Visscher, Lodders & Fegley 2010), gaseous iron is expected to be the dominant proportion of the Fe, and thus a good estimate of the total Fe content of the atmosphere. The Fe measurement of an exoplanet atmosphere also provides a direct comparison against the stellar abundance (e.g. Thiabaud et al. 2015). Our retrieved value lies in between the stellar metallicity and the trend from the Solar system planets. In addition to WASP-76 b, we also show the Fe abundance measured for WASP-121 b (Gibson et al. 2020), an UHJ with an equilibrium temperature of ∼2400 K, and the CO abundance constraints from other high-resolution studies (Gandhi et al. 2019; Line et al. 2021; Pelletier et al. 2021). These show good agreement with our result, as well as the previous work of Welbanks et al. (2019), which found species other than oxygen at stellar- or superstellar abundances.

Figure 10.

Mass–metallicity relation for known exoplanets and the Solar System, corrected for the stellar metallicity of each system. The constraint from this work is shown in black and uses our most precise value of |$\log (\mathrm{Fe}) = -4.03^{+0.28}_{-0.31}$|⁠, which arises from the evening side for the last quarter of the transit. The orange region indicates the Solar System trend of decreasing metallicity with mass. We also show the Fe abundance measurement in the UHJ WASP-121 b (Gibson et al. 2022) in red, and H₂O/H (blue) and CO/H (green) constraints for a range of exoplanets obtained with low- and high-resolution observations, respectively (Welbanks et al. 2019; Guillot et al. 2022).

The Fe abundance is of particular importance given that it can be used to determine the refractory-to-volatile abundance ratio. This provides a unique window into the formation and migration scenarios (Lothringer et al. 2021). As the Fe accreted on to the planet must be in the form of solids, it allows us to break the degeneracy often seen with volatile species, where the carbon and oxygen bearing species may be accreted as either solids or gases (e.g. Öberg, Murray-Clay & Bergin 2011; Mordasini et al. 2016). Further observations of UHJs and constraints on their refractory species will help us to discern the trends and allow us to compare with the volatile abundances measured from observations in the infrared.

3.3 Variation in the thermal structure

Our retrieved temperature and its constraints do show variation between the evening and morning sides as well as with phase, as shown in Fig. 8. These show that the temperature is best constrained for the evening side towards the end of the transit, where we retrieve a hot photosphere of |$T_\mathrm{0.1mb} = 2950^{+111}_{-156}$| K. On the other hand, we obtain a lower median temperature with a wider uncertainty of |$T_\mathrm{0.1\,mbar} = 2615^{+266}_{-275}$| K for the morning side for the |$\phi = 0.02\!-\!0.04$| phase range. The constraints for both phase ranges show that the parts of the atmosphere that are irradiated show the highest temperatures, in line with expectations from general circulation models (Wardenier et al. 2021; Savel et al. 2022). For all of our temperature constraints, the profiles also appear largely isothermal, and generally show a widening uncertainty at pressures higher than the opacity deck where the observations are unable to probe. GCMs indicate thermal inversions for the upper photosphere for the day side, although these inversions are generally consistent with the retrieved uncertainty in our temperature profiles for all cases except the evening side for the |$\phi = 0.02\!-\!0.04$| range, where our error is more tightly constrained, given the higher signal-to-noise ratio.

We compare our temperature constraints to previous analyses and other observations of the primary transit of WASP-76 b. High-resolution observations using the CARMENES spectrograph constrained a higher temperature of ∼2700–3700 K range, consistent with our temperature profiles (Landman et al. 2021). Recent analyses with the ESPRESSO and HARPS observations in the optical detected a temperature of ∼3300–3400 K (Seidel et al. 2019, 2021), which is higher than our temperature for even the hottest part of the atmosphere, namely the evening side for the last quarter of the transit. However, Von Essen et al. (2020) used low-resolution HST observations in the optical and found temperature constraints of |$2300^{+412}_{-392}$| K, while Edwards et al. (2020) analysed HST WFC3 observations in the infrared and obtained |$2231^{+265}_{-283}$| K. Fu et al. (2021) combined the WFC3 and Spitzer observations and found that the best-fitting temperature was 2000 K from forward models. These low-resolution constraints point to a temperature that is generally lower than we observe, but still consistent at ∼1.5σ with our constraint for the morning side 0.1-mbar temperature for both phase ranges. The differences may be reconciled, given that the high-resolution observations in our work are probing a higher altitude than the low-resolution observations and thus may obtain a higher temperature if there is a stratosphere present, as indicated from day-side analysis of HST observations (Fu et al. 2021). In future, combining the high-resolution and low-resolution observations may help to constrain the thermal profile to a greater extent (e.g. Brogi & Line 2019; Gandhi et al. 2019) and shed light on the potential differences seen in the constrained temperature of the terminator.

We have performed retrievals assuming an H₂/He- and an H/He-rich atmosphere, and find that the H/He-rich case (where the H₂ is dissociated) offers a substantially better fit to the observations. This is because the atmospheric composition is intricately linked to the overall constrained temperature. Such a dissociated atmosphere, consisting primarily of H and He, would be a lower mean molecular weight over an H₂- and He-rich atmosphere by a factor of almost 2. Therefore, the scaleheight, given by H = k_bT/μg, is almost double for an H/He-rich atmosphere over an H₂/He-rich atmosphere. Modelling of circulation on UHJs has shown that the majority of the day-side atmosphere is H₂-dissociated at pressures ∼0.1 bar (Bell & Cowan 2018; Tan & Komacek 2019). In the upper atmosphere (∼10⁻⁴ bar) where our observations are more sensitive, we expect a significant majority of the hydrogen in the atmosphere to be thermally dissociated, as shown in Fig. 11. Therefore, our retrievals with an H₂/He-rich atmosphere constrained unphysically high temperatures of ≳ 3800 K in order to match the strength of the features in the observations. Furthermore, we also find that the H₂/He-rich retrievals were statistically disfavoured, and thus find evidence for an atmosphere that is dissociated. The recombination of H into H₂ on the morning side of the last quarter of the transit may also be a potential reason for its lack of significant signal compared to the evening side, as its scaleheight would be decreased by a factor of ∼2 over the evening in addition to any difference in temperature and Fe abundance.

Figure 11.

Fraction of H₂ to (H + H₂) for typical photospheric pressures and temperatures for UHJs. These were derived using the analytic prescription in Parmentier et al. (2018) (see also Gandhi, Madhusudhan & Mandell 2020a). The lowest pressures with the highest temperatures have the greatest degree of H₂ dissociation, resulting in atmospheres that are H-rich.

We do observe a weak thermal inversion for the evening side in the |$\phi = -0.04\ {\rm to} \ -0.02$| range. This may be driven by a potential thermal inversion from the day side remaining present at the terminator, and detectable, given the deeper opacity deck for this part of the atmosphere (see Section 3.4). Numerous strong optical absorbers are capable of generating thermal inversions in UHJs (e.g. Lothringer, Barman & Koskinen 2018; Gandhi & Madhusudhan 2019). We note however that the retrieved inversion may also be caused by the recombination of H₂ in the deep atmosphere. At high pressures, the H₂ is present at significant abundance as shown in Fig. 11. This increases the mean molecular weight and thus reduces the scaleheight. Hence our retrieval with the H/He-rich atmosphere may act to reduce the apparent scaleheight in the lower atmosphere by lowering the temperature at higher pressures, and thus appear as a thermal inversion. We tested our temperature profile parametrization was not the root cause of this inversion by testing various pressures where the temperature was set instead of 0.1 mbar, and also varying the lower prior on P₁ and P₂, but these did not impact our results. We do note however that the inversion is weak (∼300 K) and consistent with an isotherm as these high-resolution terminator observations are only indirectly sensitive to the temperature profile from its effect on the scaleheight. We performed an additional retrieval that did not allow thermal inversions in the atmosphere and observed no significant change to our parameter estimates, further indicating that the inversion not too significant.

3.4 Opacity deck

The constraints for the opacity deck are given in Fig. 8. The range in pressures constrained for the top of the cloud deck in Fig. 8 represents the uncertainty, and the spectrum is not sensitive to pressures below this as the opacity deck is assumed to be optically thick. The retrievals show that there is a high-opacity deck/continuum present for the sides of the atmosphere that are probing most strongly the day side of the atmosphere, and our signal for Fe is probing the upper atmosphere above the opacity deck. The source of this opacity likely arises from species such as H- and/or other atomic species that may be prevalent in the atmosphere, as condensation of cloud species is unlikely on the day side for such a hot day side. This is backed up by the numerous species that have also been detected in the atmosphere with significant optical opacity (Tabernero et al. 2021; Kesseli et al. 2022). We note that species such as TiO or VO also have strong optical opacity and are capable of leading to thermal inversions (Hubeny, Burrows & Sudarsky 2003; Fortney et al. 2008; Piette et al. 2020), but there is no evidence for these species from the high-resolution observations (Tabernero et al. 2021). On the other hand, our results for the non-irradiated regions of the atmosphere show a deeper opacity deck, indicating that the night sides may be relatively free of any strong optical absorbers. Cloud formation modelling has shown that the night side should possess significant cloud opacity due to silicate clouds (Gao & Powell 2021). Such clouds form at pressures of 10⁻²–10⁻³ bar, depending on the temperature, consistent with our opacity decks for the night side of |$\log (P_\mathrm{cl}/\mathrm{bar}) = -2.08^{+0.44}_{-0.53}$| (evening side, |$\phi = -0.04\ {\rm to} \ -0.02$|⁠) and |$\log (P_\mathrm{cl}/\mathrm{bar}) = -3.33^{+0.80}_{-0.91}$| (morning side, |$\phi = 0.02\!-\!0.04$|⁠). Recent work has also shown that patchy clouds may be present on the night side of UHJs but that they only weakly affect terminator observations (Komacek et al. 2022).

3.4.1 Abundance-opacity deck degeneracy

Our constrained Fe abundance shows a degeneracy with the pressure level of the opacity deck, as shown in Fig. 12. This shows that a higher altitude opacity deck requires a higher abundance in the atmosphere, and arises because the column depth of a given species increases to match the observed feature strength when the cloud deck altitude is raised (e.g. Bétrémieux & Swain 2017; Welbanks & Madhusudhan 2019). Recent work on characterizing atmospheres with cloud/opacity decks at high resolution (Gandhi, Brogi & Webb 2020c) indicated that abundances of chemical species are still constrainable even in the presence of high-altitude opacity decks. Our R = 140 000 observations with EPSRESSO confirm this is the case and we obtain relatively tight constraints for Fe despite our limited spectral coverage of just ∼0.4 |$\mu$|m.

Figure 12.

Posterior distribution of our retrievals of WASP-76 b highlighting the degeneracy between Fe abundance and opacity deck. The left-hand panels in red show the degeneracy for the evening side of the first quarter retrieval (⁠|$\phi = -0.04\ {\rm to} \ -0.02$|⁠), and the right-hand panels shows the constraints for the morning side for the last quarter (⁠|$\phi = 0.02\!-\!0.04$|⁠).

We find that the Fe abundance retrieved for the morning side for |$\phi = 0.02\!-\!0.04$| has a wider uncertainty due to the abundance-opacity deck degeneracy. Ehrenreich et al. (2020) proposed that Fe rains out of the atmosphere and thus we are unable to observe it for the morning side later in the transit. On the other hand, GCMs have shown that the lack of Fe signal may be due to the presence of a cloud deck obscuring the signal from Fe (Savel et al. 2022). Our retrievals indicate a slight preference for the latter, but the constraints are weak due to its weaker signal compared to the much hotter evening side. We also see no preference for our spatially resolved retrieval over a 1D spatially homogeneous retrieval for the |$\phi = 0.02\!-\!0.04$| range (see Table 2), indicating that the Fe abundance of the morning may not be significantly different to that from the evening.

3.5 Constraints on winds

The constraints on the wind parameters from HyDRA-2D are shown in Fig. 13 and the corresponding constraints on the wind speed distributions are plotted in Fig. 14. In our case, the wind speed is given by ΔV_sys as we have shifted the spectra to the rest frame of the planet during our analysis, and the uncertainty in the known value of V_sys is negligible. The ΔV_sys values are dependent on the orbital period and transit time centre, but by propagating the uncertainties reported in Ehrenreich et al. (2020), we find that ΔV_sys changes by at most ∼0.5 km s⁻¹. By using other recently published precise orbital values (Ivshina & Winn 2022; Kokori et al. 2022), we also find consistent results to within the uncertainties and ∼1 km s⁻¹. Therefore, any redshift or blueshift of the spectrum indicates the presence of a wind and cannot be explained by variations in orbital parameters.

Figure 13.

Posterior distributions of ΔK_p and the wind parameters, ΔV_sys and δV_wind, for each HyDRA-2D retrieval of WASP-76 b. The red posteriors show the constraints on the |$\phi = -0.04\ {\rm to} \ -0.02$| range and the blue posteriors show the |$\phi = 0.02\!-\!0.04$| range. We also provide the median and ±1σ values for each retrieval in the table on the bottom left-hand side.

Figure 14.

Normalized wind speed distribution as derived from the ΔV_sys and δV_wind parameters given in Fig. 13. This is the day-night wind broadening kernel applied to the model as shown in Fig. 2. The constraints on the |$\phi = -0.04\ {\rm to} \ -0.02$| range are shown in red and the constraints on the |$\phi = 0.02\!-\!0.04$| range in blue. The dashed lines indicate the median values of ΔV_sys from our retrieval. The negative values indicate a day-to-night wind, which results in a net blueshift of the overall spectrum.

For both phase ranges, we see a strong day-night wind with a significant blue shift to the overall spectrum, consistent with previous observations and analyses (Ehrenreich et al. 2020; Kesseli & Snellen 2021). This is a result of the day-night wind that travels toward the observer at the terminator. For the first quarter of the transit we retrieve a wind speed of |$-5.9^{+1.5}_{-1.1}$| km s⁻¹, with the negative sign indicating a blueshift (i.e. a wind travelling from the day-to-night side). This value is consistent with previously derived wind speeds from the detection of Na in the upper atmosphere of WASP-76 b (Seidel et al. 2021). However, we constrain a stronger wind speed of |$-9.8^{+1.2}_{-1.1}$| km s⁻¹ for the |$\phi = 0.02\!-\!0.04$| range, but Na is not as strongly detected in this phase range (Kesseli et al. 2022) and thus could explain why the wind speeds are not in as good an agreement.

As well as the difference in the wind speeds, we also see a difference in the δV_wind parameter, which measures the FWHM of the spread in wind speed. We constrain a smaller δV_wind for the last quarter of the transit, as shown in Fig. 14. The differences in the parameters during the transit could potentially be due to a combination of a more dominant day-night wind and a weaker east–west jet. As the signal from the |$\phi = 0.02\!-\!0.04$| range is dominated by the evening side, any jet present would add to the overall profile and result in a greater blue shift of the signal, resulting in a lower ΔV_sys. In addition, the spread in the δV_wind would also be lower as the morning side contributes very little to the overall signal. On the other hand, the |$\phi = -0.04\ {\rm to} \ -0.02$| range has a significant contribution from both sides, and may therefore result in a lower overall wind speed and greater spread as any contribution from an east–west jet will be redshifted for the morning. In our constraints we do retrieve a wind profile that has a small redshifted contribution for the first quarter of the transit. A day-night wind of ∼6 km s⁻¹ combined with a jet of ∼4 km s⁻¹ could thus explain the relative shift in ΔV_sys and spread in δV_wind between the first and last quarter of the transit. We have also carried out retrievals with an east–west jet instead of a day-night wind and did not see strong evidence for such a case, thus any jet present would be weak compared to the day-night wind. This is consistent with the GCMs from Wardenier et al. (2021), who showed that such a weak-drag model best describes the observations.

On the other hand, the higher opacity decks constrained for both the morning and evening terminators for the |$\phi = 0.02\!-\!0.04$| range may also result in such a difference in wind parameters between the first and last quarter of the transit. If the observations in the latter half of the transit probe higher altitudes, there may be more of a blue shift in the spectrum as the day-night wind is expected to be highest at lower pressures. This would also reduce the overall spread in the wind speed, δV_wind, as the range in pressures probed reduces. Hence it is currently unclear which of these two potential causes of the difference in wind speed is dominant, and further observations of WASP-76 b at high resolution will help to shed light on the various dynamical processes occurring in the atmospheres of UHJs.

We note that there is a strong degeneracy between ΔK_p and ΔV_sys in Fig. 13, as expected given the limited phase range of the observations that would act to break this. This degeneracy results in the two ΔK_p–ΔV_sys distributions for the first and last quarter, which tend towards each other at lower K_p values. However, the true K_p value is constrained to within ∼1 km s⁻¹ (Ehrenreich et al. 2020), and so we believe the ΔV_sys offset between the first and last quarter of the transit is a result of the observations and not uncertainty in the K_p. The ΔK_p is consistent with expected value of 0 for the |$\phi = 0.02\!-\!0.04$| range, but the |$\phi = -0.04\ {\rm to} \ -0.02$| range does indicate a lower ΔK_p, albeit still consistent with 0 at ≲ 2σ. The deviation from the expected value is be due to the signal changing its velocity with phase for the first part of the transit, as seen in previous studies with other data sets and GCMs (e.g. Kesseli & Snellen 2021; Wardenier et al. 2021). This phase dependence is not modelled in our retrievals due to the complexity of including multi-dimensional phase-dependent radiative transfer into our calculations, and thus results in a slight change in ΔK_p to compensate.

3.6 Comparison to GCMs

We compare the results from HyDRA-2D to general circulation models of WASP-76 b in Fig. 15. These show the asymmetric weak-drag Monte Carlo radiative transfer model of Wardenier et al. (2021) (see ‘modification 1’ in their figs 12 and 13) at phases of ϕ = −0.031 and 0.031, along with the best-fitting retrieval model for the |$\phi = -0.04\ {\rm to} \ -0.02$| and |$0.02\!-\!0.04$| range for a small region near ∼0.366 |$\mu$|m. For both of the best-fitting models, we have removed the opacity deck as Fe is expected to be the dominant source of opacity in this spectral range, and we do not expect significant cloud condensation or hazes given the high temperature. We have also mean subtracted the data sets to remove the effect of varying planetary radius and/or reference pressure, as these dependencies are removed in our retrieval in the high-pass filter as part of the analysis.

Figure 15.

Comparison of HyDRA-2D to general circulation models of WASP-76 b. The highest likelihood model from the retrieval of the |$\phi = -0.04\ {\rm to} \ -0.02$| range was used to generate a transit spectrum as shown in the left-hand panel (in red), and the right-hand panel shows the the corresponding best-fitting model for the |$\phi = 0.02\!-\!0.04$| range. The general circulation model assumes an atmosphere with weak drag, and is calculated through Monte Carlo radiative transfer (see Wardenier et al. 2021, for further details).

Fig. 15 shows that for the first quarter of the transit, the two models agree well. This is crucial as 3D models have shown better fits to high-resolution observations (Flowers et al. 2019; Beltz et al. 2021), and thus HyDRA-2D is able to capture the relevant physics to fit the dynamics and variability in the observations. The line positions agree well and have similar strengths between our model and the GCM, indicating that the constrained wind speeds, abundances and temperatures are in line with expectations. We do however find that our retrieved best-fitting models are somewhat smoother than the predictions made from GCMs. This difference is because we model an atmosphere with a latitudinally homogeneous day-night wind contribution, but the GCM with weak drag includes latitude-dependent effects, depth-dependent wind speeds, as well as a temperature structure that varies along the line of sight. These additional effects result in additional structure to the GCM spectra. We do also see some small differences in the strengths of some of the Fe lines, but the overall agreement between the retrieval and GCM is encouraging.

The latter half of the transit shows that the retrieved features are weaker than the GCM predictions. The GCM indicates that the strength of the Fe signal should increase, but the retrieved best-fitting model from HyDRA-2D for ϕ = 0.031 does not increase in strength over the |$\phi = -0.04\ {\rm to} \ -0.02$| range. A potential factor driving this may be the temperature profile. The retrieval constrains temperatures nearer ∼3000 K for the evening side, but the GCM shows some regions of the atmosphere reaching temperatures nearer to ∼3500 K. This difference results in larger spectral features in the GCM spectra by ∼15 per cent due to the changing scaleheights of the atmosphere. There will also be an additional and potentially stronger effect of varying line strengths with temperature, and the weaker Fe lines being further reduced in the retrieval due to their strong temperature dependence as shown in Fig. 15. We note that as the retrieval is a quantitative fit to the data, the GCM may not be completely capturing the true state of the atmosphere and overestimating the temperature. This temperature difference between the retrieval and GCM may be be driven by various effects such as a supersolar Fe abundance, other optical absorbers as well as disequilibrium chemistry in the atmosphere. Further observations of the transit of WASP-76 b will help shed more light on these differences and help improve our 3D modelling of such hot planets.

3.7 Robustness tests

To verify the robustness of our results we performed numerous checks with our model. We carried out retrievals with a scalefactor that rescales the model spectrum to account for uncertainties in the measured value of the radius and mass and any non-hydrostatic effects that may be present in the upper atmosphere. Such a parameter has often been included in previous high-resolution retrievals (e.g. Gandhi et al. 2019; Gibson et al. 2020; Line et al. 2021). The scale factor did show some degeneracy with the temperature as both act to increase/decrease the overall scaleheight of the atmosphere, but there were no differences in the retrieved parameters and the overall relative evidences remained similar. We also performed separate retrievals on the |$\phi = 0.02\!-\!0.035$| and |$0.025\!-\!0.04$| range, and |$\phi = -0.035\ {\rm to} \ -0.02$| and |$-0.04\ {\rm to} \ -0.025$| ranges. The top pressure that we modelled our atmosphere out to was also varied from 10⁻⁷ to 10⁻⁸ bar. As well as these tests, we performed retrievals with additional free parameters for the radius and log (g) of the planet. None of the retrievals showed any significant difference to any parameter, and our retrievals with log (g) constrained estimates consistent with the known values (Ehrenreich et al. 2020). We were unable to constrain the planetary radius to any significant degree and thus simply recovered the prior given that the high-pass filtering of the data analysis removes the radial dependence of the spectrum. We also performed retrievals varying the number of live points with MultiNest and did not see any significant change in the convergence of any of the parameters.

We also performed additional retrievals parametrizing the temperature profile as an isotherm and with no super-Rayleigh haze prescription, and found no significant changes to our constraints on the photospheric temperature, Fe abundance, opacity deck or wind parameters. The evidences for the isothermal retrievals was also similar, with the largest difference coming from the HyDRA-2D retrieval for the first quarter, given that a weak thermal inversion was constrained for the evening side. However, this was <2σ deviant from the retrievals with the fully flexible P – T profile. These tests indicate that the observations are more sensitive to the overall temperature than to the gradients in temperature with height, in line with expectations for primary eclipse observations.

The only change we saw in the retrieved constraints was when we included the very bluest end of the data between 0.38 and 0.4 |$\mu$|m. For these tests we retrieved a higher temperature for the atmosphere and a negative velocity shift in ΔK_p. Given the very strong Fe opacity, these wavelengths probe regions of the atmosphere ≲ 10⁻⁷ bar where non-local thermodynamic effects are likely to be considerable, and our assumptions of ideal gas and hydrostatic equilibrium breakdown. Hence, it is unsurprising that our retrieval obtains a lower overall evidence for this wavelength region. We also varied the wavelength range we use for the redder end of the observations, with retrievals in the 0.4–0.65 and 0.4–0.7 |$\mu$|m range as well as our fiducial 0.4–0.8 |$\mu$|m. These showed no change in the retrieved constraints for any parameter, despite some telluric absorption from O₂ and H₂O at wavelengths ≳ 0.7 |$\mu$|m. This confirms that the spectral cleaning steps were able to effectively remove the tellurics in this wavelength range (see Section 2.5).

4 CONCLUSIONS

We have explored variations in Fe abundance, temperature profile, and opacity/haze deck in the atmosphere of the UHJ WASP-76 b through a new spatially resolved high-resolution retrieval of the terminator. This was motivated by recent analyses on ESPRESSO and HARPS observations of the primary transit that indicated the presence of an asymmetric signal in Fe during the transit. With our new retrieval model, HyDRA-2D, we explored the variability and constrained a day-night wind in the photosphere. We find that Fe has the tightest constraints in the observations on the trailing (evening) limb of the last quarter of the transit (⁠|$\phi = 0.02\!-\!0.04$|⁠), with an abundance of |$\log (X_\mathrm{Fe}) = -4.03^{+0.28}_{-0.31}$|⁠. The abundance estimates of Fe show that the atmosphere is consistent with the stellar metallicity and the CH₄/H measurement of Jupiter, as well as with previous constraints of Fe on the UHJ WASP-121 b (Gibson et al. 2022).

Our HyDRA-2D models are statistically preferred by 4.9σ over more traditional spatially homogeneous and combined-phase models. The different regions of the atmosphere probed are distinguishable given the high precision and resolution of the observations made with the ESPRESSO spectrograph. The statistical preference is driven primarily by differences in the temperature and winds, with a weaker effect also caused by the chemical abundance differences between the different regions of the atmosphere. This is consistent with previous 3D modelling of high-resolution emission spectra in the infrared (Beltz et al. 2021). Furthermore, for our separated-phase retrievals, we find that HyDRA-2D is preferred by 3.6σ over a spatially homogeneous approach for the |$\phi = -0.04\ {\rm to} \ -0.02$| range. Hence, adopting a spatially homogeneous model is likely to miss signals from some regions of the atmosphere, particularly when both sides have a high signal-to-noise ratio and there is significant asymmetry in the physical parameters. On the other hand, we find that our 1D-retrieval for the |$\phi = 0.02\!-\!0.04$| range has nearly identical evidence to the spatially resolved HyDRA-2D model, as a result of the evening side dominating over the morning side for the latter quarter of the transit.

Generally, we observe Fe abundances that are consistent between the different regions of the terminator, but we do find that the Fe abundance is highest where the region of the planet probed is being irradiated, namely the morning side for |$\phi = -0.04\ {\rm to} \ -0.02$| and the evening side for |$\phi = 0.02\!-\!0.04$|⁠. We also find that our spatially homogeneous retrievals for the first and last quarter of the transit converge towards the constraints of the evening side. This suggests that the evening side is the dominant source of the Fe signal. This remains true for the first quarter of the transit despite the higher constrained abundance for the morning side, as this higher abundance is offset by the higher altitude opacity deck in the atmosphere, resulting in a weaker overall contribution.

We do not constrain Fe as significantly on the morning side for the |$\phi = 0.02\!-\!0.04$| range, with the retrieved abundance, |$\log (X_\mathrm{Fe}) = -4.59^{+0.85}_{-1.0}$|⁠, having a wider uncertainty than for any other part of the atmosphere. Ehrenreich et al. (2020) proposed that the Fe signal is most prominent on the evening, and the lack of significant signal from the morning side may be due to Fe condensing out of the atmosphere. On the other hand, Savel et al. (2022) proposed that the signal was diminished due to a high-altitude cloud deck. Our retrievals show that either scenario can fit the observations given that there is a degeneracy between abundance and the opacity deck, with a slight preference for a high abundance and opacity/cloud deck. The degeneracy is only partial as recent work has shown that abundance constraints from high-resolution observations are attainable even with high-altitude opacity/cloud decks (Gandhi et al. 2020c), and we confirm this with observational data. However, our Fe uncertainty for the morning for the |$\phi = 0.02\!-\!0.04$| range is much wider than for other regions of the atmosphere because the signal from the evening side dominates the overall spectrum, as suggested from GCMs (Wardenier et al. 2021).

Our constraints on the temperature profile indicate an atmosphere that is between |$2950^{+111}_{-156}$| K (⁠|$\phi = 0.02\!-\!0.04$|⁠, evening side) and |$2615^{+266}_{-275}$| K (⁠|$\phi = 0.02\!-\!0.04$|⁠, morning side) at 0.1 mbar. We generally find a temperature that is hottest for the parts of the atmosphere that are irradiated, i.e. the evening side for the last quarter of the transit and the morning side for the first quarter. Our temperature profiles do not significantly constrain the temperature gradient of the atmosphere, but this is expected given that we do not directly detect photons from the planet in primary transit thus making constraints on the profile more difficult. We do however see a strong dependence of the temperature on thermal dissociation, a dissociated atmosphere consisting of H and He was statistically preferred and gave more physically plausible temperatures over an undissociated atmosphere of H₂ and He. This is because an undissociated atmosphere has a smaller scaleheight and thus the temperature must substantially increase in order to compensate and match the features in the spectrum. Hence we infer the presence of H₂ dissociation in the atmosphere of WASP-76 b, consistent with predictions made for UHJs (e.g. Parmentier et al. 2018).

We constrain a much stronger day-night wind at the end of the transit at |$-9.8^{+1.2}_{-1.1}$| km s⁻¹, where the negative sign denotes that the wind is blueshifted and hence travelling towards the observer at the terminator. For the last quarter of the transit, the signal is dominated by the evening side. On the other hand, the first quarter of the transit includes strong contributions from both sides of the atmosphere, and we find that the wind speed is reduced to |$-5.9^{+1.5}_{-1.1}$| km s⁻¹. The FWHM of the wind speed distribution, δV_wind, is also more spread for the |$\phi = -0.04\ {\rm to} \ -0.02$| range at |$10.2^{+1.0}_{-1.0}$| km s⁻¹ than for the |$\phi = 0.02\!-\!0.04$| range at |$5.62^{+0.62}_{-0.64}$| km s⁻¹. These constraints could indicate the presence of a weak equatorial jet in the photosphere, as this would induce a stronger blue shift of the evening side and combine with the day-night wind to drive the overall wind speed to higher values for the |$\phi = 0.02\!-\!0.04$| range. However, the stronger wind speed and lower δV_wind could also be a result of the higher altitudes probed in this phase range, given that the day-night wind is expected to be stronger at higher altitudes (e.g. Wardenier et al. 2021).

This work is a significant step into multi-dimensional spatially resolved retrieval analyses of exoplanetary atmospheres. Key to these retrievals is incorporating the 3D dynamical processes of GCMs with the flexibility and computational efficiency required to explore millions of models over a wide ranging parameter space. Potential future work could make the wind profiles more realistic with more complex physics (e.g. Seidel et al. 2020), or implement a varying temperature profile along the line of sight for each side of the atmosphere (Dobbs-Dixon & Blecic 2022; Nixon & Madhusudhan 2022), more akin to the thermal profiles from GCMs. Another potential improvement may be to use the likelihood method of Gibson et al. (2020). Such models are also ideally suited to simultaneous analyses combining high- and low-resolution observations (Gandhi et al. 2019), and will play a key role as we begin characterizing exoplanets with the next generation of facilities.

ACKNOWLEDGEMENTS

SG is grateful to Leiden Observatory at Leiden University for the award of the Oort Fellowship. This work was performed using the compute resources from the Academic Leiden Interdisciplinary Cluster Environment (ALICE) provided by Leiden University. We also utilize the Avon HPC cluster managed by the Scientific Computing Research Technology Platform (SCRTP) at the University of Warwick. IS acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement No. 694513. MB acknowledges support from from the UK Science and Technology Facilities Council (STFC) research grant ST/T000406/1. JPW sincerely acknowledges support from the Wolfson Harrison UK Research Council Physics Scholarship and the Science and Technology Facilities Council (STFC). LW thanks support provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51496.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. ABS acknowledges support from the Heising-Simons Foundation. We also thank the anonymous referee for a careful review of this paper.

DATA AVAILABILITY

The models underlying this paper will be shared on reasonable request to the corresponding author.

Footnotes

1

https://dace.unige.ch/dashboard/index.html.

Author notes

NHFP Sagan Fellow