Smoothed particle magnetohydrodynamics

Smoothed particle magnetohydrodynamics has reached a level of maturity that enables the study of a wide range of astrophysical problems. In this review, the numerical details of the modern SPMHD method are described. The three fundamental components of SPMHD are methods to evolve the magnetic field in time, calculate accelerations from the magnetic field, and maintain the divergence-free constraint on the magnetic field (no monopoles). The connection between these three requirements in SPMHD will be highlighted throughout. The focus of this review is on the methods that work well in practice, with discussion on why they work well and other approaches do not. Numerical instabilities will be discussed, as well as strategies to overcome them. The inclusion of non-ideal MHD effects will be presented. A prospective outlook on possible avenues for further improvements will be discussed.


INTRODUCTION
Smoothed particle magnetohydrodynamics (SPMHD) is a robust numerical method for solving the equations of magnetohydrodynamics (MHD).It is a Lagrangian, mesh-free method that builds upon the smoothed particle hydrodynamics (SPH) framework (Lucy, 1977;Gingold and Monaghan, 1977).The general picture of SPH is to solve the equations of hydrodynamics by discretising a fluid into a collection of particles that mimic fluid behaviour.Recent reviews of the fundamentals of SPH include Rosswog (2009), Springel (2010) and Price (2012).SPH, and by extension SPMHD, has many advantages for astrophysics.One, the resolution is tied to the mass.Regions of higher mass have more particles, thus more resolution, which is advantageous as the densest areas are typically the most interesting (e.g., stars forming in a molecular cloud).Two, it is trivial to incorporate gravitational N-body methods since SPH is a particle-based scheme.Three, advection is done perfectly, that is, without any dissipation, since it is a Lagrangian method.Four, it can easily handle complex geometries.Five, the Courant timestep does not depend upon the fluid velocity, thus allowing larger timesteps.And six, perhaps its strongest attribute, it has exact simultaneous conservation of mass, momentum, angular momentum, energy, and entropy to the precision of the time-stepping algorithm.This makes SPH significantly robust and stable since it reflects the conservation properties of nature.
The minimum requirements for incorporating MHD into SPH are three-fold: (1) an approach for evolving the magnetic field forward in time, (2) calculation of the accelerations deriving from the Lorentz force, and (3) an approach to uphold the divergence-free constraint of the magnetic field.These need not be independent.Adhering to the divergence-free constraint is strongly connected to the choice of how the magnetic field is numerically evolved, for instance.
The modern SPMHD method solves the set of continuum equations given by d dt where the stress S ij is defined as and ρ is the density, v is the velocity, B is the magnetic field, P is the pressure, and µ 0 is the permeability of free space.This system of equations is closed by a suitable equation of state.Note that the continuum equations are written with Lagrangian derivatives, d/dt ≡ ∂/∂t + v j ∂/∂x j , since SPMHD is a Lagrangian method.
The first attempts to include magnetic fields in SPH were performed by Gingold and Monaghan (1977) who considered magnetic polytropes, though in a form which did not conserve momentum or angular momentum.The modern SPMHD method has its roots in the work by Phillips and Monaghan (1985), who formulated equations of motion that conserve momentum by using the stress tensor, and applied the method to three-dimensional simulations of gravitationally collapsing gas clouds (Phillips, 1986a,b).
The challenge with the conservative approach is that it is numerically unstable when the magnetic pressure exceeds the thermodynamic pressure, that is, for plasma β < 1 where β ≡ P/(B 2 /2µ 0 ).Though the onset of instability is described by this criterion, the fundamental nature of the instability is linked to magnetic accelerations which are not perpendicular to the magnetic field, that is, that originate from magnetic monopoles.Solutions to this instability invariably incur some penalty to the conservation of momentum (Meglicki et al., 1995;Morris, 1996;Børve et al., 2001), though at a cost related to the magnitude of divergence errors in the magnetic field.Again, the linkage between calculating the accelerations deriving from the magnetic field and the approach used to uphold the divergence-free constraint of the magnetic field should be noted.
Carefully evolving the magnetic field in such a way as to avoid divergence errors altogether would be ideal.For grid codes, this can be accomplished with, for example, constrained transport (Evans and Hawley, 1988), but a significant impediment to do this in SPMHD is the lack of any welldefined surfaces.To date, no satisfactory approach to ensure a truly divergence-free magnetic field has yet to be found for SPMHD (Price, 2010;Brandenburg, 2010).Instead, modern SPMHD Tricco straightforwardly discretises the induction equation, with the recognition that divergence errors may be generated as there is no intrinsic divergence control.The magnetic field must be corrected through some other procedure.In practice, "divergence cleaning" (Tricco and Price, 2012;Tricco et al., 2016a) works well to remove divergence errors and is the method used in many modern SPMHD calculations.
Other pieces of the SPMHD puzzle that had to be solved include construction of the fully conservative equations incorporating varying resolution and magnetic discontinuity capturing terms (Price andMonaghan, 2004a,b, 2005).It worthy to note that modern SPMHD is equivalent in formulation to the eight-wave approach of Powell (1994); Powell et al. (1999), in that the SPMHD equations have source terms related to the divergence of the magnetic field, such that divergence errors are advected with the fluid flow.Additionally, non-ideal MHD formulations for Ohmic dissipation, ambipolar diffusion, and the Hall effect have been constructed (Wurster et al., 2014(Wurster et al., , 2016)).SPMHD has also been adapted for 2D axisymmetric geometry, promising computational efficiency for phenomena involving magnetic fields with axial geometries (García-Senz et al., 2023).
The outline of this review is as follows.A brief reminder of the fundamental SPH equations is provided first (section 2).The ingredients comprising modern SPMHD are subsequently discussed, focusing on the three main requirements of approaches to evolve the magnetic field (section 3), to calculate accelerations from the magnetic field (section 4), and to maintain the divergence-free constraint of the magnetic field (section 5).Non-ideal extensions are described in section 6, and a prospective outlook is given in section 7.

