Mesh-Free Anisotropic Diffusion

Anisotropic Diffusion in Mesh-Free Numerical Magnetohydrodynamics

Philip F. Hopkins
TAPIR & The Walter Burke Institute for Theoretical Physics, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
Submitted to MNRAS, January, 2016

We extend recently-developed mesh-free Lagrangian methods for numerical magnetohydrodynamics (MHD) to arbitrary anisotropic diffusion equations, including: passive scalar diffusion, Spitzer-Braginskii conduction and viscosity, cosmic ray diffusion/streaming, anisotropic radiation transport, non-ideal MHD (Ohmic resistivity, ambipolar diffusion, the Hall effect), and turbulent “eddy diffusion.” We study these as implemented in the code GIZMO for both new meshless finite-volume Godunov schemes (MFM/MFV). We show the MFM/MFV methods are accurate and stable even with noisy fields and irregular particle arrangements, and recover the correct behavior even in arbitrarily anisotropic cases. They are competitive with state-of-the-art AMR/moving-mesh methods, and can correctly treat anisotropic diffusion-driven instabilities (e.g. the MTI and HBI, Hall MRI). We also develop a new scheme for stabilizing anisotropic tensor-valued fluxes with high-order gradient estimators and non-linear flux limiters, which is trivially generalized to AMR/moving-mesh codes. We also present applications of some of these improvements for SPH, in the form of a new integral-Godunov SPH formulation that adopts a moving-least squares gradient estimator and introduces a flux-limited Riemann problem between particles.

methods: numerical — hydrodynamics — MHD — instabilities — conduction — diffusion

1 Introduction

Diffusion operators are ubiquitous in physical systems, manifest in e.g. conduction, viscosity, cosmic ray diffusion/streaming, photon diffusion (in optically thick media), sub-grid “turbulent eddy” diffusion, passive scalar (e.g. metal) diffusion, Ohmic resistivity, ambipolar diffusion, the Hall effect, and more. Most of these are fundamentally anisotropic. Magnetic fields, fluid flow, or radiation flux can all break isotropy and introduce a (local) preferred direction. The effects are large – for the same initial conditions but different magnetic field configurations, effects like viscosity or conduction might dominate the hydrodynamics, or be completely suppressed. Moreover, anisotropy can introduce qualitatively new behaviors and physical instabilities (e.g. Balbus, 2000; Quataert, 2008).

Clearly, it is important to capture these phenomena accurately in numerical simulations. In regular, non-moving mesh methods (including fixed-grid and adaptive mesh-refinement [AMR]), the properties of diffusion equations are reasonably well-understood, but even in these cases it is highly non-trivial to formulate anisotropic operators in a numerically stable fashion (see, e.g. Sharma & Hammett, 2007). But for some classes of problems, these methods are sub-optimal. They tend to produce excessive diffusion when a fluid is moving rapidly relative to the grid, especially across contact discontinuities (Tasker & Bryan, 2008; Springel, 2010; Hopkins, 2015); have difficulty coupling to N-body gravity methods and handling self-gravitating hydrostatic equilibrium (Müller & Steinmetz, 1995; LeVeque, 1998; Zingale et al., 2002); introduce low-order errors around refinement boundaries (O’Shea et al., 2005; Heitmann et al., 2008); and feature inherently preferred directions which can introduce systematic errors even at high resolution if physical anisotropies are not aligned with the grid axes (e.g. Peery & Imlay, 1988; Hahn et al., 2010; Hopkins, 2015; Hopkins & Raives, 2016).

Mesh-free methods can avoid these sources of error, and so are popular for many problems in astrophysics. However, the most popular mesh-free method, smoothed-particle hydrodynamics (SPH) has its own difficulties. For example, SPH requires the use of artificial diffusion operators to maintain numerical stability, which can produce greater numerical viscosity and over-smoothing of shocks and discontinuities compared to Riemann-based shock solvers (Cullen & Dehnen, 2010); it tends to suppress fluid-mixing instabilities (Morris, 1996; Ritchie & Thomas, 2001; Agertz et al., 2007); and suffers from zeroth-order errors (“E0” errors; Morris, 1996; Dilts, 1999; Read et al., 2010) which produce systematic errors that do not converge with resolution alone. Although there have been dramatic improvements in all of these (see Hopkins, 2013; Rosswog, 2014; Hu et al., 2014), the E0 errors cannot be completely eliminated from SPH without making the method numerically unstable (Price, 2012b).

Recently, however, Lanson & Vila (2008a, b); Gaburov & Nitadori (2011); Hopkins (2015) have developed a class of new, mesh-free Lagrangian finite-volume methods which are both high-order consistent and fully conservative. Similar to moving-mesh codes (Springel, 2010; Duffell & MacFadyen, 2011; Gaburov et al., 2012), these new methods appear to capture many of the advantages of AMR and SPH, while avoiding the disadvantages above (although, of course, they feature their own sources of error such as enhanced “grid noise” and volume “partition noise”; see Hopkins 2015). In Hopkins 2015 (hereafter Paper I) and Hopkins & Raives (2016) (Paper II), these are developed for magnetohydrodynamics (MHD) in the multi-method magnetohydrodynamics+gravity code GIZMO,111A public version of the code is available at built on the -body gravity and domain decomposition algorithms from GADGET-3 (Springel, 2005).

In this paper, therefore, we extend these new methods to include arbitrary anisotropic diffusion operators. We consider a systematic comparison of tests in order to determine the degree of anisotropy that can be reliably treated. In the Appendices, we show how some parts of the new methods here can also be applied to SPH.

2 General Diffusion Operators in Meshless Finite-Volume Godunov Methods

The general diffusion equation in conservative form is:


where is a Lagrangian derivative, and a diffusive flux given by a linear combination of the tensor field and the gradients of the field . We denote the inner product “” and outer product “.”

In the mesh-free finite-volume or finite-mass (MFV or MFM, respectively) formulations of the Lagrangian methods in Gaburov & Nitadori (2011) and Paper I, we begin from a general conservation equation of the form and use this to derive a second-order accurate Godunov-type numerical expression of the hydrodynamic equations. It is therefore trivial to apply the same here, giving


