Magnetic eddy viscosity of mean shear flows in two-dimensional magnetohydrodynamics

TL;DR

This paper introduces the concept of magnetic eddy viscosity in high magnetic Reynolds number regimes of 2D magnetohydrodynamics, showing it causes angular momentum transport and influences zonal flows in astrophysical objects.

Contribution

It derives and confirms the magnetic eddy viscosity in 2D MHD flows, extending the understanding of magnetic effects in high Rm regimes and generalizing Ferraro's law of isorotation.

Findings

Magnetic eddy viscosity arises in high Rm, hydrodynamic-dominated MHD flows.

Numerical simulations confirm the theoretical magnetic viscosity.

The results suggest magnetic viscosity influences angular momentum transport in astrophysical contexts.

Abstract

Magnetic induction in magnetohydrodynamic fluids at magnetic Reynolds number (Rm) less than~1 has long been known to cause magnetic drag. Here, we show that when $\mathrm{Rm} \gg 1$ and the fluid is in a hydrodynamic-dominated regime in which the magnetic energy is much smaller than the kinetic energy, induction due to a mean shear flow leads to a magnetic eddy viscosity. The magnetic viscosity is derived from simple physical arguments, where a coherent response due to shear flow builds up in the magnetic field until decorrelated by turbulent motion. The dynamic viscosity coefficient is approximately $(B_p^2/2\mu_0) \tau_{\rm corr}$, the poloidal magnetic energy density multiplied by the correlation time. We confirm the magnetic eddy viscosity through numerical simulations of two-dimensional incompressible magnetohydrodynamics. We also consider the three-dimensional case, and in cylindrical or spherical geometry, theoretical considerations similarly point to a nonzero viscosity whenever there is differential rotation. Hence, these results serve as a dynamical generalization of Ferraro's law of isorotation. The magnetic eddy viscosity leads to transport of angular momentum and may be of importance to zonal flows in astrophysical domains such as the interior of some gas giants.

Figures

Equations

$$\overline{f}\stackrel{\mathrm{def}}{=}\frac{1}{L_x}{\int }_0^{L_x}dxf,$$

$$f^{\prime }\stackrel{\mathrm{def}}{=}f-\overline{f}.$$

$$\mathsf{T}=\frac{\mathrm{B}\mathrm{B}}{{\mu }_0}-\frac{B^2}{2{\mu }_0}\mathsf{I},$$

$$F_x=\frac{\partial T_{xx}}{\partial x}+\frac{\partial T_{yx}}{\partial y}+\frac{\partial T_{zx}}{\partial z}.$$

$${\overline{F}}_x=\frac{\partial {\overline{T}}_{yx}}{\partial y}=\frac{\partial }{\partial y}\frac{\overline{B_x{B}_y}}{{\mu }_0}.$$

$$\frac{\partial \mathrm{B}}{\partial t}=\nabla \times (\mathrm{v}\times \mathrm{B})=\mathrm{B}\cdot \nabla \mathrm{v}-\mathrm{v}\cdot \nabla \mathrm{B},$$

$$\frac{\partial \mathrm{B}}{\partial t}=B_y\frac{\partial U}{\partial y}\hat{\mathrm{x}}-U\frac{\partial \mathrm{B}}{\partial x}.$$

$$\mathrm{B}(\Delta t)$$

$$\displaystyle\quad+\biggl{(}B_{y}-\Delta tU\frac{\partial B_{y}}{\partial x}\biggr{)}{\bm{\mathrm{\hat{y}}}}+O(\Delta t^{2}),$$

$$\overline{B_x(\Delta t)B_y(\Delta t)}=\overline{B_x{B}_y}+\Delta t\overline{B_y^2}\frac{\partial U}{\partial y}$$

$$\displaystyle\qquad\qquad-\Delta t\overline{\biggl{(}UB_{y}\frac{\partial B_{x}}{\partial x}+UB_{x}\frac{\partial B_{y}}{\partial x}\biggr{)}}+O(\Delta t^{2}).$$

$$\Delta \overline{B_x{B}_y}=\Delta t\overline{B_y^2}\frac{\partial U}{\partial y}+O(\Delta t^2),$$

$$\overline{B_x{B}_y}=\alpha \overline{B_y^2}{\tau }_{\mathrm{corr}}\frac{\partial U}{\partial y},$$

$${\mu }_m=\alpha \frac{\overline{B_y^2}}{{\mu }_0}{\tau }_{\mathrm{corr}},$$

$$\rho \overline{v_x^{\prime }{v}_y^{\prime }}=\gamma \rho \overline{v_y^2}{\tau }_{\mathrm{corr}}\frac{\partial U}{\partial y},$$

$${\mu }_e=-\gamma \rho \overline{v_y^2}{\tau }_{\mathrm{corr}}.$$

$$\frac{\partial}{\partial t}(\overline{\rho v_{x}})\bigr{|}_{\rm turb}=\frac{\partial}{\partial y}\left[\left(\alpha\frac{\overline{B_{y}^{2}}}{\mu_{0}}-\gamma\rho\overline{v_{y}^{2}}\right)\tau_{\rm corr}\frac{\partial U}{\partial y}\right].$$

$$\rho \overline{v_x^{\prime }{v}_y^{\prime }}={\gamma }_2\rho \overline{v_y^2}{\tau }_{\mathrm{corr}}\hat{L}\left(\frac{\partial U}{\partial y}\right).$$

$$\frac{\partial \zeta }{\partial t}+\mathrm{v}\cdot \nabla (\zeta +\beta y)=-\mathrm{B}\cdot \nabla {\nabla }^{2}A-\kappa \zeta +\nu {\nabla }^2\zeta +\xi ,$$

$$\frac{\partial A}{\partial t}+\mathrm{v}\cdot \nabla A=\eta {\nabla }^{2}A.$$

$${\mathrm{Rm}}_0\stackrel{\mathrm{def}}{=}\frac{L_0{V}_0}{\eta },$$

$$\Delta \overline{B_x{B}_y}$$

$$\Delta \overline{B_x{B}_z}$$

$$\Delta \overline{B_s{B}_{\varphi }}$$

$$\Delta \overline{B_z{B}_{\varphi }}$$

$$\Delta \overline{B_r{B}_{\varphi }}$$

$$\Delta \overline{B_{\theta }{B}_{\varphi }}$$

$$\Delta \overline{{\mathrm{B}}_p\cdot \hat{\mathrm{t}}{B}_{\varphi }}=s\Delta t\overline{{\mathrm{B}}_p\cdot \hat{\mathrm{t}}{\mathrm{B}}_p\cdot \nabla \Omega }+O(\Delta t^2),$$

$${\mu }_0{\mathrm{J}}_0$$

$${\mu }_0{\mathrm{J}}_1$$

$$\displaystyle\qquad+\frac{\partial}{\partial y}\biggl{(}B_{0y}\frac{\partial U}{\partial y}-U\frac{\partial B_{0x}}{\partial x}\biggr{)}\biggr{]}.$$

$$F_x=-J_z{B}_y.$$

$$F_{0x}$$

$$F_{1x}$$

$${\overline{F}}_{1x}=\Delta t\frac{\partial }{\partial y}\left(\frac{\overline{B_y^2}}{{\mu }_0}\frac{\partial U}{\partial y}\right)+O(\Delta t^2).$$

Full text

I Introduction

The combination of rotation and magnetic fields is common in astrophysical fluid dynamics, found in stellar interiors, gas giant atmospheres, and astrophysical disks. In the outer atmosphere of planets, magnetic fields do not play a significant role; the flow is dominated by rotation, which leads to anisotropic flow in which turbulence, waves, and coherent structures interact. A common feature is the presence of long-lived mean shear flows, called zonal flow, which alternate in latitude. In Jupiter’s atmosphere, for instance, zonal flows are remarkably persistent Vasavada and Showman (2005). However, deeper into gas giant interiors, or in exoplanets, magnetic fields may play an important role and influence the flow.