SMOOTHED PARTICLE HYDRODYNAMICS
In this review, a standard form of SPH is assumed whereby the density of an SPH particle is calculated by mass-weighted summation.A number of SPH variants and implementations have been proposed over time (Ritchie and Thomas, 2001;Hopkins, 2013;Saitoh and Makino, 2013;García-Senz et al., 2022), but here we focus on traditional mass-weighted SPH.The mass-weighted density summation is given by with h the smoothing length, m the mass, and summation is over neighbouring particles within the support radius of the smoothing kernel, , where r is the particle coordinates.
One could discretise Equation ( 1) directly, though a benefit of Equation ( 6) is that the density can then be computed at any location at any time.There is no need to time integrate the density.This also avoids assumptions about the differentiability of the density.See also Price (2012) for a nuanced presentation on the trickle down effects of using this SPH density estimate.
Smoothing lengths that are individual per particle can be obtained by simultaneously solving the density summation (Equation 6) with an expression for density given by where n dim is the number of dimensions and n h is a dimensionless quantity specifying the ratio of smoothing length to particle spacing.The smoothing length that yields agreement between Equations ( 6) and ( 7) can be obtained via a rootfinding procedure.
The conservative equations of motion derived from the Lagrangian for the discretised system need only specify the density summation in Equation ( 6) (Springel and Hernquist, 2002;Monaghan, 2002).

Doing so yields
where factors are present to account for spatially varying smoothing lengths.Obtaining ∂h a /∂ρ a may be done through Equation ( 7).See Monaghan (2005); Springel (2010); Price (2012) for detailed derivations of the SPH equations of motion from the Lagrangian.
If the pressure is a function of internal energy, u, then a suitable equation must be used to evolve the internal energy forward in time.This can be derived from the first law of thermodynamics, yielding The time derivative of ρ can be obtained by taking the time derivative of the density summation, Equation (6).Using this, the discretised internal energy equation is Note that the total specific energy or entropy of each particle could be evolved instead.The differences between these choices are minimal (see Price, 2012).
The choice of smoothing kernel has many considerations.In practice, the two families of kernels that are most widely used for astrophysical SPMHD calculations are the bell-shaped B-splines and the Wendland kernels.The kernel function can be written in functional form according to where σ is the normalisation and q = r ab /h a .The cubic spline is given by with σ = 1/π in 3D, the quintic spline by (14) with σ = 1/(120π) in 3D and Wendland C4 kernel, scaled to a radius of 2h, defined as with σ = 495/(256π).The cubic spline has a long history within SPH since it is the lowest order spline that has a continuous first derivative.The quintic spline (or other higher order splines) may be used when higher accuracy is required and the kernel bias needs to be reduced.The B-splines cannot be scaled to arbitrarily large radii because doing so will lead to the formation of close particle pairs.The Wendland kernels are attractive because they are stable against particle pairing for all neighbour numbers owing to their positive definite Fourier transform (Dehnen and Aly, 2012).The price they pay for this is that they yield larger density errors than the B-splines at lower neighbour numbers.The choice of kernel plays a role in numerical convergence and the reduction of 'E0' errors from the pressure gradient (Read et al., 2010;McNally et al., 2012;Dehnen and Aly, 2012;Hopkins, 2015), but this is dependent upon the particular details of a calculation.In the context of SPMHD, both the cubic spline and Wendland C4 kernel have been used successfully within a variety of dynamo and astrophysical applications, so the choice of kernel seems less important than other numerical details.
The approach that works best for evolving the magnetic field forward in time is to directly discretise the induction equation.This is also the most straightforward.Once the time derivative of the magnetic field has been calculated, it can be used within standard time integration schemes.A popular option is a kick-drift-kick leapfrog scheme owing to its efficiency, simplicity and conservation properties (Springel, 2005;Wadsley et al., 2017;Price et al., 2018).

Induction equation
Discretising the induction equation is as simple as applying standard differencing operators to Equation (3).Each particle evolves B/ρ according to where v ab = v a − v b .The magnetic field, when needed, can simply be reconstructed by multiplying the evolved quantity by the density.
One could evolve B directly instead of B/ρ.Expanding the left-hand side of Equation (3) yields and by making use of Equation (3) and the continuity equation (Equation 1), The corresponding discretisation is Numerically, any difference between Equations ( 16) and ( 19) should arise solely by the addition of the dρ/dt term when evolving B directly.The differential form of the continuity equation is baked into the evolution of B, whereas evolving B/ρ uses the integral formulation of the continuity equation (i.e., the density summation) to reconstruct B. In theory, evolving B/ρ should thus be preferable over B in situations where a discontinuity is present in the density as it avoids assumptions about the differentiability of the density (see also Price, 2008).In practice, however, no substantive evidence has been found of any meaningful difference between evolving B/ρ and B (e.g., Morris, 1996) An advantageous property of the above discretisations, aside from their simplicity, is that they are Galilean invariant.That is, the addition of any constant background velocity does not introduce numerical error.
The most significant disadvantage of these discretisations is that they do not place any guarantee on the divergence of the magnetic field.Even with an initially divergence-free magnetic field, numerical errors will lead to divergence errors that can grow over time.A secondary approach is required to treat these errors.The most effective approach to accomplish this at present seems to be mixed hyperbolic/parabolic divergence cleaning (Tricco and Price, 2012;Tricco et al., 2016a), as is discussed further in section 5.
Additionally, it can be seen that the induction equation solved in SPMHD mimics the source term approach of of Powell (1994); Powell et al. (1999).Consider the conservative formulation of the induction equation given by The vector calculus identity 20) in terms of a Lagrangian derivative, yielding This is identical to Equation (18) except for the term proportional to ∇ • B. In SPMHD, this term Tricco is ignored, in essence evolving the magnetic field under the assumption that the divergence of the magnetic field is zero.This treatment of monopoles is equivalent to how monopoles are treated by the eight-wave formulation.
The inclusion or absence of the ∇ • B term consequently means that divergence errors are either, respectively, advected or dispersed.Taking the divergence of Equation ( 21) yields, whereby the constraint for the divergence of the magnetic field only appears as an initial condition.This preserves the volume integral of B. By contrast, the divergence of Equation ( 18) yields In this case, the divergence of the magnetic field is treated in the same manner as density via the continuity equation, in that the volume integral of ∇ • B is conserved.The implication is that divergence errors are advected with the fluid flow.For the previous case, divergence errors are dispersed, and it has been shown that including the term proportional to ∇ • B leads to a poorer treatment of divergence errors in SPMHD (Price and Monaghan, 2005).

Example: Source term in the induction equation
The advection of a divergence 'blob' illustrates the difference between solving Equations ( 18) and (21), that is, whether the term proportional to the divergence of the magnetic field is included or not.These results mimic those presented by Price and Monaghan (2005), with the divergence advection test originating from Dedner et al. (2002).The test is set up using 50×50 particles on a square lattice set within x, y ∈ [−0.5, 1.5].The density is uniformly ρ = 1, with P = 6 and γ = 5/3.The initial velocity field is v = [1,1], such that the fluid flows from the bottom left towards the top right.The initial magnetic field is B x = 1/ √ 4π, with a perturbation B x = [(r/r 0 ) 8 − 2(r/r 0 ) 4 + 1]/ √ 4π for r < r 0 , with r 0 = 1/ √ 8.The perturbation to the magnetic field artificially introduces a non-zero divergence to the magnetic field.
Figure 1 shows the transport of the divergence error between the source term (Equation 18) and volume conservative (Equation 21) approaches.When the term related to the divergence of the magnetic field is excluded (equivalent to the source term approach of the eight-wave solver), the divergence error is passively advected with the fluid flow.The formulation remains consistent in the presence of divergence errors.By contrast, with this term included, the volume conservative form spreads the divergence error throughout the domain.

