Local simulations of MRI turbulence with meshless methods

TL;DR

This paper explores the application of meshless MHD methods, specifically SPH and MFM, to local MRI turbulence simulations, revealing their capabilities and limitations in capturing turbulence and dynamo processes.

Contribution

It demonstrates the feasibility of using meshless MHD schemes for MRI turbulence simulations and compares their performance to traditional grid-based methods.

Findings

Both schemes can simulate MRI turbulence in unstratified boxes with net flux.

MFM reproduces MRI dynamo and butterfly diagram in stratified simulations.

SPH MHD develops unphysical strong toroidal fields due to numerical issues.

Abstract

The magneto-rotational instability (MRI) is one of the most important processes in sufficiently ionized astrophysical disks. Grid-based simulations, especially those using the local shearing box approximation, provide a powerful tool to study the ensuing nonlinear turbulence. On the other hand, while meshless methods have been widely used in both cosmology, galactic dynamics, and planet formation they have not been fully deployed on the MRI problem. We present local unstratified and vertically stratified MRI simulations with two meshless MHD schemes: a recent implementation of SPH MHD (Price2012), and a MFM MHD scheme with a constrained gradient divergence cleaning scheme, as implemented in the GIZMO code \citep{Hopkins2017}. Concerning variants of the SPH hydro force formulation we consider both the "vanilla" SPH and the PSPH variant included in GIZMO. We find, as expected, that the numerical noise inherent in these schemes affects turbulence significantly. A high order kernel, free of the pairing instability, is necessary. Both schemes can adequately simulate MRI turbulence in unstratified shearing boxes with net vertical flux. The turbulence, however, dies out in zero-net-flux unstratified boxes, probably due to excessive and numerical dissipation. In zero-net-flux vertically stratified simulations, MFM can reproduce the MRI dynamo and its characteristic butterfly diagram for several tens of orbits before ultimately decaying. In contrast, extremely strong toroidal fields, as opposed to sustained turbulence, develop in equivalent simulations using SPH MHD. This unphysical state in SPH MHD is likely caused by a combination of excessive artificial viscosity, numerical resistivity, and the relatively large residual errors in the divergence of the magnetic field remaining even after cleaning procedures are applied.

Figures

Tables

Table 1. Table 1: Simulations, results and comments

Simulations	Initial fields	Boxsize	Resolution	MHD-methods	EOS	Sections/Figures	Ref
Unstratified	NZ	$H\times 6.28H\times H$	$64\times 360\times 64$	MFM	adiabatic	Sec 3.1/Fig 3	1
		$H\times 4H\times H$	48 elements per H	MFM/PSPH
	ZNZ	$H\times \pi H\times H$	$64\times 200\times 64$	MFM	isothermal	Sec 3.2/Fig 4	2
		$H\times 4H\times 4H$	48/64 elements per H	MFM	adiabatic(c)	Sec 3.2.1/Fig 5	3
Stratified	$B_y$	$\sqrt{2}H\times 4\sqrt{2}\times 24H$	1.5M elements	TSPH/PSPH	adiabatic	Sec 4.2/Fig 7	4
				MFM	ad/iso	Sec 4.3/Fig 9
			3M elements	MFM	adiabatic(c)	Sec 4.4/Fig 11

Equations

$$\frac{\partial \rho }{\partial t}+$$

$$\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}$$

$$+2q{\Omega }^{2}x\hat{\mathbf{x}}-{\Omega }^{2}z\hat{\mathbf{z}}-2\mathbf{\Omega }\times \mathbf{v},$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\frac{\partial \rho u}{\partial t}$$

$$P=\rho u(\gamma -1).$$

$$f(x,y,z)$$

$${\sigma }^2=\frac{1}{3}\int d\mathbf{x}{x}^{2}W(\mathbf{x},H_1).$$

$$Q_z=\frac{{\lambda }_{MRI}}{h}=\frac{2\pi V_{az}}{\Omega h},$$

$$\frac{Q^{\ast }}{Q_{\text{sepr}}}=(\frac{4\pi }{3N_{ngb}}{)}^{\frac{1}{3}}\kappa ,$$

$$N_{ngb}$$

$$n({\mathbf{x}}_i)$$

$$\text{max}\{(s_i-0.75)/0.75\}$$

$${\alpha }_M=\frac{⟨-2B_x{B}_y⟩}{⟨B^2⟩},$$

$$\mathbf{B}=B_0\hat{\mathbf{z}}\sin (2\pi x/H).$$

$$\rho (z)={\rho }_0\text{exp}(-\frac{z^2}{2H^2}),$$

$$divB=\frac{h|\nabla \cdot \mathbf{B}|}{B}.$$

$$\frac{d{u}_{cool}}{dt}=-\frac{u-u_{init}}{{\tau }_{cool}}$$

Full text

accretion, accretion disks — magnetohydrodynamics (MHD) — turbulence — methods: numerical

††journal: ApJS††software: GIZMO code (Hopkins, 2015), VisIt(Childs et al., 2012)

1 Introduction

The turbulence instigated by the magneto-rotational instability (MRI) can transport angular momentum outwards thus enabling accretion in several sources, such as dwarf novae, low mass X-ray binaries, and Active Galactic Nuclei (AGNs). Numerical MHD simulations are necessary to study this highly non-linear problem. Simulations of MRI range from local shearing box simulations, unstratified (e.g. Hawley et al., 1995, 1996; Sano et al., 2004; Simon & Hawley, 2009) and stratified (e.g. Brandenburg et al., 1995; Stone et al., 1996; Miller & Stone, 2000; Davis et al., 2010; Simon et al., 2011), to global simulations (e.g. Armitage, 1998; Hawley, 2000; Steinacker & Papaloizou, 2002; Fromang & Nelson, 2006; Parkin & Bicknell, 2013; Zhu & Stone, 2018). Three-dimensional simulations carried out with different grid-based codes, such as ZEUS (Hawley et al., 1995), Pencil (Brandenburg & Dobler, 2002), RAMSES (Teyssier, 2002; Fromang & Nelson, 2006), ATHENA (Stone et al., 2008), and the spectral code Snoopy (Lesur & Longaretti, 2007) report similar statistics for the turbulence.

Local MRI simulations are especially challenging because the saturated state appears to depend on the small-scale diffusion, be it physical or numerical. For instance, zero-net-flux simulations in unstratified boxes do not converge with increasing resolution, as the turbulent motions reach their peak amplitude near the smallest resolvable scales. But if physical sources of diffusivity are incorporated and resolved, turbulence can die out when the magnetic Prandtl number is too small (Fromang & Papaloizou, 2007; Fromang et al., 2007). The latter dissipation is also sensitive to the vertical aspect ratio of the computational domain ($L_{z}/L_{x}$) (Shi et al., 2015). On the other hand, in net vertical flux simulations angular momentum transport depends on the magnetic Prandtl number yet again, at least when the latter takes values of order unity (Meheut et al., 2015). Vertically stratified shearing box simulations without a net flux also suffer convergence problems (Bodo et al., 2014; Ryan et al., 2017). Adding a net vertical flux, however, can radically change the character of MRI turbulence as, for example, magnetic winds, may be launched. Some of the properties of these winds also suffer from non-convergence (Fromang et al., 2013; Bai & Stone, 2013; Lesur et al., 2013). Finally, non-ideal MHD effects can suppress or radically alter the nature and strength of turbulence (see, e.g. Fleming et al., 2000; Turner et al., 2007; Bai & Stone, 2011; Lesur et al., 2014; Bai, 2014; Simon et al., 2015). Such effects are still under investigation.

