Shearing-box simulations of MRI-driven turbulence in weakly collisional accretion discs

TL;DR

This study investigates MRI-driven turbulence in weakly collisional accretion disks using shearing-box simulations with anisotropic viscosity, revealing that turbulence levels are similar to MHD results but flow structures and heating are significantly affected.

Contribution

It is the first systematic shearing-box study including anisotropic (Braginskii) viscosity in MRI turbulence, exploring a broad parameter space and its effects on turbulence and flow structure.

Findings

Turbulence levels are similar to ideal MHD despite large anisotropic viscosities.

Anisotropic viscosity alters flow structure and redistributes turbulence.

At high Reynolds numbers, anisotropic transport becomes negligible, but dominates plasma heating.

Abstract

We present a systematic shearing-box investigation of MRI-driven turbulence in a weakly collisional plasma by including the effects of an anisotropic pressure stress, i.e. anisotropic (Braginskii) viscosity. We constrain the pressure anisotropy ($\Delta p$) to lie within the stability bounds that would be otherwise imposed by kinetic microinstabilities. We explore a broad region of parameter space by considering different Reynolds numbers and magnetic-field configurations, including net vertical flux, net toroidal-vertical flux and zero net flux. Remarkably, we find that the level of turbulence and angular-momentum transport are not greatly affected by large anisotropic viscosities: the Maxwell and Reynolds stresses do not differ much from the MHD result. Angular-momentum transport in Braginskii MHD still depends strongly on isotropic dissipation, e.g., the isotropic magnetic Prandtl number, even when the anisotropic viscosity is orders of magnitude larger than the isotropic diffusivities. Braginskii viscosity nevertheless changes the flow structure, rearranging the turbulence to largely counter the parallel rate of strain from the background shear. We also show that the volume-averaged pressure anisotropy and anisotropic viscous transport decrease with increasing isotropic Reynolds number (${\rm Re}$); e.g., in simulations with net vertical field, the ratio of anisotropic to Maxwell stress ($\alpha_{\rm A} / \alpha_{\rm M}$) decreases from $\sim 0.5$ to $\sim 0.1$ as we move from ${\rm Re} \sim 10^3$ to ${\rm Re} \sim 10^4$. Anisotropic transport may thus become negligible at high ${\rm Re}$. Anisotropic viscosity nevertheless becomes the dominant source of heating at large ${\rm Re}$, accounting for $\gtrsim 50 \%$ of the plasma heating. We conclude by briefly discussing the implications of our results for RIAFs onto black holes.

Tables

Table 1. Table 2: Summary of Braginskii MHD simulations with net vertical magnetic field and fourth-order isotropic hyperdiffusion. Full simulations were evolved for a time t ​ Ω = 100 𝑡 Ω 100 t\Omega=100 and the values of 4 ​ π ​ ⟨ Δ ​ p ⟩ / ⟨ B 2 ⟩ 4 𝜋 delimited-⟨⟩ Δ 𝑝 delimited-⟨⟩ superscript 𝐵 2 4\pi\langle\Delta p\rangle/\langle B^{2}\rangle , ⟨ 4 ​ π ​ Δ ​ p / B 2 ⟩ delimited-⟨⟩ 4 𝜋 Δ 𝑝 superscript 𝐵 2 \langle 4\pi\Delta p/B^{2}\rangle and /are temporal averages over t ​ Ω = 70 − 100 𝑡 Ω 70 100 t\Omega=70-100 . In Composite simulations, MHD fields are restarted at t ​ Ω = 100 𝑡 Ω 100 t\Omega=100 using Braginskii MHD. 4 ​ π ​ ⟨ Δ ​ p ⟩ / ⟨ B 2 ⟩ 4 𝜋 delimited-⟨⟩ Δ 𝑝 delimited-⟨⟩ superscript 𝐵 2 4\pi\langle\Delta p\rangle/\langle B^{2}\rangle , ⟨ 4 ​ π ​ Δ ​ p / B 2 ⟩ delimited-⟨⟩ 4 𝜋 Δ 𝑝 superscript 𝐵 2 \langle 4\pi\Delta p/B^{2}\rangle and /were averaged over t ​ Ω = 101 − 110 𝑡 Ω 101 110 t\Omega=101-110 , except for the LABEL:= ​ 7.5 × 10 8 LABEL:= 7.5 superscript 10 8 \ref{=}7.5\times 10^{8} simulation, for which the average is over t ​ Ω = 101 − 102 𝑡 Ω 101 102 t\Omega=101-102 .

Resolution	Sim. Type		${\mathrm{Re}}_{M,4}$	${\nu }_4/{\eta }_4$	${\mathrm{Re}}_B$	$4\pi ⟨\Delta p⟩/⟨B^2⟩$	$⟨4\pi \Delta p/B^2⟩$	/
$(256,128,64)$	Full	$7.5\times {10}^5$	$7.5\times {10}^5$	1	0.75	0.30	0.20	0.38
$(256,128,64)$	Composite	$7.5\times {10}^6$	$7.5\times {10}^6$	1	0.75	0.19	0.13	0.26
$(256,128,64)$	Full	$3\times {10}^7$	$3\times {10}^7$	1	0.75	0.071	0.026	0.13
$(384,192,96)$	Composite	$1.5\times {10}^8$	$1.5\times {10}^8$	1	0.75	0.069	0.012	0.13
$(576,288,144)$	Composite	$7.5\times {10}^8$	$7.5\times {10}^8$	1	0.75	0.034	-0.0096	0.085

Table 2. Table 3: Summary of simulations with net toroidal and vertical magnetic field ( ⟨ B y ⟩ = ⟨ B z ⟩ delimited-⟨⟩ subscript 𝐵 𝑦 delimited-⟨⟩ subscript 𝐵 𝑧 \langle B_{y}\rangle=\langle B_{z}\rangle ). The transport coefficients and mean pressure anisotropies were averaged over t ​ Ω = 100 − 200 𝑡 Ω 100 200 t\Omega=100-200 .