Artificial resistivity
The discretised induction equation assumes that the magnetic field is differentiable, which would not be true at discontinuities since the magnetic field would be multi-valued.An artificial resistivity is applied to the magnetic field to smooth discontinuities over the resolution scale so that the magnetic field remains single valued.Developed by Price andMonaghan (2004a, 2005), artificial resistivity adds dissipation to the magnetic field according to , represents the speed of information propagation between two particles.Energy dissipated from the magnetic field may be added to the internal energy, u, through (25) Importantly, the deposition of dissipated magnetic energy into the internal energy is guaranteed to be positive definite.Unlike artificial viscosity, artificial resistivity is applied to both approaching and receding particles, since discontinuities in the magnetic field can occur during both compression and rarefaction, and to all components of the magnetic field (rather than just along the line of sight like artificial viscosity), since magnetic discontinuities can occur oblique to the motion (Price andMonaghan, 2004a, 2005).
One option for the signal velocity is which uses the averaged fast magnetosonic wave speed, v mhd .A dimensionless parameter, α B , may be included per particle to switch dissipation off in regions that are not discontinuous, with such that dissipation increases in the presence of strong gradients of the magnetic field.Barnes et al. (2012) found that α B,0 = 0 leads to satisfactory results in their shock tube and other two-dimensional MHD tests.
However, Tricco and Price (2013) noted that the Price and Monaghan (2005) switch will not apply sufficient dissipation to capture discontinuities if strong shocks are present and the magnetic field is very weak (i.e., plasma β ∼ 10 6 -10 10 ), leading to numerical noise in the magnetic field and spurious growth of magnetic energy.In such a regime, using the Alfven speed instead of fast magnetosonic speed No tunable dissipation parameter in the form of α B is required.It is noteworthy to realize that this switch is not dependent upon the magnetic field at all.Yet, despite this, the switch in Equation ( 29) appears to provide less dissipation and work robustly in a variety of real applications, for example, within simulations of the MRI (Wissing et al., 2022) or the collapse of molecular cloud cores to stellar densities (Wurster et al., 2022).
Though numerical stability is the purpose of including artificial resistivity in simulations, its implementation is equivalent to a (resolutiondependent) physical dissipation term η∇ 2 B with magnetic diffusivity This is analogous to the translation of artificial viscosity to a physical viscosity (Artymowicz and Lubow, 1994;Lodato and Price, 2010).Keep in mind that advection in SPMHD is dissipationless since it is a Lagrangian scheme.This means that, for ideal MHD calculations with only numerical dissipation, it is possible to estimate the magnetic Reynolds number or magnetic Prandtl number directly from the artificial dissipation terms (e.g., Wissing et al., 2022).
Measurements of the magnetic Prandtl number, Pm = ν/η, where ν is the kinematic viscosity, typically range between Pm ∼ 0.1-2 when the only sources of dissipation are artificial viscosity and resistivity.There is slight dependence upon the numerical resolution as given by the particle smoothing lengths, but of more importance seems to be the specific problem under consideration.Overall, artificial resistivity works well to treat magnetic discontinuities.However, it is a crude tool for this purpose as it is difficult to reliably detect magnetic discontinuities, and it does not distinguish between Alfvén waves and compressive waves.Iwasaki and Inutsuka (2011) use a Godunov scheme instead of artificial resistivity to address these issues.This offers reduced numerical dissipation, but brings additional computational expense, whereas artificial resistivity is computationally inexpensive as it does not introduce additional time constraints or extra loops over the particles.

Example: Brio-Wu shock tube
The accuracy of the artificial resistivity switch given by Equation ( 29) can be demonstrated with the shocktube from (Brio and Wu, 1988), which also corresponds to shocktube 5A from Ryu and Jones (1995).This test is a magnetic extension of the Sod (1978) shocktube.The results of this particular shocktube configuration have been presented multiple times for SPMHD, for example, amongst others, Børve et al. (2001); Price and Monaghan (2004a,b); Barnes et al. (2012); Tricco and Price (2013); Hopkins and Raives (2016); Tricco et al. (2016a); Price et al. (2018); Wissing and Shen (2020).Other shocktube configurations have been tested with SPMHD, some of which are in the preceding references.The initial conditions are given by [ρ, P ] = [1,1] for the left state, and [ρ, P ] = [0.125,0.1] for the right state.The equation of state uses γ = 5/3, as in Ryu and Jones (1995), which is in contrast to γ = 2 in Brio and Wu (1988).The test is calculated in full 3D, with 800×20×20 particles for the left state and 400×10×10 particles for the right state, initially arranged on cubic lattices.The velocity is zero for both the left and right states.The magnetic field vector is There is no smoothing of the initial conditions.
Figure 2 shows the results at t = 0.1.The black circles represent the values on the particles, with the red line a reference solution given by the grid-based code ATHENA (Stone et al., 2008) using 10 4 grid cells.This shocktube contains a slow compound structure, which is notably absent from the Riemann solver solution given by Ryu and Jones (1995).SPMHD captures all of the shock structures.For the magnetic field components in particular, B y agrees with the reference solution, with the discontinuities at x ∼ −0.03 and x ∼ 0.13 represented over 4-8 particle spacings.In the analytic solution, the other magnetic field components remain constant over time, with B x = 0.75 and B z = 0.This is maintained in the SPMHD solution, except for small deviations due to the initial contact discontinuity near x ∼ 0.06 that introduce noise into the particle arrangement.The absolute max deviations in B x and B z are on the order of ∼ 10 −3 .

Euler potentials
Alternative formulations to evolve the magnetic field have been investigated, though with little success.One approach that appeared promising on the surface was defining the magnetic field in terms of the Euler potentials, B = ∇α E ×∇β E , where α E

