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

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}$ K 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$\sigma$ over traditional 1D-retrievals, and thus demonstrate the power of such modelling for robust constraints with current and future facilities.


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 ionisation of atomic species. In recent years, a number of UHJs have been characterised, 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., Hoeĳmakers et al. 2019;Seidel et al. 2019;Bourrier et al. ★ E-mail: gandhi@strw.leidenuniv.nl † NHFP Sagan Fellow 2020; Pino et al. 2020;Cabot 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-characterised atmospheres ultra-hot Jupiters to date with a mass of 0.894 J , radius of 1.854 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;Žá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 3-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;Rauscher & Menou 2013;Dobbs-Dixon & Agol 2013;Mayne et al. 2014), as well as ultra-hot Jupiters (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 spatiallyhomogeneous 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;Mikal-Evans et al. 2019;Benneke 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;Pelletier et al. 2021;Line 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, highresolution 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 characterisation. Furthermore, GCMs have also been shown to provide better fits to high-resolution observations than 1D models which 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 which 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 which 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 multidimensional retrievals on high-resolution data. HyDRA-2D allows for a full exploration of the parameter space using Bayesian algorithms whilst capturing the dynamical effects which 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 ) 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 = −0.04 -− 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 et al. 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.

METHODS
In this section we discuss the setup for HyDRA-2D. We retrieve 2 nights of observations of the primary eclipse of WASP-76 b with the ESPRESSO spectrograph on the VLT . The observations span the optical, where Fe has prominent opacity. As these observations show strong variability in both the strength and

WASP-76
Morning Evening X Fe m , T 0.1mb m , 1 m , 2 m , P 1 m , P 2 m P 3 m , haze m , haze m , P cl m X Fe e , T 0.1mb e , 1 e , 2 e , P 1 e , P 2 e P 3 e , haze e , haze e , P cl e K p , V sys , V wind 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 Fe , the temperature profile is set through the six parameters mb , 1 , 2 , 1 , 2 and 3 using the parametrisation in Madhusudhan & Seager (2009), and the opacity and haze are parametrised with haze , haze and cl . We also retrieve Δ p , Δ sys and wind to constrain the wind profile of the atmosphere (see section 2.3). 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 parametrised approach where 3-dimensional 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 firstly 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.

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( ) value for each transition, where is the degeneracy and refers to the oscillator strength. For our work, we firstly compute the line strength in cm −1 /(atom cm −2 ) for each temperature on our grid, through where 0 = upper − lower , with the upper and lower state energies of the transition upper and lower given in cm −1 . The charge of an electron, , the speed of light , and the mass of the electron are all given in c.g.s units, and 2 = 1.4387769 cm K. In addition, the partition function is defined as for a state with degeneracy . 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 −1 ), Here, is defined as (Sharp & Burrows 2007) = 1 4 10, 000 where and are the van der Waals and natural broadening coefficients respectively and 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 ionised.
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 et al. 2013), where is the Boltzmann constant and 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), Thus the cross section as a function of frequency for each transition line is then computed as 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 −1 spacing, corresponding to a resolution of = 2.5 × 10 6 at 0.4 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, Fe , which we leave as a free parameter as listed in Table 1. The extinction coefficient (in m −1 ) is given by where Fe is the total cross section of Fe and is the total number -20 → 20 wind / kms −1 1 → 20 Table 1. Parameters and uniform prior ranges for HyDRA-2D retrieval. We retrieve separate Fe abundances, temperature profile, and opacity/haze decks for each of the morning and evening sides of the terminator. The Δ p , Δ sys and wind parameters are the same for both sides (see Figure 1). density of gas. The contribution to the optical depth d arising from an element of the atmosphere of length d is then Note that we include two free parameters for the Fe abundance, one for the morning and one for the evening side, as discussed below.