where is the “effective face area” (defined in Paper I; it depends on the inter-element spacing and kernel shape, and reduces to the geometric faces of a Voronoi tesselation in the limit of a delta-function kernel). The same would obtain for Cartesian grid methods and moving-mesh codes, with the usual inter-face area. Here is the interface value of the flux.

The solution then follows our usual MHD method: we calculate the coefficients, perform a reconstruction step to estimate quantities “at the face,” replace the flux with the solution to a Riemann problem (RP), and use this to update the conserved quantities. We will explain this in more detail below.

Note that Eq. 3 is manifestly antisymmetric, so conserved quantities are manifestly conserved as desired.

Additionally, we note that the MFM and MFV methods from Paper I differ only at second order in how they handle advection between particles. In diffusion-only problems, they are manifestly identical. Therefore, we do not show both (although we have confirmed their identical results), but will, for simplicity, adopt MFM as our reference method.

2.1 Gradient Estimation & Reconstruction

We require gradients for all quantities; to do this, we adopt the standard gradient estimator in GIZMO, a moving, second-order accurate and consistent least-squares estimator. For a scalar , this is


here we assume an Einstein summation convention over the indices corresponding to the spatial dimensions, and is an (arbitrary) weight function defined in Paper I. This estimator is second-order accurate for an arbitrary mesh/particle configuration, minimizes the (weighted) least-squares deviation , and has been applied in a wide range of different numerical methods (see e.g. Oñate et al., 1996; Kuhnert, 2003; Maron & Howes, 2003; Luo et al., 2008; Lanson & Vila, 2008a). Note that Eq. 4 applies to a scalar field . For a general tensor , we apply Eq. 4 separately to every component to determine all partial derivatives needed to construct and . The gradients are slope-limited as described in Paper I, such that they do not create new local extrema within the interacting kernel.

In the reconstruction step we must extrapolate the values from and to the left and right sides of the face. For hydrodynamics and MHD, we perform a linear reconstruction of all the MHD quantities based on their slope-limited gradients (see Paper I). For additional quantities needed for diffusion, our default approach is a first-order reconstruction, e.g. . This is easy to implement, and most stable, but comes at the cost of greater numerical diffusion. We have therefore also considered limited tests using the “double linear” reconstruction of Muñoz et al. (2013). This amounts to treating like any other primitive variable. In a first loop, we calculate with our standard method; in a second loop, we calculate , and use this to linearly reconstruct (or any component of it) like any other primitive variable. This gives a simple, second-order Taylor-series representation which implicitly includes an appropriately larger stencil (since the second pass sums over the values , themselves constructed from the in all neighbors of ). To ensure smoothness, the reconstructed gradients in the double-linear method are slope-limited with respect to the particle-centered gradients using the same limiter we adopt for all the usual MHD quantities (see Paper I).

2.2 The Riemann Problem

We treat diffusion in operator-split fashion from the pure-MHD RP. We have experimented with several Riemann solvers for the diffusion RP, and find the best compromise between numerical diffusion, stability, and flexibility using the Harten et al. (1983) (HLL) solver in the Lagrangian frame. In this case:


where the maximum/minimum wavespeeds / are determined appropriate to the problem. Here , are the appropriate right/left states reconstructed at the face following Paper I (extrapolated from the particle-centered values to the face with the derivatives of , with a MINMOD slope limiter).

Consider the case where we are using our second-order reconstruction in a Lagrangian frame. We adopt the wavespeed estimate from Paper I, , . Here is a fastest wavespeed determined by the 1D RP. Then, because in our default GIZMO method our frame is moving with , in the Lagrangian frame we simply obtain , in which case the HLL solution reduces to the Global Lax Friedrich (GLF) function: .

Although using the GLF flux in non-Lagrangian methods tends to be excessively diffusive, here it gives nearly indistinguishable results from using the HLL solution with another wavespeed estimate (e.g. Roe wavespeeds or the exact eigenvalues of the Jacobian ). This is because the first-order term owing to frame motion is automatically accounted for by the Lagrangian nature of the method. Even if we use a double-linear reconstruction, we see only a tiny (percent-level) increase in diffusion using this particular form for the RP solution.

Using the above with gives an effective numerical diffusivity , which can easily exceed the physical by large factors at low resolution. We prevent this by replacing Eq. 6 with the flux-limited equation:


where above is the flux we would obtain in the simplest 2nd-order accurate (but numerically unstable) formulation. is the diffusive flux from the Riemann problem, appropriately limited by the function . Here refers to the Frobenius norm of the tensor , and refers to the interface value of the primitive variable ; where possible we adopt these from the appropriate solution to the pure-MHD RP, otherwise we approximate .

The first term in (in and ) limits for the anisotropic case, vanishing as does where there is full anisotropic suppression (even at low resolution). Note that the form of in Eq. 10 follows directly from the form of . For the cases where is represented by a different linear combination of gradients, the term in Eq. 10 should be modified accordingly.

The second term in ensures at small and at large (the functional form is motivated by Rosdahl & Teyssier, 2015); together with the MINMOD application this prevents the numerical diffusivity from exceeding physical values by more than some tolerance parameter . As usual this parameter represents some tradeoff: increasing gives smoother solutions at the expense of numerical diffusivity. We find our qualitative results are robust for all , and use as our default.

Finally, we also compute a “direct” flux based on pairwise direct-difference gradients, and use this to restrict via:


with the tolerance parameter . We find stable (albeit slightly more noisy) results in all our problems using values as large as , in fact, and larger values do give improved performance at low resolution on multi-dimensional tests (e.g. the diffusing ring), but we adopt as our default here.

We note that this is not the only way to stabilize the anisotropic diffusion equations. For example, Sharma & Hammett (2007) propose an elegant slope-limiting method; however, it is not obvious how to extend this to unstructured meshes. Recently, Kannan et al. (2015) implemented the method of Gao & Wu (2013) for moving-mesh codes. The advantage of their method is that it is extremum-preserving and generalizes relatively easily to implicit integrators. However, we choose to explore alternatives for three reasons. First, the Gao & Wu (2013) method implicitly uses the lower-order “direct” gradients, as opposed to our (in principle arbitrarily high-order) matrix-based gradient estimator, necessarily making the method lower-order. Moreover in flows where there is a clear mean-field gradient but large resolution-scale noise, matrix-based estimators are significantly more robust (see García-Senz et al., 2012; Mocz et al., 2014; Pakmor et al., 2016). Second, our method here allows for non-linear flux-limiting terms (e.g.  above) which allow us to limit numerical diffusion to physical values even at arbitrarily low resolution (potentially important in multi-physics problems where the diffusion may, in some places, dominate only on un-resolved scales). And third, our method here, unlike most in the literature, allows for any (arbitrarily complex) tensor , and/or linear combinations of and the elements of in the fluxes (relevant for e.g. radiation transport, Braginskii viscosity, and the Hall effect).

