Magnetic Drag and 3D Effects in Theoretical High-resolution Emission Spectra of Ultrahot Jupiters: the Case of WASP-76b

Hayley Beltz; Emily Rauscher; Eliza M.-R. Kempton; Isaac Malsky; Grace Ochs; Mireya Arora; Arjun Savel

doi:10.3847/1538-3881/ac897b

1. Introduction

High-resolution spectroscopy (HRS; R ≳ 15,000) has opened up exciting new pathways into exoplanet atmospheres. By taking advantage of the tens of thousands of spectral lines separable at this higher resolution, astronomers have been able to detect scores of atmospheric species, despite the orders-of-magnitude flux contrasts between the planet and the star. In order to unlock this information, the HRS method uses a template spectrum for the planet and cross correlates this across wavelengths (which is equivalent to varying the Doppler velocity shift of the spectrum) with the data to search for a peak signal. The planet's orbital motion will cause its signal to shift with wavelength during observations. Since the orbital motion of the planet is known, this method can reject spurious low-significance signals and reliably pinpoint the planet's atmospheric signature. The atmospheric conditions used to calculate the template spectrum can be adjusted until they most closely match reality, producing the largest detected signal. HRS can also be used to detect atmospheric motion on the planet in the form of net redshifts or blueshifts (such as the first-ever HRS detection from Snellen et al. 2010), as well as constrain rotational velocities from Doppler broadening (Snellen et al. 2014; Brogi et al. 2016; Schwarz et al. 2016; Flowers et al. 2019) and even obliquity (Bryan et al. 2020). Recently, the most precise C/O ratio in a hot Jupiter (HJ) was measured with HRS techniques (Line et al. 2021) using the framework described in Brogi & Line (2019). Excellent reviews of HRS techniques can be found in Birkby (2018) and Brogi & Birkby (2021).

Crucial to the successful use of this method is having a template spectrum to use in the cross correlation that accurately matches the planet's true spectrum. For this work, we will focus on the two classes of model spectra: those generated from 1D atmospheric models and those generated from more physically consistent 3D models, each of which have their advantages and disadvantages. On the one hand, 1D models are feasible for computing over large parameter spaces in a reasonable amount of computing time. They have also been used with high-resolution emission spectra to detect a variety of species, including (but certainly not limited to) Fe, Ti, Mg, TiO, CO, OH, H₂O, and CH₄ (see Brogi et al. 2014; Gandhi & Madhusudhan 2018; Guilluy et al. 2019; Nugroho et al. 2020; Giacobbe et al. 2021; Kasper et al. 2021; Nugroho et al. 2021; Pelletier et al. 2021, for some recent examples). However, planets are 3D objects, and using a 1D template can lead to unintended consequences. At low and medium resolution, using a 1D model can bias interpretations, such as by incorrectly constraining abundances and temperature structures (Feng et al. 2016; Blecic et al. 2017; Caldas et al. 2019; Pluriel et al. 2020; Taylor et al. 2020). These biases should become more pronounced when planetary properties vary strongly as a function of longitude. In HRS observations, spectral lines are more fully resolved, allowing the 3D atmospheric structure to directly impact the shapes of the observed lines, calling into question the applicability of the 1D assumption.

3D models, on the other hand, are much more computationally expensive to calculate (generally requiring at least a few weeks to run to completion). 3D models are also by nature more complicated to calculate simulated spectra from, as the geometry and radiative transfer calculations are more complex than the 1D case (Zhang et al. 2017). (However, recent open-source tools, such as the ones presented in Lee et al. 2022, are making this barrier less steep.) Importantly, if the 3D model better matches reality, its spectrum will produce a stronger cross-correlation signal with the high-resolution data. In Beltz et al. (2021), we found that our 3D HJ model outperformed hundreds of 1D models when used in the analysis of high-resolution emission spectra. 3D models are also able to capitalize on the 3D effects present in HRS, such as the Doppler shifting of lines, temperature gradients, and chemical inhomogenieties (Flowers et al. 2019; MacDonald & Lewis 2022; Prinoth et al. 2021). Other recent work has also been investigating other routes for retrieving 3D structure from high-resolution emission spectra (Herman et al. 2022; van Sluijs et al. 2022)

3D models are especially necessary for ultrahot Jupiters (UHJs), given the large spatial variations across these planets, due to their extreme irradiation regime. These planets, with equilibrium temperatures that exceed ∼2200 K, orbit so close to their host star that they are expected to have tidally synchronized rotation states (Showman & Guillot 2002). This means that these planets have a constant dayside that is continually bombarded by stellar irradiation and a permanent nightside. An emerging ubiquitous property of UHJ atmospheres is that they exhibit dayside temperature inversions—that is, regions of decreasing temperature with increasing pressure—in contrast to their cooler HJ cousins (Baxter et al. 2020; Mansfield et al. 2021). The exact cause of these inversions is still debated, but high-altitude optical and UV absorbers such as TiO or VO (Hubeny et al. 2003; Fortney et al. 2008) or Fe, SiO, and metal hydrides (Lothringer et al. 2018; Gibson et al. 2020; Lothringer et al. 2022) are likely candidates. Alternative potential causes for the inversions include high C/O ratios (Mollière et al. 2015) or inefficient IR cooling (Gandhi & Madhusudhan 2019). These temperature inversions manifest themselves in the emitted spectra of the planet, where dayside spectra show emission features. The nightsides of these planets, due to their mostly monotonic temperature structure, will only show absorption features.

