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.

1 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 (Gingold and Monaghan, 1977; Lucy, 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.

Over the past decade, SPMHD has been used to simulate the evolution and impact of magnetic fields in a wide variety of astrophysical problems, such as the study of single and binary star formation (Bürzle et al., 2011a; Bürzle et al., 2011b; Price et al., 2012; Bate et al., 2014; Tsukamoto et al., 2015a; Tsukamoto et al., 2015b; Lewis et al., 2015; Wurster et al., 2016; Wurster et al., 2017; Lewis and Bate, 2017; Wurster et al., 2018a; Wurster et al., 2018b; Tsukamoto et al., 2018; Tsukamoto et al., 2020; Wurster et al., 2021; Wurster et al., 2022), star cluster formation (Wurster et al., 2019; Dobbs and Wurster, 2021), star formation rates in spiral galaxies (Herrington et al., 2023), accretion discs (Forgan et al., 2017), tidal disruption events (Bonnerot et al., 2017), the magnetic field structure of spiral galaxies (Dobbs et al., 2016; Wissing and Shen, 2023), and galaxy cluster formation (Barnes et al., 2012; Barnes et al., 2018). It has been shown to yield correct behaviour for the small-scale dynamo amplification of magnetic fields (Tricco et al., 2016b), and can sustain turbulence incited by the magnetorotational instability (MRI) (Deng et al., 2019; Wissing et al., 2022). SPMHD has also found use outside of astrophysics, where it has been employed to study pinch plasmas and fusion (Vela Vela et al., 2019; Thompson and Cassibry, 2020; Park et al., 2023).

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

$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \frac{\partial v^i}{\partial x^i},(1)$

$\frac{\mathrm{d}{v}^i}{\mathrm{d}t}=\frac{1}{\rho }\frac{\partial S^{ij}}{\partial x^j},(2)$

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{B^i}{\rho }\right)=\frac{B^j}{\rho }\frac{\partial v^i}{\partial x^j},(3)$

$\frac{\partial B^i}{\partial x^i}=0,(4)$

where the stress S^ij is defined as

$S^{ij}=-P{\delta }^{ij}+\frac{1}{{\mu }_0}\left(B^i{B}^j-\frac{B^2}{2}{\delta }^{ij}\right),(5)$

and ρ is the density, v is the velocity, B is the magnetic field, P is the pressure, and μ₀ 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; Phillips, 1986b).

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μ₀). 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 well-defined surfaces. To date, no satisfactory approach to ensure a truly divergence-free magnetic field has yet to be found for SPMHD (Brandenburg, 2010; Price, 2010). Instead, modern SPMHD 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 and Monaghan, 2004a; Price and Monaghan, 2004b; Price and Monaghan, 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.

2 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

${\rho }_a={\displaystyle \sum _b}{m}_b{W}_{ab}\left(h_a\right),(6)$

with h the smoothing length, m the mass, and summation is over neighbouring particles within the support radius of the smoothing kernel, W_ab(h_a) ≡ W(|r_a − r_b|, h_a), where r is the particle coordinates. One could discretise Eq. 1 directly, though a benefit of Eq. 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 (Eq. 6) with an expression for density given by

${\rho }_a=m_a{\left(\frac{h_a}{n_{\text{h}}}\right)}^{-{\mathrm{n}}_{\mathrm{d}\mathrm{i}\mathrm{m}}},(7)$

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 Eqs 6, 7 can be obtained via a root-finding procedure.

The conservative equations of motion derived from the Lagrangian for the discretised system need only specify the density summation in Eq. 6 (Monaghan, 2002; Springel and Hernquist, 2002). Doing so yields

$\frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}=-{\displaystyle \sum _b}{m}_b\left[\frac{P_a}{{\mathrm{\Omega }}_a{\rho }_a^2}{\nabla }_a{W}_{ab}\left(h_a\right)+\frac{P_b}{{\mathrm{\Omega }}_b{\rho }_b^2}{\nabla }_a{W}_{ab}\left(h_b\right)\right],(8)$

where factors

${\mathrm{\Omega }}_a=1-{\displaystyle \sum _b}{m}_b\frac{\partial W_{ab}\left(h_a\right)}{\partial h_a}\frac{\partial h_a}{\partial {\rho }_a}(9)$

are present to account for spatially varying smoothing lengths. Obtaining ∂h_a/∂ρ_a may be done through Eq. 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

$\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{P}{{\rho }^2}\frac{\mathrm{d}\rho }{\mathrm{d}t}.(10)$

The time derivative of ρ can be obtained by taking the time derivative of the density summation, Eq. 6. Using this, the discretised internal energy equation is

$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{P_a}{{\mathrm{\Omega }}_a{\rho }_a^2}{\displaystyle \sum _b}{m}_b{\mathbf{v}}_{ab}\cdot {\nabla }_a{W}_{ab}\left(h_a\right).(11)$

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

$W_{ab}\left(h_a\right)=\frac{\sigma }{h^{{\mathrm{n}}_{\mathrm{d}\mathrm{i}\mathrm{m}}}}w\left(q\right),(12)$

where σ is the normalisation and q = r_ab/h_a. The cubic spline is given by

$w\left(q\right)=\left\{\begin{array}{cc}{\left(2-q\right)}^3-4{\left(1-q\right)}^3\hfill & q<1,\hfill \\ {\left(2-q\right)}^3\hfill & 1\le q<2,\hfill \\ 0\hfill & q\ge 2,\hfill \end{array}\right.(13)$

with σ = 1/π in 3D, the quintic spline by

$w\left(q\right)=\left\{\begin{array}{cc}{\left(3-q\right)}^5-6{\left(2-q\right)}^5+15{\left(1-q\right)}^5\hfill & q<1,\hfill \\ {\left(3-q\right)}^5-6{\left(2-q\right)}^5\hfill & 1\le q<2,\hfill \\ {\left(3-q\right)}^5\hfill & 2\le q<3,\hfill \\ 0\hfill & q\ge 3,\hfill \end{array}\right.(14)$

with σ = 1/(120π) in 3D and Wendland C4 kernel, scaled to a radius of 2h, defined as

$w\left(q\right)=\left\{\begin{array}{cc}{\left(1-q/2\right)}^6\left(35q^2/12+3q+1\right)\hfill & q<2,\hfill \\ 0\hfill & q\ge 2,\hfill \end{array}\right.(15)$

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; Dehnen and Aly, 2012; McNally et al., 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.

3 Evolving the magnetic field

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).

3.1 Induction equation

Discretising the induction equation is as simple as applying standard differencing operators to Eq. 3. Each particle evolves B/ρ according to

$\frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathbf{B}}{\rho }\right)}_a=-\frac{1}{{\mathrm{\Omega }}_a{\rho }_a^2}{\displaystyle \sum _b}{m}_b{\mathbf{v}}_{ab}\left[{\mathbf{B}}_a\cdot {\nabla }_a{W}_{ab}\left(h_a\right)\right],(16)$

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 Eq. 3 yields

$\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}=\rho \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\mathbf{B}}{\rho }\right)+\frac{\mathbf{B}}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}t},(17)$

and by making use of Eq. 3 and the continuity equation (Eq. 1),

$\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}=\left(\mathbf{B}\cdot \nabla \right)\mathbf{v}-\left(\nabla \cdot \mathbf{v}\right)\mathbf{B}.(18)$

The corresponding discretisation is

$\begin{array}{ll}\hfill \frac{\mathrm{d}{\mathbf{B}}_a}{\mathrm{d}t}& =-\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b{\mathbf{v}}_{ab}\left[{\mathbf{B}}_a\cdot {\nabla }_a{W}_{ab}\left(h_a\right)\right]\hfill \\ \hfill & +\frac{1}{{\rho }_a}{\displaystyle \sum _b}{m}_b{\mathbf{B}}_a\left[{\mathbf{v}}_{ab}\cdot {\nabla }_a{W}_{ab}\left(h_a\right)\right].\hfill \end{array}(19)$

Numerically, any difference between Eqs 16, 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 Powell (1994); Powell et al. (1999). Consider the conservative formulation of the induction equation given by

$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{v}\times \mathbf{B}\right).(20)$

The vector calculus identity ∇ × (v ×B) ≡ (B ⋅∇)v − (∇ ⋅v)B − (v ⋅∇)B + (∇ ⋅B)v can be used to rewrite Eq. 20 in terms of a Lagrangian derivative, yielding