Tricco
and β E are scalar potentials.This upholds ∇ • B = 0 by construction.The scalar potentials themselves are simply advected with the fluid flow.This is trivial to implement into a Lagrangian scheme, as the scalars become just constant, passive values attached to each particle.The re-construction of the magnetic field is thus dependent solely upon the configuration of the potentials.
Calculating the gradient of each potential may be done with standard SPH derivative operators, such as the differencing derivative operator, ) Plugging the re-constructed magnetic field into the conservative, Lagrangian equations of motion (see section 4) yields a numerically stable scheme, and numerical formulations of this type have been used in several astrophysical applications, such as star formation (Price and Bate, 2007, 2008, 2009), neutron star mergers (Price and Rosswog, 2006) and the magnetic fields of galaxies (Dobbs and Price, 2008;Kotarba et al., 2009).
However, the Euler potentials are incapable of correctly capturing the full range of magnetic field topologies (Brandenburg, 2010).Since the potentials are passively advected, there is a oneto-one mapping from the initial magnetic field topology to any evolved field.This is problematic for rotating discs, for example, as the magnetic field would be essentially reset with each turn once the particles return to their initial positions.Even worse is that the Euler potentials are incapable of modelling dynamos, and are incompatible with linear diffusion operators (Brandenburg, 2010).

Vector potential
Defining the magnetic field in terms of the vector potential, B = ∇ × A, would seem an ideal choice.Doing so would guarantee a divergencefree magnetic field, without any of the topological restrictions associated with the Euler potentials.However, SPMHD implementations of the vector potential appear to be plagued with numerical instabilities (Price, 2010).
An evolution equation for the vector potential can be obtained by combining B = ∇ × A with the induction equation (i.e., Equation 20), which yields where ∇ϕ represents the gauge choice, arising as a constant of integration.Price (2010) used the gauge This has the effect of moving the derivative onto v instead of A, which makes the discretisation Galilean invariant.In spite of this, using the same approach to calculate the Lorentz force as employed with the Euler potentials, whereby the re-constructed magnetic field is used within the equations of motion (see section 4), will lead to numerical instability and runaway energy growth (Price, 2010).That this occurs for the vector potential and not the Euler potentials is understandable because the Euler potentials are held constant.

Tricco
Numerical instability and runaway energy growth occur due to the mismatch of variables between the evolution of the magnetic energy (evolving A) and kinetic energy (equations of motion using B), such that the total energy in this hybrid formulism is not conserved.In theory, this is solved by deriving the consistent, conservative equations of motion from the Lagrangian that directly use the vector potential.Unfortunately, this only seems to lead to a different set of numerical issues.The resultant equations of motion contain three terms, with the first term acting like a negative magnetic pressure, the second term involving second derivatives of the kernel (which are known to be highly sensitive to particle disorder), and the third term nested low-order derivative estimates (Price, 2010).These equations do fundamentally represent the MHD equations of motion in the continuum limit, such that, for example, the second and third terms partially encode the magnetic pressure to offset the negative pressure of the first term, but the inherent numerical issues mean that this is not realized in practice.Stasyszyn and Elstner (2015) explored options to revive the hybrid formulism, such as smoothing the re-constructed magnetic field, adding numerical dissipation of the vector potential, and choosing the Coulomb gauge, ϕ = ∇ • A, which is upheld through mixed hyperbolic / parabolic divergence cleaning (more on this in section 5.1).Tu et al. (2022) also smoothed the re-constructed magnetic field, with the Weyl gauge, ϕ = 0, albeit within a meshless finite mass (MFM) scheme.Further testing is required on the robustness of these approaches.Tricco and Price (2023) derived a novel discretisation for the evolution equation of A in terms of a volume integral, but which also failed to produce a valid numerical solution.At present, it remains to be demonstrated that there is a formulation of the vector potential in SPMHD that is viable for general astrophysical application.

ACCELERATIONS FROM THE MAGNETIC FIELD
The equations of motion for SPMHD are derived from the Lagrangian.The implication of doing this is that the discretised equations of motion are the physical equations governing the system of discrete particles.This imparts a number of desirable properties, such as exact conservation of momentum, energy and entropy.Some early work directly discretised the Lorentz force, J × B, where J = ∇ × B/µ 0 is the current density, but these formulations do not provide the conservation properties of the Lagrangian approach (Gingold and Monaghan, 1977;Meglicki et al., 1995) and are poor at capturing shocks (Morris, 1996).

Conservative equations of motion
The SPMHD Lagrangian is The SPMHD equations of motion can be derived by using a variational approach (Price and Monaghan, 2004b;Price, 2012).The only ingredients necessary are the prescription for density (via the density summation in Equation 6) and the specification on how the magnetic field is evolved.The induction equation given by Equation ( 16) is assumed (or equivalently Equation 19).

Tricco
The resultant equations of motion from the Lagrangian are The second term in Equation ( 36) represents an isotropic magnetic pressure and the third term an anistropic magnetic tension.The correspondence between the SPMHD momentum equation and the stress tensor is clear upon examination of Equation (36) with Equations ( 2) and ( 5), and, indeed, Equation (36) could be written in terms of the stress tensor as Note, however, that while these equations of motion exactly conserve momentum, energy and entropy (provided they are used in conjunction with the density summation and the induction equation prescribed in their derivation), they do not conserve angular momentum.The anisotropic magnetic tension is not invariant to rotation.Notably, this term is derived solely from the numerical choice for the induction equation.
Wissing and Shen (2020) investigated switching the gradient operator in the momentum equation to one based on geometric density averaging, given by This type of gradient operator is used for the (thermodynamic) pressure gradient in the GASOLINE2 code (Wadsley et al., 2017), and promises to reduce artificial surface tension type effects (Agertz et al., 2007).For SPMHD, Wissing and Shen (2020) find that the geometric density average formulism has a lower numerical resolution requirement to produce magneticallydriven outflows from the gravitational collapse of magnetised molecular cloud cores.This is promising and warrants further examination.
Writing the momentum equation in terms of the stress tensor creates a complication in that the anisotropic tension contains a component due to monopole moments.The third term in Equation ( 36) is equivalent to (B • ∇)B/µ 0 ρ + (∇ • B)B/µ 0 ρ in the continuum limit.The force contributions proportional to ∇ • B are present in order to be momentum conserving in the presence of monopoles.The issue of monopole forces in SPMHD is complicated, in that even for a magnetic field that is constant and uniform (i.e., ∇ • B = 0), the discretisation used in the momentum equation may produce monopole forces.This is related to the 'E0' errors often discussed in relation to the pressure gradient (Read et al., 2010;McNally et al., 2012;Dehnen and Aly, 2012;Hopkins, 2015).
The real challenge with monopole accelerations is not related to the discretisation, however, but arise from any non-zero divergence errors that may be present in the magnetic field.The induction equation in SPMHD makes no guarantee of a divergence-free magnetic field.The detrimental effect of (∇ • B)B/µ 0 ρ is most strongly connected to the severity of non-zero divergence errors.
Compounding the considerations of monopole accelerations is that the conservative equations of motion themselves are tensile unstable.The anisotropic magnetic tension is an attractive force, and if it exceeds the isotropic pressure, then the particles will unphysically clump.For Equation (36), this occurs for plasma β < 1 (Phillips and Monaghan, 1985;Børve et al., 2001Børve et al., , 2004)).Importantly, this criterion is only based upon the magnitude of the magnetic field.Fortunately, removing the tensile instability is not difficult.