Due to the extreme irradiation received, the daysides of UHJs are hot enough for the thermal dissociation of important molecules. This can significantly change the main opacity sources (e.g., removing water and adding H⁻, such that the spectral features become more muted; Parmentier et al. 2018) or, in the most extreme cases, change the dominant atmospheric gas from molecular to atomic hydrogen (Bell & Cowan 2018). The dissociation of molecular hydrogen on the dayside and its recombination on the nightside should decrease the temperature gradient across the planet and decrease the expected wind speeds (Tan & Komacek 2019).

For this work, we focus our investigation on magnetism and its role in UHJ atmospheres. Magnetic field generation strength still holds uncertainties, but the work from Christensen et al. (2009) predicts that given a sufficient rotation rate (like that of tidally locked UHJs), the magnetic field strengths of planets may scale with energy flux. For the most inflated HJs, Yadav & Thorngren (2017) estimate the maximum surface field strengths to be between 50 and 100 G, with more massive planets having the potential for even stronger fields. While there have been no direct measurements of exoplanet magnetic field strength yet, several indirect methods can offer constraints for particular planets. Recently, Ben-Jaffel et al. (2021) detected the magnetosphere of the warm Neptune HAT-P-11b and estimated the equatorial magnetic field strength to be between 1 and 5 G. Magnetic planet star interactions, specifically Calcium II K emission, have also been used to indirectly estimate the surface magnetic fields of HJs in the range of 10–100 G (Cauley et al. 2019). Variability in hotspot location has been found in 3D nonideal MHD models, leading Rogers (2017) to estimate a minimum magnetic field strength of 6 Gauss for the UHJ HAP-P-7b, in order to explain the observed variability on that planet (it should be noted, however, that recent work from Lally & Vanderburg 2022 finds that the observed variation could be due to nonatmospheric sources). In our previous work, Beltz et al. (2022), we were able to estimate a minimum field strength of 3G for WASP-76b, based on comparisons with the phase curves published by May et al. (2021).

In our modeling, we remain agnostic to the intricacies of the planet's internal dynamo, and instead assume that a deep-seated global magnetic field of some magnitude exists, in order to focus on the effects of magnetic drag in the atmosphere of the representative UHJ WASP-76b and how that drag influences its emission spectra. Notably, our magnetic drag treatment does not solve the full nonideal MHD equations that, due to their added computational complexity, require other simplifying assumptions (see the discussion of this in Beltz et al. 2022). Instead, we use a parameterized approximation of magnetic drag first introduced in Perna et al. (2010a). A physical explanation of this parameterization can be thought of as follows: due to the partial thermal ionization of species on the daysides of UHJs, the strong winds on the planet will blow these charged particles (embedded in a mostly neutral flow) around the planet, crossing magnetic field lines and generating current. These currents, should they travel into the planet's interior and deposit heat, in a process known as ohmic dissipation, could play a role in the inflated radii of many HJs (Batygin & Stevenson 2010; Perna et al. 2010b; Thorngren & Fortney 2018). In the atmosphere, the bulk Lorentz force that ionized particles experience can reduce circulation efficiencies and therefore increase the day–night contrast (Perna et al. 2010a; Menou 2012; Batygin et al. 2013; Rogers & Showman 2014). Previous work has showed this in models of the HJs HD209458b and HD189733b (Rauscher & Menou 2013) and for the case of the UHJ WASP-76b (Beltz et al. 2022). Now, in this work, we take the models from Beltz et al. (2022) and postprocess their results to study the resulting high-resolution emission spectra.

Our models are all based on parameters from the UHJ WASP-76b, an inflated planet orbiting its host star every 1.81 days. This planet has been the target of many high-resolution transit analyses, and a multitude of different species have been detected in its atmosphere, including Fe I, Na, Li I, Ca II (Seidel et al. 2019; Casasayas-Barris et al. 2021; Deibert et al. 2021; Kawauchi et al. 2021; Tabernero et al. 2021), and recently molecules like OH (Landman et al. 2021). Analyses of high-resolution emission spectra have yet to be published. The work by Kesseli et al. (2021) explores an entire atomic spectral survey of the atoms and ions that have been detected on this UHJ, and finds evidence for temperature and chemical gradients on the planet. This planet has also garnered interest from astronomers due to an asymmetric detection of iron in its transmission spectrum throughout transit, first found by Ehrenreich et al. (2020) and later confirmed by Kesseli & Snellen (2021). Ehrenreich et al. (2020) proposed that this asymmetric absorption signal is a result of iron condensing on the nightside of the planet. Other proposed mechanisms for this absorption signal include temperature asymmetry between the trailing and leading limbs (Wardenier et al. 2021) or clouds (Savel et al. 2022). 3D atmospheric models of this planet have also been published in May et al. (2021), Savel et al. (2022), Wardenier et al. (2021), and Beltz et al. (2022).

In this work, we generate high-resolution emission spectra from previously published models of the UHJ WASP-76b, to explore how different magnetic drag treatments affect the resulting spectra and how these differences change as a function of phase. This particular planet is of high observational interest, but it also provides a route to making predictions regarding the emission spectra from UHJs more generally. The structure of this paper is as follows. In Section 2, we briefly describe the models used in this work, the variety of drag treatments in our models, and the process of postprocessing our model output to create simulated high-resolution emission spectra. In Section 3, we explore the features that show up in our 3D emission spectra, and how they differ between our models. In Section 4, we discuss important caveats to our models and assumptions. Finally, we summarize our main findings in Section 5.

2. Methods

2.1. General Circulation Model