$\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}=\left(\mathbf{B}\cdot \nabla \right)\mathbf{v}-\left(\nabla \cdot \mathbf{v}\right)\mathbf{B}+\left(\nabla \cdot \mathbf{B}\right)\mathbf{v}.(21)$

This is identical to Eq. 18 except for the term proportional to ∇ ⋅B. In SPMHD, this term 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 Eq. 21 yields,

$\frac{\partial \left(\nabla \cdot \mathbf{B}\right)}{\partial t}=0,(22)$

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 Eq. 18 yields

$\frac{\partial \left(\nabla \cdot \mathbf{B}\right)}{\partial t}+\nabla \cdot \left(\mathbf{v}\nabla \cdot \mathbf{B}\right)=0.(23)$

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).

3.1.1 Example: Source term in the induction equation

The advection of a divergence ‘blob’ illustrates the difference between solving Eqs 18, 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/\sqrt{4\pi }$, with a perturbation $B_x=[{(r/r_0)}^8-2{(r/r_0)}^4+1]/\sqrt{4\pi }$ for r < r₀, with $r_0=1/\sqrt{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 (Eq. 18) and volume conservative (Eq. 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.

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.

3.2 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 and Monaghan (2004a), Price and Monaghan (2005), artificial resistivity adds dissipation to the magnetic field according to

$\frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathbf{B}}{\rho }\right)}_a={\displaystyle \sum _b}{m}_b\frac{v_{\text{sig,B}}}{{\bar{\rho }}_{ab}^2}\left({\mathbf{B}}_a-{\mathbf{B}}_b\right){\widehat{\mathbf{r}}}_{ab}\cdot \overline{{\nabla }_a{W}_{ab}},(24)$

where ${\bar{\rho }}_{ab}\equiv ({\rho }_a+{\rho }_b)/2$ and $\overline{{\nabla }_a{W}_{ab}}\equiv [{\nabla }_a{W}_{ab}(h_a)/{\mathrm{\Omega }}_a+{\nabla }_a{W}_{ab}(h_b)/{\mathrm{\Omega }}_b]/2$. The signal velocity, v_sig,B, represents the speed of information propagation between two particles.

Energy dissipated from the magnetic field may be added to the internal energy, u, through

$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=-{\displaystyle \sum _b}{m}_b\frac{v_{\text{sig,B}}}{{\bar{\rho }}_{ab}^2}{\left({\mathbf{B}}_a-{\mathbf{B}}_b\right)}^2{\widehat{\mathbf{r}}}_{ab}\cdot \overline{{\nabla }_a{W}_{ab}}.(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 and Monaghan, 2004a; Price and Monaghan, 2005).

One option for the signal velocity is

$v_{\text{sig,B}}=\frac{1}{2}\overline{{\alpha }_{\text{B},ab}}\left(v_{\text{mhd},a}+v_{\text{mhd},b}\right),(26)$

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 $\overline{{\alpha }_{\text{B},ab}}=({\alpha }_{\text{B},a}+{\alpha }_{\text{B},b})/2$. Price and Monaghan (2005) created a switch for α_B based on analogy to the Morris and Monaghan (1997) switch for artificial viscosity. In this case,

$\frac{\mathrm{d}{\alpha }_{\text{B},a}}{\mathrm{d}t}=\max \left(\frac{|\nabla \times {\mathbf{B}}_a|}{\sqrt{{\mu }_0{\rho }_a}},\frac{|\nabla \cdot {\mathbf{B}}_a|}{\sqrt{{\mu }_0{\rho }_a}}\right)-\frac{{\alpha }_{\text{B},a}-{\alpha }_{\text{B},0}}{\tau },(27)$

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⁶–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 in the signal velocity is similarly problematic if shocks are present. They recommend a new switch, setting

${\alpha }_{\text{B},a}=\frac{h_a|\nabla {\mathbf{B}}_a|}{|{\mathbf{B}}_a|}.(28)$

Note that this does not require the time evolution of α_B. On test problems, this switch led to less overall magnetic energy dissipation than the Price 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.

Price et al. (2018) advocated setting the signal velocity to

$v_{\text{sig},B}=|{\mathbf{v}}_{ab}\times {\widehat{\mathbf{r}}}_{ab}|.(29)$

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 Eq. 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 (resolution-dependent) physical dissipation term η∇²B with magnetic diffusivity

${\eta }_{\text{AR}}\sim \frac{1}{2}{v}_{\text{sig,B}}h.(30)$

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. Tricco et al. (2016b) found Pm ≈ 1 in simulations of supersonic, magnetized turbulence, and Wissing 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.

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.

3.2.1 Example: Brio-Wu shock tube

The accuracy of the artificial resistivity switch given by Eq. 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),Price and Monaghan (2004b); 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 [B_x, B_y, B_z]_L = [0.75, 1, 0] for the left state and [B_x, B_y, B_z]_R = [0.75, −1, 0] for the right state. 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⁴ 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 $\sim 10^{-3}$.

FIGURE 2

FIGURE 2. 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 Eq. 29.

3.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 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,

$\nabla {\alpha }_{\text{E}}=-\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b\left({\alpha }_{\text{E},a}-{\alpha }_{\text{E},b}\right){\nabla }_a{W}_{ab}\left(h_a\right),(31)$

$\nabla {\beta }_{\text{E}}=-\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b\left({\beta }_{\text{E},a}-{\beta }_{\text{E},b}\right){\nabla }_a{W}_{ab}\left(h_a\right).(32)$

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; Price and Bate, 2008; Price and Bate, 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 one-to-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).

Star formation calculations highlight the preceding limitations. SPMHD calculations that evolve the magnetic field via the induction equation produce magnetically-driven jets and outflows (e.g., Bürzle et al., 2011b; Price et al., 2012; Bate et al., 2014; Tsukamoto et al., 2015b; Lewis et al., 2015; Hopkins and Raives, 2016; Wurster et al., 2018a; Tsukamoto et al., 2018; Wissing and Shen, 2020; García-Senz et al., 2023). However, no such outflows are produced when similar calculations are performed using the Euler potentials (Price and Bate, 2007). At present, it seems there is little use for the Euler potentials in astrophysical simulation, except perhaps as a secondary check for numerical artefacts (Dolag and Stasyszyn, 2009).

3.4 Vector potential

Defining the magnetic field in terms of the vector potential, B = ∇ ×A, would seem an ideal choice. Doing so would guarantee a divergence-free 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., Eq. 22), which yields

$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}=\mathbf{v}\times \left(\nabla \times \mathbf{A}\right)+\left(\mathbf{v}\cdot \nabla \right)\mathbf{A}+\nabla \varphi ,(33)$

where ∇ϕ represents the gauge choice, arising as a constant of integration. Price (2010) used the gauge ϕ = −v ⋅A, yielding

$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}=-\mathbf{A}\times \left(\nabla \times \mathbf{v}\right)-\left(\mathbf{A}\cdot \nabla \right)\mathbf{v}.(34)$

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.

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.

4 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/μ₀ 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).

4.1 Conservative equations of motion

The SPMHD Lagrangian is

$L_{\text{SPMHD}}={\displaystyle \sum _a}{m}_a\left(\frac{v_a^2}{2}-u_a-\frac{B_a^2}{2{\mu }_0{\rho }_a}\right).(35)$

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 Eq. 6) and the specification on how the magnetic field is evolved. The induction equation given by Eq. 16 is assumed (or equivalently Eq. 19).

The resultant equations of motion from the Lagrangian are

$\begin{array}{ll}\hfill \frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}& =-{\displaystyle \sum _b}{m}_b\left[\frac{P_a}{{\mathrm{\Omega }}_a{\rho }_a^2}{\nabla }_a{W}_{ab}\left(h_a\right)+\frac{P_b}{{\mathrm{\Omega }}_b{\rho }_b^2}{\nabla }_a{W}_{ab}\left(h_b\right)\right]\hfill \\ \hfill & -\frac{1}{2{\mu }_0}{\displaystyle \sum _b}{m}_b\left[\frac{B_a^2}{{\mathrm{\Omega }}_a{\rho }_a^2}{\nabla }_a{W}_{ab}\left(h_a\right)+\frac{B_b^2}{{\mathrm{\Omega }}_b{\rho }_b^2}{\nabla }_a{W}_{ab}\left(h_b\right)\right]\hfill \\ \hfill & +\frac{1}{{\mu }_0}{\displaystyle \sum _b}{m}_b\left[\frac{{\mathbf{B}}_a}{{\mathrm{\Omega }}_a{\rho }_a^2}{\mathbf{B}}_a\cdot {\nabla }_a{W}_{ab}\left(h_a\right)\right.\hfill \\ \hfill & \left.+,\frac{{\mathbf{B}}_b}{{\mathrm{\Omega }}_b{\rho }_b^2},{\mathbf{B}}_b,\cdot {\nabla }_a{W}_{ab}\left(h_b\right)\right].\hfill \end{array}(36)$

