Effective drag in rotating, poorly conducting plasma turbulence

Despite the increasing sophistication of numerical models of hot Jupiter atmospheres, the large time-scale separation required in simulating the wide range in electrical conductivity between the dayside and nightside has made it difficult to run fully consistent magnetohydrodynamic (MHD) models. This has led to many studies that resort to drag parametrizations of MHD. In this study, we revisit the question of the Lorentz force as an effective drag by running a series of direct numerical simulations of a weakly rotating, poorly conducting flow in the presence of a misaligned, strong background magnetic field. We find that the drag parametrization fails once the time-scale associated with the Lorentz force becomes shorter than the dynamical time-scale in the system, beyond which the effective drag coefficient remains roughly constant, despite orders-of-magnitude variation in the Lorentz (magnetic) time-scale. We offer an improvement to the drag parametrization by considering the relevant asymptotic limit of low conductivity and strong background magnetic field, known as the quasi-static MHD approximation of the Lorentz force. This approximation removes the fast time-scale associated with magnetic diffusion, but retains a more complex version of the Lorentz force, which could be utilized in future numerical models of hot Jupiter atmospheric circulation.


INTRODUCTION
Hot Jupiters (HJs) are gas giant exoplanets, with masses similar to that of Jupiter, which orbit close enough to their host star that they are generally expected to be tidally locked (Seager 2010;Heng & Showman 2015). The proximity to their host stars is also expected to partially ionize the upper atmospheres of HJs, leading to the interaction between the atmospheric flows and any present magnetic fields (Batygin & Stevenson 2010;Perna et al. 2010a,b;Koskinen et al. 2010;Menou 2012;Koskinen et al. 2014). Their relatively large profile when obscuring their host star, as well as their short orbital periods, make them ideal candidates for transiting observations. In the last two decades, these obser-vations have given us access to a great amount of information about their atmospheres, including some insight into their atmospheric dynamics -indirectly via hot spot migrations (Knutson et al. 2007(Knutson et al. , 2008Zellem et al. 2014;Bell et al. 2021) and more directly using the blue-shifting of spectra observed at the terminators (Louden & Wheatley 2015;Ehrenreich et al. 2020). Their discovery, and the subsequent observation of atmospheric dynamics, prompted the creation of a whole sub-field devoted to the numerical modeling of the atmospheres of HJs and their close relatives. These models range in sophistication and intent, from quasi-two-dimensional shallow water models (Cho et al. 2003;Langton & Laughlin 2007;Cho 2008;Showman & Polvani 2011;Perez-Becker & Showman 2013;Heng & Workman 2014;Hindle et al. 2019) to three-dimensional general circulation models (GCMs) (Showman & Guillot 2002;Dobbs-Dixon & Lin 2008;Showman et al. 2009;Menou & Rauscher 2009;Perna et al. 2010a;Rauscher & Menou 2010;Heng et al. 2011;Rauscher & Menou 2013;Batygin et al. 2013;Rogers & Komacek 2014;Rogers & Showman 2014;Rogers & McElwaine 2017).
One of the largest obstacles in modeling these atmospheres is the large conductivity contrast between the dayside and the nightside, due to the large differences in temperature (Perna et al. 2010a;Rogers & Komacek 2014;Heng & Showman 2015). The time scale of diffusion for induced magnetic fields is proportional to this conductivity, resulting in numerical models needing to resolve small time scales in the nightside, along with large time scales to capture large-scale structures in the atmosphere. Very large time-scale separations can be impractical for numerical simulations, and, as a result, many modelers resort to a parametrization of magnetic effects that doesn't directly resolve the magnetic diffusion time-scale. The most common approach begins by assuming that any induced magnetic field is a small, rapidly diffusing perturbation around a strong background magnetic field. This leads to a time-scale associated with the Lorentz force that is equal to ρσ −1 e B −2 0 , where σ e is the electrical conductivity, B 0 is the strength of the background magnetic field, and ρ is the fluid density (Davidson 1995(Davidson , 2013Knaepen & Moreau 2008). Modelers then substitute the Lorentz force with a drag term with an associated time-scale τ drag = ρσ −1 e B −2 0 (Perna et al. 2010a;Menou 2012;Rauscher & Menou 2012;Komacek & Showman 2016;Koll & Komacek 2018;Kreidberg et al. 2018;Arcangeli et al. 2019), which can also vary in space (Rauscher & Menou 2013;Beltz et al. 2021). For small values of σ e the resulting drag time-scale is not restrictive for the numerics. Studies that have implemented what is dubbed as 'MHD drag' have found that the structure of the atmospheric circulation significantly changes when the drag time-scale is similar to or smaller than the relevant dynamical timescale. There are, however, reasons to question the validity of MHD drag. Recent attempts at modeling the full MHD equations in 3D GCMs, albeit with some simplifications, have shown that, although magnetic effects do reduce the strength of atmospheric jets (as a drag would), they also cause different morphological changes in the flow which make the authors question whether the correct prescription is a drag (Rogers & Komacek 2014;Rogers 2017). Furthermore, the authors find that the Ohmic dissipation measured in their MHD models is more than an order of magnitude smaller than would be predicted by a drag term. In an independent study, Heng & Workman (2014) arrive at a similar conclusion about MHD drag by studying the shallow water MHD model. See also Pothérat & Klein (2017) for a similar discussion in the context of MHD experiments. Indeed, as we'll see in the next section, when considering the approximation of the Lorentz force in the relevant limits of low conductivity and large background magnetic field, one can show that a drag-like term appears in twodimensional flows, but in three-dimensional flows this is not the case (Davidson 1995(Davidson , 2013. Whether this approximation can be further reduced to a drag is unclear. The uncertainty of both the drag time-scale as well as the validity of the drag prescription itself, combined with their significant effects on the atmospheric circulation, makes this a central issue in the modeling of HJ atmospheres and impacts our understanding and expectation of their atmospheric circulation. In this study we introduce a simplification of MHD in the low conductivity limit and revisit the question of the Lorentz force as an effective drag. In section 2 we consider an approximation to the full MHD equations in the relevant limits of low conductivity and strong background magnetic field, called quasi-static MHD (QMHD). In particular, we focus on the form of the approximate Lorentz force and discuss its properties and potential relation to a drag-like term, whose validity we quantify using an effective drag coefficient. In section 3 we describe the numerical setup and introduce an integral length-scale along the magnetic field, as a measure of the anisotropy in the flow. Then, in section 4 we measure the anisotropy and effective drag coefficient in a series of direct numerical simulations of QMHD turbulence in an idealized setup. We find that the drag parametrization with associated time scale τ drag = ρσ −1 e B −2 0 works well for runs in which the dynamical time-scale is shorter than that which is associated with the Lorentz force. Beyond this, when the Lorentz force is sufficiently strong, the flow becomes anisotropic and the effective drag coefficient levels off. Finally, in section 5 we summarize our results and propose QMHD as an intermediate model, bridging the gap between the simplicity of a drag and the complexity of the full MHD equations, to be used in future GCMs of HJs.

ROTATING, WEAKLY CONDUCTING MHD
We aim to keep our fluid description of HJ atmospheres as simple as possible in an effort to focus purely on the dynamical effects of the Lorentz force. This means we will be ignoring many realistic features of HJ atmospheres, including stratification, radiation, daynight forcing contrast due to tidal locking, compressibility, and kinetic plasma effects. Given their moderate temperatures, HJ atmospheres are also likely partially ionized (Batygin & Stevenson 2010;Perna et al. 2010a;Koskinen et al. 2014). However, except at very low pressures, the ions and neutrals are expected to be highly coupled due to collisions, meaning that a single-fluid description is an appropriate characterization (Perna et al. 2010a;Benavides & Flierl 2020). We thus begin by considering incompressible magnetohydrodynamics (MHD) with uniform density (and no buoyancy variations), subject to rotation Ω and a uniform, steady background magnetic field B 0 . We are effectively considering a three-dimensional volume of a HJ atmosphere that is smaller than both the stratification scale height set by the entropy gradient and also the length scale of variation in the background magnetic field. Further simplifications will be achieved by considering two relevant limits, low electrical conductivity, and a strong background magnetic field.
The effect of electrical conductivity on the dynamics is quantified by the magnetic Reynolds number Re m , a dimensionless parameter comparing the magnetic diffusion time-scale to the time-scale associated with the evolving flow (Davidson 2013). We define Re m = u/η, where η = (µ 0 σ e ) −1 is the magnetic diffusivity, µ 0 is the magnetic permeability, is a dominant length-scale of the flow, and u is a velocity scale. If Re m 1, the diffusion and dissipation of induced magnetic fields is significant. Hot Jupiters with daysides cooler than roughly 1800 K are expected to have magnetic Reynolds number smaller than one throughout most of their atmospheres (Perna et al. 2010a;Hindle et al. 2021), although this assumption could break down on the dayside of the hotter HJs at lower pressures (Menou 2012;Rogers & Komacek 2014). Importantly, dynamo instabilities are not present in flows with Re m 1, and thus do not convert kinetic energy to magnetic energy, resulting in a decaying induced magnetic field b and a negligible Lorentz force, j × b, where j = µ −1 0 ∇ × b, in the absence of a magnetic field generated elsewhere.
However, in the presence of a background magnetic field (from a deeper dynamo region or from the host star), the flow can act to exchange kinetic for magnetic energy by shearing this field, inducing currents and magnetic fields. In a flow with low conductivity, the strength of this induced magnetic field scales like b ∼ Re m B 0 , and thus the Lorentz force scales like (Davidson 1995(Davidson , 2013Knaepen & Moreau 2008). This is the origin of τ drag discussed in section 1. We estimate the relevance of the Lorentz force in the dynamics by comparing the strength of the Lorentz force to the nonlinear advection term in the momentum equation, giving us our main control parameter in this study, known as the interac-tion parameter : is large enough such that N 1, then the Lorentz force can significantly affect the dynamics. Some recent studies have estimated magnetic field strengths for HJs and found magnitudes similar to that of Jupiter, but possibly up to 50 times greater for larger radius HJs, suggesting that this limit could be relevant for some HJs (Reiners & Christensen 2010;Yadav & Thorngren 2017;Rogers 2017;Hindle et al. 2021).
Flows with Re m 1 and N ∼ O(1) have the distinct property that the induced magnetic field is quickly diffused away, yet the Lorentz force is not negligible. This limit is referred to as the quasi-static approximation to MHD (which we call 'QMHD' henceforth) (Moffatt 1967; Sommeria & Moreau 1982;Davidson 1995Davidson , 2013Knaepen & Moreau 2008), and has been studied mainly in metallurgy and in MHD experiments due to the typically low conductivity of liquid metals (Alemany et al. 1979;Sommeria 1988;Gallet et al. 2009;Klein & Pothérat 2010;Pothérat & Klein 2014;Baker et al. 2018), although recent numerical studies on its turbulent properties and anisotropy have been done as well (Zikanov & Thess 1998;Burattini et al. 2008;Favier et al. 2010Favier et al. , 2011Reddy & Verma 2014;Verma 2017). After nondimensionalizing the equations of MHD using the uniform density ρ, and u, and taking the limits above, one is left with a single dynamical equation for the velocity: where p * is the total pressure modified by rotation and magnetic pressure, Ro −1 ≡ 2Ω /u is the inverse Rossby number (quantifying the relative strength of the Coriolis force),x Ω andx B0 are unit vectors in the direction of rotation and the background magnetic field, respectively, and F is a generic forcing term that can include dissipation such as viscosity and a body force (to be specified in section 3). The background magnetic field is fixed in time and is uniform in space, such that ∇×B 0 = B 0 ∇×x B0 = 0. Care must be taken if considering a spatially-dependent background magnetic field B 0 (x), as the resulting equation will not be the same. See the discussion in section 5. This equation is accompanied with the incompressibility condition ∇ · v = 0. The induced magnetic field can be found using a diag- In QMHD, the Lorentz force operator, acts to dissipate kinetic energy from any motion which varies along the direction of the background magnetic field. In two-dimensional flows, with an in-plane B 0 , it can be shown that ⊥ is the projection of the velocity perpendicular to the background field andṽ ≡ v − L −1 v dx B0 (Davidson 1995(Davidson , 2013. These two properties would, at first glance, seem to justify the use of a drag parametrization. However, most HJs are expected to be tidally locked, resulting in order one Rossby numbers (Showman et al. 2010), which is not low enough for strong two-dimensionalization of these flows (and we are not considering the effects of strong stratification or large aspect ratio). In three dimensions, there are various reasons to question the use of a drag term as an approximation to the Lorentz force in QMHD. First of all, the Lorentz force L(v) acts to create anisotropy in the flow by removing energy from motions that vary along the magnetic field, but does not affect motions that are invariant along that direction. This is in sharp contrast to a drag in which any motion is affected equally and in an isotropic way. Second, a drag term removes linear momentum from the system, whereas the Lorentz force in the QMHD approximation does not (in the absence of no-slip boundary conditions). This can have profound effects on the time evolution of the system. Finally, a drag term removes energy anywhere in space where v = 0, whereas L(v), despite removing energy in a spatially-averaged sense, may locally inject energy into the flow (Figure 1). Despite these differences, L(v) does remove energy from the flow, and is linear in v, so there is hope that the much simpler prescription of L as a drag could be approximately valid in certain regimes. One can arrive at a potential justification for this by expressing L(v) in Fourier space in terms of wavenumbers parallel (k ) and perpendicular (k ⊥ ) to the background magnetic field, resulting in , wherev is the Fourier transform of v. We see that, if max k |v|(k , k ⊥ ) happens for k ≥ |k ⊥ |, then we can approximate L(v) ≈ L drag (v) = −c 0 N v, with c 0 being an order one constant ranging from c 0 ≈ 1/2, when the max occurs near k ∼ |k ⊥ |, to c 0 ≈ 1, when the max occurs for k |k ⊥ |. The latter might be relevant in a thin atmosphere where the background magnetic field projects significantly onto the thin direction and the flow contains large horizontal structures perpendicular to the background magnetic field.
In order to investigate the validity of the Lorentz force as a drag in a quantitative way, we introduce an effective drag coefficient D eff , where · denotes a temporal and spatial average at steady state. If L(v) does indeed act like a drag, at least in a volume-averaged sense, then we would expect D eff ∝ N . However, if this is not the case, then D eff will deviate from N . v ·L(v) is the Ohmic dissipation rate, so that D eff represents the ratio of Ohmic dissipation to (twice) the kinetic energy in the flow.
In the next section we introduce a series of direct numerical simulations of QMHD turbulence which we use to investigate how L acts to create anisotropy, and how that, in turn, affects the validity of L approximated as a drag. In particular, we use these simulations to determine how the effective drag coefficient depends on the interaction parameter N .

METHODS
We performed direct numerical simulations of the QMHD system, equation (2), in a triply-periodic domain using a modified version of the Geophysical High-Order Suite for Turbulence (GHOST) (Mininni et al. 2011;Benavides 2021), a pseudo-spectral code with a fourthorder Runge-Kutta scheme for time integration and a two-thirds dealiasing rule. The generic forcing term F comprised of a 'hyper'-viscous term, −ν∇ 4 v, which acts to dissipate kinetic energy at the smallest scales, and a body forcing term f , which is random (white-in-time) and injects energy into the flow at a constant rate ε and at a single length-scale f , both of which are input parameters (Chan et al. 2012). The hyper-viscosity lets us use lower resolutions while maintaining numerical accuracy, and has been shown to have no significant effect on the turbulent properties of 3D turbulence as long as the power of the gradient is small enough, as is the case here (Agrawal et al. 2020).
The axis of rotation was chosen to be in the directionx Ω =ẑ, whereas the direction of the background magnetic field was chosen to be perpendicular to it, in the x-direction,x B0 =x. The misalignment between the rotation axis and background field is supposed to reflect a generic case, since there is no reason to expect alignment between the two outside the dynamo Figure 1. Snapshots from the Lx = 8πL runs (Table 1) of the field-aligned vorticity ω =x · (∇ × v) (lower left, in blue and red) and the Ohmic dissipation v · L(v) (upper right, in black and red) for increasing values of the interaction parameter N (equation (1)). Figures 1(b)-(d ) represent runs with approximately equal values of volume-averaged Ohmic dissipation and D eff (equation (5)). The red colors represent positive values whereas the blue and black represent negative values. All snapshots use the same colorbar scale for each given field. region (e.g., a dipole field from the interior dynamo region). In some cases misaligned rotation and background field can have significant consequences on the dynamics (Benavides et al. 2022), so we performed runs at other misalignment angles, θ, defined to be the angle the background field makes with the rotation axis (Table 1). These runs show that our results do not depend strongly on the misalignment angle, as long as it's not zero or small.
Since L dissipates motion that varies along x, we expect anisotropy to develop in our domain, manifested by structures along the x-direction which are larger than in the perpendicular directions. In order to accommodate the form of L, and the resulting anisotropy, we make a few specific choices in our implementation. First of all, we perform the same set of runs for various val-ues of L x , the domain size in the x-direction. The domain size in the perpendicular direction is fixed at L y = L z = 2πL, whereas we perform sets of runs with L x = 2πL, 4πL and 8πL. Second, the forcing function f 'stirs' the fluid in a manner that does not vary along the x-direction, so that the forcing does not project onto L, which would immediately dissipate the energy being forced. Since such two-dimensional forcing results in three-dimensional instabilities, the resulting flow is still approximately isotropic when N 1. The numerical model is nondimensionalized by L and ε, such that the domain size in the directions perpendicular to B 0 is 2π and the forcing function has an injection rate equal to 1. In all of our runs, we have ν = 2 × 10 −6 , Ω = 2, and f = 2π/k f , where k f = 9 (we randomly force modes k such that 8 < |k| < 10, making  Table 1. A summary of the runs performed for this work. All runs have (in simulations units, nondimensionalized by L and ε) ν = 2 × 10 −6 , Ω = 2, and f = 2π/k f , where k f = 9. This corresponds to an inverse Rossby number of Ro −1 = 1.98, and an 'effective' Reynolds number of Re = 14973. Nx represents the number of grid points in the x-direction. All runs have 256 grid points in both y and z directions. k f = 9). The dominant length-scale in the problem is = f /2, if we consider the size of a typical vortex produced by the forcing. The velocity scale is defined using u = (ε ) 1/3 = 1/3 = (π/9) 1/3 . These values of and u result in an inverse Rossby number of Ro −1 = 1.98, and an 'effective' Reynolds number of Re = 14973, where Re ≡ u 3 /ν is based on the hyper-dissipation used in our model. Although we do not know typical values for viscosity on hot Jupiters, given the high temperatures present, we expect it to be very small. Furthermore, with scales of motion on the order of thousands of kilometers and velocities on the order of meters per second or more, we believe that the Reynolds numbers will be much larger than one. We thus expect small scale instabilities (e.g. shear instabilities) as well as buoyancybased instabilities to be present, leading to turbulence. All of our runs have a resolution of 256 in each direction perpendicular to the background field. See Table  1 for a list of the runs performed in our study. Note that, although we are varying N over many orders of magnitude (including possibly N 1), QMHD is the limit of MHD with vanishing magnetic Reynolds number, which always holds true. Therefore, N 1 requires B 0 /( √ µ 0 ρu) to be much larger than one.
Each simulation is run until a statistical steady-state is reached, at which point the time-averages are taken until the error estimate (considering covariance) levels off. For each run we calculate the effective drag coefficient D eff , as well as an integral length-scale in the x-direction, defined as: where E(k x ) is the time-averaged, one-dimensional energy spectrum in the x-direction. Note that these integrals also include contributions from the k x = 0 mode.
x gives an estimate of the dominant length-scale in the x-direction (in units of L), and will be used as a quanti-,QWHUDFWLRQSDUDPHWHUN ,QWHJUDOOHQJWKVFDOH x L x = 2πL ℓ = • ℓ = • L x = 4πL L x = 8πL N 1/2 1.33 × f 2π 4π 8π Figure 2. The integral length-scale x (equation (6)), in units of L, as a function of the interaction parameter N (equation (1)), for various box sizes and misalignment angles. In gray we show 1.33 times the forcing length-scale range. Anisotropy develops once N > 1, and this does not depend on the domain size or misalignment angle. For N > 1 and sufficiently large domains which do not suffer from finite size effects (when x < Lx, see axis on right), we find that x ∼ N 1/2 (red, dashed line).
tative measure of anisotropy developing in the domain. In an isotropic system we would expect x ∼ f .

Anisotropy
Our simulations show that significant anisotropy develops in the flow once N 1 (Figures 1 and 2). For N < 1 the Lorentz force is negligible and does not affect the dynamics, resulting in approximately isotropic flow with x ≈ 1.33 f . The larger-than-one prefactor likely comes from the fact that the x-invariant structures formed by the forcing are unstable to three-dimensional perturbations at many length-scales. Another thing to note is that the onset of anisotropy does not depend on domain size or misalignment angle.
As N increases beyond one, the Lorentz force dissipates structures that vary along the x direction, acting more strongly on those with small length-scales. This results in more elongated structures as N increases (Figure 1). x grows until all k x (= k ) > 0 modes are stable and the flow becomes exactly two-dimensional at a critical value N 2D , where N 2D ∝ L 2 x (Zikanov & Thess 1998;Thess & Zikanov 2007;Favier et al. 2010;Gallet & Doering 2015). The effects of N 2D are seen for N < N 2Dfinite domain size effects appear when x L x ( Figure  2). The exact two-dimensionalization will also depend on the Reynolds number, with N 2D increasing with Re (Gallet & Doering 2015). Given the large physical extent and Reynolds numbers of astrophysical flows, we don't expect the two-dimensionalization to be a relevant physical phenomenon. Therefore, to avoid these effects which we believe to be irrelevant for our motivation, we consider larger domains, which allow us to push N 2D to larger values, and therefore begin to approach the astrophysically-relevant regimes.
The runs on larger domains reveal a power-law dependence of x with N , with an observed anisotropy scaling of x ∼ N 1/2 (Figure 2). This agrees with previous scaling predictions, such as that by Sommeria & Moreau (1982) who considered L as an along-field diffusion L(v) ≈ L diff (v) = κ∂ 2 x v, with κ ∼ σB 2 0 2 ⊥ /ρ and ⊥ ∼ f based on our forcing. In an eddy turnover time τ eddy , motions with horizontal extent ⊥ would diffuse vertically with a diffusion length of x ∼ √ κτ eddy ∝ B 0 ⊥ ∝ N 1/2 ⊥ . This scaling can also be arrived at by considering the wavenumbers whose turbulent timescales match that of the Lorentz force, and are therefore damped away. This gives k u(k) ∼ k 2/3 ∼ N k 2 x /k 2 , resulting in k −1 x k 4/3 The slightly different scaling for ⊥ comes from the dependence of τ eddy on ⊥ , based on 3D homogeneous and isotropic turbulence assumptions (Frisch 1995), which is not considered in the diffusivity argument.
For 1 N < N 2D the anisotropy is such that the flow is almost two-dimensional (e.g., Figure 1 (d )). Previous studies looking at turbulent energy cascades have found that inverse energy cascades (associated with twodimensional hydrodynamics) appear before exact twodimensionalization (Smith et al. 1996;Celani et al. 2010;Alexakis 2011;Deusebio et al. 2014;Sozza et al. 2015;Benavides & Alexakis 2017;Alexakis & Biferale 2018;Pouquet et al. 2019). However, in our runs we don't see any scale coarsening in the directions perpendicular to the background field. This is due to the rotation in the z-direction, coupling the horizontal and out-of-plane velocities which results in a system with a forward cascade of energy (Benavides et al. 2022). If rotation were to be weaker, we would expect the formation of an inverse cascade. Indeed, this seems to be occurring for the θ = 30 • run, where the projection of the rotation perpendicular to the background field is smaller, resulting in weaker in-plane rotation rate. The inverse cascade for this case results in larger horizontal scales ⊥ , which we believe pushes N 2D to lower values (Figure 2, light blue). On the other hand, for cases with very fast rotation, the Taylor-Proudman theorem would manifest itself as flow becoming invariant along the z-direction, resulting in a series of shear layers varying in the third direction, y, and a suppressed energy cascade (Benavides et al. 2022). This latter case might be more relevant for the transition regions of gas giant planets like Jupiter and Saturn, where the conditions for QMHD are also likely satisfied, but where rotation rates are significantly larger than those of HJs.

Effective drag
An increase in x implies a decrease in the x-derivative found in L, thereby effectively lowering the Ohmic dissipation. However, the decrease in the x-derivative occurs as we increase N , which also appears in L. What is the combined effect on D eff of N increasing but the xderivative decreasing? Figure 3 shows the effective drag coefficient D eff as we vary the control parameter N .
For N < 1, while the flow is approximately isotropic, we see a good agreement with D eff ∝ N , suggesting that a drag-like parametrization could correctly capture the dynamics and Ohmic dissipation in this regime, at least in a volume-averaged sense. Indeed, D eff ≈ N/2, which seems to validate the arguments made in the end of section 2 for a max energy occurring when k ∼ |k ⊥ |, as expected for the 3D instability of the forced 2D structures. However, as anisotropy develops for N 1, the structures that dissipate the most energy appear at larger scales ( Figure 1) and D eff begins to deviate from the one-to-one line. Much like the anisotropy, the deviation from the one-to-one line near N ∼ 1 does not depend on the domain size or misalignment angle.
It is not clear a priori what the behavior of D eff should be beyond this point. For the L x = 2πL runs, D eff begins to decrease beyond N ∼ 1. However, this is a result of the finite domain size and proximity to N 2D . By looking at successively larger L x runs, we probe what would happen in a more realistic setting. Figure 3 suggests that, for large L x , the effective drag coefficient D eff levels off and remains roughly constant, despite orders of magnitude increase in N . We found that this behavior and value of D eff does not depend strongly on Re. Figure 1 (b)-(d ) shows snapshots of runs which have approximately the same effective drag coefficient, while representing three orders of magnitude for N . Structures change in such a way so as to keep D eff roughly constant, given the increase in N .
Given our findings from Figure 2, we can see why this behavior is a result of the anisotropy scaling x ∼ N 1/2 . Combining equations (4) and (5), we can re-frame D eff as a velocity-weighted average of wavenumbers: where k 2 = k 2 x + k 2 y + k 2 z . When N > 1, we would expect k 2 x < k 2 ⊥ = k 2 y + k 2 z , so that we can replace k 2 with k 2 ⊥ in equation (7). Assuming k ⊥ doesn't vary significantly, since the forcing remains the same and no large-scale ,QWHUDFWLRQSDUDPHWHUN Figure 3. The effective drag coefficient D eff (equation (5)) versus interaction parameter N (equation (1)), for various box sizes and misalignment angles. For N < 1, the effective drag coefficient is proportional to N (black, dashed line), and seems to follow N/2 (black, dot-dashed line), suggesting that a drag formulation could be valid, with c0 ≈ 1/2, as expected for an energy maximum near k ≈ |k ⊥ |. However, the curve levels off and deviates significantly from the drag prediction by orders of magnitude when N > 1, independent of domain size and misalignment angle θ. For N 1, the effective drag drops as the flow becomes exactly two-dimensionalized at N2D; this limit does depend on the domain size (see section 4).
structures form in the flow, we can approximate it with k ⊥ ∼ 2π −1 f . Finally, substituting k 2 x = (2π x ) −2 , we end up with: In other words, the effective drag coefficient is N divided by the anisotropy of the flow squared. Since we know how the anisotropy scales with N , we can get D eff as a function of N , which in the anisotropic limit of 1 N N 2D turns out to be a constant. Recall that D eff is proportional to the Ohmic dissipation rate, so that a levelling-off of D eff also represents the same behavior for the Ohmic dissipation 1 . A similar plateauing behavior has been observed for the Ohmic dissipation in HJ GCM runs with increasing temperatures (Rogers & Komacek 2014). The authors attribute this to a change in dynamics as Re m increases past one at lower pressures. However, they note that at higher pressures, and for some of the lower temperature runs 1 The kinetic energy is approximately the same for all runs here, partly due to the presence of a forward cascade. This might change in the case of weaker rotation and subsequent formation of an inverse cascade.  Figure 1, as well as some runs with lower N . Values of a < 0 act to locally dissipate energy, whereas a > 0 act to add energy. A drag term would result in a = −1 at all points in space (blue, dot-dashed line). Although the effective drag coefficient would have us believe that the arguments used to justify L as a drag are valid, this figure shows that there are still aspects of this approximation which do not hold even for N < 1.
where the plateau begins, Re m is still low, so it's possible that the QMHD effects seen here might be present and partially responsible for what is observed. Although Figure 3 seems to suggest that a drag prescription for L is valid for N < 1, we emphasize that the effective drag coefficient is only a volume-averaged measure of this validity. As mentioned in section 2, unlike a true drag term, the Lorentz force operator only dissipates the volume-averaged energy and can locally cause energy increase in certain regions, similar to a diffusion term. This fact is visible in Figure 1 and is quantified in Figure 4, which shows, for a single snapshot in time, the probability density function (PDF) in space of the alignment of velocity and the Lorentz force in terms of the dot product a ≡ v · L(v)/(|v||L(v)|) for the four snapshots seen in Figure 1, as well as for some runs with lower N . A true drag term would result in a delta function distribution around a = −1 (blue, dot-dashed line). However, Figure 4 shows values of a between −1 and 1 even for the runs with N < 1, suggesting both positive and negative contributions to the local energy balance due to the Lorentz force, while clearly showing overall energy dissipation given by the peaking PDFs for a = −1. All of this suggests that, although a drag might capture the energy balance correctly (albeit, only for N < 1), this approximation will not result in the same spatial structures -a drag is not the same as an along-field diffusion.
That said, large-scale HJ atmosphere simulations might not resolve along-field scales very well. If one considers under-resolving as an effective averaging of such scales, then our results could justify the implementation of an effective drag.

CONCLUSIONS
In this study, we considered the turbulent dynamics of rotating MHD in the presence of a background magnetic field, in the combined limit of Re m 1 and B 0 /(u √ µ 0 ρ) 1, termed quasi-static MHD (QMHD). Motivated by approaches used in the study of hot Jupiter (HJ) atmospheres, we have shown that a drag parametrization of the Lorentz force operator L fails once the ratio of the dynamical timescale to the Lorentz timescale, quantified by the interaction parameter N , is larger than one. This happens because the Lorentz force dissipates structures that vary along the background field, creating anisotropy in the flow, which in turn acts to reduce the Ohmic dissipation, thereby reducing the effective drag. The development of anisotropy with increasing N is such that the effective drag coefficient remains constant for N > 1, despite N varying by orders of magnitude. The levelling off of D eff for N > 1 has significant implications for simulations parametrizing L as a drag, since we see values of D eff deviating by orders of magnitude from what would be predicted if one assumes L(v) = L drag (v) = −c 0 N v. This could also result in severely overestimating the amount of Ohmic dissipation, as well as misrepresenting the true dynamics of HJ atmospheres. Although we have shown more generally that the drag prescription does not fully represent the form of the Lorentz force at low conductivities, we are aware of the various difficulties in performing a GCM run with full MHD, motivating the use of MHD drag. It would be of interest to see if our results on effective drag carry through to full MHD or QMHD GCMs, and if the implementation of an MHD drag of the form D eff (x) = min{N (x), N 0 }, with N 0 being a constant of order one, is valid.
The main motivation for a drag parametrization of the Lorentz force comes from the restrictively small time scales associated with very large magnetic diffusivities (low electrical conductivity). The drag time-scale is much larger than the time-scale associated with magnetic diffusion, allowing modellers to bypass this problem. While we have found that a drag parametrization fails for N > 1, we want to emphasize that the same time-scale advantage exists in the QMHD limit, despite the more complicated form of the operator associated with the Lorentz force. This will hopefully motivate the use of the QMHD approximation in models of HJ atmospheres. Even for the case of a spatially-dependent background magnetic field or conductivity, its implementation would be straight forward if one considers separately the Lorentz force µ −1 0 (∇ × b) × B 0 (x) and the induced magnetic field b = −∇ −2 (η(x) −1 ∇×(v×B 0 (x))). An alternative which might be easier to implement would be to approximate −∇ −2 with some horizontal length-scale 2 ⊥ , similar to what was done when considering L(v) as an along-field diffusivity, L(v) ≈ L diff (v) = κ∂ 2 x v, with κ ∼ σB 2 0 2 ⊥ /ρ (Sommeria & Moreau 1982). Although the expression for L (equation (4)) would be modified in the presence of a spatially-dependent background magnetic field, we expect our results to hold for those cases, as well.
We thank Thaddeus D.
Komacek for insightful discussions and helpful suggestions. This research was carried out in part during the 2019 Summer School at the Center for Computational Astrophysics, Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. SJB acknowledges funding from the National Aeronautics and Space Administration (Award Number: 80NSSC20K1367) issued through the Future Investigators in NASA Earth and Space Science and Technology (NNH19ZDA001N-FINESST) within the NASA Research Announcement (NRA): Research Opportunities in Space and Earth Sciences (ROSES-2019).