Currently there are very few published studies of the MRI undertaken with mesh-free methods (see,e.g., Gaburov & Nitadori, 2011; Pakmor & Springel, 2013; Hopkins & Raives, 2015), and no systematic exploration of MRI properties over the various standard flow and magnetic field configurations routinely examined with grid-based codes. This is despite the frequent use of smoothed-particle magnetohydrodynamics (SPH MHD) (Springel, 2010a; Price, 2012), in particular, to probe galaxy, star, and planet formation (Dobbs et al., 2016; Dolag & Stasyszyn, 2009; Steinwandel et al., 2018; Price & Bate, 2007, 2008, 2009; Forgan et al., 2016). The neglect is perhaps connected to several numerical problems specific to meshless methods, which we now briefly discuss.

In codes without meshes or with arbitrary mesh geometries, minimizing the divergence of magnetic fields (Tricco & Price, 2012; Hopkins, 2016a) is a major challenge. Because proper minimisation of the divergence is hard to achieve, small “magnetic monopoles” can arise, leading to spurious magnetic field reconfiguration, re-connection, and artificial dissipation in neighboring domains. On fixed, rectilinear, regular, non-moving grids, the Constrained Transport (CT) scheme (Evans & Hawley, 1988) can maintain zero divergence to machine precision. Until recently, CT schemes had only been implemented for regular, non-moving meshes, but recently Mocz et al. (2014, 2016) successfully generalized the CT method to moving meshes that adopt a Voronoi tesselation as their volume partition (e.g. those in AREPO, Springel 2010b). However, most Lagrangian or quasi-Lagrangian methods, including moving-mesh as well as particle-based methods or mesh-free finite-volume methods, use “divergence cleaning” schemes to keep $\bm{\nabla}\cdot\bm{B}$ minimal (Powell et al., 1999; Dedner et al., 2002). Tricco & Price (2012) developed improved divergence-cleaning implementations in SPH (adapting the hyperbolic cleaning scheme from Dedner et al. 2002), and showed this could successfully reproduce some standard MHD tests, for example the Orszag-Tang vortex. But in non-linear MRI simulations, and in fact in any regime of MHD turbulence, effective divergence cleaning is especially difficult due to the complex, multi-scale field geometry, hence the latter methods are not guaranteed to work satisfactorily.

SPH also suffers from known numerical dissipation sourced by various terms, including the E0 error (Read et al., 2010), pairing instability (Rosswog, 2015; Dehnen & Aly, 2012) and incorrectly-triggered artificial viscosity (Deng et al., 2017a). It is well-known that this additional numerical dissipation impedes SPH’s capability to model subsonic turbulence, even without magnetic fields (Bauer & Springel, 2012; Hopkins, 2015; Deng et al., 2017b).

The Lagrangian mesh-less finite-volume (MFV) method was developed two decades ago (see, e.g., Vila, 1999; Hietel et al., 2000) and has been significantly improved since (see, e.g. Lanson & Vila, 2008a, b). Recently it has stimulated a growing interest in the astrophysical community (see, e.g., Gaburov & Nitadori, 2011; Hubber et al., 2017). Hopkins (2015) generalized the method in Gaburov & Nitadori (2011) to other mesh-free finite-volume Godunov schemes, including the closely-related “meshless finite-mass” (MFM) method. These methods, similarly to moving mesh methods, attempt to combine advantages of grid-based and particle-based codes. In particular, they can describe subsonic hydrodynamical turbulence reasonably well (with a quality comparable to regular mesh grid-codes; Hopkins 2015), though at greater computational cost; they avoid advection problems in complex flow geometries that are better modeled in the Lagrangian frame, respecting, for example, Galilean invariance, and finally they can naturally extend to self-gravitating flows by exploiting accurate state-of-the-art gravity solvers, such as treecodes, which have native implementations in particle codes. Mesh-less methods have been generalized to MHD (Hopkins, 2016a), and in subsequent work by Hopkins (2016b) a constrained-gradient (CG) divergence cleaning scheme has been developed that can maintain a much smaller $\bm{\nabla}\cdot\bm{B}$ (by $\sim$2 orders of magnitude) compared to hyperbolic divergence cleaning. These methods, as implemented in the public code GIZMO,111The public version of the code, containing all the algorithms used here, is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html have in fact already been used to simulate the MRI in two-dimensional unstratified shearing sheets (Hopkins & Raives, 2015), and these tests have demonstrated that it recovers the correct linear growth rates and behaves similarly to well-tested grid codes (e.g. ATHENA). However, how these methods perform in three dimensions, in stratified configurations, and/or during non-linear saturation, remain untested.

In this paper, we carried out MRI simulations in both unstratified and vertically stratified shearing boxes, with both SPH and MFM MHD implementations as they are implemented in the multi-method GIZMO code, in order to explore the numerical requirements for these methods to treat the MRI in the non-linear regime. We focus on MFM as opposed to MFV or more general moving-mesh schemes (several of which are also implemented in GIZMO and can, in principle, use the same constrained-gradient divergence “cleaning” method) because MFM is designed, such as SPH, to conserve exactly the mass of fluid elements (i.e. there is identically zero advection), so the method is “purely” Lagrangian. This is perhaps the most challenging case for our purpose, since hybrid moving-mesh or MFV-type methods, in which the grid moves but mass fluxes are also allowed, effectively act as a smoothing of grid motion, thus interpolating between the “pure Lagrangian” (constant mesh-motion) and the “pure Eulerian” (fixed-grid) representation of a fluid.

We will explore both the traditional “density-energy” formulation of SPH (named hereafter ‘TSPH’) (Springel, 2005) and the more recently-developed “pressure-energy” formulation (‘PSPH’) (Saitoh & Makino, 2013; Hopkins, 2012). The SPH MHD used here represents the state-of-the-art implementation described by Price (2012) with the advanced artificial viscosity/resistivity switches developed in (Cullen & Dehnen, 2010; Tricco & Price, 2013) and divergence cleaning following (Tricco & Price, 2012). In unstratified shearing box simulations, no significant density contrast is present and we expect TSPH and PSPH to perform similarly. Therefore, we did not run TSPH and PSPH comparisons for this particular setup. In the GIZMO MFM runs we adopt the constrained-gradient divergence cleaning of Hopkins (2016b).

We note that the effect of the initial noise in MFV, which appears to depend on the regularity of the initial particle distribution (Gaburov & Nitadori, 2011), is poorly understood. In general, the dependence of MRI properties on the numerical setup of the initial condition should be expected since MRI is extremely sensitive to numerical dissipation. Due to this added complexity of the initial condition design, although, in principle, MFV-type methods should be less prone to the effect of particle discretization noise during slope limiting as well as in the divergence cleaning step, we defer their scrutiny in the context of MRI to future work.

We start, in Section 2, with a discussion of our shearing box implementation and the role of the smoothing kernel function. We also tested the resolution needed for accurate MRI eigenmode growth. In Section 3, we present unstratified shearing box simulations with and without net vertical flux and in short and tall boxes. Stratified shearing box simulations are described in Section 4 where we compare different simulation setups and the two methods, MFM and SPH. A discussion and conclusion follow in Section 5 and 6.

2 The Shearing Box Approximation

The shearing box is a local model of the equations of motion widely used in MRI simulations to achieve high resolution (Goldreich & Lynden-Bell, 1965; Hawley et al., 1995; Latter & Papaloizou, 2017). One considers a small patch of a disk centered at a radius $R$ and rotating at the angular velocity $\Omega(R)$. In the corotating frame one installs a Cartesian geometry at the box centre, using $x$ and $y$ to represent the radial and the azimuthal directions respectively. In compressible ideal MHD, the governing equations are

[TABLE]