In resistive magnetohydrodynamics (MHD), it is well known that when the magnetic Reynolds number is less than one, fluid flow across a magnetic field results in magnetic drag Davidson (1995). The magnetic Reynolds number $\mathrm{Rm}=LV/\eta$ measures the relative strength of induction and magnetic diffusion, where $L$ is a characteristic length scale of the system, $V$ is a characteristic fluid velocity, and $\eta$ is the magnetic diffusivity.

The regime $\mathrm{Rm}\gg 1$ is more dynamical and complex. Within this regime we examine a subset of parameter space: we assume large Reynolds number and consequently turbulent flow. We also assume a hydrodynamic-dominated regime where the influence of the magnetic field is small. We characterize this latter restriction with the dimensionless parameter $\mathcal{A}\ll 1$, where $\mathcal{A}=|\bm{\mathrm{J}}\times\bm{\mathrm{B}}|/|\rho\bm{\mathrm{v}}\cdot\bm{\mathrm{\nabla}}\bm{\mathrm{v}}|$ is the ratio of the Lorentz force to the inertial force. If the energy-containing scale lengths of the magnetic field and velocity field are comparable, then equivalently $\mathcal{A}=(B^{2}/\mu_{0})/(\rho v^{2})$, i.e., the ratio of magnetic energy density to kinetic energy density. Here, $\bm{\mathrm{J}}$ is the current density, $\bm{\mathrm{B}}$ the magnetic field, $\bm{\mathrm{v}}$ the fluid velocity, $\rho$ the fluid mass density, and $\mu_{0}$ the permeability of free space.

In this paper, we suggest that within the above regime, the Lorentz force on mean shear flow acts not as a magnetic drag but rather as an effective magnetic eddy viscosity. We demonstrate this result in two-dimensional (2D) incompressible MHD simulations. Under conditions where $\mathcal{A}$ approaches unity, the Lorentz force becomes appreciable. If the mean shear flow is a turbulence-driven zonal flow, the magnetic eddy viscosity can be strong enough to suppress it. Figure 1 offers a schematic diagram showing the parameter space within which the magnetic eddy viscosity is expected in a turbulent flow.

This paper is primarily a study of MHD physics in an idealized setup; however, a potential example where the regime $\mathrm{Rm}\gg 1$ and $\mathcal{A}\ll 1$ could be realized is in the interiors of gas giants. In the outer atmospheres of Jupiter and Saturn, the bulk fluid is electrically insulating, but the conductivity increases rapidly with depth, giving way to a conducting fluid in the interior French et al. (2012); Liu et al. (2008). Recent spacecraft measurements have shed light on long-standing questions regarding the depth of zonal flow in planetary interiors. Gravitometric observations indicate that the zonal flows in Jupiter and Saturn terminate at a depth of approximately 3,000 and 8,500 km, respectively Kaspi et al. (2018); Guillot et al. (2018); Iess et al. (2019); Galanti et al. (2019). These depths are near the transition zone where the electrical conductivity increases rapidly. The greatly enhanced conductivity in the interior leads to a stronger influence of the magnetic field on bulk fluid motion. The coincident depths of the suppression of zonal flow and the rise in conductivity in Jupiter and Saturn have led to the hypothesis that the magnetic field is responsible for terminating differential rotation. Moreover, simulations of planetary dynamo models provide some support for this possibility Heimpel and Gómez Pérez (2011); Duarte et al. (2013); Jones (2014); Dietrich and Jones (2018). While magnetic drag at $\mathrm{Rm}<1$ has been studied for gas giants and hot Jupiters as the mechanism to suppress differential rotation Cao and Stevenson (2017); Perna et al. (2010), there remains a great deal of uncertainty about parameter values. Intermediate regimes may exist in which $\mathrm{Rm}\gg 1$ and $\mathcal{A}\ll 1$.

Idealized settings to study individual effects prove fruitful. For example, certain aspects of fluid physics under the influence of rotation and magnetization can be studied in the 2D $\beta$ plane or 2D spherical surface Tobias et al. (2007, 2011); Constantinou and Parker (2018); Durston and Gilbert (2016); Zinyakov and Petrosyan (2018). A $\beta$ plane is a Cartesian geometrical simplification of a rotating sphere that retains the physics associated with rotation and the latitudinal variation of the Coriolis effect (Pedlosky, 1992). The unmagnetized $\beta$ plane is one of the simplest settings in which coherent zonal flows emerge spontaneously from turbulence and is often used for idealized studies Rhines (1975); Vallis and Maltrud (1993). Such idealizations eliminate many sources of complexity, including thermal convection, variation in the density and electrical conductivity, and true 3D dynamics. Two-dimensional MHD studies have sufficed, however, to demonstrate magnetic suppression of zonal flow Tobias et al. (2007, 2011); Constantinou and Parker (2018). Here, our numerical simulations continue with the 2D magnetized $\beta$ plane, although our physical considerations are quite general. The magnetic eddy viscosity also sheds light on earlier results on the magnetized $\beta$ plane and magnetized 2D spherical surface.

An outline for the rest of this paper is as follows. In Sec. II, we derive the magnetic eddy viscosity of mean shear flow. In Sec. III, we recall the negative eddy viscosity for comparison. In Sec. IV, we describe the setup for our 2D simulations. The simulation results are shown in Sect. V. Section VI provides the 3D generalization of the magnetic eddy viscosity. Section VII applies these results to consider zonal flows deep in planetary interiors. We summarize our results in Sec. VIII.

II Magnetic eddy viscosity

In this section we consider a mean shear flow $\bm{\mathrm{v}}=U(y){\bm{\mathrm{\hat{x}}}}$ and investigate the zonally averaged magnetic force acting back on the flow due to the correlations induced by the flow. The zonal average of a quantity $f$ is defined as

[TABLE]

where $L_{x}$ is the length of the domain in the $x$ direction, assumed to be periodic. The deviation from the zonal mean is denoted by a prime, i.e.,

[TABLE]

The Lorentz force $\bm{\mathrm{F}}=\bm{\mathrm{J}}\times\bm{\mathrm{B}}$ acting on a fluid can be expressed as $\bm{\mathrm{\nabla}}\cdot\bm{\mathsf{T}}$, where the Maxwell stress tensor $\bm{\mathsf{T}}$ is

[TABLE]

and the stress terms associated with the electric field are neglected within the MHD approximation. This formulation of the Lorentz force will be most convenient for our purposes, though we return to consideration of the current in Appendix A. The force per unit volume in the $x$ direction is

[TABLE]

The first term on the right-hand side (RHS) vanishes after a zonal average. In the remainder of this section, we assume a quasi-2D system and set $B_{z}=0$ and $\partial_{z}=0$, so

[TABLE]

We evaluate the coherent effect that a shear flow has on inducing correlations in the magnetic field and in the corresponding backreaction on the shear flow. We assume a regime $\mathrm{Rm}\gg 1$, in which the magnetic diffusivity has minimal effect on energy-containing length scales over short timescales and the magnetic field is frozen into the fluid. The induction equation is then

[TABLE]

where the second equality applies for incompressible flow $\bm{\mathrm{\nabla}}\cdot\bm{\mathrm{v}}=0$, and $\bm{\mathrm{\nabla}}\cdot\bm{\mathrm{B}}=0$ has been used. Then

[TABLE]

Consider the evolution of the magnetic field from some arbitrary initial time $t=0$ over a short increment $\Delta t$. The magnetic field at time $\Delta t$ is