Resolution	$\mathrm{Re}$	${\mathrm{Re}}_M$	$\mathrm{Pm}$	${\mathrm{Re}}_B$	${\alpha }_{\mathrm{Re}}$	${\alpha }_M$	${\alpha }_A$	$\alpha$	$4\pi ⟨\Delta p⟩/⟨B^2⟩$	$⟨4\pi \Delta p/B^2⟩$	/
$(256,128,64)$	$3000$	1500	0.5	–	0.0056	0.022	–	0.028	–	–	–
$(256,128,64)$	$3000$	1500	0.5	0.75	0.0064	0.019	0.0076	0.033	0.32	0.24	0.40
$(256,128,64)$	$750$	750	1	–	0.0059	0.023	–	0.029	–	–	–
$(256,128,64)$	$750$	750	1	0.75	0.0064	0.020	0.0094	0.036	0.41	0.30	0.47
$(256,128,64)$	$4500$	4500	1	–	0.0067	0.033	–	0.040	–	–	–
$(256,128,64)$	$4500$	4500	1	0.75	0.011	0.042	0.011	0.063	0.18	0.11	0.25
$(192,96,48)$	$750$	3000	4	–	0.011	0.067	–	0.078	–	–	–
$(192,96,48)$	$750$	3000	4	0.3	0.014	0.063	0.022	0.099	0.28	0.19	0.36

Equations

$$\nabla \cdot \mathbf{U}=0,$$

$$\begin{array}{rl}\frac{D\mathbf{U}}{Dt}=& -2\mathbf{\Omega }\times \mathbf{U}+2\Omega Sx\hat{\mathbf{x}}-\nabla \left(p_{\perp }+\frac{B^2}{8\pi }\right)\\ & +\nabla \cdot \left[\hat{\mathbf{b}}\hat{\mathbf{b}}\left(\frac{B^2}{4\pi }+\Delta p\right)\right]+\nu {\nabla }^2\mathbf{U},\end{array}$$

$$\frac{D\mathbf{B}}{Dt}=\mathbf{B}\cdot \nabla \mathbf{U}+\eta {\nabla }^2\mathbf{B},$$

$$\Delta p=p_{\perp }-p_{\parallel }$$

$$\text{centering}\Delta p\gtrsim \frac{B^2}{8\pi }\text{@add@centering}$$

$$\Delta p\lesssim -\frac{B^2}{4\pi }.$$

$$\mathrm{Re}=\frac{S{L}_z^2}{\nu };$$

$$\mathrm{R}{\mathrm{e}}_{\mathrm{M}}=\frac{S{L}_z^2}{\eta };$$

$$\mathrm{Pm}=\frac{\nu }{\eta }.$$

$$\mathrm{R}{\mathrm{e}}_{\mathrm{B}}=\frac{S{L}_z^2}{{\mu }_{\mathrm{B}}}.$$

$$\alpha ={\alpha }_{\mathrm{Re}}+{\alpha }_{\mathrm{M}}+{\alpha }_{\mathrm{A}},$$

$$\alpha_{\rm Re}=\langle v_{x}v_{y}\rangle\ /\ \big{(}S^{2}L_{z}^{2}\big{)},$$

$$\alpha_{\rm M}=-\Big{\langle}\frac{B_{x}B_{y}}{4\pi}\Big{\rangle}\ /\ \big{(}S^{2}L_{z}^{2}\big{)},$$

$$\alpha_{\rm A}=-\Big{\langle}\frac{\Delta p}{B^{2}}B_{x}B_{y}\Big{\rangle}\ /\ \big{(}S^{2}L_{z}^{2}\big{)}$$

$$⟨\mathbf{B}⟩=B_0\hat{\mathbf{z}}.$$

$$D=D_{\eta }+D_{\nu }+D_{{\mu }_{\mathrm{B}}},$$

$$D_{\eta }=-\frac{\eta }{4\pi }\int d^3\mathbf{x}\mathbf{B}\cdot {\nabla }^2\mathbf{B},$$

$$D_{\nu }=-\nu \int d^3\mathbf{x}\mathbf{u}\cdot {\nabla }^2\mathbf{u},$$

$$D_{{\mu }_{\mathrm{B}}}=\int d^3\mathbf{x}\Delta p\hat{\mathbf{b}}\hat{\mathbf{b}}:\nabla \mathbf{U}.$$

$$\text{ref}=S{L}_z^4/{\nu }_4.$$

$$⟨\mathbf{B}⟩=B_0(\hat{\mathbf{y}}+\hat{\mathbf{z}}),$$

$$\bm{B}(t=0)=B_{0}\sin\Big{(}\frac{2\pi x}{L_{x}}\Big{)}\hat{\bm{z}}.$$

Full text

††pubyear: 2019††pagerange: Shearing-box simulations of MRI-driven turbulence in weakly collisional accretion discs–B

1 Introduction

Magnetohydrodynamic (MHD) turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1991) is widely considered to be one of the key engines powering angular-momentum transport in accretion discs. As a result, the growth of the MRI and the subsequent MHD turbulence it produces have been studied extensively over the years (see e.g. Hawley et al. 1995; Balbus & Hawley 1998; Hawley et al. 2001).

One uncertainty in the application of the MRI is that the MHD fluid approximation is not well justified in a number of astrophysical systems. This includes radiatively inefficient accretion flows (RIAFs), in which the Coulomb mean free path is larger than the typical system size (Mahadevan & Quataert 1997). Departures from the ideal-MHD framework are therefore required and at first glance it may seem necessary to model the system as a fully collisionless plasma in six-dimensional phase space. The goal of this paper is to better understand the nonlinear evolution of the MRI under such conditions.

While it has recently become possible to run kinetic simulations of the MRI using particle-in-cell codes (Riquelme et al. 2015; Hoshino 2015; Kunz et al. 2016; Inchingolo et al. 2018), such simulations remain far too expensive (at least in three dimensions) to explore parameter space. However, a variety of recent theory (Schekochihin et al. 2008; Kunz et al. 2014; Sironi & Narayan 2015; Riquelme et al. 2015) suggests that ion-Larmor scale kinetic instabilities such as the mirror and firehose instabilities, which grow readily in low collisionality, weakly magnetized plasmas, act to increase the effective collision rate via wave-particle interactions. This result is of great utility, as it at least partially motivates modeling the system as a weakly collisional plasma. The advantage of this framework is that non-ideal effects are simply introduced as additional terms in the ideal fluid equations, which is much simpler than evolving a six-dimensional distribution function of the plasma particles.