Separating the morning and evening limb
For our work we assume separate physical properties for the leading limb (or morning side) and trailing limb (or evening side) of the terminator (see Figure 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. 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.

Temperature profile
We parametrise 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 ultra-hot Jupiters (e.g., Fortney et al. 2008). We thus obtain the following parametrisation for as a function of pressure , and where we enforce 1 ≤ 3 . We additionally fix 0 = 10 −7 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 2 to vary between 0 ≤ 2 ≤ 3 . 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 − profile without the requirement that the temperature be monotonically decreasing with altitude. By enforcing continuity of the temperature between the 3 layers, we can reduce this set of equations to just 6 free parameters, a temperature value at any given pressure, 1 , 2 , 1 , 2 and 3 . For our retrieval we retrieve the temperature at 0.1 mbar, 0.1mb , 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 parametrisation of the − profile can be found in Madhusudhan & Seager (2009).

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-dependant super-Rayleigh haze. For our prescription we adopt a one sector approach with homogeneous opacity and haze parameters (e.g., Mac-Donald & Madhusudhan 2017; Welbanks & Madhusudhan 2021b) for each of the morning and evening limbs. We add this source of opacity into the overall extinction coefficient of the atmosphere, , as where is the total number density of gas, 0 = 5.31 × 10 −31 m 2 and 0 = 0.35 m. Hence, the free parameters are haze , haze and cl , whose prior ranges are given in Table 1. This therefore results in 6 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 (2021b).

Computing the unbroadened spectrum
We compute the unbroadened spectrum using numerical radiative transfer for each side of the atmosphere (e.g., Welbanks & Madhusudhan 2021b  Schematic of the broadening prescription used for HyDRA-2D. We firstly 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 blue-and red-shifted respectively. We then shift and broaden both sides by a day-night wind with Δ sys and 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. and Δ for the evening and morning sides, and is given by where p and s refer to the radius of the planet and star respectively, refers to the optical depth, and the impact parameter is given by . The subscripts and 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 p = 1.854 J at 10 −2 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.

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 Figure 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.

Planetary rotation
We firstly 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 blue-and red-shifted 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. 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 Δ sys = −10 km/s and wind = 8 km/s as shown in red. We also show an unbroadened spectrum shifted by -10 km/s in blue. rotation velocity, rot , the full rotation kernel is given by We only apply half of the kernel for each half of the terminator, as shown in Figure 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 Figure 3 for a small region near ∼0.4 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.

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 blue-shift 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 parametrise the wind with a velocity shift from the known systemic velocity, Δ sys , and a spread in its value given by the full-width half-maximum (FWHM) wind (see Figure 2). The Δ 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 which we may be sensitive to. However, in our case, the sys parameter is very precisely known at -1.17±0.02 km/s ) and thus we can assume that any deviation Δ 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, Δ p , as this often has a degeneracy with the Δ sys value. We therefore have three additional parameters in HyDRA-2D to account for the wind and its uncertainty, Δ p , Δ sys and wind as shown in Table 1. Figure 3 shows the effect of applying the wind broadening. We now see a net blue-shift in the spectrum as a result of the wind Δ sys , and an additional spread in the line profile as a result of 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 parametrisation can address the diversity of potential line profiles.

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. Figure 4 shows the setup to compute the limb darkening for each side of the atmosphere. We first compute the value of for a given phase of the observations , normalised by the stellar radius s , given by where refers to the semi-major axis of the orbit. We must then add and subtract p / 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 Figure 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 with the morning side given by m and the evening side given by e . We then compute the radial distance of the points from the centre of the star, The impact parameter is assumed to be 0.027 . We then compute the overall limb darkening function given the quadratic limb darkening coefficients 1 and 2 , Similarly for the morning side we obtain For our work we have taken 1 = 0.393 and 2 = 0.219 . We obtain the overall spectrum by simply scaling the model for each side by the value of and summing (see equation 25 below). 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.

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 summarise the steps here (see Figure 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).
We downloaded the pipeline reduced ESPRESSO spectra from the Data and Analysis Center for Exoplanets (DACE) database 1 . Before any other steps, we correct the spectra for telluric absorption lines due to H 2 O and O 2 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%. We then correct for velocity shifts due to the reflex motion of the star ( * ), the barycentric velocity ( ), and the system velocity ( ). Next, we perform a 5-clipping to remove any spurious flux spikes due to cosmic rays (top panel of Figure 5). To remove any broadband 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 −1 ; middle panel of Figure 5). Finally, we interpolate each spectrum onto 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 Figure 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, p . At a given phase ( ) the shift is given by where p = 196 km/s is the planet's rotation velocity. The Δ sys parameter is the day-night wind velocity and represents an overall shift in the planet's spectrum. We also include Δ p to account for uncertainties in the measured value of p and because there is often a degeneracy between the p and 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 where 0 is the wavelength of the unshifted spectrum. Therefore, the final combined spectrum for the terminator is given by where and refer to the stellar limb darkening factors, and Δ and Δ refer to the evening and morning side transit depths which 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 which 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 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-10-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.
Here, and are the data and model respectively for a model with spectral points, and represents the wavelength offset. This form of the likelihood accounts for the overall line shape as well as position; models which poorly match the observations or are scaled incorrectly will be penalised by the 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).