The second term in Eq. 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 Eq. 36 with Eqs 2, 5, and, indeed, Eq. 36 could be written in terms of the stress tensor as

$\frac{\mathrm{d}{v}_a^i}{\mathrm{d}t}={\displaystyle \sum _b}{m}_b\left[\frac{S_a^{ij}}{{\mathrm{\Omega }}_a{\rho }_a^2}\frac{\partial W_{ab}\left(h_a\right)}{\partial x_a^j}+\frac{S_b^{ij}}{{\mathrm{\Omega }}_b{\rho }_b^2}\frac{\partial W_{ab}\left(h_b\right)}{\partial x_a^j}\right].(37)$

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

$\frac{\mathrm{d}{v}_a^i}{\mathrm{d}t}={\displaystyle \sum _b}{m}_b\frac{S_a^{ij}+S_b^{ij}}{{\rho }_a{\rho }_b}\overline{{\nabla }_a^j{W}_{ab}}.(38)$

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 magnetically-driven 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 Eq. 36 is equivalent to (B ⋅∇)B/μ₀ρ + (∇ ⋅B)B/μ₀ρ 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; Dehnen and Aly, 2012; McNally et al., 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/μ₀ρ 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 Eq. 36, this occurs for plasma β < 1 (Phillips and Monaghan, 1985; Børve et al., 2001; Bø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.

4.2 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 eight-wave 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,

$\begin{array}{ll}\hfill {\left(\frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}\right)}_{\text{divB}}=-{\widehat{\mathit{\beta }}}_a{\displaystyle \sum _b}{m}_b& \left[\frac{{\mathbf{B}}_a}{{\mathrm{\Omega }}_a{\rho }_a^2}\cdot {\nabla }_a{W}_{ab}\left(h_a\right)\right.\hfill \\ \hfill & \left.+\frac{{\mathbf{B}}_b}{{\mathrm{\Omega }}_b{\rho }_b^2}\cdot {\nabla }_a{W}_{ab}\left(h_b\right)\right],\hfill \end{array}(39)$

where ${\widehat{\mathit{\beta }}}_a$ used to regulate the strength of the applied correction, with ${\widehat{\mathit{\beta }}}_a={\mathbf{B}}_a$ corresponding to a full correction of the force term. Adding Eq. 39 to the conservative momentum equation (Eq. 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 $\widehat{\mathit{\beta }}=\frac{1}{2}\mathbf{B}$, thereby subtracting only $\frac{1}{2}(\nabla \cdot \mathbf{B})\mathbf{B}/{\mu }_0\rho$. 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 $({\widehat{\mathit{\beta }}}_a={\mathbf{B}}_a)$ causes dispersive errors in slow magnetosonic waves. The half strength term $\left({\widehat{\mathit{\beta }}}_a=\frac{1}{2}{\mathbf{B}}_a\right)$ 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,

${\widehat{\mathit{\beta }}}_a=\left\{\begin{array}{cc}{\mathbf{B}}_a\hfill & \beta <2\hfill \\ {\mathbf{B}}_a\left(10-\beta \right)/8\hfill & 2<\beta <10\hfill \\ 0\hfill & \beta \ge 10\hfill \end{array}\right.(40)$

and

${\widehat{\mathit{\beta }}}_a=\left\{\begin{array}{cc}{\mathbf{B}}_a\hfill & \beta <1\hfill \\ {\mathbf{B}}_a\left(2-\beta \right)\hfill & 1<\beta <2\hfill \\ 0\hfill & \beta \ge 2\hfill \end{array}\right.(41)$

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 Eq. 36), but replaces the anisotropic force (third term in Eq. 36) with

${\left(\frac{\mathrm{d}{v}_a^i}{\mathrm{d}t}\right)}_{\text{aniso}}=\frac{1}{{\mu }_0}{\displaystyle \sum _b}{m}_b\frac{B_b^i{B}_b^j-B_a^i{B}_a^j}{{\rho }_a{\rho }_b}\overline{{\nabla }_a^j{W}_{ab}\left(h_a\right)}.(42)$

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).

4.2.1 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 $\widehat{\mathit{\beta }}=\mathbf{B}$ (top panel) and $\widehat{\mathit{\beta }}=\frac{1}{2}\mathbf{B}$ (bottom panel), where $\widehat{\mathit{\beta }}$ is the leading term in the tensile instability correction term given by Eq. 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 $\widehat{\mathit{\beta }}=\mathbf{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 $\widehat{\mathit{\beta }}=\frac{1}{2}\mathbf{B}$, there are no dispersive errors and the wave propagates (mostly) as expected.

FIGURE 3

FIGURE 3. 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, $\widehat{\mathit{\beta }}=\mathbf{B}$ and $\widehat{\mathit{\beta }}=\frac{1}{2}\mathbf{B}$ as used in Eq. 39. Applying the full correction term leads to supersonic dispersive errors, whereas the half correction term yields results consistent with the expected solution.

5 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 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 Eq. 39, whereby

$\nabla \cdot {\mathbf{B}}_a={\rho }_a{\displaystyle \sum _b}{m}_b\left[\frac{{\mathbf{B}}_a}{{\mathrm{\Omega }}_a{\rho }_a^2}\cdot {\nabla }_a{W}_{ab}\left(h_a\right)+\frac{{\mathbf{B}}_b}{{\mathrm{\Omega }}_b{\rho }_b^2}\cdot {\nabla }_a{W}_{ab}\left(h_b\right)\right].(43)$

Instead of this symmetric gradient estimate, one could instead measure the divergence of the magnetic field with a difference gradient estimate according to

$\nabla \cdot {\mathbf{B}}_a=-\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b\left({\mathbf{B}}_a-{\mathbf{B}}_b\right)\cdot {\nabla }_a{W}_{ab}\left(h_a\right).(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 Section 3.3 and Section 3.4. Instead, it would appear that the best option currently is to only approximately uphold this constraint.

5.1 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

$\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}=-\nabla \psi ,(45)$

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\psi }{c_{\text{h}}}\right)=-c_{\text{h}}\left(\nabla \cdot \mathbf{B}\right)-\frac{1}{\tau }\left(\frac{\psi }{c_{\text{h}}}\right)-\frac{1}{2}\left(\frac{\psi }{c_{\text{h}}}\right)\left(\nabla \cdot \mathbf{v}\right).(46)$

This formulation differs from Dedner et al. (2002), in that Eqs 45, 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

$\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{\sqrt{\rho }{c}_{\text{h}}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\psi }{\sqrt{\rho }{c}_{\text{h}}}\right)\right]-\frac{{\nabla }^2\psi }{\rho }+\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{\sqrt{\rho }{c}_{\text{h}}}\left(\frac{\psi }{\sqrt{\rho }{c}_{\text{h}}\tau }\right)\right]=0.(47)$

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,

$\frac{{\partial }^2\psi }{\partial t^2}-c_{\text{h}}^2{\nabla }^2\psi +\frac{1}{\tau }\frac{\partial \psi }{\partial t}=0.(48)$

The SPMHD discretised cleaning equations are

$\frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathbf{B}}{\rho }\right)}_a=-{\displaystyle \sum _b}{m}_b\left[\frac{{\psi }_a}{{\mathrm{\Omega }}_a{\rho }_a^2}{\nabla }_a{W}_{ab}\left(h_a\right)+\frac{{\psi }_b}{{\mathrm{\Omega }}_b{\rho }_b^2}{\nabla }_a{W}_{ab}\left(h_b\right)\right],(49)$