[TABLE]

where the magnetic-field quantities on the RHS are evaluated at $t=0$. At time $\Delta t$, we can evaluate the zonally averaged Maxwell stress as

[TABLE]

The third and fourth terms on the RHS combine to become $-\Delta tU\overline{\partial_{x}(B_{x}B_{y})}$, which is annihilated by the zonal average. We obtain

[TABLE]

where $\Delta\overline{B_{x}B_{y}}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\overline{B_{x}(\Delta t)B_{y}(\Delta t)}-\overline{B_{x}(0)B_{y}(0)}$. We have thus derived a relation between the Maxwell stress and the shear. The periodicity constraint is essential in producing this exact result; the physical systems in which zonal flow is known to occur all have such periodicity. The derivation of Eq. (10) depends on nothing but the standard MHD induction equation at large $\mathrm{Rm}$.

We would like our theoretical description to involve statistically observable quantities, but Eq. (10) involves the quantities $\Delta\overline{B_{x}B_{y}}$ and $\overline{B_{y}^{2}}$ evaluated at arbitrarily selected points in time, and the increment $\Delta t$, which is not a physically observable quantity associated with the system. What we do is neglect terms of $O(\Delta t^{2})$ and then replace $\Delta t$ by the shortest relevant correlation time $\tau_{\rm corr}$. The intuition guiding this choice is that the coherent effect of the velocity shear builds up in the magnetic field and in the Maxwell stress until some process stops it. Various processes can scramble the coherent effect of the shear, including the magnetic diffusion, turbulent inertial effects, and the shearing itself. In turbulent flow, particularly in hydrodynamic regimes where the magnetic forces are too weak to have significant effect, the turbulent eddy turnover time and perhaps the shearing time are natural correlation times. Whether the higher-order terms in $\Delta t$ are negligible in Eq. (10) depends on the specific regime and the various decorrelation mechanisms.

In light of the preceding paragraph, we make a conceptual switch, from interpreting our result in Eq. (10) as a literal initial-value calculation to instead viewing it through the lens of time-averaged statistical observables. With the replacement of $\Delta t$ by $\tau_{\rm corr}$, we interpret both the left- and right-hand sides of Eq. (10) as statistical, time-averaged quantities. The relation between the shear and the magnetic stress becomes

[TABLE]

where we have introduced an order-unity constant $\alpha$. Hence, the net result is a magnetically originated dynamic viscosity of the form

[TABLE]

such that in the zonally averaged momentum equation, the mean shear results in a coherent contribution to the zonally averaged Lorentz force of $\overline{F}_{x}=\partial_{y}(\mu_{m}\partial_{y}U)$.

We stress that the derivation of Eq. (10) requires no assumptions on the spatial structure of the magnetic field. The magnetic field could be straight (i.e., a poloidal field), but it could also have arbitrary spatial structure. $\overline{B_{y}^{2}}$ incorporates the mean magnetic field as well as all magnetic field perturbations.

Figure 2 shows the induction of a magnetic field by a shear flow over a short time $\Delta t$ for two initial field configurations. Also depicted is the decomposition of the $x$ component of the zonally averaged Lorentz force into the tension and pressure force $\bm{\mathrm{J}}\times\bm{\mathrm{B}}=\bm{\mathrm{F}}_{T}+\bm{\mathrm{F}}_{P}$, where $\bm{\mathrm{F}}_{T}=(B^{2}/\mu_{0}){\bm{\mathrm{\hat{b}}}}\cdot\bm{\mathrm{\nabla}}{\bm{\mathrm{\hat{b}}}}$, $\bm{\mathrm{F}}_{P}=-\bm{\mathrm{\nabla}}_{\perp}B^{2}/2\mu_{0}$, ${\bm{\mathrm{\hat{b}}}}=\bm{\mathrm{B}}/B$, and $\bm{\mathrm{\nabla}}_{\perp}=(\bm{\mathsf{I}}-{\bm{\mathrm{\hat{b}}}}{\bm{\mathrm{\hat{b}}}})\cdot\bm{\mathrm{\nabla}}$. When the initial magnetic field lines are straight, $\overline{F}_{x}$ is due to magnetic tension, akin to the behavior of a shear Alfvén wave. However, in the general case with nonstraight field lines, as might be expected to occur in the tangle of magnetic fields when induction is important and $\mathrm{Rm}\gg 1$, $\overline{F}_{x}$ includes contributions from both tension and pressure. In the alternative definition of pressure and tension, where the pressure is isotropic and the tension is given by $\bm{\mathrm{B}}\cdot\bm{\mathrm{\nabla}}\bm{\mathrm{B}}/\mu_{0}$, $\overline{F}_{x}$ is entirely tension. Regardless of the initial field configuration, $\overline{F}_{x}$ still acts as a viscosity, in accordance with Eq. (10).

In the derivation of Eq. (10), we have only imposed the effect of a mean shear flow. Let us consider what would happen if we included the full fluid velocity field in the induction equation, in which the mean shear flow $U\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\overline{\bm{\mathrm{v}}}\cdot{\bm{\mathrm{\hat{x}}}}$ is embedded. In a short time increment, there would be additional contributions to $\Delta\bm{\mathrm{B}}$. Those contributions not correlated across the zonal domain would mostly cancel out when the zonal average of the Maxwell stress is taken. Hence, it is the mean shear flow that gives rise to a coherent zonally averaged effect.

In Sec. V, we compare the prediction Eq. (11) for the Maxwell stress with direct numerical simulations.

III Negative eddy viscosity

For comparison with the magnetic eddy viscosity, it is worthwhile to recall the so-called negative eddy viscosity. The notion of a hydrodynamic negative eddy viscosity in 2D flow has a long history Starr (1953); Kraichnan (1976); Gama et al. (1994); Krommes and Kim (2000). A negative eddy viscosity emerges in a variety of mathematical approaches, which can be seen as a reflection of the general tendency for energy to tend towards large scales in 2D Fjørtoft (1953). The early theoretical literature was primarily concerned with the context of isotropic turbulent flows. Recent work, however, has revived the negative eddy viscosity to help explain and interpret the physics of the formation and maintenance of zonal flow.

For example, the zonostrophic instability supposes a statistically homogeneous turbulent background and considers the growth of a perturbing zonal flow Farrell and Ioannou (2007); Srinivasan and Young (2012); Parker and Krommes (2013, 2014); Constantinou et al. (2014); Fitzgerald and Farrell (2018). For a weak zonal flow of wavenumber $q$ and turbulent background with spectrum $\Phi(\bm{\mathrm{k}})$, the flow reinforces the perturbation at a rate that schematically looks like $\int d\bm{\mathrm{k}}\,q^{2}(1-q^{2}/k^{2})g(\bm{\mathrm{k}})\Phi(\bm{\mathrm{k}})$, where $g$ contains anisotropic wave-vector dependence. At small $q$, this goes like $q^{2}$. At larger $q$, there is a cutoff of the form $1-q^{2}/k^{2}$, which can be traced back to the $\bm{\mathrm{v}}^{\prime}\cdot\bm{\mathrm{\nabla}}\overline{\zeta}$ term, where $\zeta$ is the vorticity Bakas and Ioannou (2013); Parker (2016). Many treatments have neglected this term compared to $\overline{\bm{\mathrm{v}}}\cdot\bm{\mathrm{\nabla}}\zeta^{\prime}$ under the assumption of long-wavelength mean flows. We return to discussion of this cutoff at the end of the section.

