2. Rotating MHD in the presence of a background magnetic field

The equations for rotating MHD in the presence of a uniform background magnetic field are (Shebalin Reference Shebalin2006; Galtier Reference Galtier2014)

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{v} ={-} \boldsymbol{\nabla} p^* - 2 \boldsymbol{\varOmega}\times \boldsymbol{v} + ( \boldsymbol{\nabla} \times \boldsymbol{b} ) \times (\boldsymbol{B_0}+\boldsymbol{b}) + \boldsymbol{D}_v + \boldsymbol{f}, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{b}}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{b} = (\boldsymbol{B_0} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{v} + (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{v} + \boldsymbol{D}_b, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0, \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{b} = 0, \end{gather}$$

where $\boldsymbol {v} = (v_x,v_y,v_z)$ is the velocity field and $\boldsymbol {b}$ is the induced magnetic field, making up the two dynamical variables in this system. The two control parameters are $\boldsymbol {\varOmega }$, the global rotation vector (with magnitude $\varOmega$) and $\boldsymbol {B_0}$, the uniform background field (with magnitude $B_0$). Other definitions include the total pressure modified by rotation $p^*$, which is normalized by the constant density $\rho _0$ and the dissipation terms, $\boldsymbol {D}_v$ and $\boldsymbol {D}_b$, which could be regular viscosity and magnetic diffusion, respectively, but might also include other forms of dissipation such as drag or hypodiffusion. The exact form of these terms will be described in § 3, when the simulations are introduced. Magnetic fields are in Alfvén units, being normalized by $\sqrt {\rho _0 \mu _0}$, where $\mu _0$ is the magnetic permeability. Finally, $\boldsymbol {f}$ is a body force, which will be used to inject energy into the velocity field.

The inviscid and perfectly conducting system conserves the total energy,

(2.4)\begin{equation} E = \tfrac{1}{2} \int(\boldsymbol{v}^2+\boldsymbol{b}^2) \, \textrm{d}^3x. \end{equation}

However, when $\boldsymbol {\varOmega }$ and $\boldsymbol {B_0}$ are collinear, this system also conserves what is known as the parallel helicity (Shebalin Reference Shebalin2006) or hybrid helicity (Galtier Reference Galtier2014; Menu et al. Reference Menu, Galtier and Petitdemange2019). The collinear system has received considerable attention – favoured over the misaligned case in part due to its extra conserved quantity and the potential relevance of its cascade for dynamo action (Shebalin Reference Shebalin2006; Menu et al. Reference Menu, Galtier and Petitdemange2019). It also possesses simplified linear wave solutions which have been used to develop a weak wave turbulence theory (Galtier Reference Galtier2014; Bell & Nazarenko Reference Bell and Nazarenko2019). Here, we will not be considering the collinear case, and so only the total energy will be conserved in our study of $B\varOmega$-MHD in § 3. Although waves are certainly present in our system, our work concerns the strongly turbulent dynamics of energy cascades (present partly in the zero frequency modes of the system). See Appendix A for the dispersion relation of waves in the misaligned case.

Most studies, with rotation and background magnetic field aligned or not, have focused on how rotation and a moderate background field affect the decay of kinetic and magnetic energies in unforced simulations (Lehnert Reference Lehnert1955; Favier, Godeferd & Cambon Reference Favier, Godeferd and Cambon2012; Baklouti et al. Reference Baklouti, Khlifi, Salhi, Godeferd, Cambon and Lehner2019; Bell & Nazarenko Reference Bell and Nazarenko2019). Menu et al. (Reference Menu, Galtier and Petitdemange2019) investigated the sensitivity of the cascade of hybrid helicity for various rotation and guide field alignments in forced-dissipative simulations. We consider the effects of rotation and a misaligned background magnetic field on the two-dimensionalization of the flow and the energy cascade, including the limits of strong rotation and strong background magnetic field.

In our study, the rotation and background magnetic field vectors are perpendicular to each other, namely, we have chosen $\boldsymbol {\varOmega } = \varOmega \boldsymbol {\hat {z}}$ and $\boldsymbol {B_0} = B_0 \boldsymbol {\hat {x}}$, the extension to an arbitrary angle between $\boldsymbol {\varOmega }$ and $\boldsymbol {B_0}$ being discussed in § 5. The turbulence is maintained at a statistically steady state by a forcing which inputs energy at a mean rate $I$ at a length scale $1/k_f$ (see details in § 3). As a result, there is an emergent velocity scale $U$ defined to be $U^3 \equiv I k_f^{-1}$, that we compare with the background field as a measure of its strength, the inverse Alfvén Mach number

(2.5)\begin{equation} M^{{-}1} \equiv \frac{B_0}{U}. \end{equation}

This dimensionless number can also be thought of as a measure of how the third term on the right-hand side of (2.1) (the Lorentz force) and the first term on the right-hand side of (2.2) compare with the advection terms in each respective equation, which would determine whether or not the background field affects the dominant dynamics. When $M^{-1} \gg 1$ the Lorentz force acts to constrain the velocity and induced magnetic fields so that they do not vary along the $x$-direction and most of the energy lies in the $k_x = 0$ modes, resulting in 2-D MHD (Montgomery & Turner Reference Montgomery and Turner1981; Nazarenko Reference Nazarenko2007; Alexakis Reference Alexakis2011; Bigot & Galtier Reference Bigot and Galtier2011; Sujovolsky & Mininni Reference Sujovolsky and Mininni2016). It is important to note that, while the dynamics depends on $y, z$ and $t$, all vector components can be non-zero in periodic domains. This is called two-dimensional, three-component (2D3C) dynamics (Biferale, Buzzicotti & Linkmann Reference Biferale, Buzzicotti and Linkmann2017). If the induced magnetic field is not directly forced (as is the case in our study), this results in 2D3C HD and an inverse cascade of horizontal kinetic energy (Alexakis Reference Alexakis2011; Sujovolsky & Mininni Reference Sujovolsky and Mininni2016). All of our simulations lie in the regime of strong background magnetic field, $M^{-1} \gg 1$, making the rotation rate the main control parameter in our study. Does this asymptotic regime survive in the presence of rotation?

The relative strength of rotation is measured by the inverse Rossby number

(2.6)\begin{equation} Ro^{{-}1} \equiv \frac{2 \varOmega}{k_f U}. \end{equation}

This number measures the relative importance of the second term on the right-hand side of (2.1) (the Coriolis force) to the advection term, which would determine whether or not the rotation affects the dynamics. Unlike the background magnetic field, the Coriolis force only directly affects the velocity field. For non-stratified rapidly rotating HD in the absence of any magnetic field, $Ro^{-1} \gg 1$, the strong Coriolis force acts to constrain the flow such that it does not vary along the $z$-direction and most of the energy lies in the $k_z = 0$ modes (Smith & Waleffe Reference Smith and Waleffe1999; Mininni & Pouquet Reference Mininni and Pouquet2010; Gallet Reference Gallet2015; Vallis Reference Vallis2017; Buzzicotti et al. Reference Buzzicotti, Aluie, Biferale and Linkmann2018), which results in 2D3C HD where the dynamics depends only on $x, y$, and $t$. If the fluid is ionized and initialized with a non-zero seed magnetic field, rapid rotation does not necessarily result in 2D3C HD because there is no direct constraint on the induced magnetic field. Instead, if the transverse velocity component does not vanish, rapidly rotating dynamos are formed with $z$-dependent induced magnetic fields (Otani Reference Otani1993; Smith & Tobias Reference Smith and Tobias2004; Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016b; Seshasayanan, Gallet & Alexakis Reference Seshasayanan, Gallet and Alexakis2017; Tobias Reference Tobias2021). However, since our base state is the $x$-independent 2D3C HD regime found when $M^{-1} \gg 1$, rapid rotation is expected to act to constrain the flow and prevent it from varying in the $z$-direction. Note that, in this configuration, rotation is in the plane of the 2-D dynamics, not out of the plane as is often the case when it itself is the cause of the bidimensionalization. Since rotation is now in the plane of the 2-D velocities, the Coriolis force is expected to deflect horizontal velocities out of the plane, as will be discussed in § 4 when we introduce a reduced model for this system following Montgomery & Turner (Reference Montgomery and Turner1981).

Our goal in this study is to investigate the effects that in-plane rotation has on the 2-D flow caused by a strong background magnetic field. In the next section we will describe results from direct numerical simulations of the $B\varOmega$-MHD system for various rotation rates, paying particular attention to the resulting energy cascade and morphology of the flow field.

3. Strong background field limit: 3-D $B\varOmega$-MHD simulations

Equations (2.1)–(2.3) were solved numerically in a triply periodic domain of side length $2 {\rm \pi}L$ using the Geophysical High-Order Suite for Turbulence (GHOST) code (Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011). The dissipation terms, $\boldsymbol {D}_v$ and $\boldsymbol {D}_b$, each consist of a ‘hyperviscosity’ and a large-scale dissipation term called ‘hypoviscosity’. The hyperviscosity replaces the regular viscous and magnetic diffusion terms with a Laplacian of a higher order, in our case $\nabla ^2 \rightarrow -\nabla ^4$. This higher order allows for the possibility of forcing at smaller length scales while still properly resolving the smallest scales at moderate resolutions. As long as the order of the Laplacian is not very large, hyperviscosity has been shown to have no significant effect on the turbulent properties of 3-D turbulence, and we expect the same to be the case for our work (Agrawal et al. Reference Agrawal, Alexakis, Brachet and Tuckerman2020). The hypoviscosity, which would appear as $\nu _- \nabla ^{-2} \boldsymbol {v}$ on the right-hand side of (2.1) and as $\eta _- \nabla ^{-2}\boldsymbol {b}$ on the right-hand side of (2.2), acts as a large-scale dissipation term. The resulting expressions for $\boldsymbol {D}_v$ and $\boldsymbol {D}_b$ are,

(3.1)\begin{equation} \left. \begin{gathered} \boldsymbol{D}_v ={-}\nu \nabla^4 \boldsymbol{v} + \nu_- \nabla^{{-}2}\boldsymbol{v}, \\ \boldsymbol{D}_b ={-}\eta \nabla^4 \boldsymbol{b} + \eta_- \nabla^{{-}2}\boldsymbol{b}, \end{gathered}\right\} \end{equation}

where $\nu$ is the kinematic ‘hyper’-viscosity, $\eta = (\mu _0 \sigma )^{-1}$ is the magnetic ‘hyper’- diffusivity, $\sigma$ is the electrical conductivity, $\nu _-$ is the ‘hypo’-viscosity and $\eta _-$ is the magnetic ‘hypo’-diffusivity. Should an inverse cascade of a conserved quantity occur, this term ensures that no condensate forms, which would otherwise affect the cascades and inertial ranges (Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007; Xia et al. Reference Xia, Punzmann, Falkovich and Shats2008; Gallet & Young Reference Gallet and Young2013; Alexakis & Biferale Reference Alexakis and Biferale2018; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2018; van Kan & Alexakis Reference van Kan and Alexakis2019). This is done by choosing the coefficients $\nu _-$ and $\eta _-$ such that the kinetic and magnetic energy at the largest scales is smaller than that of the next largest scales. The modified GHOST code which includes these alternative dissipative terms can be found in a public Github repository (Benavides Reference Benavides2019). It is a standard parallel pseudo-spectral code with a fourth-order Runge–Kutta scheme for time integration and a two-thirds dealiasing rule. The numerical model is non-dimensionalized by $L$ and the forcing amplitude $f_0$, so that the wavenumbers $\boldsymbol {k}$ correspond to mode numbers of the domain and the forcing amplitude is one. The 3-D forcing $\boldsymbol {f}$ is isotropic and constant in time, comprising of a summation of cosines with wavenumbers in the range $8 < |\boldsymbol {k}| < 10$, making $k_f$ = 9, and random phases. The forcing wavenumber range is chosen in an attempt to properly resolve both an inverse cascade and a forward cascade. $I \equiv \langle \boldsymbol {f} \boldsymbol {\cdot } \boldsymbol {v} \rangle$ is the space- and time-averaged energy injection rate, where $\langle \boldsymbol {\cdot } \rangle$ represents a space and time average. We do not force the induced magnetic field.

All runs, unless otherwise stated, are in the large background field regime with $M^{-1} \approx 84$. We find this value to be large enough to produce the expected two-dimensionalization in the absence of rotation (figure 1a). Larger background magnetic field values result in significant restrictions in the time step which would limit our ability to perform the same parameter sweep. The Reynolds and magnetic Reynolds numbers, defined, respectively, as $Re \equiv U / k_f^3 \nu$ and $Re_m \equiv U / k_f^3 \eta$ when considering hyperviscous and hyperdiffusive terms as we do, measure the relative strength of the advection terms compared with the hyperviscous and magnetic hyperdiffusion terms. For the simulations we performed, the Reynolds and magnetic Reynolds numbers were large (approximately 300) and equal to each other, i.e. the magnetic Prandtl number is set to one. We performed 14 runs at $M^{-1} \approx 84$ but at different values of $Ro^{-1}$, ranging from $Ro^{-1} = 0$ to $Ro^{-1} = 27$. All averages and snapshots were taken in statistically steady states. See table 1 for details of the simulations and a description of how we measured the non-dimensional numbers.

Figure 1. Snapshots of the field-aligned vorticity $\omega = \hat {x}\boldsymbol {\cdot }(\boldsymbol {\nabla } \times \boldsymbol {v} )$, representing, from top left to bottom right, Regimes I (a), II (b), III (c) and IV (d), as rotation rate is increased. The red colours represent positive vorticity whereas the blue represents negative vorticity. Regime I is characterized by a bidirectional cascade, Regime II a purely forward cascade, Regime III the formation of strong shear layers (seen here in the middle of the domain) and Regime IV magnetically active shear layers. Regimes I–III have a negligible induced magnetic energy, unlike Regime IV whose magnetic energy dominates the dynamics (figure 2).

Table 1. A summary of the runs performed for this work. All runs have hyper- and hypo-viscosity of the same order (section 3). For runs with a magnetic field, $\mu = \nu$ and $\mu _- = \nu _-$. The simulations used Alfvénic units so that $B_0 / \sqrt {\rho _0 \mu _0} \rightarrow B_0$ and the other $\rho _0$ was absorbed into the pressure. The values are non-dimensionalized by $L$ and the forcing amplitude $f_0$ (or $I_0$ for random forcing), so that $k_f$ is the forcing mode number and the forcing amplitude (for constant amplitude forcing) or energy injection rate (for the random forcing) are both 1 in these units. The typical velocity, $U$, was calculated after the fact for each run using $U^3 \equiv Ik_f^{-1}$, where $I \equiv \langle \boldsymbol {f} \boldsymbol {\cdot } \boldsymbol {v} \rangle$ is the time- and space-averaged energy injection rate. The count is the number of runs in that set.

In this study, we are partly concerned with the behaviour of the energy cascade as rotation is varied. We expect the presence of a bidirectional cascade, where a fraction of the energy input by the forcing goes to large scales and the rest goes to small scales. As such, we define a measure for the fraction of energy that goes to large scales in the form of kinetic energy, $\varepsilon _-$, and that which goes to small scales in the form of kinetic energy $\varepsilon$ and magnetic energy $\varepsilon _b$. Since the large-scale magnetic energy dissipation is practically zero for every simulation performed, we ignore it from our analysis, as it plays no role. These measures are based on the dissipation rates from each of the three dissipation terms, and are defined in the following way:

(3.2a–c)\begin{equation} \varepsilon_-{\equiv} \nu_- \langle | \nabla^{{-}1} \boldsymbol{v} |^2 \rangle/I, \quad \varepsilon \equiv \nu \langle | \nabla^2 \boldsymbol{v} |^2 \rangle/I, \quad \varepsilon_b \equiv \eta \langle | \nabla^2 \boldsymbol{b} |^2 \rangle/I.\end{equation}

Energy balance at steady state tells us that $\varepsilon _- + \varepsilon + \varepsilon _b = 1$. In the limit of large Reynolds number and large forcing wavenumber, none of the energy injected is dissipated at the forcing scale and proper inertial ranges are formed. In this case, the dissipation rate at large scales represents the fraction of energy cascading to large scales, and similarly for the dissipation rate at small scales. Our runs do not reach these idealized limits. The lack of scale separation between the forcing and large-scale dissipation will manifest itself in two related ways in this paper: (i) the large-scale dissipation rate will remain non-zero despite zero average inverse cascade, because some energy that is being injected at $k_f$ will be dissipated by the large-scale dissipation mechanism ($\varepsilon _- \leq 0.1$ for the 3-D runs), and (ii) when layers form in Regime III, both the 3-D runs and 2D3C runs show an increase in large-scale dissipation rate, not because of the presence of an inverse cascade, but because the layers form coherent structures near the forcing scale, their energy grows and hence a stronger large-scale dissipation rate is achieved. These jumps in $\varepsilon _-$ denote the presence of shear layers, as discussed in § 4. To complement these estimates for energy cascades, we will look at the normalized spectral energy flux

(3.3)$$\begin{gather} \varPi_{KE}(k) \equiv \langle \boldsymbol{v}^{{<}k} \boldsymbol{\cdot} (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} )\rangle/I, \end{gather}$$

(3.4)$$\begin{gather}\varPi_{ME}(k) \equiv{-}\langle \boldsymbol{v}^{{<}k} \boldsymbol{\cdot} ((\boldsymbol{B_0} + \boldsymbol{b}) \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{b} )\rangle/I + \langle \boldsymbol{b}^{{<}k} \boldsymbol{\cdot} (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{b} - (\boldsymbol{B_0}+\boldsymbol{b})\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v})\rangle/I, \end{gather}$$

where $\boldsymbol {v}^{< k}$ stands for a filtering of the velocity $\boldsymbol {v}$ in Fourier space so that only the wavenumbers with modulus smaller than $k$ are kept. The flux $\varPi (k)$ expresses the rate at which energy is flowing out of scales larger than $2{\rm \pi} /k$ due to nonlinear interactions, normalized by the energy injection rate. Therefore, if energy is going from large to small scales, the energy flux will be positive, and vice versa. Finally, to quantify the amount and type of energy at each scale, we will also look at the energy spectra

(3.5a,b)\begin{equation} E_{KE}(k) \equiv \tfrac{1}{2}\sum_{|\boldsymbol{k}|=k}|\hat{\boldsymbol{v}}|^2(\boldsymbol{k}), \quad E_{ME}(k) \equiv \tfrac{1}{2}\sum_{|\boldsymbol{k}|=k}|\hat{\boldsymbol{b}}|^2(\boldsymbol{k}), \end{equation}

where $\hat {\boldsymbol {v}}$ denotes the Fourier transform of $\boldsymbol {v}$.

Beginning from quasi-2-D turbulence on the $y$–$z$ plane at zero rotation, we find four distinct regimes as we increased rotation (figure 1). Although not so apparent in the 3-D simulations, these regimes are separated by seemingly sharp transitions, whose boundaries are determined in § 4.

Regime I (figure 1a), defined for runs with $Ro^{-1} < 0.6$, is characterized by the presence of a bidirectional cascade. This can be seen in figure 2 as a non-zero large-scale dissipation rate as well as in figure 3(e), where the spectral energy transfers show that about half of the energy injected by the forcing goes to large scales (negative $\varPi (k)$) and the other half goes to small scales (positive $\varPi (k)$). The fraction of energy that goes to larger scales decreases with increasing rotation (figure 2). At zero rotation we do not have a purely inverse cascade ($\varepsilon _- \approx 1$) due to a combination of finite background magnetic field strength and, as we will see in § 4, the fact that we are forcing the out-of-plane velocity which acts as a passive scalar in the two-dimensionalized dynamics, thus contributing to a forward energy flux (Campagne et al. Reference Campagne, Gallet, Moisy and Cortet2014; Biferale et al. Reference Biferale, Buzzicotti and Linkmann2017). Therefore, at zero rotation rate, the system is undergoing two independent cascades: an inverse energy cascade of horizontal kinetic energy and a forward cascade of the out-of-plane kinetic energy. If we were to force only the horizontal velocity components in the $k_x = 0$ wavenumber plane, we would expect to see $\varepsilon _- \approx 1$ at zero rotation. Figure 3(a) shows the kinetic and magnetic energy spectra, which demonstrates that the magnetic energy is orders of magnitude smaller than the kinetic energy (particularly at large scales) and that the largest scales have the most energy, providing further confirmation of the presence of an inverse cascade. The spike of magnetic energy at the forcing scale is due to the excitation of Alfvén waves from the isotropic forcing. The eddy length scales seen in figure 1(a) are set by a combination of the energy injection and the large-scale hypoviscosity coefficient. Regime II (figure 1b), defined for runs with $0.6 < Ro^{-1} < 2.1$, is characterized by a purely forward cascade of energy (figures 2 and 3f). This may come as a surprise, given that the dynamics is two-dimensional. The reason for this seemingly contradictory state is that, while two-dimensional, all three velocity components are active in the dynamics and, furthermore, are coupled together with rotation. This results in a set of reduced 2D3C equations which no longer conserve enstrophy, making a forward cascade of energy possible. The rotating 2D3C system will be discussed and explored numerically in § 4.

Figure 2. Dissipation rates normalized by the energy injection rate as a function of rotation rate measured by the inverse Rossby number $Ro^{-1}$. The blue solid line shows the ohmic dissipation rate due to the magnetic diffusion term, $\varepsilon _b$, the green dashed line shows the viscous dissipation rate, $\varepsilon$, and the red dash-dotted line shows the large-scale dissipation rate due to the hypoviscosity, $\varepsilon _-$. Each regime is labelled at the top, and the vertical dashed lines represent boundaries between regimes, chosen based on the 2-D runs in § 4. Stars represent runs whose snapshots are shown in figure 1.

Figure 3. The time-averaged energy spectra (a-d) and spectral energy flux (e-h) for each regime found in our simulations. The blue dashed line shows the magnetic components, either $E_{ME}$ or $\varPi _{ME}$, the orange dash-dotted line shows the kinetic components, $E_{KE}$ or $\varPi _{KE}$, and the solid black line shows their sum. The grey box represents the forcing range. These are from the same simulations shown in figure 1 and which are starred in figure 2.

Regime III (figure 1c), defined for runs with $2.1 < Ro^{-1} < 7.5$, is characterized by the formation of strong shear layers along the $y$-direction, consisting of uniform velocity in the $x$–$z$ plane. The shear layers form when the rotational constraint on the dynamics at large scales becomes sufficiently large, requiring that $\partial _z \boldsymbol {v} \approx 0$ at those scales. The combination of $\partial _z = \partial _x = 0$ and incompressibility implies that $v_y = 0$ (since we are in a periodic domain), and thus that the last remaining component of the nonlinear advection term $v_y \partial _y = 0$ and nonlinearities are suppressed at large scales. Because of the suppressed nonlinearity at large scales, these shear layers form coherent structures that are fed directly by the forcing but that do not transfer that energy away, causing a build up of energy (not shown). The energy in the layers builds until a combination of the large-scale dissipation (figure 2) and the nonlinear term (figure 3g) are able to remove energy from those scales. Regimes I–III have negligible induced magnetic energy, as is observed in simulations of MHD with a strong background field (Alexakis Reference Alexakis2011; Sujovolsky & Mininni Reference Sujovolsky and Mininni2016), and so the induced magnetic field plays an insignificant role in the dynamics. The magnetic fluctuations are also much smaller than $B_0$ – less than 0.5 % of $B_0$ in Regimes I–III.

This changes, however, in Regime IV (figure 1d), defined for runs with $Ro^{-1} > 7.5$, where we have found the activation and growth of the induced magnetic field, which dominates both the energy as well as the nonlinear energy transfers (figure 2). The nonlinear advection term in the momentum equation is suppressed for practically all scales (figure 3h), leading to laminar-like shear-layer structures (figure 1d) and a turbulent magnetic field which is responsible for the nonlinear transfers of energy across scales, via the Lorentz force and the magnetic induction equation. The shear-layer spacing in figure 1(d) is set by the forcing scale. In this regime, significant induced magnetic field fluctuations occur both parallel and perpendicular to the background magnetic field, with a magnitude of approximately 3 % of $B_0$.

We expect the boundaries between Regimes I–III to be independent of $M^{-1}$, as they are part of the asymptotic 2D3C HD, whose sole parameter is the rotation rate. We confirm this in the next section, which deals specifically with this asymptotic set of equations, by showing that the regime transitions happen for the same values of $Ro^{-1}$. The transition from Regime III to IV is of a different nature and represents a breakdown of the hydrodynamic behaviour found for lower rotation rates. This transition is found to be $M^{-1}$-dependent, and will be discussed briefly in § 5.

4. Comparison with rotating 2D3C model

Regimes I–III can be better understood by considering the asymptotic limit of (2.1)–(2.3) when taking $M^{-1} \rightarrow \infty$ and keeping $Ro^{-1} \sim {O}(1)$ and the domain size fixed. The choice of keeping the domain size fixed is based on the fact that we are motivated mostly by astrophysical settings in confined geometries, in the presence of a strong background magnetic field. We acknowledge, however, that in many other astrophysical settings, such as the extended atmosphere of stars or the interstellar medium, a confined geometry may not be the best representative system to study. In such systems, a more appropriate limit might include taking the domain size to infinity at the same time as the $M^{-1}$, so as to prevent the exact two-dimensionalization of the flow (Thess & Zikanov Reference Thess and Zikanov2007; Gallet & Doering Reference Gallet and Doering2015). The limiting equations in this case would resemble more the Reduced MHD system, derived for tokamaks but used also to study some astrophysical systems, in which the flow is highly anisotropic, yet still three-dimensional (Strauss Reference Strauss1976; Oughton, Matthaeus & Dmitruk Reference Oughton, Matthaeus and Dmitruk2017).

Our limiting procedure, with fixed domain size, is similar to that done in Montgomery & Turner (Reference Montgomery and Turner1981), with the exception that we include the Coriolis term, and so we only briefly discuss it here. Assuming a background magnetic field in the $x$ direction and a rotation axis along the $z$ direction, the process results in a set of three dynamical equations and one nonlinear constraint for the three variables: $\psi (y,z,t)$ the streamfunction for the in-plane velocities, $v_x(y,z,t)$ the out-of-plane velocity and $A(y,z,t)$ the potential for the in-plane magnetic field. This novel constraint, which results from the presence of the Coriolis term, states that either $\delta A / \delta v_x =0$ or $v_y = \partial _z \psi = 0$, where the former is the functional derivative of $A$ with respect to $v_x$. Our 3-D simulations from § 3 seem to be consistent with these constraints, where, in Regimes I–III for $Ro^{-1} < 7.5$, we have $A\approx 0$ but $v_y \neq 0$ and, in Regime IV for $Ro^{-1} > 7.5$, we have $A\neq 0$ but $v_y = \partial _z \psi \approx 0$. For the purposes of studying the reduced dynamics of Regimes I–III, we assume $A = 0$, knowing that it would not capture the transition to Regime IV. The resulting equations form the 2D3C system with in-plane rotation

(4.1)$$\begin{gather} \frac{\partial v_x}{\partial t} + [v_x,\psi] = 2 \varOmega \frac{\partial \psi}{\partial z} + \nu \nabla_\perp^2 v_x + f_x, \end{gather}$$

(4.2)$$\begin{gather}\frac{\partial \omega}{\partial t} + [\omega,\psi] = 2 \varOmega \frac{\partial v_x}{\partial z} + \nu \nabla_\perp^2 \omega + f_\omega, \end{gather}$$

where $[F,G] \equiv \partial _y F \partial _z G - \partial _y G \partial _z F = 0, \nabla _\perp = (0,\partial _y,\partial _z), \omega = \hat {x}\boldsymbol {\cdot } (\boldsymbol {\nabla } \times \boldsymbol {v} ) = -\nabla ^2_\perp \psi$ is the out-of-plane vorticity (of the in-plane velocities), $\perp$ implies the directions perpendicular to the background magnetic field and $f_\omega = \hat {x}\cdot (\nabla _\perp \times \boldsymbol {f}_\perp )$. One could equivalently arrive at (4.1) and (4.2) by taking rotating 3-D HD and requiring that the velocity field does not depend on $x$. If considering an arbitrary angle between the background field and rotation, only the perpendicular projection of the rotation vector on the background field enters the model. For example, supposing without loss of generality that $\boldsymbol {B_0} = B_0 \boldsymbol {\hat {x}}$ and $\boldsymbol {\varOmega } = \varOmega (\sin (\theta ) \boldsymbol {\hat {z}} + \cos (\theta ) \boldsymbol {\hat {x}})$, then the Coriolis terms on the right-hand side of (4.1) and (4.2) will be multiplied by $\sin (\theta )$. This asymptotic model is in the same spirit as some of the magnetized quasigeostrophic models used in astrophysical applications (Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015; Maffei et al. Reference Maffei, Calkins, Julien and Marti2019), but it is important to note that here we have taken $M^{-1} \rightarrow \infty$ while keeping $Ro^{-1} \sim {O}(1)$, whereas the magnetized quasigeostrophic models take the rapidly rotating limit first. As is discussed in § 5, these limits are not expected to be the same.

The Coriolis force now couples the two equations together, making what would otherwise be a passive tracer into an active one. In fact, for non-zero rotation, it can be shown that the 2D3C rotating system conserves kinetic energy and helicity

(4.3)$$\begin{gather} KE = \tfrac{1}{2}\int v_x^2 + |\boldsymbol{\nabla} \psi|^2 \, \textrm{d}^2 x, \end{gather}$$

(4.4)$$\begin{gather}H = \int v_x \omega \, \textrm{d}^2 x. \end{gather}$$

These are the same conserved quantities as in 3-D HD, but we emphasize that the dynamics is two-dimensional and occurs on the $y$–$z$ plane. This is in contrast to the case of zero rotation, where the system conserves (separately) the in-plane kinetic energy $\int |\boldsymbol {\nabla } \psi |^2 \, \textrm {d}^2x$ and the out-of-plane kinetic energy $\int v_x^2 \, \textrm {d}^2x$, as well as the enstrophy, $\int \omega ^2 \, \textrm {d}^2 x$. The conservation of enstrophy can be shown to prevent the existence of a forward cascade of in-plane kinetic energy (Fjortoft Reference Fjortoft1953; Kraichnan Reference Kraichnan1967; Alexakis & Biferale Reference Alexakis and Biferale2018). Without the restriction of enstrophy conservation, though, the kinetic energy may go downscale in a forward cascade, even if one does not force the out-of-plane component.

We solve (4.1) and (4.2), with modified hyper- and hypo-viscosities as in the 3-D simulations, numerically in a doubly periodic domain of side length $2 {\rm \pi}L$ using the 2-D predecessor of GHOST. The code can be found in a public Github repository (Benavides Reference Benavides2020). Unlike the 3-D runs, whose forcing function had a constant amplitude in time, the 2D3C runs have random (white-in-time) forcing. At each time step, a wavenumber $\boldsymbol {k}_r$ of magnitude $k_f$ is chosen at random and $\hat {f}_\omega (\boldsymbol {k})$ (Fourier transform of $f_\omega$) is set to zero everywhere except for at $\boldsymbol {k}_r$, where it had a magnitude $k_f \sqrt {2 I_0 / \Delta t}$ (Chan, Mitra & Brandenburg Reference Chan, Mitra and Brandenburg2012). This has the effect of setting the energy injection rate for the in-plane flow to be $I = \langle \psi f_\omega \rangle = I_0$ on average, with $I_0$ being an input parameter of the simulation. The same forcing is applied for $f_x$, but with an amplitude of $\sqrt {2 I_0 / \Delta t}$ instead, giving the same results. Therefore, half of the energy is injected into the in-plane flow and the other half in the out-of-plane velocity. We non-dimensionalize all dynamical variables as before, using $L$ and now the energy injection rate parameter $I_0$. For all of the runs reported $k_f = 12$. See table 1 for details on the runs.