The linear growth of the MRI differs from its ideal-MHD counterpart in both a collisionless (Quataert et al. 2002) and a weakly collisional (Balbus 2004) plasma. However, despite some work with simplified fluid models and kinetic simulations, we lack detailed understanding of how kinetic physics affects the saturated MRI turbulence. Some insight has been gained by the work of Squire et al. (2017a), who focused on the nonlinear growth phase of the MRI in high-$\beta$ low-collisionality plasmas. They argued that due to the onset of the aforementioned microinstabilities, the nonlinear growth phase of the kinetic MRI (KMRI) always returns to MHD-like evolution. Similarly, the saturation into turbulence appeared to be unaffected by non-ideal physics. These results provided insight into the physics behind earlier work on collisionless accretion discs by Sharma et al. (2006), who found that the properties of KMRI-induced turbulence were not too different from MHD. This resemblance has also been found in global general-relativistic simulations (Foucart et al. 2016; Foucart et al. 2017), which employed an extended-MHD framework with anisotropic viscosity and conductivity.

In this paper, we carry out a systematic shearing-box study of MRI-induced turbulence in a low-collisionality plasma with explicit resistivity and viscosity. To model non-ideal effects we use Braginskii’s closure for magnetized, weakly collisional plasmas (Braginskii 1965), commonly referred to as “Braginskii MHD”. As explained below, the anisotropic viscosity in Braginskii’s closure is equivalent to including an anisotropic pressure stress in the MHD equations. The closure is the simplest, well motivated model to capture key aspects of kinetic physics on large scales.

There are a number of questions that motivate such a parameter exploration. Perhaps most importantly, it will clarify the relevance of non-ideal physics for angular-momentum transport and plasma heating, including the additional contribution to the total stress tensor that comes directly from the pressure anisotropy.

While the Maxwell and Reynolds stresses have been studied extensively in MHD, significantly less is known about the non-ideal, anisotropic viscous stress. Most simulations to date found that its contribution to angular-momentum transport is smaller than, but comparable to, the Maxwell stress. However, given that the pressure anisotropy is driven by gradients in the velocity field, we may expect isotropic dissipation to influence the anisotropic transport. It is therefore instructive to look at the relationship between anisotropic stress and the dimensionless isotropic Reynolds numbers.

Exploring a range of isotropic viscosities and resistivities is vital for a second reason. One of the most striking results of previous work on MRI-generated turbulence is the dependence on isotropic dissipation. Lesur & Longaretti (2007) and Fromang et al. (2007) showed that the MRI saturation amplitude is very sensitive to the choice of viscosity and resistivity, with a particularly strong dependence on their ratio, the magnetic Prandtl number. It is plausible to speculate that an additional large anisotropic viscosity may alter the effective Prandtl number. This claim is further motivated by the kinetic MRI simulations of Kunz et al. (2016), who showed that a high-$\beta$ collisionless plasma behaved in some ways like a high-magnetic-Prandtl-number fluid. It is therefore unclear to what extent we should expect to recover the usual Prandtl-number dependence in Braginskii MHD.

Another important factor to consider is how the very building blocks of MHD turbulence are modified in low-collisionality plasmas and whether this may change the large-scale turbulent state in low-collisionality accretion discs. Squire et al. (2016) and Squire et al. (2017b) showed that collisionless and weakly collisional plasmas cannot support linearly polarized shear-Alfvén waves above a critical amplitude due to a cancellation between the Lorentz force and the anisotropic-pressure force. It is unclear if and how this might affect the large-scale turbulent properties in accretion discs.

The layout of this paper is as follows. In Section 2 we discuss the method and setup for our study of turbulence in Braginskii MHD. The main focus is on boxes threaded by a net vertical magnetic field, with the corresponding results described in Section 3. We consider other initial field configurations in Section 4. Finally, Section 5 summarizes our key results and discusses current limitations and future directions.

2 Method

2.1 Equations

We use the pseudo-spectral code SNOOPY (Lesur & Longaretti 2007) to evolve the incompressible MHD equations with anisotropic pressure in a shearing box:

[TABLE]

where $\bm{\hat{b}}=\bm{B}/B$ is the unit vector along the magnetic field $\bm{B}$ and the density has been set to unity. Because our model is incompressible, we choose $p_{\perp}$ at each timestep so as to satisfy $\bm{\nabla\cdot U}=0$. We include an explicit isotropic viscosity $\nu$ and resistivity $\eta$. The velocity field $\bm{U}$ consists of a background shear $\bm{U_{0}}$ and perturbations $\bm{u}$: $\bm{U}=\bm{U_{0}}+\bm{u}$. We adopt an equilibrium Keplerian background profile $\bm{U_{0}}=-Sx\hat{\bm{y}}$, with $S=\frac{3}{2}\Omega$, and we use the code to compute the evolution of $\bm{u}$ and $\bm{B}$. We use the $2/3$ de-aliasing rule to prevent spurious modes originating from the nonlinear terms in (1)–(3).

At any timestep, the pressure anisotropy $\Delta p$ entering the momentum equation (LABEL:eq:momentum) is calculated via:

[TABLE]

where $p_{\perp}$ and $p_{\parallel}$ are the thermal pressures in the directions perpendicular and parallel to the local magnetic field (Braginskii 1965). Equation (4) can be obtained from the kinetic evolution equations of the plasma (Chew et al. 1956; Kulsrud 1983; Schekochihin et al. 2010) in the weakly collisional regime $\nu_{c}/|\bm{\nabla u}|\gg 1$, where $\nu_{c}$ is the collision rate of the plasma. In equation (4) we also assume that the effect of heat fluxes on the pressure anisotropy can be neglected, an assumption that is formally valid when $\nu_{c}/|\bm{\nabla u}|\gg\beta^{1/2}$ (Mikhailovskii & Tsypin 1971), where $\beta=8\pi p/B^{2}$ is the ratio of thermal to magnetic pressure (see Squire et al. 2017a for more discussion of the different regimes).

Upon substitution into (LABEL:eq:momentum), $\Delta p$ can be shown to behave as a diffusion operator acting along $\bm{\hat{b}}$. Its role is therefore to damp the component of the velocity along the local magnetic field that has gradients along the local magnetic field. Due to its diffusive contribution, throughout this work we will refer to the coefficient $\mu_{\rm B}$ as anisotropic (or Braginskii) viscosity.

