The Magnetic Mechanism for Hotspot Reversals in Hot Jupiter Atmospheres

Magnetically-driven hotspot variations (which are tied to atmospheric wind variations) in hot Jupiters are studied using non-linear numerical simulations of a shallow-water magnetohydrodynamic (SWMHD) system and a linear analysis of equatorial SWMHD waves. In hydrodynamic models, mid-to-high latitude geostrophic circulations are known to cause a net west-to-east equatorial thermal energy transfer, which drives hotspot offsets eastward. We find that a strong toroidal magnetic field can obstruct these energy transporting circulations. This results in winds aligning with the magnetic field and generates westward Lorentz force accelerations in hotspot regions, ultimately causing westward hotspot offsets. In the subsequent linear analysis we find that this reversal mechanism has an equatorial wave analogy in terms of the planetary scale equatorial magneto-Rossby waves. We compare our findings to three-dimensional MHD simulations, both quantitively and qualitatively, identifying the link between the mechanics of magnetically-driven hotspot and wind reversals. We use the developed theory to identify physically-motivated reversal criteria, which can be used to place constraints on the magnetic fields of ultra-hot Jupiters with observed westward hotspots.


INTRODUCTION
In recent years the field of exoplanetary research has greatly developed its understanding of exoplanet characterisation both observationally and theoretically. The field has now reached the point where light curves, infrared photometry, and spectra from spaced-based telescopes can be used to test, inform, and update our understanding of the atmospheric dynamics of these closely-orbiting gas giants.
However, recent observations suggest eastward hotspots may not be found ubiquitously, particularly on the hottest hot Jupiters (ultra-hot Jupiters). Continuous optical Kepler measurements find east-west brightspot oscillations on the ultra-hot Jupiters HAT-P-7b (Armstrong et al. 2016) and Kepler-76b (Jackson et al. 2019); optical phase curve measurements from TESS find westward brightspot offsets on the ultra-hot Jupiter WASP-33b (von Essen et al. 2020) 1 ; while thermal phase curve measurements from Spitzer find westward hotspots on the ultra-hot Jupiter WASP-12b (Bell et al. 2019) and the cooler hot Jupiter CoRoT-2b (Dang et al. 2018). There are three main explanations for these observations: reflections from cloud asymmetries confounding optical measurements (Demory et al. 2013;Lee et al. 2016;Parmentier et al. 2016;Roman & Rauscher 2017), asynchronous rotation (Rauscher & Kempton 2014), and magnetism (Rogers & Komacek 2014;Rogers 2017;Hindle et al. 2019). Ultra-hot Jupiters generally have nearzero eccentricities and are thought to be tidally-locked, so are expected to be synchronously rotating. They are also expected to have cloud-free daysides, where their atmospheres are too hot for condensates to form. Helling et al. (2019) recently ruled out cloud asymmetries as the explanation for westward brightspots on HAT-P-7b.
Using three-dimensional (3D) magnetohydrodynamic (MHD) studies, Rogers & Komacek (2014) predicted that magnetic fields could cause wind variations that drive east-west hotspot oscillations. Rogers (2017) then showed that the westward venturing hotspot displacements on the ultra-hot Jupiter HAT-P-7b can be well explained by the moderate deep-seated dipolar magnetic field strengths that are expected to be generated in the convective interior of such planets. In Hindle et al. (2019) we used a shallow-water MHD (SWMHD) model to show, firstly, that the magnetically-driven hotspot reversal mechanism is a shallow phenomenon that is driven by the flow's interaction with the planet's atmospheric toroidal magnetic field; and secondly, that the SWMHD model also requires a moderate planetary dipolar magnetic field strength to drive westward hotspot displacements on HAT-P-7b but that an excessively strong deep-seated dipolar magnetic field is required to reverse flows within the cooler (and hence less thermally-ionised) atmosphere of CoRoT-2b. The westward hotspot offsets on CoRoT-2b are therefore more plausibly explained by non-magnetic phenomena. Interestingly, the hot Jupiters Kepler-76b, WASP-12b, and WASP-33b are of the ultra-hot type so are more akin to HAT-P-7b than CoRoT-2b, making magnetically-driven reversals plausible for these observations. While 3D MHD simulations have proved crucial for identifying that magnetism can drive hotspot reversals in ultra-hot Jupiters, their dynamics is often too subtile and complex to glean physical understanding from. The aim of this study is to use a reduced physics model, alongside known features of 3D MHD simulations, to identify the mechanism by which magnetism can reverse hotspots in ultra-hot Jupiters.
In Sections 2 and 3, we use the numerical two-layer Cartesian SWMHD model of Hindle et al. (2019), with an equatorial beta-plane treatment of the Coriolis effect and different purely azimuthal equatoriallyantisymmetric initial magnetic field treatments, to study the westward transition of hotspots. In Section 4, we examine the link between magnetically-driven wind reversals and equatorial wave dynamics. Finally, in Section 5 we collate our findings, compare them to results of 3D MHD simulations, and present a physically-motivated hotspot reversal criterion for ultra-hot Jupiters.

NON-LINEAR SHALLOW-WATER MODEL
Three-dimensional models are fundamental to understanding the general features and flow behaviours of planetary atmospheres. However, with so many physical processes in play, it can be difficult to isolate the mechanisms responsible for driving a given flow pattern. In such instances, simplified models can be used to reduce the number of physical processes involved, exposing the underlying physics responsible for specific dynamical features. In this section, we present a detailed description of the reduced-gravity SWMHD model briefly described in Hindle et al. (2019), which we will use to explore the physics of wind reversals.