In this work, we use a subset of the models produced and described in Beltz et al. (2022), and will briefly summarize the details here. General circulation models (GCMs) solve the primitive equations of meteorology in three dimensions, to simulate the physics and circulation patterns of planetary atmospheres. We use the GCM from Rauscher & Menou (2012), which underwent a radiative transfer upgrade in Roman & Rauscher (2017; based on Toon et al. 1989). The absorption coefficients for the WASP-76b models were informed by the 1D temperature–pressure profiles modeled in Fu et al. (2021), chosen such that the double-gray analytic profile (Guillot 2010) roughly matched those results. Every model was run with 65 vertical layers, spread out evenly in log space over seven orders of magnitude in pressure from 100 to 10⁻⁵ bar, and with a horizontal spectral resolution of T31, corresponding to ∼4° at the equator. Each simulation was calculated for a total of 2000 planetary days. The relevant global parameters are shown in Table 1.

Table 1. WASP-76b Model Parameters

Parameter	Value
Planet radius, R_p	1.31 × 10⁸ m
Gravitational acceleration, g	6.825 m s⁻²
Orbital period	1.81 days
Rotation rate,^a ω_orb	4.018 × 10⁻⁵ s⁻¹
Substellar irradiation, F_irr	5.14 × 10⁶ W m⁻²
Planet internal heat flux, F_int	3500 W m⁻²
Optical absorption coefficient, κ_vis	2.4 × 10⁻² cm² g⁻¹
Infrared absorption coefficient, κ_IR	1 × 10⁻² cm ² g⁻¹
Specific gas constant, R	3523 J kg⁻¹ K⁻¹
Ratio of gas constant to heat capacity, R/c_p	0.286

Note.

^aAssumed to be in a synchronous rotation state.

Download table as: ASCII Typeset image

Other commonalities of the models include the addition of "sponge layers" in the top three layers of the model. These layers apply extra drag near the top boundary of our model, to reduce the buildup of artificial numerical noise. These layers have little effect on the atmospheric flow below them and can be found in many GCMs (see Mayne et al. 2014; Deitrick et al. 2020; Wang & Wordsworth 2020, for examples).

2.2. Drag Treatment in Our Three Models

The three models of WASP-76b from which we have chosen to generate emission spectra have different treatments of drag: Drag-free/0G, Uniform/10⁴s, and Active Drag/3G. Here, we will summarize the main distinctions of these treatments. For further information about the detailed differences in circulation patterns that result from these drag treatments, see Beltz et al. (2022). We review the main observable differences in Section 3.1.

1.
Drag-free/0G: this model serves as a base comparison to the other two models, as there are no additional sources of drag aside from the sponge layers restricted to the top three layers and the numerical hyperdissipation that prevents the buildup of noise on the smallest scales (Thrastarson & Cho 2011). Both of these artificial sources of drag also exist in the other models and are standard features of many GCMs.
2.
Uniform/10⁴s: this model uses a simplified treatment of drag, in which the same drag timescale (in this case, 10⁴s) is applied globally at every level to remove momentum from the east–west and north–south components. When serving as a proxy for magnetic effects, however, this method has the unfortunate implication of assuming a stronger surface magnetic field (or stronger coupling) on the nightside of the planet than the dayside, contrary to what is expected physically (Beltz et al. 2022). Uniform drag is used in many GCMs, and it can also be used to model the effects of other large-scale atmospheric events, like hydrodynamic shocks or turbulence (Li & Goodman 2010; Perez-Becker & Showman 2013).
3.
Active Drag/3G: this treatment of drag is more physically consistent than the uniform model, and falls under the "kinematic" MHD umbrella. Instead of using a single global timescale, we locally calculate our active drag timescale based on the following expression from Perna et al. (2010a):
$\begin{eqnarray}&&{\tau }_{{mag}}(B,\rho ,T,\phi )=\displaystyle \frac{4\pi \rho \ \eta (\rho ,T)}{{B}^{2}| {\sin }(\phi )| },\end{eqnarray} \tag{ 1 }$
where B is the chosen global magnetic field strength (in this case, 3 G), ϕ is the latitude, ρ is the density, and the magnetic resistivity (η) is calculated in the same way as in Menou (2012):
$\begin{eqnarray}&&\eta =230\sqrt{T}/{x}_{e}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 2 }$
where the ionization fraction, x_e, is calculated from the Saha equation, taking into account the first ionization potential of all elements from hydrogen to nickel (as in Rauscher & Menou 2013). The derivation of this timescale assumes that the planet's interior generates a dipole field, aligned with the planetary rotation axis, that is significantly stronger than any induced magnetic field in the modeled atmosphere (Perna et al. 2010a). As a result, the active drag is only applied in the east–west direction. This active drag timescale was first used in GCMs of the HJs HD 189733b and HD 209458b (Rauscher & Menou 2013). In the UHJ regime, this timescale was used in Beltz et al. (2022), where we showed that the strong day–night temperature contrast of WASP-76b would result in the timescale varying by many orders of magnitude. For example, at the 10⁻³ bar level, the timescale varies from ∼ 10¹⁸ s on the cold nightside to ∼ 10³ s on the dayside (see the right panel of Figure 1 in Beltz et al. (2022) for the global distribution of timescales for our 3G model). This timescale is the most complex treatment of magnetic drag in the models presented here, but it still does not reach the complexity of implementing the full set of nonideal MHD equations.

**Figure 1.** Abundance contours of water (left) and CO (right) as a function of temperature and pressure, used to generate emission spectra, with dayside (red) and nightside (blue) equatorial profiles overplotted. The strong temperature difference between the dayside and nightside profiles results in a large difference in water abundances. CO, on the other hand, varies less as a function of temperature, and thus the CO abundances for the dayside and the nightside are similar.
Download figure:
Standard image High-resolution image