Manifestations of a negative eddy viscosity are visible at cloud level in both Jupiter and Saturn. Careful measurements have demonstrated that the Reynolds stress of observed cloud-level velocity fluctuations are positively correlated with the mean shear Salyk et al. (2006); Del Genio et al. (2007); Del Genio and Barbara (2012). Hence, the turbulent fluctuations appear to transfer energy to the mean flow via Reynolds stress. Measurements of the turbulent Reynolds stress in laboratory plasma devices lead to similar conclusions Holland et al. (2006); Tynan et al. (2006).

It is instructive to calculate the Reynolds stress in a similar manner as we did the Maxwell stress. The calculation also serves as yet another simple demonstration of a negative eddy viscosity. The details are given in Appendix B. We obtain a prediction for the Reynolds stress,

[TABLE]

where $\gamma$ is a positive constant of order unity. This corresponds to a negative dynamic eddy viscosity,

[TABLE]

The dynamic eddy viscosity has magnitude approximately equal to the kinetic energy density multiplied by a correlation time. In pure hydrodynamic flow, the eddy turnover time at the energy-containing scale, $\tau_{\rm turb}$, is often the appropriate correlation time. Using $v_{\rm turb}\sim L_{\rm turb}/\tau_{\rm turb}$ leads to a mixing-length estimate for the kinematic negative eddy viscosity acting on mean shear flow, $\nu_{e}=-\gamma L_{\rm turb}^{2}/\tau_{\rm turb}$.

Assuming the relevant correlation times in Eqs. (11) and (13) are equal, the zonally averaged momentum flux due to the turbulent Reynolds and Maxwell stresses is

[TABLE]

Hence, the magnetic eddy viscosity becomes comparable to the negative eddy viscosity when the fluctuating magnetic energy has reached a level roughly similar to that of the fluctuating kinetic energy. However, in Sec. V we will see from our 2D simulations that the numerical prefactor $\gamma$ can be an order of magnitude smaller than $\alpha$. Therefore, if the turbulent Reynolds stress is the primary driver of zonal flow, the magnetic viscosity can have a significant influence even when the fluctuating magnetic energy is much smaller.

A refinement of Eq. (13) is possible. The aforementioned $1-q^{2}/k^{2}$ cutoff in the growth rate motivates a simple modification of Eq. (13),

[TABLE]

Here, $\hat{L}$ is an operator that multiplies the Fourier amplitude by $1-q^{2}/k_{c}^{2}$ for $q\leq k_{c}$, and suppresses the Fourier amplitude for $q>k_{c}$, where $k_{c}$ is a characteristic eddy wave number. This refinement is additionally motivated by the fact that modeling attempts that use only the local shear and neglect higher-order derivatives lead to divergent growth of small scales, but proper inclusion of the cutoff leads to well-behaved models Parker (2016).

IV Simulation Setup

We perform a series of numerical simulations to assess our prediction for the Maxwell stress and Reynolds stress. We study the behavior of a 2D incompressible MHD fluid on a magnetized $\beta$ plane Tobias et al. (2007).

The coordinates on the $\beta$ plane are $\bm{\mathrm{x}}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}(x,y)$, where $x$ is the azimuthal direction (longitude) and $y$ the meridional direction (latitude). It is convenient to represent both the velocity and magnetic field by potentials. The velocity in the rotating frame is written as $\bm{\mathrm{v}}={\bm{\mathrm{\hat{z}}}}\times\bm{\mathrm{\nabla}}\psi$, where $\psi$ is the stream function. The vorticity normal to the plane is $\zeta\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}{\bm{\mathrm{\hat{z}}}}\cdot(\bm{\mathrm{\nabla}}\times\bm{\mathrm{v}})=\nabla^{2}\psi$. The magnetic field $\bm{\mathrm{B}}=\bm{\mathrm{B}}_{0}+\bm{\mathrm{\tilde{B}}}$ consists of a constant, uniform background field $\bm{\mathrm{B}}_{0}$ and a time-varying component $\bm{\mathrm{\tilde{B}}}(\bm{\mathrm{x}},t)$. We define a corresponding vector potential $\bm{\mathrm{A}}=A{\bm{\mathrm{\hat{z}}}}$ such that $\bm{\mathrm{B}}=\bm{\mathrm{\nabla}}\times\bm{\mathrm{A}}$, where $A=A_{0}+\tilde{A}$ and $A_{0}=-B_{0y}x+B_{0x}y$. As we have imposed the presence of a background magnetic field, this 2D system is not subject to antidynamo theorems, and magnetic amplification can result such that $B\gg B_{0}$. We perform simulations in which the background magnetic field is either purely toroidal (i.e., aligned in the azimuthal direction, $B_{0y}=0$) or purely poloidal (i.e., aligned in the latitudinal direction, $B_{0x}=0$).

The evolution of the system can be described by a formulation involving vorticity and magnetic potential,

[TABLE]

Here, $\beta$ is the latitudinal gradient of the Coriolis parameter, $\kappa$ is the linear drag, $\nu$ is the viscosity, $\eta$ is the magnetic diffusivity, and $\xi(\bm{\mathrm{x}},t)$ is a mechanical forcing to excite hydrodynamic fluctuations. The first term on the RHS of Eq. (17a) stems from the curl of the Lorentz force. Equation (17b) is the magnetic induction equation written for the potential $A$. For mathematical convenience, we have set the permeability $\mu_{0}=1$ and the mass density $\rho=1$.

We numerically solve Eq. (17) using the code Dedalus Burns et al. (2019). We employ periodic boundary conditions in both directions. During time stepping, we use a high-wave-number spectral filter to remove enstrophy from accumulating at the grid scale. Unless noted, our domain size is $L_{x}\times L_{y}=4\pi\times 4\pi$ at a resolution of $512^{2}$ gridpoints. The forcing $\xi$ is statistically homogeneous and isotropic white noise. It is also small scale, localized near wave vectors with magnitude $k_{f}=12$, and injects energy per unit area at a rate $10^{-3}$ (see Ref. Constantinou and Parker (2018) for more details).

The majority of our simulations use nonzero drag $\kappa$ on the hydrodynamic flow in order to saturate without filling the largest scales of the box. However, we also show a case with only viscosity, $\kappa=0$, to confirm our results are not sensitive to the presence of drag.

V Simulation Results

We perform two series of simulations: one with a background toroidal magnetic field $B_{0}=10^{-2}$, and one with a background poloidal magnetic field $B_{0}=10^{-3}$. These runs use $\kappa=10^{-2}$, $\nu=10^{-4}$, and $\beta=10$. For each simulation series, we vary $\eta$. It is then convenient to define another magnetic Reynolds number that is referenced to the hydrodynamic limit,

[TABLE]

where $L_{0}$ and $V_{0}$ are a characteristic length scale and velocity scale in a regime where the magnetic field has no influence, i.e., large $\eta$. For $L_{0}$ we use a typical eddy length scale computed from the eddy energy spectrum as $L_{\rm turb}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}[\int d\bm{\mathrm{k}}\,k^{2}E(\bm{\mathrm{k}})/\int d\bm{\mathrm{k}}\,E(\bm{\mathrm{k}})]^{-1/2}$. For $V_{0}$ we use the rms eddy velocity $v_{\rm turb}$. Our simulations yield $L_{\rm turb}=0.14$ and $v_{\rm turb}=0.21$, giving $\mathrm{Rm}_{0}=0.03/\eta$. Our reported values of $\mathcal{A}$ are computed as the total magnetic energy divided by the total kinetic energy.

As we decrease $\eta$ from large values, we move through the three regimes depicted in Figure 1. First, is the regime of large $\eta$, $\mathrm{Rm}_{0}<1$, and magnetic energy much smaller than kinetic energy, $\mathcal{A}\ll 1$. As $\eta$ is reduced, a second regime is entered, where $\mathrm{Rm}_{0}>1$, and still $\mathcal{A}\ll 1$. Then, as $\eta$ is decreased even further, $\mathrm{Rm}_{0}\gg 1$ and $\mathcal{A}\approx 1$. The intermediate regime is the one where the assumptions behind the magnetic eddy viscosity are valid.