2.3 Timestepping

In addition to the usual timestep limiters (e.g., the CFL condition, gravitational acceleration-based limiters) which always apply, ensuring numerical stability in explicit methods for diffusion equations requires an additional timestep criterion:


For the simplest diffusion example, and an appropriately slope-limited gradient, this reduces to the common expression . We also require that the inter-particle flux between a particle and any neighbors of the conserved cannot exceed half the minimum of the in the pair, in a single timestep, but find this criterion is always satisfied if the above (more strict) criterion is as well.

3 Examples Implemented

We implement this method in the code GIZMO. A number of specific physical cases have been implemented; our focus here is not on the microphysics but on the numerical methods. Still, it is useful to list the relevant examples, both to explicitly see how they correspond to the general formulation above, and to give examples of the different forms of anisotropy that are physically expected.

3.0.1 Isotropic Passive-Scalar Diffusion

For this simple case, we have


where is the (arbitrary) diffusion coefficient, the identity matrix, and the scalar density.

3.0.2 Anisotropic Thermal Conduction

Isotropic and anisotropic thermal conduction are represented by:


where is the gas temperature, the gas density and the internal energy per unit mass, the isotropic diffusion coefficient, the anisotropic (parallel) diffusion coefficient, and the direction of the local magnetic field. Here we allow either arbitrarily specified , or default to the Spitzer conductivity (for , with ) with a limiter for cases where the implied heat transport rate exceeds the free electron streaming rate.

3.0.3 Anisotropic Cosmic Ray Diffusion & Streaming

Cosmic rays diffuse along magnetic field lines. Moreover, as shown by Uhlig et al. (2012), cosmic-ray streaming can be represented numerically via a diffusion operator. Therefore we have, for the simple case of a single-species cosmic-ray population


where is the cosmic ray pressure (proportional for a single-species model to the CR number density), is the cosmic ray energy density, and and are the appropriate “effective” coefficients for diffusion and streaming processes, respectively.

3.0.4 Turbulent Eddy Diffusion Models

In popular “sub-grid” models for turbulence, the effect of unresolved eddies is treated as a diffusion process. Following the Smagorinsky (1963) model, we model this for e.g. diffusion of a scalar field via


where is the resolution scale (for our meshless methods, defined to be equal to the rms inter-element spacing inside the kernel), is a constant calibrated to numerical simulations in Smagorinsky (1963), and is the symmetric shear tensor.

3.0.5 Anisotropic Radiation Transport in the Diffusion Limit

In the optically thick limit, the radiative transfer moment equations are commonly expressed as a diffusion equation. It is convenient to represent this in the form


where is the speed of light, the opacity at frequency , the photon number (or energy) density, an optional “radiative flux limiter,” and the dimensionless Eddington tensor. Note the instead of has no effect on our methodology, it simple re-orders the gradient terms in . Various numerical methods (e.g. flux-limited diffusion and “variable Eddington tensor” methods) rely on this description.

3.0.6 Anisotropic Viscosity

For a viscous fluid, we have


where is the velocity and the density. In this case the form of depends on the viscous parameterization. In the case of a magnetized fluid, the Braginskii viscosity can be written


while for the un-magnetized case, it is common practice to decompose the viscosity into shear () and bulk () terms, following


Here “” denotes the double-dot-product. As above, the particular arrangement of and within the double-dot-product (for the magnetized case) or decomposition of into shear/bulk terms have no effect on our methodology; they simply re-order the gradients of within . As with conduction we allow either freely-specified , , , or can calculate the coefficients accoding to Spitzer-Braginskii theory with the appropriate limiters. We also add the corresponding viscous term to the energy equation.

3.0.7 Non-Ideal MHD

Astrophysical non-ideal MHD effects are typically parameterized as Ohmic dissipation (controlled by the resistivity ), the Hall effect () and ambipolar diffusion (). All appear as diffusion operators in the induction equation; if we operator-split the ideal MHD term (already solved in GIZMO), we have


where . We can, with some elaborate algebra, write this as an equation in , but it is easier to cast this directly into the alternative Godunov form:


where , , are the components of and the coefficients are calculated from the plasma properties. The appropriate magnetic terms are added to the energy equation.

Figure 1: Diffusing sheet test (§ 4.1). A 3D box ( elements across plotted -range) is initialized with a step-function discontinuity () in diffusing quantity . The solution is scale-free in units shown, and independent of the physics (we confirm identical behavior for implementations of passive-scale diffusion, conduction, cosmic rays, turbulent eddy diffusion, radiation transport in the diffusion limit, viscosity, and Ohmic resistivity; § 3). We consider isotropic and anisotropic cases with diffusion tensor , with parallel, partially-aligned, and perpendicular fields . In each case we compare our mesh-free finite-element methods (MFM, our MFV method is manifestly identical on these problems; § 2) to exact solutions.

4 Numerical Tests

4.1 Diffusing Sheet

For a first test in Fig. 1, we study a simple discontinuity in one dimension, for , for , with all other properties fixed. To focus on the diffusion equations, we turn off all MHD forces (for now). We also simulate this in a 3D box (with an arbitrarily large periodic domain), because although it is essentially a 1D problem, sensitivity to particle arrangement makes accurate solution much more challenging in 2D or 3D. Finally, we also add noise (which is critical to test for numerical instability): we add a uniform random value to between , for all particles.

4.1.1 Isotropic Case

In the isotropic case, after the noise has been (quickly) damped away,222We have considered both a perfect step-function initial condition, and (alternatively) initializing the profile with the analytic solution corresponding to time such that . These are indistinguishable after a short time (). this has the trivial analytic solution