2.3. Generating Emission Spectra from the GCM

In order to determine how different magnetic drag parameterizations in our GCM of the UHJ WASP-76b manifest themselves in the high-resolution emission spectra, we take our model output and postprocess it with a detailed radiative transfer code. First, we regrid the GCM's 65 vertical layers in pressure to 250 layers at constant altitude, by assuming vertical hydrostatic equilibrium, consistent with the GCM.³ This is necessary in order to calculate line-of-sight columns and to increase the spatial resolution of the radiative transfer simulations. Based on resolution tests performed by Malsky et al. (2021), this number of altitude layers should be sufficient to accurately calculate the emergent spectra.

We perform the radiative transfer calculations for the postprocessed spectra following the methods in Zhang et al. (2017). We implement the two-stream approximation for inhomogeneous multiple scattering atmospheres from Toon et al. (1989). The 3D planet model is divided into 1D line-of-sight columns (48 latitude and 96 longitude points), and the outgoing intensity from each column is calculated independently. For each column, we calculate the total optical depth of each of the 250 layers. The model has the capability of simulating different cloud species (e.g., Harada et al. 2021), but no clouds were included in this work. Therefore, the total optical depth at each wavelength is simply

$\begin{eqnarray}&&{\tau }_{\lambda }=\int {\kappa }_{\mathrm{gas}}\ {dl},\end{eqnarray} \tag{ 3 }$

where κ_gas is the local gas opacity at each point along the line-of-sight path. The opacities are evaluated at their Doppler-shifted wavelengths, according to the line-of-sight velocity at each location:

$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{LOS}} & = & u\sin \theta +v\cos \theta \sin \phi -w\cos \theta \cos \phi \\ & & +\,{\rm{\Omega }}({R}_{p}+z)\sin \theta \cos \phi ,\end{array}\end{eqnarray} \tag{ 4 }$

where ϕ and θ correspond to the latitude and longitude, respectively, Ω is the planetary rotation rate, and u, v, and w are the wind speeds in the east–west, north–south, and vertical directions, respectively.

In calculating our spectra, we assumed local chemical equilibrium (LCE) and solar elemental abundances, which result in the amount of any particular molecular species varying as a function of temperature and pressure. In Figure 1, we show the abundance contours for water and CO (the two main sources of opacity at the wavelength ranges presented in this work), as well as equatorial dayside and nightside temperature–pressure profiles. Because of the strong temperature difference between the two hemispheres, we also expect a strong chemical gradient for some species as well. We can see from Figure 1 that the water abundances differ more strongly between the dayside and the nightside compared to CO. These differences in abundances will have implications for the line strengths in the corresponding emission spectra.

For this work, we calculated the spectra at two different wavelength ranges. The first of these spans from 1.14 to 1.35 μm, which is inside the wavelength range covered by the upcoming WINERED instrument (Ikeda et al. 2016). This wavelength range has features from H₂O, TiO, VO, K, and Na, and most of the spectra and results that we show here are from this wavelength range. The other wavelength range that we explored was 2.3–2.35 μm, which contains opacities from CO, H₂O, and TiO, and sits inside the range covered by IGRINS (Yuk et al. 2010; Park et al. 2014; Mace et al. 2016), CRIRES, and CRIRES+ (Kaeufl et al. 2004; Follert et al. 2014). We explore the differences between these two wavelength ranges in Section 3.4. The spectra were calculated from planetary snapshots at phases of every 11 fdg 25 for the entire orbit. We use the GCM outputs for the cases of 0 G, 3 G, and the shortest uniform drag timescale (10⁴ s) from Beltz et al. (2022). We calculate the spectra for two different conditions: one in which the spectral lines are broadened and shifted due to winds and rotation ("Doppler on" spectra) and one without these effects ("Doppler off"). All spectra were calculated at a resolution of R = 125,000.

3. Results

3.1. 3D Atmospheric Structure

Before diving into the spectra, we will briefly summarize the atmospheric structures of our models, as the temperature and wind structures will play a dominant role in shaping the resulting emission spectra. In Figure 2, we show the temperature and wind structures on the hemisphere facing an observer, at four phases throughout the planet's orbit, and for each of the three models examined in this work. The temperature and wind maps (shown at 10⁻⁴ bars, around where spectral line cores form) help to visualize the differences in the circulation patterns between the three models. The maps of where strong temperature inversions exist (calculated as per Harada et al. 2021) and the line-of-sight velocities (from winds and rotation) help to visualize which parts of the planet will emit spectra with emission or absorption features, and the net Doppler shift of those features.

**Figure 2.** Orthographic projections for the three models presented in this paper (in pairs of rows, from top to bottom: 0 G, 3 G, and uniform) at four orbital phases (from left to right: 0, 0.25, 0.50, and 0.75; transit would be at 0). The first row of each model shows the temperature structure at 10⁻⁴ bars (within the region probed by spectral line cores), with the wind directions plotted as arrows. The second row for each model shows the maximum vertical temperature inversion at each location. The blue and red contours show constant line-of-sight velocities in increments of ±2 km s⁻¹ at the same pressure level as the temperature plot immediately above it. The differences in temperature structures and wind patterns will influence the high-resolution emission spectra generated from these models.
Download figure:
Standard image High-resolution image