where $\hat{\bm{x}}$ and $\hat{\bm{z}}$ are the unit vectors in the $x$ and $z$ directions, and $\rho$, $u$, $P$, $c_{s}$, $\bm{v}$ represent the density, specific internal energy, gas pressure, sound speed and velocity respectively. The tidal acceleration $2q\Omega^{2}x\hat{\bm{x}}$ in equation 2 comes from the expansion of the effective potential (gravitational plus centrifugal). The constant $q\equiv-d\ln\Omega/d\ln R$, and for a Keplerian disk $q=1.5$. The vertical component of the star’s gravity is represented by $-\Omega^{2}z\hat{\bm{z}}$ which, if included, results in a vertical density stratification with a scale height of $H=c_{s}/\Omega$, where $c_{s}$ is the initial sound speed. In simplified models examining motions confined near the disk midplane this term can be dropped. The ratio between the gas pressure and magnetic energy, $\beta\equiv P/(B^{2}/8\pi)$, is a dimensionless measure of the magnetic field strength.

We assume an ideal gas equation of state (EOS),

[TABLE]

We choose $\gamma=5/3$ except when we set $\gamma=1.001$ to mimic an isothermal EOS. In particular, we have one stratified simulation with $\gamma=1.001$ to show how such a soft EOS exacerbates long-term numerical dissipation.

2.1 Shearing Box Boundary Conditions

The computation domain is a rectangular prism with sides of length $L_{x}$, $L_{y}$ and $L_{z}$. In unstratified boxes, the domain is periodic in $y$ and $z$, and shear periodic in $x$. These boundary conditions can be expressed mathematically for a fluid variable $f$ as

[TABLE]

They apply to the thermodynamic variables and the perturbed fluid velocity, where the background is $\bm{v}_{0}=-q\Omega x\hat{\bm{y}}$. The implementation of the shearing periodic boundary conditions in Lagrangian codes is relatively easy since we do not need to extrapolate fluid quantities to ghost zones, as in grid-codes. When a fluid element (“particle”) moves across the radial boundary it reappears at the other radial boundary with a velocity offset added to its azimuthal velocity.

In vertically stratified simulations we apply outflow boundary condition in the $z$ direction by removing any element whose smoothing length is larger than 1.2H. This yields a density floor of about 0.0002 in code units (see below).

2.2 Equilibrium Tests and the Kernel Function

The shearing box admits the following simple equilibrium: $\bm{v}=-\bm{v}_{0}$, $\rho=$ constant. To test whether the code properly describes this state, in addition to the shearing periodic boundary conditions, we conduct an MFM simulation using this as an initial condition. We employ a cubic box of one disk scale height per side. In the calculation we choose units so that $\Omega=1$, $c_{s}=1$, and $\rho=1.$ At a resolution of $48\times 48\times 48$ elements with the Wendland C4 kernel (200 neighbours, $N_{ngb}=200$), the equilibrium can be maintained to machine precision for the duration of the simulation ($\sim 200\Omega^{-1}$) with no sign of breakdown.

We next reran the simulation using the cubic spline kernel (55 neighbours) and found that the radial velocity becomes non-zero and the perfect lattice breaks into a glass configuration. The velocity errors are a few percent of the sound speed. In this case, elements form pairs as shown in figure 1. In SPH this pairing (or clumping) instability (Springel, 2010a; Price, 2012) happens with any kernel whose Fourier transform is negative for some wave vectors, at sufficiently large neighbour number (Dehnen & Aly, 2012). It would appear then that MFM also suffers a similar instability if “too many elements” are included in the kernel radius of compact support (i.e. if one does not, as one should, use higher-order kernels with a larger number of neighbors).222The interpretation of the pairing instability in MFM is slightly different to that in SPH – the kernel function is used in MFM to define the volume partition between neighboring resolution elements. If one takes a low-order kernel, say, the cubic spline, and forces its radius of compact support to enclose too many elements (such that the mean inter-element separation is much smaller than the kernel function width), the effective faces between elements essentially “overlap” into a single face (which becomes ill-defined). The Wendland C4 kernel does not suffer from these issues at this “enclosed neighbor number” (see e.g. Dehnen & Aly, 2012). It also helps to keep elements well ordered, which is crucial for accurate gradient estimation in any unstructured method (Rosswog, 2015). MRI turbulence is generally subsonic (except in the disk corona of stratified box), so we always use the Wendland C4 kernel to minimize numerical noise/dissipation (except when noted).

2.3 Resolution

Ideally we would like to estimate the shortest lengthscale adequately resolved in a simulations $h$, the number of resolution elements required to furnish that level of accuracy, and ultimately the numerical resources required by a given code to achieve this.

Unfortunately, it is not straightforward to determine the shortest resolvable lengthscale, even for grid codes (which possess an unambiguous grid spacing $\delta z$). Particle codes exhibit several lengthscales that might be used as references in this context. Perhaps the most obvious choice is the mean particle separation $\overline{\delta}$. Of course, this is rather crude, because particle codes are intrinsically adaptive, with denser regions better resolved than less dense regions. But even if we put this concern to one side, we must emphasise that resolution is not so much determined by the size of $\overline{\delta}$ as it is by the kernel volume, since this is where the interpolation occurs (which is instrumental in defining a volume element in MFM). We may then be tempted to use, instead of $\overline{\delta}$, the kernel support radius $H_{1}$, because the particles within a kernel form a ‘resolution unit’ together. But more accurate estimates exist, as we now discuss.

In SPH, an alternative definition for the resolution scale is $h=2\sigma$ where $\sigma$ is the standard deviation of the weighting kernel function $W(\bm{x},H_{1})$, as defined in Dehnen & Aly (2012),

[TABLE]

We hence will denote this lengthscale by $h_{\text{SPH}}$.

In MFM, for a well-chosen kernel (i.e. $H_{1}$ is within a factor of a few rms inter-particle separation, so that faces are well-defined), an alternative estimate of $h$ is the face-area weighted inter-neighbor separation (where the face areas are themselves determined by the volume partition from the kernel function; Hopkins 2015). Note that, consequently, h so defined is the scale upon which the Riemann problem is solved, diferent to both $H_{1}$ and $\overline{\delta}$. For the Wendland C4 kernel parameters adopted here, this yields $h\approx 0.45H_{1}$, so we will use this value throughout. We denote this by $h_{\text{MFM}}$.

2.4 MRI Quality Factor

In grid-code MRI simulations, the number of cells per fastest growing mode’s wavelength is an often used resolution metric (Hawley et al., 2011; Parkin & Bicknell, 2013). It can potentially tell us where and when in a simulation the MRI is adequately or inadeqately resolved.

We may define a quality parameter (Noble et al., 2010; Hawley et al., 2011) via

[TABLE]

where $V_{az}$ is the $z$ component of the Alfven velocity and $h$ is the (vertical) resolution length. $\lambda_{MRI}$ is close to but not exactly the fastest growing linear mode’s wavelength, $\lambda_{fastest}=\sqrt{16/15}\lambda_{MRI}$, in the presence of a net vertical flux. Perturbations with wavelengths smaller than $\lambda_{MRI}/\sqrt{3}$ are stable for the same configuration. $Q_{z}$ can be measured and averaged over the disk body during the saturated state. It can also be defined with respect to other coordinate directions, most often in the $y$ direction. Regions where the plasma beta is high are of lower quality factor, and in these regions there may be insufficient resolution.

We must state from the outset that a quality factor, so defined, is a rather crude measure of how well the turbulence is resolved. First, it is based on the linear theory of the net-flux MRI set-up and hence may not be generally applicable during the ensuing nonlinear saturation; certainly its relevance for zero-net flux simulations is debatable. Second, it only describes whether the input scale of the turbulence is resolved and has nothing to say about the ensuing turbulent cascade on smaller scales. If $Q_{z}\gtrsim 1$ then there is no inertial range to consider, and, in addition, the input and dissipative scales are directly adjacent: strictly, there is no real turbulence, but rather a monoscale chaotic flow. Nonetheless, vertically stratified shearing box simulations indicate that $Q_{z}>10,Q_{y}>20$ ensures the convergence with resolution of some large-scale average flow quantities (Hawley et al., 2011).