Removing the tensile instability
The approach that seems to work best to counteract the tensile instability is to explicitly subtract the non-physical force arising from the monopole contribution (Børve et al., 2001).This aligns with the source-term approach of the eightwave solver (Powell, 1994;Powell et al., 1999).Importantly, the removal of monopole accelerations must use the same discretisation for ∇ • B as in the momentum equation, that is, where βa used to regulate the strength of the applied correction, with βa = B a corresponding to a full correction of the force term.Adding Equation (39) to the conservative momentum equation (Equation 36) yields a numerically stable solution.
The consequence of removing monopole accelerations is that the conservative equations of momentum no longer conserve momentum.The severity of this non-conservation of momentum is dependent upon the magnitude of divergence errors.In the worst case, this can corrupt the obtained solution (for a dramatic example of this, see Tricco and Price, 2012).Thus, controlling the magnitude and growth of divergence errors will improve momentum conservation.Børve et al. (2004) showed that the tensile instability can be corrected with β = 1 2 B, thereby subtracting only 1 2 (∇ • B)B/µ 0 ρ.In theory, this provides a factor of two improvement on the non-conservation of momentum.Barnes et al. (2012) recommended this for general SPMHD calculations, though Tricco and Price (2012) found that this might cause numerical artefacts as, while technically sufficient to prevent instability, it may leave particles in a near-pressureless state.
Compounding this is that subtracting the full strength term ( βa = B a ) causes dispersive errors in slow magnetosonic waves.The half strength term ( βa = 1 2 B a ) does not introduce these type of errors (see section 4.2.1, c.f. Iwasaki, 2015).
Since the tensile instability only manifests for plasma β < 1, Børve et al. (2006), Price et al. (2018) and Wissing and Shen (2020) have proposed schemes where the tensile instability correction is switched off in pressure-dominated regimes, suggesting, respectively, and Both provide a linear transition from full correction to no correction once plasma β is above a critical threshold.This provides exact momentum conservation when possible.
An alternative approach to deal with the tensile instability is to use a more accurate derivative estimate for the anisotropic term Morris (1996).This uses the conservative form for the isotropic hydrodynamic and magnetic pressure (first and second terms in Equation 36), but replaces the anisotropic force (third term in Equation 36) with

Tricco
This is not momentum conserving, but yields a numerically stable solution.A disadvantage is that the Morris (1996) approach cannot be switched off, and also has dispersive errors in slow magnetosonic waves (Iwasaki, 2015).

Example: Propagation of an isolated wave
The tensile instability correction term, depending upon its implementation, can introduce dispersive errors to slow MHD waves.Here, this is demonstrated using the test from Iwasaki (2015) of the propagation of an isolated wave.The initial density is ρ = 1, with 320×160 particles arranged on a square lattice in the domain x ∈ [−2, 2] and y ∈ [−1, 1].An isothermal equation of state is used, with P = 1 and sound speed, c s = 1.The magnetic field is uniform along the x-direction with plasma β = 0.1.The test is performed using the cubic and quintic spline kernels.
Figure 3 shows the rightward propagating wave at t = 0.5 for β = B (top panel) and β = 1 2 B (bottom panel), where β is the leading term in the tensile instability correction term given by Equation ( 39).The location of the expected peak of the wave at t = 0.5 is given by the vertical dashed line (x = 0.5).For β = B, dispersive errors cause supersonic velocities preceding the wave.In this case, the quality of the kernel affects the result, with the cubic spline showing a larger phase velocity than the quintic spline.For β = 1 2 B, there are no dispersive errors and the wave propagates (mostly) as expected.

DIVERGENCE-FREE CONSTRAINT OF THE MAGNETIC FIELD
Upholding the divergence-free constraint on the magnetic field is critical.In one respect, real magnetic fields are purely solenoidal, so modelling magnetic fields that preserve this topological constraint is only logical.From a theoretical perspective, magnetic fields within simulation should avoid unphysical configurations.In another respect, there are numerical justifications for  39).Applying the full correction term leads to supersonic dispersive errors, whereas the half correction term yields results consistent with the expected solution.minimizing divergence errors.Doing so avoids numerical artefacts in the momentum equation and the overall conservation of momentum and energy.
An important question is: how does one define a "divergence-free magnetic field" numerically?One such definition has already been used in Tricco Equation ( 39), whereby Instead of this symmetric gradient estimate, one could instead measure the divergence of the magnetic field with a difference gradient estimate according to (44) Measuring zero divergence of the magnetic field in one metric does not guarantee zero in another.As a case in point, neither of the above discretisations are zero when the magnetic field is specified in terms of the Euler potentials, even though it intrinsically guarantees a divergence-free magnetic field (Tricco and Price, 2012).Nonetheless, minimizing divergence error in one discretisation will often reduce divergence error in another discretisation.This is important for SPMHD, as it means the detrimental side effects of the tensile instability correction will benefit from divergence control in a different discretisation.
At present, there are no robust methods in SPMHD for exactly preserving the divergence-free constraint on the magnetic field.The challenges faced with the Euler potentials and vector potential implementations are discussed in sections 3.3 and 3.4.Instead, it would appear that the best option currently is to only approximately uphold this constraint.