2.2 Modeling Kinetic Microinstabilities

The background shear and MRI-induced magnetic-field growth naturally lead to a finite pressure anisotropy. However, once $\Delta p$ becomes comparable to the magnetic pressure, plasma instabilities are excited that are not fully captured by the set of equations (1)–(3). The two main microinstabilities that need to be accounted for are the mirror instability (Barnes 1966; Hasegawa 1969), excited if

[TABLE]

and the firehose instability (Rosenbluth 1956; Chandrasekhar et al. 1958; Parker 1958), which is excited when

[TABLE]

Wave-particle interactions induced by these instabilities increase the effective collisionality of the system, which in turn acts to isotropize the pressure tensor. As a result, the mirror/firehose instability, excited by a growing/declining pressure anisotropy, acts to halt further growth/decline. Kinetic simulations have shown that these instabilities tend to pin the anisotropy near the instability thresholds (Kunz et al. 2016).

We include this kinetic result by imposing hard-wall limits on $\Delta p$. If $\Delta p$ is driven outside of the stability limits given by $-B^{2}/4\pi<\Delta p<B^{2}/8\pi$, then it is pinned to either $\Delta p=-B^{2}/4\pi$ or $\Delta p=B^{2}/8\pi$, depending on which boundary is crossed. Otherwise, $\Delta p$ is determined by equation $\eqref{eq:delp}$.

By using instantaneous bounds on $\Delta p$, we essentially assume that the only effect of microinstabilities is to halt the growth of $\Delta p$, with no direct change to other fluid quantities. For the firehose instability, this limiting behaviour arises due to particle scattering, while for the mirror instability, there is a long phase where small-scale magnetic fluctuations grow secularly in time. Although the limiter model can, in principle, capture the effect of either scattering or secular growth if we consider $\bm{u}$ and $\bm{B}$ in equations (1)–(3) to be large-scale averages, there may be other poorly understood effects that are not captured. In addition, the assumption that $\Delta p$ limiters act instantaneously is incorrect for motions with timescales approaching the ion gyro time. Because we find that the cascade continues almost unaffected to scales well below the Braginskii viscous scale, this could mean that the smallest-scale motions will be more strongly affected by $\Delta p$ forces than we assume here. It is unclear if and how these additional effects could be incorporated in a fluid model, and a kinetic description is likely necessary to fully capture all the underlying physics (see Schekochihin et al. 2008, Kunz et al. 2014, Squire et al. 2017a for more discussion).

We have also run two simulations without limiters at the firehose instability threshold. These no-firehose-limiter simulations are partially justified by the fact that the parallel firehose instability is already present in the Braginskii equations (by contrast, the mirror instability is not). However, the kinetic simulations in Kunz et al. (2014) showed that the oblique firehose instability is probably more important for maintaining $\Delta p$ near the instability threshold. For this reason, most of our simulations have both firehose and mirror limiters included.

2.3 Setup

Throughout this work we set $\Omega=1$. Since we want to explore any pressure-anisotropy-induced differences, for every Braginskii MHD simulation we also carry out a corresponding MHD simulation in which $\Delta p=0$. Given the sensitivity of MRI turbulence to isotropic dissipation, we test a number isotropic viscosities $\nu$ and resistivities $\eta$. We define the associated dimensionless Reynolds number,

[TABLE]

the magnetic Reynolds number,

[TABLE]

and their ratio, the magnetic Prandtl number,

[TABLE]

We also define the analogous Braginskii Reynolds number,

[TABLE]

To quantify the turbulent angular-momentum transport, we define the dimensionless transport coefficient

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

are the contributions from the volume-averaged ($\langle...\rangle$) Reynolds stress, Maxwell stress and anisotropic viscous stress respectively, normalized by $S^{2}L_{z}^{2}$. This is the incompressible version of the compressible transport parameter, which is usually normalized using the initial pressure (see e.g. Hawley et al. 1995; Sharma et al. 2006). While $\alpha_{\rm M}$ and $\alpha_{\rm Re}$ are present both in ordinary MHD and in Braginskii MHD, $\alpha_{\rm A}$ requires a pressure anisotropy and is therefore only nonzero in Braginskii MHD.

To capture the most important parasitic modes that break up the MRI “channel” modes into turbulence, most of our simulations are in horizontally elongated boxes of size $L_{x}=4$, $L_{y}=4$ and $L_{z}=1$ (Bodo et al. 2008; Pessah & Goodman 2009; Longaretti & Lesur 2010). However, we did also explore other aspect ratios. We find that Braginskii MHD results are particularly sensitive to box size, with dramatically different results in horizontally narrow boxes with net vertical flux (see Appendix A.1).

In order to satisfy the Courant condition at large $\mu_{\rm B}$, we are limited to rather modest resolutions by current standards, despite sub-cycling over the $\Delta p$ term in equation (LABEL:eq:momentum) in simulations with large $\mu_{\rm B}$ (where we used 5 or 8 as the maximum number of sub-cycles per main MHD timestep). Most of our full simulations in $4\times 4\times 1$ boxes have resolution $256\times 128\times 64$. To test very large $\mu_{\rm B}$, we also ran a number of lower resolution simulations. In each case, the initial magnetic field and background shear are perturbed with small-amplitude white noise.

The isotropic viscosities and resistivities are chosen such that the dissipative scales are properly resolved. This places an upper bound on the $\mathrm{Re}$ and $\mathrm{Re_{M}}$ that we can explore in a full Braginskii simulation, given the attainable resolutions. To explore the relationship between anisotropic stress and the isotropic diffusivities, ideally we would like to cover a broad range of viscosities and resistivites, and explore the limit $\mathrm{Re}\rightarrow\infty$, $\mathrm{Re_{M}}\rightarrow\infty$ at fixed $\mathrm{Pm}$. However, the required resolutions are computationally unfeasible with the numerical methods that we use here.

To test larger Reynolds numbers in higher resolution, we perform a number of “Composite” MHD–Braginskii MHD simulations. In these simulations, we first evolve the equations of MHD for a time $t\Omega=100$. The turbulent MHD flow fields are then restarted with anisotropic viscosity included and evolved further for several $\Omega^{-1}$. These composite simulations are motivated by our observation that anisotropic viscosity transforms MHD flow fields into Braginskii-like flow fields on timescales shorter than the orbital time. We have tested that, by restarting MHD turbulence with anisotropic viscosity, we are able to recover the typical $\Delta p$ and $\alpha_{\rm A}$ of the corresponding full Braginskii simulation in a fraction of $\Omega^{-1}$. We show an example of this behavior in Appendix B.