The differences in the temperature structures and atmospheric flow patterns between these models are discussed in detail in Beltz et al. (2022). The main changes that we saw once our active magnetic drag was turned on were an increase in the day–night temperature contrast, a decrease in the hotspot offset, and a change in the flow pattern. Under this active drag treatment, the dayside flow pattern switches from a day-to-night flow over most of the terminator to a flow that is channeled over the poles. Some interesting features not discussed previously are found in the line-of-sight velocities for our active magnetic drag model. Especially at the phase of 0.25, one can see redshifted regions near the terminator on the dayside blow up and over the pole, appearing as blueshifted regions on the nightside. In the uniform drag model, we see that significant changes to the circulation over the terminator to the west of the substellar point results in some changes in the line-of-sight velocity patterns, particularly at phases of 0 and 0.75.

3.2. Spectra as a Function of Orbital Phase

In Figure 3, we compare spectra from our three different models at three orbital phases: 0 (where the observer only sees the nightside of the planet), 0.25 (the first quadrature, where part of dayside and part of the nightside are in view), and 0.50 (where the viewer sees only the dayside of the planet). Due to the higher temperature of the dayside, the spectra corresponding to a phase of 0.50 have a much higher continuum flux than the nightside spectra. The continuum flux level for each model is set by the disk-integrated photospheric brightness temperature at each phase. We can see the largest day–night temperature difference in the 3 G model, followed by the uniform model, and the 0 G model, consistent with the simulated orbital phase curves from these models (Beltz et al. 2022). Similarly, the 0 G model has a higher continuum flux at an orbital phase of 0.25, compared to the dragged models, since it more efficiently advects hot gas from the dayside toward the eastern terminator, which is the region in view at this phase.

In Figure 4, we show spectra from our 3 G model at a variety of phases, both with and without the Doppler effects from winds and rotation.⁴ The spectra show an obvious switch between absorption features and emission features as different parts of the planet come into view. While the spectral lines begin in absorption, when the nightside of the planet faces us (at a phase of 0), by a phase of 0.375, enough of the dayside has come into view that the feature switches to emission, as there are temperature inversions found throughout the dayside in our model. Eventually, the features switch back to absorption at a phase of 0.875, as enough of the nightside comes back into view. Another interesting feature occurs at the phases of 0.25 and 0.75, where both absorption and emission features are present. Additionally, the emission features are blueshifted and the absorption features are redshifted at the orbital phase of 0.25, while the opposite directions are found in the shifts at 0.75.

We can understand this behavior by looking at Figure 2 and the animated Figure 5, which highlights the spectra and corresponding temperature maps at each phase. When looking at the line-of-sight velocities at a phase of 0.25, for example, we see both redshifted and blueshifted components near the limbs, which results in the emission features possessing a different net shift than the absorption features coming from the opposite limb. This effect would only be present in spectra of sufficient spectral resolution, as lower resolutions could result in flat, featureless spectra at these phases. We can roughly estimate the minimum spectral resolution needed to detect both the absorption and emission components by using $R\sim \tfrac{\lambda }{{\rm{\Delta }}\lambda }$ , where λ is the line center of the non-Doppler-shifted spectra and Δλ is the difference between the peak of the emission feature and the trough of the absorption feature. When applied to the feature shown in Figure 4, this roughly corresponds to a minimum resolution of ∼50,000. Note that in the spectra calculated without Doppler shifts, the absorption and emission features almost balance out, resulting in largely featureless spectra. By taking Doppler shifts at sufficient resolution into account, we do not lose the spectral line information. The high-resolution emission spectra that we present have inherently 3D effects, which would not be produced from a 1D model. These effects are especially important and become more computationally complex when considering nontransiting planets (as shown in Malsky et al. 2021). In nontransiting cases, even at phase of 0.5, part of the nightside hemisphere will be visible, which influences the resulting high-resolution spectra. As more high-resolution observations are taken, it will be critical to take into account multidimensional effects.

Figure 5. A comparison between the temperature and line-of-sight velocity structure visible on the planet (right) and the disk-integrated emission spectrum (left, highlighted) at each phase from our 3 G model. As the planet rotates and different regions come into view, the corresponding spectra will vary accordingly. Spectra with emission features are emitted from the hot dayside regions, where temperature inversions exist, while the cooler nightside regions emit spectra with absorption features. The line-of-sight velocities, due to winds and rotation, modulate the net Doppler shift of these different spectral components. This figure is available as an animation in the HTML version of the final article. The 10 s animation shows the planet rotating from a phase of 0 (shown in the static image) throughout its entire orbit. As the planet rotates, the temperature and wind fields change, and the corresponding emission spectra is highlighted on the left. During the first half of the animation, the hotter dayside comes into view, increasing the continuum level of the emission spectra. At a phase near 0.5, the dayside of the planet is in view, and the corresponding emission spectrum has the highest flux continuum level. Additionally, as the dayside of the planet comes into view, the spectral features switch from absorption lines to emission lines. From phases 0.5 to 1.0, more of the nightside starts to come into view, and the continuum level of the flux decreases and the spectral features switch back to absorption lines.

(An animation of this figure is available.)

Download figure:

Video Standard image High-resolution image

While Figures 3 and 4 allow us to view by eye the change from absorption features to emission features as a function of phase, it is helpful to examine this trend quantitatively. By generating spectra with and without water, we can subtract one from the other to isolate the water features, both in emission (positive) and absorption (negative). Using numerical integration, the absorption and emission components are summed at each phase. Spectra without Doppler effects were used to generate this plot, so all differences are a result of differing temperature structures between the models. In Figure 6, we show the integrated line flux (positive for emission, negative for absorption) for water features (the most prominent source of opacity within the 1.14 μm to 1.35 μm wavelength range) from our three models. The shapes of all three curves match the basic predictions: emission peaks at an orbital phase of 0.50 and absorption peaks at a phase of 0. Interestingly, the uniform drag model consistently shows the strongest integrated emission, while the 0 G model has the strongest integrated absorption. Another detail to note from this figure is that the models do not transition from being emission-dominated to absorption-dominated at the same phase. These differences in the intensity and orbital phase dependence of the spectral features are a result of complexities in the models' corresponding temperature structures, which vary based on the strength and type of magnetic drag applied.