Figure 3 shows results from three representative simulations. The left column shows space-time diagrams of the zonal flow $U(y,t)=\overline{v}_{x}$. The right column shows time series of the different contributions to the domain-averaged energy density: zonal-flow kinetic energy (ZKE) $\langle U^{2}/2\rangle$, eddy kinetic energy (EKE) $\langle(v_{x}^{\prime 2}+v_{y}^{\prime 2})/2\rangle$, zonal magnetic energy (ZME) \bigl{\langle}\overline{\tilde{B}}_{x}^{2}/2\bigr{\rangle}, eddy magnetic energy (EME) $\langle(B_{x}^{\prime 2}+B_{y}^{\prime 2})/2\rangle$, and the background magnetic energy (BME) $B_{0}^{2}/2$. Here, $\langle f\rangle\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\int d\bm{\mathrm{x}}\,f/L_{x}L_{y}$ denotes a domain average of quantity $f$.

Figures 3(a) and 3(b) show a case with $\mathrm{Rm}_{0}=0.3$. Figure 3(a) shows that zonal jets form spontaneously. After an initial transient time of about $5/\kappa=500$, seven jets equilibrate at finite amplitude and remain coherent throughout the rest of the simulation. The eddy magnetic energy is about 4 orders of magnitude smaller than the kinetic energy of the flow [Fig. 3(b)]. Figures 3(c) and 3(d) show a case with toroidal background field at $\mathrm{Rm}_{0}=300$. The magnetic energy is somewhat smaller than the kinetic energy, $\mathcal{A}=0.05$. The zonal flow is mostly suppressed; it has smaller amplitude and is less coherent. Figures 3(3) and 3(f) show a case with poloidal background field at $\mathrm{Rm}_{0}=30$. Here, the zonal kinetic energy is not reduced from the hydrodynamic case as in Fig. 3(a), but the coherence of the zonal flow is clearly reduced in Fig. 3(e). There is a continual cycle of jet merging and formation of new jets, and the zonal flow never settles down to a steady state. Interestingly, this behavior occurs even when $\mathcal{A}=0.006$.

V.1 Maxwell stress

We compare the Maxwell stress predicted by Eq. (11) with the actual stress observed in the simulations. To facilitate the comparison, we perform a least-squares best fit for the single free parameter $\alpha$, which serves merely as an overall constant of proportionality and does not affect the spatial structure. The fit is performed on the Maxwell stress itself, not its derivative.

Figure 4 examines a simulation with toroidal background field. For this simulation, $\mathrm{Rm}_{0}=30$ and $\mathcal{A}=0.005$. The curves labeled “predicted $\overline{B_{x}^{\prime}B_{y}^{\prime}}$” are computed using correlation time $\tau_{\rm corr}=L_{\rm turb}/v_{\rm turb}$. Both short-duration averages (50 time units, left column) and long-duration averages (500 time units, right column) are shown for this comparison to illustrate that there is no substantial difference. The proportionality constant $\alpha=1.02$ for the short-time average and 1.04 for the long-time average. The predictions for the Maxwell stress from Eq. (11) are in very good agreement with those from the simulations. $\overline{B_{y}^{\prime 2}}$ does not vary much with $y$ [Figs. 4(a) and 4(b)], and therefore the Maxwell stress closely resembles the mean flow shear; compare Figs. 4(c) and 4(d) with Figs. 4(e) and 4(f). Additionally, by comparing the spatial dependence of $\partial_{y}\overline{B_{x}^{\prime}B_{y}^{\prime}}$ [Figs. 4(i) and 4(j)] with $-U$ and $\partial_{y}^{2}U$ [Figs. 4(g) and 4(h)], we see that $\partial_{y}\overline{B^{\prime}_{x}B^{\prime}_{y}}$ more closely resembles $\partial_{y}^{2}U$. Thus, the Lorentz force acts like viscosity rather than like drag.

Figure 5 is much the same, but for a poloidal background field, with $\mathrm{Rm}_{0}=30$ and $\mathcal{A}=0.006$. There is an additional subtlety here regarding the calculation of the Maxwell stress. The prediction in Sec. II uses the total magnetic field $\bm{\mathrm{B}}$. For a poloidal background field, $\overline{B_{x}B_{y}}=\overline{B}_{x}B_{0y}+\overline{B_{x}^{\prime}B_{y}^{\prime}}$ (whereas for a toroidal background field, $\overline{B_{x}B_{y}}=\overline{B_{x}^{\prime}B_{y}^{\prime}}$). However, Eq. (11) can be used as a prediction for $\overline{B_{x}^{\prime}B_{y}^{\prime}}$ if we replace $\overline{B_{y}^{2}}$ with $\overline{B_{y}^{\prime 2}}$. The figure shows the prediction for $\overline{B_{x}^{\prime}B_{y}^{\prime}}$ rather than $\overline{B_{x}B_{y}}$, though this distinction has little quantitative significance here because $B_{0}\ll\tilde{B}$ for $\mathrm{Rm}\gg 1$.

In Fig. 5, the prediction for the Maxwell stress agrees quite well with the actual values. The proportionality constant $\alpha=1.05$ for the short-time average and 1.17 for the long-time average. It is interesting to note that unlike the toroidal case, $\overline{B_{y}^{\prime 2}}$ has significant variation in $y$. This dependence leaves a noticeable imprint in the Maxwell stress, showing that $\overline{B_{x}^{\prime}B_{y}^{\prime}}$ is not merely proportional to the flow shear $\partial_{y}U$. Figures 5(i) and 5(j) show, in addition to the derivative of the fluctuating stresses, the term $\partial_{y}(\overline{B}_{x}B_{0y})$. This term, which is responsible for magnetic drag at $\mathrm{Rm}\ll 1$, is much smaller than $\partial_{y}\overline{B^{\prime}_{x}B^{\prime}_{y}}$ here.

Figure 6 shows simulations results at different values of $\mathrm{Rm}_{0}$. The left column shows results in the statistically steady state for a toroidal background field. The regime of magnetic eddy viscosity is valid for $3\lesssim\mathrm{Rm}_{0}\lesssim 300$. The upper bound is determined where the magnetic field has a significant influence on the flow and the condition $\mathcal{A}\ll 1$ no longer holds. Figure 6(e) shows the correlation coefficient between the derivative of the Maxwell stress and the prediction, which is close to 1 in the magnetic eddy viscosity regime. Figure 6(g) shows the rms values of the Reynolds stress and Maxwell stress, which exhibit a similar dependence as the kinetic energy and magnetic energy. The proportionality constant $\alpha$ for both the toroidal and poloidal cases is approximately equal to 1, as seen in Figs. 6(i) and 6(j).

The right column depicts the case with a poloidal background field. The details are similar, with the regime of magnetic eddy viscosity existing for $10\lesssim\mathrm{Rm}_{0}\lesssim 100$. Figure 6 also demonstrates that the Lorentz force transitions from a magnetic drag at low $\mathrm{Rm}_{0}$ to magnetic viscosity at high $\mathrm{Rm}_{0}$. Figure 6(f) shows that at low $\mathrm{Rm}_{0}$, the correlation of $\partial_{y}(\overline{B}_{x}B_{0y})$ with $-U$ is equal to 1, indicating drag. A straightforward calculation shows the amplitude of the magnetic drag at low $\mathrm{Rm}_{0}$ is expected to be $\overline{F}_{x}=-(B_{0}^{2}/\eta\mu_{0})U$ (momentarily returning to dimensional quantities); this amplitude is quantitatively recovered in the simulation. As $\mathrm{Rm}_{0}$ increases, the Lorentz force due to fluctuating magnetic fields, $\partial_{y}\overline{B_{x}^{\prime}B_{y}^{\prime}}$, correlates well with the predicted magnetic viscosity. When $\mathrm{Rm}_{0}\gtrsim 10$, the viscouslike fluctuating force begins to dominate, as seen in Fig. 6(h).