This method offers insight into the statistics of Braginskii MHD turbulence at large Reynolds numbers, even when the Braginskii equations are evolved for a rather short amount of time. As a result, the composite simulations enable us to explore isotropic Reynolds numbers ($\mathrm{Re}\sim 10^{4}$) and resolutions ($768\times 384\times 192$ and $384\times 192\times 96$) that are otherwise unattainable for a full Braginskii simulation with large $\mu_{\rm B}$.

We also perform a number of simulations in which the viscosity and resistivity are replaced by hyperdiffusion operators, $\nu_{4}\nabla^{4}\bm{U}$ and $\eta_{4}\nabla^{4}\bm{B}$. We define the associated dimensionless quantities $\ref{=}SL_{z}^{4}/\nu_{4}$ and $\mathrm{Re_{M,4}}=SL_{z}^{4}/\eta_{4}$. Using hyperdiffusion serves as an alternative method to probe larger effective Reynolds numbers, by potentially increasing the size of the inertial range without increasing the resolution. The $k^{4}$ dependence of hyperdiffusion allows us to dissipate energy above the grid scale even for very small $\nu_{4}$, without constraining the turbulence at intermediate wavenumbers, thus mimicking turbulence at somewhat larger $\mathrm{Re}$ than would be possible with standard diffusion operators (note, however, that the \ref{}definedaboveis\textit{not}theeffectiveReynoldsnumberofourhyperdiffusionsimulations).\par Whileourmainfocusisonsimulationdomainswithanetverticalfield(Section\ref{sec:vertical}),wealsodiscussotherinitialmagnetic-fieldconfigurations:netverticalandtoroidalfield(\ref{sec:toroid}),aswellaszeronetflux(\ref{sec:zero}).\par\par\par$

3 Net Vertical Field

The main part of our study concerns boxes initially threaded by a purely vertical magnetic field,

[TABLE]

We choose $B_{0}=\sqrt{8\pi/1348}\ \Omega L_{z}$, so that the fastest-growing MRI mode in MHD has wavelength $\lambda_{\rm MRI}=0.25L_{z}$.

Table 1 gives a summary of our different choices of $\mathrm{Re}$, $\mathrm{Re_{M}}$ and $\mathrm{Re_{B}}$. The main result of this section concerns anisotropic stress and its dependence on the values of $\mathrm{Re}$ and $\mathrm{Re_{M}}$. However, we postpone our discussion of this result until section 3.2 and first look at the overall Braginskii MHD evolution, and its similarities and differences relative to MHD.

3.1 Comparison to MHD

In each of our simulations the qualitative evolution follows the same pattern. The initial small-amplitude perturbations are amplified by the MRI. The amplification continues until the MRI modes reach large amplitudes and become unstable to parasitic instabilities (Goodman & Xu 1994). These are secondary instabilities that grow on the large field and flow gradients in the MRI mode, causing it to break up into turbulence. Squire et al. (2017a) argued that the dominant parasitic modes are not too different in Braginskii MHD. This is broadly consistent with our numerical results summarized below, although we do often see somewhat larger $\alpha$ in the Braginskii case during a short initial transient phase. In addition, we find a very significant box-size dependence in Braginskii MHD (see Appendix A.1).

Figure 1 shows the evolution of the simulation with $\mathrm{Re}=1500$ and $\mathrm{Pm}=2$. The dotted lines track the MHD evolution, while the solid lines correspond to Braginskii MHD with $\mathrm{Re_{B}}=0.75$. The two models show similar behavior and the resultant turbulent energy densities and angular-momentum transport are close to identical. The main qualitative differences appear at early times, during the MRI growth phase. Braginskii viscosity delays the growth along the initial field direction ($\langle B_{z}^{2}\rangle$ and $\langle v_{z}^{2}\rangle$) through stronger damping of the initial white-noise perturbations. In addition, in the Braginskii MHD simulation the MRI is able to grow to larger amplitudes, before it eventually saturates.

For all of the cases we simulated, the transport is not greatly affected by the anisotropic viscosity. This includes our simulations with both limiters included, as well as our full simulation without a firehose limiter. The main difference is the presence of the anisotropic viscous stress. The Maxwell and Reynolds stresses are very similar in both MHD and Braginskii MHD.

The remarkable similarities between the MHD and weakly collisional solutions draw us to an interesting conclusion: even with large anisotropic viscosity, the turbulent amplitudes are still set primarily by the isotropic Reynolds and magnetic Reynolds numbers. This is illustrated in Figure  2, where we show the time-averaged transport coefficients of the MHD (black) and Braginskii MHD (red) runs at different $\mathrm{Pm}$. For the Braginskii runs we also plot the transport due to just the Reynolds and Maxwell stresses as empty red diamonds. The error bars shown in Figure 2 are estimates for the standard deviations of the average transport coefficients, which use a binning time of $t=10\Omega^{-1}$ (due to the short Braginskii timesteps, our simulations were not evolved long enough for an error analysis similar to Longaretti & Lesur 2010). Nevertheless, they illustrate that the $\alpha_{\rm Re}+\alpha_{\rm M}$ of the Braginskii calculation agrees well with MHD. In Braginskii MHD we recover the usual Prandtl-number dependence found in MHD, thus showing that turbulence is still strongly influenced by the isotropic diffusivities, even when $\mathrm{Re_{B}}\ll\mathrm{Re}$.

