Exploring the Effects of Active Magnetic Drag in a GCM of the Ultra-Hot Jupiter WASP-76b

Ultra-hot Jupiters represent an exciting avenue for testing extreme physics and observing atmospheric circulation regimes not found in our solar system. Their high temperatures result in thermally ionized particles embedded in atmospheric winds interacting with the planet's interior magnetic field by generating current and experiencing bulk Lorentz force drag. Previous treatments of magnetic drag in 3D General Circulation Models (GCMs) of ultra-hot Jupiters have mostly been uniform drag timescales applied evenly throughout the planet, which neglects the strong spatial dependence of these magnetic effects. In this work, we apply our locally calculated active magnetic drag treatment in a GCM of the planet WASP-76b. We find the effects of this treatment to be most pronounced in the planet's upper atmosphere, where strong differences between the day and night side circulation are present. These circulation effects alter the resulting phase curves by reducing the hotspot offset and increasing the day-night flux contrast. We compare our models to Spitzer phase curves which imply a magnetic field of at least 3 G for the planet. We additionally contrast our results to uniform drag timescale models. This work highlights the need for more careful treatment of magnetic effects in atmospheric models of hot gas giants.


INTRODUCTION
Gas giant planets orbiting extremely close to their host stars make excellent targets for observers due to their favorable planet-star flux ratios and offer avenues for testing prescriptions of high-temperature physical processes in theoretical models. The quintessential category of these planets is ultra-hot Jupiters (UHJs), which orbit so close to their host star that their equilibrium temperatures exceed ∼2200 K. The extreme temperatures of these planets warrant a careful consideration of the relevant physical processes included in the models, as the increased irradiation results in different dominant mechanisms than those of their cooler cousins, "normal" hot Jupiters. While magnetic effects may begin to alter the circulation patterns of gas giant planets that reach temperatures ≈ 1500 K (Menou 2012;Rogers & Komacek 2014), it is not until planets reach an equilibrium temperature 2000 K that non-ideal MHD effects are predicted to become very strong due to the coupling between the circulation and atmospheric magnetic field, potentially resulting in effects such as hot spot reversals (Hindle et al. 2021a). Since no solar system analogue Corresponding author: Hayley Beltz hbeltz@umich.edu exists for these planets, intricate multi-dimensional atmospheric modeling is critical for understanding and interpreting observations of their atmospheres.
The consequences of magnetism manifest themselves on hot gas giants in a multitude of ways. The interior of the planet is expected to host a magnetic field of comparable or even greater strength than solar system Jovian planets (Yadav & Thorngren 2017). In the planet's upper atmosphere, high energy stellar photons can photo-ionize species and drive evaporative winds, whose outflow can be shaped by the planet's magnetic field (Owen & Adams 2014). Here, we focus on the atmospheric effects of thermal ionization of species near the planet's photosphere interacting with the planet's magnetic field. Due to the extreme temperatures on the daysides of these planets, species undergo thermal ionization while remaining embedded in the mostly neutral atmosphere; these ions then are advected around the planet via strong winds. The currents generated by this interaction could travel into the planet's interior and deposit heat there via Ohmic dissipation, perhaps explaining the inflated radii of many hot Jupiters (Batygin & Stevenson 2010;Perna et al. 2010a;Thorngren & Fortney 2018). The coupling between the ions and the mostly neutral winds results in a bulk Lorentz force drag, potentially reducing circulation efficiencies and increasing the day-night contrast (Perna et al. 2010b;Menou 2012;Batygin et al. 2013;Rauscher & Menou 2013;Rogers & Showman 2014). Current treatment of magnetic effects in 3D atmospheric models vary in complexity, from computationally expensive non-ideal MHD simulations to the use of a universal drag timescale applied throughout the planetary atmosphere, with this range in modeling complexity serving a variety of purposes. Perna et al. (2010b) introduced the idea of parameterizing the effects of the bulk Lorentz force on atmospheric winds with a drag timescale. The formulation of this timescale was derived from order of magnitude approximations of the terms in the non-ideal MHD induction equation in the hot Jupiter regime, thus effectively operating as a "kinematic" MHD model. This timescale was used in Rauscher & Menou (2013) in their models of two different hot Jupiters. This work found that for the hotter planet modeled, (HD 209458b, T eq ≈ 1500 K), the drag prescription resulted in slower wind speeds and a reduction in the strength of the equatorial jet often seen in hot Jupiter atmospheric models.
One simpler treatment for magnetic drag is applying a uniform drag timescale to the entire modeled atmosphere. Komacek & Showman (2016) explored the use of various uniform values for their model of a HD 209458blike hot Jupiter to approximate the effect of Lorentz forces. So long at the drag timescale was short enough ( 10 4 s), the zonal jet was eliminated and day-night temperature differences were large. Additionally, they noted that their drag timescale had only a secondary effect on the day-night temperature contrast when it was longer than the rotation rate of the planet. Uniform drag timescales are found in many works involving GCMs of hot (or ultra hot) Jupiters Koll & Komacek 2018;Kreidberg et al. 2018;Mansfield et al. 2018;Arcangeli et al. 2019). These universal drag timescales can be converted to global magnetic field strengths with order-of-magnitude estimates, such as in Kreidberg et al. (2018), which estimated a magnetic field stronger than ∼ 1 Gauss for the UHJ WASP-103b. This uniform drag treatment will also have an effect on the predicted phase curves of the planet, as shown in Tan & Komacek (2019). The resulting phase curves of their modeled planets had smaller offsets and larger amplitudes when the atmospheric drag was stronger. This means that increasing the strength of the drag resulted in hotter daysides and cooler nightsides on the planet as well as reduced the eastward shift of the planet's hotspot, as should be expected.
A limitation of universal drag timescales is that they do not allow for spatial variation in drag strength between the hot day and cold nightside of the planet: given the hundreds of Kelvin temperature contrasts, the magnetic resistivity should vary by many orders of magnitude, so we should expect effective magnetic timescales to also vary by orders of magnitude for a single pressure level (Rauscher & Menou 2013). Uniform drag timescales also do not account for the directional dependence of magnetic drag: assuming the planet's deepseated magnetic field is a dipole, only winds in the eastwest direction should experience Lorentz drag (Perna et al. 2010b).
The most complex treatments of magnetism in the atmosphere of a hot gas giant involve the use of magnetohydrodynamic (MHD) models. Rogers & Showman (2014) introduced the first non-ideal MHD simulation of a hot Jupiter. The authors compared the Lorentz force from their MHD simulations to the active magnetic drag timescale prescription from Rauscher & Menou (2013). Their Lorentz force peak value was within an order of magnitude of the timescale prescription, but the spatial extent of their Lorentz force was more localized. In a follow up paper, Rogers & Komacek (2014) ran a grid of MHD hot Jupiter models and found that time variability becomes important at the higher temperatures modeled (1400-1800 K) pointing to the intrinsic importance of feedback between magnetism and the planet's thermal structure. In the UHJ regime, Rogers (2017) offers the closest instance of a physically consistent nonideal magnetohydrodynamic treatment of an UHJ. This work, focused on HAT-P-7b, highlighted the complexity of dynamic magnetic field lines and estimated the minimum global magnetic field strength of the planet to be 6 Gauss, based on comparison to observed phase curve variability for this planet (Armstrong et al. 2016). This minimum field strength of 6 Gauss is also inferred from shallow-water magnetohydrodynamic models from Hindle et al. (2019).
Magnetism is not the only high-temperature physical process of note in these planets. If high enough temperatures are reached on the dayside of UHJs, H 2 is expected to dissociate, which results in a local cooling effect as dissociation requires an input of energy. As winds transport the gas to the nightside of the planet, the temperatures can drop enough that recombination occurs, locally heating the atmosphere and reducing the day-night temperature contrast (Bell & Cowan 2018). In Tan & Komacek (2019), the authors investigated the effects of molecular hydrogen dissociation and recombination, in addition to the presence of uniform drag, across a range of planetary equilibrium temperatures. Their work found that beginning near T eq = 2200 K, the dissociation and recombination of hydrogen begins to play a significant role in the circulation of the atmosphere.
Including hydrogen dissociation decreased the strength of the equatorial jet, which disappeared at the highest equilibrium temperature examined (T eq = 3600). Overall, including molecular hydrogen dissociation and recombination reduced the day-night temperature contrast, reducing the predicted phase curve amplitude.
Beyond H 2 , there are other molecular species whose dissociation has significant impact on the atmospheric structure of UHJs by changing the main sources of opacity. Parmentier et al. (2018) describes the influence of thermal dissociation and ionization of a variety of species including H 2 O and H − on the thermal and spectral properties of UHJs, using four UHJ models made from SPARC/MITgcm. (The only main molecular opacity source that did not dissociate anywhere in the UHJ atmospheres was CO.) These models showed stark differences between day and nightside temperatures and abundances. They found dissociation of water and presence of H − opacity may result in muted of spectral features. Additionally, when included as an opacity source, H − increases the optical depth across all wavelengths, causing the photosphere to move to a lower pressure (Lothringer et al. 2018).
Many UHJs exhibit what are known as temperature inversions-portions of the atmosphere that increase in temperature with decreasing pressure, opposite to standard behavior. These hot upper atmosphere inversions are likely due to high altitude optical/UV absorbers, such as TiO/VO (Hubeny et al. 2003;Fortney et al. 2008) or Fe, SiO, and metal hydrides (Lothringer et al. 2018;Gibson et al. 2020), but high C/O ratios (Mollière et al. 2015) or inefficient IR cooling (Gandhi & Madhusudhan 2019) could also play a role. The strength of this inversion may also scale with stellar host type (Baxter et al. 2020). Additionally, temperature inversions on UHJs are spatially inhomogeneous, as the stratospheric heating on the dayside will be absent from the shadowed nightside (Kreidberg et al. 2018).
Due to the high temperatures of the daysides of UHJs, clouds are expected to be mostly constrained to the nightside of the planets. Recently, Roman et al. (2021) explored cloud effects in our GCM for a variety of irradiation temperatures, including UHJs. This work found that in the UHJ regime (T irr > 3250K), 1 clouds were restricted to the nightside at high latitudes and were absent near the equator globally. Using a more complex cloud microphysics model, the GCM from Mansfield et al. (2018) was post-processed in Helling et al. (2019b), which concluded that some clouds would exist on the nightside of the UHJ HAT-P-7b and the dayside could host some very optically thin clouds, away from the equator. Their similar analysis of WASP-18b in Helling et al. (2019a) also found a cloud-free dayside and heterogeneous clouds on the nightside. As far as hazes are concerned, Helling et al. (2020) showed that the dayside of all UHJs are too hot for hydrocarbon hazes to be stable and that these hazes should play no role in the aerosol opacities on the nightside and terminator region.
One particularly interesting UHJ is WASP-76b. First detected by West et al. (2016), WASP-76b is an inflated ultra-hot Jupiter orbiting every 1.81 days with an equilibrium temperature just over 2200 K. Since its detection, a broadened sodium feature was detected in transmission spectra by Seidel et al. (2019), hinting at the super-rotation of the upper atmosphere. The sodium detection was confirmed shortly after inŽák et al. (2019). Ehrenreich et al. (2020) detected Fe absorption in transmission spectra from one side of the planet's terminator but not the other, which was interpreted as evidence for nightside condensation of the species. This detection was later confirmed by Kesseli & Snellen (2021). Recent work from Wardenier et al. (2021) suggests that this differential absorption could instead be explained by a strong temperature difference on the trailing and leading limb of the planet, without needing to invoke iron condensation. Fu et al. (2020) additionally detected TiO and H 2 O in transmission spectra form HST and Spitzer. From the emission spectra of the planet, CO emission features are present and their models suggested a temperature inversion of ∼ 500K. HST emission and transmission spectra also suggests the presence of TiO, H 2 O, and thermal inversions (Edwards et al. 2020). Tentative detections (4σ) of VO also exist (Tsiaras et al. 2018).
Perhaps due to their difficulty to model, (as a result of their extremely short radiative and dynamic timescales, in addition to the complicating physics discussed above) only a handful of GCMs for UHJs have been published. Because atmospheric dynamics will manifest in observables, 3D GCMs are useful in interpreting the spectra, (such as HAT-P-7b in Mansfield et al. 2018) phase curves, (such as the case of WASP-103b in Kreidberg et al. 2018) or both (see WASP-18b from Arcangeli et al. 2019). A common theme that arises from these works is that the GCM had difficulty reproducing UHJs with very low heat redistribution, which could potentially be due to an underestimation of magnetic drag strength, uncertainties surrounding dissociation effects, or warming from nightside clouds. When additional sources of drag are included, the resulting phase curves produced by the GCM in the works referenced above are a better match to observed values.
In this work, we explore the effects of active magnetic drag in our GCM of the UHJ, WASP-76b. Although all of the physical processes mentioned above are no doubt present to some degree on WASP-76b, we choose to not include clouds or H 2 dissociation effects at this time to focus solely on the influence of active magnetic drag in the planet's atmosphere. Future work is necessary to characterize the mutual interactions between these physical processes for UHJs. Based on observational constraints on the strength of the temperature inversion in this planet's atmosphere (Fu et al. 2020), we set up our model to include this stratosphere. In Section 2, we describe our GCM and the active magnetic drag treatment. In Section 3, we examine the atmospheric structures of our models with and without the active magnetic drag treatment. We additionally compute models featuring a universal drag timescale and compare the atmospheric effects of the two different treatments of magnetic drag. In Section 4, we contextualize our results and discuss the limitations of our models. Finally, in Section 5 we summarize the main points of this work.