**Figure 6.** The integrated line flux for the three models in this paper, as a function of orbital phase. From this plot, we see that the emission lines from the uniform drag model are consistently the strongest in the emission-dominated regime and the 0 G model has the strongest absorption features. Another interesting aspect is that the model spectra switch between being dominated by emission vs. absorption lines at different orbital phases.
Download figure:
Standard image High-resolution image

3.3. Cross-correlation Results

Another way of combining information from multiple lines within the planet's spectra is to use cross correlation to effectively integrate over all spectral lines, which is typically how signals are extracted from HRS observations. As described previously, two types of spectra are generated at each phase: Doppler on (which contains broadening and line shifting from winds and planetary rotation) and Doppler off (which contains none of this broadening). By cross correlating the Doppler on version of the spectrum with the Doppler off version (in velocity space), the peak of the cross-correlation function (CCF) will tell us the net Doppler shift of the signal, while the width of the CCF quantifies the total amount of broadening in the spectrum (see Figure 9 in Beltz et al. 2021, for example).

Figure 7 shows a series of CCFs from each of our three GCMs. At each phase, different parts of the planet's temperature and velocity fields are within the observable hemisphere, which influences the broadening and net shift of the spectra. Zhang et al. (2017) first calculated simulated high-resolution emission spectra from 3D models of three HJs (without any added drag) and found some similarities in the patterns of net Doppler shift versus orbital phase between those planets (their Figure 7). We see that our Doppler shifts reach higher absolute values, likely due to the more extreme winds and faster rotation of our UHJ; our 0 G model is the most similar to their results for WASP-43b, the fastest rotator (and one of the more highly irradiated) of their set. Similar to Zhang et al. (2017), the net Doppler shift varies in sign and magnitude throughout orbit, as different parts of the planet rotate into view. However, the spectra presented here have their peak redshifts slightly later in phase (∼0.75) compared to the ones from Zhang et al. (2017), which peak in redshift closer to a phase near 0.61.

These net Doppler shifts are informed by the temperature structure of the planet, which we see in Figure 2 and the animated Figure 5. For example, looking at a phase of 0, we see that the the drag-free model has a warmer, brighter region on the left, because the equatorial jet has advected that gas around from the dayside, which aligns with the blueshifted velocity structure, resulting in a net blueshifted CCF value at this phase. In contrast, the temperature patterns on the nightsides of the dragged models are more complex, with cooler regions somewhat associated with blueshifted structure, ultimately resulting in net redshifts. At phases where much of the dayside is in view (∼0.25–0.75), the spectra from the drag-free and uniform drag models show similar behavior in their net Doppler shifts, with a slow trend of increasing net shift with orbital phase (i.e., the spectra move from blueshifted to redshifted). Our active drag model, however, displays more complex behavior, where this overall trend is disrupted by a superimposed net redshift to blueshift pattern right around 0.50. At phases where much of the nightside is in view (<0.25 and >0.75), the drag-free model shows opposite net shifts than the 3 G and uniform models. This is a significant result: the net Doppler shift in the nightside spectra may be indicative of the presence of magnetic drag (if redshifted) or not (if blueshifted). A thorough analysis of the velocity precision for HRS emission spectra has yet to be done, but the theoretical work on simulated CRIRES data in Brogi & Line (2019) constrains the corresponding velocities of water and CO with errors as small as ±0.93 km s⁻¹ and ±0.26 km s⁻¹, respectively. At this level of precision, the differences between the Doppler shifts of our models would be feasible to detect. However, we do caution that there may be some wavelength dependence to these predictions (see Section 3.4). This is complementary to the evidence for magnetic drag that could appear in net Doppler shift versus line strength trends in the transmission spectra (Miller-Ricci Kempton & Rauscher 2012).

Figure 8 highlights the danger of ignoring the multidimensionality of planets, as expressed in their emission spectra. The teal curve, which reproduces the peak location of the CCF at each phase, as shown in Figure 7, shows the "true" function of the peak velocity shift as a function of orbital phase. The orange curve was generated by cross correlating a single dayside template spectrum—calculated from the 3D model at a phase of 0.50—with the spectra from our 3D model at each phase. The purple curve similarly does this with a single nightside spectrum—corresponding to a phase of 0. At phases near the template spectrum's origin phase, the peak velocity shifts match well with the teal curve. However, the farther away from the origin phase the single template is used, the larger the discrepancy from the true peak velocity. Most notably, at particular phases, such as 0.25 and 0.75, the dayside and nightside templates return net velocity shifts of different signs at very high net velocity, as these template spectra are just picking out the most redshifted or blueshifted component of the observed hemisphere, as the dayside or nightside region rotates in or out of view (see Figure 2).

This figure serves as a warning to observers to ensure that the template spectra that they use in their data analysis are generated at phases close to the observations. Longer observations over large ranges of orbital phase should use multiple templates for best results and to avoid misleading retrieved net Doppler shifts. While here we find different velocity components due to the spatial variation in the planet's thermal structure, a similar effect resulting from the spatial variation in chemical abundances may have been seen in Cont et al. (2021), who detected two atmospheric species with significant differences in associated velocities in the UHJ WASP-33b. They interpret this large velocity offset as a result of the depletion of TiO (but not Fe) near the substellar point, while both species exist on the terminator of the planet. Thus, their cross-correlation technique may have detected TiO only near the terminator with large line-of-sight velocity components, leading to a higher velocity offset than the more evenly distributed Fe.