The result that angular-momentum transport is not significantly affected by large anisotropic viscosities is partly due to the anisotropic-pressure limiters. These limit $\Delta p$ to being comparable to the local field strength, which implies that the anisotropic stress cannot become significantly larger than the Maxwell stress. However, we see here that the box-averaged anisotropic stress can be substantially smaller than the Maxwell stress (e.g., $\alpha_{\rm A}=0.25\alpha_{\rm M}$ in the full Braginskii run with $\mathrm{Re}=4500$ and $\mathrm{Pm}=1$; we also find that $\alpha_{\rm A}\ll\alpha_{\rm M}$ in high-$\mathrm{Re}$ composite simulations, as discussed in Section 3.2), which is less obvious, and surprising in light of previous results (e.g., Sharma et al. 2006 and Kunz et al. 2016, where $\alpha_{\rm A}$ is comparable to $\alpha_{\rm M}$). Perhaps even more surprising is that we find little dependence of the angular-momentum transport on anisotropic viscosity even though the effective ($\Delta p$-limited) Braginskii viscosities considered here are much larger than the isotropic diffusivities. This is in contrast to the strong dependence of angular-momentum transport in a shearing box on isotropic diffusivities (Fromang et al. 2007; Lesur & Longaretti 2007). In addition, we show in Section 3.2 that angular-momentum transport is similar to MHD despite the fact that the flow structure in Braginskii MHD is quite different from MHD (e.g., Figure 9).

3.2 Anisotropic Transport in Braginskii MHD

In this section we explore aspects of angular-momentum transport that are specific to Braginskii MHD. Having demonstrated that the Maxwell and Reynolds stresses tend to track their values in the complementary MHD simulations, our main focus is on the evolution of the anisotropic stress and pressure anisotropy.

In Figure 3a we show the evolution of the Maxwell, Reynolds and anisotropic viscous stresses for the simulation with $\mathrm{Re}=1500$, $\mathrm{Pm}=2$ and $\mathrm{Re_{B}}=0.75$. The overall transport is dominated by the Maxwell stress, followed by the anisotropic and Reynolds stresses. This is similar to the results of the kinetic MRI simulations in Kunz et al. (2016) and the global extended-MHD simulation in Foucart et al. (2017). In all of our simulations we find that the Maxwell stress dominates.

Figure 3b shows the evolution of the box-averaged pressure anisotropy. In the initial growth phase of the MRI, $\Delta p$ grows steadily until all cells are pinned at the mirror boundary. In the turbulent phase it then shows small oscillations around a roughly constant value. The background coloring shows the underlying $4\pi\Delta p/B^{2}$ distribution of cells over time. The inset shows this distribution across the simulation domain at a selected time ($t\Omega=150$). Most cells are pinned at the microinstability limits, the majority being on the mirror side, giving an overall positive pressure anisotropy.

Figure 3c shows the evolution of heating due to isotropic and anisotropic diffusion, normalized by the total dissipation,

[TABLE]

where $D_{\eta}$ is the resistive heating,

[TABLE]

$D_{\nu}$ is the isotropic viscous heating,

[TABLE]

and $D_{\mu_{\rm B}}$ is the anisotropic viscous heating,

[TABLE]

The pressure anisotropy $\Delta p$ in equation (19) is computed with mirror and firehose limiters included. In this simulation with $\mathrm{Re}=1500$, $\mathrm{Pm}=2$ and $\mathrm{Re_{B}}=0.75$, resistive heating dominates, followed by anisotropic viscous heating and isotropic viscous heating.

3.2.1 Dependence on $\mathrm{Re_{B}}$

One might expect that anisotropic viscous transport and anisotropic viscous heating will depend primarily on our choice of anisotropic viscosity $\mu_{\rm B}$. We show the dependence on Braginskii Reynolds number for simulations with $\mathrm{Re}=1500$ and $\mathrm{Pm}=2$ in Figure 4. Figure 4a shows that /increases with increasing anisotropic viscosity at large $\mathrm{Re_{B}}$, reaching a plateau at small $\mathrm{Re_{B}}$. Figure 4b shows the dependence of the time-averaged heating fractions on $\mathrm{Re_{B}}$. $D_{\mu_{\rm B}}/D$ has a dependence similar to /, increasing with increasing $\mu_{\rm B}$ at large $\mathrm{Re_{B}}$ and approaching an approximately constant value at small $\mathrm{Re_{B}}$. This is qualitatively similar to the results found by Squire et al. (2018) for strong Alfvénic turbulence. Note, however, that the final values of /and $D_{\mu_{\rm B}}/D$, reached at $\mathrm{Re_{B}}\lesssim 1$, depend on $\mathrm{Re}$ and $\mathrm{Re_{M}}$ (see Figures 10 and 11).

The plateaus in Figure 4 at small $\mathrm{Re_{B}}$ (high $\mu_{\rm B}$) are related to the fact that, due to the presence of limiters, the effective $\mu_{\rm B}$ saturates. Anisotropic transport and heating are most sensitive to the choice of anisotropic viscosity at small $\mu_{\rm B}$, when most fluid cells have $\Delta p$ within the microinstability limits (eq. 5 & 6). But once $\mu_{\rm B}$ is sufficiently large such that most cells lie outside of the limiter region, /and $D_{\mu_{\rm B}}/D$ reach a plateau. We illustrate this in Figure 5a, where we show the distributions of $3\mu_{\rm B}\bm{\hat{b}\hat{b}:\nabla U}4\pi/B^{2}$ for the different choices of $\mathrm{Re_{B}}$. This is the pre-limiter $\Delta p$ distribution divided by (twice) the magnetic energy, before any hard-wall limits are applied to $\Delta p$ (the distribution with limiters is shown in the inset of Figure 3b). The anisotropic stress and anisotropic heating fraction reach an almost constant value once the pressure anisotropy distribution lies mostly outside of the dashed vertical lines denoting the mirror and firehose limits. We find that typically $\mathrm{Re_{B}}\lesssim 3$ is enough to approach the asymptotic value. This can be explained as follows: the presence of limiters causes the effective $\mu_{\rm B}$ to saturate when $3\mu_{\rm B}\bm{\hat{b}\hat{b}:\nabla U_{0}}=-3\mu_{\rm B}b_{x}b_{y}S\sim B^{2}/4\pi$, or equivalently when $\mathrm{Re_{B}}\sim 4\pi S^{2}L_{z}^{2}/B^{2}$. This gives $\mathrm{Re_{B}}\sim$ a few for typical turbulent energy densities, which explains the plateaus in Figure 4.

Figure 5b shows that while the $3\mu_{\rm B}\bm{\hat{b}\hat{b}:\nabla U}4\pi/B^{2}$ distribution becomes broader, the distribution of $\bm{\hat{b}\hat{b}:\nabla U}4\pi/\langle B^{2}\rangle$ becomes narrower with increasing $\mu_{\rm B}$. In addition, large negative $\bm{\hat{b}\hat{b}:\nabla U}$ is more suppressed in simulations without firehose limiter. This can be explained by the work of Squire et al. (2018), who demonstrated that anisotropic viscosity acts to minimize field-line stretching $\bm{\hat{b}\hat{b}:\nabla U}$ to resist changes in magnetic-field strength and make the flow “magneto-immutable”.