$\begin{array}{ll}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\psi }{c_{\text{h}}}\right)}_a& =\frac{c_{\text{h},a}}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b\left({\mathbf{B}}_a-{\mathbf{B}}_b\right)\cdot {\nabla }_a{W}_{ab}\left(h_a\right)-\frac{1}{\tau }{\left(\frac{\psi }{c_{\text{h}}}\right)}_a\hfill \\ \hfill & +\frac{1}{2}{\left(\frac{\psi }{c_{\text{h}}}\right)}_a{\displaystyle \sum _b}{m}_b{\mathbf{v}}_{ab}\cdot {\nabla }_a{W}_{ab}\left(h_a\right).\hfill \end{array}(50)$

Note that the difference estimate is used for ∇ ⋅B (Eq. 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_{\psi }={\psi }^2/2{\mu }_0\rho c_{\text{h}}^2$, 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

$\tau =\frac{h}{\sigma c_{\text{h}}},(51)$

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

${ϵ}_{\text{divB}}=\frac{h|\nabla \cdot \mathbf{B}|}{|\mathbf{B}|}.(52)$

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.

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 5

FIGURE 5. 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).

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 non-conservative 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 Eq. 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 (Eq. 43). Though numerically stable, this low-order estimate was found to be overly dissipative and could introduce artefacts into the physical portions of the magnetic field since this estimate is sensitive to particle disorder.

5.1.1 Example: 3D Orszag-Tang vortex

The Orszag-Tang vortex (Orszag and Tang, 1979) is a widely used test that consists of interacting vortices that incite magnetized turbulence, and is used here to demonstrate the efficacy of mixed hyperbolic/parabolic divergence cleaning. Results for this test have been shown for SPMHD many times, for example, Price and Monaghan (2005); Børve et al. (2006); Rosswog and Price (2007); Dolag and Stasyszyn (2009); Barnes et al. (2012); Tricco and Price (2012), Tricco and Price (2013); Price et al. (2018); Wissing and Shen (2020). Originally a 2D test, it is extended to 3D for this test by adding a small thickness in z. The simulation domain is x, y ∈ [0, 1] and z ∈ [0, 3/128], with 512 × 512×12 particles arranged on a cubic lattice. The initial conditions are ρ = 25/(36 π), p = 5/(12π), v = [-sin(2πy), sin(2πx), 0] and B = [-sin(2πy), sin(4πx), 0] with γ = 5/3. The calculations use the artificial resistivity switch given by Eq. 29.

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 25th to 75th percentiles. There is a significant reduction in divergence error when divergence cleaning is applied. The median ϵ_divB with divergence cleaning is $\sim 0.1\%$ (the average ϵ_divB is also $\sim 0.1\%$). This is below the 25th percentile without divergence cleaning.

FIGURE 6

FIGURE 6. Density cross-sections of the 3D Orszag-Tang vortex.

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 25th and 75th percentiles.

5.2 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, η_AR∇²B = η_AR∇(∇ ⋅B) − η_AR∇ × (∇ ×B). 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 Section 3.3 and Section 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 the field (since the divergence of the curl is zero). From this, it can be stated that ∇ ⋅B* = ∇²ϕ, 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.

6 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), Braiding and Wardle (2012b); 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., 2016; Wurster 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

$\begin{array}{ll}\hfill \frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}& =-\nabla \times \left[{\eta }_{\text{O}}\left(\nabla \times \mathbf{B}\right)\right]\hfill \\ \hfill & +\nabla \times \left\{{\eta }_{\text{A}}\left[\left(\nabla \times \mathbf{B}\right)\times \widehat{\mathbf{B}}\right]\times \widehat{\mathbf{B}}\right\}\hfill \\ \hfill & -\nabla \times \left[{\eta }_{\text{H}}\left(\nabla \times \mathbf{B}\right)\times \widehat{\mathbf{B}}\right],\hfill \end{array}(53)$

where $\widehat{\mathbf{B}}=\mathbf{B}/|\mathbf{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

$\delta t=\frac{C_{\eta }{h}^2}{\eta },(54)$

which is calculated for each non-ideal effect. Price et al. (2018) use C_η = 1/(2π). This timestep constraint can become dominant for large η, and additionally note that it scales ∝ h², whereas the Courant condition is only ∝ h. Super-time-stepping may be of some benefit to alleviate pressure from the non-ideal MHD timestep constraints (Alexiades et al., 1996; Tsukamoto et al., 2013; Wurster et al., 2016), though this is only applicable for Ohmic dissipation and ambipolar diffusion.

6.1 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

$\nabla \times {\mathbf{B}}_a=\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}{\displaystyle \sum _b}{m}_b\left({\mathbf{B}}_a-{\mathbf{B}}_b\right)\times \nabla W_{ab}\left(h_a\right),(55)$

$\begin{array}{ll}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathbf{B}}{\rho }\right)}_a=-{\displaystyle \sum _b}& \left[\frac{{\eta }_{\text{O},a}{\left(\nabla \times \mathbf{B}\right)}_a}{{\mathrm{\Omega }}_a{\rho }_a^2}\times {\nabla }_a{W}_{ab}\left(h_a\right)\right.\hfill \\ \hfill & \left.+\frac{{\eta }_{\text{O},b}{\left(\nabla \times \mathbf{B}\right)}_b}{{\mathrm{\Omega }}_b{\rho }_b^2}\times {\nabla }_a{W}_{ab}\left(h_b\right)\right].\hfill \end{array}(56)$

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

$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{{\eta }_{\text{O},a}{\left(\nabla \times \mathbf{B}\right)}_a^2}{{\rho }_a}(57)$

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).

6.2 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 (Eq. 55). The second derivative is estimated with a symmetric derivative, according to

$\begin{array}{ll}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathbf{B}}{\rho }\right)}_a=-{\displaystyle \sum _b}& \left[\frac{{\mathbf{D}}_{\text{A},a}}{{\mathrm{\Omega }}_a{\rho }_a^2}\times {\nabla }_a{W}_{ab}\left(h_a\right)\right.\hfill \\ \hfill & \left.+\frac{{\mathbf{D}}_{\text{A},b}}{{\mathrm{\Omega }}_b{\rho }_b^2}\times {\nabla }_a{W}_{ab}\left(h_b\right)\right],\hfill \end{array}(58)$

with

${\mathbf{D}}_{\text{A},a}=-{\eta }_{\text{A},a}\left[{\left(\nabla \times \mathbf{B}\right)}_a\times {\widehat{\mathbf{B}}}_a\right]\times {\widehat{\mathbf{B}}}_a.(59)$

A positive-definite increase in thermal energy is ensured (Wurster et al., 2014), with

$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{{\eta }_{\text{A},a}}{{\rho }_a}\left\{{\left(\nabla \times \mathbf{B}\right)}_a^2-{\left[{\left(\nabla \times \mathbf{B}\right)}_a\cdot {\widehat{\mathbf{B}}}_a\right]}^2\right\}.(60)$

6.3 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 Eq. 55. The outer derivative is calculated using a symmetric derivative estimate according to Eq. 58, but with

${\mathbf{D}}_{\text{H},a}={\eta }_{\text{H},a}\left(\nabla \times {\mathbf{B}}_a\right)\times {\widehat{\mathbf{B}}}_a(61)$

in place of D_A.

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 asymmetric in the presence of the Hall effect, that is, they have different phase velocities. The right-handed 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.

7 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 ($\sim 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 zero-net 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.

FIGURE 8

FIGURE 8. Spacetime diagrams of magneto-rotational 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).

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, particle-based 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 Godunov-type 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 constrained-gradient 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., 2014; Hopkins et al., 2018; Hopkins 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 gas-dust mixtures in SPH (Laibe and Price, 2012a; Laibe and Price, 2012b; 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. (2021), Tsukamoto 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.

Author contributions

TT: Writing–original draft.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. TT was supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. This research was enabled in part by support provided by ACENET and the Digital Research Alliance of Canada.

Acknowledgments

TT thanks Daniel Price for many conversations over the years that has refined his understanding of, and appreciation for, SPH. TT thanks Robert Wissing for permission to reproduce the MRI 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).

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Agertz, O., Moore, B., Stadel, J., Potter, D., Miniati, F., Read, J., et al. (2007). Fundamental differences between SPH and grid methods. MNRAS 380, 963–978. doi:10.1111/j.1365-2966.2007.12183.x

CrossRef Full Text | Google Scholar

Alexiades, V., Amiez, G., and Gremaud, P.-A. (1996). Super-time-stepping acceleration of explicit schemes for parabolic problems. Commun. Numer. Methods Eng. 12, 31–42. doi:10.1002/(SICI)1099-0887(199601)12:1⟨31::AID-CNM950⟩3.0.CO;2-5