Hyperbolic / parabolic divergence cleaning
Dedner et al. (2002) introduced a divergence cleaning approach that transports divergence errors away from their source and damps them through a set of coupled hyperbolic and parabolic equations.In essence, divergence errors are propagated through a damped wave equation.A new, scalar field, ψ, is coupled to the magnetic field to facilitate this.Tricco and Price (2012) adapted the Dedner et al. (2002) divergence cleaning scheme for SPMHD, with later improvements by Tricco et al. (2016a).The cleaning equations, in continuum form, are given by dB dt = −∇ψ, (45) (46) This formulation differs from Dedner et al. (2002), in that Equations ( 45) and ( 46) use Lagrangian derivatives, use ψ/c h for the evolved variable, and account for compression and rarefaction of the fluid.Also accounted for are wave cleaning speed, c h , and parabolic damping, τ , that vary in time.Tricco et al. (2016a) showed that the above system of equations can be combined to create a generalised wave equation of the form If c h , τ , ρ and the fluid velocity are held constant, this reduces to the usual damped wave equation, specified in the original Dedner et al. (2002) formulation, The SPMHD discretised cleaning equations are Note that the difference estimate is used for ∇ • B (Equation 44).These 'constrained' hyperbolic / parabolic divergence cleaning equations have formal guarantees about numerical stability.The discretised cleaning equations were derived by considering the energy content of the ψ field, e ψ = ψ 2 /2µ 0 ρc 2 h , ensuring that any change in the magnetic energy is conserved by an equal but opposite change in ψ energy so that total energy is conserved.Imposing this constraint further ensures that the parabolic damping will only ever siphon energy from the magnetic field.
The hyperbolic wave speed, c h , is typically taken to be the fast magnetosonic wave speed, so that propagation is at the fastest rate permissible within the Courant timestep condition.This is given, per particle, by ∆t a = Ch a /c h,a where C ∼ 0.3 is the Courant factor.Since the hydrodynamical timestep is already constrained by the fast magnetosonic speed, this does not impose any additional timestep constraint.It is possible to use a faster speed, so long as the timestep is commensurately decreased (Dobbs and Wurster, 2021).The rate of damping is set by with σ = 1 an optimal choice (Tricco and Price, 2012;Barnes et al., 2018).
The divergence error within a calculation can be measured by the dimensionless quantity Typically, 'constrained' divergence cleaning keeps the mean ϵ divB in a simulation around ∼1%.This is sufficient for practical applications, and, generally, reduces average divergence error by about an order of magnitude for most test problems and astrophysical applications.Divergence cleaning can subsequently yield multiple orders of magnitude improvement in momentum conservation through the connection of the tensile instability correction with the divergence of the magnetic field.An example of this is shown in Figure 4, which showcases the gravitational collapse of a molecular cloud core from Tricco and Price (2012) (their Figure 21).Without divergence cleaning, material is ejected out of the plane of the disc due to spurious momentum gain from the tensile instability correction arising from divergence errors.With divergence cleaning, the disc remains stable over long-term evolution.Figure 5 shows the corresponding total momentum change.A sink particle is introduced around one free-fall time.Divergence cleaning improves momentum conservation by approximately two orders of magnitude, which is sufficient to keep the sink particle within the disc.
Importantly, the 'constrained' divergence cleaning scheme is stable in the presence of free boundaries or sharp density contrasts (Tricco and Price, 2012).Price and Monaghan (2005) tested a nonconservative scheme on the Orszag-Tang vortex (Orszag and Tang, 1979), finding that it gives, at best, a factor two reduction in average divergence error, but could potentially increase overall divergence error through spurious energy creation (a 'divergence creating' scheme).Stasyszyn et al. (2013) demonstrate that a non-conservative scheme can corrupt shocktube solutions with sharp density jumps.
Although divergence cleaning uses the difference derivative estimate for ∇ • B, as given by Equation (44), it is possible in principle to use other discretisations.As an example, Tricco and Price (2012) constructed 'constrained' divergence cleaning equations that used the symmetric derivative estimate (Equation 43).Though numerically stable, this low-order estimate was found to be overly dissipative and could .Change in total linear momentum over 1.4 free-fall times for the gravitational collapse of a molecular cloud core (see Figure 4).A sharp injection of momentum occurs when a sink particle is introduced around one free-fall time.Divergence cleaning improves the conservation of momentum by two orders of magnitude, which is sufficient to keep the sink particle within the disc.Reproduced from Figure 22 in Tricco and Price (2012).
introduce artefacts into the physical portions of the magnetic field since this estimate is sensitive to particle disorder.
Figure 6 shows the density cross-sections at t = 0.5 and t = 1 for the calculation with divergence cleaning.The divergence error can be measured by the dimensionless quantity ϵ divB = h|∇ • B|/|B|, which provides a standard metric for comparison.Figure 7 shows the median ϵ divB over time for calculations of the 3D Orszag-Tang vortex with and without divergence cleaning.The shaded regions represent the inter-quartile range, that is, the 25 th to 75 th percentiles.There is a significant reduction in divergence error when divergence cleaning is applied.The median ϵ divB with divergence cleaning is ∼ 0.1% (the average ϵ divB is also ∼ 0.1%).This is below the 25 th percentile without divergence cleaning.

Alternative divergence control approaches
SPMHD calculations enjoy a minor level of divergence control through artificial resistivity.As discussed in section 3.2, the explicitly added dissipation from artificial resistivity is equivalent to a (resolution-dependent) physical dissipation, Though artificial resistivity is intended to treat discontinuities in the magnetic field, it may help dissipate magnetic monopoles.The caveat is that dissipation is also applied to physical portions of the field, so using artificial resistivity primarily as a divergence control measure has its drawback.
Periodically smoothing the magnetic field to remove small-scale noise has occasionally been suggested (Børve et al., 2001;Dolag and Stasyszyn, 2009;Stasyszyn and Elstner, 2015;Tu et al., 2022).However, the frequency of the smoothing procedure needs to remove small-scale noise without overly smoothing the field.Smoothing procedures are adhoc without any formal guarantees.
Formulations that yield a truly divergence-free magnetic field would be ideal.SPMHD appears to be far from achieving this currently (see discussion on the Euler potentials and vector potential in sections 3.3 and 3.4).One avenue that is under-explored are projection methods (Brackbill and Barnes, 1980;Tóth, 2000).Considering an "unclean" magnetic field, B * , it can be written in terms of its physical and unphysical components according to B * = ∇ × A + ∇ϕ, where A is the vector potential and is the physical portion of Tricco the field (since the divergence of the curl is zero).From this, it can be stated that ∇ • B * = ∇ 2 ϕ, and then by solving for ϕ, the divergence-free magnetic field can be obtained from B = B * − ∇ϕ.It is noteworthy to consider that projection methods are commonly used to obtain divergence-free velocity fields for the simulation of incompressible fluids in SPH (Cummins and Rudman, 1999;Hu and Adams, 2007;Lind et al., 2020).Price and Monaghan (2005) tested projection methods in SPMHD, finding that they worked well for simple test problems.The most significant disadvantage is the computational cost in solving the elliptic set of equations.Using a tree may help improve efficiency, but would still be similar in cost to tree-based gravity solvers.Individual particle timesteps add further complication.It may be worth revisiting projection methods and testing further.