In our small number of grid-based simulations we substitute $h=\delta z$, the grid spacing. In our particle simulations we use both $\overline{\delta}$ and $h_{\text{SPH}}$ or $h_{\text{MFM}}$ for comparison. We denote by $Q_{\text{sepr}}$ the quality factor when $\overline{\delta}$ is used for h, and denote by $Q^{*}$ the quality factor when either $h_{\text{SPH}}$ or $h_{\text{MFM}}$ is used.

For MFM the two quality factors are related by

[TABLE]

where $\kappa\equiv H_{1}/h_{\text{MFM}}$ and is related to the type of kernel used in the simulation. The cubic spline kernel can achieve good density estimation with a small number of neighbours $N_{ngb}$ (around 42), where we define

[TABLE]

where $m_{i}$ is the mass of the $i$’th particle. The cubic spline yields $Q^{*}$ and $Q_{\text{sepr}}$ that are relatively close (the latter only a factor 0.9 of the other). But the Wendland C4 kernel ($N_{ngb}=200$) gives $Q^{*}=0.6Q_{\text{sepr}}$.

2.5 Channel flow growth

The linear MRI eigenmodes in a net vertical flux are called channel flows. Being nonlinear solutions, these eigenmodes will grow to large and nonlinear amplitudes before being destroyed, in the first instance, by parasitic instabilities (Goodman & Xu, 1994). A robust turbulence then ensues, which sometimes exhibits the recurrent generation and destruction of the channels. In this section we measure the growth rates of the simulated channel flow in SPH and MFMand compare it with the theoretical value. In particular, we calculate the error as a function of resolution (in fact by, $Q$), to see if certain scalings can be discerned. For comparison we also obtain error estimates for the finite volume Godunov code, ATHENA (Stone & Gardiner, 2010), with second order reconstruction and either the Roe and HLLD solvers. Our primary aim here is to compare the performance of the codes on this linear problem not with respect to computational effort directly, which we address in detail in Section 5.1, but with respect to $h$. Is a plausible estimate of the minimum resolvable lengthscale a useful diagnostic for measuring the accuracy of MRI growth, or are other things going on? If the importance of $h$ is overwhelming we may in fact expect all three codes to perform similarly on this test.

We initialize a box of size $H\times H\times H$, threaded by uniform vertical background fields of magnitude $B_{0}$. We set $\beta=84$, so that the fastest growing channel mode just fits into the box, and $\gamma=1.001$ so the gas is effectively isothermal. The initial amplitude of the channel mode is $0.001B_{0}$. The theoretical growth rate of the fastest channel mode is $0.75\Omega^{-1}$. The simulations are run for $8\Omega^{-1}$ so that the channel mode is smaller than $0.4B_{0}$ at the end of the simulation. We calculate the growth rate using the magnitude of the magnetic field at two consecutive snapshots taken every $0.5\Omega^{-1}$. The growth rate relative error is defined as

[TABLE]

where $s_{i}$ is the $i$th measure of the growth rate.

In figure 2, we plot the growth rate error as a function of $Q_{\text{sepr}}$, i.e. the (mean) number of resolution elements per $\lambda_{MRI}$, be they grid cells (in ATHENA) or the average particle spacing (in SPH and MFM). In MFM the Wendland C4 kernel is employed in order to avoid pair instability. As is clear, and to be expected, the errors decline with increasing resolution. They are less than 1% when $Q_{\text{sepr}}>32$ in MFM. TSPH captures the MRI better than MFM in the low resolution simulations but it converges slower. It is known that SPH has zeroth order errors that only vanish when both $Q_{sepr}$ and $N_{ngb}$ approach infinity ($N_{ngb}$ is fixed here) (Read et al., 2010; Zhu et al., 2015).

ATHENA performs better than both SPH and MFM, if we are to directly compare number of grid points against number of particles. But, as explained in the previous two subsections, even if such a comparison may reflect relative numerical effort it doesn’t do justice to the particle codes, for whom the mean particle spacing is not equal to the shortest resolvable lengthscales. In figure 2 we also plot the relative growth rate error for both SPH and MFM versus $Q^{*}$, the number of resolution lengths $h$ per MRI wavelength (less than the number of particles). These dashed curves now sit on top or adjacent to each other. If we are to consider resolution this way, actually the codes perform comparably; ATHENA and MFM (which share the same Riemann solver) agree exceptionally well. In contrast, SPH possesses a different and adverse scaling in comparison, indicative of its sensitivity not only to the number of particles but to the number of near neighbours. This poorer performance we will see again later in nonlinear tests of MRI saturation.

Finally, we say a few words dealing with the cubic spline kernel in MFM, which for smaller $N_{ngb}$ yields a larger effective quality parameter $Q^{*}$ (see section 2.3). Indeed, it outperforms the Wendland C4 kernel when the resolution is low. However, when $Q_{sepr}>20$ the growth rate errors increase due to the pairing instability (which itself is resolution dependent, see Dehnen & Aly (2012)). This numerical noise (see section 2.2) can dominate over the weak channel mode in the early stage. The channel modes does eventually outcompete this noise but the errors in gradient estimation lead to extra dissipation (see section 4.2) which is hard to quantify.

3 Unstratified Shearing Box Simulations

We summarize all the simulations we undertook in both unstratified and stratified boxes in table 1, which includes key parameters, physical and numerical configurations, comments, and references. Further details can be found in the referenced subsections.

3.1 Net-Vertical-Flux Simulations

We first ran an unstratified shearing box simulation with net vertical flux similar to the fiducial model of Hawley et al. (1995). This is the simplest 3D MRI setup. For such configuration we are able to reproduce the main features of previous grid-based simulations using high resolution MFM simulations. The box is of size $H\times 6.28H\times H$ and threaded by vertical fields with $\beta=400$. We used a resolution of $64\times 360\times 64$ elements which corresponds to 28 elements per $\lambda_{MRI}$, and hence roughly 17 $h_{\text{MFM}}$ per $\lambda_{MRI}$ (or roughly 9 $H_{1}$ per $\lambda_{MRI}$). In comparison, Hawley et al. (1995) employed about 16 grid points per MRI wavelength. The box size also affects the simulated turbulence: smaller box tends to have stronger outbursts in the turbulent state (Bodo et al., 2008; Lesaffre et al., 2009). We set two other runs in a box of $H\times 4H\times H$ with $\beta=330$ ($H=2\lambda_{MRI}$), and with either PSPH and MFM using 48 elements per H. We added random velocity perturbations (5% of the sound speed) to the shear flow at initialisation. The simulations are run for 11 orbits with $\gamma=5/3$.

To characterize the saturated turbulence we plot several density weighted-averaged quantities in figure 3. We take the arithmetic average of fluid variables at all the MFM fluid elements so the average is naturally density-weighted because of the adaptive nature of GIZMO. In unstratified turbulence, the density fluctuations are small so the density-weighted average should be close to the volume average in previous studies. This should be the case also for stratified turbulence, another situation in which we will apply this averaging method (see next section 4), because the stress is almost independent of the density when $|z|<2\sqrt{2}H$ (Simon et al., 2011). In the $L_{y}=6.28H$ simulation (red curves), both the magnetic energy and stresses are in agreement with the results of Hawley et al. (1995), which were obtained with an Eulerian code. The ratio of the Maxwell stress to the magnetic pressure,

[TABLE]

is about $0.5$, similar to the aforementioned previous adiabatic calculations. The agreement between the older ZEUS runs and MFM may have been expected because the resolution is similar in the two (about 16 resolution lengths, either $\delta z$ or $h_{\text{MFM}}$, per MRI wavelength). It should be noted, however, that the ZEUS code is more diffusive than ATHENA - our point of comparison in the linear growth tests.

The simulations in smaller boxes, $L_{y}=4H$, show stronger bursts in stresses and magnetic energy because fewer active (non-axisymmetric) modes can fit in the box leading to an artificial truncation of the participating modes in the nonlinear dynamics. As a consequence, the system is nearer criticality and single channel modes can intermittently dominate.