CrossRef Full Text | Google Scholar

Artymowicz, P., and Lubow, S. H. (1994). Dynamics of binary-disk interaction. 1: resonances and disk gap sizes. ApJ 421, 651. doi:10.1086/173679

CrossRef Full Text | Google Scholar

Bai, X.-N. (2014). Hall-effect-Controlled gas dynamics in protoplanetary disks. I. Wind solutions at the inner disk. ApJ 791, 137. doi:10.1088/0004-637X/791/2/137

CrossRef Full Text | Google Scholar

Barnes, D. J., Kawata, D., and Wu, K. (2012). Cosmological simulations using GCMHD+. MNRAS 420, 3195–3212. doi:10.1111/j.1365-2966.2011.20247.x

CrossRef Full Text | Google Scholar

Barnes, D. J., On, A. Y. L., Wu, K., and Kawata, D. (2018). SPMHD simulations of structure formation. MNRAS 476, 2890–2904. doi:10.1093/mnras/sty400

CrossRef Full Text | Google Scholar

Bate, M. R., Tricco, T. S., and Price, D. J. (2014). Collapse of a molecular cloud core to stellar densities: stellar-core and outflow formation in radiation magnetohydrodynamic simulations. MNRAS 437, 77–95. doi:10.1093/mnras/stt1865

CrossRef Full Text | Google Scholar

Biriukov, S., and Price, D. J. (2019). Stable anisotropic heat conduction in smoothed particle hydrodynamics. MNRAS 483, 4901–4909. doi:10.1093/mnras/sty3413

CrossRef Full Text | Google Scholar

Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., and Borgani, S. (2011). A non-ideal magnetohydrodynamic GADGET: simulating massive galaxy clusters. MNRAS 418, 2234–2250. doi:10.1111/j.1365-2966.2011.19523.x

CrossRef Full Text | Google Scholar

Bonnerot, C., Price, D. J., Lodato, G., and Rossi, E. M. (2017). Magnetic field evolution in tidal disruption events. MNRAS 469, 4879–4888. doi:10.1093/mnras/stx1210

CrossRef Full Text | Google Scholar

Børve, S., Omang, M., and Trulsen, J. (2001). Regularized smoothed particle hydrodynamics: a new approach to simulating magnetohydrodynamic shocks. ApJ 561, 82–93. doi:10.1086/323228

CrossRef Full Text | Google Scholar

Børve, S., Omang, M., and Trulsen, J. (2004). Two-dimensional MHD smoothed particle hydrodynamics stability analysis. ApJS 153, 447–462. doi:10.1086/421520

CrossRef Full Text | Google Scholar

Børve, S., Omang, M., and Trulsen, J. (2006). Multidimensional MHD shock tests of regularized smoothed particle hydrodynamics. ApJ 652, 1306–1317. doi:10.1086/508454

CrossRef Full Text | Google Scholar

Brackbill, J. U., and Barnes, D. C. (1980). The Effect of Nonzero ∇ ⋅ B on the numerical solution of the magnetohydrodynamic equations. J. Comput. Phys. 35, 426–430. doi:10.1016/0021-9991(80)90079-0

CrossRef Full Text | Google Scholar

Braiding, C. R., and Wardle, M. (2012a). The Hall effect in accretion flows. MNRAS 427, 3188–3195. doi:10.1111/j.1365-2966.2012.22001.x

CrossRef Full Text | Google Scholar

Braiding, C. R., and Wardle, M. (2012b). The Hall effect in star formation. MNRAS 422, 261–281. doi:10.1111/j.1365-2966.2012.20601.x

CrossRef Full Text | Google Scholar

Brandenburg, A. (2010). Magnetic field evolution in simulations with Euler potentials. MNRAS 401, 347–354. doi:10.1111/j.1365-2966.2009.15640.x

CrossRef Full Text | Google Scholar

Brio, M., and Wu, C. C. (1988). An upwind differencing scheme for the equations of ideal magnetohydrodynamics. J. Comput. Phys. 75, 400–422. doi:10.1016/0021-9991(88)90120-9

CrossRef Full Text | Google Scholar

Brookshaw, L. (1985). A method of calculating radiative heat diffusion in particle simulations. PASA 6, 207–210. doi:10.1017/S1323358000018117

CrossRef Full Text | Google Scholar

Bürzle, F., Clark, P. C., Stasyszyn, F., Dolag, K., and Klessen, R. S. (2011a). Protostellar outflows with smoothed particle magnetohydrodynamics. MNRAS 417, L61–L65. doi:10.1111/j.1745-3933.2011.01120.x

CrossRef Full Text | Google Scholar

Bürzle, F., Clark, P. C., Stasyszyn, F., Greif, T., Dolag, K., Klessen, R. S., et al. (2011b). Protostellar collapse and fragmentation using an MHD GADGET. MNRAS 412, 171–186. doi:10.1111/j.1365-2966.2010.17896.x

CrossRef Full Text | Google Scholar

Cleary, P. W., and Monaghan, J. J. (1999). Conduction modelling using smoothed particle hydrodynamics. J. Comput. Phys. 148, 227–264. doi:10.1006/jcph.1998.6118

CrossRef Full Text | Google Scholar

Cummins, S. J., and Rudman, M. (1999). An SPH projection method. J. Comput. Phys. 152, 584–607. doi:10.1006/jcph.1999.6246

CrossRef Full Text | Google Scholar

Dedner, A., Kemm, F., Kröner, D., Munz, C. D., Schnitzer, T., and Wesenberg, M. (2002). Hyperbolic divergence cleaning for the MHD equations. J. Comput. Phys. 175, 645–673. doi:10.1006/jcph.2001.6961

CrossRef Full Text | Google Scholar

Dehnen, W., and Aly, H. (2012). Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. MNRAS 425, 1068–1082. doi:10.1111/j.1365-2966.2012.21439.x

CrossRef Full Text | Google Scholar

Deng, H., Mayer, L., Latter, H., Hopkins, P. F., and Bai, X.-N. (2019). Local simulations of MRI turbulence with meshless methods. ApJS 241, 26. doi:10.3847/1538-4365/ab0957

CrossRef Full Text | Google Scholar

Dobbs, C. L., and Price, D. J. (2008). Magnetic fields and the dynamics of spiral galaxies. MNRAS 383, 497–512. doi:10.1111/j.1365-2966.2007.12591.x

CrossRef Full Text | Google Scholar

Dobbs, C. L., Price, D. J., Pettitt, A. R., Bate, M. R., and Tricco, T. S. (2016). Magnetic field evolution and reversals in spiral galaxies. MNRAS 461, 4482–4495. doi:10.1093/mnras/stw1625

CrossRef Full Text | Google Scholar

Dobbs, C. L., and Wurster, J. (2021). The properties of clusters, and the orientation of magnetic fields relative to filaments, in magnetohydrodynamic simulations of colliding clouds. MNRAS 502, 2285–2295. doi:10.1093/mnras/stab150

CrossRef Full Text | Google Scholar

Dolag, K., and Stasyszyn, F. (2009). An MHD GADGET for cosmological simulations. MNRAS 398, 1678–1697. doi:10.1111/j.1365-2966.2009.15181.x

CrossRef Full Text | Google Scholar

Evans, C. R., and Hawley, J. F. (1988). Simulation of magnetohydrodynamic flows: a constrained transport model. ApJ 332, 659. doi:10.1086/166684

CrossRef Full Text | Google Scholar

Forgan, D., Price, D. J., and Bonnell, I. (2017). On the fragmentation boundary in magnetized self-gravitating discs. MNRAS 466, 3406–3416. doi:10.1093/mnras/stw3314

CrossRef Full Text | Google Scholar

García-Senz, D., Cabezón, R. M., and Escartín, J. A. (2012). Improving smoothed particle hydrodynamics with an integral approach to calculating gradients. A&A 538, A9. doi:10.1051/0004-6361/201117939

CrossRef Full Text | Google Scholar

García-Senz, D., Cabezón, R. M., and Escartín, J. A. (2022). Conservative, density-based smoothed particle hydrodynamics with improved partition of the unity and better estimation of gradients. A&A 659, A175. doi:10.1051/0004-6361/202141877

CrossRef Full Text | Google Scholar