Even though the kinetic energy is much greater than the magnetic energy for $\mathrm{Rm}_{0}<300$, Fig. 6(d) shows how a poloidal background field affects the coherence of the zonal flow even starting at $\mathrm{Rm}_{0}=3$. Figure 6(d) is computed using the time average $\langle U\rangle_{t}^{2}/2$, which is reduced compared to the zonal kinetic energy when the spatial structure of the zonal flow is not steady in time [see Figs. 3(a) and 3(e)]. The longer the duration of the time average, the more $\langle U\rangle_{t}^{2}$ will be reduced by the decoherence.

To check whether the prediction for the Maxwell stress is sensitive to the presence of frictional drag, we also perform a simulation with $\kappa=0$ but maintaining nonzero viscosity. The results are shown in Fig. 7. This simulation uses domain size $L_{x}\times L_{y}=2\pi\times 2\pi$, numerical resolution of $256^{2}$, and parameter values $\beta=2$, $\nu=10^{-4}$, and $\eta=10^{-4}$. The background magnetic field is toroidal, with $B_{0}=10^{-2}$. Without drag, the structures tend to fill the box size, and only a single jet exists at late times. Time-averaging around $t=2500$, we find $L_{\rm turb}=0.56$, $v_{\rm turb}=0.46$, and hence $\mathrm{Rm}\approx 2600$. The magnetic energy is much smaller than the kinetic energy, with $\mathcal{A}\approx 6\times 10^{-4}$.

Equation (11) accurately predicts the Maxwell stress in this simulation as well. This is true both at early times $200<t<350$, when there are two jets, and at late times, $3000<t<3250$, after the two jets have merged so that a single wavelength fills the domain. Moreover, at late times, the spatial structure of $U$ and $\partial_{y}^{2}U$ are significantly different, and $\partial_{y}\overline{B_{x}^{\prime}B_{y}^{\prime}}$ is clearly closer to $\partial_{y}^{2}U$. We also observe that Eq. (11) holds even though the kinetic energy has not yet saturated under viscous dissipation, because the timescale setting the magnetic quasi-equilibrium is much shorter than the viscous time.

V.2 Reynolds stress

We also examine the negative eddy viscosity in our simulations to determine whether the relation in Eq. (13) holds in our simulations. Figure 8 shows results for a simulation with a background toroidal field of $B_{0}=10^{-2}$. The system has $\mathcal{A}\ll 1$ and so the magnetic field has little influence on the flow. Correlations induced in the fluctuating Reynolds stress by the zonal flow are picked out by zonal averaging. The zonally averaged Reynolds stress and its divergence are in overall good agreement with the predictions from Eq. (13), with a best-fit $\gamma=0.08$. Observe that the proportionality constant for the Reynolds stress $\gamma$ is about 10 times smaller than the proportionality constant for the Maxwell stress $\alpha$ [see Figs. 6(i) and 6(j)].

Equation (13) produces a noticeable over-prediction of the jet forcing $-\partial_{y}\overline{v_{x}^{\prime}v_{y}^{\prime}}$ at the eastward peaks of the jets. Figure 8 also shows the refined prediction of Eq. (16). The refined prediction removes the over-estimate of $-\partial_{y}\overline{v_{x}^{\prime}v_{y}^{\prime}}$ at the eastward jets and shows much better agreement with the actual value. We have used $k_{c}=1/L_{\rm turb}=7.1$, so the only free parameter is the overall proportionality factor which, for this case, is $\gamma_{2}=0.13$.

VI Generalization to three dimensions

The simple derivation of magnetic eddy viscosity in Sec. II assumed a 2D Cartesian system. Here, we follow a similar procedure to generalize Eq. (10) to three dimensions in Cartesian, cylindrical, and spherical coordinates. We find the relevant components of the Maxwell stress tensor that are required to compute the zonally averaged Lorentz force. As in Sec. II, the following equations are exact to first order in $\Delta t$ and depend only on the large-$\mathrm{Rm}$ limit of the MHD induction equation.

In Cartesian coordinates, let $\bm{\mathrm{v}}=U(y,z){\bm{\mathrm{\hat{x}}}}$. A zonal average is over $x$, and periodicity in $x$ assumed. We obtain

[TABLE]

In cylindrical coordinates $(s,\phi,z)$, let $\bm{\mathrm{v}}=U(s,z)\bm{\hat{\phi}}$, and a zonal average is over the azimuthal angle $\phi$. We find

[TABLE]

In spherical coordinates ($r,\theta,\phi)$, let $\bm{\mathrm{v}}=U(r,\theta)\bm{\hat{\phi}}$, and again a zonal average is over $\phi$. We obtain

[TABLE]

We make several observations. First, in cylindrical or spherical geometry, the first-order coherent effect in the Maxwell stress due to $U$ vanishes if $U\sim s\sim r\sin\theta$, i.e., solid body rotation. Second, one can straightforwardly compute that in the cross terms $\overline{B_{y}B_{z}}$, $\overline{B_{s}B_{z}}$, or $\overline{B_{r}B_{\theta}}$, there is no first-order coherent effect due to $U$. One might expect that the fluctuating field’s contribution to this term averages to approximately zero. A mean magnetic field could lead to nonzero correlation such as $\overline{B_{r}B_{\theta}}$, although in a regime $\mathrm{Rm}\gg 1$ the $\overline{B_{r}^{2}}$ and $\overline{B_{\theta}^{2}}$ terms likely dominate. Third, the computation of the zonally averaged Lorentz force in non-Cartesian geometries, $\overline{F}_{\phi}=\overline{\bm{\mathrm{\nabla}}\cdot\bm{\mathsf{T}}\cdot\bm{\hat{\phi}}}$, introduces additional geometric factors in the tensor divergence. Fourth, for rotation about an axis, $U=s\Omega$, Eqs. (20) and (21) may be written compactly as

[TABLE]

where $\bm{\mathrm{B}}_{p}$ is the poloidal (nonazimuthal) part of the magnetic field and ${\bm{\mathrm{\hat{t}}}}$ is a unit vector in the plane orthogonal to ${\bm{\mathrm{\hat{\phi}}}}$, e.g., ${\bm{\mathrm{\hat{t}}}}={\bm{\mathrm{\hat{s}}}},{\bm{\mathrm{\hat{z}}}},{\bm{\mathrm{\hat{r}}}},\bm{\hat{\theta}}$.

Equation (22) functions as a dynamical generalization of Ferraro’s law of isorotation, which states that for an axisymmetric magnetic field and axisymmetric azimuthal flow with angular velocity $\Omega(s,z)$, a steady state is possible only when $\bm{\mathrm{B}}_{p}\cdot\bm{\mathrm{\nabla}}\Omega=0$ Ferraro (1937). Hence, Ferraro’s law has been cited for the notion that magnetic fields tend to suppress differential rotation (although this is not required; the law merely requires that $\bm{\mathrm{B}}_{p}$ lies on surfaces of constant $\Omega$). Our derivation of a magnetic eddy viscosity generalizes the isorotation law because it (i) allows for turbulent dynamics, (ii) does not assume axisymmetry of the magnetic field, and (iii) provides an estimate for the magnitude of the resulting force that counteracts differential rotation. The isorotation law is recovered from Eq. (22) if $\bm{\mathrm{B}}_{p}$ is axisymmetric.