The solution is completely self-similar so the absolute values of all quantities are irrelevant; the only numerically important quantity is the number of resolution elements over which the contact discontinuity has been diffused (as we evolve the solution, it becomes better-resolved). We therefore plot results at fixed .

Our results converge well to the analytic solution with both methods even at low resolution (). Because the initial noise is grid-scale, it should be damped away on a small timescale ; we confirm this.

Note that the solution here is independent of the physics; we have explicitly verified that our implementations of passive-scalar diffusion ( above), conduction ( with ), cosmic rays (, with ), eddy diffusion (), isotropic radiation transport ( with ), viscosity ( with , , and ), and Ohmic resistivity () all give identical results on this test (given the same diffusivity), as they should.

4.1.2 Anisotropic Case

Next we consider the anisotropic case. Here we take with and constant and . In this case, the solution is identical to the isotropic case but with , so it is entirely specified by the absolute value of the projection of in the gradient () direction, .

We have considered values between to check for pathological behavior; these are summarized with three representative cases shown here: (perpendicular fields, which should completely suppress diffusion), (parallel fields; the solution should be identical to the isotropic case), and . For all scalar diffusion cases (passive scalars, cosmic rays, eddy diffusion, conduction) this gives identical results. For our radiation diffusion, we make the system anisotropic by instead taking , with constant. We confirm this produces exactly identical solutions to the cases above.

Our MFM method is able to handle all three cases accurately; we confirm that there is zero diffusivity for the perpendicular case (up to machine precision, if our initial depends on alone) and that the parallel case exactly matches the isotropic case. Convergence is again good even at low resolution , and (we show below) the results are insensitive to noise or particle order.

Figure 2: Anisotropic sheet (as Fig. 1) with varied ICs (§ 4.1.3). We compare (1) A “default” case with resolution particles across the plotted range in the -direction ( particles across the range of diffusion at this time, in MFM), and effective neighbor number within the kernel partition function. Particles are initially in a cubic lattice, with noise in equal to of the jump value. (2) High-resolution () and neighbor number (, with a quintic kernel function). (3) Initial noise of jump value. (4) Initial particles laid down randomly (according to a Poisson distribution). Our MFM results are indistinguishable from each other and the exact solution in each case.

4.1.3 Dependence on Particle Order, Resolution, Noise, and Neighbor Number

Fig. 2 considers variations of the anisotropic sheet. In order to test whether our methods are sensitive to the local arrangement of particles, we have considered (in both 2D and 3D tests) an initial particle distribution following (1) a regular square lattice, (2) uniformly randomly-distributed particles over the volume (a Poisson distribution), (3) a glass (generated from the random distribution), and (4) a densest sphere packing. In 2D we have also considered a regular triangular and hexagonal grid. For our MFM/MFV method, the results are almost indistinguishable (after the initial noise is damped) in every one of these cases (even the “worst-case” Poisson distribution, shown in Fig. 2), clearly demonstrating that the method does not depend sensitively on particle order.

We consider resolution tests, varying both the absolute resolution and number of kernel neighbors. Our MFM/MFV results are well-converged even with just resolution elements across the jump (i.e. ). The convergence rate of our method is also independent of neighbor number, as shown in Paper I and Paper II (going to larger neighbor number simply trades a reduction in noise for increased numerical diffusivity).

We have also tested for sensitivity to initial noise, varying the seed noise level from of the jump level. In all cases our results are stable and do not qualitatively change with respect to this. Likewise, the steepness of the initial discontinuity has no effect on our results (whether we begin with the full solution at a resolved scale, or a perfectly steep discontinuity across a single particle).

Figure 3: Anisotropic diffusing sheet problem, as Fig. 1 (), but here focusing on two special cases where the solution is not the same as in Fig. 1 (see § 4.1.4). There is no analytic solution for the cases shown, so, we compare a converged solution from ATHENA. Left: Braginskii viscosity (§ 3.0.6). The initial discontinuity is in only, but this form of viscosity generates at later times. We consider two cases: (upper; this should produce non-zero diffusion) and (lower; this should have vanishing diffusion). Right: Hall effect (§ 3.0.7). The initial discontinuity is in ( constant, vanishes). We consider a case which should diffuse (; upper) and one which should be suppressed (; lower). MFM and ATHENA agree well in all cases.

4.1.4 Braginskii Viscosity and the Hall Effect

Two cases in § 3 are more complicated, even in a simple diffusing sheet setup. These are Braginskii viscosity (§ 3.0.6) and the Hall effect (§ 3.0.7). In both, is a vector ( or , respectively), and even if we initialize a gradient only in one element of that vector and set the other elements to vanish everywhere, the form of the diffusion operator leads to growth of other components of . Mathematically, the anisotropic part of the diffusion operator in most cases is a projection operator; here, it also includes a rotation operator. This leads to different non-linear solutions and makes it more challenging to achieve stability.

Therefore we consider these cases specifically in Fig. 3. Since the non-linear solutions do not have analytic forms even for this simple test, we compare to a converged, high-resolution ( across the plotted domain), one-dimensional solution from the well-tested grid-based code ATHENA (Stone et al., 2008). We run ATHENA in its most accurate constrained-transport, PPM-CTU (highest-order) mode. For simplicity, we do not disable the other hydrodynamic forces, but these are not dominant. To avoid certain boundary condition effects in ATHENA we initialize the test problem in slightly different fashion: we take


For the Braginskii problem, we take initial and consider both (which should produce zero diffusion) and (which should diffuse); we take , and to be small so it has no effect on the dynamics except to control the anisotropy. For the Hall problem, we take vanishing initial velocities and (no diffusion) (diffusion). Our MFM simulations are run in 3D boxes with resolution across the plotted domain.

Fig. 3 shows these develop non-zero and , despite these vector components initially vanishing. In both cases our MFM/MFV methods produce solutions in excellent agreement with ATHENA, even well into non-linear evolution. Both cases with diffusion, and cases with full anisotropic suppression, are captured.