García-Senz, D., Wissing, R., Cabezón, R. M., Vurgun, E., and Linares, M. (2023). Axisymmetric smoothed particle magnetohydrodynamics. MNRAS 518, 4115–4131. doi:10.1093/mnras/stac3328

CrossRef Full Text | Google Scholar

Gingold, R. A., and Monaghan, J. J. (1977). Smoothed particle hydrodynamics: theory and application to non-spherical stars. MNRAS 181, 375–389. doi:10.1093/mnras/181.3.375

CrossRef Full Text | Google Scholar

Harris, A., and Tricco, T. (2023). Sarracen: a Python package for analysis and visualization of smoothed particle hydrodynamics data. J. Open Source Softw. 8, 5263. doi:10.21105/joss.05263

CrossRef Full Text | Google Scholar

Hennebelle, P., and Inutsuka, S.-i. (2019). The role of magnetic field in molecular cloud formation and evolution. Front. Astronomy Space Sci. 6, 5. doi:10.3389/fspas.2019.00005

CrossRef Full Text | Google Scholar

Herrington, N. P., Dobbs, C. L., and Bending, T. J. R. (2023). The role of previous generations of stars in triggering star formation and driving gas dynamics. MNRAS 521, 5712–5723. doi:10.1093/mnras/stad923

CrossRef Full Text | Google Scholar

Hopkins, P. F. (2013). A general class of Lagrangian smoothed particle hydrodynamics methods and implications for fluid mixing problems. MNRAS 428, 2840–2856. doi:10.1093/mnras/sts210

CrossRef Full Text | Google Scholar

Hopkins, P. F. (2015). A new class of accurate, mesh-free hydrodynamic simulation methods. MNRAS 450, 53–110. doi:10.1093/mnras/stv195

CrossRef Full Text | Google Scholar

Hopkins, P. F. (2016). A constrained-gradient method to control divergence errors in numerical MHD. MNRAS 462, 576–587. doi:10.1093/mnras/stw1578

CrossRef Full Text | Google Scholar

Hopkins, P. F., Kereš, D., Oñorbe, J., Faucher-Giguère, C.-A., Quataert, E., Murray, N., et al. (2014). Galaxies on FIRE (Feedback in Realistic Environments): stellar feedback explains cosmologically inefficient star formation. MNRAS 445, 581–603. doi:10.1093/mnras/stu1738

CrossRef Full Text | Google Scholar

Hopkins, P. F., and Raives, M. J. (2016). Accurate, meshless methods for magnetohydrodynamics. MNRAS 455, 51–88. doi:10.1093/mnras/stv2180

CrossRef Full Text | Google Scholar

Hopkins, P. F., Wetzel, A., Kereš, D., Faucher-Giguère, C.-A., Quataert, E., Boylan-Kolchin, M., et al. (2018). FIRE-2 simulations: physics versus numerics in galaxy formation. MNRAS 480, 800–863. doi:10.1093/mnras/sty1690

CrossRef Full Text | Google Scholar

Hopkins, P. F., Wetzel, A., Wheeler, C., Sanderson, R., Grudić, M. Y., Sameie, O., et al. (2023). FIRE-3: updated stellar evolution models, yields, and microphysics and fitting functions for applications in galaxy simulations. MNRAS 519, 3154–3181. doi:10.1093/mnras/stac3489

CrossRef Full Text | Google Scholar

Hosking, J. G., and Whitworth, A. P. (2004). Modelling ambipolar diffusion with two-fluid smoothed particle hydrodynamics. MNRAS 347, 994–1000. doi:10.1111/j.1365-2966.2004.07273.x

CrossRef Full Text | Google Scholar

Hu, X. Y., and Adams, N. A. (2007). An incompressible multi-phase SPH method. J. Comput. Phys. 227, 264–278. doi:10.1016/j.jcp.2007.07.013

CrossRef Full Text | Google Scholar

Iwasaki, K. (2015). Minimizing dispersive errors in smoothed particle magnetohydrodynamics for strongly magnetized medium. J. Comput. Phys. 302, 359–373. doi:10.1016/j.jcp.2015.09.022

CrossRef Full Text | Google Scholar

Iwasaki, K., and Inutsuka, S.-I. (2011). Smoothed particle magnetohydrodynamics with a Riemann solver and the method of characteristics. MNRAS 418, 1668–1688. doi:10.1111/j.1365-2966.2011.19588.x

CrossRef Full Text | Google Scholar

Kotarba, H., Lesch, H., Dolag, K., Naab, T., Johansson, P. H., and Stasyszyn, F. A. (2009). Magnetic field structure due to the global velocity field in spiral galaxies. MNRAS 397, 733–747. doi:10.1111/j.1365-2966.2009.15030.x

CrossRef Full Text | Google Scholar

Laibe, G., and Price, D. J. (2012a). Dusty gas with smoothed particle hydrodynamics - I. Algorithm and test suite. MNRAS 420, 2345–2364. doi:10.1111/j.1365-2966.2011.20202.x

CrossRef Full Text | Google Scholar

Laibe, G., and Price, D. J. (2012b). Dusty gas with smoothed particle hydrodynamics - II. Implicit timestepping and astrophysical drag regimes. MNRAS 420, 2365–2376. doi:10.1111/j.1365-2966.2011.20201.x

CrossRef Full Text | Google Scholar

Laibe, G., and Price, D. J. (2014). Dusty gas with one fluid in smoothed particle hydrodynamics. MNRAS 440, 2147–2163. doi:10.1093/mnras/stu359

CrossRef Full Text | Google Scholar

Lewis, B. T., and Bate, M. R. (2017). The dependence of protostar formation on the geometry and strength of the initial magnetic field. MNRAS 467, 3324–3337. doi:10.1093/mnras/stx271

CrossRef Full Text | Google Scholar

Lewis, B. T., Bate, M. R., and Price, D. J. (2015). Smoothed particle magnetohydrodynamic simulations of protostellar outflows with misaligned magnetic field and rotation axes. MNRAS 451, 288–299. doi:10.1093/mnras/stv957

CrossRef Full Text | Google Scholar

Lind, S. J., Rogers, B. D., and Stansby, P. K. (2020). Review of smoothed particle hydrodynamics: towards converged Lagrangian flow modelling. Proc. R. Soc. Lond. Ser. A 476, 20190801. doi:10.1098/rspa.2019.0801

PubMed Abstract | CrossRef Full Text | Google Scholar

Lodato, G., and Price, D. J. (2010). On the diffusive propagation of warps in thin accretion discs. MNRAS 405, 1212–1226. doi:10.1111/j.1365-2966.2010.16526.x

CrossRef Full Text | Google Scholar

Lucy, L. B. (1977). A numerical approach to the testing of the fission hypothesis. AJ 82, 1013–1024. doi:10.1086/112164

CrossRef Full Text | Google Scholar

Maury, A., Hennebelle, P., and Girart, J. M. (2022). Recent progress with observations and models to characterize the magnetic fields from star-forming cores to protostellar disks. Front. Astronomy Space Sci. 9, 949223. doi:10.3389/fspas.2022.949223

CrossRef Full Text | Google Scholar

McNally, C. P., Lyra, W., and Passy, J.-C. (2012). A well-posed kelvin-helmholtz instability test and comparison. ApJS 201, 18. doi:10.1088/0067-0049/201/2/18

CrossRef Full Text | Google Scholar

Meglicki, Z., Wickramasinghe, D., and Dewar, R. L. (1995). Gravitational collapse of a magnetized vortex: application to the Galactic Centre. MNRAS 272, 717–729. doi:10.1093/mnras/272.4.717

CrossRef Full Text | Google Scholar

Monaghan, J. J. (2002). SPH compressible turbulence. MNRAS 335, 843–852. doi:10.1046/j.1365-8711.2002.05678.x

CrossRef Full Text | Google Scholar

Monaghan, J. J. (2005). Smoothed particle hydrodynamics. Rep. Prog. Phys. 68, 1703–1759. doi:10.1088/0034-4885/68/8/R01

CrossRef Full Text | Google Scholar

Morris, J. P. (1996). Analysis of smoothed particle hydrodynamics with applications. Ph.D. thesis. Melbourne, Australia: Monash University.

Google Scholar

Morris, J. P., and Monaghan, J. J. (1997). A switch to reduce SPH viscosity. J. Comput. Phys. 136, 41–50. doi:10.1006/jcph.1997.5690

CrossRef Full Text | Google Scholar