METHODS
General Circulation Models (GCMs) are threedimensional numerical tools that simulate the underlying physics and circulation patterns of planetary atmospheres. To do this, GCMs solve the simplified set of fluid dynamics equations known as the primitive equations of meteorology. We use the GCM from Rauscher & Menou (2012) with the updated radiative transfer scheme in Roman & Rauscher (2017) (based on Toon et al. 1989). Our GCM solves the radiative transfer with a double-gray treatment, meaning that two absorption coefficients are used: one in the visible wavelength regime to account for absorption from the host star and one in the infrared regime for the planet's thermal emissionThe nuances of atmospheric implications for doubley-gray radiative transfer versus more complex treatments, such as correlated-k or picket fence methods, are beyond the scope of this work. Interested readers are directed to the recent paper from Lee et al. (2021) that explores in depth the consequences of these different radiative transfer schemes. We modeled the planet with 65 vertical layers evenly spaced in log pressure over 7 orders of magnitude from the bottom boundary of 100 bars, at a horizontal spectral T31 resolution (roughly 4 degrees at the equator) and the parameters listed in Table 1. Our absorption coefficients were informed by Fu et al. (2020); by changing the ratio of the infrared 4.018 × 10 −5 s −1 Substellar irradiation, Firr 5.14 × 10 6 W m −2 Planet internal heat flux, Fint 3500 W m −2 Optical absorption coefficient, κvis 2.4 × 10 −2 cm 2 g −1 Infrared absorption coefficient, κIR 1 × 10 −2 cm 2 g −1 Specific gas constant, R 3523 J kg −1 K −1 Ratio of gas constant to heat capacity, R/cp 0.286 and visible coefficients, our resultant atmosphere will have an inverted temperature profile at locations that are highly irradiated. Due to the increased difficulty of modeling planets in the ultra-hot regime, where heating rates are stronger and winds can be faster, we added sponge layers to the top three layers of all the models presented for numerical stability purposes. Sponge layers act as a damping mechanism in the top of atmospheric models to reduce the buildup of artificial noise brought about by atmospheric waves reflecting off of the top boundary of the model (Forget et al. 1999;Wills & Schneider 2016). The strength of these sponge layers decreases with height, with the strongest drag being applied in the top layer and then linearly decreasing with log pressure. By restricting the sponge layers to the top few levels, we leave the rest of the atmosphere unaffected. In our initial investigations, we found that after 1000 orbits were completed, the overall kinetic energy in the sponge layers was decreased by a factor of a few percent, and below the sponge layers, the change was even smaller. This is the first time sponge layers have been used in our GCM, but sponge layers are relatively common and have been implemented in other exoplanet GCMs (see Mayne et al. 2014;Deitrick et al. 2020;Wang & Wordsworth 2020, for examples).

