No Umbrella Needed: Confronting the Hypothesis of Iron Rain on WASP-76b with Post-processed General Circulation Models

For our study, we make use of a line-of-sight ray-striking radiative transfer model that includes Doppler effects, inherited from Kempton & Rauscher (2012), Kempton et al. (2014), Rauscher & Kempton (2014), and Flowers et al. (2019). We begin with the GCM output, which we interpolate from a grid evenly spaced in log-pressure onto a grid evenly spaced in altitude—thus ensuring that our simulated rays propagate in straight lines. We further only consider the GCM output at pressures generally less than 1 bar (the cut is made in altitude, not pressure, so the bottommost pressure value is weakly latitude-/longitude-dependent), as the atmosphere is fully optically thick in transmission geometry at pressures greater than this value (Figure 3). We convert between pressure and altitude assuming hydrostatic equilibrium, accounting for the nonuniform mean molecular weight between GCM grid cells due to the thermal dissociation of H₂ in hotter regions of the atmosphere. Our calculation assumes a nominal atmosphere in cooler regions that is 83.6% H₂, 16% He, and 0.4% metal-bearing molecules (the latter having an average mean molecular weight of μ = 12) by number.

Zoom In Zoom Out Reset image size

Accounting for the 3D geometry of our GCM structures, we calculate the slant optical depth,

$$\begin{eqnarray}&&{\tau }_{\lambda }=\int {\kappa }_{\lambda }{ds},\end{eqnarray} \tag{ 1 }$$

along a given trajectory after determining the opacity κ_λ in each grid cell as a local function of temperature and pressure (see Section 2.2.2 for a discussion of our opacity sources), where ds is the line-of-sight path length through an individual grid cell. The wavelength at which κ_λ is evaluated is adjusted depending on the line-of-sight velocity of the grid cell v_LOS, as per

$$\begin{eqnarray}\begin{array}{rcl}\lambda & = & {\lambda }_{0}\left(1-\displaystyle \frac{{v}_{\mathrm{LOS}}}{c}\right)\\ {v}_{\mathrm{LOS}} & = & -u\sin (\phi )-v\cos (\phi )\sin (\theta )\\ & & -{\rm{\Omega }}({R}_{{\rm{P}}}+z)\sin (\phi )\cos (\theta )\end{array}\end{eqnarray} \tag{ 2 }$$

for rest wavelength λ₀, east–west velocity u, north–south velocity v, latitude θ, longitude ϕ, altitude z, speed of light c, and planetary angular rotation speed Ω. These Doppler shifts assume that the orbital velocity of the planet has already been entirely accounted for; we address potential consequences of this assumption with respect to WASP-76b in Section 5.2.1. As described in Kempton & Rauscher (2012), higher-order line-shifting effects such as gravitational redshifting and microturbulent broadening are expected to be negligible.

After calculating the above τ_λ , we can calculate the intensity I_λ along each line of sight, assuming absorption-only radiative transfer:

$$\begin{eqnarray}&&{I}_{\lambda }={I}_{\lambda ,0}{e}^{-{\tau }_{\lambda }}\end{eqnarray} \tag{ 3 }$$

for stellar intensity incident upon a grid cell, I_λ,0. We assume that the stellar spectrum follows a blackbody distribution for T = 6329 K, the effective temperature of WASP-76 (Ehrenreich et al. 2020). The total flux transmitted through the atmosphere is arrived at by integrating I_λ over the sky-projected solid angle of each respective grid cell through which a beam emerges. The wavelength-dependent transit depth D_λ is then given by

$$\begin{eqnarray}&&{D}_{\lambda }=\displaystyle \frac{{F}_{\mathrm{out}}-{F}_{\mathrm{in}}}{{F}_{\mathrm{out}}},\end{eqnarray} \tag{ 4 }$$

where F_out is the out-of-transit flux (i.e., the stellar flux) and F_in is the in-transit flux, accounting for both the blockage of stellar light by the optically thick core of the planet and the transmitted flux through the atmosphere.

Our spectral simulations are computed over the 379–789 nm range to coincide with the ESPRESSO data from Ehrenreich et al. (2020), with a total of 298,766 individual wavelength points corresponding to a resolution of R ≈ 4 × 10⁵. After computation, our spectra are convolved to the native ESPRESSO resolution in the singleHR mode (R ≈ 138,000; Pepe et al. 2013). Because of the very large file sizes and RAM requirements associated with the substantial number of wavelength points, we generate the spectra in discrete wavelength "chunks" on a high-performance computing cluster and then stitch the entire spectrum together once the full calculation is complete.