Orszag, S. A., and Tang, C. M. (1979). Small-scale structure of two-dimensional magnetohydrodynamic turbulence. J. Fluid Mech. 90, 129–143. doi:10.1017/S002211207900210X

CrossRef Full Text | Google Scholar

Pandey, B. P., and Wardle, M. (2008). Hall magnetohydrodynamics of partially ionized plasmas. MNRAS 385, 2269–2278. doi:10.1111/j.1365-2966.2008.12998.x

CrossRef Full Text | Google Scholar

Park, S.-S., Kim, D.-K., Kim, J.-H., and Kim, E. S. (2023). Smoothed particle hydrodynamics method for pinch plasma simulation with non-ideal MHD model. Phys. Plasmas 30, 053901. doi:10.1063/5.0138221

CrossRef Full Text | Google Scholar

Phillips, G. J. (1986a). Three-dimensional numerical simulations of collapsing, isothermal magnetic gas clouds. MNRAS 221, 571–587. doi:10.1093/mnras/221.3.571

CrossRef Full Text | Google Scholar

Phillips, G. J. (1986b). Three-dimensional numerical simulations of collapsing isothermal magnetic gas clouds - non-uniform initial fields. MNRAS 222, 111–119. doi:10.1093/mnras/222.1.111

CrossRef Full Text | Google Scholar

Phillips, G. J., and Monaghan, J. J. (1985). A numerical method for three-dimensional simulations of collapsing, isothermal, magnetic gas clouds. MNRAS 216, 883–895. doi:10.1093/mnras/216.4.883

CrossRef Full Text | Google Scholar

Powell, K. G. (1994). Approximate Riemann solver for magnetohydrodynamics (that works in more than one dimension).

Google Scholar

Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., and De Zeeuw, D. L. (1999). A solution-adaptive upwind scheme for ideal magnetohydrodynamics. J. Comput. Phys. 154, 284–309. doi:10.1006/jcph.1999.6299

CrossRef Full Text | Google Scholar

Price, D. J. (2008). Modelling discontinuities and kelvin helmholtz instabilities in SPH. J. Comput. Phys. 227, 10040–10057. doi:10.1016/j.jcp.2008.08.011

CrossRef Full Text | Google Scholar

Price, D. J. (2010). Smoothed particle magnetohydrodynamics - IV. Using the vector potential. MNRAS 401, 1475–1499. doi:10.1111/j.1365-2966.2009.15763.x

CrossRef Full Text | Google Scholar

Price, D. J. (2012). Smoothed particle hydrodynamics and magnetohydrodynamics. J. Comput. Phys. 231, 759–794. doi:10.1016/j.jcp.2010.12.011

CrossRef Full Text | Google Scholar

Price, D. J., and Bate, M. R. (2007). The impact of magnetic fields on single and binary star formation. MNRAS 377, 77–90. doi:10.1111/j.1365-2966.2007.11621.x

CrossRef Full Text | Google Scholar

Price, D. J., and Bate, M. R. (2008). The effect of magnetic fields on star cluster formation. MNRAS 385, 1820–1834. doi:10.1111/j.1365-2966.2008.12976.x

CrossRef Full Text | Google Scholar

Price, D. J., and Bate, M. R. (2009). Inefficient star formation: the combined effects of magnetic fields and radiative feedback. MNRAS 398, 33–46. doi:10.1111/j.1365-2966.2009.14969.x

CrossRef Full Text | Google Scholar

Price, D. J., and Laibe, G. (2015). A fast and explicit algorithm for simulating the dynamics of small dust grains with smoothed particle hydrodynamics. MNRAS 451, 813–826. doi:10.1093/mnras/stv996

CrossRef Full Text | Google Scholar

Price, D. J., and Monaghan, J. J. (2004a). Smoothed Particle Magnetohydrodynamics - I. Algorithm and tests in one dimension. MNRAS 348, 123–138. doi:10.1111/j.1365-2966.2004.07345.x

CrossRef Full Text | Google Scholar

Price, D. J., and Monaghan, J. J. (2004b). Smoothed Particle Magnetohydrodynamics - II. Variational principles and variable smoothing-length terms. MNRAS 348, 139–152. doi:10.1111/j.1365-2966.2004.07346.x

CrossRef Full Text | Google Scholar

Price, D. J., and Monaghan, J. J. (2005). Smoothed particle magnetohydrodynamics — III. Multidimensional tests and the ∇·B= 0 constraint. MNRAS 364, 384–406. doi:10.1111/j.1365-2966.2005.09576.x

CrossRef Full Text | Google Scholar

Price, D. J., and Rosswog, S. (2006). Producing ultrastrong magnetic fields in neutron star mergers. Science 312, 719–722. doi:10.1126/science.1125201

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, D. J., Tricco, T. S., and Bate, M. R. (2012). Collimated jets from the first core. MNRAS 423, L45–L49. doi:10.1111/j.1745-3933.2012.01254.x

CrossRef Full Text | Google Scholar

Price, D. J., Wurster, J., Tricco, T. S., Nixon, C., Toupin, S., Pettitt, A., et al. (2018). Phantom: a smoothed particle hydrodynamics and magnetohydrodynamics code for astrophysics. PASA 35, e031. doi:10.1017/pasa.2018.25

CrossRef Full Text | Google Scholar

Pudritz, R. E., and Ray, T. P. (2019). The role of magnetic fields in protostellar outflows and star formation. Front. Astronomy Space Sci. 6, 54. doi:10.3389/fspas.2019.00054

CrossRef Full Text | Google Scholar

Read, J. I., Hayfield, T., and Agertz, O. (2010). Resolving mixing in smoothed particle hydrodynamics. MNRAS 405, 1513–1530. doi:10.1111/j.1365-2966.2010.16577.x

CrossRef Full Text | Google Scholar

Ritchie, B. W., and Thomas, P. A. (2001). Multiphase smoothed-particle hydrodynamics. MNRAS 323, 743–756. doi:10.1046/j.1365-8711.2001.04268.x

CrossRef Full Text | Google Scholar

Rosswog, S. (2009). Astrophysical smooth particle hydrodynamics. New A Rev. 53, 78–104. doi:10.1016/j.newar.2009.08.007

CrossRef Full Text | Google Scholar

Rosswog, S., and Price, D. (2007). MAGMA: a three-dimensional, Lagrangian magnetohydrodynamics code for merger applications. MNRAS 379, 915–931. doi:10.1111/j.1365-2966.2007.11984.x

CrossRef Full Text | Google Scholar