Our Magnetic Drag Treatment
Our model's implementation of magnetic drag is unique among GCMs as it is calculated based on local atmospheric conditions and applied in a geometrically and energetically consistent way. The physical origin of this drag comes from thermally ionized particles interacting with the planet's magnetic field due to strong (mostly neutral) winds advecting the particles across the planet. From the atmosphere's perspective, the wind feels a bulk Lorentz drag force. We choose to model the case of a dipole field aligned with the planet's rotation axis as a simplifying assumption, because detailed information about the shape of exoplanet's magnetic field lines is not currently known. Additionally, we fix the shape of the magnetic field lines and assume the local field strength does not vary as a function of radius in the model. Our prescription for magnetism stands out from other GCMs because we use an active magnetic drag timescale; that is, one that is locally calculated and updated throughout the simulation to model the effects of magnetism in the planet's atmosphere. We can refer to this framework as a "kinematic" MHD framework. The critical assumption with this treatment is that we are assuming the dipolar magnetic field, generated in the planet's interior, is much stronger than any field induced in the atmosphere. As a result of this assumption, we apply the drag only in the zonal (east-west) direction. This assumption remains reasonable when the magnetic Reynolds number is less than 1 (Menou 2012;Hindle et al. 2021b,a). The magnetic Reynolds number can be approximated with R m ≈ U H η (Hindle et al. 2021a) where U is the zonal wind speed, H is the pressure scale height, and η is the magnetic resistivity. In this work, we calculate resistivity in the same way as Menou (2012): We calculate the ionization fraction, x e , using the Saha equation, taking into account the first ionization potential of all elements from hydrogen to nickel (as in Rauscher & Menou 2013). As shown on the left in Figure  1, R m < 1 for the majority of temperatures and pressures found in our model, implying that our assumptions are reasonable for much of the atmosphere, although we may be missing more complex behavior in the hottest dayside regions of the upper atmosphere. Although this is a simplification of non-ideal atmospheric MHD effects, it nevertheless represents a reasonable starting point for exploring the effects of active magnetic drag and is more complex than the uniform drag timescales applied in other GCMs.
Our model also has unique advantages compared to current non-ideal MHD simulations in the low magnetic Reynolds number regime. In models such as those found in Rogers & Komacek (2014), magnetic resistivity is calculated from a static, horizontally homogenous (no latitude or longitude dependence) initial temperature profile. Other work, such as Rogers (2017), calculates resistivity based on a static temperature profile with horizontal variations, but does not change as the simulation runs and local temperature conditions change. Our GCM calculates this resitivity locally and often, allow-ing the resitivity to evolve temporally in a more self consistent manner. By doing this, our resistivity varies by many orders of magnitude from the dayside to the nightside for a single pressure level and updates as the model progress, allowing for feedback with the atmospheric thermal structure. Additionally, our model is coupled with double-gray radiative transfer equations. Current non-ideal MHD models lack this radiative transfer coupling and instead employ a Newtonian relaxation for radiative forcing, a simplified scheme where the temperatures relax torward a prescribed profile at a chosen timescale (Rogers & Showman 2014;Rogers & Komacek 2014). Finally, our model's dynamical core does not have any explicit viscosity in the momentum equations, as atmospheres are generally highly inviscid. Because current non-ideal MHD models such as those discussed above have heritage in dynamo/interior models with explicit viscosity in the solved equations, their resulting wind speeds are up to an order of magnitude weaker than expected from inviscid GCMs in the hot Jupiter regime.
As described in detail in Rauscher & Menou (2013) and shown in Equation 2, the magnetic drag is applied by subtracting U/τ mag from the east-west momentum equation, where the magnetic timescale is calculated locally as: where B is the chosen global magnetic field strength, φ is the latitude, ρ is the density. The kinetic energy lost as a result of this drag is then returned to the atmosphere as localized ohmic heating in the energy equation as U 2 τmag , for both the active drag discussed here and uniform drag models discussed below. This local treatment of magnetic drag allows the magnetic timescale to vary by many orders of magnitude throughout the entire planet without being as computationally expensive as solving the full non-ideal MHD equations. Additionally, this treatment updates the timescales as the model progresses, allowing feedback between the magnetic effects and the atmosphere's thermal structure. Figure  1 (right shows the magnetic drag timescales across the atmosphere's temperature and pressure space from our 3 G model of WASP-76b. Across a single pressure level in the upper atmosphere, our magnetic drag timescale can vary by nearly 15 orders of magnitude. The shortest magnetic timescales are located on the dayside of the upper atmosphere, due to the high temperatures and low densities. Deeper in the atmosphere, where the density is greater, the timescales are longer and the effect of our magnetic drag on the circulation pattern will be smaller. This large variation in magnetic drag strength highlights the versatility of our active, locally calculated drag treatment.
Both τ mag and R −1 m scale linearly with the magnetic resistivity and so in Figure 1 we see some correlation between the two. However, the other physical variables that factor into the magnetic Reynolds number and drag timescale mean that they do not trace each other exactly. For example, almost all of the regions with R m > 1 also have τ mag < 10 5 s, but there are locations with τ mag = 10 5 s that have R m values as low as 10 −10 . Even the shortest drag timescales, τ mag ∼ 1 s, can exist in regions with R m ∼ 10 −2 , although generally those timescales are found in regions with R m ≥ 1.
We ran our model of WASP-76b for variety of different magnetic field strengths: 0 G, 0.3 G, 3 G, and 30 G to encompass the field strengths measured in our solar system and to explore the range of effects our magnetic treatment will have on our modeled atmospheres. We additionally ran two instances of a uniform drag with timescales equal to 10 4 and 10 7 seconds following the GCMs in Tan & Komacek (2019). All of these models were calculated for 2000 orbits, to allow the atmosphere enough time for winds to accelerate (starting from rest) and reach a steady state.