Under appropriate physical conditions where second-order effects can be neglected, one can apply the same arguments as in Sec. II to relate the coherent part of the statistically averaged Maxwell stresses to the mean shear flow. If the first-order truncation is applicable, then it is clear that a differentially rotating body will induce coherent Maxwell stresses that act to redistribute angular momentum in a way similar to a viscosity.

While the prediction of the Maxwell stress was verified in 2D MHD $\beta$-plane simulations in Sec. V, comparing the equations of this section with 3D numerical simulations is out of scope of the present work. The expressions in Eqs. (19)–(21) could be directly compared with the Maxwell stress observed in 3D simulations to determine the applicability of the magnetic eddy viscosity in a particular regime.

VII Estimates for Jupiter and Saturn

Could the magnetic eddy viscosity regime apply in some planetary interiors and perhaps be responsible for suppression of differential rotation? We consider this question along with some Jupiter and Saturn data. We explore the possibility that $\mathrm{Rm}\gg 1$ occurs in regions with significant mean flow (implicitly with $\mathcal{A}\ll 1$). This analysis is restricted to the upper atmosphere as the deeper planetary dynamo region likely has $\mathcal{A}\sim O(1)$.

There is currently a great deal of uncertainty about the dynamics below the cloud-top level in these two gas giants. While the flow speed is expected to decrease with depth due to the increasing mass density, for this estimate we will use the observed outer-atmosphere values.

Cloud tracking of Jupiter’s and Saturn’s upper-level atmosphere showed that eddies act to reinforce the zonal jets with an effective eddy viscosity $\nu_{e,J}\approx-10^{6}\,\rm m^{2}\,s^{-1}$ for Jupiter Salyk et al. (2006) and $\nu_{e,S}\approx-2\times 10^{5}\,\rm m^{2}\,s^{-1}$ for Saturn Del Genio et al. (2007); Del Genio and Barbara (2012). We observe that $\nu_{e,J}$ is roughly equal to the mixing-length estimate $\nu_{e}=-\gamma L_{\rm turb}^{2}/\tau_{\rm turb}$, using $L_{\rm turb}=1000$ km and $\tau_{\rm turb}=2.5\times 10^{5}$ s observed at the Jupiter cloud tops. Supposing the negative eddy viscosity is the sole driving force of the zonal flow, we ask: if the negative eddy viscosity is to be counterbalanced solely by magnetic eddy viscosity, then (i) what is the strength of the magnetic field needed to do so and (ii) how deep below the cloud tops of the two gas giants do we expect such strong magnetic fields?

The mass density increases rapidly beneath the cloud tops. We use a number intermediate within the expected density range for both Jupiter and Saturn, $\rho\approx 200$ kg m*-3* French et al. (2012); Anderson and Schubert (2007); Mankovich et al. (2019). For the correlation time, we use the turbulent eddy turnover time $\tau_{\rm turb}=2.5\times 10^{5}$ s. This time is not that different from the typical shear time associated with the surface-level zonal jets $t_{\textrm{shear},J}\approx 2\times 10^{5}\,\rm s$ and $t_{\textrm{shear},S}\approx 5\times 10^{4}\,\rm s$ Vasavada and Showman (2005); Del Genio et al. (2007). We also use typical mean magnetic field values $B_{0,J}=4\times 10^{-3}\,\rm T$ Moore et al. (2018) and $B_{0,S}=2\times 10^{-5}\,\rm T$ Gombosi et al. (2009). With these in hand, the kinematic magnetic eddy viscosity $\nu_{m}=\mu_{m}/\rho$ computed with the mean magnetic field values is 2–4 orders of magnitude smaller than the corresponding negative eddy viscosity. Thus, magnetic viscosity can only start becoming important if we take into account the magnetic fluctuations.

For convenience, we assume the scaling $B^{2}\sim B_{0}^{2}\text{Rm}$ at $\mathrm{Rm}\gg 1$. Accounting for the factor of $\mathrm{Rm}$, we estimate for Jupiter that $\nu_{m}\approx\nu_{e}$ requires $\mathrm{Rm}_{J}\approx 60$ and for Saturn requires $\mathrm{Rm}_{S}\approx 700$. Hence these results are consistent with the assumption $\mathrm{Rm}\gg 1$. Then, using the outer-atmosphere typical eddy length scale $L_{\rm turb}$ and speed $v_{\rm turb}$, we compute that the required $\mathrm{Rm}$ values correspond to magnetic diffusivities of $\eta_{J}\approx 6\times 10^{4}\,\rm m^{2}\,s^{-1}$ and $\eta_{S}\approx 6\times 10^{3}\,\rm m^{2}\,s^{-1}$. Finally, from Refs. French et al. (2012); Liu et al. (2008); Cao and Stevenson (2017) we find that such magnetic diffusivity values are achieved at approximately $0.96R_{J}$ and $0.8R_{S}$. Therefore, under the assumptions we have made, we expect zonal jets would be magnetically terminated by these depths. These estimates are not inconsistent with the recent gravitometric estimates of depths of $0.96R_{J}$ and $0.85R_{S}$ Kaspi et al. (2018); Guillot et al. (2018); Galanti et al. (2019); Iess et al. (2019).

The decrease of flow velocity with depth will modify these calculations. Suppose that at the critical depth, the characteristic velocity and turbulent length scale have decreased by a factor of 10 relative to the observed surface values. The negative eddy viscosity would be smaller by a factor of 100, using the scaling of Sec. III. Hence, the required $\mathrm{Rm}_{J}$ falls below 1, making it marginal, while the required $\mathrm{Rm}_{S}$ is about 7, still within the regime of magnetic eddy viscosity. The required $\eta_{S}=VL/\mathrm{Rm}_{S}$ is unchanged from the previous estimate.

This discussion is not intended to assert that magnetic eddy viscosity is the mechanism responsible for terminating the zonal flows in Jupiter and Saturn, but rather serves as a consideration of whether it is a plausible mechanism for gas giants at all. There are substantial complexities not included in this discussion, such as the variation of conductivity and density with depth, meridional flow, thermal convection, and realistic physics of planetary dynamos. Still, these estimates suggest that magnetic eddy viscosity is a potential candidate for inhibiting differential rotation in the interior of some gas giants.

VIII Conclusion

In 2D incompressible MHD, under the assumptions of $\mathrm{Rm}\gg 1$ and a turbulent flow regime where the magnetic energy is smaller than the kinetic energy, we have shown that a mean shear flow induces correlations in the magnetic field that result in a magnetic eddy viscosity. Numerical simulations of a 2D magnetized $\beta$ plane, in which zonal flows arise naturally, confirm this result.

Two key ingredients lead directly to the prediction of magnetic eddy viscosity, given in Eqs. (11) and (12). First is the large-$\mathrm{Rm}$ limit of the MHD induction equation. Second is the assumption of turbulent flow, which acts to scramble any coherent effect built up in the magnetic field, by enforcing a short correlation time.

The expressions for the 3D generalization of the components of the Maxwell stress tensor are given in Sec. VI. They show directly how differential rotation leads to zonally averaged Maxwell stress that transports angular momentum. Additionally, these expressions function as a dynamical generalization of Ferraro’s law of isorotation, because they (i) allow for turbulent dynamics, (ii) do not assume axisymmetry of the magnetic field, and (iii) provide an estimate for the magnitude of the resulting force that counteracts differential rotation.

In 2D simulations with a mean poloidal magnetic field, we observe the Lorentz force transition from magnetic drag to magnetic eddy viscosity as $\mathrm{Rm}$ increases beyond about 10. Magnetic eddy viscosity also occurs in simulations where the mean magnetic field is toroidal, in which case the $\omega$ effect is absent and there is no associated magnetic drag; the magnetic eddy viscosity is due solely to magnetic fluctuations.