3.2.2 Dependence on $\mathrm{Re}$ and $\mathrm{Re_{M}}$: composite simulations

Is isotropic dissipation important for the level of anisotropic transport? It is well known that the Maxwell and Reynolds stresses show a strong dependence on isotropic viscosity and resistivity. Given the interdependence of the pressure anisotropy and velocity-field gradients, it is plausible that $\alpha_{\rm A}$ will also be sensitive to the choice of isotropic Reynolds numbers. Our fully evolved Braginskii simulations (Table 1) tentatively suggest that the relative importance of $\alpha_{\rm A}$ is greater at low $\mathrm{Re}$ and $\mathrm{Re_{M}}$. This is best seen from the two simulations with $\mathrm{Pm}=1$: /decreases from 0.47 for $\mathrm{Re}=750$ to 0.25 for $\mathrm{Re}=4500$, a change significantly larger than the characteristic temporal fluctuations. To explore this dependence in detail and cover a broad range of viscosities and resistivities, we first make use of our composite MHD–Braginskii MHD simulations (see Section 2.3) and then simulations with hyperdiffusion.

Our composite simulations are summarized in the bottom part of Table 1. While we focus on simulations with $\mathrm{Pm}=1$, we also include a simulation with a larger Prandtl number, for which we used $\mathrm{Re}=10500$, $\mathrm{Re_{M}}=42000$. The presented values of $4\pi\langle\Delta p\rangle/\langle B^{2}\rangle$, $\langle 4\pi\Delta p/B^{2}\rangle$ and /are temporal averages, where the averaging is started at time $t\Omega=101$, a time $\Omega^{-1}$ after anisotropic viscosity is added to the system.

We show the evolution of the composite simulations in Figure 6, using four different $\mathrm{Re}$ with fixed Prandtl number $\mathrm{Pm}=1$. The $\mathrm{Re}=750$, 3000, 10500 simulations were performed at resolution $384\times 192\times 96$, while for $\mathrm{Re}=21000$ we went up in resolution to $768\times 384\times 192$. We checked that our $384\times 192\times 96$ simulations are converged by also running $\mathrm{Re}=10500$ at resolution $576\times 288\times 144$. Figure 6a shows the evolution of the box-averaged pressure anisotropy; the anisotropic stress evolution is given in Figure 6b. The evolution of the anisotropic heating fraction $D_{\mu_{\rm B}}/D$ is shown in Figure 6c. The dashed vertical line indicates the time when MHD snapshots were restarted using Braginskii MHD with $\mathrm{Re_{B}}=0.75$. The MHD $\Delta p$, $\alpha_{\rm A}$ and $D_{\mu_{\rm B}}/D$ are computed from the MHD flow fields using the same $\mu_{\rm B}$ that is used for the Braginskii runs, even though $\Delta p$ is not dynamically present in MHD.

The rapid changes in $\Delta p$, $\alpha_{\rm A}$, $D_{\mu_{\rm B}}/D$ and subsequent plateau are a convincing demonstration that we reach the Braginskii state very quickly. What is driving the abrupt transition? To understand this, it is instructive to look at $4\pi\Delta p/B^{2}$ snapshots, before and after anisotropic viscosity is introduced. The upper panels of Figure 7 show this for the $\mathrm{Re}=\mathrm{Re_{M}}=21000$ simulation. In the MHD snapshot in Figure 7a, the vast majority of cells are pinned at the mirror/firehose limit. Moreover, there are many “opposite” cells in close proximity of each other. In these regions anisotropic viscosity operates on a short timescale, trying to eliminate the strong gradients. Doing so produces the smoother distribution of $\Delta p$ depicted in Figure 7b and causes a rapid change in the volume-averaged $\Delta p$. The jump is less pronounced at low $\mathrm{Re}$, because the MHD flow field has already been smoothed by isotropic dissipation, thus diminishing the dynamical importance of anisotropic viscosity. In the bottom panels of Figure 7 we show how in the process of changing the statistics of $\Delta p$, Braginskii viscosity also reduces small-scale gradients in the velocity field.

The damping of high-$k$ modes in the Braginskii velocity field leaves an imprint on the power spectrum. We show the kinetic and magnetic energy spectra of the composite simulation with $\mathrm{Re}=21000$, $\mathrm{Pm}=1$ in Figure 8, both for the MHD part and the Braginskii MHD ($\mathrm{Re_{B}}=0.75$) part of the run. Figure 8a shows that the kinetic-energy spectra in both models have spectral slopes close to $k^{-3/2}$, the MHD case being slightly shallower. There is also extra suppression of high-$k$ velocity fluctuations in the Braginskii case, which is consistent with the snapshot in Figure 7d. Figure 8b demonstrates that the magnetic-energy spectra are hardly modified in Braginskii MHD.

3.2.3 Dependence on $\mathrm{Re}$ and $\mathrm{Re_{M}}$: magneto-immutable turbulence

In Figure 9a we show the pre-limiter anisotropic pressure distribution, before (MHD, dotted line) and after (solid line) the transition. Braginskii viscosity drives more cells into the region inside of the microinstability limits (dashed vertical lines), while damping the tails of the distribution (note that it is the projection of $\bm{\nabla u}$ onto the magnetic-field direction, $\bm{\hat{b}\hat{b}:\nabla u}$, that is strongly modified by anisotropic viscosity, $\bm{\hat{b}\hat{b}}$ and $\bm{\nabla u}$ alone are only mildly affected). In addition, Figure 9a clearly shows that at low $\mathrm{Re}$ we get a distribution that is heavily skewed towards positive $\bm{\hat{b}\hat{b}:\nabla U}$. At higher $\mathrm{Re}$, the distribution becomes more symmetric, leading to a smaller average pressure anisotropy, which explains the results in Figure 6.