The Base Case: No Magnetic Effects
To understand the effects of magnetic drag on our models of WASP-76b, we begin by examining the "base" case atmosphere, free of magnetic drag through a snapshot of the atmospheric structure taken from the end of the simulation. Also, since this is the first planet with a forced inverted temperature profile we have published with this GCM, we will examine the extent and effects of the temperature inversion. Overall, our model of WASP-76b exhibits the typically expected hot Jupiter flow patterns seen in our previous works (e.g., Rauscher & Menou 2013;Beltz et al. 2021;Roman et al. 2021). These features include a strong, eastward equatorial jet present in the deeper atmosphere ( Figure 4) and an eastward advected hotspot. Figure 2 shows the temperature-pressure profiles for the planet as well as the 1-D profile we use globally to initialize the simulation and theoretical 1-D averages (based on Guillot 2010). From this figure, a few things are apparent. First, our model shows a large (>1500 K) day-to-night temperature contrast along the equator. This is in broad agreement with the ∼ 1440K brightness temperature contrast measured from Spitzer phase curves at 4.5 µm (May et al. 2021).
Second, while the temperature inversion must disappear on the nightside, as there is no stellar flux to be absorbed at high altitude, the east and west terminators are markedly different in the temperature inversions at those locations (roughly yellow and blue lines in Figure  2, respectively). Because low pressures heat and cool rapidly, the upper atmosphere temperatures are more similar at each terminator. The dominant eastward direction of the winds means that at deeper pressures the east terminator remains warm from gas heated on the dayside, while the west terminator has colder air advected around from the nightside.
Aside from the inverted temperature profile itself, another notable difference between this model and previous models with our GCM of non-inverted hot Jupiters is the decreased frequency and strength of chevron-shaped temperature features. These features are a result of vertical motions caused byconverging winds, which heats the atmosphere in the chevron-shaped pattern (e.g., Rauscher & Menou 2010) and can can result in the hottest spot on particular pressure levels to be east or west of the substellar point (e.g., Beltz et al. 2021). In our model of WASP-76b however, these chevron features are weak enough that the hottest point stays east of the substellar point at all pressures ( Figure 3) and the chevron features are less prominent than in our noninverted models.