Figure 4: Diffusion across a super-sonically moving, rotated, shearing contact discontinuity (§ 4.2), as Fig. 1, with full MHD physics enabled. Top: Illustration of the initial problem setup. We measure properties along the axis, perpendicular to the initial discontinuity in density and . We consider an isotropic case (middle) and an anisotropic, perpendicular (fully-suppressed; bottom) case. The contact discontinuity means the particle arrangement is different on either side of , and the shear field means that the arrangement of neighbor particles is constantly changing. These do not destroy our solution.
Figure 5: Diffusion of a sinusoidal pulse (; § 4.3), in a 1D periodic box with with elements across the -range. We show anisotropic cases with perpendicular (full suppression) and parallel (no suppression) fields. In the perpendicular case, the initial pulse does not evolve; in the parallel case, the amplitude of the pulse decays as (we show results at ).

4.2 Diffusion Across a Moving, Rotated, Shearing Contact Discontinuity

We next consider a more challenging variation of the diffusion sheet, illustrated in Fig. 4. We (1) re-enable our normal MHD physics. (2) Insert a contact discontinuity at (the location of the jump in ) with a density jump of a factor of , so the diffusion (of some passive scalar) is across the discontinuity. This implies a different particle arrangement (since our particles are equal-mass) across . (3) Make the discontinuity shearing, with . This means the particle geometry around the diffusing interface is constantly being re-arranged. (4) Uniformly boost the system by . (5) Rotate the entire system by .

None of these changes the physical solution for the diffusion. However they are all, in principle, numerically challenging. Because our methods are fully Lagrangian, our solutions are trivially invariant to the boost (4) and rotation (5) operations. This is not the case, however, for Eulerian methods. Our MFM method is also invariant to (1) and (2), i.e. the evolution of a stable contact discontinuity produces vanishing fluxes, and the gradient estimator is explicitly insensitive to particle arrangement within the kernel for linear gradients (this is not always the case in other Lagrangian methods such as SPH, where solutions around contact discontinuities can depend on particle arrangement; see references in § 1).

Fig. 4 shows that, despite these complications, our MFM solution still agrees very well with the exact result. There is some noise around , introduced by the shear (3), particularly around the contact discontinuity, but it is not visible on the scale plotted. This constant re-arrangement of the particles (hence constant re-arrangement of the effective faces and implicit mesh) introduces “grid noise” (for detailed discussion of this noise term see Duffell & MacFadyen, 2014; Hopkins, 2015; Mocz et al., 2015; Pakmor et al., 2016).

A detailed comparison with grid-based methods is outside the scope of this paper; however, this problem is very challenging for such codes. A moving contact discontinuity, especially one not aligned with the grid, produces large numerical diffusivity at low resolution (see Hopkins, 2015; Springel, 2010). Using ATHENA to run a version of this problem with anisotropic conduction, we find the numerical diffusivity exceeds the physical diffusion even in the isotropic case unless we use more than elements across the plotted range (compared to here); even higher resolution is required for the anisotropic case.

Figure 6: Diffusion with a variable diffusivity and diffused quantity (§ 4.4). Both vary linearly in . The box is 3D with elements across the -range, but the results are converged (up to tested ). We consider anisotropic cases with parallel (left; no suppression) and perpendicular (right; full suppression) fields.

4.3 Sinusoidal Temperature Distribution

Fig. 5 considers a one-dimensional test problem in which a scalar follows a sinusoidal distribution from Arth et al. (2014). This tests the same physics as the later stages of the diffusing sheet, but is much “easier” (since it is 1D and there is never a steep gradient). No other MHD physics are active, and we take a periodic box of unit length , with unit density and sound speed and , with the physical solution


initialized at . We have confirmed (as expected) that all the conclusions from our diffusing sheet tests apply in this test.

4.4 Variable Diffusivity

We previously took to be constant. Here, consider the case , ; this produces the analytic solution . In the simple anisotropic case with this becomes . We consider this in the same 3D box setup as before (here with large enough distance in the -direction so that the boundary conditions do not enter the considered domain).

The results are shown in Fig. 6. Our methods perform well in both the anisotropic case and (more trivially, therefore not shown), the isotropic case.

Figure 7: Three-dimensional Gaussian pulse test (§ 4.5). A Gaussian distribution of the diffused quantity is injected at the center with initial spherically symmetric width . We plot a slice through the plane, at time for diffusivity , in a box of unit size with elements. Color denotes , scaled linearly, in arbitrary units (so we simply set ). Left: Isotropic diffusion. The Gaussian simply expands; all methods accurately capture the exact solution (bottom). Right: Anisotropic diffusion with and . The distribution should expand only in the -direction. MFM captures this, with minor artifacts from the initial particle arrangement which converge away at higher resolution.

4.5 Gaussian Pulse

We now consider a multi-dimensional problem – the diffusion of a quantity injected as a -function instantaneously into a homogenous background. In a periodic box of unit size, we initialize a 3D, spherically-symmetric Gaussian for centered on the origin, where the diffused quantity is treated as a passive scalar with all other background properties constant (so no other MHD effects appear). In the isotropic case, this evolves as:


where is an arbitrary normalization, is the diffusivity, the time since the problem was initialized, and defines the initial width of the distribution ( becomes a -function; larger correspond to starting from an already-evolved solution). We take , comparable to our inter-particle spacing, so that there is a well-defined gradient in our initial condition.

In the anisotropic case, if we assume with constant (in space and time) , we can always define our axes so that ; then this evolves as:

i.e.  diffuses normally along , and not perpendicular.

Fig. 7 shows the results. The isotropic case is easily recovered accurately (even at lower resolution than shown). MFM is also able to recover the anisotropic case; however at low resolution (a box) the agreement is not perfect, as some artifacts from the grid structure (here particles were initially laid in a Cartesian grid) are present. These are invisible by-eye if we go to resolution, although still measurable in the L1 error norm. In multi-dimensional problems such as this, the perpendicular width of structures must, in general, be a few particles across before complete anisotropy can be fully captured. This is required so that a reliable gradient in the relevant direction can be determined (similar to the requirement in grid-based codes).