Because WASP-76b is an ultrahot Jupiter observed at optical wavelengths, we compile a set of opacities that is appropriate for these conditions. The atomic species that we utilize are Fe i, Fe ii, Mg i, Mn i, Na i, Ca ii, Ti, Li i, and K i; and our molecular species are TiO, VO, OH, and H₂O. Rayleigh scattering from H₂, H, and He, as well as H⁻ bound–free and free–free absorption are included as continuum opacity sources. Finally, collision-induced opacity from collisional pairs of H₂–H₂, H₂–H, H₂–He, and H–He is considered. This list of opacity sources includes all of the species with reported detections in WASP-76b optical spectra from Ehrenreich et al. (2020), Tabernero et al. (2021), and Fu et al. (2021). We do not include ZrO, which was searched for in WASP-76b by Tabernero et al. (2021) but resulted in a nondetection. In modeling WASP-76b, we generate two sets of transmission spectra, one with Fe (i and ii) and continuum opacity sources only, and another with the full set of molecular, atomic, and continuum absorbers (see Figure 4 and Section 4.1.1).

Zoom In Zoom Out Reset image size

All opacities are generated on a temperature–pressure grid spanning 500–5000 K and 1 × 10⁻⁶–1 × 10³ bar, covering the range of the WASP-76b GCM output. We initially calculate opacities on an evenly spaced wavenumber grid, from a wavenumber of 0 to 30,000 cm⁻¹ with a spacing of 0.01 cm⁻¹. These results are then downsampled with the k-table feature of HELIOS (Malik et al. 2017).

We use the newest available open-source line lists in constructing our molecular opacities: VO (VOMYT; McKemmish et al. 2016), TiO (TOTO; McKemmish et al. 2019), OH (MoLLIST; Tennyson et al. 2016), and H₂O (POKAZATEL; Polyansky et al. 2018). Atomic opacities are drawn from the Kurucz database (Kurucz & Bell 1995). In selecting our opacities, we make use of the most abundant isotopologue for each species. Our atomic sources do not include pressure broadening—only Doppler (thermal) and natural broadening. ⁸ From these line lists, we calculate the relevant opacities with the GPU-enabled HELIOS-K code (Grimm & Heng 2015).

Opacities from individual species are combined by a weighting of their respective mixing ratios, under the assumption that the atmosphere of WASP-76b resides in a state of chemical equilibrium. To determine the chemical equilibrium mixing ratios, we employ the FastChem model (Stock et al. 2018), which accounts for gas-phase chemistry only. In Section 2.2.4, we describe how we treat condensation for the purposes of our modeling.

As in Flowers et al. (2019), to compare our model transmission spectra directly against the Ehrenreich et al. (2020) results, we must calculate the transmission spectra as a function of orbital phase throughout the duration of transit. To account for orbital phase dependencies, we apply the following procedure:

• 1.  
  Account for phase-dependent backlighting of the planet (i.e., stellar limb-darkening effects). At different points of its transit, the planet will occult regions of its host star of varying brightness. Furthermore, at a fixed orbital phase, different regions of the planet's limb will be backlit by varying intensities of stellar light. Similar to Flowers et al. (2019), we calculate the normalized stellar intensity at the center of each cell of the 2D projected planetary grid produced by our GCM at each modeled orbital phase of the planet. We use the quadratic limb-darkening coefficients reported by Ehrenreich et al. (2020) to establish the stellar center-to-limb intensity profile, and we take into account the 89623 orbital inclination of WASP-76b (Ehrenreich et al. 2020) to determine where the planet resides on the stellar disk as a function of its orbital phase. We make the assumption of constant impact parameter b over the course of transit. ⁹ This procedure allows us to calculate a backlighting factor f, which ranges from 0 to 1, effectively replacing the constant I_λ,0 from Equation (3) with a variable ${I}_{\lambda ,0}\times f(\theta ^{\prime} ,z,\varphi )$, for a given orbital phase φ and 2D projected polar angle $\theta ^{\prime}$.
• 2.  
  Account for the decreasing of the continuum by interpolating a light curve produced by the batman code (Kreidberg 2015). Step 1 ensures that less light is transmitted through the planet's atmosphere than a uniform stellar disk would emit. Step 2 further enforces that the inner, optically thick core of the planet is simulated crossing a limb-darkened star, as opposed to a star of uniform brightness.
• 3.  
  Account for the planet's rotation over the course of transit. Because the planet is continually rotating as it travels across the face of its host star, we must transform the GCM coordinate system so that the correct observer-facing hemisphere is modeled at each instance during transit. For simplicity, we assume zero obliquity, which allows us to calculate the coordinate transform simply by assigning a linear offset to each planetary longitude; i.e., ϕ_rotated = ϕ + φ.

In choosing the phases at which to evaluate our models, we more densely sample ingress and egress because the transmission signal quickly varies in strength and velocity during these times. The sampling that we employ is