The Magnetic Models
We can now begin to look at the differences between our active magnetic drag models and our drag-free models. Figure 3 shows horizontal temperature maps for these scenarios at various depths in the atmosphere. The wind fields are overlaid in white stream lines. The differences between the non-magnetic models and magnetic models are most stark at the lowest pressures. The wind fields on the dayside of the magnetic models only flow in the north-south direction. This is because our magnetic drag is applied only to the east-west momentum equation (see Section 2.1). At this low pressure and high temperature regime, the magnetic timescales are very short, even for the smallest magnetic field strength tested. This causes the elimination of east-west momentum, leaving only the weaker north-south wind vectors. As the pressure increases and we get deeper into the atmosphere, the choice of magnetic field strength results in more disparate atmospheric structures. At the second pressure level shown, P ∼ 1 mbar, the drag timescales become slightly longer, due to the increase in density. At this level, the weakest magnetic field, 0.3 G, shows some wind flowing in the east-west direction, but a majority of the flow remains in the north-south direction. The Figure 1. Magnetic Reynolds number for the 3 G model of WASP-76b (left) and the magnetic drag timescales (right) based on the atmospheric structure at the end our simulation. Our kinematic MHD framework is appropriate where Rm < 1, which is the case for the majority of the atmosphere. The hottest regions on the dayside upper atmosphere, where Rm is the largest, is where atmospheric circulation can induce a magnetic field component comparable to the global field strength. On the right, we show the range in values for our active drag timescale. Across a single pressure level, the magnetic drag timescale can vary by many orders of magnitude. Since our magnetic drag varies with latitude, a single temperature and pressure will have a range of timescales which explains the non-monotonic behavior in portions of the plot.
two strongest magnetic field cases on the other hand still show exclusively north-south flow. Going deeper into the atmosphere still, this pattern continues. At P ∼ 0.01 bar, the 0 G and 0.3 G cases are very similar because the magnetic timescales at this pressure level have become so large that there is only a weak influence on the flow. In the 0 G case, the equatorial jet is slightly stronger and broader. Additionally, the nightside vortices near ±50 • latitude are more pronounced than the 0.3 G model. The 3 G and 30 G cases still predominately have north-south winds on their dayside because the increased field strength allows the timescales to remain short enough to strongly affect the flow. At P ∼ 0.1 bar, the magnetic timescales become so long that only the strongest magnetic field tested shows a substantially different flow pattern than the non-magnetic case. At this pressure level, the 0 G case hosts the strongest equatorial jet, which becomes weaker as the magnetic field strength increases.
Increasing the magnetic field strength from 3G to 30 G switches the atmospheric circulation pattern from a global equatorial super-rotation found in all models < 30 G, to a dayside dominated by day-night flow and nightside eastward equatorial flow. It is interesting that this eastward flow is maintained on the nightside, even though there is no longer eastward flow on the directly forced dayside. While we leave a detailed analysis of this flow pattern to future work, it is worth noting that such behavior is not found in models with uniform drag timescales.  Guillot (2010). The dashed black line is the profile used to initialize conditions globally. A strong temperature inversion is present on the day side, which becomes non-existent on the un-illuminated nightside. The large day-to-night temperature contrast is also apparent.
The equatorial jet is also affected by the magnetic drag, as seen in Figure 4. The non-magnetic case has the strongest and deepest jet. When active drag is present, the jet is weakened, especially at low pressures. This can be explained by the low pressures hosting the shortest drag timescales, which inhibit east-west flow and instead directs day-to-night flow meridionally up and over the poles. The stronger the magnetic field, the weaker the corresponding jet. In the strongest magnetic field we tested, the jet almost completely disappears. The reduction of the equatorial jet will result in less of an eastward shift of the planet's hotspot, meaning that the resulting phase curves will experience a smaller peak offset.
In Figure 5, we show the temperature-pressure profiles for all of our active magnetic drag models. To first order, these profiles share many properties: the hottest parts of the dayside exceed 3000K, temperature inversions occur on the dayside but not the nightside, and the inversion varies in strength between the two extremes, most notably near the terminators. We can see the influence of the active magnetic drag in the nightside equatorial profiles; as we increase the strength of the magnetic field, these profiles become colder. The models with magnetic drag also show more similarity between their east and west terminator profiles, in constrast to the non-or weakly-magnetic models where there is still significant eastward advection of hot gas on the dayside and cool gas from the nightside, leading to a warmer east terminator and cooler west side. Additionally, we see that at higher magnetic field strengths, the temperature profiles in the lower atmosphere become more uniform.
We note that the increase in magnetic field strength also increases the amount of large temperature jumps between layers (zig-zag like behavior). This may be due to variability in heating and cooling at smaller spatial/temporal scales than our model resolvescould result from the numerical dissipation being too weak in these uppermost layers (Heng et al. 2011), although since the best numerical dissipation strength to use should likely vary over the conditions throughout the atmosphere (Thrastarson & Cho 2011), increasing the dissipation in these upper layers would likely overdamp the higher pressure regions. This variability is especially present with the dark blue curves, corresponding to near the terminator west of the substellar region. Given the fact that this region lies at the boundary between stellar irradiation and no irradiation, it is not surprising that these profiles exhibit these characteristics, nor are they particularly worrisome.