3.4. Cross Correlation at Different Wavelength Ranges

As discussed in Section 2, we calculated spectra at two different wavelength ranges: 1.135–1.355 μm (wavelength 1, λ₁), used in all analyses up until this point, and 2.3–2.35 μm (wavelength 2, λ₂). In Figure 9, we show the Doppler shifts of the three different models at the two different wavelength ranges. Notably, the velocity shifts differ by some amount at each orbital phase tested, in some cases altering the trends identified above. This is a result of the different dominant atmospheric species in the two wavelengths ranges. In λ₁, the dominant absorbing species is water, while in λ₂ it is CO. These two species may be probing different atmospheric pressure levels, which can have different corresponding wind speeds and will result in different net Doppler shifts. In addition, the water abundance varies strongly as a function of temperature, while the CO abundance does not (see Figure 1). Since the water opacity relative to the CO opacity will therefore be location-dependent, this means that the difference between the pressures probed by each wavelength range should also vary with the spatial location around the planet. This is a good demonstration that the 3D structure of the planet's atmosphere influences its emission spectrum in subtle and complex ways.

**Figure 9.** Net Doppler shifts vs. orbital phase for emission spectra over the two wavelength ranges modeled. The colored lines show the Doppler shifts for wavelength range one, where water is the dominant absorber. The gray lines correspond to the shifts of wavelength range two, where CO is the dominant absorber. The different wavelength ranges probe different pressures in the atmosphere and, due to the temperature dependence of the water abundance, the disparity between the wavelength ranges should be spatially dependent. As such, the net Doppler shifts differ slightly. The bottom right panel shows this difference in net Doppler shift as a function of phase.
Download figure:
Standard image High-resolution image

4. Discussion

Important caveats to our GCM are already discussed in great detail in Beltz et al. (2022). We will briefly mention these points here, and invite the reader to examine that work for a more detailed explanation.

Because we have chosen to focus on the effects of magnetic drag in our modeling and postprocessing, we are unable to include all the physical processes expected to be present on UHJs. Potentially, the most consequential of these is molecular hydrogen dissociation and recombination. This would reduce one of the main effects that we see from our active magnetic drag, in that it would result in a reduction of the day–night temperature contrast (Bell & Cowan 2018; Tan & Komacek 2019). We also chose not to include clouds in these models, but work from Roman et al. (2021) suggests that if present, cloud coverage would be very minimal and confined to the nightside upper latitudes, due to the high temperatures on the planet. The inclusion of these clouds could alter the nightside temperature structure, potentially warming and cooling at different atmospheric heights. Additionally, depending on the height and opacity of the cloud deck, the strength of the absorption features coming from this region of the planet could be muted. The strength of these effects would vary as a function of phase, but would likely be most influential near the quadrature phases (0.25 and 0.75), where emission features and absorption features are comparable in magnitude and potentially lead to these spectra becoming more emission-dominated.

An important caveat to keep in mind when interpreting these results is that our magnetic drag prescription is a physically motivated approximation of the nonideal MHD equations. Our work assumes that the toroidal component of the magnetic field—that is, magnetic field lines in the latitudinal direction—is negligible. Instead, based on our dipole assumption, our magnetic field lines are solely in the poloidal (longitudinal) direction. Models devoted to solving MHD equations, like the ones presented in Rogers & Komacek (2014) and Hindle et al. (2021), treat both of these magnetic field components, but have their own limitations, such as precalculated local conductivities (that do not self-consistently update as the simulation runs), using a simplified treatment for heating and cooling due to radiative fluxes, and the inclusion of explicit viscosity (resulting in significantly slower wind speeds). Our models and postprocessed spectra presented here straddle the bridge between GCMs using simplified uniform drag and dedicated nonideal MHD models. Our work highlights the ongoing need for model-to-model comparison while still presenting testable observables to gauge the validity of our approximations. Specifically, the results from Figures 7, 8, and 9 are examples of the type of analysis that, if compared to observational data, can test our magnetic drag scheme's ability to approximate the physical environment of UHJs.

One of the crucial assumptions in our drag treatment is that any magnetic field induced in the atmosphere is much smaller than the global magnetic field, a relationship that holds when the magnetic Reynolds number R_m < 1. As we show in Figure 1 of Beltz et al. (2022), this assumption holds true for the vast majority of the atmosphere, except for some parts of the dayside near the low-pressure boundary of our model. To more physically consistently describe the magnetic effects in this region requires full nonideal MHD. Some 3D nonideal MHD models exist for (U)HJ atmospheres (Rogers & Komacek 2014; Rogers 2017; Rogers & McElwaine 2017; Hindle et al. 2021), and show that variability in the circulation pattern (including the westward flow near the substellar point) and the strength of the induced magnetic field may occur in the regime of high magnetic Reynolds numbers.