Rotating, magnetized fluids are commonly found in astrophysical domains. The magnetic eddy viscosity may, in certain parameter regimes, provide a simple characterization of the transport of angular momentum and the magnetic braking force due to differential rotation. Our main numerical example here is a case where the mean shear flow consists of a turbulence-driven zonal flow, with a simple idealized model that could provide insight into the interiors of gas giants. Mean shear flow also occurs in astrophysical disks, and magnetic eddy viscosity could be explored in that context as well.

Acknowledgements.

We gratefully acknowledge useful discussions with Bill Dorland, Eli Galanti, Tristan Guillot, Yohai Kaspi, Brad Marston, Felix Parra, David Stevenson, and Steve Tobias. The work of J.B.P. was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. N.C.C. thanks Sean Haney for his hospitality in Barber Tract in December 2018, where part of this paper was written. Numerical simulations were conducted on the Australian National Computational Infrastructure at ANU, which is supported by the Commonwealth of Australia.

Appendix A Calculation of current and force directly from $\bm{\mathrm{J}}\times\bm{\mathrm{B}}$

Although the Maxwell stress tensor is the most direct route to the force and momentum flux derived in Sec. II, we can also calculate the force by finding the current and using $\bm{\mathrm{J}}\times\bm{\mathrm{B}}$. As in Sec. II, we stick to two dimensions for simplicity. The current is given by $\mu_{0}\bm{\mathrm{J}}=\bm{\mathrm{\nabla}}\times\bm{\mathrm{B}}$. For convenience, we use somewhat different notation than before; here we refer to quantities at $t=0$ as $\bm{\mathrm{J}}_{0}$ and $\bm{\mathrm{B}}_{0}$ and refer to their first-order perturbations at $t=\Delta t$ as $\bm{\mathrm{J}}_{1}$ and $\bm{\mathrm{B}}_{1}$. From the magnetic field in Eq. (II), we can compute the current,

[TABLE]

The Lorentz force per unit volume $\bm{\mathrm{F}}=\bm{\mathrm{J}}\times\bm{\mathrm{B}}$ has an $x$ component

[TABLE]

Thus, we have

[TABLE]

We assume that $\overline{F}_{0x}=0$, which is satisfied if $\overline{B_{0x}B_{0y}}=0$.

Let us now consider $F_{1x}$, and substitute the previous expressions to write everything in terms of $U$ and $\bm{\mathrm{B}}_{0}$. From now on, we drop the 0 subscript for simplicity. After some algebra, we obtain for the zonally averaged force,

[TABLE]

This expression is equivalent to Eq. (10).

To summarize, we have recovered the same net force using $\bm{\mathrm{J}}\times\bm{\mathrm{B}}$ directly as we did earlier with the Maxwell stress tensor. In the 2D MHD system where $\bm{\mathrm{B}}$ is in the $xy$ plane, both $\bm{\mathrm{J}}_{0}$ and $\bm{\mathrm{J}}_{1}$ are along ${\bm{\mathrm{\hat{z}}}}$, implying currents must be able to flow in the $z$ direction.

Appendix B Derivation of negative eddy viscosity

In this appendix we derive for 2D hydrodynamic flow the negative eddy viscosity discussed in Sec. III using the same technique as we used to derive the magnetic eddy viscosity.

The zonally averaged force per unit volume in the $x$ direction is given by the Reynolds stress

[TABLE]

We use the vorticity formulation of the equation of motion. We suppose at $t=0$ there is a blob of vorticity $\zeta_{0}(x,y)=z_{0}\cos(\bm{\mathrm{k}}\cdot\bm{\mathrm{x}})$. One could multiply the spatial dependence by an envelope function to localize the spatial extent, but that is unnecessary for the effect of interest. Consider the evolution of this vorticity blob due to advection by an imposed shear $\bm{\mathrm{v}}=U(y){\bm{\mathrm{\hat{x}}}}=U_{0}\sin(qy){\bm{\mathrm{\hat{x}}}}$. We neglect both the Coriolis and the Lorentz force in this derivation. We have

[TABLE]

To include the high-wave-number cutoff discussed in Sec. III, one should retain the $\bm{\mathrm{v}}^{\prime}\cdot\bm{\mathrm{\nabla}}\overline{\zeta}$ term, but we neglect it here for simplicity.

We evolve $\zeta$ for a short time increment $\Delta t$. The vorticity becomes

[TABLE]

where $\bm{\mathrm{k}}_{\pm}=\bm{\mathrm{k}}\pm\bm{\mathrm{q}}$ and $\bm{\mathrm{q}}=q{\bm{\mathrm{\hat{y}}}}$. From $\zeta=\nabla^{2}\psi$, we obtain

[TABLE]

Calculating $\bm{\mathrm{v}}(\Delta t)={\bm{\mathrm{\hat{z}}}}\times\bm{\mathrm{\nabla}}\psi(\Delta t)$ and then constructing $\overline{v_{x}(\Delta t)v_{y}(\Delta t)}$, we obtain

[TABLE]

We identify $qU_{0}\cos(qy)=\partial_{y}U$, as well as $k_{x}^{2}z_{0}^{2}/k^{4}=k_{x}^{2}\psi_{0}^{2}\approx v_{y}^{2}$ and $\overline{v_{y}^{2}}=v_{y}^{2}/2$. Furthermore, we observe that the factor $1-4k_{y}^{2}/k^{2}$ is only positive for $k_{x}^{2}>3k_{y}^{2}$. This $\bm{\mathrm{k}}$-dependent prefactor does not arise in the magnetic calculation. We collapse this factor into a $\bm{\mathrm{k}}$-independent positive constant $\gamma$, recognizing that $\gamma$ may be small in magnitude because $1-4k_{y}^{2}/k^{2}$ is only positive in a small fraction of $\bm{\mathrm{k}}$ space. Finally, as in Sec. II, we replace $\Delta t$ with $\tau_{\rm corr}$, the maximum time the shear can act coherently on a structure before being decorrelated. Again, we interpret our results in a statistically averaged sense. We obtain a Reynolds stress

[TABLE]

This expression corresponds to a negative dynamic eddy viscosity

[TABLE]

This calculation does not directly generalize to 3D. Negative eddy viscosities have in the literature been limited to 2D or quasi-2D flow, and the quasi-2D nature of flow is a consequence of rotation. In contrast, the positive magnetic eddy viscosity is derived in 3D in Section VI without such complications.

We remark on the relation between the wave-number dependence here and that found in a related Kelvin–Orr calculation. In that calculation with an initial wave shearing in a velocity profile $U=Sy{\bm{\mathrm{\hat{x}}}}$, the energy transfer to the mean flow is found to have a factor $1-4k_{y}^{2}/k^{2}$, as in Eq. (33) Bakas and Ioannou (2013); Bakas et al. (2015); Constantinou and Parker (2018). This wave-number dependence arose when a $\bm{\mathrm{k}}$-independent friction was used as the decorrelation mechanism. When a viscosity instead acts as the decorrelation mechanism, the length-scale dependence of the viscosity modifies the factor in the energy transfer from $1-4k_{y}^{2}/k^{2}$ to $1-6k_{y}^{2}/k^{2}$. Returning to the calculation here, we observe that we did not include any explicit friction or viscosity. Instead, we posited that decorrelation occurred due to turbulent inertial motion, and we implicitly assumed it was independent of wave number. Hence, we recovered the same factor $1-4k_{y}^{2}/k^{2}$ as found using a $\bm{\mathrm{k}}$-independent friction.