Uniform Timescale Drag Cases
We additionally wanted to compare the effects of our active magnetic drag treatment with uniform drag timescale models. Following Tan & Komacek (2019), we ran our model of WASP-76b with uniform drag of timescales 10 4 and 10 7 seconds. As in our active magnetic drag models, we included sponge layers for the top three layers in these models for numerical stability. The resulting temperature-pressure profiles for uniform drag are shown in Figure 6. Overall, these profiles are relatively similar to each other, especially in the upper atmosphere. The lower atmosphere (below ∼ 1 bar) becomes increasingly uniform in temperature for shorter drag timescales. In fact, due to the uniform application of magnetic drag, the shortest drag model is essentially a single temperature at each pressure level below ∼10 bars. We saw a similar effect in our active magnetic models, but to a lesser degree. Notably, for the model with the shorter uniform drag timescale (10 4 s), the thermal inversion is slightly stronger and extends deeper in the atmosphere. The longer uniform timescale (10 7 s) exhibits these effects to a smaller degree. This effect was not seen with our active magnetic drag models, where the timescales increase with increasing pressure. Figure 3. Maps of temperature and wind structure at several pressure levels (top to bottom: lowest to highest pressure) for each of the active magnetic drag cases studied (from left to right: increasing magnetic field strength). Each plot has the substellar point centered and the thickness of the wind vectors scale with the relative wind speeds; maximum wind speeds for each map are reported in Table 2. The effects of our active magnetic drag are most readily seen at the lowest pressure, where all the magnetic cases show exclusively north-south flow on the dayside, as the magnetic treatment has removed all the east-west momentum. As the pressure increases, the magnetic drag timescale also increases, reducing its effects. The stronger magnetic fields will have the shortest timescales and thus have the effects of magnetic drag present in deeper levels. . Longitudinally averaged east-west winds for the active magnetic drag cases studied. The black contours denote the boundary between positive (eastward) and negative (westward) winds. When the active magnetic drag is incorporated, the jet strength immediately diminishes, especially in the upper atmosphere. The stronger the magnetic field strength, the weaker the jet and the more it is constrained only to deeper pressure levels. Figure 5. Temperature-Pressure profiles for our active magnetic models. The rainbow hues denote equatorial profiles, while the grey curves show profiles for the rest of the planet. All models show strong day-night contrasts and and temperature inversions on the daysides. The increase in magnetic field strength causes the nightside temperatures to decrease and the lower atmospheres to become more uniform.
In Figure 7, we show a similar plot to Figure 3, but for the case of uniform drag timescale. The differences between the drag-free model and the uniform drag models are more subtle than for our active magnetic drag models with different assumed field strengths. Because the uniform drag is applied evenly throughout the atmosphere, the resulting wind structure is similar in shape but weaker in strength than the drag-free case for the longer timescale. For the shorter timescale, the flow pattern west of the substellar point actually shows some westward motion at nearly all latitudes at and above 0.1 bar. This behavior is seen in other GCMs with uniform drag and is a result of the frictional drag becoming the dominant term in the momentum equation (Komacek & Showman 2016;May et al. 2021). Figure 8 shows the longitudinally averaged eastwest winds (akin to Figure 4) for the uniform drag timescale cases. We can see that the uniform drag timescale prescription destroyed the equatorial jet for the shortest, 10 4 s timescale. Interestingly, some weaker, higher latitude eastward jets are present for most of atmosphere(P 1 bar). For the longer uniform drag timescale, we see that the jet has been weakenedto a point where there is actually a slight net westward flow Figure 6. Temperature-Pressure profiles for different uniform drag timescales. The shortest uniform drag timescale case (10 4 s) exhibits a stronger thermal inversion ending deeper in the atmosphere and a more uniform lower atmosphere compared to the drag-free case. Additionally, the shortest uniform drag timescale causes cooler nightside temperatures, as it diminishes all of the temperature-homogenizing flow instead of just the east-west component, as in our active magnetic drag models. The 10 7 s uniform drag timescale case shows these same effects, but to a lesser degree.
at very low pressures (P 0.1 mbar). However, from examining Figure 7, we see that the wind structure contains both significant east and westward flow, largely canceling each other out in the zonal average. .
We can make a more direct comparison of our uniform drag timescale models to the "cold interior" models of WASP-76b in May et al. (2021). For the weak drag timescale (10 7 s), the flow patterns and temperature structures between our two works (our Figure 7 and their Figure 5) are very similar. One difference to note is that our upper atmospheres tend to be hotter and our lower atmospheres are colder than those in May et al. (2021), likely due to their inclusion of hydrogen dissociation and recombination. We also note that the longitudinally averaged east-west wind speeds (our Figure 8 and their Figure 7) are broadly similar in overall wind pattern in the case of weak drag. Both models exhibit a strong equatorial jet of similar strength that extends to ∼ 1 bar. The model in May et al. (2021) exhibits westward flow at high latitudes on order of a few km/s that are not present in our model. Comparing the strong drag cases shows larger differences between our models. Although both models display the disruption of the equatorial jet at all pressure levels modeled, our model shows the emergence of net eastward flow at high latitudes extending to the top of the modeled atmosphere. The corresponding models in May et al. (2021) show westward flow for high latitudes at the top of the atmosphere (P < 10 −2 bars) and eastward flow of a smaller magni-tude at pressures lower than this., The general qualitative agreement between these uniform drag models (which differ in other parameter choices) suggests that the differences we find between our active and uniform drag models may be contextualized within the broader work of the community, where uniform drag timescales have been more commonly used. Figure 9 shows our calculated phase curves from all of the models. As we increase the magnetic field strength in our active magnetic drag models, the peak of the phase curve shifts nearer to 0.5 (secondary eclipse) and the amplitude is increased. This is the expected result of our magnetic drag disrupting the eastward advection of hot gas away from the substellar point and generally reducing the transport of heat to the nightside. The difference in the phase of peak flux between our 3 G and 30 G models is small, as in both cases the hot spot is essentially unmoved from the substellar point. The larger amplitude of the 30 G case can be attributed to a cooler nightside and a stronger suppression of eastward flow in the upper atmosphere. A similar pattern exists for the uniform drag timescale models. The longer uniform drag timescale peaks at a similar phase to our model with no magnetic (or other) drag applied, meaning that its hot spot is still efficiently advected east of the substellar point despite the drag applied. The phase curve for the shorter uniform drag timescale model peaks between the curves for our 0.3 and 3 G models, reflective of Figure 7. Maps of temperature and wind structure at the same pressure levels as Figure 3, but comparing the uniform drag timescale cases. The differences in atmospheric circulation caused by the uniform drag timescale are more subtle but most prominently seen in the weakening of the equatorial jet in the lower atmosphere.