The results in Figure 6 and Figure 9a can be explained qualitatively as follows. At low $\mathrm{Re}$, high wavenumber modes of the flow field are efficiently damped by isotropic dissipation and so the fluctuating part of the $\bm{\hat{b}\hat{b}:\nabla U}$ distribution is narrow. Because the shear part of $\bm{\hat{b}\hat{b}:\nabla U}$ is strongly skewed towards positive values (since $\langle\bm{\hat{b}\hat{b}:\nabla U_{0}}\rangle=\langle-\frac{3}{2}b_{x}b_{y}\Omega\rangle>0$), the sum of the two is also going to be biased towards positive values. This leads to a positive and appreciable $\Delta p$ approaching the mirror threshold.

Increasing $\mathrm{Re}$ and $\mathrm{Re_{M}}$ means that there are higher wavenumber modes in the turbulent velocity field. In particular, note that $\bm{\nabla u}\sim k^{3/4}$ for a $k^{-3/2}$ spectrum (see Figure 8a), so it increases with increasing resolution. As a result, the fluctuating part of $\bm{\hat{b}\hat{b}:\nabla U}$ becomes more important and can drive more cells towards the firehose limit. Moreover, in Figure 9b we show that the $\bm{\hat{b}\hat{b}:\nabla u}$ distribution actually has a negative skew for large-$\mathrm{Re}$ turbulence in Braginskii MHD, so that it cancels to a large extent the positive $\bm{\hat{b}\hat{b}:\nabla U_{0}}$ from the mean shear.

Figure 9c offers insight into the physics behind the negative skew of the $\bm{\hat{b}\hat{b}:\nabla u}$ distribution. The 2D histogram shows how in the $\mathrm{Re}=21000$, $\mathrm{Re_{B}}=0.75$ simulation the plasma rearranges itself to produce a $\bm{\hat{b}\hat{b}:\nabla u}$ that locally balances the shear. The fact that the turbulence counters the largely positive $\bm{\hat{b}\hat{b}:\nabla U_{0}}$ can be attributed to anisotropic viscosity causing the rearrangement of the flow field so as to resist changes in $B$ by minimizing $\bm{\hat{b}\hat{b}:\nabla U}$ (Squire et al. 2018). We interpret the results of Figure 9 as a consequence of this magneto-immutability: the turbulence can more effectively cancel the field-line stretching of the mean shear ($\bm{\hat{b}\hat{b}:\nabla U_{0}}$) at large $\mathrm{Re}$, when the plasma is less constrained by isotropic dissipation. Anisotropic viscosity does also rearrange the turbulence at low $\mathrm{Re}$ and $\bm{\hat{b}\hat{b}:\nabla u}$ is able to locally cancel $\bm{\hat{b}\hat{b}:\nabla U_{0}}$ to an appreciable extent (see Figure 9a for $\mathrm{Re}=750$), but the effect is more pronounced at high $\mathrm{Re}$.

As the $\bm{\hat{b}\hat{b}:\nabla U}$ distribution becomes broader and more symmetric with increasing $\mathrm{Re}$, comparable numbers of cells land on the mirror and firehose sides. As a result, $\langle 4\pi\Delta p/B^{2}\rangle$ and /decrease in value, as in Figure 6. /can nevertheless remain more positive in comparison, primarily due to $\langle-B_{x}B_{y}\rangle$ typically being larger when averaged over cells at the mirror limit than cells at the firehose limit. It is not entirely surprising that the Maxwell stress is different at the two microinstability boundaries. For example, where the Maxwell stress is negative, the mean shear drives cells towards the firehose side (as $\bm{\hat{b}\hat{b}:\nabla U_{0}}<0$), whereas $\bm{\hat{b}\hat{b}:\nabla U_{0}}$ is skewed towards the mirror side where $-B_{x}B_{y}>0$. The microinstabilities also significantly affect the dynamical effects of the Maxwell stress: at the firehose limit there is effectively no magnetic tension, and the Maxwell and anisotropic stresses cancel, while at the mirror boundary the effective magnetic tension is enhanced by a factor $(1+4\pi\Delta p/B^{2})$.

The values of $\langle 4\pi\Delta p/B^{2}\rangle$ and /as a function of isotropic Reynolds number are plotted in Figure 10a and Figure 10b respectively. In addition to $\langle 4\pi\Delta p/B^{2}\rangle$, Figure 10a also shows the values of $4\pi\langle\Delta p\rangle/\langle B^{2}\rangle$ as filled semi-transparent diamonds to demonstrate that the exact choice of averaging does not affect our conclusions. We show all our simulations with $\mathrm{Pm}=1$, which includes both full Braginskii (black) and composite MHD–Braginskii MHD simulations (red). The two procedures give consistent results in the shared range of $\mathrm{Re}$, which supports the plausibility of composite simulation predictions at large Reynolds numbers (see also Appendix B).

Figure 10 shows that anisotropic transport and anisotropic pressure decrease significantly as we go to higher $\mathrm{Re}$ and $\mathrm{Re_{M}}$111The trend that $\alpha_{\rm A}/\alpha_{\rm M}$ and $4\pi\langle\Delta p\rangle/\langle B^{2}\rangle$ decrease with increasing Reynolds numbers also seems to be present in our two simulations without firehose limiter (see Table 1; note, however, that the two simulations have different $\mathrm{Pm}$). Quite notably, in spite of a large negative pressure anisotropy in the simulation with $\mathrm{Re}=10500$ and $\mathrm{Pm}=1$, anisotropic viscous transport remains positive.. It remains unclear how they will behave as we let $\mathrm{Re}$, $\mathrm{Re_{M}}$ $\rightarrow\infty$. Unfortunately probing larger isotropic Reynolds numbers is beyond our current computational capabilities.

In spite of $4\pi\langle\Delta p\rangle/\langle B^{2}\rangle$, $\langle 4\pi\Delta p/B^{2}\rangle\rightarrow 0$ at large $\mathrm{Re}$, anisotropic pressure is an important source of dissipation. This is because regions at the firehose and mirror boundaries contribute positively to the anisotropic viscous heating rate (eq. 19), as $\Delta p$ and $\bm{\hat{b}\hat{b}:\nabla U}$ have the same sign (see eq. 4). In Figure 11a we show that anisotropic viscosity becomes the dominant source of dissipation at large $\mathrm{Re}$, even though the volume-averaged pressure anisotropy is a steadily decreasing function of $\mathrm{Re}$ (Figure 10).

3.2.4 Simulations with hyperdiffusion