The internal energy increases due to the turbulent dissipation. This is most significant in the two $L_{y}=4H$ simulations as they exhibit the strongest bursts from the channel flows: since these flows achieve large amplitudes, when they break down a great deal of energy is dissipated into heat. The PSPH simulation, in particular, is some four times ‘hotter’ than the large box MFM simulation. The PSPH run undergoes also a much greater increase of internal energy compared to the equivalent MFM run (see Figure 3, upper right panel) due to stronger channel activity near the beginning of the run (signaled by the very large initial spike in the various diagnostics shown in Figure 3). The dominance of channels in the PSPH run early on suggests that the system is closer to marginal stability than in MFM; this is probably due to the additional numerical diffusivity in PSPH. Note that the plasma $\beta$ increases substantially as the gas is heated up, and the boxes will ultimately approach the incompressible zero-net-flux regime. This explains why, in general, the bursts become less powerful as the simulations continue. Finally, for the net-flux unstratified setup we only ran PSPH, and not TSPH; since there are no steep density gradients, we expect no significant difference due to the actual formulation of the SPH hydro force.

3.2 Zero-Net-Vertical-Flux Simulations in a ‘Standard’ Box

We run a standard zero-net-flux unstratified box simulation with MFM with no explicit physical dissipation terms (Stone et al., 1996; Fromang & Papaloizou, 2007). We initialize a box of size $H\times\pi H\times H$ with $64\times 200\times 64$ elements, threaded by magnetic fields,

[TABLE]

The field strength $B_{0}$ is chosen so that the volume averaged $\beta$ equals 400. We use an isothermal EOS ($\gamma=1.001$) to align with Fromang & Papaloizou (2007). MRI turbulence is sensitive to the nature of both physical and numerical dissipation: without physical viscosity and resistivity, Fromang & Papaloizou (2007) found that zero-net-flux MRI turbulence was driven to smaller scales as resolution increased with no sign of convergence. Fromang et al. (2007) showed that the saturated state depended on the magnetic Prandtl number when a source of diffusivity, physical or numerical, is present; if this was too low, turbulence would die after some period of time. With these problematic results in mind, we will now assess how well a zero-net-flux MRI setup can be modelled by a meshless code.

We plot the averaged magnetic energy, Maxwell stress, and $\alpha_{M}$ in figure 4. In contrast to MRI runs with grid codes, after an initial burst the magnetic fields and magnetic stress rapidly decay. There is no sustained turbulence, as in (Fromang & Papaloizou, 2007) nor is there some period of MRI turbulence before decay, as in Fromang et al. (2007). It is true that MFM smooths fluid variables within a kernel, so the resolution (with respect to $h_{\text{MFM}}$) is roughly half that of the standard simulation with $64\times 200\times 64$ grid cells in Fromang et al. (2007) (see section 2.3). A simulation with $>$128 particles per scale height is prohibitively expensive with MFM (see section 5.1). But even an ‘ideal MHD’ run undertaken with low resolution in a grid code can sustain MRI turbulence (Stone et al., 1996). Thus our result is disappointing.

One way to interpret it is to consider the relative sizes of the numerical resistivity and viscosity. At low resolutions, MFM should have a moderate numerical viscosity, $\nu$ (see appendix A) and relatively large numerical resistivity, $\eta$ (see appendix B), as a consequence the effective numerical Prandtl number $P_{m}=\nu/\eta$ must be small (smaller than 1). It is then perhaps not surprising that the turbulence decays (Fromang et al., 2007). However, the fact that it decays so abruptly might point to a simpler reason: the numerical resistivity is just very large and prohibits turbulence of any kind past the initial spike. Indeed, our MFM run resembles in some respects the Fleming et al. (2000) run with a magnetic Reynolds number of 13000, which after an initial burst abruptly dies off.

3.2.1 Zero-Net-Vertical-Flux Simulations in a Tall Box

It has been shown that when the numerical domain is reshaped, so it exhibits a large vertical aspect ratio ($L_{z}/L_{x}\geq 2.5$), a new more vigorous and cyclical MRI dynamo emerges (Shi et al., 2015) (see also Lesur & Ogilvie (2008)). Importantly, its saturated stress is independent of resolution. To test the effect of a taller box, we redo the simulations in Section 3.2 in a box of size $H\times 4H\times 4H$. We present two simulations with 48 and 64 elements per scale height; their details are shown in figure 5. We set $\gamma=5/3$ and add a cooling term in order to keep the internal energy roughly constant (see equation 19). The numerical dissipation is EOS-related and we expect very fast decay with the isothermal EOS (see section 4.3). This does mean, however, that we can’t make direct quantitative comparison with Shi et al. (2015).

Initially our simulations exhibit turbulence as shown in Figure 5. Moreover they reproduce the averaged toroidal field patterns produced by Shi et al. (2015). However, while the turbulence is sustained for much longer than in the standard box, after some 30-40 orbits the activity ultimately dies out. During the turbulent phase, $\alpha_{M}\sim 0.44$ in the higher resolution tall box simulation, but the averaged magnetic energy $\langle B^{2}/8\pi\rangle/P_{0}\sim 0.01$ is much lower than the values ($>0.1$) obtained with the ATHENA code (Shi et al., 2015). As a result, the stress ($\sim 0.006P_{0}$) is also much smaller.

In tall box simulations, when the magnetic Prandtl number, $P_{m}\geq 4$ the saturated stress is independent of $P_{m}$ while the turbulence vanishes for $P_{m}=1$ with even 128 cells per scale height (Shi et al., 2015). While $P_{m}$ effects might be at play in our simulations it should be noted that our 64 particle simulation is quite low resolution (elements are further smoothed within a kernel). Worse resolution leads to more rapid decay. It is likely that both meshless simulations possess too great a numerical resistivity to support a sustained MRI dynamo.

To summarize, the MFM simulations – at least with the current implementation of the method – do not allow us to sustain the MRI in unstratified zero-net flux simulations, either in short or tall boxes. This is probably the result of either too large a numerical resistivity, or, more generally, too low an effective $P_{m}$, at least at the resolutions we were able to access (see Section 5.1). It will be particularly interesting to explore both (a) higher-resolution MFM simulations (where the resistivity should be lower and $P_{m}$ larger), and (b) simulations using other, closely-related schemes which are not completely “fixed mass” schemes but closer to moving mesh schemes (e.g. the MFV scheme or arbitrarily shearing-mesh schemes with the constrained-gradient divergence cleaning). We did not run zero-net-flux simulations with SPH, because of evidence suggesting that it has an intrinsically larger numerical viscosity, and thus would yield higher Prandtl numbers, quite irrespective of the specific implementation of artificial viscosity (Bauer & Springel, 2012; Deng et al., 2017a, b), and produces substantially larger magnetic field divergence with available divergence cleaning methods (see Section 4.1).

4 Stratified Shearing Box Simulations

In this section we undertake simulations in the vertically stratified shearing box, in which the (leading order) vertical acceleration from the central star’s gravity is incorporated. We initialise the simulations with a Gaussian density profile with uniform temperature

[TABLE]

where $H=c_{s}/\Omega$ is a factor different from $\sqrt{2}c_{s}/\Omega$ in some previous work (Davis et al., 2010; Simon et al., 2011). Note that $c_{s}$ is the initial sound speed; in adiabatic runs the sound speed (and hence scale height) will change. We adopt units so that $H=1,\Omega=1,c_{s}^{2}/\gamma=1$ and use the adiabatic EOS with $\gamma=5/3$. The density profile is sampled using the Monte Carlo method and then relaxed to a glassy configuration. The density errors in the disk body ($-3H<z<3H$) are below the 1% level. We initialize azimuthal magnetic fields with $\beta=25$ in a box of size $\sqrt{2}H\times 4\sqrt{2}H\times 24H$ (the box is extremely tall but no element has $\vline z\vline>6$ in our simulations, see figure 12). Outflow boundary condition are applied but in fact there is no significant outflow and few elements are clipped. Random velocity perturbations $\sim 0.01c_{s}$ are added to seed the instability.

