MAGNETICALLY CONTROLLED CIRCULATION ON HOT EXTRASOLAR PLANETS

Dynamic meteorology is a phenomenologically rich subject because of the lack of separation of physical scales. In other words, differences in the microphysical nature of a given system can have profound effects on its macroscopic state. As a result, it comes as no surprise that circulation patterns on hot Jupiters generally do not resemble those on solar system gas giants (Showman et al. 2008; Menou & Rauscher 2009).

From a hydrodynamical point of view, the circulational modes of typical hot Jupiter atmospheres differ in two principal ways compared to solar system gas giants. The first and most obvious difference is the energetics. Unlike the outer solar system, hot Jupiters reside in an environment where the incoming stellar irradiation completely dominates over the intrinsic planetary heat-flux. As a result, circulation patterns on hot Jupiters are driven primarily by the congenital dayside-to-nightside temperature differences (Showman & Polvani 2011). Furthermore, concurrent with the cooling of the planetary interior, the top-down heating of the atmosphere ensures the onset of stable stratification in the observable (P > 100 bars) atmospheric region (Guillot & Showman 2002; Burrows et al. 2007a).

The second difference lies in the extent to which the atmospheres are rotationally dominated. While solar system gas giants rotate rapidly (i.e., $\mathcal {T}_{\rm {Jup}} \simeq \mathcal {T}_{\rm {Sat}} \simeq 10$ hr), hot Jupiters are thought to rotate pseudo-synchronously with their orbital periods (i.e., $\mathcal {T}_{\rm {HJ}} \sim 3\hbox{--}5$ days) as a result of tidal despinning (see Hut 1981). This implies that although still dynamically significant, rotation alone does not exhibit commanding control over the atmospheric flow.

Since the discovery of the first transiting extrasolar gas giant HD209458b (Charbonneau et al. 2000; Henry et al. 2000), numerous authors have explored the atmospheric dynamics of hot Jupiters with a variety of numerical techniques. Because of the inherent differences in the frameworks of the simulations, today there exists a hierarchical collection of results that correspond to variable degrees of sophistication. On the simpler end of the spectrum are two-dimensional (2D) shallow-water simulations (Cho et al. 2008, 2003; Langton & Laughlin 2008, 2007), while the more intricate global circulation models (GCMs) include solvers of the three-dimensional (3D) "primitive" equations (Cooper & Showman 2005; Showman et al. 2008; Menou & Rauscher 2009; Heng et al. 2011a) as well as the 3D fully compressible Navier–Stokes equations (Dobbs-Dixon & Lin 2008; Dobbs-Dixon & Agol 2013). Simultaneously, various groups have gone to different lengths in their treatment of radiative transfer, with exploited models ranging from simple prescriptions such as Newtonian cooling (Showman et al. 2008) to double-gray (Rauscher & Menou 2012) and non-gray (Showman et al. 2009) schemes. An important step toward delineating the correspondence among results obtained with different solvers has been recently performed by Heng et al. (2011b).

Although there are quantitative differences in the results generated by different GCMs, there is general agreement on the qualitative features of the circulation. Specifically, there are three aspects of interest. First, super-rotating zonal jets exist in all simulations. Their number (and naturally, the widths) ranges between 1 and 4, depending on the model (see Showman et al. 2009), but the relative sparsity of the jets compared to Jupiter and Saturn is understood to be a result of diminished rotation rate (Showman & Guillot 2002). Moreover, in a recent study, Showman & Polvani (2011) showed that the formation of jets is ordained by the interaction of the atmospheric flow with standing Rossby waves which in turn result from the strong difference in the radiative forcing between the planetary dayside and the nightside.

Second, the characteristic wind speeds produced by different models are consistent within a factor of a few, and are generally in the ∼km s⁻¹ range. This is likely a direct result of the overall similarity in the force-balance setup within the models. Specifically, Showman et al. (2011) argue that near the equator, the horizontal pressure-gradient acceleration caused by the asymmetric irradiation is balanced by the advective acceleration. Meanwhile, Coriolis force takes the place of advective acceleration as the primary balancing term in the mid-latitudes. Both force-balances yield ∼km s⁻¹ as the characteristic wind speeds, in agreement with the numerical models.

Finally, in GCMs that resolve the vertical structure of the atmosphere (e.g., Showman et al. 2008; Menou & Rauscher 2009; Heng et al. 2011a), eastward jets consistently dominate the lower atmosphere while the upper atmosphere is characterized by more or less symmetric dayside-to-nightside circulation. In other words, winds originate at the sub-solar point and flow to the anti-solar point over the terminator in the upper atmosphere. The transition between the circulation patterns takes place at P ∼ 0.1–0.01 bars and is a consequence of the substantial reduction of the radiative time constant with diminishing pressure (Iro et al. 2005). In particular, because the radiative adjustment timescale is much shorter than the advective timescale in the upper atmosphere, the flow is unable to perturb the temperature structure away from radiative equilibrium significantly. Figure 1 depicts a schematic representation of the characteristic features of atmospheric circulation on hot Jupiters.

Zoom In Zoom Out Reset image size

There exists another important, distinctive feature of hot Jupiter atmospheres, namely their non-negligible electrical conductivity (see Figure 2). Electrical conductivity in hot Jupiter atmospheres does not originate from the ionization of H or He, but rather from the stripping of the valence electrons belonging to alkali metals such as K and Na (Batygin & Stevenson 2010; Perna et al. 2010). While these elements are thought to be present in trace abundances (e.g., [K]/[H] ∼ 10^−6.5, [Na]/[H] ∼ 10^−5.5; Lodders 1999), temperatures of ∼2000 K at upper atmospheric pressures lead to total and partial ionization of K and Na, respectively. In fact, at mbar levels, the conductivity can reach values as high as σ ∼ 1 S m⁻¹ (Batygin et al. 2011; Rauscher & Menou 2012; Heng 2012). Furthermore, it is generally expected that much like solar system gas giants, hot Jupiters posses interior dynamos, which produce surface fields comparable to, or in slight excess of, Jupiter's field⁴ (e.g., B ∼ 3–30 G; Stevenson 2003; Christensen et al. 2009). Consequently, there is a distinct possibility that atmospheric circulation on hot Jupiters may be in part magnetically controlled. That is to say, highly irradiated atmospheres may be sufficiently conductive for the Lorentz force to play an appreciable, if not dominant, role in the force balance.

Zoom In Zoom Out Reset image size