Figure 8: The diffusing ring test problem (§ 4.6; style as Fig. 7). A “hot spot” (high-) within a small radial annulus is initialized with pure azimuthal fields; this should diffuse into a ring without additional radial diffusion. We compare MFM at two different resolutions (labeled), and both early (top) and late (bottom) times. For early times there is an approximate analytic solution, shown; for late times there is no analytic solution, but the high- material should remained confined to the same radial annulus as at early times, gradually diffusing around the ring until it is isothermal. In MFM/MFV we see good agreement with the expected behaviors at both times, even at low resolution. There is some numerical radial diffusion, but this gradually converges away as we go towards higher resolution.
Figure 9: Convergence in the diffusing ring test (Fig. 8), with our MFM method. We plot the L1 error norm in averaged over the domain as a function of the number of elements on a side , at a fixed time . Convergence is close to the ideal linear scaling (dotted red).
Figure 10: Diffusing ring test for radiative diffusion (§ 4.6.2; style as Fig. 7), where the diffused quantity is an anisotropic tensor, with (i.e. a purely azimuthal Eddington tensor). We initialize a small (), thin (), azimuthally symmetric ring (in a box of unit length); the ring should expand with speed and remain thin. The resolution here is , and we show the result at early (; top) and late (; bottom) times. MFM captures the physically correct behavior.

4.6 Diffusing Ring

A more challenging version of this problem is diffusion with azimuthal anisotropy, following Parrish & Stone (2005); Sharma & Hammett (2007, 2011). In a periodic box of unit size centered on the origin and cylindrical coordiates, , , and purely azimuthal magnetic fields , we initialize . Here and are an arbitrary background and normalization (we take and ), and define a Gaussian ring at radius of width , and is an initial Gaussian spread about in the direction. We assume .

Because the fields are purely azimuthal, the quantity should diffuse in the purely azimuthal direction, “around the ring,” rather than in the radial direction. At early times, this has an exact solution of the same functional form: where (with normalization ). This assumes we can neglect the periodic boundary conditions around the ring – i.e. is valid for (hence early times). At late times, the diffusion from both directions self-intersects on the opposite side of the ring, and there is no simple exact solution. Eventually, though, as , the system becomes isothermal within each azimuthal annulus. This is a challenging problem even in high-order grid codes (see Parrish & Stone, 2005).

Fig. 8 compares the results at early and late times, at two different resolution levels. As before, MFM is able to capture the azimuthal anisotropy. Even at low resolution (), there are only weak grid artifacts, but these and the amount of perpendicular diffusion improve at higher resolution. At late times, for comparison, in the fixed-grid code ATHENA on a Cartesian mesh (where the preferred direction of the grid is not the azimuthal direction), it requires going to resolution before the diffusion properly “wraps” into a ring at all (see e.g. Sharma & Hammett, 2011); our case resembles a case with ATHENA. And we note that we have not aligned the particles with the anisotropy (the particles are in a regular triangular lattice).

Note that we have tested both the 2D ring version of this problem and the 3D version, where the ring becomes a cylinder elongated in the direction and the box is periodic in that direction. The results are very similar in both cases.

4.6.1 Convergence

Note that our diffusing ring setup is slightly different from that in Sharma & Hammett (2007, 2011) and Kannan et al. (2015), who used this problem to measure the convergence properties of their method. We have therefore also compared initial conditions which match their choice: we initialize a step-function discontinuity with and inside or outside (respectively) of an annulus and . The qualitative results with all methods are identical to those shown in Fig. 8.

Fig. 9 quantifies the convergence (in 2D tests using this initial condition) by measuring the L1 norm of at time relative to a high-resolution solution () interpolated to the particle positions. We find a convergence rate , close to the ideal . This is competitive with and in some cases superior to the implementations studied in fixed-grid and moving-mesh codes in Sharma & Hammett (2007, 2011); Kannan et al. (2015). However we caution that we are comparing to our own high-resolution solution, not an exact solution (because none exists); so systematic errors which may converge more slowly do not appear.

4.6.2 Tensor Diffusion: The Radiative Diffusion Case

Interestingly, if we consider the radiative diffusion version of this problem, the behavior is qualitative different. Take and (define ); because the anisotropy is inside the first gradient in the diffusion equation (i.e. we have as opposed to ), the solution is distinct. For an azimuthally symmetric and , the diffusion equation for scalar reduces to , i.e. the ring expands radially with a speed . Fig. 10 shows the results of this test. In our MFM/MFV methods, the ring expands as expected, and the correct qualitative behavior is captured even at extremely low resolution ().