In our fiducial model we use 1.5M elements in total, leading to $\overline{\delta}\sim 0.03$ at the disk midplane. However, due to the adaptive feature of our method, the resolution is lower the further away from the midplane. This helps in order to save some computational resource because high resolution is not needed in the MRI-stable disk corona with strong fields (see figure 12 and also Miller & Stone, 2000). Yet, the nearly zero-flux MRI turbulence in the disk body still requires high resolution and is computationally demanding.

4.1 Divergence Cleaning of Magnetic Fields

Both SPH and MFM are not able to strictly maintain exactly solenoidal magnetic fields naturally, and thus must employ cleaning schemes to keep their divergences minimal. We try to quantify the efficacy of this procedure in this section before showing our main results.

We define the dimensionless divergence of magnetic fields as

[TABLE]

In our unstratified box simulations of section 3, the $divB$ diagnostic is smaller than $10^{-3}$ at the location of most fluid elements in MFM. Divergence control in the stratified shearing box MRI is more challenging, however. We run the fiducial model to compare the level of non-zero divergence in MFM with the CG cleaning (Hopkins & Raives, 2015; Hopkins, 2016a) and TSPH with the hyperbolic divergence cleaning (Tricco & Price, 2012). In figure 6 the hyperbolic cleaning keeps $divB\sim 0.1$ in TSPH while the CG cleaning keeps $divB$ two orders of magnitude lower in MFM. Large $divB$ only occurs at the vertical boundaries and in the weak field regions in MFM; the vertical boundaries are poorly resolved because MFM fluid elements, which are built from particles, are fewer, but there the divergence should have negligible influence on the turbulent disk body since the correlation length of magnetic fields is smaller than H (Davis et al., 2010; Bai & Stone, 2013). We also stress that the maps shown in Figure 6 are quite representative of the differences between SPH and MFM in our tests.

Summarizing, the CG cleaning method in MFM significantly outperforms its competitors here, and we shall see how this is important in the following subsection. As a word of caution, we note that, SPH methods that are recast at least partially in a finite-volume formulation, such as Godunov-SPH (Inutsuka, 2002) , might be amenable to implementations of the GC cleaning method. It would be interesting to explore the latter avenue in order to find out how much better an SPH method can perform once it is equipped with a superior divergence cleaning scheme.

4.2 Unphysical Behaviour in SPH Simulations

We present three SPH simulations, two of which are run with the Wendland C4 kernel and adiabatic EOS but using different SPH formulations. In figure 7 we plot $\alpha_{M}$ and the scaled magnetic energy versus time, in addition to space time diagrams of the horizontally averaged toroidal field. As is clear, the TSPH and PSPH simulations provide similar results. At first the MRI grows and expels the initial azimuthal fields to the disk corona where strong fields accumulate and are amplified (at about $\sim 60\Omega^{-1}$). Domains dominated by magnetic energy propagate from the corona to the disk midplane and ultimately the entire box is dominated by strong growing azimuthal fields ($\beta\sim 1$). Note that $\alpha_{M}$ is negligible from some 10-20 orbits, indicative that the MRI is quenched. At the end of the simulation the magnetic fields are at equipartition with the gas pressure, and almost entirely azimuthal with no turbulent activity. Simultaneously the disk expands vertically as it becomes magnetically supported. As is obvious, these simulations bear little resemblance to previous grid-code stratified shearing box simulations, which report robust subsonic turbulence in the disk body (Davis et al., 2010; Simon et al., 2011). To test if the strong fields persist in a more numerically dissipative setup, we reran the TSPH simulation with the quartic spline kernel (which has a noisier element distribution) and an isothermal EOS (see figure 9). We find that this TSPH variant in fact damps the turbulence faster and reaches the $\beta\sim 1$ state earlier.

We should note, the result that SPH grows strong toroidal fields is not unique to our simulation setup or code. Dobbs et al. (2016), using the SPHNG SPH code to simulate global galactic disk models, also reported unaccountable growth of magnetic fields. Likewise, similar behavior has also been seen in MHD-SPH simulations of disk formation in tidal disruption events with the code PHANTOM (Bonnerot, private communication). Both of these codes also implement the hyperbolic cleaning method from Tricco & Price (2012) in SPH. Stasyszyn & Elstner (2015) presents a detailed study that discusses the likely numerical issues: they consider 3D, global simulations of a differentially rotating disk with an initially pure-toroidal field, designed so the system is stable and should exhibit no field growth. Using more accurate (CT or vector-potential based) schemes they show that they recover this solution. But using SPH with similar hyperbolic divergence cleaning, they show that discretization error produces small radial field components, which couples to the rotational shear and amplifies this and in turn the toroidal field exponentially. They specifically show that the form of the SPH MHD induction equation leads (in essentially any internally-consistent SPH based cleaning scheme) to the divergence-cleaning amplifying the vertical field, instead of damping the radial field, in order to locally restore $\nabla\cdot{\bf B}=0$.

This demonstrates a few key ingredients that interact here: the particularly virulent form of this instability in SPH requires shear/differential rotation (either in global disk simulations or shearing boxes), non-zero radial, azimuthal, and vertical field components where there is a vertical gradient present that can offset the radial gradient (hence 3D, stratified simulations), and relatively-large $\nabla\cdot{\bf B}$ errors (note these are large here, with $divB\sim 0.01-0.1$).

We should also note that it is possible to construct divergence-cleaning schemes such as that in Tricco & Price (2012) which are total-energy conserving. In highly-idealized test problems, this will serve to limit the non-linear magnitude of any erroneous magnetic field amplification. However, in a shearing box or global thin disk simulation, there is an essentially infinite source of energy from shear, so this does not “rescue” the simulations from excessive numerical dissipation.

More generally, it is well-known that, without any divergence-cleaning, the $\nabla\cdot{\bf B}$ errors are violently numerically unstable: magnetic monopoles grow explosively and the amplitude of ${\bf B}$ is correspondingly rapidly-amplified. It is also well-established that this artificial, explosive field growth can occur even with divergence-cleaning, if the cleaning is not sufficiently accurate, or if it acts “too slowly” to respond to the growth rate. For example, Mocz et al. (2016) showed that using just Powell et al. (1999)-type (considerably less-sophisticated) divergence cleaning in even ordered meshes produces large artificial magnetic field growth (on essentially the Courant timescale) and much larger magnetic field strength, in idealized tests compared to CT methods.

Regarding the damping of turbulence, we should of course note that SPH requires artificial viscosity and resistivity to capture MHD shocks (Cullen & Dehnen, 2010; Tricco & Price, 2013, and references therein). It is well known that SPH tends to over-damp subsonic turbulence due to imperfectly triggered artificial viscosity (Bauer & Springel, 2012; Hopkins, 2015; Deng et al., 2017b). The GIZMO code applies an artificial viscosity switch similar to that described in Cullen & Dehnen (2010) with $\alpha_{min}=0.05$ and $\alpha_{max}=2$ (see Hopkins, 2015, Appendix F2 for details) to suppress unwanted artificial viscosity. This switch works most efficiently in regions away from shocks and may not be effective at regions with large velocity derivatives (Deng et al., 2017a). In our SPH simulations, relatively large artificial viscosity with $\alpha_{sph}>0.2$ is still triggered (see figure 8). Artificial viscosity certainly helps the turbulence dissipate. With damped velocity fluctuations, the MRI and its parasitic modes (Latter et al., 2009; Pessah & Goodman, 2009) cannot grow efficiently. We therefore also explore what happens if we revert to the more dissipative artificial viscosity in GADGET2 (Springel, 2005) and restart the TSPH simulation from $t=50\Omega^{-1}$ (see figure 12); the turbulence does decay faster (as expected) and the strong toroidal fields develop more quickly. Thus, as expected, the MRI turbulence damping owes significantly to the artificial viscosity.