Realizing the potential importance of the coupling between the mean flow and the planetary magnetic field, Perna et al. (2010) modeled the Lorentz force as a Rayleigh drag (a velocity-dependent force that opposes the flow) and incorporated it into the GCM previously utilized by Menou & Rauscher (2010). This effort was later amended by Rauscher & Menou (2013), who also modeled the Lorentz force as a Rayleigh drag but self-consistently accounted for spatial variability in the electrical conductivity (by extension the drag timescale) in the weather layer. The results obtained with dragged GCMs exhibit significant differences in the obtained flow velocities relative to the standard GCMs. Namely, Perna et al. (2010) found a factor of ∼3 decrease in the peak wind speeds as the background dipole magnetic field increased from B_dip = 3 G to B_dip = 30 G, while Rauscher & Menou (2012) found a similar decline in the jet speeds as the field increased from B_dip = 0 G to B_dip = 10 G. The magnetic limitation of the peak wind speeds is of considerable importance as it may prevent the global circulation from approaching a super-sonic state (note that the characteristic sound speed is of the order of $c_s \sim \sqrt{k_{\rm {B}} T / \mu } \sim 3$ km s⁻¹, where k_B is Boltzmann's constant, T is the temperature, and μ is the mean molecular weight), thereby inhibiting the formation of shocks.

Although dragged 3D GCMs clearly highlight the quantitative importance of the magnetic effects in hot Jupiter atmospheres, they fail to accentuate significant qualitative differences in the obtained flows. Specifically, much like conventional GCMs, magnetically dragged GCMs still show deep-seated zonal jets, overlaid by complex flow patterns that intersect the poles of the planets. This lack of qualitative differences may arise from two distinct possibilities. The first is that beyond diminishing the peak wind speeds, the background magnetic field has little effect on the global circulation. In actuality, this may very well be true for pressure levels where the circulation is dominated by zonal jets, because of the geometrical simplicity of the flow-field interactions. Indeed, the coupling between the zonal flow and the pole-aligned background dipole field is azimuthally symmetric: differentially rotating jets convert the poloidal field into toroidal field (Liu et al. 2008; Batygin & Stevenson 2010). As will be discussed in detail below, beyond the reduction of velocities, this conversion poses few dynamical ramifications for the jets. Furthermore, owing to higher pressure and somewhat diminished temperatures compared with the upper atmosphere (and the associated decrease in conductivity), the zonal jets may reside in the kinematic regime, where the effects of the Lorentz force are modest (Batygin et al. 2011; Menou 2012a).

The second possibility is that although in reality the interactions with the background field impel the circulation to strongly deviate from its purely hydrodynamical counterpart, the procedure of modeling the Lorentz force as a Rayleigh drag does not capture the essential features of the dynamics. This is likely true in the upper atmosphere, where azimuthal symmetry is broken and the flow takes on a more topologically complex form. After all, in such a setting there is no requirement for the Lorentz force to simply oppose the flow everywhere, as is done by Rayleigh drag. Thus, there is a distinct possibility that previous modeling efforts have consistently misrepresented the circulation patterns of the upper atmospheres of hot Jupiters. Accordingly, the investigation of this possibility is the primary purpose of this work.

A statistically sound comparison between theoretical models and observations requires the incremental decrease in the goodness of fit to outweigh the cost of introducing new degrees of freedom into the model (see Rodgers 2000 for an in-depth review). Within the context of extrasolar planets, the limitations in observational capabilities and the quality of the data render the construction of highly sophisticated models unjustified (Line et al. 2012). Although a rigorous comparison with observational data is not the focus of this paper, our modeling efforts will lie in the same rudimentary spirit. In other words, here we shall focus on understanding the qualitative, rather than quantitative, nature of the circulation. Numerous simplifying assumptions will be made and the representation of the flow (including flow velocities, dayside-to-nightside temperature differences, etc.) should only be viewed as approximate. However, unlike all previous works on the subject, the model we shall utilize will remain self-consistently magnetohydrodynamic (MHD). In taking this approach, we hope to successfully capture the essential features of magnetic effects in highly irradiated planetary atmospheres.

The paper is organized as follows. In Section 2, we describe the equations inherent to our numerical GCM and reproduce the main features of hot Jupiter atmospheric flows in the purely hydrodynamic regime. In Section 3, we discuss the qualitative features of the atmospheric flows, treating the Lorentz force as a hydrodynamic drag. Specifically, we develop an analytical theory for magnetically dragged circulation patterns in the upper atmosphere and test the resulting scaling law against numerical simulations with enhanced viscosity and explore the effects of varying the radiative timescale. In Section 4, we introduce a pole-aligned background magnetic field and demonstrate the transition of the upper atmosphere's dayside-to-nightside circulation into a globally zonal state with the onset of the background field. We conclude and discuss our results in Section 5.

Momentum. The circulation model solves the Navier–Stokes equation for an electrically conductive Boussinesq fluid

$$\begin{equation} \frac{D \mathbf {v}}{D t} = - 2 \mathbf {\Omega } \times \mathbf {v} - \frac{\nabla \mathcal {P}}{\bar{ \rho } } + \frac{\delta \rho }{\bar{\rho }} \mathbf {g} + \frac{\mathbf {J} \times \mathbf {B} }{\bar{\rho }} + \nu \nabla ^2 \mathbf {v} \end{equation} \tag{ 1 }$$

in a rotating spherical shell of finite thickness. Here, D/Dt = ∂/∂t + v · ∇ is the material derivative, v is the velocity vector, $\mathbf { \Omega } = (2 \pi / \mathcal {T}) \mathbf {\hat{z}}$ is the rotation vector, $\mathcal {P}$ is the modified pressure, ρ is the density, J is the current density and ν is the kinematic viscosity. The bar denotes an average, whereas δ denotes the perturbation away from the background state. The density and temperature are related to the pressure through the ideal gas equation of state⁵

$$\begin{equation} P = \frac{\rho }{\mu } k_{\rm {B}} T, \end{equation} \tag{ 2 }$$

where P is the total pressure. The dynamic domain where the equation is solved is confined above a rigidly rotating spherical shell. We denote the inner and outer radii of the atmosphere as r₁ < r₂, respectively.

Continuity. The model is formally 3D and the vertical component of the motion enters into the continuity equation

$$\begin{equation} \nabla \cdot \mathbf {v} = 0. \end{equation} \tag{ 3 }$$

However, the nearly constant density, incompressible fluid approximation prevents us from self-consistently modeling a radially extensive atmosphere. Indeed, the atmospheric density does not change, except by thermal expansion/contraction⁶

$$\begin{equation} \frac{\delta \rho }{\bar{\rho } } = - \frac{\delta T}{\bar{T}}. \end{equation} \tag{ 4 }$$

Consequently, we limit the thickness of the atmosphere to a single scale-height in our simulations: $r_2 - r_1 = \mathcal {H} = k_{\rm {B}} \bar{T} / \mu g$, where g is the acceleration due to gravity. Additionally, we set $r_1 = \mathcal {R}$, the radius of the planet.

Generally, because of the above-stated reasons, our model should be viewed as more closely related to the 2D shallow water GCMs (Cho et al. 2003; Langton & Laughlin 2008) rather than the family of 3D models (Showman et al. 2009; Perna et al. 2010; Rauscher & Menou 2012). However, proper treatment of the induction equation (see below) in absence of predescribed symmetry requires the model to remain 3D.

Energy. The energy equation, that governs the temperature, inherent to the model reads

$$\begin{equation} \frac{D T}{D t} = \kappa \nabla ^2 T, \end{equation} \tag{ 5 }$$

where κ is the coefficient of thermal diffusivity (kept constant throughout the computational domain). Strictly speaking, this equation governs diffusive heat flux and (in direct interpretation) is unsuitable for modeling a medium where energy transport is accomplished mainly by radiation. This is because the above equation rests on the approximations of short photon mean free path and the neglect of the temperature-dependence of the opacity (in the context of the Boussinesq treatment employed here, the latter makes sense but the former breaks down at pressure levels corresponding to optical depth of the order of unity, allowing for only a crude approximation to reality). However, if it is possible to relate κ to radiative properties of the gas, the above energy equation can still be used to effectively mimic the appropriate heat transport.

In a radiatively dominated atmosphere, the correct energy equation reads (Peixoto & Oort 1992)

$$\begin{equation} \frac{D T}{D t} = \frac{1}{\bar{\rho } c_p} \nabla \cdot \mathcal {F}, \end{equation} \tag{ 6 }$$

where c_p is the specific heat capacity at constant pressure and $\mathcal {F}$ is the radiative heat flux. The expression for the radiative heat flux reads (Clayton 1968)

$$\begin{equation} \mathcal {F} = \frac{4 \sigma _{\rm {sb}} \bar{T}^3}{3 \bar{\rho } \psi } \nabla T, \end{equation} \tag{ 7 }$$

where σ_sb is the Stefan–Boltzmann constant and ψ is the opacity. At this point, the relationship between κ and the atmospheric temperature, density, opacity, and heat capacity is obvious. However, before proceeding further, let us recall that, to an order of magnitude, the Newtonian cooling timescale τ_N is given by the ratio of the atmosphere's excess heat content to its excess radiative flux:

$$\begin{equation} \tau _{\rm {N}} \simeq \frac{\bar{\rho } c_p \mathcal {H}}{4 \sigma _{\rm {sb}} \bar{T}^3}. \end{equation} \tag{ 8 }$$

Consequently, we may express

$$\begin{equation} \kappa = \frac{4 \sigma _{\rm {sb}} \bar{T}^3}{3 \bar{\rho }^2 c_P \psi } \simeq \frac{\mathcal {H}^2}{\tau _{\rm {N}}}, \end{equation} \tag{ 9 }$$

where we have implicitly assumed an infrared optical depth of the order of unity at the pressure-level of interest. The relationship between κ and τ_N is convenient, as it can be related to previous works. In particular, Showman et al. (2008) have calculated τ_N using a state of the art radiative transfer model and tabulated the results on a pressure-temperature grid. Here, we utilize their computations as a guide in estimating the thermal diffusivity.

Note that we could have arrived at the relationship (9) more intuitively by dimensional analysis. Specifically, noting that the radial extent of the atmosphere is much smaller than the lateral extent, the relevant length scale is the vertical scale-height, $\mathcal {H}$. Meanwhile, because the heat transport is primarily radiative, τ_N is clearly the relevant timescale. Bearing in mind the units of diffusivities (i.e., m² s⁻¹), Equation (9) naturally emerges as an estimate.

Magnetic field. The evolution of the magnetic field is governed by the induction equation (Moffatt 1978)

$$\begin{equation} \frac{\partial {\mathbf {B}}}{\partial t} = \eta \nabla ^2 \mathbf {B} + \nabla \times (\mathbf {v} \times \mathbf {B}), \end{equation} \tag{ 10 }$$

where η = 1/μ₀σ is the magnetic diffusivity (kept constant throughout the computational domain) and μ₀ is the permeability of free space. Meanwhile, the absence of magnetic monopoles implies a divergence-free magnetic field:

$$\begin{equation} \nabla \cdot \mathbf {B} = 0. \end{equation} \tag{ 11 }$$

Once the structure of B is known, the current density (within the MHD approximation) is given by

$$\begin{equation} \mathbf {J} = \frac{1}{\mu _0} \nabla \times \mathbf {B}. \end{equation} \tag{ 12 }$$

At this point, the full set of governing differential equations is presented. Paired with a matching set of boundary and initial conditions, the system can be integrated forward in time self-consistently.

The equations are solved using a mixed spectral–finite difference algorithm and, following Kuang & Bloxham (1999), the spherical harmonic decomposition is taken up to degree ℓ_max = 33 in the latitude m_max = 21 in the azimuthal angle. The computational domain is broken up into 64 radial shells. The model is integrated forward in time until equilibration in the thermal, kinetic, and magnetic energies are attained.

The physical parameters employed in the numerical experiments we report are loosely based on the planet HD209458b (Charbonneau et al. 2000). Aside from being the first extrasolar planet found to transit its host star, it has become a canonical example used in the studies of hot Jupiter atmospheres (Burrows et al. 2007b; Snellen et al. 2010). To this day (along with HD189733b; Knutson et al. 2009) it remains the best characterized extrasolar planet. The object has a mass $\mathcal {M} = 0.69 M_{\rm {Jup}}$, a radius $\mathcal {R} = 1.35\, R_{\rm {Jup}}$, and a rotational (assumed synchronous with orbit) period of $\mathcal {T}_{\rm {HD209458b}} = 3.5$ days.

Momentum. Because we are integrating a quasi-2D atmospheric shell, intended to be representative of a single pressure level, we apply impenetrable, stress-free boundary conditions to the velocity field. Thus, flows are free to develop without friction along the boundaries, although at initialization the planet is cast into solid-body rotation meaning v = 0 at t = 0.

Energy. The temperature is initialized at $\bar{T} = 2000$ K with $\delta T/\bar{T} = \delta \rho /\bar{\rho } = 0$. The atmosphere is kept stably stratified by the imposition of a stable background temperature gradient (set to the adiabatic lapse rate h = g/c_p) similar to the implementation in other Boussinesq dynamo models with stable layers (e.g., Stanley & Mohammadi 2008; Stanley 2010; Christensen & Wicht 2008). Additionally, an azimuthally variable heat-flux is applied at the outer boundary, r₂, to account for the variable stellar irradiation (in practice it does not matter whether the variable heat flux is applied at the outer or the inner boundaries, since it is the horizontal temperature gradient that controls the flow). The functional form of the heat-flux is chosen to be that of the $\mathcal {F} \propto Y_1^1$ spherical harmonic, while the amplitude is taken to be a free parameter (see the discussion in the next section).

Magnetic field. Nominally, an electrical conductivity of σ = 1 S m⁻¹, characteristic of P ≃ 1 mbar, T ≃ 2000 K is prescribed to the computational domain. Simultaneously, negligible electrical conductivity is assigned outside the computational domain. That is, σ ≃ 0 at r < r₁ and r > r₂, meaning that all of the current generated by the atmospheric flow is confined to the weather layer. This is in agreement with the approach of Perna et al. (2010) and Rauscher & Menou (2013), but in some contrast with the analytical work of Batygin & Stevenson (2010) and Batygin et al. (2011), where the current is allowed to penetrate the convective interior of the planet.

Because this work is primarily aimed at studying the upper atmosphere, the pressure-level of interest may reside above the atmospheric temperature inversion, characteristic of some hot Jupiter atmospheres (Burrows et al. 2007b; Spiegel et al. 2009). The presence of such an inversion may provide an electrically insulating layer,⁷ justifying σ ≃ 0 at r < r₁. The capability of electrical currents to penetrate the upper atmosphere and close in the magnetosphere is determined by the upper atmosphere's temperature structure and the abundance of alkali metals at P ≲ 1 mbar altitudes. Although a clear possibility, the physics of atmosphere–magnetosphere coupling is no doubt complex and is beyond the scope of the present study. Consequently, we neglect it for the sake of simplicity.

In all of our simulations, we initialize J = 0 in the weather layer. However, the weather layer is still permeated by the background magnetic field, B_dip, presumed to be generated by dynamo action in the convective interior of the planet. This zero current field reads

$$\begin{equation} \mathbf {B}_{\rm {dip}} = \nabla \times \left(k_m \frac{\sin (\theta)}{r^2} \hat{\mathbf {\phi }} \right), \end{equation} \tag{ 13 }$$

where k_m is a constant that sets the surface field strength, i.e., $|B_{\rm {dip}}| = k_m / \mathcal {R}^3$ at the equator. Within the context of the model, the background field is implemented by modifying the governing equations to account for a time-independent axially dipolar external field. This occurs specifically in two terms. In the momentum Equation (1), the Lorenz force $(\mathbf {J} \times \mathbf {B})/\bar{\rho }$ is replaced by $\mathbf {J} \times (\mathbf {B}+\mathbf {B}_{\rm {dip}})/\bar{\rho }$. In the magnetic induction Equation (10), the induction term ∇ × (v × B) is replaced by ∇ × (v × (B + B_dip)). This method is similar to that used by Sarson et al. (1997) and Cebron et al. (2012) to implement external fields.

Upon non-dimensionalization of the governing equations (Kuang & Bloxham 1999), we obtain six dimensionless numbers that completely describe the system. They are the Ekman number E, the magnetic Rossby number R_o, the magnetic Prandtl number q, the magnetic Rayleigh number R_th, and the aspect ratio χ, as well as an additional value, Γ, which parameterizes the incoming stellar flux. In terms of physical parameters, these quantities read

$$\begin{eqnarray} &E& \equiv \frac{\nu }{2 \Omega \mathcal {R}^2}, \ \ \ R_{\rm {o}} \equiv \frac{\eta }{2 \Omega \mathcal {R}^2}, \ \ \ \ \ \Gamma \equiv \frac{\mathcal {F}}{\kappa \bar{\rho } c_p h}, \nonumber \\ &q& \equiv \frac{\kappa }{\eta }, \ \ \ \ \ \ \ \ R_{\rm {th}} \equiv \frac{g h \mathcal {R}^2}{2 \bar{T} \Omega \eta }, \ \ \ \ \ \chi \equiv \frac{\mathcal {H}}{\mathcal {R}}. \end{eqnarray} \tag{ 14 }$$

The aspect ratio of the atmosphere we consider is χ = 7 × 10⁻³. For numerical stability, the model requires the Ekman number to exceed a critical value (which we adopt in the simulations) of the order of E_crit ∼ 10⁻⁵. If we consider the molecular viscosity of hydrogen $\nu = \bar{n} a \sqrt{3 k_{\rm {B}} \bar{T}/m_{\rm {H_2}}}/2$, where $\bar{n}$ is the number density and a is the molecular cross-section, we obtain a hopelessly small E ∼ 10⁻²⁰ ≪ E_crit at mbar pressure. We note that the actual Ekman number is orders of magnitude higher, since on the global planetary scale, small-scale turbulence is a far more relevant source of viscosity than transfer of momentum at the molecular level (Peixoto & Oort 1992). Still, the true value of E is probably much smaller than E_crit.

For a given electrical conductivity (equivalently diffusivity), R_o and R_th simply encompass the physical units of the model. Specifically, for σ = 1 S m⁻¹, R_o = 1.79 × 10⁻⁶ and R_th = 1.38 × 10⁹. Meanwhile, all of the information regarding the considered pressure level is provided by q. Taking σ = 1 S m⁻¹ as before, we obtain q ≃ 6, 60 and 600, corresponding to τ_N = 10⁵, 10⁴, and 10³ s, appropriate for P = 1, 0.1 and 0.001 bars, respectively. Finally, as already mentioned above, we take Γ to be an adjustable parameter, tuned to obtain the desired temperature gradients. It is important to note that our model differs somewhat from typical dynamo models in a sense that Γ, rather than R_th, parameterizes the extent to which the system is driven.

Let us conclude the description of the numerical model with a brief discussion of its shortcomings as a guide for future work. First and foremost, the limitation of the atmosphere to a single scale-height may be prejudicial, as previous work has shown that including an extended vertical extent of the atmosphere is important to capture circulation features (Heng et al. 2011b). Second, while we have kept thermal and magnetic diffusivities uniform throughout the computational domain, it should be understood that in reality these values vary with temperature and pressure. Third, an implicit assumption of an infrared optical depth of the order of unity is rather crude at ∼mbar pressure levels and should be lifted in more rigorous treatments of radiative transfer. Finally, as already mentioned above, the artificially enhanced viscosity inherent to our model almost certainly smoothes out smaller scale flow to an unphysical extent.

Prior to performing unabridged MHD simulations, it is useful to first compare our model to previously published results. Accordingly, we begin by setting the strength of the background field to B_dip = 0. Because the computational domain is initialized with a null current density, no magnetic fields are generated yielding purely hydrodynamic simulations.

Horizontal slices of the atmosphere (through the center of the computational domain) in the cylindrical projection are shown in Figure 3. The figures are centered on the substellar point and the background color shows the temperature distribution. The arrows denote the circulation vector field. Peak wind speeds as well as the maximal temperature deviations from $\bar{T} = 2000$ K are labeled.

Zoom In Zoom Out Reset image size

The typically observed transition from zonal flows to dayside–nightside flows is observed as the pressure is decreased sequentially from P = 1 bars (panel (C) of Figure 3) to P = 0.001 bars (panel (A) of Figure 3). However, important differences exist in our results, contrasted against, say, the results of Showman et al. (2008). The first important distinction is that in the zonally dominated parameter regime, rather than developing a single broad jet, our model shows the development of three counter-rotating jets.

This is a simple consequence of angular momentum conservation in the computational domain, and is not uncharacteristic of 2D models (Showman et al. 2011). Because of free-slip boundary conditions employed in our model, the atmosphere is not allowed to exchange angular momentum with the interior. As such, because the atmosphere is initialized in solid-body rotation, the development of any prograde jets must be accompanied by the development of retrograde jets. Because the retrograde jets reside at a high latitude θ_ret, whereas the prograde jets are essentially equatorial, angular momentum conservation requires them to be faster by a factor of |v_ret/v_pro| ∼ 1/cos θ_ret, as observed in the simulations. The angular momentum conserving 3D simulations of Heng et al. (2011b) exhibit similar behavior, although in their model the counter-rotating jets develop below the prograde ones, and are considerably slower because of the associated density enhancement.

The second distinction of interest is the direction of the equatorial jet. While 3D GCMs consistently produce eastward equatorial jets, the equatorial jets in our hydrodynamical simulations are westward. This is not too surprising, as shallow water and equivalent barotropic models are known to produce both eastward and westward equatorial jets, depending on the details of the simulation setup (Showman et al. 2011; Heng et al. 2011b). Consequently, this difference can probably be attributed to the limited vertical extent of our model.

The third important difference is the fact that flow velocities are not consistent along the pressure levels. In particular, for the same values of δT, zonal flows have ∼few km s⁻¹ peak wind speeds, while the dayside-to-nightside flows are more than an order of magnitude slower. This is largely a consequence of the difference in the geometry of the circulation and the fact that viscosity enters as a significant member in the force balance for dayside-to-nightside flows. This can be seen by approximating $\nu \nabla ^2 \sim \nu / \mathcal {L}^2$ (Holton 1992), where $\mathcal {L}$ is a characteristic length scale associated with the curvature of the circulation. For zonal jets, $\mathcal {L}_{\rm {z}} \sim \mathcal {R}$, while for dayside-to-nightside circulation, the return flow is in part radial implying $\mathcal {L}_{\rm {dn}} \sim \mathcal {H}$. As will be shown in the following section, faster velocities can be attained by either increasing the aspect ratio of the atmosphere or by artificially enhancing the radiative heat flux as a result of a linear proportionality between peak wind speeds and δT. An additional point of importance is that in typical 3D simulations, dayside-to-nightside circulation is nearly uniform over the terminator with the return flow residing at greater depth, while our results depict a partial return flow over the poles. This is again a consequence of the quasi-2D geometry of our model.

Although these quantitative distinctions are certainly worthy of attention, the typical features of the flow are approximately captured by our simplified model. Consequently, while being mindful of the model's limitation, we do not view the quantitative dissimilarities as critical, as they are not central to the argument of the paper. After all, recall that we are primarily concerned with the possibility of a qualitative alteration of the dayside-to-nightside flow by magnetic effects.

Let us begin with the estimation of characteristic timescales. In order to accomplish this, we first simplify the Lorentz and viscous forces to resemble the functional form of Rayleigh drag. Utilizing Ohm's law, we have⁸

$$\begin{eqnarray} \frac{\mathbf {J} \times \mathbf {B}}{\bar{\rho }} &\sim & \frac{\sigma (\mathbf {v} \times \mathbf {B}_{\rm {dip}}) \times \mathbf {B}_{\rm {dip}} }{\bar{\rho }} \sim - \left(\frac{\sigma k_m^2}{\mathcal {R}^6\bar{\rho }} \right) \mathbf {v} = - \frac{\mathbf {v}}{\tau _{\rm {L}}} \nonumber \\ \nu \nabla ^2 \mathbf {v} &\sim & - \left(3 \frac{\nu }{\mathcal {R}^2} \right) \mathbf {v} = - \frac{\mathbf {v}}{\tau _{\nu }}. \end{eqnarray} \tag{ 15 }$$

This allows us to rewrite Equation (1) in a simpler form:

$$\begin{eqnarray} &&\frac{D \mathbf {v}}{D t} = - 2 \mathbf {\Omega } \times \mathbf {v} - \frac{\nabla \mathcal {P}}{\bar{ \rho } } + \frac{\delta \rho }{\bar{\rho }} \mathbf {g} - \frac{\mathbf {v}}{\tau _{f}}, \end{eqnarray} \tag{ 16 }$$

where τ_f = (1/τ_L + 1/τ_ν)⁻¹. Taking |B_dip| = 1 G, the characteristic timescales are: τ_L ∼ 10³ s and τ_ν ∼ 10⁵ s. Other relevant timescales in the problem are the rotational (Coriolis) timescale τ_Ω ∼ 2π/Ω ∼ 10⁵ s, radiative timescale τ_N ∼ 10³ s, and the advective timescale $\tau _{\rm {adv}} \sim \mathcal {R}/v \sim 10^{5}$ s.

There exists a clear separation of timescales in the system. As a result, upon including the parameterized Lorentz force into the equations of motion, the inertial and Coriolis terms can be dropped (the viscous term can be dropped as well, although this simplification is unnecessary). The removal of the inertial terms implies a steady-state solution. The removal of the Coriolis term creates a symmetry characterized by an axis that intersects the sub-solar and anti-solar points. Taking advantage of this symmetry, we orient the polar axis of the coordinate system such that it intersects the sub-solar point. Upon doing so, we can specify a null azimuthal velocity and drop all azimuthal derivatives in the equations of motion. In a local Cartesian reference frame, this leaves us with horizontal ($\hat{\mathbf {y}}$) and vertical ($\hat{\mathbf {z}}$) momentum equations, where the latter simplifies to the equation of hydrostatic balance:

$$\begin{eqnarray} &&\frac{1}{\bar{\rho }}\frac{\partial \mathcal {P}}{\partial y} = - \frac{v_y}{\tau _f} \nonumber \\ &&\frac{1}{\bar{\rho }}\frac{\partial \mathcal {P}}{\partial z} = g \frac{\delta T}{\bar{T}}. \end{eqnarray} \tag{ 17 }$$

Following Schneider & Lindzen (1977) and Held & Hou (1980), we shall adopt a Newtonian energy equation with implicit stable stratification (recall that we have set the potential temperature gradient to h = g/c_p):

$$\begin{equation} \xi v_z \left(\frac{g}{c_p} \right) = - \frac{\delta T - \delta T_{\rm {rad}}}{\tau _{\rm {N}}}. \end{equation} \tag{ 18 }$$

Here, ξ ⩾ 1 is a constant that parameterizes lateral heat advection and δT_rad is a purely radiative perturbation to the background state, $\bar{T}$. In other words, as the damping of the circulation is strengthened and v → 0, δT → δT_rad.

Retaining the incompressibility condition (3), we introduce a stream-function Ψ, defined through v = ∇ × Ψ (Landau & Lifshitz 1959). Taking a partial derivative of the y-momentum equation with respect to z and of the z-momentum equation with respect to y, we obtain

$$\begin{equation} -\frac{1}{\tau _{f}} \frac{\partial ^2 \Psi }{\partial z^2} = \frac{g}{\bar{T}} \frac{\partial (\delta T)}{\partial y}. \end{equation} \tag{ 19 }$$

Taking a derivative of Equation (19) with respect to y and switching the order of partial differentiation yields

$$\begin{equation} - \frac{1}{\tau _{f}} \frac{\partial }{\partial z} \left(\frac{\partial ^2 \Psi }{\partial y \partial z} \right)= \frac{g}{\bar{T}} \frac{\partial ^2 (\delta T)}{\partial y^2}. \end{equation} \tag{ 20 }$$

Meanwhile, differentiating the Newtonian cooling Equation (18) with respect to z gives

$$\begin{equation} \left(\frac{\partial ^2 \Psi }{\partial y \partial z} \right) = \frac{c_p}{g \xi } \frac{\partial }{\partial z} \left(\frac{\delta T - \delta T_{\rm {rad}}}{\tau _{\rm {N}}} \right), \end{equation} \tag{ 21 }$$

allowing us to eliminate ∂²Ψ/∂y∂z:

$$\begin{equation} \frac{\partial ^2 (\delta T_{\rm {rad}})}{\partial z^2} - \frac{\partial ^2 (\delta T)}{\partial z^2} = \xi \frac{ \tau _{f} \tau _{\rm {N}} g^2}{c_p \bar{T}} \frac{\partial ^2 (\delta T)}{\partial y^2}. \end{equation} \tag{ 22 }$$

Note that $g^2/c_p \bar{T}$ is simply the square of the Brunt–Vaisala frequency for an isothermal atmosphere.

Equation (22) admits the trial solutions

$$\begin{eqnarray} &\delta & T_{\rm {rad}} = (\delta T_{\rm {rad}}^0) \cos \left(\frac{y}{\mathcal {R}} \right) \sin \left(\frac{\pi z }{\mathcal {\mathcal {H}}} \right) \nonumber \\ &\delta & T = (\delta T^0) \cos \left(\frac{y}{\mathcal {R}} \right) \sin \left(\frac{\pi z }{\mathcal {\mathcal {H}}} \right), \end{eqnarray} \tag{ 23 }$$

where $(\delta T_{\rm {rad}}^0)$ and (δT⁰) are constants that represent the maximal deviations in the respective quantities from the background state. Note that $(\delta T_{\rm {rad}}^0)$ is a parameter inherent to the model rather than a variable. Upon substitution of the above solutions (23) into Equation (22), we obtain a relationship between (δT⁰) and $(\delta T_{\rm {rad}}^0)$:

$$\begin{equation} (\delta T^0) = (\delta T_{\rm {rad}}^0) \left[ 1 + \xi \zeta \frac{ \tau _{f} \tau _{\rm {N}} g^2}{\pi ^2 c_p \bar{T}} \left(\frac{\mathcal {H}}{\mathcal {R}} \right)^2 \right]^{-1}. \end{equation} \tag{ 24 }$$

In the above equation, ζ is an empirical factor that has been introduced to account for the approximations inherent to Equations (15).

With an analytical solution for the temperature perturbation in hand, we substitute Equations (23) into Equation (19) and integrate twice to obtain

$$\begin{equation} \Psi = \left[ \frac{g \zeta \tau _{f} \mathcal {H}^2 (\delta T^0) }{\pi ^2 \bar{T} \mathcal {R}} \right] \sin \left(\frac{y}{\mathcal {R}} \right) \sin \left(\frac{\pi z}{\mathcal {H}} \right) \hat{\mathbf {x}}. \end{equation} \tag{ 25 }$$

This solution implies the same functional form for laterally averaged heat transport in the vertical and horizontal advection terms in the energy Equation (6), lending some support for the approximation inherent to Equation (18). Once the stream function is obtained, the maximal horizontal and vertical velocities are given by

$$\begin{equation} v_y^{\rm {max}} = \Psi _0 \frac{\pi }{\mathcal {H}} \ \ \ \ \ v_z^{\rm {max}} = \Psi _0 \frac{1}{\mathcal {R}}, \end{equation} \tag{ 26 }$$

where Ψ₀ is the term in square brackets in (25). Note that the above theory automatically implies a quasi-2D flow since $v_z^{\rm {max}}/v_y^{\rm {max}} \sim (\mathcal {H}/\mathcal {R}) \ll 1.$

With a simple analytical theory at hand, we performed a series of numerical simulations, varying the radiative and drag timescales in the ranges 10² < τ_N < 10³ s and 10³ < τ_f < 10⁵ s. Although we observed a considerable variability in the wind speeds and dayside-to-nightside temperature differences in our simulations, the nature of the flow was largely the same as that seen in panel (A) of Figure 3 across the runs. The peak wind speeds obtained in the simulations as functions of τ_N and τ_f are presented as black dots in Figure 4.

Zoom In Zoom Out Reset image size

In addition to the simulation results, Figure 4 shows $v_y^{\rm {max}}$ given by Equation (26) for the same parameters. As can be seen from the figures, the scaling law inherent to Equation (25) matches the numerical experiments quite well. Extrapolating toward larger values of τ_f, it can be inferred that our simulations would have produced peak wind speeds of the order of ∼km s⁻¹ if not for the numerical requirements of enhanced viscosity. However, it can also be expected that the character of the flow would also change qualitatively with diminishing viscosity, as the force balance shifts away from that implied by Equation (15).⁹

It should be kept in mind that the adjustable parameters ξ and ζ were fit to the data. Moreover, the value of $(\delta T_{\rm {rad}}^0) = 3360$ K¹⁰ was chosen by running a simulation where viscous forces completely dominated the force balance ensuring v = 0. In other words, the quantitative agreement seen in Figure 4 is a consequence of the fact that the adjustable parameters of the analytical solution have been fit to the numerical data, but the fact that the functional form of the analytical model conforms with numerical experiments suggest that the stream-function (25) captures the main features of dragged upper-atmospheric circulation on hot Jupiters.

With the exception of a rather limited number of problems, self-consistent MHD solutions can only be attained with the aid of numerical simulations. However, for the system at hand, the qualitative effect of magnetic induction can be understood from simple theoretical considerations. As we already argued, the two characteristic states of hot Jupiter atmospheric circulation are a zonally dominated state and a meridionally dominated state (Showman et al. 2011). Whether or not a given configuration will be significantly affected by the introduction of the magnetic field can be established by analyzing its stability. More specifically, we can work within a purely kinematic (rather than dynamic) framework to understand if the Lorentz force acts to perturb the flow away from its hydrodynamic counterpart or simply damps the circulation.

Zonal flows. Although not directly applicable, recent studies of Ohmic dissipation that arises from zonal flows performed by Liu et al. (2008; within the context of solar system gas giants) and by Batygin & Stevenson (2010) as well as Menou (2012b; within the context of hot Jupiters) have already produced some results on a related problem. Here we work in the same spirit as these studies and prescribe the following functional form to the zonal flow to approximately represent three jets, such as those shown in panel (C) of Figure 3:

$$\begin{equation} \tilde{\mathbf {v}} = \tilde{v}_0 \sin (3 \theta) \hat{\mathbf {\phi }}, \end{equation} \tag{ 27 }$$

where $\tilde{v}_0$ is a negative constant whose magnitude corresponds to the peak wind speed. This prescription trivially satisfies the continuity Equation (3), although we note that a more realistic zonal flow should also exhibit differential rotation.

The interaction between this flow and the background magnetic field (13) will induce a field B_ind in the atmosphere. Because B_dip is entirely poloidal, and $\tilde{\mathbf {v}}$ is strictly toroidal, B_ind will also be strictly toroidal (Moffatt 1978). As can be readily deduced from Equation (10), this means that B_ind cannot interact with $\tilde{\mathbf {v}}$ to further induce new field unless $\tilde{\mathbf {v}}$ deviates from a purely zonal flow. As a result, in steady state, the induction equation reads

$$\begin{equation} -\eta \nabla ^2 \mathbf {B}_{\rm {ind}} = \nabla \times (\tilde{v} \times \mathbf {B}_{\rm {dip}}) = -\frac{6 k_m \tilde{v}_0 }{r^4} \cos (\theta) \sin (\theta) \hat{\mathbf {\phi }}. \end{equation} \tag{ 28 }$$

It is noteworthy that had we chosen to represent a single broad jet (such as that seen in most 3D simulations; Showman et al. 2008; Menou & Rauscher 2010; Rauscher & Menou 2013) by setting $\tilde{v} \propto \sin (\theta)$ (in this case $\tilde{v}_0$ is positive), Equation (28) would have looked the same, with the exception of the coefficient on the right-hand side, which would have been 2 instead of 6. As a result, it should be kept in mind that the following kinematic solution applies to the case of a single jet as well.

The angular part of Equation (28) is satisfied by the expression

$$\begin{equation} \mathbf {B}_{\rm {ind}} = \mathcal {A}(r) \cos (\theta) \sin (\theta) \hat{\mathbf {\phi }}, \end{equation} \tag{ 29 }$$

where $\mathcal {A}(r)$ is a yet undefined function. This form ensures that the meridional component of the induced current vanishes at the poles. Meanwhile, the radial impenetrability of the boundaries requires $\mathcal {A}(\mathcal {R}) = \mathcal {A}(\mathcal {R}+ \mathcal {H}) = 0$ as dictated by Equation (12). With these boundary conditions, Equation (28) can be solved to yield

$$\begin{eqnarray} \mathcal {A}(r) &=& 3 k_m \tilde{v}_0 (\mathcal {R}-r) (\mathcal {H}^4+5 \mathcal {H}^3 \mathcal {R}+10 \mathcal {H}^2 \mathcal {R}^2 \nonumber \\ &&+10 \mathcal {H} \mathcal {R}^3+4 \mathcal {R}^4-\mathcal {R}^3 r-\mathcal {R}^2 r^2-\mathcal {R} r^3-r^4) /(2 \eta r^3 \nonumber \\ &&\times (\mathcal {H}^4+5 \mathcal {H}^3 \mathcal {R}+10 \mathcal {H}^2 \mathcal {R}^2+10 \mathcal {H} \mathcal {R}^3+5 \mathcal {R}^4)). \end{eqnarray} \tag{ 30 }$$

The induced field and the associated electrical current are shown in panel (A) of Figure 5.

Zoom In Zoom Out Reset image size

The Lorentz force that arises from the interactions between B_ind and B_dip takes the form

$$\begin{eqnarray} \mathbf {F}_{\rm {L}} &=& \frac{(\nabla \times \mathbf {B}_{\rm {ind}}) \times \mathbf {B}_{\rm {dip}}}{\rho \mu _0} = \left(\frac{ \sigma k_m^2}{r^6 \rho } \right) \tilde{v}_0 \hat{\mathbf {\phi }} \nonumber \\ &&\times (3 \sin (\theta) (-11 \mathcal {R} (\mathcal {H}+\mathcal {R}) (\mathcal {H}+2 \mathcal {R}) (\mathcal {H}^2+2 \mathcal {H} \mathcal {R} \nonumber \\ &&+ 2 \mathcal {R}^2)+9 r (\mathcal {H}^4+5 \mathcal {H}^3 \mathcal {R}+10 \mathcal {H}^2 \mathcal {R}^2+10 \mathcal {H} \mathcal {R}^3 \nonumber \\ &&+ 5 \mathcal {R}^4)-r^5)+\sin (3 \theta) (-7 \mathcal {R} (\mathcal {H}+\mathcal {R}) (\mathcal {H}+2 \mathcal {R}) \nonumber \\ &&\times (\mathcal {H}^2 + 2 \mathcal {H} \mathcal {R}+2 \mathcal {R}^2)+5 r (\mathcal {H}^4+5 \mathcal {H}^3 \mathcal {R} \nonumber \\ &&+10 \mathcal {H}^2 \mathcal {R}^2+10 \mathcal {H} \mathcal {R}^3+5 \mathcal {R}^4)+3 r^5))/(8 r \nonumber \\ &&\times (\mathcal {H}^4+5 \mathcal {H}^3 \mathcal {R}+10 \mathcal {H}^2 \mathcal {R}^2+10 \mathcal {H} \mathcal {R}^3+5 \mathcal {R}^4)). \end{eqnarray} \tag{ 31 }$$

Because the Lorentz force acts in the same sense as the flow itself (that is, $\mathbf {F}_{\rm {L}} \times \tilde{\mathbf {v}} = 0$), it can only act to accelerate/decelerate the jets but not alter their directions. Indeed, the functional form of F_L is that of a Rayleigh drag (Equation (15)), however, the characteristic timescale is non-uniform in latitude and radius, i.e., τ_L = f(r, θ). The non-uniformity we derive here should not be confused with the variability in F_L that can arise from the spatial dependence of the electrical conductivity (see Rauscher & Menou 2013).

It is noteworthy that the radial dependence of F_L can give rise to differential rotation. However, this is not particularly important, since in some similarity with the above discussion, differential rotation will only induce toroidal fields through the ω-effect (Moffatt 1978) and will therefore only change the solution obtained here on a detailed level (i.e., the added dependence of $\tilde{\mathbf {B}}$ on r will subtly modify the function $\mathcal {A}(r)$). In other words, a differentially rotating zonal flow still results in a purely toroidal induced field.

For a sensible comparison with previous works (e.g., Perna et al. 2010; Menou 2012a; Rauscher & Menou 2012) and the simple theory presented in the previous section, it is instructive to evaluate the maximal magnitude of F_L, which corresponds to the upper and lower boundaries of the domain in question,¹¹ i.e., $r = \mathcal {R}, r = \mathcal {R} + \mathcal {H}$. To leading order in χ, the expression reads

$$\begin{equation} \rm {max}(\mathbf {F}_{\mathrm{L}}) \simeq - 6 \left(\frac{\sigma \mathit {k}_m^2}{\mathcal {R}^6 \rho } \right) \tilde{v}_0 \chi \cos (\theta)^2 \sin (\theta) \hat{\mathbf {\phi }}. \end{equation} \tag{ 32 }$$

From this expression, it is clear that F_L acts primarily in the mid-latitudes rather than the equator. As a result, the damping of the jets is latitudinally differential, meaning that even if the flow is initially composed of multiple bands (as we consider here), it will approach a single equatorial jet as the conductivity and/or the magnetic field is increased. Furthermore, recall that the functional form of Equation (32) is also valid in the case of a single jet. Qualitatively, this seems to imply that the Lorentz force acts to collimate the jet toward the equator. Such an effect is sensible given that the radial component of the field is stronger as one approaches the pole for a simple dipole. However, it should also be kept in mind that a true planetary magnetic field might be more complicated, leading to further lack of triviality in the circulation.

Dayside-to-nightside flows. Let us now consider the more topologically complex interaction between meridional flows and a spin-pole aligned dipole magnetic field. As in Section 3.1, we shall work in a coordinate frame where the polar axis intersects the sub-solar point and is directed at the host star. Unlike the case of zonal circulation, this configuration has no exploitable symmetry. Consequently, a simple solution to the steady-state induction Equation (28) is difficult, if not impossible, to obtain. We shall therefore make substantial simplifications.

In our prescription for the velocity field, we neglect radial flow altogether (thereby violating continuity) and adopt an expression similar to Equation (27):

$$\begin{equation} \tilde{\mathbf {v}} = \tilde{v}_0 \sin (\theta) \hat{\mathbf {\theta }}. \end{equation} \tag{ 33 }$$

Because of our choice of coordinate system, Equation (13) cannot be used directly. However, keeping in mind that the background dipole field originates in a much deeper region of the planet than the atmosphere, we can write down the magnetic field in a current-free representation (Jackson 1998):

$$\begin{equation} \mathbf {B}_{\rm {dip}} = - \nabla \Upsilon = - \nabla \left(k_m \frac{\sin (\theta) \cos (\phi)}{r^2} \right). \end{equation} \tag{ 34 }$$

As in Batygin & Stevenson (2010), we assume that the induction term is dominated by the interaction with the background field, rather than the induced field: $(\tilde{\mathbf {v}} \times \mathbf {B}) \simeq (\tilde{\mathbf {v}} \times \mathbf {B}_{\rm {dip}})$. Upon making this simplification and uncurling Equation (10), the steady state induction equation reduces to Ohm's law:

$$\begin{equation} \mathbf {J} = \sigma \left(\tilde{\mathbf {v}} \times \nabla \Upsilon - \nabla \Phi \right), \end{equation} \tag{ 35 }$$

where ∇Φ is the electric field.

Because the current is necessarily divergence-free, the scalar potential Φ can be obtained from the following equation:

$$\begin{equation} \nabla ^2 \Phi = \nabla \cdot \left(\tilde{\mathbf {v}} \times \nabla \Upsilon \right) = \frac{k_m \tilde{v}_0}{r^3} \sin (\theta) \sin (\phi). \end{equation} \tag{ 36 }$$

It can be easily checked that the angular part of this relationship is satisfied by

$$\begin{equation} \Phi = \mathcal {A}(r) \sin (\theta) \sin (\phi). \end{equation} \tag{ 37 }$$

As before, confining the current to the atmosphere implies the boundary conditions: $k_m \tilde{v}_0 = \mathcal {R}^3 \mathcal {A}(\mathcal {R}) = (\mathcal {R}+\mathcal {H})^3 \mathcal {A}(\mathcal {R}+\mathcal {H})$. In turn, the radial part of the solution reads

$$\begin{eqnarray} \mathcal {A}(r) &=& k_m \tilde{v}_0((-\mathcal {H} (\mathcal {H}^2+3 \mathcal {H} \mathcal {R}+3 \mathcal {R}^2) (\log (r)+2) \nonumber \\ &&+ \log (\mathcal {R}) (\mathcal {H}^3+3 \mathcal {H}^2 \mathcal {R}+3 \mathcal {H} \mathcal {R}^2+\mathcal {R}^3+2 r^3) \nonumber \\ &&-(\mathcal {R}^3+2 r^3) \log (\mathcal {H}+\mathcal {R})))/(3 \mathcal {H} r^2 \nonumber \\ &&\times (\mathcal {H}^2+3 \mathcal {H} \mathcal {R}+3 \mathcal {R}^2)). \end{eqnarray} \tag{ 38 }$$

The induced current density can now be obtained through Ohm's law.

The Lorentz force can be approximated as originating from the interactions between the induced current and the background magnetic field. The resulting expression is quite cumbersome. However, all of the important features of F_L can be seen by evaluating it at the center of the dynamic domain. To leading order in χ, the expression takes the form:

$$\begin{eqnarray} \mathbf {F}_{\rm {L}}|_{(r = \mathcal {R}+\mathcal {H}/2)}& \simeq & \left(\frac{\sigma k_m^2 }{\mathcal {R}^6 \rho } \right) \tilde{v}_0 \nonumber \\ &&\times \left [ \hat{r} \cos (\theta) (1 -2 \sin ^2(\theta) \cos (2 \phi)+ \cos (2 \theta))/2 \right. \nonumber \\ &&+ \hat{\mathbf {\theta }} (\sin (3 \theta)-\sin (\theta)) \cos ^2(\phi) \nonumber \\ &&\left. - 2 \hat{\mathbf {\phi }} \sin (\theta) \cos (\theta) \sin (\phi) \cos (\phi) \right ]. \end{eqnarray} \tag{ 39 }$$

A similar evaluation of F_L at $r = \mathcal {R}$ and $r = \mathcal {R}+\mathcal {H}$ shows that the $\hat{r}$ and $\hat{\phi }$ components of the force do not change significantly with radius, although the $\hat{\theta }$ component does.

Indeed, the Lorentz force that arises from the interactions between the dayside-to-nightside circulation and the background magnetic field does not only oppose the flow. Instead, it acts to introduce both radial and zonal components to the circulation. Importantly, the typical magnitude of the zonal component of F_L is commensurate with the meridional component (although of course their spatial dependence is different). As argued in Section 3.1, the characteristic timescale associated with the Lorentz force is comparable to the radiative timescale at mbar pressures and is generally shorter than that, corresponding to other relevant forces. This means that the force-balance implied by Equations (17) is in essence not relevant to circulation on hot planetary atmospheres.

In summary, we conclude that dayside-to-nightside flow is unstable to perturbations arising from the Lorentz force. Consequently, we expect that the upper atmospheric circulation will change qualitatively once a substantial magnetic field is introduced into the system. We now turn our attention to numerical MHD simulations with the aim of testing this presumption and quantifying the dynamical state of magnetized upper atmospheres of hot Jupiters.

The hydrodynamic simulation parameters are chosen as described in Section 2, corresponding to the P = 1 mbar pressure level (i.e., τ_L = 10³ s). We start out with the equilibrated hydrodynamic flow shown in panel (A) of Figure 3 and introduce a weak pole-aligned dipole magnetic field into the system. Upon equilibration, we take the approach of sequentially increasing the magnitude of B_dip. At each step, we allow the flow to reach a steady state before increasing the field strength further. We have checked that the flows obtained by successive enhancement of B_dip are identical to those obtained by initializing the atmosphere in solid-body rotation with a given value of B_dip. Consequently, in agreement with Liu & Showman (2013), we conclude that the obtained flows are insensitive to initial conditions.

The panels of Figure 6 show the upper atmospheric circulation for a series of magnetic field strengths. From this series of results, a clear pattern emerges: as the magnitude of B_dip is increased, the flow takes on an exclusively zonal character. Specifically, it is clear that the circulation patterns characteristic of |B_dip| = 0.33 G (panel (B)) are already markedly different from the |B_dip| = 0.025 G case (panel (A)), which clearly resembles the unmagnetized circulation. The flow is in essence entirely azimuthal once the field is increased to |B_dip| = 0.5 G (panel (C)). This is in contrast to the non-uniformly dragged simulations of Rauscher & Menou (2013), who find the flow to become less zonally dominated with enhanced field strength.

Zoom In Zoom Out Reset image size

It is noteworthy that the flow speeds up once it takes on a zonal nature. This is almost certainly due to the fact that viscosity acting on vertical motion more strongly affects the divergent flow, and is therefore not a physically significant result. Increasing the field strength further diminished the flow velocities but did not alter the qualitative nature of the solution, although somewhat higher values of the Ekman number were required to ensure numerical stability.

At the expense of a great inflation in the required computational time, we have extended the simulations presented in Figure 6 to higher resolution. Namely, we prolonged the spherical harmonic decomposition up to degree ℓ_max = 34 and m_max = 29 while eliminating hyperviscosity from our runs entirely. Aside from a mild (i.e., few percent) increase in the velocities, the results observed in these simulations were largely unchanged from the nominal simulations. This implies that the presented solutions do not depend sensitively on small-scale flows. In other words, the transition of the atmosphere to a state dominated by zonal jets is a result of interactions between global circulation and the large-scale magnetic field.

Provided the zonal nature of the flow observed in the magnetized simulation, we can expect that the induced field will be almost entirely toroidal and will approximately be described by Equations (29) and (30). As shown in Figure 5, this indeed appears to be the case. The numerically obtained azimuthal component of the field (panel (B)) is qualitatively similar to its analytically computed counterpart (panel (A)), although the field lines are concentrated toward the vicinity of the equator in the numerical solution (this is simply a consequence of the fact that the circulation is not exactly given by the expression (27)). The magnitude of the induced field is also in good agreement with the analytical theory. For $\tilde{v}_0 = 440$ m s⁻¹, and $k_m/\mathcal {R}^3 = 0.5$ G, Equation (29) yields max(B_ind) = 0.64 G, where as the numerically computed value is max(B_ind) = 0.52 G.

Unlike the case considered in the previous section (where the Lorentz force is treated as a drag), within the framework of MHD simulations, the relationship between the peak wind speed and the temperature perturbation is not necessarily simple. Consequently, in order to preliminarily explore the sensitivity of our results on irradiation, we performed an additional suite of simulations where the applied heat flux was enhanced by a factor of three, compared to the simulations shown in Figure 6. In these overdriven simulations, we found the peak wind speeds to be a factor of ∼2–2.5 higher. However, the characteristic flow patterns closely resembled those, shown in Figure 6, and specifically, the circulation with |B_dip| = 0.5 G remained dominated by azimuthal jets. Consequently, we conclude that the transition of the circulation to a zonal regime with increasing magnetic field strength is robust within the context of our model.

It is interesting to note that the zonal nature of the circulation is ensured at a comparatively low magnetic field strength. If we adopt a scaling based on an Elsasser number of the order of unity (Stevenson 2003), typical hot Jupiter magnetic fields should exceed that of Jupiter by a factor of a few, e.g., |B_dip| ∼ 10 G. Moreover, the arguably more physically sensible scaling based on the intrinsic energy flux (Christensen et al. 2009) suggests that typical hot Jupiter fields may be still higher by another factor of ∼5. Cumulatively, this places the critical magnetic field needed for the onset of zonal flows a factor of ∼10–100 below the typical field strengths. As a result, it would be surprising if a more sophisticated treatment of the hydrodynamics and radiative transfer proved the critical field strength to be higher than the typical one. Nevertheless, such simulations should no doubt be performed.