Since the double-gray radiative transfer in our GCM is agnostic to the specifics of the atmospheric composition, it was necessary to make assumptions regarding the abundances of species in the atmosphere for our calculation of the emission spectra. For simplicity's sake, we assumed solar metallicity with abundances set by LCE. Measuring the bulk metallicity for an UHJ is difficult, but the measurements of water composition (Kreidberg et al. 2014) and CO (Line et al. 2021), as well as the results from interior models (Thorngren et al. 2016; Cridland et al. 2019; Thorngren & Fortney 2019), suggest that our choices of solar abundances are reasonable. The works on this planet by Deibert et al. (2021), Landman et al. (2021), and Wardenier et al. (2021) all assume solar metallicity for their analysis. May et al. (2021) suggest that a lower bulk metallicity could explain the low internal heat flux, which is more consistent with the observed phase curve. On the other hand, Tabernero et al. (2021) retrieve a slightly higher than solar metallicity for the host star, which Casasayas-Barris et al. (2021) use as justification to also test supersolar composition models. With this uncertainty in mind, more work is necessary to truly understand the metallicity of this planet. As it pertains to our results, increasing the metallicity could result in moving the photosphere to slightly lower pressure in our GCM, perhaps with minor changes to the circulation pattern, and would lead to stronger spectral features in the emission spectra. These changes would likely not alter the main results and conclusions presented here.

Another important aspect to consider is that in our calculations we have assumed local thermal equilibrium (LTE). This assumption is true for a majority of the atmosphere modeled, except perhaps near the top boundary of our model. Deibert et al. (2021) found their NLTE model of WASP-76b to be a slightly better fit than their LTE model, but both models significantly underestimated the line depths. Its important to note, however, that their work used transmission spectra, which probe higher up in the atmosphere than emission, and therefore we would expect fewer NLTE effects in the emission spectra presented here. Other work on NLTE effects in UHJs has focused on the hottest known exoplanet to date, KELT-9b. For pressures greater than ∼ 10⁻⁶ bar (which encapsulates all of the pressures probed in our GCM), LTE and NLTE models were nearly identical in predicting the hydrogen level populations. In Fossati et al. (2021), they use NLTE as a framework for pressures less than 10⁻⁴ bars when modeling temperature–pressure profiles. For pressures greater than 10⁻⁴ bars, they note that the retrieved temperature profiles are very similar. Based on these results, we expect that LTE is a valid approximation for the pressures modeled in our GCM. It is tempting to extrapolate the trends in net Doppler shifts for the different models presented here to uniquely identify evidence of magnetic effects or estimates of magnetic field strengths. However, even if we take the assumption that the planet's magnetic field is a dipole aligned with its rotational axis to be completely true, these are Doppler shifts for one planet at one wavelength range. The additional physics not included in our model, such as clouds, hydrogen dissociation, or more complex magnetic feedback mechanisms, could alter the strength of the predicted Doppler shift. Future work exploring a variety of field strengths and sweeping over a range of planetary parameters would be necessary to draw more quantitative conclusions about the minimum field strengths. We can note, however, that our 3G model is the closest match to the phase curve published in May et al. (2021). One feasible potential observational test from this work would be the case where, if the trend in net Doppler shift versus orbital phase can be measured to sufficient precision, we could differentiate between whether the drag can be approximated with a uniform timescale or whether it is better matched by our more complex treatment. Future observational and theoretical work is therefore needed to determine how robust these predictions are.

5. Conclusion

UHJs are planets with strong temperature and compositional gradients, due to their extremely short orbital periods, and thus are excellent examples of the critical need for 3D models to fully understand the planet's emission spectra. In this work, we explore these 3D effects and illuminate how different forms of drag manifest in high-resolution emission spectra. Our main points are as follows.

1.
The strong temperature and chemical gradient on WASP-76b results in high-resolution emission spectra that vary strongly as a function of phase, showing at times absorption features, emission features, or both. These are inherently 3D features, in particular, as the emission and absorption components may have different net Doppler shifts, due to their emergence from regions of the planet with different line-of-sight velocities.
2.
The thermally ionized atmosphere of a UHJ should experience a magnetic drag acting on its winds; this influences the flow pattern in the upper atmosphere and the Doppler shifts in the planet's emission spectra. We predict that the presence of this drag should result in a slight net redshift in the nightside emission spectra, compared to the slight net blueshift seen in the drag-free model (both shifts ∼1 km s⁻¹). In addition, the complex form of magnetic drag may be apparent, over a simpler uniform-timescale drag, in the trend of net Doppler shift versus phase when the dayside of the planet is in view (showing a reversal in behavior around the secondary eclipse). We caution that these predictions are dependent on the wavelength of the observations.
3.
Using a single atmospheric template for cross correlation at multiple phases may provide poor matches to the planet, as the spectral features change inherently with the viewing geometry. This can bias the recovered net Doppler shifts, especially at phases where parts of both the dayside and the nightside are visible. In the most egregious cases, the amplitude of the calculated net Doppler shift can misleadingly exceed several km s⁻¹, when the inverted or noninverted region of the planet (producing emission or absorption features, respectively) only appears at the edge of the planet disk, where the line-of-sight velocities from winds and rotation are maximized.

The extreme spatial variations in UHJ atmospheres, exacerbated by the complication of magnetic drag, influence their high-resolution emission spectra. While these signatures can be subtle and highly dependent on the detailed 3D structure of the atmosphere, they may in turn offer us an avenue for empirically constraining those detailed physical properties and processes, possibly with current instruments, but likely within the reach of future measurements with Extremely Large Telescopes.

This work was generously supported by a grant from the Heising-Simons Foundation. Many of the calculations in this paper made use of the Great Lakes High Performance Computing Cluster maintained by the University of Michigan. The authors also thank the reviewer for their helpful comments that improved the quality of this work.

Magnetic Drag and 3D Effects in Theoretical High-resolution Emission Spectra of Ultrahot Jupiters: the Case of WASP-76b

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction