Magnetorotational Dynamo Action in the Shearing Box

TL;DR

This study investigates how the geometry of the shearing box influences magnetic dynamo action driven by the magnetorotational instability, highlighting the importance of box aspect ratio and turbulent magnetic flux distribution.

Contribution

It demonstrates that the dynamo's efficiency depends on the box aspect ratio and provides a physical explanation involving magnetic flux spreading and turbulence.

Findings

Dynamo action is easier in tall boxes than in long boxes.

Large magnetic Reynolds number is necessary for dynamo in both geometries.

Intermittency of azimuthal magnetic flux correlates with dynamo efficiency.

Abstract

Magnetic dynamo action caused by the magnetorotational instability is studied in the shearing-box approximation with no imposed net magnetic flux. Consistent with recent studies, the dynamo action is found to be sensitive to the aspect ratio of the box: it is much easier to obtain in tall boxes (stretched in the direction normal to the disk plane) than in long boxes (stretched in the radial direction). Our direct numerical simulations indicate that the dynamo is possible in both cases, given a large enough magnetic Reynolds number. To explain the relatively larger effort required to obtain the dynamo action in a long box, we propose that the turbulent eddies caused by the instability most efficiently fold and mix the magnetic field lines in the radial direction. As a result, in the long box the scale of the generated strong azimuthal (stream-wise directed) magnetic field is always comparable to the scale of the turbulent eddies. In contrast, in the tall box the azimuthal magnetic flux spreads in the vertical direction over a distance exceeding the scale of the turbulent eddies. As a result, different vertical sections of the tall box are permeated by large-scale nonzero azimuthal magnetic fluxes, facilitating the instability. In agreement with this picture, the cases when the dynamo is efficient are characterized by a strong intermittency of the local azimuthal magnetic fluxes.

Tables

Table 1. Table 1: Setups examined in this study. The final case (III) is decaying until t ≈ 600 𝑡 600 t\approx 600 , and the numbers given for E 𝐸 E and α 𝛼 \alpha are the values averaged over short time intervals around the indicated times.

Case	$L_x\times L_y\times L_z$	${\nu }^{-1}$	$E$	$\alpha (\times {10}^{-2})$
I	$1\times 4\times 4$	5000	0.015	0.015
II	$1\times 4\times 16$	5000	0.15	0.0081
IIIa (t=222)	$2\times 4\times 1$	45000	0.017	0.28
IIIb ($t>800$)	$2\times 4\times 1$	45000	0.0015	0.0037

Equations

$${\partial }_t\mathbf{v}=q{\Omega }_{0}x{\partial }_y\mathbf{v}-(\mathbf{v}\cdot \nabla )\mathbf{v}-\nabla P+\mathbf{b}\cdot \nabla \mathbf{b}+\nu {\nabla }^2\mathbf{v}-$$

$$-2{\mathbf{\Omega }}_0\times \mathbf{v}+q{\Omega }_0{v}_x{\mathbf{e}}_y,$$

$${\partial }_t\mathbf{b}=q{\Omega }_{0}x{\partial }_y\mathbf{b}+\nabla \times (\mathbf{v}\times \mathbf{b})+\eta {\nabla }^2\mathbf{b}-q{\Omega }_0{b}_x{\mathbf{e}}_y,$$

$$\frac{d}{dt}⟨\frac{v^2}{2}+\frac{b^2}{2}⟩=\tilde{\alpha }-\nu ⟨{\mathbf{\omega }}^2⟩-\eta ⟨{\mathbf{j}}^2⟩,$$

$$\frac{d}{dt}⟨\mathbf{v}\cdot \mathbf{b}+(2-q){\mathbf{\Omega }}_0\cdot \mathbf{A}⟩=-(\nu +\eta )⟨\mathbf{\omega }\cdot \mathbf{j}⟩,$$

$$\frac{d}{dt}⟨\mathbf{b}\cdot \mathbf{A}⟩=-2\eta ⟨{\mathbf{j}}^2⟩.$$

Full text

††pubyear: 2017

1 Introduction

Magnetorotational instability (MRI), the instability of a differentially rotating, conducting fluid, permeated by a weak magnetic field Velikhov (1959); Chandrasekhar (1960); Balbus & Hawley (1991), is believed to be the leading mechanism that renders astrophysical accretion flows turbulent and allows them to efficiently lose their angular momentum Balbus & Hawley (1998). The local properties of the instability and the resulting turbulence may be studied in the shearing-box approximation, which models a small box in the midplane of a disc so that the large-scale effects such as stratification and stream-line curvature can be neglected Goldreich & Lynden-Bell (1965); Hawley et al. (1995); Umurhan & Regev (2004); Fromang et al. (2007); Longaretti & Lesur (2010); Lesur & Longaretti (2011); Riols et al. (2015). The dynamics studied in the shearing box are therefore expected to be universal, that is, applicable to other astrophysical and laboratory systems possessing rotating, shearing flows (e.g., Kagan & Wheeler, 2014; Petitdemange et al., 2008; Ji & Balbus, 2013).

The properties of magnetic turbulence driven by the MRI in the shearing box are highly nontrivial. In order to discuss them we use the following standard orthogonal coordinate frame. Assuming that the box is positioned inside a disc, we chose the vertical coordinate ($z$) along the direction of the angular velocity, the coordinate ($x$) along the radial direction, and the azimuthal coordinate ($y$) along the unperturbed rotating flow. When the box is permeated by a nonzero magnetic flux in the vertical direction, the linear instability develops, eventually leading to a strongly nonlinear steady-state turbulent regime. Although the resulting turbulence is still not fully understood, it has been recently established that it resembles the standard magnetohydrodynamic turbulence Walker et al. (2016). The azimuthal magnetic field is concentrated at large scales and it plays the role of the guide field for small-scale fluctuations, whose energy spectrum and other statistical characteristic are close to that of homogeneous Alfvénic MHD turbulence Zhdankin et al. (2017). A similar result has subsequently been obtained using hybrid (kinetic ions) plasma simulations Kunz et al. (2016), which also found that sub-proton-scale MRI turbulence is consistent with the kinetic-Alfvén turbulence. When the net vertical magnetic flux is zero, strong transient amplification is still possible if the net azimuthal magnetic flux is nonzero (e.g., Balbus & Hawley, 1992; Vishniac & Diamond, 1992).

In this work we concentrate on the case when the net magnetic flux through any cross-section of the box is zero. In this zero-net-flux case, a linear instability is not possible, but the flow may become unstable for a finite initial perturbation (subcritical instability) Rincon et al. (2007); Lesur & Ogilvie (2008); Herault et al. (2011); Riols et al. (2013); Squire & Bhattacharjee (2014); Riols et al. (2016). Our recent high-resolution numerical simulations Walker et al. (2016) concentrated on the case of unit magnetic Prandtl number ($\mathrm{Pm}$) and used the “long” shearing box ($L_{x}:L_{y}:L_{z}=2:4:1$). They did not find dynamo action for the magnetic Reynolds numbers several times larger than those required for the dynamo action in a non-shearing, non-rotating box. This observation cast doubt on the existence of MRI dynamo action for $\mathrm{Pm}\leq 1$, in a stark contrast with the non-shearing, non-rotating case where the dynamo action is expected to exist for any given $\mathrm{Pm}$ as long as the magnetic Reynolds number is large enough (e.g., Boldyrev & Cattaneo, 2004; Boldyrev et al., 2005; Iskakov et al., 2007; Malyshkin & Boldyrev, 2010). A different outcome of the shearing-box dynamo simulations was however reported in Shi et al. (2016); Nauman & Pessah (2016) who used the “tall” shearing boxes, that is, boxes that vertical sizes exceed the radial sizes. Quite interestingly, they observed that the critical Prandtl number for dynamo action111We use the term “dynamo action” to refer to the self-sustained turbulence that results from the subcritical MRI instability that, in turn, reinforces the magnetic field. is reduced as long as $L_{z}/L_{x}\gtrsim 2.5$.

In order to understand the different outcomes of these studies, we have performed a series of direct numerical simulations of incompressible MHD dynamo action for varying aspect ratios of the shearing box. We found that, in agreement with Shi et al. (2016) and Nauman & Pessah (2016), the dynamo action is sensitive to the aspect ratio of the box. In our case of $\mathrm{Pm}=1$, tall boxes (stretched in the $z$-direction) rapidly exhibit dynamo action. To explain these results we propose that the turbulent eddies caused by the dynamo, efficiently fold the magnetic field lines in the radial ($x$) direction. As a consequence, in the long box the $x$-scale of the generated $B_{y}$ component of the magnetic field is always comparable to the scale of the turbulent eddies. In contrast, in the tall box the flux of $B_{y}$ can spread in the vertical direction over the distances exceeding the scales of the turbulent eddies, which are constrained by the short $x$ dimension of the box. The vertical mixing of the $B_{y}$ field is thus suppressed in tall boxes. As a result, different vertical sections of the tall box are permeated by large-scale nonzero fluxes of the azimuthal field $B_{y}$, leading to the instability.

This phenomenological picture motivated us to review one of the results from our previous work (Walker et al., 2016, case IV). In this work, the dynamo action was not observed in the long box ($L_{x}:L_{y}:L_{z}=2:4:1$) for the long time during which we integrated the solution. The initial large-scale fluctuations kept decaying for more than 200 shearing times. During this decay, the scale of the turbulence kept decreasing as well, as to maintain the balance between the linear and turbulent shear rates. For the large Reynolds number that we used, the turbulence would decay to progressively smaller scales. In this case, however, the scale of the turbulent fluctuations should eventually become sufficiently smaller than the vertical extent of the box $L_{z}$, so that according to our phenomenological picture, the dynamo should become possible. To check this hypothesis, we significantly extended the running time of simulations of (Walker et al., 2016, case IV), and after about 600 shearing times did observe the possible onset of the dynamo action.

This may reconcile the available numerical results. Based on our findings, we suggest that the dynamo action is always possible, no matter what the aspect ratio of the box is. However, the long boxes require significantly larger Reynolds number and (assuming large scale of the initial fluctuations) significantly longer running time in order to observe the dynamo action. We have also established that the dynamo action leads to a quite intermittent distribution of the azimuthal magnetic fluxes, in agreement with the proposed phenomenological picture.

2 Numerical setup

The shearing-box approximation captures the local small-scale dynamics of an accretion disc (Goldreich & Lynden-Bell, 1965). A small box centered at a radial distance $r_{0}$ and orbiting with the angular velocity $\Omega(r)$ has a shearing rate $q\equiv-(d\ln\Omega/d\ln r)_{r_{0}}$. For the Keplerian orbital flow, we have $q=3/2$, and we denote $\Omega_{0}=\Omega(r_{0})$. At the center of the box $r_{0}$ the gravitational force balances the centrifugal force. Assuming an incompressible flow, one writes the fluid momentum equation and the induction equation for the magnetic field in the vicinity of $r_{0}$, a detailed discussion of the derivation is given in (e.g., Umurhan & Regev, 2004). The local orthogonal coordinate frame in the box is chosen such that the direction of $z$ coincides with the direction of $\Omega_{0}$, $x$ with the radial direction, and $y$ with the orbiting direction of the flow. The equilibrium background flow in this box then corresponds to the constant linear shear $\bm{v}_{0}(x)=-q\Omega_{0}x{\hat{y}}$. By representing the velocity field as $\bm{u}=\bm{v}+\bm{v}_{0}$, we formally eliminate the background shear. The shearing-box equations for the velocity and magnetic fields then take the form:222Because the mean field is zero in the setups we study here, we will refer only to the fluctuating part of the field $\bm{b}$.

[TABLE]

In this system, the magnetic field is measured in the Alfvénic units, $v_{A}=b/\sqrt{4\pi\rho_{0}}$, the density $\rho_{0}$ is constant, and we use the dimensionless variables: the time is normalized to $t_{0}=\Omega_{0}^{-1}$, the spatial variables to $L^{*}$, and the velocity to $\Omega_{0}L^{*}$, where $L^{*}\equiv\min(L_{x},L_{y},L_{z})$. The dimensionless pressure $P$ is not an independent field, but rather is defined as to ensure the incompressibility of the flow. In our numerical simulations we chose $\mathrm{Pm}=\nu/\eta=1$.

We will solve Equations (1)-(2) in a rectangular box $(L_{x},L_{y},L_{z})$, assuming periodic boundary conditions in the $z$- and $y$-directions, and shear-periodic boundary conditions in the $x$-direction Umurhan & Regev (2004). It is useful to discuss the conservation laws of Equations (1)-(2). In the non-rotating system, the MHD equations conserve the quadratic integrals of energy, cross helicity, and magnetic helicity if the external energy supply and energy dissipation are absent (e.g., Biskamp, 2003; Tobias et al., 2011). In the MRI case, energy is injected into the system by the instability. From Eqs. 1 and 2 one then derives for the energy Longaretti & Lesur (2010):

[TABLE]

where ${\tilde{\alpha}}\nobreak\hskip 1.66672pt=\nobreak\hskip 1.66672ptq\Omega_{0}\langle v_{x}v_{y}-b_{x}b_{y}\rangle$ is the energy injection rate, and the angular brackets denote an average over the box. The conservation law of cross helicity has to be modified in the shearing box as:

[TABLE]

and the conservation law for magnetic helicity remains unchanged:

[TABLE]

In these equations we introduce the vorticity $\bm{\omega}=\nabla{\times}{\bm{v}}$, the current density $\bm{j}=\nabla{\times}{\bm{b}}$, and the vector potential $\bm{A}$, where $\bm{b}=\nabla\times\bm{A}$. In the presence of the magnetorotational instability and in the absence of dissipation, the energy grows while the cross helicity and magnetic helicity do not.

The system (1)-(2) is solved using the pseudo-spectral, incompressible version of the code Snoopy (Lesur & Longaretti, 2007). Snoopy includes physical dissipation, while it uses a Fourier transform (FFTW 3 library) and a low-storage third-order Runge-Kutta (RK3) scheme for time evolution of all the fields. To employ the shear-periodic boundary conditions, Snoopy remaps the Fourier components of the fields periodically every $\Delta t_{\mathrm{remap}}=|L_{y}/(q\Omega_{0}L_{x})|$ (see the mathematical details of the procedure in  Umurhan & Regev (2004)).

In our work, we consider three simulation boxes with the following dimensions $(L_{x},L_{y},L_{z})\in\{(1,4,4),(1,4,16),(2,4,1)\}$. The numerical resolution in the first two cases is 128 points per unit length of $L_{x}$ and $L_{z}$, and 64 points per unit length of $L_{y}$. In the third case, there are 512 points per unit length of $L_{x}$ and $L_{z}$, and 256 points per unit length of $L_{y}$. Thus, the first and last cases have the same number of points in the vertical direction and are comparable after rescaling the unit of length. The first two cases have a viscosity of $\nu=1/5000$. The third case has a viscosity of $\nu=1/45000$. A summary of the cases studied is given in Table 1. Note that the energy and transport parameter reported are volume densities of these quantities. The transport parameter is defined in the standard way: $\alpha=\langle v_{x}v_{y}-b_{x}b_{y}\rangle/(q\Omega_{0}L_{z})^{2}$. The values shown in the table are averages over some time intervals in the steady states (cases I, II, and IIIb), and over a short time interval in the decaying case IIIa, where the energy is approximately constant. Cases I and II were initialized with $\bm{B}_{0}=(0,0,B_{0}\cos(2\pi x/L_{x}))$ and a large-scale, white-noise velocity configuration, and then evolved until a steady state is reached. Case III is the continuation of the zero-net-flux simulations discussed in our paper (Walker et al., 2016, case IV).

3 Results

Consistent with previous simulations (Nauman & Pessah, 2016), we find that our boxes with the aspect ratio $L_{z}/L_{x}\geq 4$ (cases I and II) are turbulent. They rapidly reach a steady state and maintain it for the extent of the simulation, while $E$ and $\alpha$ fluctuate around a mean value. Case III, which is a continuation of (Walker et al., 2016, case IV), exhibits very slowly decaying turbulence until $t\approx 600$, and appears to reach a steady state at later times (see Fig. 1).

In an attempt to explain these results, we find it useful to study the behavior of the azimuthal magnetic field $b_{y}$. As was noted in our previous work Walker et al. (2016), this field is concentrated at relatively large scales, and plays the role of the guide field for the remaining fluctuations. We found, however, an important difference in the distribution of this field in tall and long boxes. The net flux of the $b_{y}$ field is zero in the dynamo case. Given some initial perturbation of the magnetic field, the instability sets in, which leads to stretching and folding of the magnetic field lines in the $x-y$ direction. The shearing flow then increases the strength of the $b_{y}$ field. An important observation, however, is that the $x$-scale of the resulting folds of the $b_{y}$ field is always on the $x$-scale of the fluctuating velocity field. In the box extended in the $x$-direction (case III), this means that both the $b_{y}$ and the $v_{x}$ fields are still concentrated at comparable scales. This is seen in the spectra shown in Fig. 2.

The situation is qualitatively different in a box extended in the vertical direction (case II). It turns out that in this case the flux of $b_{y}$ can spread in the vertical direction over the scales much larger than the scales of the velocity field. This is seen in Fig. 3, where the $v_{z}$ component of the velocity field is concentrated at the scale comparable to the small horizontal box size $L_{x}$, while the $b_{y}$ magnetic field is concentrated at the larger scale $L_{z}$. The structure of the $b_{y}$ field is shown in Fig. 4, which also reveals separation of the regions of positive and negative azimuthal magnetic flux in the vertical direction. This is in contrast to Fig. 5 where $b_{y}$ is more homogeneous, due to large-scale turbulent mixing. Fig. 6 shows case I, an intermediate between cases II and III.

In order to quantify this flux intermittency and to compare it with the other cases, we subdivided $y=\rm{const}$ cross-sections of the simulation boxes into small $1/16\times 1/16$ squares and calculated azimuthal magnetic fluxes for each of these squares. We then plotted histograms of the obtained fluxes $\Phi_{y}$. The results are shown in Fig. 7. The flux intermittency is smallest in case IIIa, larger in case I, and the largest in the tallest box of case II. The increased flux intermittency implies more favorable conditions for the dynamo action.

We are now in a position to answer why our case III did not exhibit dynamo action until the turbulence decayed to very small amplitudes (see Fig. 1). The matter is that when the energy of fluctuations declines, so does the scale of the turbulence (see Fig. 5). In our work Walker et al. (2016) this is explained by the necessity to maintain the balance between the linear and nonlinear shearing rates of the turbulence. When the $x$-scale of the turbulence was relatively large compared to the box size $L_{z}$, the dynamo did not operate. When the $x$-scale of the turbulence decreased as to become smaller than the $L_{z}$ size of the box, in analogy with cases II and I the dynamo action became possible. Consistent with our discussion above, the intermittency of the magnetic flux also increased in this regime (see Fig. 7).

Although the dynamo action seems to operate in case IIIb, we should exercise caution in claiming that a true steady state is observed. As pointed out in (Rempel et al., 2010) it may be possible in such a case that the system is described by a supertransient state and may ultimately decay if integrated long enough 333Although, with $\mathrm{Rm}=45000$, the expected lifetime is $\mathcal{O}(10^{32}$).. Also, unlike the other dynamo cases in this study, IIIb has developed, in addition to the small-scale fluctuations, a large-scale $v_{y}(x)$ zonal flow resulting from geostrophic balance (Johansen et al., 2009) (see Fig. 8). Comparing the spectra of case III from early times in Fig. 2 with those from late times in Fig. 9, we see that this zonal flow is long-lived, or at least steadily reinforced. The possibility of such a flow may be related to the long radial extent of the box. This flow however is practically uncorrelated with the small-scale fluctuations; for instance, it does not contribute to the transport coefficient $\alpha_{v}$ in Fig. (1).

4 Conclusions

We have found in incompressible MHD simulations of the shearing box that the properties of the zero-net-flux dynamo action depend intimately on the relation between the scale of the turbulence and the vertical size of the box. For a given set of physical parameters, transport is sustained and a dynamo action rapidly sets in in systems with high $L_{z}/L_{x}$ aspect ratio, while in the opposite case of small aspect ratio, the dynamo action is not observed until the scale of the turbulence significantly decreases. This is the despite the size of the box being much greater than the dissipative scale of turbulence in both cases. Based on our results we suggest the following explanation for this phenomenon. We propose that an inherent property of the shearing-box turbulence is its tendency to build up local regions of nonzero $b_{y}$ flux, separated in the vertical direction $z$. This may be related to the fact that the configurations ${\bm{b}}=(0,b_{y}(z),0)$, ${\bm{v}}=(v_{x}(z,t),v_{y}(z,t),0)$ are exact solutions of the ideal equations (e.g., Goodman & Xu, 1994). Those sub-regions with nonzero $b_{y}$ fluxes are MRI-unstable. The resulting turbulent eddies, on the other hand, are trying to mix and homogenize the $b_{y}$ field. The size of the turbulent eddies is always comparable to the size of the folds of the $b_{y}$ field in the $x$-direction. Indeed, such eddies are caused by the MRI instability of the $b_{y}$ field. In “long” boxes, the eddies caused by an initially unstable large-scale perturbation have a typical length-scale larger than the vertical extent of the box, and they are able to effectively mix the field in that direction, thus reducing the $b_{y}$ fluxes and effectively reducing the rate of the instability as compared to the rate of turbulent energy dissipation. The resulting turbulence will therefore decay, until the scale of the turbulence becomes smaller than the vertical size of the box, at which point the vertical flux separation enhances the efficiency of the instability. In “tall” boxes, on the other hand, the eddies are always smaller than the vertical extension of the box. This allows oppositely-directed $b_{y}$ fluxes to build up in vertically separated regions, allowing the instability to win.

This poses a broader question to what extent the shearing-box results describe the natural MRI turbulence. The shearing-box model suggests that the scale of the MRI-dynamo-driven turbulence should be smaller than the thickness of the disc and that the resulting transport coefficients are very small. We, however, notice that the vertical limitations in real accretion disks are never as severe as the periodic boundary conditions of the shearing-box simulations. Allowing the flux of the $b_{y}$ magnetic field to partially escape the box in the vertical direction, for instance, will facilitate the dynamo action rather than impede it. This indicates that the shearing-box model with periodic boundary conditions may not allow one, at least at the present level of understanding, to realistically simulate the scale and the transport properties of the MRI-dynamo-driven turbulence at low magnetic Prandtl numbers. What seems to be a robust feature of the MRI-dynamo-driven turbulence, however, is its tendency to separate the regions of positive and negative $b_{y}$ fluxes, leading to highly intermittent azimuthal magnetic flux patches. In a stratified disk, these magnetic flux tubes can then buoy above the surface, forming an intermittently magnetized corona (e.g. Galeev et al., 1979). Non-periodic large-scale conditions (e.g., stratification, flux escape, etc.) need to be incorporated in the shearing-box model in order to provide a more realistic description.

Acknowledgments

We thank Fausto Cattaneo and Geoffroy Lesur for useful comments. JW is supported by the Wisconsin Alumni Research Foundation at the University of Wisconsin-Madison. SB is partly supported by the National Science Foundation under the grant NSF AGS-1261659 and by the Vilas Associates Award from the University of Wisconsin-Madison. Simulations were performed at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin under the NSF-Teragrid Project TG-PHY110016.