Governing equations
The reduced-gravity SWMHD model, is an adaptation of the SWMHD model of Gilman (2000) and is the MHD analogue of its hydrodynamic namesake (e.g., Vallis 2006), which has been used extensively in hydrodynamic studies of hot Jupiters (Langton & Laughlin 2007;Showman & Polvani 2010;Showman et al. 2013;Perez-Becker & Showman 2013). It is made up of two constant density fluid layers: a shallow active upper layer and an infinitely-deep inactive lower layer, which has no pressure gradients, velocity fields or induced magnetic fields in the horizontal direction (see Hindle et al. 2019, for a model schematic). Physically, the upper layer represents the meteorologically active upper atmosphere and the lower layer represents the deep atmosphere and deep interior of a hot Jupiter. The interface between the two layers is a material surface over which no magnetic flux is permitted to cross. When the system's length scales approach the shallow-water limit (i.e., if typical active layer horizontal scales, L, are much larger than the active layer's thickness, H), the vertical momentum equation of the full 3D system approaches magnetohydrostatic balance. This limiting approximation may be used together with the model's interface constraints to vertically integrate the 3D MHD equations over the vertical coordinate, z, to yield a shallow-water model with vertically independent variables (see Gilman 2000;Hindle et al. 2019, for further discussion). Using Cartesian horizontal spatial coordinates, (x, y), the dynamical behaviour of the active layer can be described by the following governing equations: where u(x, y, t) ≡ (u, v), is the horizontal active layer fluid velocity, h(x, y, t) is the active layer thickness which is used as the model's temperature proxy (see below), B(x, y, t) ≡ (B x , B y ) is the horizontal active layer magnetic field (in velocity units), and A(x, y, t) is the magnetic flux function of the active layer. We comment that the magnetic flux function definition differs from its two-dimensional definition through the inclusion of h in Equation (4). This arises as the magnetic flux function describes the vertically-integrated horizontal magnetic field over the whole fluid column, rather than simply the horizontal magnetic field at a specific vertical level. We use ∇ ≡ (∂ x , ∂ y ) to define the horizontal gradient operator, d/dt ≡ ∂/∂t + u · ∇ to define the Lagrangian time derivative operator, and ∇ × A z ≡ (∂ y A, −∂ x A) is the horizontal curl of the scalar field A about the vertical coordinate.
Defining the system in terms of the magnetic flux function guarantees that the SWMHD divergence-free condition, ∇ · (hB) = 0, remains satisfied throughout the domain at all times. This is the shallow-water analogue of Gauss' law of magnetism, which excludes magnetic monopoles. This shallow-water divergence-free condition is obtained by integrating the full 3D form of Gauss' law over the vertical coordinate, while imposing zero magnetic flux constraints across our model's layer interfaces. Using this formulation also highlights that in the absence of magnetic diffusion (D η = 0), A is a materially conserved quantity (see Equation (3) with D η = 0).
For numerical stability, we apply the following explicit diffusion prescriptions (Gilbert et al. 2014, A. D. Gilbert et al. 2021: where ν is the kinematic viscosity and η is the magnetic diffusivity. Geometrically, we fix a local Cartesian coordinate system about the equator, with −Rπ ≤ x < Rπ and −Rπ/2 < y < Rπ/2. We centre the system about the planet's substellar point, so x/R approximately corresponds to the azimuthal coordinate and y/R approximately corresponds to the latitudinal coordinate. Rotational effects are included via the so-called equatorial beta-plane approximation of Rossby (1939). Specifically, the only effects of sphericity the equatorial beta-plane approximation captures are the dynamical effects caused by latitudinal variations in the planetary rotation vector's vertical component. The approximation also uses the fact that in equatorial regions the Coriolis parameter, f , is approximately linear to set f = βy, where the constant β = 2Ω/R is the local latitudinal variation of the Coriolis parameter at the equator, Ω is the planetary rotation rate and R is the planetary radius.
The system is driven by a Newtonian cooling treatment, Q, in the continuity equation (Equation (2)), which relaxes the system towards the prescribed radiative equilibrium thickness profile, h eq , over a radiative timescale, τ rad . The Newtonian cooling is implemented with where H is the system's reference active layer thickness at radiative equilibrium and ∆h eq is the difference in h eq between this reference thickness and the radiative equilibrium layer thickness at the substellar point. This profile is similar to the spherical forcing prescriptions used in comparable hydrodynamic models (e.g., Shell & Held 2004;Langton & Laughlin 2007;Showman & Polvani 2010Showman et al. 2012;Perez-Becker & Showman 2013). The transfer of mass caused by Q generates horizontal pressure gradients, which drive recirculation via the generation of planetary scale shallow-water waves. Similarly, in three dimensional models, pressure gradients caused by heating drive recirculation via internal gravity waves. Using this analogy, mass sources and sinks represent heating and cooling respectively. This connection has been used extensively in hydrodynamic models of hot Jupiters, with active layer geopotential, gh, used as a proxy for specific thermal energy (Langton & Laughlin 2007;Showman & Polvani 2010;Showman et al. 2013;Perez-Becker & Showman 2013). Using this physical link, we equate the model's active layer reference geopotential, gH, to the reference thermal energy, RT eq , of the modelled planet's atmosphere, where R and T eq respectively denote the specific gas constant and the equilibrium reference temperature. With the addition of Q, a vertical mass transport term, R, needs to be introduced to enforce specific momentum conservation. In "cooling" regions (Q < 0) mass sinks from the active layer to the quiescent layer and causes no active layer accelerations 2 . However, in "heating" regions (Q > 0) mass transport causes decel-eration of the active layer as motionless fluid is transferred upwards. This deceleration due to heating is calculated by requiring specific momentum conservation in the active layer, yielding which has also been used in the hydrodynamic version of this model (e.g., Shell & Held 2004;Showman & Polvani 2010Showman et al. 2012;Perez-Becker & Showman 2013). We parameterise atmospheric drag with a linear Rayleigh drag treatment, −u/τ drag , where τ drag is the timescale of the dominant horizontal drag process in the thin active layer. Previous hydrodynamic studies use this Rayleigh drag to parameterise Lorentz forces (e.g., Rauscher & Menou 2013) or basal drag at the bottom of the radiative zone (e.g., Held & Suarez 1994;Liu & Showman 2013;Komacek & Showman 2016). In our study, we include Lorentz forces explicitly. However, due to the geometry of the SWMHD model, we only explicitly include the Lorentz forces caused by the atmospheric toroidal magnetic field (see Section 2.2 for a discussion of the magnetic field geometry in the atmosphere). We hence use the Rayleigh drag treatment to parameterise the Lorentz forces caused by planet's deep-seated poloidal magnetic field, which are not included explicitly. This is consistent with the treatment proposed by , whose τ drag parameterisation was based on estimating the direct influence that the planet's deep-seated poloidal magnetic field has on zonal flows. Though one could argue that in this setting the Rayleigh drag should have no meridional component, for comparison with past hydrodynamic results, we follow the commonly applied treatment of using Rayleigh drag in both horizontal directions (e.g., Showman & Polvani 2011;Rauscher & Menou 2013;Perez-Becker & Showman 2013). 3 We also comment that Rogers & Komacek (2014) found that magnetically driven wind variations emerge in the upper radiative atmosphere (where basal drags are negligible), so we do not consider basal drag in this work.

Magnetic field profile
The extension of planetary dynamo theory into the hot Jupiter regime is not well understood. That said, from current dynamo theory one would expect hot Jupiters to have planetary dynamos that are sustained within the convective deep interior, generating deepseated poloidal magnetic fields. The hottest hot Jupiters also have weakly-ionised atmospheres. If the atmospheres are sufficiently ionised, the zonally-dominated atmospheric flows become sufficiently connected to the planet's deep-seated poloidal magnetic field to induce a strong toroidal field that dominates the atmospheric magnetic field geometry (Menou 2012). Assuming this picture, and the planet's deep-seated magnetic field's geometry is dominated by an axial dipole, the induction of the toroidal component of the magnetic field can be approximated by where is the 3D gradient operator, and the electric currents generating the dipolar planetary field are implicitly assumed to be located far below the atmospheric region of interest Batygin et al. 2011;Menou 2012). Therefore, if toroidal field induction dominates toroidal field diffusion, the atmospheric toroidal field profile is expected to be equatoriallyantisymmetric, as found in the simulations of Rogers & Komacek (2014). There are not enough degrees of freedom in the SWMHD induction equation to simultaneously model the planetary dipolar field and the atmospheric toroidal field, so we only model the dominant atmospheric toroidal field self-consistently. We choose to enforce the simple equatorially-antisymmetric, purely azimuthal, initial magnetic field: where V A is the constant parameter that sets the magnitude of the azimuthal magnetic field. This profile may appear an unintuitive choice at first, but London (2017) noted that it has the useful properties for wave dynamics, which we shall exploit in Section 4. It is monotonic, behaves linearly in the equatorial region, and is bounded as y/L eq → ∞. The approximately linear latitudinal dependence of B 0 in the equatorial region means one can choose V A in accordance with the first order Taylor expansion of non-monotonic equatorially-antisymmetric profiles. Upon comparing to other field profiles, we generally find that doing so reproduces similar equatorial dynamics. To illustrate this, in Section 3 we compare some basic results to the profile B 0 = V A (y/L eq ) exp(1/2 − y 2 /2L 2 eq ), which is the equatorially-antisymmetric profile used in Hindle et al. (2019). This has the same first order Taylor expansion as Equation (10), has the maximum B 0 = V A at y = L eq (i.e., V A is the maximal initial Alfvén speed), and can be motivated from both Equation (9) and the simulations of Rogers & Komacek (2014). We implement the initial magnetic field profile of Equation (10) across an initially flat layer (h(x, y, 0) = H, everywhere), using the initial magnetic flux function, A 0 (y) = HV A L eq e 1/2 ln(cosh(y/L eq ))).

Numerical method and parameter choices
Numerical solutions are obtained by evolving Equations (1) to (4) from an initial uniformly-flat rest state (i.e., h(x, 0) = H, u(x, 0) = 0), in the presence of a purely azimuthal magnetic field (A(x, 0) = A 0 (y)). For hydrodynamic solutions we evolve until steady-state is achieved and for MHD solutions we run for a magnetic diffusion timescale. The system is solved on a 256 × 511 x-y grid, using an adaptive third-order Adam-Bashforth time-stepping scheme (Cattaneo et al. 2003), with spatial derivatives taken pseudo-spectrally in x and using a fourth-order finite difference scheme in y. We use periodic boundary conditions on u, h, and A in the x direction. On the y boundaries we impose v = 0 (impermeability), ∂u/∂y = 0 (stress-free), ∂A/∂x = 0 (no normal magnetic flux), and maintain the total columnar horizontal magnetic flux of the system. These conditions do not fix values of h on the y boundaries, which are up-dated to satisfy a consistency condition that results from mass conservation and our other boundary conditions. 4 We choose simulation parameters based on the planetary parameters of HAT-P-7b, an ultra-hot Jupiter with observed east-west brightspot variations (Armstrong et al. 2016) that can be well explained by 3D MHD simulations (Rogers 2017). Relevant planetary parameters are presented in Table 1. As discussed above, we equate the active layer's reference geopotential with a radiative equilibrium thermal energy reference level. Therefore the gravity wave speed is set using c g ≡ √ gH = RT eq = 3.0 × 10 3 m s −1 , where we use the planet's orbit-averaged effective temperature for the equilibrium reference temperature and the specific gas constant is calculated using the solar system abundances in Lodders (2010). We assume synchronous orbits, so Ω = 2π/t orbit , where t orbit is the orbital period. We calculate β ≡ 2Ω/R = 6.6 × 10 −13 m −1 s −1 , so the equatorial Rossby deformation radius is This is a fundamental length scale over which gravitational and rotational effects balance, and is the interaction length scale of planetary scale flows that corresponds to their latitudinal widths. The characteristic wave travel timescale, τ wave , is defined by the time a shallow-water gravity wave takes to travel over the distance L eq , and is τ wave ≡ L eq c g ≈ 2.2 × 10 4 s ≈ 0.26 Earth days, We set the reference thickness of the model's active layer to the atmospheric pressure scale height, that is H ≡ RT eq R 2 /GM = 4.3 × 10 5 m, where M is the planetary mass and G is Newton's gravitational constant.
In hydrodynamic shallow-water models (e.g., Shell & Held 2004;Langton & Laughlin 2007;Showman & Polvani 2010Showman et al. 2012;Perez-Becker & Showman 2013), the forcing profile is usually set so that ∆h eq /H ∼ (T day − T eq )/T eq , where T eq is the average reference temperature (for a given atmospheric depth) and T day is the maximal dayside reference temperature (at that atmospheric depth). For comparison, applying the reference temperatures used for HAT-P-7b in Rogers (2017), this equates to ∆h eq /H ∼ 0.22, 0.19, and 0.14 at P = 10 −3 bar, 10 −2 bar, and 10 −1 bar respectively. We consider models with ∆h eq /H = {0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6} to cover forcing parameter regimes within and either side of the expected range.
The simulations presented in this paper have a viscous diffusion of ν = 4 × 10 8 m 2 s −1 . In terms of "true" physical values, this diffusion coefficient is comparatively large; yet, upon checking, we find that viscous components of Equation (1) remain negligibly small. This is to be expected as we are predominantly modelling large scale planetary flows, upon which viscous dissipation generally has little direct influence. We set the magnetic diffusivity to η = 4 × 10 8 m 2 s −1 , which within the expected η range on HAT-P-7b's nightside. These values of η and ν are both small enough to make the dynamical timescales of our system much smaller than the diffusion timescales. In 3D geometries, longitudinal variations in η are likely to play an important role in the evolution of the magnetic field, but we defer considerations of this more complicated problem to a future paper.
The timescales τ rad and τ drag respectively determine the frequency over which Newtonian cooling and magnetic drag from the deep-seated (but not atmospheric) magnetic field are allowed to occur. Studies of hydrodynamic shallow-water analytics (Showman & Polvani 2011) and simulations (Perez-Becker & Showman 2013) show that varying τ rad controls the efficiency of (geopotential) energy redistribution occurs; whereas varying τ drag adjusts the distance over which atmospheric recirculation patterns can flow before becoming significantly damped. Hence, since τ rad and τ drag adjust qualitatively similar (albeit non-identical) fundamental flow features, it can be beneficial to reduce the modelling problem by fixing τ drag = τ rad . We do so in three of the examined cases of Section 3: (a) short τ rad and strong drag, τ rad = τ drag = τ wave ; (b) moderate τ rad and moderate drag, τ rad = τ drag = 5τ wave ; and (c) long τ rad radiative and weak drag, τ rad = τ drag = 25τ wave . However, as τ rad and τ drag are not necessarily equivalent in hot Jupiter atmospheres, we also consider the additional two cases: (d) short τ rad and weak drag, τ rad = τ wave , τ drag = 25τ wave ; and (e) long τ rad and strong drag, τ rad = 25τ wave , τ drag = τ wave . Rogers & Komacek (2014) found magnetically-driven reversals to occur in the upper atmospheres of ultra-hot Jupiters, where τ rad ∼ τ wave and τ drag ∼ τ wave , the conditions are most akin to case (a), though τ rad and τ drag are not generally exactly equal.
The remaining free parameter in our system is V A , which determines the magnitude of the system's magnetic field. Our general approach is to increase V A , from V A = 0, until we find a change in the nature of the SWMHD system (i.e., hotspot reversals). Here we highlight that, for large enough V A , we always find hotspot reversals in the SWMHD model, regardless of our choices of ∆h eq /H, τ rad , and τ drag . In Section 3 we will discuss both hydrodynamic and magnetohydrodynamic solutions over a wide range of parameter choices to illustrate the magnetic mechanism that drives reversals, and its robustness to changes in parameter space.

Model validity
Here we briefly discuss validity criteria for our model in the context of our parameter choices. First, we comment that H/L eq 6 × 10 −3 1, so the shallowwater approximation is well-founded and vertical dependences in the atmosphere are not of leading order importance. Secondly, we take Ω = Ω z, which is typically known as the traditional approximation and is formally valid in the limit of strongly stable stratification (N 2 /Ω 2 1, e.g., Vallis 2006). For our parameters, N 2 /Ω 2 ∼ 4 × 10 4 1, so this approximation is also well-founded. Thirdly, our y boundaries are located at y = ±Rπ/2 ∼ ±2.3L eq , so the impermeable wall at our model's "poles" has little physical influence on our solutions and does not interact with the equatorial dynamics we wish to study. Finally, the equatorial beta-plane truncation of the Coriolis parameter is f = βy + O((y/R) 3 ), so we are careful not to draw conclusions about the polar flows (with y R), where the Coriolis parameter is overestimated and boundary effects can occur. The initial magnetic field choice of Equation (10) is based on an analogous Taylor truncation, so places no further constraint on our discussion.

NUMERICAL SOLUTIONS
In this section we discuss numerical solutions of the model presented in Section 2. First, in Sections 3.1 and 3.2 we respectively highlight the basic flow behaviours of hydrodynamic and magnetohydrodynamic solutions. Then, in Section 3.3, we discuss detailed force balances of these numerical solutions. In Sections 3.1 to 3.3, we focus on solutions with ∆h eq /H = 0.2, which lies within the expected forcing range of our fiducial planet HAT-P-7b (see Section 2). Finally, in Sections 3.4 and 3.5, we discuss the extension of the developed theory to other forcing magnitudes and toroidal field profiles.
We visualise the basic form of our numerical solutions by plotting their (non-dimensionalised) geopotential distributions in Figure 1. As discussed in Section 2, we use geopotential energy, gh, as a shallow-water proxy of thermal energy so the geopotential distributions are analogous to those of temperature perturbations. In the hydrodynamic version of our shallow-water model, solutions are known to converge upon a steady state (e.g., Langton & Laughlin 2007;Showman & Polvani 2010;Showman et al. 2013;Perez-Becker & Showman 2013) and we replicate such hydrodynamic steady state solutions in the lefthand column of Figure 1 for comparison with our MHD simulations, which we plot in the middle and righthand columns for two difference solution phases (see Section 3.2). In each row of Figure 1 (from top to bottom) we display the solutions for (a) short τ rad and strong drag, τ rad = τ drag = τ wave ; (b) moderate τ rad and moderate drag, τ rad = τ drag = 5τ wave ; (c) long τ rad and weak drag, τ rad = τ drag = 25τ wave ; (d) short τ rad and weak drag τ rad = τ wave , τ drag = 25τ wave ; and (e) long τ rad and strong drag τ rad = 25τ wave , τ drag = τ wave .

Basic hydrodynamic solutions
Generally, in the hydrodynamic steady state solutions ( Figure 1, lefthand column) there are two dominant flow features. Drag-adjusted geostrophic circulations dominate at mid-to-high latitudes; while zonal jets dominate at the equator. The drag-adjusted geostrophic circulations satisfy a three-way force balance between horizontal pressure gradients, the Coriolis force, and Rayleigh drag (see Section 3.3). In the northern hemisphere, this balance is characterised by flows that circulate clockwise about the geopotential maximum and anticlockwise about the geopotential minimum; while the converse is true in the southern hemisphere. The dominant acceleration components in the equatorial regions are horizontal pressure gradients, which are largest in the zonal direction; the Rayleigh drag, which is simply a damping force that reduces wind speeds; and an advection correction, which is of lower order importance if drags are not weak (again, see Section 3.3). Hotspots are, by definition, located at the equatorial pressure maxima so the pressure driven zonally-directed equatorial jets diverge from them.
Newtonian cooling drives a solution's geopotential distribution towards the equilibrium geopotential (see that gh → gh eq as τ rad → 0). Therefore τ rad determines two things: how far planetary flows can redistribute geopotential energy before cooling occurs; and the magnitude of pressure gradients in the system, which in-turn determine planetary flows magnitudes (see Figure 1, lefthand column and axis scales). The Rayleigh drag reduces wind speeds everywhere. At equatorial latitudes, a strong Rayleigh drag decreases the distance that the zonal jets can redistribute geopo- Figure 1. The effect of azimuthal magnetic fields on energy redistribution. Contours of the relative layer thickness deviations (rescaled geopotential energy deviations) are plotted on colour axes that are shared along rows, with (individually-normalised) velocity vectors, hotspots (cyan crosses), and lines of constant A (solid/dashed for Bx positive/negative) over-plotted. In each column, reading from left to right, we present hydrodynamic steady state solutions (VA = 0), supercritical MHD solutions moments before reversal, and supercritical MHD solutions in the reversed quasi-steady phase. We present solutions in the following parameter regimes: (a) τ rad = τ drag = τwave, with VA = 0 or VA = 1.6cg in the top row; (b) τ rad = τ drag = 5τwave, with VA = 0 or VA = 0.7cg in the second row; (c) τ rad = τ drag = 25τwave, with VA = 0 or VA = 0.2cg in the third row; (d) τ rad = τwave, τ drag = 25τwave, with VA = 0 or VA = 1.4cg in the fourth row; (e) τ rad = 25τwave, τ drag = τwave, with VA = 0 or VA = 0.5cg in the bottom row. tential energy along the equator, increasing the relative severity of zonal geopotential gradients. At midto-high latitudes the Coriolis force becomes significant and solutions satisfy the aforementioned drag-adjusted geostrophic balance. In a "true" geostrophic balance, without suppression from drags and forcing, pressure gradients are exactly balanced by the Coriolis force, which acts perpendicularly to the velocity causing flows to rotate (to their right in the northern hemisphere and to their left in the southern hemisphere). This yields large-scale mid-to-high latitude vortices that are aligned with isobars, similar to those seen in the short τ rad , weak drag, hydrodynamic solution ( Figure 1 (c), lefthand column). However, the slowing of winds from the Rayleigh drag reduces the magnitude of Coriolis deflection. Therefore in the strong drag limit largescale vortices cannot fully develop. Similarly, when τ rad is short, heating/cooling occurs before large-scale vortices fully develop. Comparing the mid-to-high latitude flows of the hydrodynamic solutions, one finds a transition between the long-τ rad /weak-drag solutions, with fully-formed geostrophic vortices, to the short-τ rad /strong-drag solutions, in which the drag-adjusted geostrophic circulations are approximately aligned with the isobars of the equilibrium geopotential (see Figure 1, lefthand column). Aside from an unphysical special case discussed in Showman & Polvani (2011) and Perez-Becker & Showman (2013), for all finite physicallyrelevant choices of τ rad and τ drag , the meridional mass transport into the equator, caused by the drag-adjusted geostrophic circulations, is maximised east of the substellar point.
These solutions always exhibit eastward hotspots. This is because the equatorward (rescaled) geopotential energy transport from the mid-to-high latitude circulations, −∂(hv)/∂y, always has its equatorial maximum located eastward of the substellar point. At the equator, the pressure gradient drives winds that diverge from hotspots, causing equatorial geopotential energy transport away from the hotspot regions (i.e., −∂(hu)/∂x < 0 in hotspot regions). Hence, by Equation (2) (geopotential energy conservation), the hotspots locate themselves at the equatorial point of maximal incoming geopotential energy flux, which is located between the equatorial maxima of −∂(hv)/∂y and Q. The Newtonian cooling (Q) attempts to return a solution to its forcing equilibrium (i.e., with its hotspot at the substellar point); whereas, as stated above, the equatorial maximum of −∂(hv)/∂y is always eastward. The degree of the hotspot's eastward offset is therefore determined by the location of the equatorial maximum of −∂(hv)/∂y and its relative magnitude compared to Q. In short, the size of the (eastward) hotspot offset is determined by the efficiency over which the drag-adjusted geostrophic circulations can redistribute thermal 5 energy from the western equatorial dayside to the eastern equatorial dayside, by circulating it to-and-from the higher latitudes.

Basic magnetohydrodynamic solutions
In the weakly-magnetic limit, shallow-water magnetohydrodynamic solutions behave much like their hydrodynamic counterparts (i.e., solutions reach a steady state that is characterised by eastward hotspots, zonal equatorial winds, and drag-adjusted geostrophic circulations at mid-to-high latitudes). However, when the azimuthal magnetic field exceeds a critical magnitude the nature of the solution changes. Supercritical magnetic solutions have three phases: an initial phase, in which winds and geopotentials resemble their hydrodynamic counterparts but their circulations induce magnetic field evolution; a transient phase, in which mid-to-high latitude winds align with the azimuthal magnetic field and dayside equatorial winds experience a net westward acceleration, driving an east-to-west hotspot transition; and a reversed quasi-steady phase, in which westward zonally-dominated dayside winds maintain westward hotspots (until, after a comparably long period of time, the magnetic field decays via magnetic diffusion). 6 We present geopotential distributions of supercritical magnetic solutions in the transient and quasisteady phases in the two righthand columns of Figure 1 (middle and right respectively). The supercritical magnetic solutions are plotted for the same drag choices as the hydrodynamic solutions that they share a row with (see Section 3.1), but now lines of constant A, which approximately correspond to field lines of the horizontal magnetic field, are also over-plotted for visualisation of the magnetic field.
After a magnetic solution's initial phase, in which it behaves similarly to its hydrodynamic counterpart, in mid-to-high latitude regions there is a competition between the drag-adjusted geostrophic balance and the magnetic tension (i.e., B · ∇B, the restorative force that acts to straighten bent horizontal magnetic field lines) that the circulating flows generate. Initially, the magnetic field is purely azimuthal, with only latitudinal gradients in its profile, so magnetic tension is zero everywhere. To understand the magnetic field's evolution we highlight that, as the magnetic diffusion timescale is large in comparison to the dynamical timescales of the system, A is approximately materially conserved. This means that lines of constant A are advected by the mid-to-high latitude circulations, bending them and causing a growth of magnetic tension. For subcritical magnetic field strengths, a drag-adjusted magnetogeostrophic balance can be supported, with winds and geopotential profiles making small adjustments to balance the magnetic contribution (before magnetic diffusion eventually returns the system to a hydrodynamic steady state). In contrast, for supercritical magnetic field strengths, magnetic tension becomes strong enough to obstruct the drag-adjusted geostrophic circulations and solutions enter into a transient phase, which ultimately results in hotspot reversals. In Section 3.3, we shall see that the reversal is driven by a westward Lorentz force acceleration in the region surrounding the hotspot, which is itself generated by this obstruction of geostrophic balance. The westward Lorentz force acceleration causes the point of zonal wind divergence on the equator to shift eastwards, so that in hotspot regions geopotential energy flux is westward (i.e., ghu < 0) rather than zero. This shifts the hotspot westward until the system rebalances into a state with a westward hotspot (again, see Section 3.3).
We find that this reversal mechanism (i.e., westward equatorial-dayside Lorentz force accelerations driven by the obstruction of geostrophic balance) always leads to hotspot reversals in the SWMHD model, regardless of our choice of ∆h eq /H, τ rad , and τ drag . However, since these parameters control pressure gradient magnitudes and recirculation efficiency, they determine the critical magnetic field strength sufficient for reversal. We present bounds on the magnetic field strength's critical magnitude, V A,crit , for various parameter choices in Figure 5. Generally, ∆h eq /H and τ rad set the magnitude of a solution's pressure gradients, and therefore the magnitude of the circulations to be overcome, so shorter τ rad and larger ∆h eq /H correspond to larger V A,crit magnitudes. Initially in long τ drag solutions the fully formed large scale geostrophic vortices advect the lines of constant A efficiently until they are resisted by magnetic tension; whereas, for short τ drag solutions, the slowing of winds from drags decreases the distance over which winds initially advect the lines of constant A. Therefore weak drag solutions generally experience a larger degree of field line bending and hence more magnetic tension (relative to the other accelerations in their solutions for a given V A ) than strong drag solutions. Put simply, strong drag solutions require larger V A,crit magnitude to reverse. We quantify dependences of V A,crit on ∆h eq /H, τ rad , and τ drag in later discussion.
In the quasi-steady phase of supercritical SWMHD solutions, the magnitudes of τ rad and τ drag determine the efficiency of the westward energy redistribution. For large τ rad and τ drag timescales, the (westward) hotspot offsets are large as the equatorial pressure-Lorentz balance is free to redistribute energy towards the point where the zonal winds converge, almost entirely without restriction; Conversely, for short τ rad and τ drag timescales, this equatorial energy redistribution is less efficient and hotspot offsets are smaller. Comparing between rows in Figure 1 (righthand column), suggests τ drag is the most influential timescale in determining westward hotspot offsets in the SWMHD system.

Force balances
In this subsection we compare the force balances of Equation (1) for hydrodynamic and supercritical MHD solutions with the parameters of regime (b) in Figure 1 (i.e., for ∆h eq /H = 0.2, τ rad = τ drag = 5τ wave , with either V A = 0 or V A = 0.7c g ). We highlight how the pres-ence of a strong equatorially-antisymmetric azimuthal magnetic field modifies the force balances of different planetary regions, and link these modifications to the more general discussions of Sections 3.1 and 3.2.
In Figures 2 and 3 we respectively plot the dominant meridional and zonal acceleration components of Equation (1), for solutions in regime (b). In the lefthand column of Figures 2 and 3, we present the acceleration components for the hydrodynamic steady state solution; whereas in the middle and righthand columns of Figures 2 and 3, we present the acceleration components of the transient and quasi-steady phases of its supercritical MHD counterpart. Along each row of Figures 2 (meridional components) and 3 (zonal components), we plot (from top downwards) the acceleration contributions due to horizontal pressure gradients (−g∇h), the Coriolis effect (−f z × u), the Lorentz force (B · ∇B), Rayleigh drag (−u/τ drag ), and advection (−u·∇u). Additionally, in the bottom row of Figure 2 we plot the total meridional acceleration (∂v/∂t) and, likewise, in the bottom row of Figure 3 we plot the total zonal acceleration (∂u/∂t). For the presented parameter choices the acceleration contributions due vertical mass transport (R) and viscous diffusion (D ν ) are much weaker so are not included in the plots.
At mid-to-high latitudes, the force balances of hydrodynamic solutions in steady state are well described by the three-way drag-adjusted geostrophic balance discussed in Section 3.1. In particular, Figures 2 and 3 (lefthand column) highlight this for regime (b), showing that in both horizontal directions the mid-to-high latitude accelerations due to horizontal pressure gradients and the Coriolis force almost exactly cancel, albeit with small Rayleigh drag adjustment and a yet smaller advection contribution. The meridional components of these accelerations remain balanced in equatorial regions, with all of them vanishing at the equator. However, in the zonal direction, the Coriolis force vanishes in equatorial regions but zonally-directed pressure gradients do not, so zonal pressure gradients are balanced by the Rayleigh drag, with an advection adjustment. Since hotspots in hydrodynamic solutions are always located where zonal equatorial jets diverge, these three acceleration components are equally zero at hotspots (see cyan markers in Figure 3). As discussed in Section 3.1, hotspots are driven eastward by the net west-to-east equatorial energy transfer that results from the mid-to-high latitude drag-adjusted geostrophic circulations.
As discussed in Section 3.2, magnetic tension (B·∇B) is initially zero everywhere so MHD solutions initially resemble their hydrodynamic counterparts. However, lines of constant A (which closely follow magnetic field lines) Figure 2. Meridional force balances. In each column, reading from left to right, we plot meridional accelerations corresponding to hydrodynamic steady state solutions, transient phase supercritical MHD solutions, and quasi-steady supercritical MHD solutions. In rows one to four, we respectively plot meridional accelerations due to horizontal pressure gradients, the Coriolis effect, the Lorentz force, and Rayleigh drag; the summed meridional accelerations are plotted in row five. The solutions are presented for ∆heq/H = 0.2, τ rad = τ drag = 5τwave, with VA = 0 (HD) or VA = 0.7cg (MHD) (i.e., parameter regime (b) in Figure 1). are advected by the mid-to-high latitude circulations that are archetypal of hydrodynamic solutions. This causes them to bend equatorward between the western and eastern dayside (where the initial circulations are poleward and equatorward respectively; see Figure 1, row (b), middle column). Consequently, a restorative Lorentz force that resists meridional winds is produced (see Figure 2, third row, middle column). For subcritical MHD solutions (not plotted) this Lorentz force resists but does not fully obstruct the mid-to-high latitude circulations, which adjust into a (drag-adjusted) magneto-geostrophic balance. However, in supercritical MHD solutions, the Lorentz force resists meridional winds strongly enough to zonally-align the mid-to-high latitude winds. Hence, supercritical MHD solutions enter into the transient phase discussed in Section 3.2.
When the magnetic field geometry is azimuthally dominated, understanding Lorentz force accelerations is less intuitive in the zonal direction than in the meridional direction (in which they simply oppose meridional flows). The zonal Lorentz accelerations, B · ∇B x , are most easily understood geometrically when considered as the directional derivative of B x along horizontal magnetic field lines, which are approximately equivalent to lines of constant A. When the magnetic field lines bend equatorward they generally move into regions of  Figure 2, we present solutions for the parameter choices ∆heq/H = 0.2, τ rad = τ drag = 5τwave, with VA = 0 (HD) or VA = 0.7cg (MHD) (i.e., parameter regime (b) in Figure 1). To aid discussion in the text, hotspot locations have been marked with cyan crosses in hydrodynamic solution panels that correspond to zonal acceleration components with a non-zero equatorial contribution.
smaller |B x |, hence the zonal Lorentz force component generally accelerates flows westward (as B · ∇B x < 0); conversely, when they bend poleward they generally move into regions of larger B x , hence the zonal Lorentz force component generally accelerates flows eastward (as B · ∇B x > 0). One can see this by comparing lines of constant A in mid-to-high latitudes of Figure 1 (row (b), middle column) with the corresponding mid-to-high latitude zonal Lorentz force accelerations in Figure 3 (third row, middle column). Since magnetic field lines bend equatorward between the western and eastern dayside at mid-to-high latitudes, the Lorentz force accelerates mid-to-high latitude dayside flows westward (and eastward on the nightside).
Similar westward dayside Lorentz force accelerations are generated along the equator by magnetic field lines bending into equatorial regions. To visualise this, in Figure 4 we plot the horizontal magnetic field geometry (top row) and the zonal component of the Lorentz force (bottom row) in the equatorial region, −π/8 < y/R < π/8, for the transient (left) and quasi-steady (right) phases of the supercritical MHD solution (again, for parameter regime (b)). In the initial phase the Lorentz force primarily acts to resist drag-adjusted geostrophic circulations (see above). Therefore, in the early transient phase the magnetic field lines bend equatorward between the western and eastern dayside (where the initial circulations are poleward and equatorward respectively). For  Figures 2 and 3, we present the transient (lefthand column) and quasi-steady (righthand column) phases of the supercritical MHD solution with ∆heq/H = 0.2, τ rad = τ drag = 5τwave, and VA = 0.7cg (i.e., parameter regime (b) in Figure 1), though we restrict this plot to the equatorial region, −π/8 < y/R < π/8. The zonal Lorentz force acceleration is the directional derivative of Bx along horizontal magnetic field lines (approximately lines of constant A, see Section 2.3).
the lowest equatorial regions (|y/R| π/32 in Figure 4) such equatorward magnetic field line bending causes the lines to move into regions smaller |B x |. Consequently, zonal Lorentz force accelerations are westward in regions surrounding the hotspot (see Figure 4, lefthand column). In fact, zonal Lorentz force accelerations are always westward in hotspot regions, regardless of radiative/drag/forcing parameter choices, because in hydrodynamic (and weak/early-phase MHD) solutions hotspots are located between the substellar point and the (eastward) maximum of equatorward flow (where lines of constant A are bent most equatorward). The resulting westward accelerations cause an equatorial imbalance in the zonal momentum equation (see Figure 3, bottom row, middle column), which drives the point of zonal equatorial wind divergence eastwards of the hotspot and, consequently, shifts the hotspot westward (see discussion in Section 3.2). Finally, as these westward accelerations cause dayside equatorial winds to become more westward, lines of constant A are swept from east to west along the equator, bending them further and thus enhancing equatorial Lorentz force accelera-tions across all equatorial latitudes (see Figure 4, righthand column) 7 .
Across radiative/drag/forcing parameter choices, when the hotspots have transitioned westwards the system rebalances into a quasi-steady state, which is characterised by westward hotspots, zonally-aligned winds, and magnetic field lines that have an equatorward bend along the line x = 0 in equatorial regions. The predominant meridional balance is between pressure gradients, the Coriolis force, and the Lorentz force (see Figure 2, righthand column); whereas the predominant zonal balance is between pressure gradients, the Lorentz force, and the Rayleigh drag (see Figure 3, righthand column). In these balances the zonally-aligned winds cause the meridional Rayleigh drag and the zonal Coriolis force to be small. We comment that as the magnetic field eventually diffuses away, the balance adjusts to the decreasing Lorentz force contribution, eventually restor-ing the drag-adjusted geostrophic/magneto-geostrophic balances associated with hydrodynamic and weaklymagnetic solutions (and hence eastward hotspots).

Forcing dependence
We find that, when one compares marginally supercritical magnetic solutions with ∆h eq /H = {0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6}, the qualitative physical behaviours and balances discussed in Sections 3.1 to 3.3 (and illustrated in Figures 1 to 4) remain highly similar (in fact, remarkably so). The only discernible changes we observe between marginally supercritical magnetic solutions, upon increasing ∆h eq /H, are an approximately linear scaling of dependent variable magnitudes and a correction from advection, which generally only provides a lower order correction. This is to be expected from the theory we have developed so far, as the process that needs to be overcome in order to trigger hotspot reversals (i.e., the drag-adjusted geostrophic balance) is a linear one. Consequently, choices of ∆h eq /H do not change the mechanics of the hotspot reversals, though they do determine quantitive features of the system (such as magnitudes and V A,crit ).
We can use our developed understanding of the reversal mechanism to predict magnitudes V A,crit , with simple scaling arguments based on the respective magnitudes of geostrophic circulations and the restorative Lorentz force. Let τ −1 geo = U/L eq be the frequency over which geostrophic flows circulate and τ −1 A = V/L A be the (Alfvén) frequency over which the azimuthal field attempts to zonally-align these circulations, where U, V, L eq and L A are the typical velocity and length scales associated with the two opposing processes. Reversals occur when τ −1 A τ −1 geo or equivalently when V UL A /L eq (i.e., when the azimuthal field is strong enough to restrict the geostrophic flows). Perez-Becker & Showman (2013) showed the velocities of geostrophic circulations in Coriolis dominated regions scale like highlighting that the reversal threshold is expected to have a linear dependence on ∆h eq /H. In Figure 5 we plot the dependence of V A,crit on ∆h eq /H from our simulations. For comparison, we overplot the lines V A,crit c g = 2πR(∆h eq /H) κL eq (τ rad /τ wave ) where, since the circulations bend field lines on the planetary azimuthal scale, we take L A = 2πR and κ is a constant of order unity based on the profile of B 0 (y). 8 We generally find reasonable agreement between this simple scaling prediction and numerical simulations, particularly in the realistic regimes of τ rad short and ∆h eq /H ∼ 0.1-0.3, but note that V A,crit approaches a minimum as ∆h eq /H → 0, which we shall consider in Section 4. This scaling law approximation deals less favourably in the (less physical) long τ rad cases, where τ drag dependencies become important. However, as we shall discuss in Sections 5 and 6, the other uncertainties in atmospheric characteristics are likely to provide much larger uncertainties than those arising from this scaling law approximation.

Linear-Gaussian magnetic field profiles
Upon comparing the discussed results to their equivalents for the initial magnetic field profile B(x, 0) = V A (y/L eq ) exp(1/2 − y 2 /2L 2 eq ), we found the same mechanical features. Namely, subcritical solutions behave similarly to their hydrodynamic counterparts; whereas, for supercritical magnetic solutions, the obstruction of geostrophic circulations by the magnetic field causes zonal wind alignment, a westward Lorentz force acceleration, and therefore reversed hotspots. The only different qualitative flow features arise at the poles, where V A (y/L eq ) exp(1/2−y 2 /2L 2 eq ) decays, but our model and aims are not directed towards the polar regions. The quantitative differences between solutions are also tend to be minor, with a second order change in V A,crit as the two profiles cause a slightly different magnitude of Lorentz force to be generated for a given V A . To make this comparison, we have marked V A,crit for Linear-Gaussian profiles on Figure 5 with starred markers. We conclude that the choice of a B 0 ∝ tanh(y/L eq ) profile is a useful simplification when considering reversals. This can be advantageous due to properties of the hyperbolic tangent function, which is both monotonic and bounded as y → ∞.

Summary of findings
In this section we have identified the mechanism responsible for driving hotspot reversals in our SWMHD model. The reversals are caused by the westward Lorentz force acceleration that is generated when strong equatorially-antisymmetric azimuthal magnetic fields 8 κ is an estimate of the relative strength of B 0 (compared to V A ) at low latitudes, y 0 , where westward Lorentz force accelerations first develop. In Figure 5, we take κ = e 1/2 tanh(y 0 /Leq) ≈ 0.47 (using y 0 ≈ Rπ/16 based on Figure 4). B0,max = |ωR,1|/k B0,max = |ωR,3|/k B0,max = |ωR,5|/k τ rad = τ wave , τ drag = τ wave τ rad = 5τ wave , τ drag = 5τ wave τ rad = 25τ wave , τ drag = 25τ wave τ rad = τ wave , τ drag = 25τ wave τ rad = 25τ wave , τ drag = τ wave B 0 = V A e 1/2 tanh(y/L eq ) B * 0 = V A (y/L eq )e 1/2−y 2 /2L 2 eq Figure 5. Quantitive dependencies of critical magnetic field amplitudes on the forcing magnitude parameter, ∆heq/H, for different choices of τ rad and τ drag . Critical magnetic field amplitudes are illustrated with marker points. These are located mid-way between the upper/lower bounds of the identified critical amplitude range, for a particular parameter set, with error bars indicating these upper/lower bounds. Lines indicating scaling law predictions (dashed) and zero-amplitude limits based on the linear theory (dotted; see Section 4) are overlaid obstruct the geostrophic circulation patterns responsible for energy redistribution in the hydrodynamic system. The understanding we have developed explains why such hotspot reversals always emerge in the SWMHD model, regardless of our choices for the free forcing/drag parameters ∆h eq /H, τ rad , and τ drag . Moreover, this developed understanding has allowed us to use simple scaling arguments to predict the reversal threshold, V A,crit , in terms of planetary parameters, finding reasonable agreement between predictions and numerical simulations in realistic forcing regimes for our fiducial planet HAT-P-7b. However, our simulations also show that V A,crit approaches a minimal threshold in the zero amplitude limit. In Section 4 we shall probe linear theory to explain this finding. For this, we shall use our finding that, when compared, equatorially-antisymmetric azimuthal magnetic field profiles with similar latitudinal dependence at equatorial and mid-latitudes behave similarly to one another.

Linearised steady state solutions
First we seek to establish the features of the reversals that linear theory can capture, and its limitations. We do so by linearising the non-diffusive versions of Equations (1) to (4) about the background state {u 0 , v 0 , h 0 , A 0 } = {u 0 (y), 0, H, A 0 (y)}, where H is the (constant) background layer thickness, A 0 is defined such that dA 0 /dy = HB 0 for the latitudinallydependent azimuthal background magnetic field, B 0 = B 0 (y) x, and u 0 (y) is to be fixed in a manner that balances the zeroth order zonal momentum equation of the hydrodynamic version of the system which we wish to investigate. To probe the system at the reversal threshold, we assume steady state perturbations exist about this background state and apply the plane wave ansatz, {u 1 , v 1 , h 1 , A 1 } = {û(y),v(y),ĥ(y),Â(y)}e ikx , where k denotes the azimuthal wavenumber and subscripts of unity denote perturbations from the background state. Such perturbations satisfy is the first order forcing contribution in the system based on the equilibrium thickness profile, h eq = H+∆h eq cos (kx) exp(−y 2 /2L 2 eq ), and based on our numerical findings we have assumed that R does not make a first order contribution to Equations (15) and (16). Before solving, we note that hydrodynamic solutions are never singular, but that Equation (18) causes the magnetic version of the system to be singular if u 0 = 0. To compare to the simulations of Section 3, we solve the system for f = βy and B 0 = V A e 1/2 tanh(y/L eq ).
For a given u 0 (y), we seek solutions of Equations (15) to (18) on −L y < y < L y , with impermeable boundaries at y = ±L y , using the shooting method outlined in Appendix A. We take L y = 5L eq (see Equation (11)), which is large enough to ensure that the outer boundary condition has a negligible influence on solutions. We solve the system for u 0 (y) = U 0 exp(−y 2 /2L 2 eq ), where U 0 is chosen so that in the hydrodynamic limit the zonallyaveraged zonal-acceleration in Equation (22) of Showman & Polvani (2011) vanishes at the equator. We plot linear solutions for ∆h eq /H = 0.01 in Figure 6, on the reduced domain −π/2 < y/R < π/2, for three τ rad and τ drag choices, comparing hydrodynamic solutions with MHD solutions at the threshold of criticality, as found by simulations.
Hydrodynamic solutions generally resemble those discussed in Showman & Polvani (2011), albeit with an adjustment due to u 0 (as discussed by Tsai et al. 2014, for u 0 constant). They are characterised by geostrophic circulations at mid-to-high latitudes and zonal pressure driven jets at equatorial latitudes. Such solutions closely resemble the non-linear hydrodynamic steady state solutions we discussed in Section 3. The characteristic flow patterns of hydrodynamic steady state solutions can also be directly linked to the forcing responses of specific standing, planetary scale, equatorial shallowwater waves (Matsuno 1966;Showman & Polvani 2011;Tsai et al. 2014). The geostrophic circulations are linked to the planetary scale equatorial Rossby waves, which are geostrophic in nature at mid-to-high latitudes; while the equatorial jets are linked the superposition of the planetary scale equatorial Rossby waves and the equatorial Kelvin wave, which travels eastward about the equator in response to pressure perturbations. The pre-sented linear hydrodynamic solutions all have eastward hotspots (located at points of zonal wind divergence), as the linearised meridional convergence of geopotential flux into the equator, −gH∂v 1 /∂y| y=0 , is maximised eastward of the substellar point (due to the form of the geostrophic circulations; further discussion in Section 3.1).
The marginally critical MHD solutions share some common characteristics with their non-linear simulated counterparts. Specifically, in these solutions the aligning influence of the meridional Lorentz force is strong enough obstruct geostrophic circulations, which are replaced by zonally-aligned winds. However, unlike their simulated non-linear counterparts, the magnetohydrodynamic solutions do not have westward hotspots. The arises because in this simple linear model one can show that the Lorentz force components, which drive hotspots reversals in non-linear simulations (see Section 3), vanish at the equator. 9 Instead, marginally critical MHD solutions approach a limit of zero hotspot offset, as the obstruction of geostrophic circulations causes −gH∂v 1 /∂y| y=0 → 0. This highlights that in simple linear models, with similar linearisations of the Lorentz force and the induction equation (i.e., without more sophisticated treatments of magnetic diffusion and non-linear effects), one can identify the obstruction of geostrophic circulations that cause hotspot reversals in non-linear simulations, but not westward hotspot offsets explicitly. This observation is useful in the remainder of this section, where we aim to link the magnetic obstruction of geostrophic circulation patterns to wave dynamics.

Wave background: Alfvén-Rossby wave coupling
Various authors have studied the linear waves present in rotating MHD systems. Early studies, which used quite general (usually uniform) flow/field geometries, focussed on the influence that these waves have on the geodynamo (Hide 1966(Hide , 1969aAcheson & Hide 1973). Since the development of SWMHD (Gilman 2000), authors have been able to utilize its reduced geometry to study waves in more specific flow/field geometries. Rotating SWMHD waves have now been studied for a variety of thin-layered astrophysical and geophysical systems including the geodynamo, the solar tachocline,   Figure 5). Solutions are plotted for (a) τ rad = τ drag = τwave (left); (b) τ rad = τ drag = 5τwave (middle); and (c) τ rad = τ drag = 25τwave (right). Solutions are calculated for −5Leq < y < 5Leq, but are cut off for −Rπ/2 < y < Rπ/2 (recall, Leq/R ≈ 0.67). The strong magnetic field aligns flows preventing geopotential recirculation between latitudes, but in the linearised model the (non-linear) equatorial Lorentz force is zero. Consequently, in the linearised model hotspot offsets of marginally critical MHD solutions tend to zero, but do not reverse like full SWMHD simulations. and neutron star atmospheres (Schecter et al. 2001;Zaqarashvili et al. 2007Zaqarashvili et al. , 2009Heng & Spitkovsky 2009;Márquez-Artavia et al. 2017;Zaqarashvili 2018). Of these, the solar tachocline, which is also expected to have an equatorially-antisymmetric toroidal dominant magnetic field geometry, can be considered as similar to the hot Jupiter system. Schecter et al. (2001) studied waves in the local regions of the solar tachocline, focusing on regions away from the equator; whereas Zaqarashvili et al. (2007,2009) studied the global dynamics of these waves for the two extreme cases ( 1 and 1) of the rotation-stratification parameter, = 4Ω 2 R 2 /c 2 g . However, the atmosphere of HAT-P-7b lies in the region of parameter space between these two extremes ( ≈ 4.8). Zaqarashvili (2018) studied equatorial SWMHD waves using an equatorial betaplane model for two purely-azimuthal magnetic field geometries: firstly, uniform and, secondly, equatoriallyantisymmetric (latitudinally-linear). London (2017) and London (2018) studied some asymptotic solutions of the beta-plane and spherical version of the system, with an equatorially-antisymmetric azimuthal field, in certain weak and strong field limits, but we wish to study the transition where magnetism becomes dynamically important. The predictions made in Hindle et al. (2019) were based on the equatorially-antisymmetric azimuthal magnetic field results of Zaqarashvili (2018). However, in this section we relax the weakly-magnetic assumptions that those analyses take.
Past works of linear waves in rotating MHD systems find that Alfvén waves and Rossby waves are coupled. To illustrate this, we highlight the known local dispersion relations of waves in the non-diffusive, unforced, drag-free SWMHD system, with the uniform azimuthal background magnetic field, B 0 = V A,0 x, and the generalised beta-plane treatment f = f 0 + βy. 10 In local regions (i.e., if |y/R| 1 and |βy| |f 0 |), this linearised system associated with this background state may be approximately solved with the plane wave ansatz: {u 1 , v 1 , h 1 , A 1 } = {û,v,ĥ,Â}e i(kx+ly−ωt) , where hatted variables are constant amplitudes of the plane 10 Here f 0 ≡ f (y 0 ) = 2Ω sin θ 0 and β ≡ df /dy|y=y 0 = (2Ω/R) cos θ 0 are respectively the Coriolis parameter and the Coriolis parameter's local latitudinal variation at a reference latitude, θ 0 = y 0 /R, about which the system is centred.
wave solutions, k is the azimuthal wavenumber, l is the latitudinal wavenumber, and ω is the oscillation frequency. Seeking solutions that are first order in the Coriolis parameter only, yields the following dispersion relation (Zaqarashvili et al. 2007;Heng & Spitkovsky 2009): for arbitrary wave amplitudes, where K ≡ (k 2 + l 2 ) 1/2 . For V A,0 > 0, Equation (19) has four solutions, which in the rotation-free limit (f 0 = β = 0) are (Schecter et al. 2001) where c m,0 = (c 2 g + V 2 A,0 ) 1/2 is the magneto-gravity wave speed for a constant background magnetic field. The first pair of solutions are Alfvén waves, which are driven by magnetic tension and travel parallel to the background magnetic field; the second pair of solutions are magneto-gravity waves, which propagate horizontally to restore pressure gradients and magnetic tension. For rotationally modified waves, we follow Schecter et al. (2001) by labelling the rotationally modified Alfvén waves as slow "Alfvén branch" solutions and the rotationally modified magneto-gravity waves as fast "magneto-gravity branch" solutions.
In the fast wave limit (|ω|/|2Ω| 1), to leading order, fast magneto-gravity branch solutions satisfy (Heng & Spitkovsky 2009, but These two magneto-gravity branch solutions travel in opposite directions in order to restore pressure gradients and magnetic tension, but with a Coriolis modification. Solutions of this kind are known as magneto-Poincaré waves (as they reduce to Poincaré waves for V A,0 = 0; e.g., Heng & Spitkovsky 2009) or magneto-inertial gravity waves (and inertial gravity waves in hydrodynamics; e.g., Márquez-Artavia et al. 2017;Zaqarashvili 2018). We choose the inertial gravity and magnetoinertial gravity nomenclature (IG and MIG hereafter). The independence of Equation (21) on β highlights that IG/MIG solutions do not generally have a leading order dependence on β (and exist on the f-plane; e.g., Vallis 2006). In the slow wave (|ω|/|2Ω| 1) limit, the dispersion relation evaluated at the equator (f 0 = 0 and β = 2Ω/R), yields the two Alfvén branch solutions (Heng & Spitkovsky 2009, but where M = 4k 2 V 2 A,0 (c 2 m,0 k 2 +c 2 g l 2 )((c 2 m,0 +V 2 A,0 )k 2 +c 2 g l 2 ) is the magnetic component of the numerator. In the limit where these solutions are dominated by the Alfvén speed, these waves Alfvénic in nature (see by taking V A,0 dominatingly large). Conversely, in the hydrodynamic limit (V A,0 = 0) the Alfvén branch solutions reduce to so the eastward Alfvén branch solution vanishes and the westward Alfvén branch solution reduces to a Rossby wave.
This Alfvén-Rossby wave coupling is a well-documented feature of MHD in systems with a latitudinally dependent planetary vorticity (Hide 1966(Hide , 1969aAcheson & Hide 1973). However, Alfvén and Rossby waves are fundamentally different in nature. Rossby waves arise due to potential vorticity conservation and the latitudinal variation of the Coriolis parameter. They behave geostrophically and are highly dispersive, so can transfer energy and angular momentum to the surrounding system (e.g., Vallis 2006;Pedlosky 2013). Conversely, pure Alfvén waves travel parallel to the dominant azimuthal magnetic field geometry and are non-dispersive, so travel as zonally-aligned solitons. Comparing the oscillation frequency of Rossby (ω R ) and Alfvén (ω A ) waves gives suggesting that, for given choices of β and V A,0 , Rossby wave characteristics dominate at large scales; whereas Alfvén wave characteristics dominate at small scales. In Section 3, we showed that reversals on hot Jupiters are closely tied to the zonal-alignment of equatorially-adjacent geostrophic circulations by equatorially-antisymmetric azimuthal magnetic fields. Therefore, to investigate reversals in the zero amplitude limit, we examine the behaviour of equatorial waves as the Alfvén oscillation frequency approaches ω R in magnitude for an antisymmetric azimuthal background magnetic field.
Applying the plane wave ansatz, {u 1 , v 1 , h 1 , A 1 } = {û(y),v(y),ĥ(y),Â(y)}e i(kx−ωt) , the evolution of the perturbations is determined by the following linearised SWMHD system: From this we eliminateû,ĥ,Â,B x , andB y to obtain the single ordinary differential equation: for the latitudinal solving domain, −L y < y < L y , with where c m (y) ≡ (c 2 g + B 2 0 ) 1/2 denotes the (rotationless) magneto-gravity wave speed. This system can the contain singular points at y = y s , if ω = ±B 0 (y s )k (Alfvén singularity) or ω = ±c m (y s )k (magneto-gravity singularity), which we label based on the ω-regions each singularity is associated with.
If one attempts to write L in Sturm-Liouville form 11 , through use of an integrating factor, it is found that the highest order functional coefficient of the Sturm-Liouville operator, p = (ω 2 − B 2 0 k 2 )/(ω 2 − c 2 m k 2 ), is not independent of the oscillation frequency. Therefore, the desirable properties of the Sturm-Liouville eigenvalue problem (e.g., real eigenvalues and orthogonality of eigenfunctions) are not generally guaranteed. Zaqarashvili (2018) studied this system in the weaklymagnetic limit where singular points do not influence the planetary scale waves. 12 In this approximation L can be re-expressed in terms of the parabolic cylinder Sturm-Liouville operator (the hydrodynamic version of L, see Matsuno 1966). Therefore, away from singular ω-regions, where the approximations of Zaqarashvili (2018) hold, one may expect solutions to conform to Sturm-Liouville properties (which we find in the following analysis).

Equatorial wave solving method
We now examine non-trivial eigenvalue-eigenfunction pairs, {ω,v(y)}, that satisfy L{v} = 0 everywhere in the latitudinal domain, −L y < y < L y , subject to impermeable boundary conditions (i.e.,v(±L y ) = 0). We use the planetary parameters discussed in Section 2, f = βy and B 0 = V A e 1/2 tanh(y/L eq ). This B 0 (y) choice is useful because it is both monotonic and bounded as y → ∞ (London 2017), so there is at most one Alfvén singularity in each hemisphere. For this B 0 (y) choice, solutions with c g k ≤ |ω| ≤ (c 2 g + V 2 A e) 1/2 k have magneto-gravity singularities; while solutions with |ω| ≤ V A e 1/2 k have Alfvén singularities. We seek wave-like solutions with the planetary scale azimuthal wavenumber, k = 1/R. We find that solving this eigenvalue problem, without further approximation on L, is an analytically intractable problem so we use a semi-analytic approach.
Since L is symmetric about the equator, homogeneous solutions will be either symmetric (v symmetric and u,ĥ,Â antisymmetric) or antisymmetric (v antisymmetric andû,ĥ,Â symmetric) about the equator. 13 Although the system we solve here is unforced, we wish to compare solutions to the numerical simulations of Section 3, which had equatorially-symmetric forcing on h. Therefore, we only consider antisymmetric homogeneous solutions and solve L{v} = 0 in the upper-half domain, 0 < y < L y , with the antisymmetric lower boundary conditionv(0) = 0, which replacesv(−L y ) = 0. Eigenfunctions are defined up to a constant factor, so a third and final normalisation boundary condition must also be included. We set dv/dy| y=0 = N , where N is a normalisation constant chosen for numerical convenience, and take L y = 5L eq to ensure boundary influences are negligible.
We use a shooting method to seek eigensolutions. The shooting method calculates successive "shots" (or test solutions,v T ) for given test frequencies, ω T , where each shot satisfies L{v T } = 0, subject to two of the three boundary conditions. The third boundary condition is then satisfied by varying ω T so that the deviation from the third boundary condition, G[ω T ], vanishes.
If the system has no singular points, shots are carried out by the inversion of the tridiagonal matrix that corresponds to Equation (29), with finite difference discretizations, such that the lower boundary conditions are satisfied. We find that magneto-gravity singularities are false singularities (i.e., L is singular but solutions are not; see Appendix B), so, for V A > 0, solutions in the magneto-gravity singularity ω-range can also be treated as regular everywhere. For solutions in the Alfvén singularity ω-range, we construct Frobenius power series solutions in the singular region (see Appendix B), fix constants of integration by shooting into, and matching with, the y = 0 boundary conditions, before finally shooting towards y = L y to obtain G[ω T ]. Solutions are then checked via back-substitution.
As discussed above, Sturm-Liouville theory only guarantees real eigenvalues in the weakly-magnetic limit. Therefore, we examine convergence for complex test frequencies, which have G = G r + iG i = 0 for G r , G i ∈ R. We find that G r G i is antisymmetric about ω i = 0, with contours G r = 0 and G i = 0 crossing exclusively on the real line, so ω ∈ R. We find the position of eigensolutions on the real line using the bracketed Newton-Raphson method discussed in Press et al. (1992).

Free wave eigensolutions
We label non-singular eigensolutions with a meridional mode number, n, based on the hydrodynamic convention. Generally, when the domain is finite and large enough, magnetic eigenfunctions for solutions without singularities are qualitatively similar to their hydrodynamic counterparts and n is the number of internal points wherev(y) = 0 in −L y < y < L y . However, hydrodynamic Kelvin solutions have the propertyv = 0 everywhere so represent a special case. They are typically labelled with the meridional mode number n = −1, with ψ −1 = 0 (Matsuno 1966). We find that solutions with c g k ≤ |ω| ≤ (c 2 g + V 2 A e) 1/2 k are the magnetic versions of Kelvin solutions, so we label them with n = −1 for consistency, although we find they have small non-zerô v (see below). For hydrodynamic and weakly-magnetic systems there are three solutions for each n ≥ 1: one equatorial Rossby/magneto-Rossby solution, one westward equatorial IG/MIG solution, and one eastward equatorial IG/MIG solution. When magnetism is included another two sets of solutions (one east; one west), with |ω| ≤ V A e 1/2 k, emerge. These solutions, which have Alfvén singularities (where ω 2 = B 0 (y s ) 2 k 2 ), differ .v (blue) and u (red) are respectively purely real and purely imaginary for the normalisation we apply. We mark asymptotes at y = ±ys with dotted black lines. The solutions are labelled with the latitudinal mode number, n, based on the latitudinal dependence ofv. The corresponding profiles for VA = 0.15cg/R are qualitatively identical.
significantly from regular equatorial wave solutions (see below). For convenience, we label these with a meridional mode number, n, determined by the scale of latitudinal variations inv (for n = 1, 3, 5,v is plotted in Figure 7). In Table 2 we present oscillation frequencies, ω, for the n = 1, n = 3, and n = −1 free wave eigensolutions, with each row representing a specific type of equatorial wave (see caption). We present the oscillation frequencies for V A = 0, V A = 0.15c g /R and V A = 0.2c g /R and, in cases where eigenfunctions are finite everywhere, we plot the corresponding free wave eigenfunctions for the equatorial n = 1 and n = −1 waves in Figure 8.
Eastward and westward equatorial IG/MIG solutions are the system's most rapidly oscillating waves (with |ω| > c g k). The azimuthal background magnetic field slightly increases the phase speed of the MIG modes (see Table 2). However, their energy redistribution patterns remain qualitatively similar to their hydrodynamic IG counterparts (see Figure 8, rows one and three).
Kelvin/magneto-Kelvin solutions are characterised by zonally-dominated winds. The are two hydrodynamic Kelvin solutions: an eastward equatorial Kelvin solution, with ω = c g k,v = 0, {û,ĥ} ∝ exp(−y 2 /2L eq ), and a westward boundary Kelvin solution, with ω = −c g k, Figure 8. The regular equatorial n = 1 (rows one to three) and n = −1 (row four) free wave eigenfunctions (geopotential contours with overlaid velocity vectors) are plotted for VA = 0, VA = 0.15cg, and VA = 0.2cg, taking k = 1/R. We label rows according to their wave types (see Table 2). Solutions are calculated for −5Leq < y < 5Leq, but are cut off for −Rπ/2 < y < Rπ/2 (Leq/R ≈ 0.67). v = 0, {û,ĥ} ∝ exp(y 2 /2L eq ). 14 These hydrodynamic solutions are special cases of magneto-Kelvin eigensolutions, which have c g k ≤ |ω| ≤ (c 2 g + V 2 A e) 1/2 k. While hydrodynamic Kelvin solutions havev = 0 everywhere, we find that magneto-Kelvin solutions acquire a non-zerô v in order to maintain latitudinally-independent oscillation frequencies. This can be understood by combining Equations (25), (27) and (28) to yield For hydrodynamic Kelvin solutions, the lefthand and righthand sides of Equation (33) are identically zero throughout the domain; whereas magneto-Kelvin solutions have c g k ≤ |ω| ≤ (c 2 g + V 2 A e) 1/2 k, {û,ĥ} similar to their hydrodynamic counterparts, and a non-zerov that ensures Equation (33) remains balanced. Like in the hydrodynamic limit, we find two magneto-Kelvin solutions: an eastward equatorial magneto-Kelvin solution and a westward boundary magneto-Kelvin solution. Magnetism causes both varieties to have a small non-14 The westward boundary Kelvin solution is removed when the condition is {û,v,ĥ} → 0 as |y| → 0 is imposed (Matsuno 1966). zero meridional velocity component (|v/û| 1) and an increased |ω|, but both are characteristically similar to their hydrodynamical counterparts. For the equatorial magneto-Kelvin solution, this is illustrated in Figure 8, which shows its energy redistribution pattern remains qualitatively similar as V A is increased.
In the hydrodynamic version of the system, equatorial Rossby solutions propagate westward and oscillate slowly (|ω| < c g k), with their azimuthal phase speeds, |ω|/k, successively decreasing for larger n solutions. In the hydrodynamic limit, the structures of equatorial Rossby solutions are characterised by midto-high latitude geostrophic vortices (see Figure 8, row two, lefthand column). For weakly-magnetic equatorial magneto-Rossby solutions, we find that the presence of the azimuthal background magnetic field has little effect on the form of the waves' eigenfunctions, which are magnetogeostrophic in nature. Weakly-magnetic solutions adjust to the contribution of magnetic tension with small increases to their azimuthal phase speeds. However, when their oscillation frequencies are exceeded by the maximal background azimuthal Alfvén frequency (i.e., when V A ≥ e −1/2 |ω|/k), equatorial magneto-Rossby solutions enter the ω-range of Alfvén singularities and are Table 2. Oscillation frequencies, ω, for the n = 1, n = 3, and n = −1 equatorial wave solutions with the planetary scale azimuthal wavenumber, k = 1/R, are tabulated for three choices of VA. In the Solution type column we use the following shorthands: R/MR denotes Rossby/magneto-Rossby solutions, WIG/WMIG denotes westward inertial gravity/magneto-inertial gravity solutions, EIG/EMIG denotes eastward inertial gravity/magneto-inertial gravity solutions, WA denotes (singular) westward Alfvén solutions, EA denotes (singular) eastward Alfvén solutions, K/MK denotes equatorial Kelvin/magneto-Kelvin solutions, and BK/BMK denotes boundary Kelvin/magneto-Kelvin solutions.
n Solution type ω/(cg/R) ω/(cg/R) ω/(cg/R) removed from the system. Higher n equatorial magneto-Rossby solutions are removed for the weakest V A values, before successively lower n solutions are removed for larger V A values (as Alfvénic properties become dynamically important at larger and larger scales). We attribute the removal of the planetary scale equatorial magneto-Rossby solutions to the breaking of potential vorticity conservation in regions of large Lorentz force The shallow-water hydrodynamic definition of potential vorticity is, q = h −1 (∂v/∂x−∂u/∂y +f ) (e.g., Vallis 2006). In the non-diffusive, unforced, drag-free version of the SWMHD model, the potential vorticity evolution satisfies where J = (∂B y /∂x − ∂B x /∂y) z. Equation (34) shows that the curl of the Lorentz force generated by the horizontal magnetic field component generally prevents potential vorticity conservation in the magnetic limit. 15 Since the material conservation of potential vorticity is essential to the propagation mechanism of Rossby waves (e.g., see Vallis 2006), in regions of large Lorentz force their generation is inhibited. In magnetic systems, two additional sets of solutions emerge. These solutions have |ω| ≤ V A e 1/2 k, so contain singularities, yet present some distinguishable properties of Alfvén waves. Specifically, they arise in both eastward and westward travelling varieties, and |ω| increases with V A and n. To assess their nature as y → y s (for ω 2 = B 0 (y s ) 2 k 2 ), we use the Frobenius solutions discussed in Appendix B. In singular regions,v = O(ln |(y − y s )/L eq |), so by Equation (33) u = O([(y − y s )/L eq ] −1 ), meaning that |û/v| → ∞ as y → y s . This highlights that Alfvén singularities cause a wave barrier to emerge at y = y s , over-which wavedriven meridional energy/momentum transport mechanisms cannot cross. Since they are not finite everywhere, equatorial wave structures with Alfvén singularities cannot determine global energy redistribution in the same way that planetary-scale equatorial waves do in hydrodynamic hot Jupiter models. Hence, in the limit where magnetism becomes significant, dissipative and non-linear effects become essential for understanding equatorial dynamics. A non-singular analogue of these solutions could be present in systems that include these extra physical processes but, since we are focused on the breakdown of geostrophic balance, we do not investigate solutions of this kind further. 16 Thus far, we have discussed magnetic free wave solutions about a flat rest state. However, Tsai et al. (2014) and Debras et al. (2020) find the redistributing properties of waves can be altered by the presence of a 15 Further, Dellar (2002) showed that potential vorticity has no materially invariant counterpart in SWMHD. 16 In the very strong field limit, London (2017) identified "outer band" solutions akin to these Alfvenic solutions that were trapped in polar regions in linear non-diffusive beta-plane systems, but concluded that they do not have a finite global (linear, nondiffusive) counterpart in London (2018). Spherical (linear, nondiffusive) SWMHD waves studies in other geometries have found additional slow magneto-Rossby (Márquez-Artavia et al. 2017) and magnetostrophic (Heng & Spitkovsky 2009) type waves at the poles of shallow-water systems, which may be useful in explaining the dynamics of the polar MHD flows. Márquez-Artavia et al. (2017) also found polar trapping of the "fast" magneto-Rossby solutions, which can plausibly be related to the removal of equatorial magneto-Rossby solutions (i.e., magneto-Rossby waves could become confined to regions of the atmosphere less influenced by magnetism).
background zonal flow, though the fundamental characteristics of these waves remain unchanged. Compared to the system we have so far explored, taking u 0 = U 0 constant (as in Tsai et al. 2014), simply manifests itself in the trivial phase translation ω → ω * − U 0 k, where ω and ω * are oscillation frequencies for a background at rest and a background with a zonal flow respectively. For this translation, Alfvén singularities emerge where B 0 (y s ) 2 k 2 = (ω * − U 0 k) 2 = ω 2 , which is the same condition as the rest case. We have also considered solutions about the latitudinally dependent background state, u 0 = u 0 (y), finding that Alfvénic singularities, with similar Frobenius solution dependencies, emerge at points where B 0 (y s ) 2 k 2 = (ω * − u 0 (y s )k) 2 .

Comparisons with non-linear simulations
Our findings concerning Alfvén-Rossby wave coupling in an equatorial beta-plane model, with an equatoriallyantisymmetric azimuthal background magnetic field, are consistent with our developed theory of hotspot reversals from the simulations of Section 3. In the hydrodynamic limit, planetary scale geostrophic circulations associated with equatorial Rossby waves are free to recirculate energy between the equatorial and mid-tohigh latitudes in a manner described by Showman & Polvani (2011). In the weakly-magnetic limit, planetary scale circulations remain largely unchanged, with equatorial magneto-Rossby waves only altering slightly to account for the magnetic contribution to their magnetogeostrophic circulations. However, at a critical threshold magnetic tension becomes large enough to inhibit the magneto-geostrophic circulations associated with equatorial magneto-Rossby waves. This is the free wave manifestation of the obstruction of geostrophic circulations, which we identified as the trigger for hotspot reversals in Section 3. Here the analogy between global circulations and the standing wave description of linear steady-state solutions described by Showman & Polvani (2011) breaks down and the force balance description used in Section 3 is preferable. In Section 3, we saw that the meridional Lorentz force responsible for obstructing geostrophic circulations always has a corresponding westward component that, ultimately, results in hotspot reversals. Together, the developed theory of Sections 3 and 4 can be used to place a zero-amplitude limit on the reversal threshold, V A,crit . In linear theory, magnetic tension inhibits the propagation of equatorial Rossby waves, with the oscillation frequency ω R,n = −βk k 2 + (2n + 1)β/c g , when ω A,max ≥ |ω R,n |, where ω A,max = B 0,max k = V A e 1/2 k is the maximal Alfvén frequency. Our find-ings suggest that, when the slowest (largest n) equatorial Rossby wave that is important for supporting the planetary scale mid-to-high latitude geostrophic balance becomes inhibited by magnetic tension, geostrophic circulations are obstructed and hotspots are driven westward by the resulting zonal Lorentz force.
In Figure 5, we have overplotted the theoretical thresholds associated with the obstruction of the n = 1, 3, 5 equatorial Rossby solutions, for comparison with the zero-amplitude (∆h eq /H → 0) limits of the simulated reversal thresholds, V A,crit . We generally find acceptable agreement between the simulations and these theoretical criteria, noting that in the most physically relatable case, where τ rad is short, reversals occur at the point where the n = 1 equatorial Rossby wave is overcome by magnetic tension. When τ rad and τ drag act over longer timescales, Figure 5 suggests that the obstruction of geostrophic circulations is associated with the loss of larger n equatorial Rossby solutions. This is somewhat consistent with the standing wave description of linear steady-state hydrodynamic solutions, as geostrophic circulations in solutions with longer τ rad and τ drag timescales are located at higher latitudes (e.g., see Figure 1), so require contributions to their energy recirculation patterns from larger n equatorial Rossby waves (e.g., Matsuno 1966;Showman & Polvani 2011;Tsai et al. 2014). While a wave analysis with non-linear effects and diffusion may be able to more precisely define these weakly-forced limits, we note that this description provides a vast improvement on scaling predictions of typical toroidal field strengths on hot Jupiters, which have order of magnitude (or larger) uncertainties (discussion in Section 5).

THE MAGNETIC REVERSAL MECHANISM
In Section 3, we identified the mechanism that drives magnetic hotspot reversals in SWMHD simulations of hot Jupiters. We provide a schematic and summarised explanation of the mechanism in Figure 9 and its caption.
This mechanism is also relevant for other, less idealised, magnetic field geometries. The reversal mechanism requires two features in the azimuthal field geometry: (1) large |B x | at mid-to-high latitudes to block/obstruct the circulation of the energy transporting geostrophic flows and (2) smaller or zero |B x | at equatorial latitudes, so that when magnetic field lines are bent into the equatorial region (by the mid-to-high latitude circulations), they pass into regions of smaller |B x |, generating a westward Lorentz force acceleration. This suggests that, as long as the profile is characterised by these two features, the developed theory does not de-  Figure 9. A schematic of the magnetic reversal mechanism, with grey temperature contours and white magnetic field lines (solid for Bx > 0; dashed for Bx < 0). (a) In hydrodynamic steady state solutions, drag-adjusted geostrophic circulations dominate at mid-to-high latitudes; whereas zonal pressure-driven jets dominate at the equator. Hotspots are shifted eastward as these circulations transport thermal energy from the western equatorial dayside to the eastern equatorial dayside, via higher latitudes. (b) In ultra-hot Jupiters, partially-ionised winds flow through the planet's deep-seated magnetic field, inducing a dominant equatorially-antisymmetric atmospheric toroidal magnetic field. When field lines are parallel to the equator, magnetic tension is zero, so flows behave hydrodynamically. (c) As the field and flow couple, the geostrophic circulations bend the magnetic field lines poleward on the western dayside and equatorward on the eastern dayside, generating a Lorentz force, (B · ∇)B. The meridional Lorentz force component acts to resist the geostrophic circulations; whereas, since |Bx| is smallest in equatorial regions, the zonal Lorentz force component, (B · ∇)Bx, is westward in hotspot regions, where field lines bend equatorward (and vice versa where field lines bend poleward). (d) Beyond a magnetic threshold, the system's nature changes. The meridional Lorentz force obstructs the circulating geostrophic winds, causing zonal wind alignment. This confines thermal structures and blocks the hydrodynamic transport mechanism. The zonal Lorentz force accelerates winds westward in the hottest dayside regions, causing a net westward dayside temperature flux. This drives the hottest thermal structures westward, until zonal pressure gradients can balance the zonal Lorentz force.
pend on exact antisymmetry in the dominant magnetic field geometry. This observation is useful when comparing to the 3D MHD simulations of Rogers & Komacek (2014) and Rogers (2017), which are characterised by antisymmetrically-dominant, but not exactly antisymmetric, toroidal magnetic field geometries.

Hotspot reversal criterion
In Sections 3 and 4, we identified two physicallymotivated reversal criteria on the Alfvén speed. The azimuthal Alfvén speed is defined as where µ 0 and ρ are the permeability of free space and the density. Taking c g = √ RT (see Section 2) and applying the ideal gas law therefore yields where T and P are the temperature and pressure at which the reversal occurs. From this, we have the following critical reversal criterion on the toroidal field mag-nitude: where n = 1 (largest scale Rossby wave) and κ ≈ 1 (B φ approaches maximal amplitudes close to the equator, as in Rogers & Komacek 2014) have been taken. This criterion quantifies the toroidal field magnitude sufficient to obstruct geostrophic circulations, with the first term in the maximum relating to when the toroidal field inhibits the propagation of the largest scale equatorial Rossby wave (in the small ∆h eq /H limit). Further, if the electric currents that generate the planet's assumed deep-seated dipolar field are located far below the atmosphere, Menou (2012) argued that the toroidal and dipolar field magnitudes should be related by the scaling law: is the magnetic Reynolds number and U φ is the magnitude of zonal wind speeds. We use the toroidal field criterion, and apply B φ ∼ R m B dip , to quantitively compare the predictions of SWMHD theory to the 3D MHD simulations of Rogers & Komacek (2014) and Rogers (2017).

Linking hotspot and wind reversals
Thus far, we have considered hotspot reversals, rather than the reversal of zonal-mean zonal winds,ū. Though time-correlated in 3D MHD models (Rogers 2017), hotspot and wind reversals are not necessarily synonymous. While thermal/wind structures and geopotential/wind structures compare well between hydrodynamic shallow-water and 3D models (e.g., Perez-Becker & Showman 2013;Komacek & Showman 2016), Debras et al. (2020) found a consistent treatment of the vertical component of the eddy-momentum flux (i.e., the vertical Reynolds stress) is critical to the development of equatorial superrotation (ū > 0). 17 In hydrodynamic models of hot Jupiters, equatorial superrotation emerges from the momentum transport mechanism of Showman & Polvani (2011). Showman & Polvani (2011) noted that the necessity for such a mechanism is a consequence of an angular momentum conservation theorem arising from Hide (1969b), which implies that equatorial superrotation can only be maintained if driven by an up-gradient angular momentum pumping mechanism. Showman & Polvani (2011) showed that this up-gradient mechanism is provided by the same geostrophic circulations that result in eastward hotspots. Therefore, since we have shown that magnetically-driven hotspot reversals are caused by the obstruction these recirculation patterns, Hide's theorem provides an anti-theorem, which implies that the magnetically-driven hotspot reversals are accompanied by a disruption of superrotation.
The realisation of this anti-theorem can be identified in 3D MHD simulations. These found that mid-tohigh latitude vortical structures zonally-align and, consequently, the transport of eastward eddy-momentum (horizontal Reynolds stress) from mid-latitudes into equatorial regions is reduced at atmospheric depths were reversals occur (compare Figures 2, 9, and 11 in Rogers & Komacek 2014). Rogers & Komacek (2014) found that, when the up-gradient horizontal Reynolds stress component diminishes, westward equatorial zonal-mean zonal accelerations are driven by the remaining downgradient momentum transport components (i.e., the vertical Reynolds stress and the Maxwell stresses). Thus, the above the application of Hide's theorem provides a meaningful connection between wind reversals and the magnetically-driven hotspot reversals mechanism we have presented.

Wave dynamics and turbulence
While we have not modelled turbulence in this work, actual planetary flows are expected to be highly turbulent. In hydrodynamic planetary systems, wave arguments have historically proved useful for developing understanding of geostrophic turbulence and how its conservational properties relate to eddies. Specifically, potential vorticity conservation is fundamental for both Rossby wave propegation and geostrophic turbulence, so Rossby wave properties can be used to understand the structures of planetary scale turbulence (e.g., Rhines 1975;Vallis 2006). Rogers & Komacek (2014) found that the relationship between zonal jets and magnetic fields in 3D MHD simulations shared intermittent features with MHD turbulence on a beta-plane that were identified by Tobias et al. (2007). Hydrodynamic geostrophic turbulence and MHD beta-plane turbulence have very different characteristics. Amongst them, the wave-wave/wave-zonal flow interactions associated with the inverse cascade of geostrophic turbulence are replaced with interactions that result in a forward MHD cascade, with MHD interactions occurring over scales on (and below) the planetary scale when the azimuthal Alfvén wave frequencies exceed the planetary scale Rossby wave frequency ). This turbulence condition is remarkably similar to the hotspot reversal criterion we identified in the weakly forced regime, which was motivated by wave dynamics and the findings of non-turbulent SWMHD simulations. We attribute this kinship to the breaking of potential vorticity conservation in MHD models in regions of large horizontal Lorentz force, which inhibits geostrophic characteristics such as Rossby wave propagation (as discussed in Section 4). We also highlight that forcing and drags generate potential vorticity sources/sinks, so potential vorticity conservation is modified when drag and forcing treatments are strong, which is why reversal thresholds deviate from this simple criterion in the strongly forced limit.

Magnetic field evolution and structure
After the initial hotspot transition, long term temporal differences between SWMHD and 3D MHD models arise because SWMHD can only model the planetary dipolar field or the atmospheric toroidal field selfconsistently (see Section 2.2), meaning that it cannot take into account toroidal field induction from reversed conducting zonal winds passing through the planetary dipolar field. If a strong toroidal field can be maintained indefinitely, the shallow-water theory predicts completely reversed winds, even in 3D models. However, at the onset of the wind reversals, the induction caused by the reversed winds flowing through the deep-seated magnetic field will result in a reduction of the atmospheric toroidal field's magnitude. Hence, while the quasi-steady magnetically-driven wind reversals of SWMHD are useful for modelling the reversal process, in reality one would expect to see oscillatory wind variations as toroidal fields successively strengthen and weaken in a wind-up-wind-down cycle of the toroidal magnetic field. Wind variations of this kind can be both observationally inferred from the oscillating peak brightness offsets of HAT-P-7b (Armstrong et al. 2016) and directly measured in 3D MHD simulations of the HAT-P-7b parameter space (Rogers 2017). This in itself has the interesting consequence that the reversal mechanism may provide a saturation process for the atmospheric toroidal magnetic field.
Due to the density dependence of the Alfvén speed, B φ,crit has a ∼ P 1/2 pressure dependence (see Equation (36)). This explains why Rogers & Komacek (2014) and Rogers (2017) found that wind reversals first onset in the upper atmosphere, but move deeper for stronger field strengths. Furthermore, if the reversal mechanism is a toroidal field saturation process (as discussed above), B φ should not greatly exceed B φ,crit . Hence, B φ should decrease above the deepest region where reversals occur (since B φ,crit decreases upwards), which is a feature of the toroidal field profiles found in Rogers & Komacek (2014), though other processes may also cause an upwards reduction in B φ . Comparing the geometry of the toroidal fields in the quasi-steady reversed SWMHD solutions with those in oscillating 3D MHD solutions is difficult. However, when the toroidal field is approaching criticality in strength, we do find similarities between our toroidal field geometries and those of Rogers & Komacek (2014). In both models the equatoriallyantisymmetric toroidal fields couple to mid-to-high latitude circulations in a manner that bends them towards the equator from west to east, which we showed is a geometry that results in westward Lorentz force accelerations (see Section 3).

Quantitive comparisons with 3D MHD
In Table 3, we compare predictions of the reversal criterion to magnetic field strengths of in three 3D MHD simulations: M7b1 and M7b2 of Rogers & Komacek (2014), and the HAT-P-7b model of Rogers (2017), all of which display wind reversals at some critical pressure depth, P crit . In these estimates, we take T eq =T , ∆T = T day −T , τ rad = τ wave , and set ∆h eq /H = ∆T /T eq . For comparisons to the simulations of Rogers & Komacek (2014), we compare B φ,crit to |B φ |, the horizontallyaveraged toroidal field component at the end of the run; whereas, for the HAT-P-7b simulation of Rogers (2017), we estimate the critical dipolar field strength at the atmospheric base, B dip,crit,base . This is calculated using the B φ ∼ R m B dip scaling law of Menou (2012). We take η = 2 × 10 6 m 2 s −1 and U φ ∼ 10 2 m s −1 from 3D simulations, to yield B dip,crit ≈ 4.3 G at P = 1 mbar, then noting that the atmospheric base is located at r = 0.15R in the simulations yields B dip,crit,base = 7. We note that the reversal criterion compares reasonably to the magnitude of the horizontally-averaged toroidal component field in the simulations of Rogers & Komacek (2014), with uncertainties in T day bracketing the true |B φ | value. This occurs both at P = P crit and above P crit , supporting the idea of reversals providing a toroidal field saturation process. The prediction of B dip,crit,base = 7 lies within the range 3 G < B dip,crit,base < 10 G identified by Rogers (2017). We note that, while B φ,crit has dependencies on ∆T /T eq and τ rad , η can vary significantly between the day and night sides of ultra-hot Jupiters (by orders of magnitude). Therefore, current understanding of the connection between toroidal and poloidal fields on hot Jupiter is constrained by large uncertainties (in Table 3. Estimates of reversal criteria compared to field strengths and reversal criteria from 3D MHD simulations (see Section 5.2.4 for definitions and an accompanying discussion).
References-a Rogers & Komacek (2014); b Rogers (2017) B φ ∼ R m B dip ), which far outweigh uncertainties in the toroidal field criterion that we have developed.

DISCUSSION
In this work we have explained the atmospheric mechanics of magnetically-driven hotspot reversals in hot Jupiters using numerical (Section 3) and semi-analytic (Section 4) analyses of a SWMHD model (Section 2), where we have applied parameters based on the ultrahot Jupiter HAT-P-7b. In Section 5 we used the theory developed throughout this study to identify a criticality criterion and discussed our findings in the context of 3D MHD simulations. This criticality criterion can be used to place physically-motivated constraints on the magnetic fields of ultra-hot Jupiters with observed westward hotspots. It also represents the point where hydrodynamic models with Lorentz force mimicking Rayleigh drag treatments should be replaced with self-consistent MHD modelling. In Section 5, we also identified the link between wind reversals and the hotspot reversal mechanism, highlighted relevant shared features between modifications to wave dynamics and atmospheric turbulence, discussed the role of reversals on the magnetic field's evolution, and made quantitive comparisons between the reversal criterion 3D MHD simulations.
Using the numerical SWMHD simulations, we demonstrated that hotspot reversals occur when equatoriallyantisymmetric azimuthal components of the magnetic field are strong enough to obstruct the geostrophic circulations that transport energy to the eastern dayside in hydrodynamic models. The magnetic field geometry that results from this obstruction always drives westward Lorentz force accelerations in hotspot regions, causing hotspots to transition from east-to-west. Using this finding we identified a reversal criterion for the toroidal field in the strong forcing regime using a simple argument based on the timescales of the two competing processes.
The recent observational drive in exoplanet meteorology provides a timely backdrop around which theories regarding the mechanism of wind/hotspot reversals can be tested and developed. Observational constraints on atmospheric properties continue to improve whilst a combination of archival data and dedicated observational missions from Kepler, Spitzer, Hubble, TESS, CHEOPS (and in the future JWST) are accelerating our understanding of the atmospheric theory of exoplanets. Since the prediction of magnetically-driven wind variations in hot Jupiters (Rogers & Komacek 2014), westward hotspots/brightspots have been inferred on the hot Jupiters HAT-P-7b (Armstrong et al. 2016), CoRoT-2b (Dang et al. 2018), Kepler-76b (Jackson et al. 2019), WASP-33b (von Essen et al. 2020), and WASP-12b (Bell et al. 2019). In Hindle et al. (2019) we inferred that the observed westward offsets of CoRoT-2b are unlikely to be driven by magnetism and Dang et al. (2018) proposed that such observations of CoRoT-2b could be explained by nonsynchronous rotation. However, as HAT-P-7b, Kepler-76b, WASP-33b, and WASP-12b are all ultra-hot Jupiters, the westward hotspots/brightspots observations on these planets are likely to be driven by magnetism. In future work, we will estimate the magnetic field strengths sufficient to explain their westward hotspots/brightspots observations. The toroidal field hotspot/wind reversal criterion we have developed is observationally motivated and appears to reproduce results of 3D MHD simulations. While this criterion does come with uncertainties due to the simplifications we have made, there are currently much larger uncertainties in the B φ ∼ R m B dip magnitude scaling of Menou (2012), which is used to connect the poloidaltoroidal magnitudes. This is because η is highly temperature dependent (e.g., Rogers & Komacek 2014), meaning that R m can vary by orders of magnitude between sides of the same hot Jupiter. Three-dimensional mod-els can further inform about connections between the poloidal-toroidal fields, which cannot be studied with SWMHD. In particular, we highlight the need to develop theoretical understanding of the effects that accompany strong day-night η dependencies in hot Jupiter atmospheres. Such comparisons offer important testcases for the extension of dynamo theory into the hot Jupiter regime. The predictions and constraints in this paper are clearly not the end of the story and, ultimately, bespoke 3D MHD simulations offer the best prospect for providing accurate constraints on the magnetic field strengths of ultra-hot Jupiters. That said, understanding the reversal mechanism is an important theoretical step and the reversal criteria we have presented en-ables modellers/observers to gain intuition into the most important atmospheric characteristics concerning reversals, particularly if one is comparing between multiple hot Jupiters in an ensemble approach.
Further details of the concepts, models, and applications discussed in this work are included in A. W. Hindle's forthcoming PhD thesis.
We acknowledge support from STFC for A. W. Hindle's studentship (ST/N504191/1). T. M. Rogers and P. J. Bushby acknowledge the Leverhulme grant RPG-2017-035. We thank Andrew Gilbert, Andrew Cumming, Natalia Gómez-Pérez, and Toby Wood for useful conversations leading to the development of this manuscript.

A. LINEARISED STEADY STATE SOLUTIONS
To solve the linearised, non-diffusive, steady state SWMHD system considered in Section 4.1 (i.e., Equations (15) to (18)), we reduce the system to a single inhomogenenous ordinary differential equation of the form where F 1 (y), F 2 (y), and F 3 (y) are latitudinally dependent coefficient functions, L is the system's second order differential operator, and Q(y) is the system's source term We have omitted the exact dependencies of F 1 (y), F 2 (y), F 3 (y), Q(y) for steady forced solutions (due to their cumbersome forms). These can be provided upon reasonable request. If S(y) and u 0 (y) are symmetric about the equator and B 0 (y) is antisymmetric about the equator, L and Q are respectively symmetric and antisymmetric about the equator. Solutions of Equation (A1) on −L y < y < L y are obtained by noting that, since L and Q are respectively symmetric and antisymmetric about the equator, inhomogeneous solutions are antisymmetric (i.e.,v antisymmetric andû,ĥ,Â symmetric). Consequently, we solve Equation (A1) in the upper-half domain, 0 < y < L y , withv(L y ) = 0 (impermeability) andv(0) = 0 (antisymmetry), before reflecting solutions. This reduced boundary value problem is solved by inverting the tridiagonal matrix that corresponds to Equation (A1) with finite difference discretizations. We fix the equatorial boundary conditionv(0) = 0 and vary dv/dy| y=0 in order to satisfyv(L y ) = 0, converging upon dv/dy| y=0 with the complex equivalent of the bracketed Newton-Raphson method discussed in Press et al. (1992).

B. SINGULAR TEST SOLUTIONS IN THE LINEAR EQUATORIAL WAVE SOLVING METHOD
For V A > 0, we examine the nature of the test solutions of Equation (29) in the upper-half domain, 0 < y < L y , about the singular points, y = y s . For |ω| ≤ B 0,max k, y s is located where B 0 (y s )k = |ω| (Alfvén singularities); whereas for c g k ≤ |ω| ≤ (c 2 g + B 2 0,max ) 1/2 k, y s is located where B 0 (y s )k = (ω 2 − c 2 g k 2 ) 1/2 (magneto-gravity singularities). The method of Frobenius giveŝ v = C 1v1 + C 2v2 ,v 1 = ∞ n=0 a nŷ n+µ1 ,v 2 = Dv 1 ln |ŷ| + ∞ n=0 b nŷ n+µ2 ; (B2) whereŷ = (y − y s )/L eq , C 1 and C 2 are the constants of integration,v 1 andv 2 are the first and second fundamental solutions, a n , b n and D are constant coefficients to be set or determined, and µ 1 ∈ Z and µ 2 ∈ Z are the roots of the indicial equation given by Equation (29). About magneto-gravity singular points, µ 1 = 2 and µ 2 = 0, sô where one is free to set a 0 = 1, b 0 = 1, b 2 = 0 (in fact, or b 2 can be set to any constant), and use Equation (29) to determine D, a n , and b n . Solutions of this kind are not singular at y = y s , so magneto-gravity singularities are in fact false singularities where the solution remains finite as y → y s . About Alfvén singular points, µ 1 = 0 and µ 2 = 0, sô where one is free to set a 0 = 1, b 1 = 1, b 0 = 0 (again, or b 0 can be set to any constant), and use Equation (29) to determine D, a n , and b n . Solutions of this kind are dominated by thev = O(ln |ŷ|) component as y → y s , so solutions with Alfvén singularities have infinite discontinuities for D = 0 (which we always find).