Observable Consequences: Phase Curves
the increased ability of the drag to prevent advection of the hot spot. Because phase curves integrate flux over the entire visible hemisphere, many of the differences between our active magnetic cases and uniform are lost in this integration. Published Spitzer 4.5 micron phase curves of this planet show a very small hotspot offset (< 0.003 in orbital phase from 0.5, May et al. (2021)), making our 3 G and 30 G models most consistent with the data. Given the double-grey nature of our GCM, comparing our bolometric phase curve amplitudes and the Spitzer amplitudes in parts per million of the stellar flux is not trivial and imprecise, as it relies on the use of multiple assumptions. Thus, we refrain from those comparisons here. We do note however that our equatorial day-night contrasts are in broad agreement with the 1440K brightness temperature contrast reported in May et al. (2021). . Longitudinally averaged east-west wind speeds for the uniform drag timescale cases. The 10 7 s uniform drag timescale weakens the equatorial jet, while the shorter, 10 4 s timescale completely eliminates it. In this shorter timescale we see instead weaker, higher latitude jets, which also begin to appear in the 10 7 s drag timescale case. Figure 9. Left: orbital phase curves of total (bolometric) thermal emission for all of the models from this paper. Right: the amplitude ((Fmax − Fmin)/Fmax) and phase of peak flux for each model's phase curve. For our active magnetic drag cases, increasing the magnetic field strength increases the amplitude of the phase curve and causes it to peak closer to a phase of 0.5 (which is secondary eclipse, not shown). This reflects the diminished day-night heat transport and disrupted eastward flow, due to the influence of the drag. The uniform drag timescale cases show similar behavior. The longer uniform drag timescale case has a similar hot spot shift to our non-dragged, 0 G case but an amplitude more similar to our 0.3 G model, while the shorter drag timescale has a phase curve amplitude and peak between the 0.3 and 3 G active drag models.
Compared to the less complex treatment of magnetism as a uniform drag, our active magnetic drag has some advantages. First, we are able to see clearly the difference in the magnitude of magnetic effects on the dayside vs the nightside. The uniform drag scenario forces all sides of the planet to be treated equally, despite the large difference in temperature around the planet. Many of the small scale atmospheric effects of our active magnetic drag disappear when the drag becomes uniformly applied. Similarly to Kreidberg et al. (2018), we can roughly estimate the magnetic field strength that corresponds to uniform drag timescales used in this work. Because the same timescale is applied to the atmosphere regardless of temperature, the corresponding field strength on the dayside and the nightside will differ dramatically. Based on equation 12 from Perna et al. (2010b), we estimate that our shorter drag timescale (10 4 s) would correspond to a field strength of ∼6 G on the dayside equator and nearly 300 G on the nightside equator. 2 For the longer uniform drag timescale (10 7 s) we estimate the corresponding field strengths to be ∼0.2 G and ∼10 G on the dayside and nightside respectively. Physically, we do not expect the nightside field strength to be nearly two orders of magnitude stronger than the dayside. Allowing our magnetic resistivity to vary spatially and temporally is also an improvement to treatments found in current non-ideal MHD models which employ temporally fixed resistivity. Because our active magnetic drag treatment allows the drag strength to vary as a function of temperature and pressure, we are able to model magnetic resistivity in a more physically consistent fashion.
Second, it is important that our active magnetic drag is only applied to the east-west momentum equation as we are assuming a planetary dipole field. As a result, there is a fundamental change in the direction of atmospheric circulation. Our dayside flow patterns were very distinct from the 0 G model in the upper atmosphere, as shown in Figure 3. Uniform drag timescales, on the other hand, are applied to east-west and northsouth momentum equations. When the uniform drag timescale is long enough that equatorial super-rotation is still present, the corresponding flow patterns are similar in shape to the drag-free models but lower in strength. The shorter uniform drag model was able to disrupt this superrotation in a way that can produce slight net westward flow near the equator.
Our active magnetic drag prescription is able to produce dis-tinct differences in flow patterns that a uniform drag timescale is not.
Many hot and ultra-hot Jupiters have inflated radii, which is likely due to energy deposited in the internal adiabat of the planet (Batygin et al. 2011;Laughlin et al. 2011;Lopez & Fortney 2016;Komacek & Youdin 2017;Thorngren et al. 2021). Energy dissipation due to currents flowing in a resistive atmosphere has been theorized as a potential mechanism to deposit this energy (Batygin & Stevenson 2010;Perna et al. 2010a;Rogers & Showman 2014;Thorngren & Fortney 2018).
In our models, we return the kinetic energy lost from our magnetic drag timescales in the energy equation as u 2 /t mag , mimicking local Ohmic dissipation. Notably, the deepest pressure modeled in this work is only 100 bars, so we cannot self-consistently predict the amount of heating that would be deposited deeper. However, using results from 1D global models of Ohmic dissipation profiles (Huang & Cumming 2012;Wu & Lithwick 2013) as guidance, we might expect the Ohmic dissipation in the convective interior to be around 1% that of the Ohmic dissipation at the base of our dynamically active atmosphere. For each of global fields tested (0.3 G, 3.0 G, and 30 G) this approximates to 10 18 W, 10 20 W, amd 10 19 W respectively deposited in the interior, which correspond to a maximum value of ∼ 0.01% of the stellar irradiation for the 3G model. The internal heating needed to inflate a planet will be dependent on model assumptions and evolutionary history, but using the general estimate from Komacek & Youdin (2017), which suggests that internal heating rates 1% of stellar irradiation at a minimum depth of 100 bars can explain inflated hot Jupiter radii, we cannot conclude that our model's Ohmic heating could explain the inflated radius of WASP-76b. Future studies, particularly those that model deeper pressures than the ones presented here, will be needed to explore this issue more thoroughly and self-consistently.
In order to first focus on the effects of our active magnetic drag treatment, we have chosen to omit molecular hydrogen dissociation/recombination, thermal ionization, and clouds. All of these processes should be present on WASP-76b, so further work is necessary to better understand the interactions of these processes with our active magnetic drag. Although we did not include these effects in our models, we can estimate how they would alter our model based on results from other GCMs. One of the main effects expected from molecular hydrogen dissociation and recombination is a reduction of the day/night temperature contrast (Bell & Cowan 2018;Tan & Komacek 2019). Cooling the dayside and warming the nightside would likely work to reduce some of the effects of magnetism on the upper atmosphere daysides of our models. By increasing the length of our magnetic timescale on the dayside of the planet, we would expect the day flow pattern of only north-south winds to not persist as deep in the atmosphere. The pressure at which east-west flows emerge again on the dayside will likely still be determined to first order by the strength of the global magnetic field.
Based on the condensation curves from Roman et al. (2021), clouds would likely be confined entirely to the nightside, particularly at upper latitudes of the planet due to the high temperature of the dayside These clouds could work to blanket part of the nightside, trapping heat and warming the area by a couple hundred Kelvin below the cloud. This would decrease the magnetic drag timescales on these nightside pressure regions, but seeing that these timescales are already many orders of magnitude longer than the dayside timescales, the effects of our magnetic drag would still be predominately on the dayside of the planet. Furthermore, above the clouds, we can expect a similar amount of cooling, which would decrease the likelihood of the overall nightside flow patterns being significantly altered from the cloudfree models we present here.
Though we did not include clouds in these models, we can use the temperature profiles from Figure 2 and the condensations curves from Roman et al. (2021) to speculate what potential cloud species could exist on the nightside of the planet. Based on the temperature structure, KCl, ZnS, and MnS clouds are only possible at high latitudes in the upper atmosphere, at pressures lower than ∼ 10 −2 bars or so. The condensation curves of Cr 2 O 3 , SiO 2 , VO, and Mg 2 SiO 4 show that these species could exist as condensates on the nightside equator at pressures near ∼ 10 −1 bars and at high latitudes for lower pressures. Ni, Fe, Ca 2 SiO 4 , CaTiO 3 , and Al 2 O 3 are also potentially present at nearly all pressures modeled, but would be confined to the nightside of the planet.
Our active magnetic drag treatment makes simplifying assumptions regarding the shape of the planet's magnetic field (a dipole), its orientation (aligned with the rotation axis), and that it lacks any time variability. To more accurately represent the magnetic field of this planet, one potential route forward would be to incorporate non-ideal MHD equations, which are extremely computationally expensive in a 3D atmosphere. Work done at this level of complexity highlights the expected variability and asymmetry of the induced magnetic field lines and strength (Rogers & Komacek 2014;Rogers & McElwaine 2017), and westward advection of the hot spot (Rogers 2017;Hindle et al. 2021b)However, it is important to note that non-ideal MHD has never been coupled to a dynamical solver without explicit viscosity. Models with explicit viscosity return wind speeds significantly reduced compared to those found in primitive equation GCMs like the one presented here, even in the purely hydrodynamic case (Rogers & Showman 2014;Rogers & Komacek 2014). Additonally, current non-ideal MHD models do not currently incorporate the large day-night magnetic resistivity variations presented in our model. A useful path forward could involve improvements to our kinematic MHD model by relaxing our strict dipole assumption and accounting for the toroidal component of the field. This would increase the validity of our model in the regions where the magnetic Reynolds number exceeds 1.