Artificial resistivity dissipates magnetic fields, and a switch to minimize artificial resistivity away from shocks was developed by Tricco & Price (2013). We here apply this artificial resistivity switch with $\alpha_{B,min}=0.005$ and $\alpha_{B,max}=0.1$. We choose this conservative $\alpha_{B,max}$ because the turbulence is subsonic (Hopkins & Raives, 2015). We have re-run the TSPH simulation in figure 7 with $\alpha_{B,max}=1$ (suggested by Tricco & Price, 2013) and obtain similar results. The numerical resistivity in SPH MHD is evidently different from that of Riemann solvers (see appendix B).

4.3 A Transient MRI dynamo in MFM Simulations

Our MFM simulations use the same initial conditions as those of the SPH simulations. We present three simulations. The fiducial model is run with the Wendland C4 kernel and adiabatic EOS. In addition, to test the effect of the kernel function and EOS, we run two simulations with the quartic spline kernel ($N_{ngb}=60$) and with an isothermal EOS. In the isothermal run we solve the energy equation instead of dropping it as done in Stone et al. (2008). To mimic the isothermal EOS we set $\gamma=1.001$ so that the thermal energy dominates the total energy and large truncation errors affect the accuracy of magnetic energy calculation. This can cause the fast dissipation of the magnetic fields.

In figure 9 various quantities are plotted as functions of time. We see here that the MRI grows faster and reaches the magnetic pressure maximum quicker in the simulation with the quartic spline kernel (aqua) compared with the other two simulations. The quartic spline kernel is more compact and has a larger quality factor, $Q^{*}$ than the Wendland C4 kernel (see section 2.3). Although Dehnen & Aly (2012) shows that the quartic spline kernel is superior to the traditional cubic spline kernel it is still vulnerable to pair instability, which introduces numerical noise in gradient estimation. It would appear that this noise provides a significant degree of numerical dissipation because we find that the MRI dies quickly after its initial spike. The isothermal and fiducial adiabatic runs are qualitatively similar: both can sustain MRI turbulence for a period of some 30-40 orbits before dying.

In the first 30 orbits the fiducial model successfully reproduces the quasi-periodic ($\sim$10 orbits) butterfly pattern of the averaged azimuthal fields. Note, however, that the butterfly diagram becomes erratic at $\sim 200\Omega^{-1}$ as seen in other thermal MRI runs (Gressel, 2013; Riols & Latter, 2018). The saturated $\alpha_{M}\sim 0.4$ and $\langle B^{2}/8\pi\rangle/P_{0}\sim 0.01$, however, are both in agreement with previous isothermal grid-codes’ results (Simon et al., 2011). During the simulation, the box expands vertically due to accretion heating leading to a decrease of resolution in the disk body. In figure 10, the density at the disk midplane drops to 0.6 at $200\Omega^{-1}$ which corresponds to 1.2 times larger mean fluid element separation. The decrease in resolution certainly must affect the sustainability of the turbulence. We turn to higher resolution simulations to assess if our results can improve.

4.4 High resolution MFM runs

In order to maintain good resolution over the course of the simulations, we avert expansion resulting from turbulent heat transport and add an ad hoc cooling term as in Noble et al. (2010); Parkin & Bicknell (2013),

[TABLE]

where $\tau_{cool}=2\pi/\Omega$, and $u_{init}$ is the initial specific internal energy constant. This fast cooling maintains the disk scale height nearly constant, thus preserving the initial resolution across the disk. In addition we increase the number of elements to 3 million which results in $\langle Q_{y}\rangle\sim 30,\langle Q_{z}\rangle\sim 10$ in the turbulent state.

In figure 11 various flow properties are plotted. The most important result is that the MRI turbulence is sustained for a longer time (as it should if the method is converging properly). During this phase the main flow diagnostics are in good agreement with those of grid code runs: $\alpha_{M}\sim 0.4$, the averaged magnetic energy is a few percent of the gas pressure, and the Maxwell stress about 4 times larger than the Reynolds stress (Hawley et al., 2011). In figure 11, the butterfly diagram is reproduced but after $300\Omega^{-1}$ the pattern becomes erratic. Comparing to the 1.5M elements MFM fiducial model in figure 9, finer magnetic field structures are captured (see figure 12) and the butterfly diagram/dynamo is better resolved. We note that even this simple cooling can introduce additional numerical noise at the kernel scale (Rice et al., 2014). However, we cannot afford higher resolution, or to run longer simulations, with this set-up (see section 5.1).

5 Discussion

5.1 Computational Cost and Possible Applications

In addition to robustness of numerical results, another worthy metric of comparison between codes is their computational cost to carry out a comparable calculation. We run a setup equivalent to our 3M elements local stratified simulation with ATHENA, and compare the computational cost directly using a fixed timestep in both. We found that the MFM simulation is $>100$ times more computationally expensive than an equivalent run with 32 cells per scale height using ATHENA with the orbital advection method (Masset, 2000; Stone & Gardiner, 2010) for optimization. The lower computationally efficiency of MFM has nothing to do with the hydro solver, rather owes to the neighbor “search tree” which needs to be updated constantly and walked to find neighbors and re-build the domain (because it allows for arbitrary particle re-configuration between timesteps). Of course, in simulations where particle order is not dramatically changing, and the only forces are local, we could in principle save considerable computational expense by simply storing the interacting neighbor lists and re-building the domain less often. Furthermore, there is room to improve significantly the neighbor search algorithm on modern massively parallel architectures coupled with accelerators, as it is being currently investigated for a range of particle-based codes (Guerrera et al., 2018).

Nevertheless, such a difference in performance is highly problem-dependent. The tree algorithm can be efficiently exploited – and the difference in performance is dramatically mitigated – if other physics which involves non-local forces is calculated. A prime example of the latter is self-gravity. The tree-based gravity solver coupled with MFM, indeed, is, generally speaking, both faster and more accurate compared to traditional gravity solvers coupled with grid-based codes, which is the reason why particle-based methods with tree-based gravity have been since long very competitive in comparison with grid-based methods in the modeling of self-gravitating protoplanetary disks (see e.g. Mayer & Gawryszczak, 2008). Furthermore, addressing self-gravitating disks ultimately requires global calculations (Durisen et al., 2007). This is by itself a natural regime for mesh-free codes since one of their major goals is to enable adaptive resolution on global problems, more akin to adaptive-mesh-refinement codes, which have similar computational and memory cost.

The most interesting applications of the mesh-free methods studied here are thus not in idealized MRI setups where accuracy of the MHD calculation over long timescales, absent other physical effects, is the prime objective. Rather, these methods may be more promising for studies of turbulence in magnetized self-gravitating disks, especially the strong dynamo action reported by Riols & Latter (2018). This spiral wave dynamo is vigorous even with large magnetic resistivity and may be responsible for the primordial magnetic field amplification in galaxy formation (Rieder & Teyssier, 2016, 2017), a field where adaptive resolution (either with Lagrangian or AMR-type codes) is essentially required. Likewise, the methods described in this paper have considerable potential for applications in other areas of astrophysics where self-gravitating magnetized disks should be relevant, such as the central regions of massive protogalaxies where self-gravitating circumnuclear gas disks could trigger the formation of supermassive black holes (Regan & Haehnelt, 2009; Choi et al., 2013; Mayer et al., 2010, 2015) or the outer regions of accretion disks around AGNs (Rafikov, 2001). In the case of the protogalactic nuclei, in particular, adaptivity is necessary to capture a wide range of spatial and temporal scales. Furthermore, understanding the interplay between the stabilizing effect of magnetic pressure, turbulence and gas inflows governed by global self-gravitating modes might be the key to understand whether a monolithic central collapse into a supermassive star occurs, which later will turn into a massive black hole, as opposed to fragmentation into stars (Latif et al., 2014). The Riols & Latter dynamo action might play an important role in this latter case as it might reveal itself as an important element to understand the process of angular momentum transport, and thus the evaluate better the possibility of a central monolithic collapse.