Retrieval setup
We retrieve the high-resolution observations of WASP-76 b, with two transits obtained on 02/09/2018 and 30/10/2018 with the ESPRESSO spectrograph on the VLT ). 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 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 which are consistent for both sides, Δ p , Δ sys and wind . Each spectral model is generated at a spectral resolution of R=500,000, assuming a planet radius of 1.854 J , planet mass of 0.894 J and stellar radius 1.756 . We use the Nested Sampling algorithm MultiNest (Feroz & Hobson 2008;Feroz et al. 2009;Buchner et al. 2014) for parameter estimation with HyDRA-2D, adopting 1000 live points.
We perform two separate retrievals, in the = −0.04 -− 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 between −0.0432 ≤ ≤ 0.0432 and the whole planet is in the transit between 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 Figure 6). This allows us to explore variations within the morning and evening side as the transit progresses and compare any differences.
To compare our results with HyDRA-2D we also perform a spatially-homogeneous retrieval across the combined = −0.04 -− 0.02 and = 0.02 -0.04 ranges, which we refer to as the combinedphase 1D-retrieval. We also perform two additional spatiallyhomogeneous retrievals for the = −0.04 -− 0.02 and = 0.02 -0.04 ranges separately (referred to as the separated-phase 1Dretrievals). Such spatially-homogeneous retrievals are well estab- Preference over combined-phase 1D 4.9 Preference over separated-phase 1D 3.6 - Table 2. Retrieved Fe abundances with HyDRA-2D for each side of the atmosphere and with the 1D retrievals for each phase range (see Figure 7). 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 = 0.02 − 0.04 range. . We assume uniform Fe abundance across the whole terminator for these 1D-retrievals, with a single pressuretemperature profile for the atmosphere and uniform opacity deck and haze. Hence for each retrieval we have 13 free parameters, including Δ p , Δ sys and 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 .

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 Figure 7 and the constraints on the temperature profile and opacity deck in Figure 8. These indicate that the median values and their errors vary significantly between the evening and morning, and between the retrievals across the = −0.04 -− 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 = −0.04 -− 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.

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( 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( Fe ) = −4.89 +0.43 −0.45 for the = −0.04 -−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 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 Figure 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 = −0.04 -− 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 Figure 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 Figure 9. We do this to investigate the potential biases which may arise by assuming spatial homogeneity. We calculate the average abundance of the terminator from HyDRA-2D through two methods. Firstly, 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 2021a). 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.
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 Figure 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 ( = 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.

Variation in Fe
The Fe is most strongly constrained from HyDRA-2D for the evening side of the atmosphere for the = 0.02 -0.04 range (see Figure 7). We see a tight well-constrained abundance of log( Fe ) = −4.03 +0.28 −0.31 , slightly above the solar value of log( 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 = −0.04 -− 0.02 range, with log( 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 which is being irradiated due to the geometry of the transit (see Figure 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 below.
Our morning side constraint for the = 0.02 -0.04 range has the lowest abundance and the weakest constraint, at log( 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 which 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 which 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. = −0.04 -− 0.02 evening side and  Figure 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.
= 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. Figure 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 et al. 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  . 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(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 ultra-hot Jupiter WASP-121 b (Gibson et al. 2022) in red, and H 2 O/H (blue) and CO/H (green) constraints for a range of exoplanets obtained with low-resolution and high-resolution observations respectively Guillot et al. 2022). 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), a UHJ with an equilibrium temperature of ∼2400 K, and the CO abundance constraints from other high resolution studies Pelletier et al. 2021;Line et al. 2021). These show good agreement with our result, as well as the previous work of  which found species other than oxygen at stellar-or super-stellar abundances.

Mass-metallicity relation
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 onto 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 et al. 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.

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 Figure 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 0.1mb = 2950 +111 −156 K. On the other hand, we obtain a lower median temperature with a wider uncertainty of 0.1mb = 2615 +266 −275 K for the morning side for the = 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 = 0.02 -0.04 range, where our error is more tightly constrained given the higher signal-to-noise.
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(Seidel et al. , 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 fit 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;) and shed light on the potential differences seen in the constrained temperature of the terminator.
We have performed retrievals assuming an H 2 /He rich and an H/He rich atmosphere, and find that the H/He rich case (where the H 2 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 2 and He rich atmosphere by a factor of almost two. Therefore, the scale height, given by = / , is almost double for an H/He rich atmosphere over an H 2 /He rich atmosphere.  Figure 11. Fraction of H 2 to (H+H 2 ) for typical photospheric pressures and temperatures for ultra-hot Jupiters. These were derived using the analytic prescription in Parmentier et al. (2018) (see also Gandhi et al. (2020a)). The lowest pressures with the highest temperatures have the greatest degree of H 2 dissociation, resulting in atmospheres which are H rich.
Modelling of circulation on UHJs has shown that the majority of the day side atmosphere is H 2 -dissociated at pressures ∼0.1 bar (Bell & Cowan 2018;Tan & Komacek 2019). In the upper atmosphere (∼ 10 −4 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 Figure 11. Therefore, our retrievals with an H 2 /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 2 /He rich retrievals were statistically disfavoured, and thus find evidence for an atmosphere that is dissociated. The recombination of H into H 2 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 scale height would be decreased by a factor of ∼2 over the evening in addition to any difference in temperature and Fe abundance.
We do observe a weak thermal inversion for the evening side in the = −0.04 -− 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 ultra-hot Jupiters (e.g., Lothringer et al. 2018;. We note however that the retrieved inversion may also be caused by the recombination of H 2 in the deep atmosphere. At high pressures, the H 2 is present at significant abundance as shown in Figure 11. This increases the mean molecular weight and thus reduces the scale height. Hence our retrieval with the H/He rich atmosphere may act to reduce the apparent scale height in the lower atmosphere by lowering the temperature at higher pressures, and thus appear as a thermal inversion. We tested our temperature profile parametrisation 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 1 and 2 , 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 scale height. We performed an additional retrieval which 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.

Opacity deck
The constraints for the opacity deck are given in Figure 8. The range in pressures constrained for the top of the cloud deck in Figure 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 which 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 which 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 which 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 et al. 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 −2 -10 −3 bar depending on the temperature, consistent with our opacity decks for the night side of log( cl /bar) = −2.08 +0.44 −0.53 (evening side, = −0.04 -− 0.02) and log( cl /bar) = −3.33 +0.80 −0.91 (morning side, = 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).

Abundance-opacity deck degeneracy
Our constrained Fe abundance shows a degeneracy with the pressure level of the opacity deck, as shown in Figure 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;. Recent work on characterising atmospheres with cloud/opacity decks at high resolution (Gandhi et al. 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 m.
We find that the Fe abundance retrieved for the morning side for = 0.02 -0.04 has a wider uncertainty due to the abundanceopacity 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 = 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.

Constraints on winds
The constraints on the wind parameters from HyDRA-2D are shown in Figure 13 and the corresponding constraints on the wind speed distributions are plotted in Figure 14. In our case the wind speed is given by Δ 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 sys is negligible. The Δ 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 Δ sys changes by at most ∼0.5 km/s. By using other recently published precise orbital values (Kokori et al. 2022;Ivshina & Winn 2022), we also find consistent results to within the uncertainties and ∼1 km/s. Therefore, any red-or blue-shift of the spectrum indicates the presence of a wind and cannot be explained by variations in orbital parameters.
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 Kesseli & Snellen 2021). This is a result of the day-night wind which 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 blue-shift (i.e. a wind travelling from 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 = 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 wind parameter, which measures the full-width half maximum of the spread in wind speed. We constrain a smaller wind for the last quarter of the transit, as shown in Figure 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 = 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 Δ sys . In addition, the spread in the wind would also be lower as the morning side contributes very little to the overall signal. On the other hand, the = −0.04 -− 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 red-shifted for the morning. In our constraints we do retrieve a wind profile that has a small red-shifted 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 Δ sys and spread in 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 = 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, 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 Δ p and Δ sys in Figure 13, as expected given the limited phase range of the observations which would act to break this. This degeneracy results in the two Δ p -Δ sys distributions for the first and last quarter, which tend towards each other at lower p values. However, the true p value is constrained to within ∼1 km/s , and so we believe the Δ sys offset between the first and last quarter of the transit is a result of the observations and not uncertainty in the p . The Δ p is consistent with expected value of 0 for the = 0.02 -0.04 range, but the = −0.04 -− 0.02 range does indicate a lower Δ 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 datasets 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 Δ p to compensate.

Comparison to GCMs
We compare the results from HyDRA-2D to general circulation models of WASP-76 b in Figure 15. These show the asymmetric weak-drag Monte Carlo radiative transfer model of Wardenier et al. (2021) (see "modification 1" in their Figures 12 and 13) at phases of = −0.031 and = 0.031, along with the best fit retrieval model for the = −0.04 -− 0.02 and = 0.02 -0.04 range for a small region near ∼0.366 m. For both of the best fit 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 datasets 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 shows that for the first quarter of the transit, the two models agree well. This is crucial as 3-dimensional 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-fit 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 fit model from HyDRA-2D for = 0.031 does not increase in strength over the = −0.04 -− 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% due to the changing scale heights 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 Figure 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 super-solar Fe abundance, other opti- cal 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.

Robustness tests
To verify the robustness of our results we performed numerous checks with our model. We carried out retrievals with a scale factor which rescales the model spectrum to account for uncertainties in the measured value of the radius and mass and any non-hydrostatic effects which may be present in the upper atmosphere. Such a parameter has often been included in previous high resolution retrievals (e.g., 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 scale height 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 = 0.02 -0.035 and = 0.025 -0.04 range, and = −0.035 -− 0.02 and = −0.04 -− 0.025 ranges. The top pressure that we modelled our atmosphere out to was also varied from 10 −7 bar to 10 −8 bar. As well as these tests, we performed retrievals with additional free parameters for the radius and log( ) of the planet. None of the retrievals showed any significant difference to any parameter, and our retrievals with log( ) constrained estimates consistent with the known values . 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 parametrising 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-0.4 m. For these tests we retrieved a higher temperature for the atmosphere and a negative velocity shift in Δ p . Given the very strong Fe opacity, these wavelengths probe regions of the atmosphere 10 −7 bar where nonlocal thermodynamic effects are likely to be considerable, and our assumptions of ideal gas and hydrostatic equilibrium break down. 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 m and 0.4-0.7 m range as well as our fiducial 0.4-0.8 m. These showed no change in the retrieved constraints for any parameter, despite some telluric absorption from O 2 and H 2 O at wavelengths 0.7 m. This confirms that the spectral cleaning steps were able to effectively remove the tellurics in this wavelength range (see section 2.5).

CONCLUSIONS
We have explored variations in Fe abundance, temperature profile and opacity/haze deck in the atmosphere of the ultra-hot Jupiter 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 which 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 ( = 0.02 -0.04), with an abundance of log( 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 4 /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 separatedphase retrievals, we find that HyDRA-2D is preferred by 3.6 over a spatially-homogeneous approach for the = −0.04 -− 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 high signal-to-noise and there is significant asymmetry in the physical parameters. On the other hand, we find that our 1Dretrieval for the = 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 which 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 = −0.04 -− 0.02 and the evening side for = 0.02 -0.04. We also find that our spatiallyhomogeneous 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 = 0.02 -0.04 range, with the retrieved abundance, log( 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 = 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 ( = 0.02 -0.04, evening side) and 2615 +266 −275 K ( = 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 2 and He. This is because an undissociated atmosphere has a smaller scale height 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 2 dissociation in the atmosphere of WASP-76 b, consistent with predictions made for ultra-hot Jupiters (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 blue-shifted 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, wind , is also more spread for the = −0.04 -−0.02 range at 10.2 +1.0 −1.0 km/s than for the = 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 daynight wind to drive the overall wind speed to higher values for the = 0.02 -0.04 range. However, the stronger wind speed and lower 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 spatiallyresolved retrieval analyses of exoplanetary atmospheres. Key to these retrievals is incorporating the 3-dimensional 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 (Nixon & Madhusudhan 2022;Dobbs-Dixon & Blecic 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 , and will play a key role as we begin characterising exoplanets with the next generation of facilities.

DATA AVAILABILITY
The models underlying this article will be shared on reasonable request to the corresponding author.