The goal of these simulations is to reproduce the parameter sweep performed in § 3, but with the added advantage of working with a 2-D code, thus allowing a larger quantity of runs, higher resolutions (larger Reynolds numbers, around 600) and longer time integration. We have performed 23 runs, with $Ro^{-1}$ ranging from 0 to approximately 5, at four times the horizontal resolution. Our results confirm the presence of Regimes I–III, going from a bidirectional cascade to a forward cascade to a shear-layer configuration (figure 4).

Figure 4. Snapshots of the out-of-plane vorticity, $\omega = \hat {x}\boldsymbol {\cdot } (\boldsymbol {\nabla } \times \boldsymbol {v} ) = -\nabla ^2_\perp \psi$, for the 2D3C rotating simulations, representing, from left to right, Regimes I (a), II (b) and III (c), as rotation is increased. We see striking similarities to figure 1, confirming that the asymptotic 2D3C model captures the different regimes found in the 3-D simulations of § 3.

At zero rotation we see a bidirectional cascade with half of the injected energy going to large scales and half going to small scales (figure 5), similar to what was found in the 3-D runs (figure 2). For the 2D3C rotating system this is the case because of the choice of forcing, which injects half of the energy to the in-plane flow and the other half to the out-of-plane velocity. Since the two flows are completely decoupled at zero rotation, they each follow the standard behaviour observed in 2-D and passive tracer turbulence, that is, an inverse and forward cascade of energy, respectively. As we increase rotation, the Coriolis force couples the two fields, enstrophy is no longer conserved, and the in-plane velocities no longer cascade all the injected energy to large scales, resulting in a bidirectional cascade with decreasing inverse energy flux. There is an approximately linear approach to zero inverse energy flux, and at $Ro^{-1} \approx 0.6$ there is a transition to a purely forward cascade. With a larger number of simulations, Regimes I and II are much more clearly separated, and their transition appears to be sharp (figure 5). This transition is qualitatively similar to other bidirectional to forward cascade transitions seen in other studies and could hint at a universal mechanism responsible for this transition (Seshasayanan, Benavides & Alexakis Reference Seshasayanan, Benavides and Alexakis2014; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016a; Benavides & Alexakis Reference Benavides and Alexakis2017; van Kan & Alexakis Reference van Kan and Alexakis2020).

Figure 5. Dissipation rates normalized by the energy injection rate as a function of rotation rate measured by the inverse Rossby number $Ro^{-1}$. The green dashed line shows the viscous dissipation rate, $\varepsilon$, the red dash-dotted line shows the large-scale dissipation rate due to the hypoviscosity, $\varepsilon _-$ and the purple dash-dotted line shows the same but for hysteresis runs initialized with layers. Each regime is labelled at the top, and the vertical dashed lines represent boundaries between regimes. These denoted boundaries are placed at the same value of $Ro^{-1}$ as those seen in figure 2. Triangles represent runs whose snapshots are shown in figure 4.

Upon further increase of the rotation, the forward cascade regime (Regime II) transitions to a shear-layer configuration (figure 4c), entering Regime III. This corresponds to the case when the Coriolis force dominates at large scales, making the dominant balance in (4.1) and (4.2) $\partial _z \psi \approx \partial _z v_x \approx 0$, hence the layers. There are a few differences in the morphology of the shear layers seen for these runs, compared with Regime III in the 3-D simulations (figure 1c). Here they take up the whole domain and also appear to equilibrate at scales larger than the forcing, through a series of mergers (not shown). Neither of these characteristics are seen in the shear layers of the 3-D simulations. We believe this is due to a few factors, including the longer integration times and the change in forcing. A surprising feature of this transition, revealed by the better-resolved parameter sweep, is that it is discontinuous (figure 5). (An increase in large-scale dissipation marks this transition not because an inverse cascade forms (a weakness of this measure), but because of a lack of separation of scales. The layers form at or near the forcing scale and remain there as coherent structures, fed directly by the forcing, resulting in a build up of energy at those scales. This, in turn, results in a larger dissipation rate from the large-scale dissipation. If we were to perform runs at a larger $k_f$, this effect would disappear. The discontinuous transition is also observed in the kinetic energy, which is not shown.) Discontinuities are a characteristic of subcritical bifurcations, which should also display hysteresis. By initializing in the layered regime, we confirmed the presence of hysteresis as we reduced the rotation rate (figure 5 inset).

Despite differences in the forcing, Reynolds numbers and values of $M^{-1}$, the regime transitions seem to occur for the same values of $Ro^{-1}$, suggesting that the rotating 2D3C system successfully describes the dynamics observed in the 3-D simulations from § 3 and that Regimes I–III are robust properties of the system. The 2-D asymptotic system has allowed us to perform a more detailed parameter sweep of this parameter space, and has revealed sharp transitions and non-trivial behaviour near those transitions which we did not anticipate from the 3-D simulations.