Up, Up, and Away: Winds and Dynamical Structure as a Function of Altitude in the Ultra-Hot Jupiter WASP-76b
Aurora Y. Kesseli Hayley Beltz Emily Rauscher I. A. G. Snellen





Aurora Y. Kesseli
aurorak@ipac.caltech.edu
0000-0002-3239-5989]Aurora Y. Kesseli
These authors contributed equally to this work.
IPAC, Mail Code 100-22, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA


0000-0002-6980-052X]Hayley Beltz
These authors contributed equally to this work.
Department of Astronomy, University of Maryland, College Park, MD 20742, USA


0000-0003-3963-9672]Emily Rauscher
Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

0000-0003-1624-3667]I.A.G. Snellen
Leiden Observatory, Leiden University, Postbus 9513, 2300 RA, Leiden, The Netherlands



§ ABSTRACT 

Due to the unprecedented signal strengths offered by the newest high-resolution spectrographs on 10-m class telescopes, exploring the 3D nature of exoplanets is possible with an unprecedented level of precision.
In this paper, we present a new technique to probe the vertical structure of exoplanetary winds and dynamics using ensembles of planet absorption lines of varying opacity, and apply it to the well-studied ultra-hot Jupiter WASP-76b. We then compare these results to state-of-the-art global circulation models (GCMs) with varying magnetic drag prescriptions. We find that the known asymmetric velocity shift in Fe I absorption during transit persists at all altitudes, and observe tentative trends for stronger blueshifts and more narrow line profiles deeper in the atmosphere.  
By comparing three different model prescriptions (a hydrodynamical model with no drag, a magnetic drag model, and a uniform drag model) we are able to rule out the uniform drag model due to inconsistencies with observed trends in the data. 
We find that the magnetic model is slightly favored over the the hydrodynamic model, and note that this 3-Gauss kinematic magnetohydrodynamical GCM is also favored when compared to low-resolution data.
Future generation high-resolution spectrographs on Extremely large telescopes (ELTs) will greatly increase signals and make methods like these possible with higher precision and for a wider range of objects. 
 











§ INTRODUCTION


Transmission spectroscopy has dominated the study of exoplanet atmospheres since the first atmospheric detection <cit.>, and has moved the study of exoplanets from bulk properties to in-depth studies of composition, chemical processes, and dynamics in exoplanet atmospheres. Transmission spectroscopy performed at high spectral resolution has proven to be a powerful technique for detecting a range of atomic and molecular species that are often difficult to distinguish at lower spectral resolution <cit.>. By taking advantage of the fine radial velocity spacing of the spectrographs, high resolution spectroscopy has also been a valuable tool to uncover winds and atmospheric dynamics <cit.>. 


Hot and ultra-hot Jupiters offer the best test beds for assessing both chemical models and dynamical global circulation models (GCMs) for exoplanetary atmospheres due to their large scale-heights and favorable star-to-planet radius ratios. Recent instrumentation advances from facilities such as JWST, and increased throughput and wavelength coverage of next generation high-resolution spectrographs like ESPRESSO <cit.>, MAROON-X <cit.>, KPF <cit.>, and IGRINS <cit.>, mean that more complex planetary models can be tested. GCMs are 3D numerical models that simulate the dynamics of planetary atmospheres. These models are computationally expensive, but can capture many of the complex 3D processes such as clouds <cit.>, scale height effects <cit.>, magnetic drag <cit.>, hydrogen dissociation <cit.>, and changing line contrasts due to a 3D temperature structure <cit.> that influence the atmosphere and resulting spectra. 


The publicly available ESPRESSO dataset consisting of two transits of WASP-76b represents the highest-signal-to-noise-ratio high-resolution transmission spectrum that is publicly available to date. This dataset has been used in many papers due to its unprecedented quality, including the novel detection of asymmetric Fe I absorption during transit <cit.>, and subsequent detection and analysis of asymmetries in many other species <cit.>. It was also used to detect single absorption lines for a range of metals <cit.>, and in retrievals to produce phase- and spatially-resolved  parameters (abundances, temperature-pressure profiles) for the first time <cit.>. Most recently, <cit.> also precisely constrained the abundances of eight different atomic species.  


In this paper, we use the same WASP-76 b data and present a novel method to probe 3D atmospheric dynamics by exploring how radial velocities and line shapes change as a function of altitude (vertically). As was shown by <cit.>, measurements of how wind speeds and net Doppler shifts vary in altitude may help differentiate between different physical conditions input into GCMs, such as magnetic field strength and drag parameters. We aim to use the added vertical dimensionality offered by this technique to put more robust constraints on the physical processes (e.g., drag, magnetism, winds) occurring within WASP-76 b. 