NON-IDEAL MHD
Non-ideal MHD concerns plasmas that are only partially ionized.A population of neutral particles can introduce important deviations from the flux freezing condition of ideal MHD.In the context of star formation, these non-ideal effects may play a role in molecular cloud formation, disc fragmentation and braking, and the seeded magnetic field strength inside protostars (see reviews by Pudritz and Ray, 2019;Hennebelle and Inutsuka, 2019;Maury et al., 2022).
Consider a partially ionized plasma consisting of electrons, ions, and neutrally charged particles.Ohmic resistivity refers to the impediment of free electron flow in conducting plasmas.Ambipolar diffusion is the process of neutral particles drifting through ions.The magnetic field is tied to the electrons and ions, thus there is a movement of mass that does not result in transport of the magnetic field.Both Ohmic resistivity and ambipolar diffusion are dissipative processes.The Hall effect, on the other hand, is not dissipative, but rather dispersive.For the Hall effect, ions or charged dust grains are collisionally coupled to the neutral species, but the electrons are not.The magnetic field is only transported by the electrons.See, for example, Wardle and Ng (1999); Braiding and Wardle (2012a,b); Pudritz and Ray (2019);Tsukamoto et al. (2022) for discussion on non-ideal effects in the context of star formation.
The three non-ideal effects of Ohmic dissipation, ambipolar diffusion, and the Hall effect have all been implemented into SPMHD (Tsukamoto et al., 2013;Wurster et al., 2014) and have been used for a variety of simulations related to star formation (e.g., Tsukamoto et al., 2015a;Wurster et al., 2016Wurster et al., , 2019)).Generally speaking, the context in which non-ideal effects would be expected to play a substantial role is when the ionization rate is sufficiently low.For this reason, Wurster et al. (2014) implement non-ideal effects using a single species of SPMHD particles using the strong coupling approximation, though it is possible to model using two species of particles, e.g., ambipolar diffusion (Hosking and Whitworth, 2004).
The continuum equations for the non-ideal MHD effects are given by where B = B/|B|, that is, the unit vector in the direction of the magnetic field, and η O , η A , and η H are the coefficients for Ohmic resistivity, ambipolar diffusion, and the Hall effect, respectively.
Non-ideal MHD effects introduce additional timestep constraints for numerical stability.This is given by ∆t which is calculated for each non-ideal effect.et al., 1996;Tsukamoto et al., 2013;Wurster et al., 2016), though this is only applicable for Ohmic dissipation and ambipolar diffusion.

Ohmic dissipation
Tsukamoto et al. ( 2013) and Wurster et al. (2014) implemented Ohmic resistivity in SPMHD using a 'two-first derivatives' approach (also for ambipolar diffusion and the Hall effect).This approach makes including non-constant resistivity coefficients straightforward.The discretised equations are given by The curl of B is calculated first using a difference derivative estimate, then the curl of the result with a symmetric derivative estimate.This conjugate pair of derivative operators appear in many contexts (e.g., Cummins and Rudman, 1999;Price, 2010;Tricco and Price, 2012).The Ohmic dissipation yields a positive-definite increase in thermal energy according to Note that artificial resistivity is equivalent to an Ohmic dissipation (see section 3.2).In this case, it is calculated using a second derivative directly in a manner equivalent to Brookshaw (1985) and (Cleary and Monaghan, 1999), rather than with two first derivatives as above.Ohmic resistivity implemented with a direct second derivative does work well, provided the resistivity coefficient is spatially constant (Bonafede et al., 2011;Biriukov and Price, 2019).

Ambipolar diffusion
Ambipolar diffusion can be implemented using a two-first derivatives approach, with the inner derivative, (∇ × B) a , calculated using a difference derivative estimate (Equation 55).The second derivative is estimated with a symmetric derivative, according to A positive-definite increase in thermal energy is ensured (Wurster et al., 2014), with (60)

Hall effect
The Hall effect calculates (∇ × B) a in the same manner as Ohmic resistivity and ambipolar diffusion, that is, according to the difference derivative estimate in Equation ( 55).The outer derivative is calculated using a symmetric derivative estimate according to Equation ( 58), but with As the Hall effect is dispersive, and not dissipative, it only leads to the re-distribution of magnetic energy.There is no deposition into thermal energy.
The Hall effect introduces an additional wave type known as whistler modes (Sano and Stone, 2002;Pandey and Wardle, 2008;Bai, 2014).The left and right polarizations of Alfvén waves become Tricco asymmetric in the presence of the Hall effect, that is, they have different phase velocities.The righthanded polarization is called the whistler wave.The introduction of an additional wave type can be generally understood in that the Hall effect is hyperbolic in nature rather than parabolic.