CONCLUSION
In this work we presented GCMs of the ultra hot Jupiter WASP-76b, focusing on the role of magnetic drag in shaping this planet's atmospheric structure. We used the "active magnetic drag" treatment from Rauscher & Menou (2013), which calculates a drag timescales based on local conditions throughout the atmosphere and updates these values as the simulation runs. Our main results are as follows: • The influence of our active magnetic drag was most strongly seen in the upper atmosphere dayside of the planet, where the drag timescales were the shortest. The winds here were solely in the north-south direction, such that hot gas from the dayside mainly flowed to the nightside over the poles, as our magnetic drag reduces momentum in the east-west direction.
• The stronger the global magnetic field strength, the deeper in the atmosphere the effects of our magnetic drag were seen. Our drag treatment can dramatically disrupt the equatorial jet seen in hot Jupiter models, reducing its strength and causing it to be confined to the lower atmosphere, where the magnetic timescales become longer. At the strongest case we examined, 30 G, the equatorial jet was reduced significantly and the atmosphere no longer displayed superrotation.
• We calculated bolometric thermal emission phase curves from our GCMs and found that the models with active magnetic drag had larger amplitudes (due to less efficient day-night heat transport) and peaked nearer to secondary eclipse. In comparison to Spitzer phase curve observations of this planet (May et al. 2021), our models would imply a magnetic field strength of at least 3 G for WASP-76b.
We do note, however, that other important physical processes are missing from these models (hydrogen dissociation and clouds) and a full comparison with the data requires modeling all of these processes together.
• We also looked at the effects of using a uniform drag timescale as an estimate for magnetic drag. Unlike our active drag timescale, the uniform treatment applied the same drag timescale globally in the east-west and north-south momentum equations. The circulation patterns in the upper atmosphere of the uniform case differed significantly from our active magnetic drag treatment, especially on the dayside.The corresponding phase curves put the longer timescale (10 7 s) model between our drag-free and 0.3 G active model and the short timescale (10 4 s) between our 0.3 G and 3.0 G active models. Order of magnitude estimates based on local conditions of the dayside and nightside of the planet indicate field strengths from ∼6-300 G and ∼0.2-10 G for the short and long timescale respectively. We have shown that our active magnetic drag case does a better job of reproducing the expected physics of UHJs.
Our active magnetic drag prescription is an improvement over the use of uniform drag timescales to model the interaction between partially thermally ionized winds and a deep-seated planetary magnetic field, in that it captures more of the physical behavior we expect. However, we re-emphasize that the physically correct solution will require complex and computationally expensive full non-ideal MHD treatments in regions where the magnetic Reynolds number is greater than 1. In this regime, the induced atmospheric field becomes comparable or larger than the dipolar field. As a result of this feedback in atmospheric circulation, hot spot reversal and time-variability can occur (Rogers & Komacek 2014;Hindle et al. 2019Hindle et al. , 2021b Future changes to our model should improve on our kinematic MHD framework and incorporate the effects of the toroidal field on meridional winds to replicate these steady state features in the high Reynolds number regime. . Given the high observational interest for UHJs, our community should invest in further development and use of those models especially for UHJs with strong thermal ionization such as WASP-76b, HAT-P-7b, WASP-18b, and KELT-9b. With future work, we will more carefully examine the influence of active magnetic drag on various types of atmospheric characterization measurements, beyond just the bolometric phase curves shown here. In particular, high-resoultion spectroscopy may be a promising avenue for empirically constraining wind speed and drag mechanisms (Miller-Ricci Kempton & Rauscher 2012;Flowers et al. 2019;Beltz et al. 2021).