Ryu, D., and Jones, T. W. (1995). Numerical magetohydrodynamics in astrophysics: algorithm and tests for one-dimensional flow`. ApJ 442, 228. doi:10.1086/175437

CrossRef Full Text | Google Scholar

Saitoh, T. R., and Makino, J. (2013). A density-independent formulation of smoothed particle hydrodynamics. ApJ 768, 44. doi:10.1088/0004-637X/768/1/44

CrossRef Full Text | Google Scholar

Sano, T., and Stone, J. M. (2002). The effect of the Hall term on the nonlinear evolution of the magnetorotational instability. I. Local axisymmetric simulations. ApJ 570, 314–328. doi:10.1086/339504

CrossRef Full Text | Google Scholar

Sod, G. A. (1978). A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J. Comput. Phys. 27, 1–31. doi:10.1016/0021-9991(78)90023-2

CrossRef Full Text | Google Scholar

Springel, V. (2005). The cosmological simulation code GADGET-2. MNRAS 364, 1105–1134. doi:10.1111/j.1365-2966.2005.09655.x

CrossRef Full Text | Google Scholar

Springel, V. (2010). Smoothed particle hydrodynamics in astrophysics. ARA&A 48, 391–430. doi:10.1146/annurev-astro-081309-130914

CrossRef Full Text | Google Scholar

Springel, V., and Hernquist, L. (2002). Cosmological smoothed particle hydrodynamics simulations: the entropy equation. MNRAS 333, 649–664. doi:10.1046/j.1365-8711.2002.05445.x

CrossRef Full Text | Google Scholar

Stasyszyn, F. A., Dolag, K., and Beck, A. M. (2013). A divergence-cleaning scheme for cosmological SPMHD simulations. MNRAS 428, 13–27. doi:10.1093/mnras/sts018

CrossRef Full Text | Google Scholar

Stasyszyn, F. A., and Elstner, D. (2015). A vector potential implementation for smoothed particle magnetohydrodynamics. J. Comput. Phys. 282, 148–156. doi:10.1016/j.jcp.2014.11.011

CrossRef Full Text | Google Scholar

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., and Simon, J. B. (2008). Athena: a new code for astrophysical mhd. ApJS 178, 137–177. doi:10.1086/588755

CrossRef Full Text | Google Scholar

Thompson, S., and Cassibry, J. (2020). 2D modeling of fusion ignition conditions for a multilayer plasma liner magneto-inertial fusion target in a cylindrical configuration. Phys. Plasmas 27, 022701. doi:10.1063/1.5132336

CrossRef Full Text | Google Scholar

Tóth, G. (2000). The ∇·B=0 constraint in shock-capturing magnetohydrodynamics codes. J. Comput. Phys. 161, 605–652. doi:10.1006/jcph.2000.6519

CrossRef Full Text | Google Scholar

Tricco, T. S. (2019). The Kelvin-Helmholtz instability and smoothed particle hydrodynamics. MNRAS 488, 5210–5224. doi:10.1093/mnras/stz2042

CrossRef Full Text | Google Scholar

Tricco, T. S., and Price, D. J. (2012). Constrained hyperbolic divergence cleaning for smoothed particle magnetohydrodynamics. J. Comput. Phys. 231, 7214–7236. doi:10.1016/j.jcp.2012.06.039

CrossRef Full Text | Google Scholar

Tricco, T. S., and Price, D. J. (2013). A switch to reduce resistivity in smoothed particle magnetohydrodynamics. MNRAS 436, 2810–2817. doi:10.1093/mnras/stt1776

CrossRef Full Text | Google Scholar

Tricco, T. S., and Price, D. J. (2023). An integral-based approach for the vector potential in smoothed particle magnetohydrodynamics. arXiv e-prints , arXiv:2306.15039. doi:10.48550/arXiv.2306.15039

CrossRef Full Text | Google Scholar

Tricco, T. S., Price, D. J., and Bate, M. R. (2016a). Constrained hyperbolic divergence cleaning in smoothed particle magnetohydrodynamics with variable cleaning speeds. J. Comput. Phys. 322, 326–344. doi:10.1016/j.jcp.2016.06.053

CrossRef Full Text | Google Scholar

Tricco, T. S., Price, D. J., and Federrath, C. (2016b). A comparison between grid and particle methods on the small-scale dynamo in magnetized supersonic turbulence. MNRAS 461, 1260–1275. doi:10.1093/mnras/stw1280

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Iwasaki, K., and Inutsuka, S.-i. (2013). An explicit scheme for ohmic dissipation with smoothed particle magnetohydrodynamics. MNRAS 434, 2593–2599. doi:10.1093/mnras/stt1205

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., and Inutsuka, S. (2015a). Bimodality of circumstellar disk evolution induced by the Hall current. ApJ 810, L26. doi:10.1088/2041-8205/810/2/L26

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., and Inutsuka, S. (2015b). Effects of Ohmic and ambipolar diffusion on formation and evolution of first cores, protostars, and circumstellar discs. MNRAS 452, 278–288. doi:10.1093/mnras/stv1290

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Machida, M. N., and Inutsuka, S. (2021). Conditions for justifying single-fluid approximation for charged and neutral dust fluids and a smoothed particle magnetohydrodynamics method for dust-gas mixture. ApJ 913, 148. doi:10.3847/1538-4357/abf5db

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Machida, M. N., and Inutsuka, S.-i. (2023). Co-evolution of dust grains and protoplanetary disks. PASJ 75, 835–852. doi:10.1093/pasj/psad040

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Machida, M. N., Susa, H., Nomura, H., and Inutsuka, S. (2020). Early evolution of disk, outflow, and magnetic field of young stellar objects: impact of dust model. ApJ 896, 158. doi:10.3847/1538-4357/ab93d0

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Maury, A., Commerçon, B., Alves, F. O., Cox, E. G., Sakai, N., et al. (2022). The role of magnetic fields in the formation of protostars, disks, and outflows. arXiv e-prints , arXiv:2209.13765. doi:10.48550/arXiv.2209.13765

CrossRef Full Text | Google Scholar

Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., and Inutsuka, S. (2018). Does misalignment between magnetic field and angular momentum enhance or suppress circumstellar disk formation? ApJ 868, 22. doi:10.3847/1538-4357/aae4dc

CrossRef Full Text | Google Scholar

Tu, X., Wang, Q., Zheng, H., and Gao, L. (2022). Meshless methods for magnetohydrodynamics with vector potential. J. Comput. Phys. 470, 111596. doi:10.1016/j.jcp.2022.111596

CrossRef Full Text | Google Scholar

Vela Vela, L., Sánchez, R., Reynolds-Barredo, J. M., and Geiger, J. (2019). Magneto-hydrodynamical nonlinear simulations of magnetically confined plasmas using smooth particle hydrodynamics (SPH). Phys. Plasmas 26, 012511. doi:10.1063/1.5078451

CrossRef Full Text | Google Scholar

Wadsley, J. W., Keller, B. W., and Quinn, T. R. (2017). Gasoline2: a modern smoothed particle hydrodynamics code. MNRAS 471, 2357–2369. doi:10.1093/mnras/stx1643

CrossRef Full Text | Google Scholar

Wardle, M., and Ng, C. (1999). The conductivity of dense molecular gas. MNRAS 303, 239–246. doi:10.1046/j.1365-8711.1999.02211.x

CrossRef Full Text | Google Scholar

Wissing, R., and Shen, S. (2020). Smoothed particle magnetohydrodynamics with the geometric density average force expression. A&A 638, A140. doi:10.1051/0004-6361/201936739

CrossRef Full Text | Google Scholar

Wissing, R., and Shen, S. (2023). Numerical dependencies of the galactic dynamo in isolated galaxies with SPH. A&A 673, A47. doi:10.1051/0004-6361/202244753

CrossRef Full Text | Google Scholar

Wissing, R., Shen, S., Wadsley, J., and Quinn, T. (2022). Magnetorotational instability with smoothed particle hydrodynamics. A&A 659, A91. doi:10.1051/0004-6361/202141206

CrossRef Full Text | Google Scholar

Wurster, J., Bate, M. R., and Bonnell, I. A. (2021). The impact of non-ideal magnetohydrodynamic processes on discs, outflows, counter-rotation, and magnetic walls during the early stages of star formation. MNRAS 507, 2354–2372. doi:10.1093/mnras/stab2296

CrossRef Full Text | Google Scholar

Wurster, J., Bate, M. R., and Price, D. J. (2018a). The collapse of a molecular cloud core to stellar densities using radiation non-ideal magnetohydrodynamics. MNRAS 475, 1859–1880. doi:10.1093/mnras/stx3339

CrossRef Full Text | Google Scholar

Wurster, J., Bate, M. R., and Price, D. J. (2018b). The effect of extreme ionization rates during the initial collapse of a molecular cloud core. MNRAS 476, 2063–2074. doi:10.1093/mnras/sty392

CrossRef Full Text | Google Scholar

Wurster, J., Bate, M. R., and Price, D. J. (2019). There is no magnetic braking catastrophe: low-mass star cluster and protostellar disc formation with non-ideal magnetohydrodynamics. MNRAS 489, 1719–1741. doi:10.1093/mnras/stz2215

CrossRef Full Text | Google Scholar

Wurster, J., Bate, M. R., Price, D. J., and Bonnell, I. A. (2022). On the origin of magnetic fields in stars - II. The effect of numerical resolution. MNRAS 511, 746–764. doi:10.1093/mnras/stac123

CrossRef Full Text | Google Scholar

Wurster, J., Price, D., and Ayliffe, B. (2014). Ambipolar diffusion in smoothed particle magnetohydrodynamics. MNRAS 444, 1104–1112. doi:10.1093/mnras/stu1524

CrossRef Full Text | Google Scholar

Wurster, J., Price, D. J., and Bate, M. R. (2016). Can non-ideal magnetohydrodynamics solve the magnetic braking catastrophe? MNRAS 457, 1037–1061. doi:10.1093/mnras/stw013

CrossRef Full Text | Google Scholar

Wurster, J., Price, D. J., and Bate, M. R. (2017). The impact of non-ideal magnetohydrodynamics on binary star formation. MNRAS 466, 1788–1804. doi:10.1093/mnras/stw3181

CrossRef Full Text | Google Scholar