Figure 11: Magneto-thermal instability (MTI) and heat-flux driven bouyancy instability (HBI) tests (§ 4.7), with our MFM method. Top: MTI: The initial condition (left) is a 2D box of length , with resolution elements, with a vertically stratified atmosphere with (temperature decreases upwards), in equilibrium with a constant vertical gravitational field and conductivity . Magnetic field lines are shown at different times ; these are initialized with trace values aligned along . Color denotes , increasing linearly from the minimum to maximum values at each time (cyan to purple to magenta, respectively; see Fig. 13 for quantitative values). With no diffusion, or isotropic diffusion (; center), the system is stable and the field configuration should be preserved. With anisotropic diffusion (; right), the system is unstable and develops convection. This re-orients the field to be initially near vertical (at e.g. the time shown), then breaks up into turbulence and the field becomes isotropic. The characteristic timescale is the bouyancy time in these units. Bottom: HBI. The resolution and conductivity are the same but the initial atmosphere now has , and initial . Again with no diffusion or isotropic diffusion, the system is stable. With anisotropic diffusion, the HBI amplifies transverse motions that re-orient the field to be perpendicular (). MFM recovers all of the expected behaviors for both instabilities with anisotropic diffusion, and correctly suppresses them with isotropic diffusion.
Figure 12: Development of the MTI in an MFM simulation at low () resolution, now with a box of side-length . Field lines are shown as Fig. 11. Initially small seed velocity perturbations grow in the vertical direction and become the convective cells. Once the cells reach the box boundaries (, the convection is sustained by the constant-temperature, reflecting boundary conditions and produces sustained turbulence that isotropizes the magnetic field. Again, MFM can capture all the important behaviors even at very low resolution, especially for the larger box studied here which produces larger Mach numbers (here ) compared to Fig. 11. The bouyancy time for this setup is .

4.7 Anisotropic Diffusion-Driven Instabilities: MTI & HBI

Figure 13: Time evolution of our low-resolution MTI (left) and HBI (right) simulations from Fig. 12 (§ 4.7.2). Top: rms Mach number (volume-averaged) in the center of the computational domain, as a function of time. Bottom: Mean squared -component of the magnetic field orientation, . In MFM, the instabilities develops and velocities increase from their seed values rapidly around . The magnetic fields are rapidly re-oriented at the same time. In the MTI, fields go from horizontal to vertical, then the sustained convection with steady-state Mach numbers (given the large box-size ) isotropizes the fields (red dotted line represents isotropic fields). The fluctuations about isotropy owe to the low resolution. In the HBI, fields go from vertical to isotropic, but the large Mach numbers (), low resolution, and boundary conditions here cause an overshoot that sustains fluctuations around isotropy until the velocities decay to , at which point the instability rapidly completes the horizontal re-orientation of the field.

We now consider two instabilities specific to anisotropically conducting plasmas: the magneto-thermal instability (MTI) and heat-flux driven bouyancy instability (HBI) (Balbus, 2000; Quataert, 2008). These have been studied in various astrophysical contexts as drivers of turbulence, convection, and mechanisms to enhance or suppress conduction; our specific problem setup is motivated by the studies in Parrish & Stone (2005); Parrish et al. (2008); Parrish & Quataert (2008); McCourt et al. (2011); Kannan et al. (2015).

Here, we are not interested in the physics of the instabilities themselves, but they are useful numerical tests for several reasons. (1) They require accurate coupling of the anisotropic conduction to the magneto-hydrodynamics of the flow (not guaranteed in operator-split methods). (2) They are specific to anisotropic conduction and are suppressed with isotropic conduction, so directly test whether isotropic numerical diffusion can overwhelm physical diffusion. (3) They test the ability of the anisotropic conduction operator to recover small-amplitude seed perturbations. (4) They lead to non-linear, sub-sonic turbulence, which is particularly challenging for mesh-free methods (historically, SPH) to treat accurately (see e.g. Bauer & Springel, 2012; Price, 2012a), and require that our operator preserves anisotropy even in a turbulent flow.

4.7.1 High-Resolution Tests

Following Parrish & Stone (2005); Parrish et al. (2008); Parrish & Quataert (2008); McCourt et al. (2011); Kannan et al. (2015), we initialize a 2D box in coordinates with an analytic constant gravitational acceleration , size , resolution , polytropic , and conductivity ( for the isotropic case, for the anisotropic case). For the MTI, we initialize with , and , with a seed velocity perturbation333In our high-resolution tests, we use a seed velocity perturbation with magnitude instead of as in Parrish & Stone (2005); Parrish et al. (2008); Parrish & Quataert (2008). This was chosen for computational convenience. Because we use an explicit integration method, the timestep for high resolution in a small spatial domain becomes very short (). Using the smaller seed velocity requires approximately an order-of-magnitude longer runtimes to develop non-linear behavior, so we adopt the larger seed for high-resolution tests. However we explicitly show in § 4.7.2 below that our MFM method can accurately capture very small seed velocities even at much lower resolution . with . For the HBI, with , and and . The density is initialized so that the initial pressure gradient balances gravity with in the box. The domain is divided vertically into three equal sub-domains of length ; the top and bottom layers have isotropic conduction (i.e. are buoyantly neutral) while the central layer has anisotropic conduction (this follows the previous studies and reduces the sensitivity to boundary conditions). The boundaries are periodic in ; and constant-temperature reflecting boundaries in . Note that a sharp, reflecting and conducting “wall” is particularly challenging to implement in mesh-free methods. We treat the reflecting boundaries as follows: a layer (3 particles deep) of boundary particles with fixed positions and temperatures is placed at and . For every interaction between a boundary particles and normal particle (gradient calculation, the hydrodynamic operations, etc), the boundary particle is assumed to have all specific properties matched to the normal particle, except the temperature (fixed to the IC value), velocity and magnetic field (which follow the usual reflection rules) and density (adapted to give equal pressure, given the different temperatures).

Both of these ICs represent a stably stratified atmosphere. With no conduction, or isotropic conduction, the system should remain in equilibrium and the seed perturbations should damp. With anisotropic conduction, provided a large enough diffusivity (as chosen here) such that the system is approximately isothermal along field lines, the instabilities should grow and eventually re-orient the field lines. In the MTI (vertical temperature profile ), the small vertical velocity grows into large, non-linear convective cells which carry the field lines and re-orient the field in the vertical direction, until the cells cross the domain and the fixed-temperature boundary conditions produce sustained turbulence. In the HBI (), the initial mixed perturbation generates growing separation/compression of field lines, which leads to horizontal motions that try to stretch the field lines horizontally, until the instability saturates when the field lines are horizontal and suppress further convection. In both the HBI and MTI, the characteristic timescale is the bouyancy time ( in code units for our MTI and HBI setups, respectively).

Fig. 11 shows the results of these tests at a few bouyancy times, using our MFM method. Recall the relatively large initial perturbation here means non-linear behavior can develop in just a couple bouyancy times, and we see that occur in the anisotropic case. The qualitative behavior of the non-linear HBI and MTI compares well to that seen in Parrish & Stone (2005); Parrish & Quataert (2008); McCourt et al. (2011). In both MTI and HBI, the cases with isotropic diffusion show (correctly) no evidence of instability, and the initial perturbations are eventually fully damped.

4.7.2 Low-Resolution Tests

To explore late-time evolution and resolution effects, we now consider an initial condition which is identical to our high-resolution tests, but with lower resolution and larger box size (these both allow larger timesteps), smaller seed velocity , and larger conductivity (to preserve the desired limit where the conduction is sufficiently fast).

Fig. 12 shows the time evolution of the MTI in one of these runs; we obtain similar (good) behavior in our low-resolution MTI runs. In both cases, even at resolution, the instabilities are (correctly) completely suppressed with isotropic diffusion. Quantitatively, the results from these runs are shown in Fig. 13.

In MFM, we see the MTI grow (despite the very low resolution), going non-linear after bouyancy times (similar to high-resolution cases with the same seed perturbation amplitude; Parrish & Stone 2005; McCourt et al. 2011; Kannan et al. 2015), at which point the convective plumes re-order the magnetic field from horizontal to vertical. At late times, the plumes break up into sustained, non-linear convection, which isotropizes the field (although this is known to be sensitive to the boundary conditions). The fluctuations around isotropy owe to the ongoing turbulence (and are smaller at higher resolution). The saturated Mach numbers are larger than those in the smaller box, as expected based on the scaling seen in McCourt et al. 2011. Both qualitative and quantitative evolution compare favorably to higher-resolution studies in fixed-grid and moving-mesh schemes (McCourt et al., 2011; Kannan et al., 2015).

For the HBI, the instability goes non-linear on a similar timescale to the MTI, but the Mach numbers slowly decay as the field re-aligns. Interestingly, for the box, there is an intermediate period of isotropic fields, before the instability completes their horizontal re-alignment; this does not appear in the box, so we suspect it owes to the larger turbulent “overshoot” through the unstable zone induces by the first phase of the instability, which must be damped before the re-alignment can complete.

Figure 14: Hall MRI tests (§ 4.8). We compare the field lines (as Fig. 11) well into non-linear evolution. The initial condition is a 2D shearing box of length (the pressure scale length) with resolution, a constant field , trace pressure/velocity perturbations, and explicit Ohmic resistivity (with magnetic Reynolds number ) and Hall effects (with Hall parameter ). Opposite signs of correspond to the identical setups with opposite signs of ; this has no effects unless the Hall term is present. Results are shown at (). Left: For , modes should grow quickly. In the non-linear state, the fastest-growing modes become bigger than the box, leading to the horizontal channel modes seen; these increase without limit. Right: For , growth is slower and smaller modes grow fastest, leading to near-isotropic MRI turbulence in the saturated state as opposed to channel modes.
Figure 15: Growth of the volume-averaged magnetic pressure (relative to the initial thermal pressure ) in the mean-field Hall MRI tests. We show results for different initial Hall parameters (as labeled), at low resolution (; thin lines) and intermediate resolution (; thick lines). We also show the expected analytic maximum linear growth rate (; dotted black line), which should be close to the simulated growth rate for . The growth for is similar to the ideal MHD case (as expected), with growth rates in good agreement with the analytic expectation for resolution, and late-time formation of channel modes (Fig. 14). Growth rates are suppressed at very low resolution, as expected from the MRI studies in Paper II. For , linear growth is suppressed and the magnetic field strength saturates when the system saturates in MHD turbulence (also as expected). For , the system should be stable against MRI growth; we confirm this and that the initial seed noise damps at the expected rate (decay rate ). All the results here are in good agreement with well-tested Eulerian codes (see e.g. Sano & Stone 2002, Fig. 5).

4.8 Hall MRI

As discussed in § 4.1.4, the Hall effect presents unique numerical challenges; therefore, we now consider the magneto-rotational instability (MRI) with a Hall term and Ohmic resistivity (“Hall MRI”). This again involves the effects of anisotropic diffusion on plasma instabilities and turbulence. The effects of the Hall term are especially important in the context of proto-stellar and proto-planetary disks (for reviews, see Balbus, 2003; Wardle, 2007), and have been studied in detail using Eulerian methods (Flock et al., 2011; Simon et al., 2011; Bai, 2011; Simon et al., 2015).

The MRI itself is astrophysically interesting in a wide range of contexts involving magnetized disks, and is numerically particularly interesting because it has historically proven challenging for SPH methods to correctly capture its growth (see e.g. Rosswog & Price, 2007; Price & Bate, 2008; Dolag & Stasyszyn, 2009). In Paper II and Hopkins (2016), we consider extensive studies of the MRI in ideal MHD (, ) in GIZMO. We showed that our MFM and MFV methods recover the correct linear growth rates and non-linear properties in good agreement with well-tested higher-order Eulerian codes such as ATHENA. We also showed that at least some SPH schemes were capable of doing the same (albeit with greater noise), if a larger neighbor number is used.

Here, we follow Sano & Stone (2002) and consider a simplified test problem. We adopt a 2D axisymmetric shearing box from Guan & Gammie (2008); this is locally Cartesian with one coordinate () representing the radial direction and the other coordinate () the vertical direction, representing a small, azimuthally symmetric “ring” about the midplane of a disk in a Keplerian potential. Details of the boundary conditions and gravitational terms, as implemented in GIZMO, are in Paper II (this is the same box used for the MRI simulations therein). The box has initially constant density , side-length (where is the scale-height), orbital frequency , mean pressure (), spatially uncorrelated random pressure and velocity fluctuations with uniform distribution and , and uniform vertical field . The MRI is then characterized by three numbers, the plasma beta (); the magnetic Reynolds number , which determines the Ohmic resistivity ; and the Hall parameter , which determines (we assume the free electron fraction is constant).

Figs. 14-15 show the resulting evolution of the magnetic fields, for . For (no Hall term), we simply have the MRI with explicit resistivity. In 2D shearing boxes, the fastest-growing modes for should have growth rates only slightly smaller than the expected in the ideal limit; once non-linear, they should form an inverse cascade until horizontal channel modes appear which grow without limit. For , the behavior should be essentially identical, with slightly faster mode growth at higher . When , the system is fully stable against the MRI, and the initial perturbations should decay. At intermediate , the MRI should grow but with a suppressed maximum growth rate (for finite ). Within this range, as becomes more negative (larger ), smaller-scale modes grow faster, until for there is no critical scale at all; because of this, the system can saturate for many orbital periods with steady-state MRI turbulence (as the growing small-scale, high- modes prevent the formation of low- channel modes). We confirm all these behaviors here (compare Fig. 5 in Sano & Stone 2002). As expected and shown in detail in Paper II, the linear growth rates are suppressed at very low resolution (), but rapidly approach the analytic solution at higher resolution. Moreover, we note that we have re-run every problem in Sano & Stone (2002) with our MFM method (varying , , , and the box size; considering zero net-field cases; and their whistler wave test problem) and confirm identical qualitative behavior.

Note that, for , given our setup, the instantaneous Hall parameter and magnetic Reynolds number , so as modes grow non-linear ( increases), the system moves closer to the ideal MHD limit. Because our problem has finite , if we begin with , there is a fastest-growing mode, with relatively large wavelengths ( for ). During the turbulent phase, these modes can increase non-linearly, which in turn increase the fastest-growing mode growth rate and wavelength, and suppress the growth of the smaller-scale modes. This can eventually trigger a runaway inverse cascade that produces the channel modes seen for