PROSPECTIVE OUTLOOK
If there is any indication on the future potential of SPMHD as a tool to study astrophysical problems, one only needs to look at the large bodies of work that have been conducted over the past decade in a variety of astrophysical contexts.It may taken decades for SPMHD to reach its current level of maturity, but all major roadblocks have been cleared.SPMHD is finally a method that is generally applicable for the study of astrophysics.This is not to say that there is no room to improve the method further.Far from it.
With respect to grid-based codes, SPMHD possesses a number of advantages, such as its conservation properties, adaptive resolution with respect to the density, and absence of dissipation due to advection.On the other hand, the numerical dissipation arising from artificial viscosity and resistivity is typically greater than that stemming from Reimann-based solvers in grid codes.SPMHD is also typically more computationally expensive than grid-based codes, owing to the number of neighbours under the kernel (∼ 60-200) and the need for a nearest neighbour search.The inherent particle 're-meshing' motions additionally incur a level of background noise in the velocity field, though this can be addressed by using more neighbours or better kernels, but at increased computational expense (e.g., Tricco, 2019).For magnetic fields, SPMHD can only approximately uphold the divergence-free constraint at present, whereas grid-based codes can utilize approaches such as constrained transport to exactly solve for a divergence-free magnetic field.Furthermore, the magnetic field can only be specified where there is mass present (i.e., particles).
One area that deserves more study is the MRI.Deng et al. (2019) found that SPMHD adequately simulates the MRI in unstratified shearing boxes with net flux, but have decaying solutions for zeronet flux.They also found that SPMHD produces unphysical behaviour for stratified shearing boxes, though the calculations of Wissing et al. (2022) were able to sustain MRI-induced turbulence in stratified shearing boxes for over 100 orbits, seemingly due to the geometric density average formulism (see Figure 8).Wissing et al. (2022) were also able to sustain turbulence in zero-nut flux unstratified shearing boxes, provided the magnetic Prandtl number was above a critical threshold.Part of the complication of studying the MRI is the breadth of initial configurations and parameters.One avenue worth exploring is the inclusion of physical dissipation to set the Reynolds and magnetic Reynolds numbers, along with the magnetic Prandtl number.Another avenue that seems ripe is global disc simulations.SPH has a rich history of simulating accretion discs owing to its conservation properties, and applying SPMHD in this direction would appear sensible.
There is also need to explore the properties of SPMHD on other fundamental MHD phenomena and to validate correctness of obtained solutions.Dynamos are one class of problem deserving more attention.To date, there has only been limited exploration of the small-scale dynamo (Tricco et al., 2016b) and galactic dynamos (Wissing and Shen, 2023).And only the surface has yet been scratched on the properties of magnetised turbulence in SPMHD.
There should be continued effort to improve the fundamental accuracy or convergence of SPMHD.There are discussions and advancements within hydrodynamical schemes for Lagrangian, particlebased methods that may be relevant for SPMHD.García-Senz et al. (2012) developed an SPH scheme that has higher-order gradient estimates through matrix inversion, and Iwasaki and Inutsuka (2011) use Riemann solvers inside their Godunov-SPH scheme.Meshless finite mass (MFM) (Hopkins, 2015;Hopkins and Raives, 2016) is another Godunovtype scheme that uses a least-squares matrix gradient operator and Riemann solvers with slope limiters.Like SPMHD, MFM uses mixed hyperbolic / parabolic divergence cleaning to uphold the divergence-free constraint on the magnetic field, but also employs a constrainedgradient method for further reduction of divergence errors in the magnetic field (Hopkins, 2016).At present, it is difficult the quantify how important is the choice of numerical scheme as there is a lack of overlap in astrophysical simulations between MFM and SPMHD, with use of MFM focused primarily on galaxy formation (e.g., the FIRE simulations Hopkins et al., 2014Hopkins et al., , 2018Hopkins et al., , 2023) ) and SPMHD on star formation calculations.Perhaps there may be lessons from these numerical methods, or even others, in regards to improvements for SPMHD.
Coupling magnetic fields with charged dust is a perfect opportunity to extend SPMHD in new directions.Non-ideal MHD effects and dust are tied together, as dust grains can adsorb electrons or ions.There are already robust solvers for gasdust mixtures in SPH (Laibe and Price, 2012a,b;Price and Laibe, 2015), and there is nothing in principle preventing creation of algorithms to couple magnetic fields and dust together in a unified solver.Tsukamoto et al. (2021Tsukamoto et al. ( , 2023) ) have made initial steps in this direction, expanding upon the one-fluid dust approach of Laibe and Price (2014) to include magnetic fields.They make a number of simplifications, and better treatment on the differences between neutral and charged dust are still required (depending upon the physical regime).
Finally, it is worth noting that adapting existing SPH codes to include magnetic fields is straightforward.At its core, this involves evolving another variable (the magnetic field), adding magnetic dissipation terms, calculating accelerations from the magnetic field, and solving the divergence cleaning equations alongside the evolution of the magnetic field.While there are improvements that can still be made to the method, recent advancements and achievements with SPMHD suggest that it is capable of general theoretical studies of astrophysics, providing a complementary approach to grid-based methods.

Tricco
figures.Figures and analysis in this review made use of SARRACEN, a Python package for the analysis and visualization of SPH data (Harris and Tricco, 2023).

Figure 1 .
Figure 1.The advection of ∇ • B with and without a source term related to the divergence of the magnetic field in the induction equation (top and bottom rows, respectively).With the source term (top row), the volume integral of ∇ • B is preserved and the divergence of the magnetic field is advected with the fluid flow.Without the source term (bottom row), the divergence of the magnetic field is dispersed throughout the domain.

Figure 2 .
Figure2.Brio-Wu shocktube at t = 0.1 with γ = 5/3.The initial configuration leads to a compound shock structure.The black circles represent the SPMHD solution, with the red line a reference solution from the ATHENA code.The magnetic field components of the SPMHD solution agree with the reference solution.Notably, B x remains constant even in the presence of density jumps and velocity discontinuities, with only small deviations in B x and B z located near the initial contact discontinuity.Results are presented using the artificial resistivity switch given by Equation (29).

Figure 3 .
Figure3.The propagation of an isolated wave in a strongly magnetized medium (plasma β = 0.1) with the full tensile instability correction term (top) and the correction term reduced by half (bottom), that is, β = B and β = 1 2 B as used in Equation (39).Applying the full correction term leads to supersonic dispersive errors, whereas the half correction term yields results consistent with the expected solution.

Figure 4 .
Figure 4.The gravitational collapse of a molecular cloud core.Without divergence control (left panel), divergence errors lead to material becoming kicked out of the plane of the disc due to momentum errors arising from the tensile instability correction.With divergence cleaning (right panel), divergence errors are kept sufficiently low to avoid catastrophic momentum injection, and the disc remains planar.Reproduced from Figure 21 in Tricco and Price (2012).

Figure 7 .
Figure 7.The median divergence error across all particles for the 3D Orszag-Tang vortex, as measured by h|∇ • B|/|B|, with and without divergence cleaning (orange dashed line and blue solid line, respectively).The shaded regions represent the inter-quartile range, that is, the 25 th and 75 th percentiles.

Figure 8 .
Figure 8. Spacetime diagrams of magnetorotational instability (MRI) calculations in a stratified net flux shearing box over 100 orbits.The top panel shows the horizontal averaged radial magnetic field, and the bottom panel the azimuthal magnetic field.The characteristic 'butterfly' diagram is produced, whereby azimuthal fields are transported out of the plane of the disc and the direction of the magnetic field periodically reverses.Reproduced with permission from Figure 15 in Wissing et al. (2022).
Barnes et al. (2018)2005)ity is similarly problematic if shocks are present.They recommend a new switch, settingα B,a = h a |∇B a | |B a | .(28)Notethatthisdoesnotrequire the time evolution of α B .On test problems, this switch led to less overall magnetic energy dissipation than thePrice and Monaghan (2005)switch.However,Barnes et al. (2018)found that the Tricco and Price (2013) switch was too dissipative in cosmological simulations of galaxy cluster formation.
Wissing and Shen (2023)und Pm ≈ 1 in simulations of supersonic, magnetized turbulence, andWissing et al. (2022)found Pm ≈ 1.5 in stratified, net-flux shearing box simulations.Wissing and Shen (2023)performed simulations of the galactic dynamo in spiral galaxies, measuring Pm ≈ 0.1 in the inner radius of the galaxy and Pm ≈ 0.7 in the outer regions.