5.2 Other Lagrangian MHD methods

We have restricted our study to just two classes of numerical methods, SPH and MFM (although we did consider a few “variants” of SPH). Furthermore, we only considered TSPH and PSPH variants of the SPH method. We should point out that caution is warranted in generalizing any of these results to other Lagrangian methods. Moving meshes or mesh-free finite-volume (MFV)-type methods with divergence-cleaning methods can arbitrarily “smooth” the mesh motion, decreasing the “mesh deformation noise” (McNally et al., 2012; Muñoz et al., 2014) and likely allowing for more accurate divergence cleaning simply because the mesh is deforming less rapidly and less irregularly (so e.g. smaller gradient errors can be ensured). As noted above, unstaggered CT schemes have now been developed (Mocz et al., 2014, 2016) for certain specific types of moving-mesh schemes, which can maintain $\nabla\cdot{\bf B}\approx 0$ at machine precision, so should perform more similarly to CT-grid schemes here, although the numerical noise/dissipation properties of moving-mesh codes (which determine the MRI damping) are often very different.

Fundamentally distinct SPH MHD methods have also been developed. Although early attempts at implementing SPH MHD based on vector potentials did not allow reconnection (e.g. Rosswog & Price, 2007), newer hybrid methods that combine vector potentials with divergence-cleaning in the vector potential space appear to avoid exactly the runaway field amplification discussed here (see Stasyszyn & Elstner, 2015). To our knowledge, however, these schemes have not yet been explored in a broader context or used for MRI simulations. Finally, in the final stages of the preparation of this paper we became aware that a simple variant of an SPH MHD solver based on the GDSPH method in the GASOLINE2 code (Wadsley et al., 2017) is currently being tested in local MRI setups similar to those described here (Robert Wissing et al, private communication). In the latter, the Lorentz force is smoothed in the same way as the hydro force, possibly helping to reduce numerical dissipation.

6 Conclusions

We presented the results of a series of MRI simulations with two meshless MHD methods, SPH and MFM, in both vertically unstratified and stratified boxes. Two variants of SPH were considered, a “vanilla” SPH method based on the density-energy formulation, and PSPH (both as implemented in the GIZMO code). The MRI, especially in its zero-net-flux configuration, is sensitive to numerical or physical dissipation which makes it challenging for any code to simulate, when physical diffusivities are omitted. This is true for both Eulerian and Lagrangian codes, and the results in the zero-net-flux case will always, to some extend, depend on the nature of their numerical dissipation. It is perhaps then no surprise that the biggest discrepancies between the codes are observed for this magnetic configuration.

Our main findings can be summarized as:

1. The use of an appropriate kernel function which does not exhibit the pairing instability and allows a relatively large radius of compact support (e.g. Wendland C4) is crucial for maintaining element or mesh-generating-point order and for accurate gradient calculation. In MFM, this is directly akin to using a larger stencil to obtain more accurate, higher-order gradient estimators in traditional regular-grid codes. Although these kernels are less compact than the traditional spline kernels and tend to over-smooth fluid variables in SPH, they help to sustain the turbulence longer.

2. A stiff adiabatic EOS can help to control the noise in solving the energy equation where the truncation errors can be significant, because the magnetic energy is much smaller than the internal energy.

3. In unstratified shearing boxes with a net vertical field, MFM exhibits a similar error scaling in the linear growth MRI rates compared to the finite volume Eulerian code ATHENA. Both SPH and MFM can adequately simulate the ensuing turbulence, though the former is more diffusive and thus the MRI is closer to criticality. Two consequences of higher diffusivity in SPH are more vigorous channel bursts and very severe heating, as a consequence.

4. In unstratified shearing boxes with zero-net vertical field, SPH and MFM exhibit decaying turbulence at the (relatively low) resolution we are able to simulate here. It is possible this decay is linked to a very low numerical magnetic Prandtl number, but it is more likely that the numerical resistivity is simply too high at this resolution for sustaining the MRI.

5. In vertically stratified shearing box simulations, SPH MHD produces radically unphysical behaviour: turbulence dies out but strong toroidal fields continue to grow to equipartition with the gas pressure. This owes to non-trivial coupling of poorly-controlled magnetic field divergence, differential rotation/shear, and vertical stratification, at least in the most common SPH form of the induction and divergence-cleaning operators.

6. In vertically stratified shearing boxes, high resolution MFM simulations produce results comparable to grid codes implementing the CT cleaning method. For several tens of orbits the classical MRI dynamo is captured, with its characteristic butterfly diagram. Nonetheless, the turbulence ultimately dies out after some 50 orbits, at the relatively low resolution studied as our “baseline” here. Going to higher resolution sustains the dynamo for longer, indicating that the decay is likely due to residual numerical resistivity.

We thank Matthew Bate, Daniel Price, Stephen Rosswog, Jim Stone, James Wadsley, Sijing Shen and Robert Wissing for useful discussions. We acknowledge support from the Swiss National Science Foundation via the NCCR PlanetS. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214.

Appendix A Numerical viscosity

To test the numerical viscosity we perturb the background shear flow by adding a radial dependent azimuthal velocity (a ‘zonal flow’), ie, $\delta{\bm{v}}=Usin(2\pi n_{x}x)\hat{\bm{y}}$. We adjust the internal energy (for the EOS used here, $\gamma=5/3$) so that the pressure perturbation is $\delta P=-\frac{\Omega U}{\pi n_{x}}cos(2\pi n_{x}x)$ and the initial setup is in equilibrium. Numerical viscosity will cause the perturbation to decay. By drawing an analogy to the Navier-Strokes equations (not necessarily true here), $U$ will decay at a rate $\nu_{num}k_{x}^{2}$, where $\nu_{num}$ is the effective numerical viscosity and $k_{x}=2\pi n_{x}$. We fit the decay of $\langle(\delta v_{y})^{2}\rangle$ to determine $\nu_{num}$ and the decay rate of $\langle(\delta v_{y})^{2}\rangle$ is shown in figure 13. We use a shearing box of size $H\times H\times H$ resolved by $32\times 32\times 32$ elements to carry out simulations with $U=0.1$ and $n_{x}=[1,2,3]$. In figure 13, the numerical viscosity in MFM is larger than TSPH when the resolution is low. However, MFM outperforms TSPH at 32 elements per wavelength which is in line with the channel flow growth rate test in section 2.5.

Appendix B Numerical resistivity

Numerical dissipation can destroy magnetic fields. In a periodic box of size $H\times H\times H$ resolved by $32\times 32\times 32$ elements, we initial vertical magnetic fields and adjust the internal energy to set the box in pressure equilibrium. Here we set $\gamma=5/3$ but tests with an isothermal EOS behave similarly. The fields take the form of $\bm{B}=B_{0}\hat{\bm{z}}sin(2\pi n_{x}x)$, where $B_{0}=\sqrt{8\pi P_{0}/\beta}$. In MFM simulations, the field structure decays due to numerical resistivity (see figure 14). $B_{0}$ should decay at a rate $\eta_{num}k_{x}^{2}$. The decay rate can be determined by fitting the decay of the averaged magnetic energy (decays twice as fast as $B_{0}$). We vary the wavelength ($n_{x}$) and field strength ($\beta$) to test how strong the numerical dissipation is. In this test, the dimensionless divergence of the magnetic fields is $\sim$0.0001 and we believe the dissipation is not caused by the non-zero divergence. However, in SPH simulations the magnetic energy doesn’t decay even with $n_{x}=3$ (see figure 15). The numerical noise from the artificial viscosity and resistivity break the perfect lattice.