$$\begin{eqnarray}d\varphi =\left\{\begin{array}{ll}0\buildrel{\circ}\over{.} 416, & -15\buildrel{\circ}\over{.} 94\leqslant \varphi \leqslant -8\buildrel{\circ}\over{.} 25\\ 1\buildrel{\circ}\over{.} 39, & -8\buildrel{\circ}\over{.} 25\leqslant \varphi \leqslant 6\buildrel{\circ}\over{.} 0\\ 0\buildrel{\circ}\over{.} 416, & 6\buildrel{\circ}\over{.} 0\leqslant \varphi \leqslant 15\buildrel{\circ}\over{.} 94\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where ${\varphi }_{\max }=15\buildrel{\circ}\over{.} 94$ represents the phase at which the planet's trailing limb no longer occults the stellar disk. This procedure results in 50 total orbital phases being explicitly modeled with our radiative transfer post-processing code.

The large blueshifts observed by Ehrenreich et al. (2020) have been attributed to iron condensation on the nightside and cooler "morning" limb of WASP-76b. To confront this hypothesis within the framework of our transmission spectrum modeling, we account for clouds and condensation in the following ways.

To incorporate the effects of condensation out of the gas phase, we utilize the iron condensation curve for a solar composition mixture as computed in Mbarek & Kempton (2016). At each cell in our GCM-produced atmosphere, we interpolate the condensation curve at the cell's pressure; if the temperature in the cell is lower than the interpolated condensation temperature, then we remove iron from the total opacity calculation at that location by setting its opacity contribution to zero. This approach maximizes the effect of Fe condensation, as in reality homogeneous nucleation of Fe is inefficient (Gao & Powell 2021). This step effectively treats the physics considered in the Ehrenreich et al. (2020) toy model, which only accounted for gas-phase removal as the process responsible for the large net blueshifts in WASP-76b's transmission spectrum.

In our full-species modeling of WASP-76b's transmission spectrum, we allow other species to condense out of the atmosphere following a similar procedure. We condense Ca (i and ii) following the Ca₂SiO₄ condensation curve, Ti and TiO following CaTiO₃, K following KCl, Na following Na₂S, Mn following MnS, and Cr following Cr₂O₃. Our prescription of complete removal implies full rainout from the gas phase once the species-limiting condensate forms.

Condensation not only causes gas-phase depletion, but can also result in the formation of an optically thick cloud made up of liquid droplets or solid particles. Such a cloud layer would block the transmission of stellar light through the atmosphere in the location where the cloud forms, and it could therefore also be the main driver of the "missing" iron absorption in the receding limb of WASP-76b's atmosphere—a scenario that has not yet been explored in the literature for this planet. To account for such cloud layers, our models further contain a post-processing implementation of gray, optically thick clouds. These clouds, which are self-consistent with our radiative transfer (but not with our GCMs), are placed such that they span up to several (1–10) scale heights above the condensation curves of iron, forsterite (Mg₂SiO₄), or aluminum oxide (Al₂O₃).

Fe is selected because it is the species that is putatively condensing in WASP-76b's atmosphere; Mg₂SiO₄ is modeled because per the microphysics model of Gao et al. (2020), silicate clouds are chemically favored over iron ones, and they appear to be a defining feature in the transmission spectra of giant planets with comparable irradiation to WASP-76b; and Al₂O₃ is considered because it is the highest-temperature condensate predicted in hot Jupiter atmospheres (Mbarek & Kempton 2016).

Notably, our implemented clouds never fully reach the limiting-case thickness of 10 scale heights. Because our temperature–pressure profiles tend to cross condensation curves twice (Figure 2), clouds cannot exist in the uppermost regions of the atmosphere—the upper boundary of the cloud is (indirectly) set by the pressure at which the atmospheric inversion or isothermal region occurs. This effect can be seen in the bottom row of Figure 5, in which the "10 scale height" clouds do not extend much past the "1 scale height" clouds. Hence, while our clouds can exist up to 10 scale heights in thickness (provided they fulfill their criterion with respect to a given condensation curve), in practice they are limited by the temperature/pressure structure of the atmosphere.

Zoom In Zoom Out Reset image size

Similar to our condensation treatment, our cloud effect is modeled post hoc, in that the radiative effects of cloud formation and the resultant changes to opacities are not accounted for natively in the GCM (see Section 5.3 for a discussion of this limitation). Furthermore, our simple cloud model does not account for the wavelength-dependent scattering and/or absorption from cloud particles—a simplification that is justified because cloud opacities are typically smoothly varying with wavelength, in contrast with the sharply wavelength-dependent atomic and molecular features that are primarily responsible for the high-resolution transmission spectrum signal in cross-correlation. Additionally, cloud particle sizes are represented by broad distributions, and no clouds in a planetary atmosphere are characterized by a single particle size (e.g., Pont et al. 2013).

For our first set of results, we perform our analysis on our iron-only transmission spectra.

We find that controlling the drag timescale over our explored parameter range has the strongest first-order effect on Doppler shift over all phases (Figure 7), with the weak-drag models providing a much better match to the Ehrenreich et al. (2020) data than the strong-drag models. Decreasing the drag timescale in an atmosphere is equivalent to making drag processes more efficient. In the Helmholtz decomposition framework of Hammond & Lewis (2021), about half the day–night heat redistribution in a fiducial hot Jupiter atmosphere is driven by substellar–antistellar flow; shorter drag timescales would serve to slow the speed of this flow, reducing heat transport. The remainder of the day–night heat redistribution, per Hammond & Lewis (2021), is due to rotational flow—e.g., the superrotating equatorial jet. The speed of this jet would similarly suffer directly from a shorter drag timescale. Additionally, per Showman et al. (2013), stronger drag also reduces the pumping of eddy angular momentum onto the dayside equator because drag damps the propagation of planet-scale wave patterns, further decreasing the jet strength. In sum, increased drag reduces longitudinal asymmetries and slows wind speeds, both of which reduce net Doppler shifts in transmission spectra.

Zoom In Zoom Out Reset image size

As a result of the aforementioned effects, our high-drag models have more homogenized limbs, whereas low-drag models have limbs that are more heterogeneous, with a visible equatorial jet structure (Figures 3, 5, and 8). These features can be seen in our simulated Doppler shifts (Figure 7). The strong-drag case has Doppler shifts that are roughly symmetric about the center of transit, with only a slight (<2 km s⁻¹) blueshift at midtransit. Moreover, the maximum and minimum of the Doppler shift are both near ±4–5 km s⁻¹, which is close to the equatorial rotational velocity of WASP-76b (±5.3 km s⁻¹). These features are expected for a solely rotationally broadened profile—which would be predicted for the high-drag models, which have weaker winds. Additionally, a clear prediction from Figure 8 is that the strong-drag case is redshifted at ingress and blueshifted at egress—as seen in Figure 7.

Zoom In Zoom Out Reset image size

Similarly, the Fe abundance map (Figure 5) combined with the wind map (Figure 8) explains the behavior of the weak-drag case in Figure 7. Namely, the Fe abundance is asymmetrically distributed, preferentially on the eastern limb, where the strongest blueshifts for this model are. Therefore, the model produces Doppler shifts that are not evenly distributed around 0 km s⁻¹ across phase and are more strongly blueshifted, which is more consistent with the Ehrenreich et al. (2020) observations.

Including iron condensation removes the gas-phase iron in the cooler regions of the atmosphere (Figure 5). As can be inferred from our temperature–pressure structure (Figure 2), condensation is primarily confined to the nightside in all our GCMs. Especially in our weak-drag models, condensation is less prominent at the equator, whereas it is more prominent at the midlatitudes and poles. The effect of condensation tends to be marginal, dependent on cloud placement, and not necessarily linear, reflective of the complexities in the 3D thermal structure of the atmosphere.

We further find that our deep RCB models are a better fit to the data, generating stronger net blueshifts. A shallow RCB is expected to create a more homogeneous atmosphere, because convective motions quickly mix entropy across large spatial scales, reducing day–night temperature contrasts on deep isobars relative to the deep RCB case. Atmospheric circulation in (ultra)hot Jupiters is driven and supported by day–night pressure gradients (e.g., Showman & Guillot 2002), so the smaller day–night contrasts in the shallow RCB case cause winds to be slower. These slower winds in turn decrease east–west limb asymmetries as well as the strength of the blueshift signal (as shown in Figures 7 and 9).

We find that the presence of optically thick clouds, as opposed to gas-phase iron condensation, also strongly and consistently increases the measured blueshift (Figure 7). This is because clouds will preferentially form on the cooler, receding (eastern) limb, hence blocking its contribution to the transmission spectrum; these clouds would lessen the redshift contribution of the receding limb, thereby boosting the measured blueshift.

As seen in Figure 2, Al₂O₃ condenses out at hotter temperatures than Fe. Hence, the placement of clouds above the Al₂O₃ condensation curve as opposed to the Fe condensation curve results in a cloud deck situated in warmer regions of the atmosphere. The effect of this change at the blueshift level is highly nonlinear, as it depends somewhat strongly on the exact placement (and thickness) of the clouds. Notably, though, the optically thick clouds produce a "jump" in blueshift near a phase of −7° that the other models cannot (Figure 7)—a distinct feature of the Ehrenreich et al. (2020) data. While this jump does not extend to the exact observed magnitude, it appears that an Al₂O₃ cloud deck (or a similarly positioned cloud of arbitrary composition) provides the type of spatial inhomogeneity required to reproduce the Ehrenreich et al. (2020) data in broad strokes.

Interestingly, although the deep RCB/strong-drag GCM would have much more condensation in its atmosphere (Figure 2), the deep RCB/weak-drag GCM has much stronger blueshifts (Figure 9). It thus appears that the importance of the drag timescale and faster winds leading to more east–west asymmetry "wins out" over the existence of condensation alone.

Zoom In Zoom Out Reset image size

Our net center-of-transit Doppler shifts are summarized in Figure 9 for all of the model permutations explored. In total, we find that the following effects produce increased blueshifts in WASP-76b's atmosphere when iron is considered alone: weak drag, deep RCBs, and optically thick clouds composed of Al₂O₃, Mg₂SiO₄, or Fe. Gas-phase iron condensation increases net blueshifts in some cases, but (contrary to the Ehrenreich et al. 2020 interpretation) actually decreases net blueshifts in others.

By toggling our backlighting effect, we determine that accurate treatment of limb-darkening (see Section 2.2.3) has a median effect over the course of transit on the order of 0.1 km s⁻¹. Considering the mean error bar size on the Ehrenreich et al. (2020) blueshift data (Figure 7), it appears that although including this effect improves physical consistency, it does not have strong (or even currently detectable) effects at the Doppler shift level.

For our next set of model runs, we include all of the opacity sources described in Section 2.2.2 and cross-correlate these "full species" models against Doppler-off templates that include one chemical species at a time. We are able to identify peaks in the CCFs at most in-transit phases for Fe, Ti, Mn, Ca ii, and Cr, indicating that these species are present at an observable level in our modeled spectra.

Similar to the Fe-only models, we find that spectrally post-processed, unaltered GCMs cannot reproduce the Ehrenreich et al. (2020) piecewise blueshift trend without added modifications such as condensation or cloud formation (Figure 10, panels (a) and (b)). Most of the detected species tend to follow the trend of Fe, especially in the weak-drag case. The clear exception is Ca ii. It exhibits much weaker blueshifts than the other species in the weak-drag case, and is mostly redshifted over the course of transit in the strong-drag case. Interestingly, its Doppler shift remains mostly constant out to egress in the weak-drag case, whereas it experiences a sharper blueshift in egress in the strong-drag case. This trend implies that in the weak-drag case, the Ca ii signature is proportionally stronger in the approaching limb, so as the receding limb exits the stellar disk in egress, the Doppler signal remains largely unchanged until the approaching limb exits as well. In contrast, in the strong-drag case, it appears that the Ca ii signal is more evenly distributed on both limbs, so when the receding limb begins to exit the stellar disk, the total contribution of the Ca ii signal shifts to the approaching limb, rapidly increasing the measured blueshift over egress due to the rotational component.

Zoom In Zoom Out Reset image size

These characteristics can be physically motivated by considering the temperature structure of the atmosphere. The abundance of Ca ii in equilibrium chemistry is strongly dependent on thermal ionization, and hence temperature. Any temperature asymmetries, then, will be reflected in asymmetries in the abundance profile of Ca ii. Per Figure 3, our weak-drag case is much hotter on the approaching limb, because the lack of strong drag is conducive to the formation of a superrotating equatorial jet that advects heat from the planet's hotspot, causing a thermal offset in the direction of the approaching limb. Therefore, in the weak-drag case, the approaching limb will be more abundant in Ca ii than the receding limb. In the strong-drag case, however, the lack of a jet and hotspot offset imposes limb homogeneity, giving each limb relatively equal weight in the Ca ii signal. See Section 5.2.2 for a further discussion of Ca ii and its potential nonhydrostatic behavior in observations.

Adding clouds to our models (Figure 10, panel (c)) has the same effect as identified in the Fe-only models: stronger blueshifts and a trend that better matches the Ehrenreich et al. (2020) data for clouds with higher condensation temperatures.

Notably, the cloud thickness seems to strongly control the jump in blueshift prior to center of transit—the greater the maximum thickness of the cloud, the earlier in phase the jump occurs. The best-fitting model in this grid includes clouds up to 10 scale heights above the Al₂O₃ condensation curve (corresponding to a 0.7 mbar cloud top on the western limb and a 1.1 mbar cloud top on the eastern limb). A cloud top up to 10 scale heights above the Fe condensation curve performs nearly, but not quite as well. Even in the Al₂O₃ cloud case, though, the fit to the data is not perfect; namely, the magnitude of our Doppler shift is a few km s⁻¹ short of the observed data at egress (Figure 10, panel (c)).

Contrary to the cloudless case, including gas-phase condensation when clouds are also included increases the measured blueshift (Figures 9 and 10). This effect is often very slight; in the 10 scale height, Al₂O₃ cloud case, including condensation increases blueshift by 6 m s⁻¹ at center of transit.

The overall best fit to the Ehrenreich et al. (2020) data from our suite of models is obtained for the case of low drag (τ_drag = 10⁷), a deep RCB, a 10 scale height, optically thick Al₂O₃ cloud, and gas-phase iron condensation. Even this model, however, is not perfect. For instance, this model fails to match the egress blueshift magnitude of the Ehrenreich et al. (2020) data (Figure 7). We discuss these points in further detail in Section 5.

As shown in Figure 8, our 3D, physically self-consistent GCM output wind field is more complex than the substellar–antistellar and rotational wind field described by the Ehrenreich et al. (2020) toy model. For instance, blueshifts from winds are not uniform across the entire limb, and there even exists a redshifted annulus on the terminator (the bottom row of Figure 8), which represents return flow from the nightside to the dayside. ¹⁰ Furthermore, our GCMs are unable to reproduce the wind speeds required by the Ehrenreich et al. (2020) toy model. Our weak-drag GCM has an eastern limb average line-of-sight velocity of −6.5 km s⁻¹; this is in excess of the −5.3 km s⁻¹ day–night flow in the toy model. When averaged with the western limb, though, the average line-of-sight velocity reduces to −2.3 km s⁻¹—significantly less than the toy model wind field.

A crucial component of the Ehrenreich et al. (2020) toy model is an asymmetry in the Fe abundance distribution in the atmosphere of WASP-76b. Broadly, our iron abundance maps (Figure 5) indicate that, regardless of drag timescale, there does exist a projected asymmetry of iron on the eastern limb by the end of transit. The eastern limb is preferentially blueshifted (Figure 8), so more gas-phase iron on the eastern limb produces a larger measured blueshift. While this asymmetry is exacerbated by gas-phase condensation of Fe (the top two rows versus the next two rows of Figure 5), our self-consistent, GCM-based forward models with Fe condensation alone do not readily match either the magnitude or the shape of the Ehrenreich et al. (2020) observed blueshifts, and we find that we must invoke additional physics (e.g., clouds) in order to achieve a better fit.

As shown in Figures 11 and 12, our forward models generally provide good agreement with the multiple species detected across the entire ESPRESSO wavelength range by Tabernero et al. (2021). However, we bring up several noteworthy exceptions below.

Tabernero et al. (2021) report a detection of Hα in one of their two analyzed transits. In our Doppler-off, single-species template models, though, we note that the hydrogen Balmer lines do not protrude above the continuum (comprised of collision-induced absorption, H⁻ opacity, and Rayleigh scattering). However, when performing a blind search in our full-species spectrum, we find that there does exist an absorption line within the Hα spectral region. When compared to the Hα rest wavelength, this line (which we identify as a blend of TiO and VO) displays a blueshift relative to rest Hα that is 2σ consistent with the reported Hα data. We propose that either Hα was mistakenly identified by Tabernero et al. (2021) or that non-LTE, nonhydrostatic, and/or nonequilibrium chemistry effects may be boosting the strength of the Hα line compared to what is predicted by our GCM post-processing scheme.

Conversely, Tabernero et al. (2021) were not able to detect Cr or Ti in the WASP-76b high-resolution transmission spectrum, whereas our chemical equilibrium forward models are readily able to identify these species in cross-correlation. This mismatch may point to unaccounted-for opacity sources that wash out the Cr and Ti features, increased continuum absorption that masks the lines of these species (perhaps excess H⁻ brought about by photoionization, for example), or nonsolar abundances of Cr and Ti.

A notable result of the Tabernero et al. (2021) work is the height of their reported Ca ii lines. The Ca ii H line was detected at a height of 1.57 R_P, and the Ca ii K line was detected at a height of 1.78 R_P—both at high significance and at much higher altitudes than in our models. The error bars on these detections are large enough such that they are still consistent with our models' line heights at 2σ. However, the consistency of the observed Ca ii H&K line heights with one another supports their fidelity, as their similar line strength implies that they should become optically thick at roughly the same altitude. Furthermore, the deep absorption of the Ca ii triplet near 850 nm detected by two other instruments is supporting evidence for abundant high-altitude Ca ii in WASP-76b (Casasayas-Barris et al. 2021; Deibert et al. 2021). Finally, high-altitude Ca ii H&K lines have been detected in other ultrahot Jupiter atmospheres and have been attributed to hydrodynamic outflows and high-altitude photoionization (Yan et al. 2019; Borsa et al. 2021). Therefore, the discrepancy with respect to our models may be linked to the hydrostatic assumption of our GCM (see Section 5.2.2), which only extends up to 1.19 R_P. Exploring the outflow hypothesis is beyond the scope of this paper, and we encourage further nonhydrostatic investigations of this and similar planets.

Similar to this work, Wardenier et al. (2021) aimed to post-process a 3D GCM with a radiative transfer code to explore the role of Fe condensation in WASP-76b. Given that they use a different GCM (the nongray SPARC/MITgcm; Showman et al. 2009) and a different radiative transfer scheme (Monte Carlo photon-tracking) from those used in this work, comparing the basic results of both works serves to validate both approaches.

Our results are generally consistent with the 3D post-processing results of Wardenier et al. (2021). Both studies favor weak-drag over strong-drag models and resort to GCM modifications to reproduce the Doppler shift behavior reported by Ehrenreich et al. (2020) and Kesseli & Snellen (2021).

Other details of the modeling approach differ between our two works. For example, Wardenier et al. (2021) are able to use their nongray GCM to explore the effect of optical opacities on the predicted thermal structure. Our work explores a much larger spectral range than Wardenier et al. (2021), enabled by our more computationally efficient ray-striking radiative transfer approach. This method allows us to avoid potential biases incurred by only modeling a small portion of the Fe i band structure. (We find this to be up to a 1.5 km s⁻¹ effect that can further alter the Doppler shift trend over phase, due to Fe lines over a narrow wavelength range only probing a commensurately narrow range of altitudes.) We further model additional opacity sources beyond those considered by Wardenier et al. (2021).

Notably, Wardenier et al. (2021) were unable to reproduce the post-ingress jump in Doppler shift exhibited by the Ehrenreich et al. (2020) and Kesseli & Snellen (2021) data without artificially restricting the longitude range of iron condensation within their model domain or removing TiO and VO opacity from their GCM opacities (thus altering the thermal structure of the planet's atmosphere). We are able to produce the same jump in our simulations—effectively muting portions of the gas-phase Fe Doppler shift signal—but by using a more physically motivated model involving optically thick clouds (Figure 10). At no point do we actually change the GCM output (and hence the underlying physics); rather, we use post hoc, painted-on effects to reproduce the Doppler shift observations.

By comparing the results presented in this paper with those from our companion paper on the Spitzer phase curve of WASP-76b (May et al. 2021), we can attempt to seek out a set of globally consistent models that describe the ensemble of high-resolution and broadband photometric data. We find that models with cold interiors (i.e., a deep RCB) are favored by both sets of observations, as such a model allows for the cold nightside and large phase curve amplitude observed at 4.5 μm, constrained by the Spitzer study. However, the phase curve data prefer the strong-drag (τ_drag ≤ 10⁴ s) GCMs, driven largely by the near-zero phase curve offset observed in both Spitzer wavelength channels. This result is in contrast with our own conclusions that the high-resolution data can only be explained by weak-drag models, which allow for a significant east–west limb asymmetry.

The mixed success of our joint modeling efforts to consistently describe all observables—hotspot offsets, phase curve amplitudes, phase-resolved Fe blueshifts, and multispecies blueshifts—hints that our physical understanding of WASP-76b is as of yet incomplete.

We have shown that our forward models can reproduce much of the observed Doppler shift trend once optically thick clouds are included (Figures 7 and 10, panel (c)). The full magnitude of the observed Doppler shift, however, is not fully accounted for by this approach. The orbit of WASP-76b itself may provide the rest of the solution.

A Kozai–Lidov mechanism is often invoked in the literature to explain the migration of hot Jupiters to their current-day short-period orbit (for reviews, see Dawson & Johnson 2018 and Fortney et al. 2021). Any Kozai-like process, though, would result in vestigial eccentricity, as the mechanism presupposes a high past eccentricity that then enables tidal interactions to reduce the semimajor axis. Other formation mechanisms also have the potential to drive large eccentricities for Jupiter-mass planets as well—for instance, migration through the disk into a magnetospheric cavity could do so (Debras et al. 2021). Finally, hot Jupiters subject to other potential formation channels (even such as in situ formation; e.g., Batygin et al. 2016) could have initially low eccentricities excited by secular interactions (e.g., Wu & Lithwick 2011) or an external perturber (e.g., Zakamska & Tremaine 2004).

As demonstrated by Montalto et al. (2011), the effect of even a small (e = 0.01) unaccounted-for eccentricity can produce velocity shifts for hot Jupiters on the order of km s⁻¹. With respect to WASP-76b, West et al. (2016) placed a 3σ upper bound on eccentricity at 0.05. Fu et al. (2021) additionally constrained the planet's eccentricity as ${0.016}_{-0.011}^{+0.018}$. A nonzero—and, for our purposes, significant—eccentricity is therefore certainly within the realm of possibility, given the extant data.

In the analyses of Ehrenreich et al. (2020) and Tabernero et al. (2021), the eccentricity of the planet is held to 0. The authors' rationale appeals to a short (on the order of tens of Myr) circularization timescale, the age of the star (∼2 Gyr), a low stellar rotation rate (∼0.03 day), and observational constraints favoring small eccentricities. If we allow for even a small eccentricity in our models, however, we find that we do not need to appeal to ad hoc removal of iron at specific longitudes (the favored explanation of Wardenier et al. 2021) to produce the observed blueshift signature. Because WASP-76b subtends more than 30° of phase during its transit, including a small eccentricity can both increase the magnitude of our blueshift at every point in transit and produce a trend consistent with the Ehrenreich et al. (2020) results.

Using the exoplanet package (Foreman-Mackey et al. 2021) to model planet-frame radial velocity deviations from a circular orbit within a χ² fitting routine (see the Appendix for the fitting details), we find that the data are best explained by an eccentricity of e = 0.018 and a longitude of periastron ω = 125° for a cloudless atmosphere. ¹¹

But eccentricity alone cannot produce the data's characteristic jump in Doppler shift post-ingress that our optically thick cloud models reproduced. Combining 10 scale height Al₂O₃ clouds and condensation with a small eccentricity produces this work's best qualitative model fit (Figure 13, with e = 0.017 and ω = 117°). The quantitative fit also improves correspondingly: when including eccentricity, the χ² value for the clear case decreases from 401 to 121, and the χ² value for the cloudy case decreases from 329 to 112. The blueshift at the correct phases is increased by the eccentricity trend, thereby fully reproducing both the magnitude and time-dependence of the Ehrenreich et al. (2020) data.

Zoom In Zoom Out Reset image size

Both our cloudless and cloudy eccentricity-inclusive modeling results are consistent within 1σ with the ω derived from Fu et al. (2021; ${62}_{-82}^{{+67}^{\circ}}$). Furthermore, both estimates are well within the 3σ limit derived by West et al. (2016).

While the above eccentricity argument is perhaps satisfying, it may not accurately represent reality. Increasing blueshift as a function of phase has also been observed in the ultrahot Jupiter WASP-121b (Bourrier et al. 2020; Borsa et al. 2021). It is perhaps more likely that there exists an underlying physical mechanism common to both WASP-121b and WASP-76b than their orbital configurations happening to be similar; the parameter estimation of WASP-121b at discovery constrains e < 0.07 at 3σ (Delrez et al. 2016), and later studies have held the planet to a circular orbit (Evans et al. 2016, 2017; Mikal-Evans et al. 2019; Borsa et al. 2021). Furthermore, Bourrier et al. 2020 place a tighter constraint on eccentricity (e < 0.0078) at 3σ from TESS data. It is yet to be determined whether a post-processed GCM could account for the blueshift trend in WASP-121b without modifications such as added eccentricity (e.g., as in the weak-drag, deep RCB case in Figure 7). Hence, the question of whether cloudless models or the addition of optically thick clouds alone could provide a population-level explanation of phase-resolved phenomena warrants further investigation.

As explored in Section 5.1.2, the presence of nonhydrostatic effects in WASP-76b's atmosphere is supported by the behavior of the Ca ii H&K lines observed in Tabernero et al. (2021). In turn, it is possible that a hydrodynamic outflow may be responsible for some portion of the anomalous blueshift signature.

By virtue of its close proximity to its host star, WASP-76b should experience strong UV radiation flux; such irradiation can result in mass loss in its upper layers (e.g., Erkaev et al. 2007; Des Etangs et al. 2010; Linsky et al. 2010; Ehrenreich et al. 2015). Indeed, Seidel et al. (2021) recently found evidence for high-altitude vertical winds in WASP-76b, which may deliver material to a hydrodynamic outflow (or serve as the base of the outflow itself) in the planet's exosphere. As proposed for WASP-121b (Bourrier et al. 2020), an anisotropic expansion of the atmosphere could lead to variable Doppler shifts over the course of transit.

Both our ray-striking radiative transfer and our GCM presuppose that WASP-76b's atmosphere is in local hydrostatic equilibrium. While our assumption may be valid strictly in the pressure regimes that we consider in our study—simulations of atmospheric loss generally begin in the nanobar regime (e.g., Murray-Clay et al. 2009; Bourrier et al. 2015, 2016), whereas the upper boundary of our atmosphere is at a pressure of ≈11 μbar—enforcing the hydrostatic assumption prevents us from exploring this hypothesis. Fully nonhydrostatic models that include winds driven by UV radiation would be required to determine whether a physically plausible outflow could become optically thick enough at the correct phases to meaningfully contribute to the WASP-76b Doppler shift trend.