We then compare these results to output from our state-of-the-art GCMs. This paper focuses on the treatment of drag in GCM models as the drag prescription is the main process that sets the wind speeds, and in turn the main influencer on measured radial velocities. Here we use the RM-GCM[https://github.com/emily-rauscher/RM-GCM] which was the first GCM to be post-processed and directly compared to high resolution emission <cit.> and transmission <cit.> spectra. These works, among others <cit.> highlight the inherent “3Dness” present in high resolution spectroscopic observations. An extremely multidimensional aspect of our models is the inclusion of a spatially varying magnetic drag. This state-of-the-art active magnetic drag prescription allows for the strength of magnetic drag to be calculated based on local conditions. Most GCMs use a single “uniform” drag timescale to simulate the effects of magnetic drag <cit.>. This assumption is particularly problematic for ultra-hot Jupiters (UHJs), as it can imply unphysically large magnetic field strengths on the nightside of the planet <cit.>. Previous works have shown the inclusion of magnetic effects can alter predicted phase curves by reducing hotspot offsets and increasing amplitude due to an increased day-night contrast <cit.>. This is in agreement with observations of ultra-hot Jupiters, which often show smaller hotspot offsets than their hot Jupiter counterparts <cit.>. Additionally, the inclusion of magnetic drag can result in net Doppler shifts that can vary on the order of several kilometers per second in high resolution emission <cit.> and transmission <cit.> spectra. We conduct our analysis on a uniform drag as well as a drag-free model for comparison to the active drag case. 

We present the data in Section <ref>, and our cross correlation analysis in Section <ref>. The GCM models are then described in Section <ref>. We present the results of the cross correlation analysis on the WASP-76 b ESPRESSO data in Section <ref> and <ref>, and compare the observed trends in altitude present in the data to state-of-the-art GCM models in Section <ref>. We discuss how our results fit in with previous work and potential future directions for GCMs in Section <ref>. Finally, we summarize our conclusions in Section <ref>. 



§ OBSERVATIONS AND DATA ANALYSIS
 



We downloaded the 1D blaze-corrected, stitched spectra from the Data and Analysis Center for Exoplanets (DACE) database[<https://dace.unige.ch/dashboard/index.html>], which were reduced with version 1.3.2 of the ESPRESSO pipeline. More details on the observing strategy are given in <cit.>. Before any cleaning or analysis steps, we removed contamination by telluric H_2O and O_2 lines using molecfit <cit.>, which has proven adept at removing the shallow absorption lines that mainly permeate the optical portion of the spectrum <cit.>. 

We outline the general cleaning steps we performed here, but refer the reader to <cit.> for more details, as the data presented here are the same as those used in that study. In order to remove the contribution from the host star so that the much weaker signal from the exoplanet can be seen, we created a uniform spectral time series along both the wavelength and flux axes. We performed the cleaning steps on each of the two nights separately, and then combined the nights after cross correlation. Using the host star's known radial velocity induced by the planet <cit.> and the system velocity <cit.> we shifted all of the spectra to rest. We then performed a simple 5-σ clipping to remove any spurious pixels due to cosmic rays. In order to preserve the relative flux values between different parts of the spectrum while removing any broadband noise or slight deviations in the shape of the blaze functions across observations, we placed all the spectra on a “common" blaze in a method similar to <cit.>. Finally, we interpolated all the spectra onto a uniform wavelength grid. 



 §.§ Cross Correlation Analysis
 




At this point the cleaned and uniform spectral time series are usually cross correlated with a model atmosphere created using the planet's parameters and one or more molecular or atomic species. This process has proven extremely successful at detecting a range of atoms and molecules in exoplanet atmospheres <cit.>. However, as discussed in detail in <cit.>, information on basic absorption line properties such as line amplitude, full width at half maximum (FWHM), or continuum dispersion are lost. Alternatively, a binary mask <cit.> can be used, which simply consists of a list of absorption line positions for the atom or molecule of interest. The main choice that affects the binary mask is the number of lines used, as some species have hundreds or even thousands of absorption lines. By including a large number of weak absorption lines that are unlikely to penetrate above the continuum of the planet, the signal will be diluted. This mask is then applied to the spectrum and any pixels that fall within the mask are averaged together (the weight of every mask pixel is treated the same and given a value of one in our analysis). In this way, the FWHM or amplitude of the resulting cross correlation function (CCF) from the binary mask is the average FWHM or amplitude of however many absorption lines were included in the mask. Hence, the amplitude of the co-added 1D CCF in the planet's rest frame can be expressed as an average excess absorption from a number of lines in parts per million (ppm) or percent, and is similar to what is reported for single lines resolved by high-resolution spectroscopy. 

Using this binary mask CCF approach, we created multiple masks from the opacity function of the atomic species Fe I. By comparing behaviors between different binary mask CCFs of the same chemical species, we attempt to isolate effects that are due to altitude from effects that may be present when CCFs of different species are compared to each other (i.e., chemical effects). The opacities were downloaded from the DACE opacity database[<https://dace.unige.ch/opacityDatabase>] and created using the open source HELIOS-K opacity calculator <cit.> and the Kurucz atomic line lists <cit.>. We created separate masks for the most opaque lines, the moderately opaque lines, and the weak lines. The most opaque lines become optically thick high up in the atmosphere and they trace the processes occurring at high altitudes, while the resulting CCF from the mask created from the lower opacity lines trace deeper in the atmosphere. This method allows us to directly probe different atmospheric depths which we compare to 3D atmospheric models. 

The opacity model of Fe I, with the absorption lines separated into their three bins is shown in Figure <ref>. The lower boundary of the strong lines bin was chosen to approximately map to a pressure of 10^-5 bars, which corresponds to the top of the model atmosphere. As detailed in section <ref>, extending the model to higher altitudes is computationally difficult, and so this strong line bin does not have a model counterpart. For the other two bins, we chose the number of lines for each which led to an even split in signal in the resulting CCFs between the final ESPRESSO spectra and the mask (see SNR values in Table <ref>). To determine the lower bound of the weak lines bin, we looked for the point where adding more Fe lines into our CCF mask no longer increased the resulting signal to noise ratio (SNR), calculated by averaging the CCFs in planet's rest frame and comparing the peak value to noise far from the expected peak. By comparing the bottom and top panel of Figure <ref>, the lower bound of the weak line bin approximately maps to a pressure of 1 millibar, which is consistent with the cloud deck location from the retrievals in <cit.>. While the binary mask CCF approach is less sensitive to parameter choices such as temperature than a model CCF approach due to line shapes not being important, opacity is still temperature sensitive. We find that we need to adjust the values that are used for the bin cutoffs (e.g., bottom panel of Figure <ref>) when using opacities of different temperatures, but when the bins are chosen so that the heights of the resulting CCFs are consistent, the derived parameters (radial velocity, FWHM, etc.) are also consistent to within the quoted uncertainties. 



We cross correlated each night of cleaned spectra with the three different binary masks. While the velocity bins for both nights were the same, the phase coverage and spacing on each night was different. We therefore interpolated the CCF grids onto a uniform phase-space grid, and then co-added the two together, weighting the contribution of each by their SNRs and the number of spectra taken during transit. 










§ GCM MODELING






In this work, we post-processed three models of WASP-76b <cit.>. These are double-grey models using the RM-GCM <cit.>, with 65 vertical layers evenly spaced in log pressure, from 100 to 10^-5 bars, and a horizontal spectral resolution of T31, corresponding to roughly 4 degree spacing at the equator. In this paper, we focus on the effect of drag and differing drag prescriptions since the main observables we are comparing to are radial velocities and line shapes, which are measurements or winds and are dominated by the treatment of drag. The three models presented differ in their treatment of magnetic effects. The simplest case, the 0G/drag-free model, represents the hydrodynamics only case. The uniform drag model applies a single global drag timescale of 10^4 s, chosen to match the strong drag case from <cit.>. 
Our final model, the 3 G case,  applies an active drag prescription, also known as a “kinematic MHD" method, and was first applied to hot Jupiters in <cit.>. This method applies a drag on the winds in the east-west direction <cit.> and with a timescale calculated based on local conditions, using the following expression from <cit.>:

    τ_mag(B,ρ,T, ϕ) = 4 πρ η (ρ, T)/B^2 |sin(ϕ) |

where B is the chosen global magnetic field strength (in this case 3 G), ϕ is the latitude, ρ is the density, and η the magnetic resistivity, which is a strong function of temperature. The 3 G model is presented here as it most closely matched the Spitzer phase curve of the planet <cit.>. Compared to the drag-free models, previous work has shown that the consideration of active drag results in a variety of behaviors not seen in drag free or uniform drag models. The differences most relevant to this paper are as follows:

    
  * The active drag model produces a different dayside upper atmosphere flow pattern. The 0 G model and uniform drag model show day to night flow centered on the dayside, with the uniform drag model showing slower wind speeds. Our active model shows flow moving mostly in the north-south direction, up and over the poles <cit.>.
    
  * Increasing the strength of our magnetic field results in a decrease of the hotspot offset (until it rests at the substellar point) and increase the day-night temperature contrast <cit.>.
    
  * Active drag models will show different Doppler shift trends in high-resolution emission <cit.> and transmission <cit.> spectra, compared to 0 G and uniform models of the same planet.   


In Figure <ref> we show line of sight (LOS) velocities including the effects of both winds and rotation for the east and west limbs during ingress, midtransit, and egress. These are plotted as a function of altitude and thus the atmospheric extent is highly influenced by the temperature structure. The two black contours indicate the pressure levels probed by the middle (outer contour, roughly 10 μbar) and weak lines (inner contour, roughly 1 mbar), respectively. The strongest lines probe pressures not captured by the GCM and are not shown. These differences in velocity structure—between types of drag as well as pressure levels—influence the net Doppler shifts found in the postprocessed spectra (see Section <ref>).  Although there are also varying temperature structures between the models, the winds will dominate the determination of the net Doppler shift. 




We then post-processed the two GCMs to generate high resolution transmission spectra at R∼ 450,000. Using ray-striking radiative transfer <cit.>, we take into account 3D effects such as temperature structure and Doppler shifts from winds and rotations. To generate Fe opacity tables, we assumed solar-abundance equilibrium models <cit.> with line lists from <cit.>. 


As the data are cleaned and processed in many steps, which masked or removed some fraction of pixels, we injected the model transmission spectra into the data at the same K_p and a systemic velocity of 100 km s^-1 (leading to an offset of 101 km s^-1 from the true planet's signal) in order to perform the same cleaning steps on the models and therefore accurately compare the models to the data. The model transmission spectra were generated for five phases: -0.04, -0.02, 0.0, 0.02, 0.04. For each model transmission spectrum, we interpolated the model onto the same wavelength grid as the data, and then follow convention and simply multiply the data by the interpolated model transmission spectrum for the corresponding phase. We then performed the exact same cleaning steps on the data that have been inject with the model as we did for the data that was not injected with the model <cit.>. 




§ RESULTS


We begin by discussing the trends we observe in Fe I in the atmosphere of WASP-76b, as Fe I shows the strongest signal. We then move on to comparing the Fe I trends to trends we detect in other atoms. Finally, we compare the trends we see in the data, with the trends we find from the injected GCM models. 



 §.§ Observed Fe I Trends with Altitude







We created three separate CCF grids from the separate strong-line, middle-line, and weak-line binary masks applied to the combined two nights of ESPRESSO data. Due to the slow rotation of the host star (v sin i=1.48 km s^-1; ) and misaligned orbit of the planet (λ=61.28 ^∘; ), contamination from Doppler shadow due to the Rossiter-McLaughlin effect is confined to a small region around ±10 km s^-1 in the rest frame of the star. Due to known difficulties in correcting for the Rossiter-McLaughlin effect <cit.>, we simply choose to mask out this region of the CCF grids. We used the planet's well-known velocity semi-amplitude (K_p) of 196.5±1 km s^-1 from <cit.>, and shifted all of the cross correlation functions to rest. These results are shown in Figure <ref>, and the presence of Fe I lines in the planet's atmosphere in all three opacity bins is clearly detected as white residuals in the planet's rest frame (along 0 km s^-1). As expected, the CCF grid created using the strong-line binary mask shows the largest residual amplitudes, meaning that these lines probe higher in the atmosphere. The middle and weak-lined masks show decreasing values in the residual amplitude as they probe deeper in the atmosphere. 

The residual amplitudes in Figure <ref> are related to the planetary radius or number of scale heights above the continuum where the absorption occurs <cit.>. We converted the residual amplitudes to planetary radii using the following equation: 


    R_λ = √(1 + h/δ) R_p ,


where R_λ is the effective planetary radius above the continuum where the average absorption occurs, h is the measured residual amplitude, and δ is the white light transit depth of the planet. We also convert this radius to scale heights byassuming a scale height of 1500 km, as reported for the dayside in <cit.>. 



We combined each different CCF grid in time by co-adding all of the phases within transit to get three 1D CCFs in the planet's rest frame (Figure <ref>). Consistent with Figure <ref>, the CCF created from the mask that used the weak Fe lines only extends to about 1.01 R_p, whereas the strong line CCF extends to 1.06 R_p. 
In the absence of any atmospheric dynamics and planetary rotation, we would expect the full phase-combined 1D CCFs to peak in the planet's rest frame, close to 0 km s^-1, but Figure <ref> shows that at all altitudes the planet's absorption is significantly blueshifted. As the signals are noisy and we do not know the underlying functional form of the CCFs, we simply fit Gaussians to each 1D CCF in the region around the observed signal (±80 km s^-1) using the python curve fitting code, , to determine the average radial velocity shift and an uncertainty on that radial velocity for each bin (see Appendix <ref> for more information on the reliability of these values and uncertainties). We also calculated SNRs by taking the peak value of the Gaussian fit and dividing by the standard deviation far from the peak (outside of ±80 km s^-1). The measured radial velocities for each bin are plotted in Figure <ref> and the SNR of each CCF is reported in the caption of Figure <ref>. These measured radial velocities for each CCF show a significant trend where the Fe lines in the lower atmosphere appear to be more blueshifted, and as we look higher in the atmosphere the lines become less blueshifted.












WASP-76 b shows a strong asymmetry between the absorption depth and the radial velocity shift measured in the first half of the transit versus the second half of the transit <cit.>, which can be seen clearly in Figure <ref>. This blueshift evolution over the course of the transit has been explained with GCMs <cit.> as being due to a combination of day-to-night winds and the hotter east limb (which is more blueshifted from tidally locked rotation) coming into view at the end of the transit. To try to isolate changes that are due to the vertical structure and not the change in viewing angle,  we also split the full transit into two phase ranges that separated the beginning and end of the transit and co-added these two parts of the CCF grid separately (Figure <ref>). When we examine Figure <ref>, we see that the significant trend of blueshift decreasing with altitude observed in Figure <ref> is no longer clearly observed. The majority of the perceived increasing blueshift is due to the fact that at lower altitudes the signal from the beginning of the transit (ϕ_-) is drastically diminished and so the combined CCF from the whole transit is influenced more by the ϕ_+ phases which are more blueshifted. 

We report the measured radial velocity shifts for the three altitude bins, separated by phase, in Table <ref>. We find highly consistent values and trends for the measured radial velocities at the beginning and end of the transit when we compare our values to those from previous studies of <cit.> and <cit.> even though we use a completely different method to either paper, demonstrating the reliability of the measured velocities and their lack of dependence on reduction method (see also Appendix <ref> for more discussion). Looking at how the measured radial velocities change with altitude, we find that the radial velocities at all altitudes are consistent within their uncertainties for the beginning of the transit (ϕ_-). For the ϕ_+ bin, we find that the radial velocities show a trend for stronger wind speeds deeper in the atmosphere (i.e. larger blueshifts are measured for the opacity bin that traces the deepest point in the atmosphere). In order to determine the significance of this trend we performed a linear regression using , which was chosen due to its handling of uncertainties in regressions, and find a low-significance trend with a p-value of 0.02 or 2.3-σ. Only observing a trend during the end of the transit (ϕ_+) can be explained given that recent GCM modeling <cit.> and phase-resolved retrieval results for WASP-76 b <cit.> suggest that during the second half of the transit the evening side completely dominates and the morning side of the atmosphere cannot be resolved, while during the beginning of the transit both the morning and evening side of the atmosphere can be observed. Therefore, the ϕ_+ CCFs originate from a single hemisphere, making it easier to observe any trends within the data. The three bins trace altitudes of ∼1.02 to 1.1 R_p, or 1 to 4 scale heights, meaning that at most we measure a 1.15 km s^-1 change in the wind between these altitudes. 



We also report the FWHM of the CCFs and the ratio of the heights between the phase bin for the beginning of the transit and the phase bin for the end of the transit in Table <ref>. These values and their uncertainties were also determined by fitting Gaussians to the CCFs shown in Figure <ref>. The CCFs for the phase bin encompassing the beginning of the transit are on average 6 km s^-1 wider than for the phase bin at the end of the transit. This is expected given the results from the GCM modeling and retrievals mentioned above which conclude that the Fe signal at the beginning of the transit originates from a combination of the morning and evening terminators, while the signal at the end of the transit is dominated by the evening side alone. Again, we only see a trend for the ϕ_+ phases, and measure an increasing FWHM as we probe higher altitudes. As previously mentioned, the majority of the perceived blueshift evolution in altitude shown in Figure <ref> is caused by the CCF signal in the beginning of the transit becoming weaker for the CCF mask that contained the least opaque lines (weak lines), and is reflected in our reported height ratio in Table <ref>, which decreases at lower altitudes. For the lowest altitude bin we find that during the beginning of the transit the Fe lines become optically thick much deeper in the atmosphere than during the end of the transit when the hotter evening side of the atmosphere is in view. 

Finally, motivated by recent modeling work that discussed how asymmetric absorption would manifest in a K_p vs. V_sys diagram <cit.> and observations of different offsets for different species in the same K_p vs. V_sys diagram <cit.>, we created separate K_p vs. V_sys diagrams for the 3 altitude bins and compare the relative offsets from the planet's known position in Appendix <ref>. We do not find any clear trend in K_p or V_sys offset for the three bins. 













 §.§ Observed Trends in Altitude from other Atoms and Ions





Fe I shows the strongest signal of any of the detected atoms or ions in WASP-76 b <cit.>, and so we focus our findings on Fe I. However, we also checked whether any of the trends seen in Fe I were also present in the absorption signals from other atoms (Cr I, V I, and Mn I) that showed significant absorption in <cit.>. For each atom we created two masks (instead of 3 as was done for Fe I) since the signals were not as strong and using more bins would result in less robust detections. One mask contained the most opaque lines while the second mask contained lines that were less opaque. We split the bins at a point where the two masks led to SNRs in the resulting CCFs that were roughly equal. The binary mask CCFs for the two opacity bins and three different atoms, co-added in the planet's rest frame are shown in Figure <ref>. 

We find the same trends in all three other species that we see in the analogous plot for Fe I (Figure <ref>), and in making the same K_p vs. V_sys plots for these species as we did for Fe I, we do not find any clear trends with altitude either (Appendix <ref>), reinforcing both conclusions drawn from Fe alone. Figure <ref> demonstrates that for these three other species the CCF made using the mask of weaker lines shows a larger blueshift than the CCF made using the mask of stronger lines. For Cr I and V I, we measure a significant (>3-σ) difference between the radial velocities of our two altitude bins. Mn I has the lowest SNR of the three trace species, and so even though we find the same blueshifting pattern, the difference between the two measured velocities is less than 2-σ. The consistency of the observed pattern in radial velocity with altitude for multiple species gives credence to the trend we observe with Fe I. It also demonstrates that the physical mechanism causing the signal to be weaker at lower altitudes during the beginning of the transit (decreasing height ratio trend) is not unique to Fe I, and therefore likely is a process that acts on a global scale (e.g., temperature or wind structure) and not on a single species (e.g., ionization, condensation). 



 §.§ Comparison with GCM models







We next aimed to use GCMs to help us interpret these trends and better understand the physical processes occurring in the atmosphere of WASP-76 b. To that end, we compared the measured radial velocities and line shapes of the injected models to the data to determine which drag prescription best fits the data. We injected the different models into the ESPRESSO data at a systemic velocity of +100 km/s so that the true signal and the injected signal would not overlap. We tested injecting the signal at different velocities and found that any differences were within the uncertainties. 

We found that all three models underestimate the strength of the Fe signal at all altitudes (see Figure <ref>), but that the mismatch with the bin that contains the strong lines was the worst (the models only reach a radius of 1.02 R_p while the data reach 1.05 R_p). The top of the atmosphere for the models occurs at 10^-5 bars, and so this significant underestimation from the models is likely due to the strong Fe lines absorbing at altitudes above 10^-5 bars in the planet's atmosphere. The top boundary of the GCM was chosen to be 10^-5 bars, as the the lower atmospheric densities result in numerical instabilities due to extremely short timescales. This upper boundary was also enforced in the post-processing, which is why the CCFs of the GCMs are much weaker for the strong lines. Additionally, non-LTE effects and photoionization — which are not modeled by the GCM — play stronger roles in the extended atmosphere <cit.>. Therefore, we focus the rest of our comparisons on the bins containing the middle and weak Fe lines, as these probe pressure levels covered by our GCMs. Even at lower altitudes, the models still underestimate the signal strength, which may be due to non-solar abundances of Fe – an assumption made in the post-processing routine – or differences in temperature structure compared to the true planet. 


Since our modeled spectra do not create a signal that is as strong as the data, and therefore have even lower SNRs (∼2.5-3.5 for the models and >5 for the data), in order to better compare the model to the data we inject the model into the data using standard injection and recovery methods <cit.> at 20 times the expected strength to create `zero noise model' CCFs. We treat these CCFs in the same way as the data and split them into different phase bins, fitting Gaussians to each CCF in order to measure the same CCF parameters for each model as we measured for the data. These newly measured CCF parameters are given in Table <ref>. More complex structure (e.g., double peaks) is now apparent in the CCFs due to the drastically increased SNRs, and demonstrate how this binary mask CCF method preserves line shapes and therefore information about atmospheric wind structures. By comparing Figure <ref>, we can see how differing wind patterns manifest in the CCFs. Despite the more complex structure, we still fit a single Gaussian to the models in order to treat the models and data in a consistent manner and derive parameters that can be compared in a one-to-one manner. The `zero noise' model CCFs are plotted along with the real data CCFs in Figure <ref>. We also plot all of the summarizing quantities from the two Tables in Figure <ref> to enable an easier visual comparison between the models and data. Each model shows unique features that we will attempt to match to the data. 



By comparing the CCF parameters for all the models, we see that all of the GCM models reproduce the observed pattern of more blueshifted CCFs in the second half of the transit, which naturally arises from the combination of tidally locked planet rotation and the hotter eastern limb coming into view during the later part of the transit (see Figure  <ref>). As the transit progresses the blueshifted contours become more dominant and are surrounded by regions of higher temperatures. 

Again looking at Table <ref> and Figure <ref>, we see that all of the models follow the clear trend in the data that the FWHMs at both phase bins increase with altitude (larger FWHMs for the middle Fe line bin than the weak line bin). We also see that all the models follow the tentative trend that we found in the data where for the ϕ_+ phases, the degree of blueshifting increases as we probe deeper in the atmosphere (weaker Fe I lines). Examining Figure <ref>, the 3D wind structures at the pressure contours vary significantly, which contribute to the differing Doppler shifts found at these two pressures. This effect is a result of the complex, multidimensional interplay between winds, temperature, atmospheric circulation. The equatorial jet is more pronounced at all longitudes at the deeper pressure level probed by the weak lines resulting in a narrower FWHM and a stronger average blueshift at post transit phases. At the higher level probed by the middle lines the wind structure is far more spatially variable, resulting in a wider FWHM and a smaller average blueshift at post transit phases (even though the winds speeds are typically larger at this altitude).   







We find that both the 0G and 3G models show significantly larger FWHMs for the first half of the transit than the second half of the transit, as is seen in the data. Conversely, the uniform model shows larger FWHMs for the second half of the transit, which is disfavored by the data. The uniform drag model also shows a stronger signal at the beginning of the transit and a weaker signal at the end of the transit (represented in Table <ref> by a height ratio greater than one). 
Both the 0G and 3G models reproduce the observed trend in the data of a weaker signal at the beginning of the transit. 
The uniform drag model applies the same drag at every spatial point on the planet, which causes less asymmetry than the other two models, and therefore a less dominant evening side of the planet. The signal is therefore weaker and has a wider FWHM due to contribution from both hemispheres. This is evidence that the uniform drag model does not accurately reproduce the observations and that a more complex treatment of drag is necessary. 








Between the 0G and 3G models, we find some instances where the 0G models seems to be preferred and some where the 3G model seems to be preferred. In the data, we find a decreasing height ratio as deeper pressures are probed, which means that the signal from the ϕ_+ bin become stronger compared to the ϕ_- bin as one looks deeper in the atmosphere. The 3G model shows a decreasing height ratio, like the data, while the 0G model has an almost constant height ratio at the two altitudes. This trend is likely influenced by multiple concurrent physical effects. In our GCMs, the 3G model has the largest day-night temperature contrast, resulting in the largest scale height difference between egress and ingress. The model also shows a jet that is more localized to the lower atmosphere than the 0G case, potentially leading to this large scale height difference being more localized to the lower atmosphere. Other processes, such as clouds or hydrogen dissociation (though not included in these models) could also likely influence this trend, and are discussion more in Section <ref>. 
On the other hand, the 0G model matches the measured RV shifts the best, in that it exhibits the strongest blueshift for the ϕ_+ phases, and the difference between the RV shift between the two phase bins is largest. This is because the 0G model has no additional drag, resulting in a stronger super-rotational jet. It is important to note that no GCM models have been able to reproduce the magnitude of the measured RV shifts found on WASP-76 b, and so there is likely some missing physics in the GCMs. 
Even though all of the models underpredict the observed blueshifts, both the 3 G and 0 G models follow the trend of increasing blueshifts as the transit progresses from ϕ_- to ϕ_+, and increasing blueshifts deeper in the atmosphere.





§ DISCUSSION



From our comparisons, we are able to show that applying drag uniformly to the atmosphere is not sufficient to fit the data and we are able to reject the uniform drag model. Previous comparisons to high resolution data found that models with weak drag were required to fit the asymmetric shape of the Fe CCF <cit.>. Conversely, low resolution Spitzer phase curves showed little phase offset, and therefore required models with significant drag, therefore setting up a tension between the two methods. As the previous GCM modeling results used a uniform drag prescription <cit.>, we suggest that the poor fits for the drag models could have been due to the spatially uniform application of drag, and find that the inclusion of spatially varying drag is necessary to best fit these high resolution data. From our analysis it is initially not obvious whether the 0 G or the 3 G model is the better fit, as both of the models are able to fit many of the trends. However, some of the first GCM modeling work for hot Jupiters <cit.> predicted that the trend of increasing blueshifts for deeper pressures—which is shown by the data—could be indicative of magnetic effects. Finally, it is notable that the 3 G model was also able to fit the low resolution Spitzer phase curves better than the 0 G model <cit.> meaning that our single kinematic MHD model fits both high resolution and low resolution data, which is a significant feat. As we have not done a full parameter sweep of varying field strengths or magnetic topologies, we are not claiming that the magnetic field strength of this planet is 3G, but take the agreement with both phase curves and high-resolution data as potential indication of magnetic effects shaping the environment.  


All of the models underpredict the wind speeds, and previous GCMs of this planet have been unable to match the magnitude of the phase-resolved Doppler shifts <cit.>. <cit.> suggested this discrepancy between observed radial velocities and those predicted by GCMs could be due to imprecision/uncertainties in the ephemeris and eccentricity, but recent precision measurements of the eccentricity and ephemeris of WASP-76 b from <cit.> have shown that orbital uncertainties are not to blame for inconsistencies with GCMs. It therefore seems that there is some missing physics that has not been taken into account for the GCMs that is needed to completely fit the WASP-76 b data. As of now it is unclear whether this problem of underpredicting planetary radial velocities is widespread, as there are few planets with this type of measurement, but initial wind measurements in a sample of 6 ultra-hot Jupiters found that WASP-76 b had the highest wind speed <cit.>.

Future work should explore in more depth why GCMs can underpredict wind speeds in HRS. Running GCMs with lower hyperdissipation and/or higher resolution may lead to higher model wind speeds. In addition, including physical processes such as hydrogen dissociation, clouds, and magnetic drag concurrently should be ran to better physically understand the atmosphere of this planet. Hydrogen dissociation would reduce the day-night temperature contrast <cit.> and the corresponding difference in scale height. Clouds, which may be present on one or both terminators, can increase the net blueshift, particularly during ingress <cit.>. Understanding how all of these processes act together and influence measured observables will shed light on which are the most important in future exoplanet modeling and the physical mechanisms occurring in planetary atmospheres. 








§ CONCLUSIONS




Using two ESPRESSO transits of WASP-76b, we explored how winds and dynamics in the planet's atmosphere changed as a function of vertical altitude, and compared these results to the output of state-of-the-art global circulation models (GCMs). To resolve the vertical structure of the atmosphere, we created binary masks containing Fe I lines in three different opacity bins. The bin with the Fe I lines that were the most opaque probes the highest layers of the atmosphere. The bins with less opaque Fe I lines probe lower in the atmosphere since light can travel through deeper layers of the atmosphere before it becomes optically thick. We then cross correlated these binary masks with the data and explored trends in the cross correlation functions. By cross correlating with a binary mask we are able to preserve information on the true line shapes, and therefore accurately extract the average height in the atmosphere where the signal arises and full width at half maximum (FWHM) of this signal. As the Fe I signal is famously also known to show an asymmetric velocity shift from the start of the transit to the end of the transit <cit.>, we also split up the cross correlation functions (CCFs) into two phase bins. 

We find that at all heights, the phase bin for the second half of the transit (ϕ_+) is more blueshifted than the first half of the transit (ϕ_-), and that the asymmetry first found by <cit.> and confirmed in <cit.> persists at all altitudes. We also see that at all altitudes the FWHM for the phase bin encompassing the beginning of the transit is significantly larger than the phase bin for the end of the transit. We take this as evidence that at the end of the transit the signal is dominated by only one side of the planet (the hotter evening side), whereas at the beginning of the transit the signal is a combination of the beginning and end of the transit, consistent with phase-resolved retrievals of WASP-76b <cit.>. When we focus on the ϕ_+ phases, which we assume come form a single hemisphere, we find that at higher altitudes the signal is both wider (has a larger FWHM) and less blueshifted. The blueshift evolution with altitude only shows a marginal trend (2.3-σ), but the FWHM trend is robust. Finally, we find that as we probe deeper in the atmosphere, the signal from the beginning of the transit is less significant (SNR=5.1 for the strong-line bin and SNR=2.84 for the weak-line bin) and peaks lower in the atmosphere (height ratios between the two half decrease from 0.76 for the strong-line bin to 0.44 for the weak-line bin). 

We next compared the data and the observed trends to GCMs with different prescriptions for magnetic drag by injecting the models into the data and performing the same analysis on the models.

All of the models show the trend of larger FWHMs at higher altitudes and a stronger blueshifted signal at lower altitudes for the ϕ_+ phase bin, as is seen in the data. We explain both of these trends as being due to the fact that at lower altitudes the models show smaller absolute velocities, but a more ordered velocity pattern due to the presence of a jet. At higher altitudes there are larger absolute velocities, but in a less ordered pattern, leading to broader CCFs with a smaller average blueshift. However, the models show key differences and we are able to rule out the uniform drag model due to the fact that the FWHMs for that model are larger at ϕ_+ phases and the signal is stronger during the first half of the transit. As both of these trends are clearly ruled out by the data, we are able to conclude that applying a spatially uniform drag prescription via a single drag coefficient cannot reproduce the data. The remaining two models both adequately reproduce the RV shifts, FWHMs, and height ratios seen in the data, but we argue that while the 3G model does not predict the degree of blueshifting as well as the 0G model, it better follows the trends seen in height ratio and FWHM. As all GCM models are unable to reproduce the degree of blueshifting <cit.>, and the 3G model also fit low resolution Spitzer phase curves better <cit.>, we take this as tentative evidence of magnetic field effects shaping the environment of WASP-76b. While more work is needed to explore if other field strengths, physical processes, or more complex radiative transfer routines  can provide a better fit to all of the data, it is promising that a single GCM model can adequately fit both the low-resolution Spitzer phase curve and the high-resolution optical transmission spectra, as previous modeling efforts found that to fit key features in the data they required significantly different drag prescriptions <cit.>. 







This study is a first step towards exploring exoplanets along new dimensions. As detection signal strengths increase, it becomes possible to better resolve the atmosphere by breaking the signal up. <cit.>, and later <cit.> and <cit.>, showed that by observing the shape and velocity evolution of planetary absorption lines during a transit, differences between the two terminators of a planet can be inferred. <cit.> performed retrievals that successfully characterized WASP-76b's atmosphere at 4 discreet locations across the planetary surface by breaking the transit up both temporally (probing different regions as the planet rotates during transit) and in velocity space (separating the oppositely rotating eastern and western hemisphere). In this work, we present a novel method to probe the vertical structure of the planet's atmosphere by taking advantage of the fact that absorption lines of different opacity become optically thick at different altitudes. With the current limits on signal, we are unable to draw extensive conclusions or distinguish between a 3-G and 0-G model, but with future facilities such as high resolution spectrographs on ELTs, this will be possible. 








We would like to thank Adam Langeveld for useful discussions that prompted inclusion of the K_p vs. V_sys diagrams. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. Many of the calculations in this paper made use of the Great Lakes High Performance Computing Cluster maintained by the University of Michigan. 


 













aasjournal





§ DETERMINING UNCERTAINTIES



As our conclusions are drawn from the measured radial velocities and their uncertainties, ensuring the accuracy of our values and uncertainties is vital for the interpretation of the data and in turn the robustness of our conclusions. In this section we try to calculate errorbars in a few different ways in order to ensure that they give consistent results and are truly representative of the data. We also compare our uncertainties to other publications of different data sets and discuss how the reduction and analysis steps affect the measured velocities to ensure that there are no hidden factors that are not represented in the errorbars. 

In the simplest method, we used  (a python package for curve fitting) to fit Gaussians to the CCFs. As the data are noisy and we do not know the underlying functional form of the data (e.g., Gaussian, double-peaked, etc.), some fitting or smoothing needs to be done to obtain measured values and Gaussians provides the best solution given the data quality and information.  reports uncertainties in the parameters it fits, one of which is the center of the Gaussian. These uncertainties come from the covariance matrices. Using the CCF from the strong line bin during the second half of the transit as an example (black line in top panel of Figure <ref>), we obtain a best-fit value of -9.565 ± 0.479 km s^-1 from . Using the same fitted Gaussians as above, but instead of simply assuming the uncertainties given by , we calculate the standard deviation of the fitted Gaussians and divide by the signal to noise ratio of the CCF <cit.>. This method gives an uncertainty of 0.348 km s^-1 for the same example as above. Finally, we also calculate uncertainties by recording the velocity shift of each individual CCF in time and then taking the standard deviation of that measurement and dividing by the square root of the number of measurements. This method was also used to derive alternative uncertainty values in <cit.>, and is similar to one of the methods used to derive radial velocity uncertainties in <cit.>. This final method assumes that the true radial velocity is same at every measurement (at least within a small phase bin), and therefore that any change in velocity is uncertainty in the measurement and not physical. This is not the case for at least the beginning of the transit, but it appears to be a reasonable assumption for the end of the transit (looking at Figure <ref>, phases from 0.02 to 0.04 have a vertical residual pattern). Using this approach, again for the same example as the above two, we are able to get 4 separate measurements of the radial velocity and find an uncertainty of 0.341 km s^-1. The fact that the three methods explored here return similar uncertainties lends credence to our uncertainty results. Throughout the paper, we quote the uncertainties as those from  as these seem to be the most conservative and are the simplest to derive. 

Next, we compare our uncertainty values to other recent papers to ensure that our values are on the same order as these.  <cit.> use a single ESPRESSO transit of WASP-121b and report a radial velocity uncertainty of 0.16 and 0.28 km s^-1 in the first and second half of the transit, respectively, when the CCF is split in two phase bins in a similar manner as we have done. Our uncertainties, which are on the order of 0.4 km s^-1, seem reasonable considering we have two transits, but we split our CCFs into 3 line-strength bins before we split by the two phase bins. 

Finally, there may be data reduction or analysis steps that cause small changes in the measured radial velocities. This uncertainty in the radial velocity would not be taken into account in the above calculations, and so we explore how some of our choices affect the measured radial velocities. We obtain consistent measurements for the radial velocity of the Fe I signal with <cit.> (ϕ_-=-3.89km s^-1, ϕ_+=-10.26km s^-1), and  <cit.>, as reported by <cit.> (averaging ϕ_- phases ∼-3.5km s^-1, and ϕ_+∼-10.5km s^-1), using a completely different method to either paper. In our own analysis, we do not perform any cleaning steps during our analysis process that significantly alters the data, unlike similar analyses in the NIR, which use PCA to clean the data. Indeed, <cit.> calculate uncertainties in their radial velocity measurements by using different numbers of PCA iterations and taking the resulting radial velocities as different ”measurements" to obtain errorbars, but as our analysis does not perform any drastic cleaning steps, an approach like this is not possible. We did test a few of our data processing choices (i.e., masking bad pixels, etc.) and found that the resulting measured radial velocity differences were smaller than the uncertainties we quote in the paper. The only choice that made a significant difference in the measured radial velocities was the choice of how many lines to include in each mask, but as long as the same choice is used during the GCM model analysis and throughout the analysis, the comparisons and conclusions will be accurate. We therefore do not think that the analysis itself imparts significant uncertainties into the measured radial velocities, and that both the values and uncertainties accurately represent the data. 

The data themselves are noisy because we have split the signal into so many different parts and due to the inherent increased noise when using the binary mask method as opposed to a traditional CCF method, but we note that every split CCF has a SNR>4.5, which is usually considered a reliable detection for high resolution observations. Because of this, we believe that the data are of high enough quality to measure reliable RVs. As data quality improves with the coming of ELTs, more robust measurements will be possible that rely on fewer underlying assumptions (e.g., Gaussian shapes).



§ TRENDS WITH ALTITUDE IN MEASURED K_P





Motivated by recent work that find differing offsets for different chemical species in the standard K_p vs V_sys diagram and suggested this difference could be due to species residing in different layers (altitudes) of a planet's atmosphere <cit.>, we tested to see whether our CCFs created with different altitude bins showed noticeable offsets in a K_p vs V_sys diagram. 

These diagrams are created in the standard method, by taking the 2D CCFs for each species and each altitude bin and shifting them to rest assuming different values of planet semi-amplitude (K_p) and different systemic velocities (V_sys). Once the 2D CCFs are shifted to rest, they are co-added in this rest frame and the SNR is calculated as described in Section <ref>. In order to compare the positions of the peaks in these diagrams, we have overlayed the SNR contours for the different altitude bins of each individual species in a single plot (Figure <ref>). By comparing the peaks for each bin (marked with crosses) and their associated SNR contours, we do not see a clear pattern emerge of differing K_p offsets at different altitudes. For the most part, the measured SNR peaks for a single species are consistent to within 1-σ, and the strong line peak is not always offset from the weak line peak in the same direction (i.e. the strong line bin peak for V is at a higher K_p value than the weak line bin, whereas for Cr and Mn it is the opposite). 

We see this lack of trend as potential indication that altitude is not the sole driver of differences in the measured peak K_p, and other processes, such as condensation or ionization, may lead to more significant differences in the measured K_p offset. It is important to note, however, that the species studied here only range in altitude from about 1.01 to 1.10 R_p (∼10^-3 - 10^-8 bar), and so species that primarily absorb in regions of the atmosphere outside of this range (e.g., H, He, etc.) may show significantly offset K_p values that are indicative of their altitudes.  






