© ESO 2020

1. Introduction

Magnetic fields are important in a wide array of different astrophysical systems. In star formation, they govern the dynamics at several stages during collapse. They are critical in the launching of jets from a broad range of sources. They also play a major role in the transport of angular momentum in ionized accretion disks due to the magnetorotational instability. Magnetic fields have been largely neglected in galaxy formation simulations, mostly due to the technical difficulties associated with them. It is only recently that researchers have begun to apply them (Wang & Abel 2009; Kotarba et al. 2011; Pakmor & Springel 2013; Rieder & Teyssier 2016; Butsky et al. 2017; Su et al. 2017; Pakmor et al. 2017; Steinwandel et al. 2019). The importance of magnetic fields in galaxy formation is clear from observations of the Milky Way and nearby galaxies, which reveal that the magnetic energy is in equipartition with the thermal and turbulent energies (Boulares & Cox 1990; Beck et al. 1996). This means that they are likely to have a large dynamical effect on the evolution of the galaxy, adding significant non-thermal pressure that can suppress star-formation (Pakmor & Springel 2013). In addition, it has been shown that magnetic fields have a strong impact on fluid instabilities (Jun et al. 1995; McCourt et al. 2015), which may affect how gas in the intergalactic medium (IGM) accretes onto galaxies and how gas in galactic outflows leaves (or cycles back to) galaxies. The strength and structure of magnetic fields in galaxies also determine the transport of cosmic rays (CRs), which has recently emerged as a promising candidate for driving galactic outflows because they have long cooling time scales (Uhlig et al. 2012; Booth et al. 2013; Pakmor et al. 2016; Butsky & Quinn 2018).

Apart from the improvements in the general method, advances have been made in the magnetohydrodynamics extension for smoothed particle hydrodynamics (SPMHD). The modern foundation of SPMHD comes largely from the work of Price & Monaghan (2004) which was built on the earlier work of Phillips & Monaghan (1985). The two main technical difficulties to overcome for SPMHD, are the handling of divergence errors and the choice of an artificial resistivity term to capture shocks and discontinuities in the magnetic field.

Artificial dissipation terms are required to smooth out discontinuities in any fluid quantity in all numerical hydrodynamics methods. In SPH, this is most commonly achieved via explicit artificial dissipation. To avoid excessive dissipation away from shocks and discontinuities, switches have been developed to limit where the artificial dissipation terms are active. For magnetic fields, newly developed artificial resistivity switches (Price et al. 2018; Tricco & Price 2013) have significantly reduced the amount of dissipation and improved the method in the weak field regime.

Unphysical divergence errors (magnetic monopoles) can arise from the discretization and numerical integration of the MHD equations. Divergence errors in SPMHD for magnetic-dominated scenarios need to be handled with care, as they can produce a negative force between particles which leads to the tensile instability (Monaghan 2000). As such, the force produced from the divergence needs to be partly removed in the strong field regime for the method to remain stable, this breaks momentum and energy conservation in proportion to the divergence error. It is therefore crucial to try to keep the divergence error as close to zero as possible. While grid codes have access to the constrained transport scheme (Evans & Hawley 1988), which ensures a divergence free field up to machine precision, it cannot easily be implemented within meshless methods, due to the absence of regular spatial grid surfaces. Generation of divergence free fields in SPMHD have been explored in detail, however, all of them suffer from problems. Generation of magnetic fields from Euler potentials (B = ∇α × ∇β) cannot wind the magnetic field and thereby not produce a dynamo (Brandenburg 2010). Price (2010) showed that vector potential implementations (B = ∇ × A) are plagued with numerical instabilities. However, Stasyszyn & Elstner (2015) recently showed that with additional diffusion, smoothing of the magnetic field, and enforcing the Coulomb gauge (∇ ⋅ A = 0), the vector potential formalism could remain stable for a handful of test cases. Additional testing would be required to determine the robustness of the method. The most popular method to deal with divergence error in meshless methods, is to evolve the magnetic field via the induction equation and then to “clean” the divergence away. In general, this is done by introducing a separate scalar field which couples to the induction equation such that it produces a damped wave equation for the divergence error, so the divergence is spread outward like a damped wave. The method was first developed in Dedner et al. (2002) and was improved by Tricco & Price (2012), who introduced a constrained version of the method. This ascertains that the magnetic energy is either conserved or dissipated. This was updated in Tricco et al. (2016) to correctly allow variable cleaning speed, which further improved the method.

These new improvements in artificial dissipation and divergence error controlling have significantly increased the accuracy and convergence of the SPMHD method. There have also been implementations of non-ideal MHD in SPMHD proposed recently (Tsukamoto et al. 2013, 2015a,b; Wurster et al. 2014, 2016; Price et al. 2018), which include Ohmic resistivity, ambipolar diffusion, and the Hall effect.

As mentioned previously, the numerical surface tension seen in traditional SPH (TSPH) can be solved by using a different gradient operator (GDSPH; Wadsley et al. 2017). This substantially improves the accuracy of pressure forces across density jumps and provides a more physical form for the internal energy equation, where it represents a direct discretization of $\frac{\mathrm{d}u}{\mathrm{d}t}=-\frac{P}{\rho }\mathrm{\nabla }·v$ from the Euler equations while retaining all the usual conservation properties. In Wadsley et al. (2017), the authors show that GDSPH, together with an explicit turbulent diffusion term on thermal energy, yields excellent results in fluid mixing test cases, such as the Kelvin-Helmholtz instability and the blob test.

In this paper, we investigate how SPMHD benefits from the use of GDSPH. As such, we have implemented MHD within the GASOLINE2 (Wadsley et al. 2017) and CHANGA (Menon et al. 2015) codes, which both utilize the GDSPH formalism. GASOLINE2 is a highly parallel, state-of-the-art code for cosmological structure formation simulations which includes all the features of modern SPH methods. CHANGA includes all the same SPH methods as GASOLINE2, but it is written in an inherently parallel language CHARM++ (Kale & Krishnan 1993) which enables more efficient parallelization. The major difference between the two codes lies in the gravity solver, which is different in CHANGA because it uses an oct-tree, rather than an arbitrary binary KD-tree as in GASOLINE2.

This paper is organized as follows. In Sect. 2, we go through the SPMHD theory and show how the equations can be formulated using the GDSPH approach. In Sect. 3, we test our implementation on a large suite of standard test cases and in Sect. 4 we apply the code to an astrophysical application: the collapse of a magnetized cloud. In Sect. 5, we discuss our results and present some concluding remarks.

2. Theory

In this section, we show how the MHD equations can be formulated in a conservative way within the GDSPH framework. The development is similar to the findings in previous work (Price & Monaghan 2004; Price 2012; Tricco et al. 2016; Price et al. 2018) and we direct the reader to these papers for additional background details.

2.1. MHD theory

The two main equations which are relevant for ideal MHD are the Lorentz force law and the induction equation. Assuming that the fluid is an ideal conductor (E = 0), the Lorentz force law can be written as:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=\frac{1}{{\mu }_0\rho }(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}=\frac{1}{{\mu }_0\rho }(-\frac{1}{2}\mathrm{\nabla }{\mathbf{B}}^2+(\mathbf{B}·\mathrm{\nabla })\mathbf{B}),\end{array}$$(1)

where v, ρ, B and μ₀ is the velocity, density, magnetic field and vacuum permeability, respectively. The first term acts like an isotropic magnetic pressure term, while the other term acts as an attractive term along magnetic field lines (tension). Going forward, we define code units such that μ₀ = 1. The conservative form of SPMHD is attained by using the stress tensor to describe the momentum equation. Assuming that the magnetic field is divergence free, the MHD stress tensor can be written as:

$$\begin{array}{c}\hfill S^{\mathit{ij}}=-{\delta }^{\mathit{ij}}(P+\frac{{\mathbf{B}}^2}{2})+B^i{B}^j,\end{array}$$(2)

where P is the thermal pressure and δ^ij is the Kronecker delta. The momentum equation can then be written as:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=\frac{\mathrm{\nabla }·\mathbf{S}}{\rho }=-\frac{1}{\rho }\mathrm{\nabla }(P+\frac{{\mathbf{B}}^2}{2})+\frac{1}{\rho }[(\mathbf{B}·\mathrm{\nabla })\mathbf{B}+\mathbf{B}(\mathrm{\nabla }·\mathbf{B})].\end{array}$$(3)

There is an extra tension force term which would normally have no effect due to the assumption ∇ ⋅ B = 0. However, as mentioned in the introduction, this constraint is usually not fully upheld in numerical codes. To avoid numerical instability within SPH, this term needs to be negated when the magnetic pressure exceeds the thermal pressure.

The change in the magnetic field is obtained from the induction equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}=\mathrm{\nabla }\times (\mathbf{v}\times \mathbf{B})=(\mathbf{B}·\mathrm{\nabla })\mathbf{v}-\mathbf{B}(\mathrm{\nabla }·\mathbf{v}),\end{array}$$(4)

where the first term affects the magnetic field through shearing motion, while the second will increase the magnetic field when undergoing compression. A combined effect of the two terms is to enhance the field due to compression perpendicular to the field direction (for example, B ∝ ρ^2/3 for spherical collapse). Compression in the direction of the field has no effect.

2.2. SPH discretization

Derivatives within SPH can be discretized in a number of ways, and a general formulation is given by Price (2012):

$$\begin{array}{cc}& \frac{\mathrm{\nabla }A}{\rho }=\frac{\varphi }{\rho }[\frac{A}{{\varphi }^2}\mathrm{\nabla }\varphi +\mathrm{\nabla }\left(\frac{A}{\varphi }\right)]\approx \sum _b\frac{m_b}{{\rho }_a{\rho }_b}(A_a\frac{{\varphi }_b}{{\varphi }_a}+A_b\frac{{\varphi }_a}{{\varphi }_b}){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}},\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{\mathrm{\nabla }A}{\rho }=\frac{1}{\varphi \rho }[\mathrm{\nabla }(\varphi A)-A\mathrm{\nabla }\varphi ]\approx \sum _b\frac{m_b}{{\rho }_a{\rho }_b}\frac{{\varphi }_b}{{\varphi }_a}(A_b-A_a){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}},\hfill \end{array}$$(6)

where ϕ can be any arbitrary, differentiable scalar quantity. The geometric density average force formulation (GDSPH) corresponds to using ϕ = 1 while traditional SPH corresponds to using ϕ = ρ. GDSPH therefore gives the following symmetric and anti-symmetric gradient operators:

$$\begin{array}{cc}& \frac{\mathrm{\nabla }A}{\rho }\approx \frac{1}{\rho }(A\mathrm{\nabla }1+\mathrm{\nabla }A)=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}(A_a+A_b){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \frac{\mathrm{\nabla }A}{\rho }\approx \frac{1}{\rho }(\mathrm{\nabla }A-A\mathrm{\nabla }1)=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}(A_b-A_a){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}.\hfill \end{array}$$(8)

Here, ${\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}$ is a symmetric gradient of the smoothing kernel:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}=\frac{1}{2}[f_a{\mathrm{\nabla }}_{a}W(r_{\mathit{ab}},h_a)+f_b{\mathrm{\nabla }}_{b}W(r_{\mathit{ab}},h_b)],\end{array}$$(9)

where W is the smoothing kernel, h_a is the smoothing length of particle a, and r_ab = |r_a − r_b| is the distance between particle a and b. Here, f_a is a correction term introduced in Wadsley et al. (2017) to ensure that internal energy and density evolve consistently, such that entropy is tightly conserved. To attain a conservative formalism for SPH, the symmetric gradient operator is applied to the equations of motion and the anti-symmetric gradient operator is applied to the internal energy equation¹. As a consequence, zeroth order errors arise in the equations of motions which will depend on the local particle distribution². A generalized error term for the zeroth order errors is given by:

$$\begin{array}{c}\hfill {\mathbf{E}}_0=\sum _b\frac{m_b}{{\rho }_b}({\mathrm{\Phi }}_{\mathit{ab}}+{\mathrm{\Phi }}_{\mathit{ab}}^{-1}){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}},\end{array}$$(10)

where $\mathrm{\Phi }=\frac{{\varphi }_a}{{\varphi }_b}$ depend on the chosen scalar quantity ϕ in Eqs. (5) and (6). As shown by Read et al. (2010), in TSPH ${\mathrm{\Phi }}_{\mathit{ab}}=\frac{{\rho }_a}{{\rho }_b}$, while in GDSPH Φ_ab = 1. It is then evident that these errors are more severe at density gradients in TSPH than in GDSPH (where they are explicitly independent of the density gradient). A similar improvement can be seen for the linear errors.

Applying the symmetric gradient operator to the momentum equation (Eq. (3)) and the anti-symmetric gradient operator to the induction equation (Eq. (4)) gives:

$$\begin{array}{cc}& \frac{\mathrm{d}{v}_a^i}{\mathrm{d}t}={\displaystyle \sum _b}\frac{m_b}{{\rho }_a{\rho }_b}(S_a^{\mathit{ij}}+S_b^{\mathit{ij}}){\mathrm{\nabla }}_a^j{\overline{W}}_{\mathit{ab}}+f_{\mathrm{div}B,a}^i,\hfill \end{array}$$(11)

$$\begin{array}{cc}& \frac{\mathrm{d}{\mathbf{B}}_a}{\mathrm{d}t}={\displaystyle \sum _b}\frac{m_b}{{\rho }_b}[{\mathbf{B}}_a({\mathbf{v}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}})-{\mathbf{v}}_{\mathit{ab}}({\mathbf{B}}_a·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}})],\hfill \end{array}$$(12)

where v_ab = v_a − v_b. The stability term $f_{\text{div}B,a}^i$ is added to avoid the tensile instability. This can occur due to divergence errors when the magnetic pressure exceeds the gas pressure ($\frac{B^2}{2}>P$) (Phillips & Monaghan 1985). The stability term is defined as:

$$\begin{array}{c}\hfill f_{\mathrm{div}B,a}^i=-\widehat{B_a^i}{\displaystyle \sum _b}\frac{m_b}{{\rho }_a{\rho }_b}({\mathbf{B}}_a+{\mathbf{B}}_b)·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}.\end{array}$$(13)

This basically removes the divergence term $(-\frac{\mathbf{B}}{\rho }\mathrm{\nabla }·\mathbf{B})$ from Eq. (3) (Børve et al. 2001; Price 2012). Removing a term from the conservative momentum equation effectively breaks momentum conservation. However, the error introduced will be proportional to the divergence. To minimize its effect in the weak field regime, we use the scheme from Børve et al. (2004) with a factor of $\widehat{B_a^i}=B_a^i$ for β <  1 as advocated by Tricco & Price (2012):

$$\begin{array}{c}\hfill \widehat{B_a^i}=\{\begin{array}{cc}{B}_a^i\hfill & \beta <1\hfill \\ B_a^i(2-\beta )\hfill & 1<\beta <2\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\end{array}$$(14)

where $\beta =\frac{2P}{B^2}$ is the plasma beta.

2.3. Treating discontinuities

When fluid quantities become discontinuous, they are no longer differentiable, which is problematic as differentiability is assumed by the SPMHD equations. Artificial resistivity is required to smooth out discontinuities in the magnetic field which can occur both along and orthogonal to the fluid motion and in both compression and rarefaction. The artificial resistivity can be represented as an isotropic diffusion:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{B}}{{\mathrm{d}t}_{\mathrm{diss}}}=\eta {\mathrm{\nabla }}^2\mathbf{B},\end{array}$$(15)

where η is a resistivity parameter. We use the Brookshaw method (Brookshaw 1985), which estimates the second derivative by using the first derivative kernel and the difference in the field divided by the particle spacing. Following the GDSPH discretization, we get:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{B}}_a}{{\mathrm{d}t}_{\mathrm{diss}}}={\displaystyle \sum _b}\frac{m_b}{{\rho }_b}\left(\frac{{\eta }_a+{\eta }_b}{|{\mathbf{r}}_{\mathit{ab}}|}\right){\mathbf{B}}_{\mathit{ab}}({\widehat{\mathbf{r}}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}),\end{array}$$(16)

where B_ab = B_a − B_b and ${\widehat{\mathbf{r}}}_{\mathit{ab}}={\mathbf{r}}_{\mathit{ab}}/|{\mathbf{r}}_{\mathit{ab}}|$. To conserve energy, the change in the internal energy becomes:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{u}_a}{{\mathrm{d}t}_{\mathrm{diss}}}=-\frac{1}{2}{\displaystyle \sum _b}\frac{m_b}{{\rho }_a{\rho }_b}\left(\frac{{\eta }_a+{\eta }_b}{|{\mathbf{r}}_{\mathit{ab}}|}\right){\mathbf{B}}_{\mathit{ab}}^2({\widehat{\mathbf{r}}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}).\end{array}$$(17)

To reduce dissipation away from shocks, we can introduce a varying and resolution dependent resistivity parameter:

$$\begin{array}{c}\hfill \eta =\frac{1}{2}{\alpha }_B{v}_{\mathrm{sig},B}|{\mathbf{r}}_{\mathit{ab}}|,\end{array}$$(18)

where α_B is a dimensionless coefficient and v_sig, B is the signal speed. Proper choice of α_B and v_sig, B makes the artificial resistivity second order accurate away from shocks (η ∝ h²). We choose to implement the resistivity from PHANTOM (Price et al. 2018) where the signal speed is activated following:

$$\begin{array}{c}\hfill v_{\mathrm{sig},B}=|{\mathbf{v}}_{\mathit{ab}}\times {\widehat{\mathbf{r}}}_{\mathit{ab}}|.\end{array}$$(19)

The dimensionless coefficient α_B is set to a constant. In PHANTOM, this coefficient is set to α_B = 1, however, from our tests we find that α_B = 0.5 provides sufficient dissipation. This switch was shown to be the least dissipative compared to previous switches, while still capturing the correct magnetic features (Wurster et al. 2017).

2.4. Divergence cleaning

As we discussed in the introduction, divergence errors are generated by the discretization and integration of the MHD equations. Apart from creating an unphysical magnetic field, it also forces us to introduce a stability term (Eq. (13)), which breaks momentum conservation in the strong field regime. This makes it crucial to reduce the divergence errors as much as possible. The best way found in SPMHD is by introducing a divergence cleaning scheme (Tricco & Price 2012). In general, this is done by introducing a separate scalar field which couples to the induction equation, such that it produces a damped wave equation for the divergence error. That is, the divergence is spread outward like a damped wave. In our implementation, we employ the constrained hyperbolic divergence cleaning from Tricco et al. (2016), an improved version of the method presented by Dedner et al. (2002). The constrained hyperbolic divergence cleaning ensures that magnetic energy is either conserved or dissipated. In this method, a scalar field ψ is coupled to the induction equation as follows:

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}\right)}_{\psi }=-\mathrm{\nabla }\psi .\end{array}$$(20)

The scalar field ψ evolves according to:

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\psi }{c_h}\right)=-c_h\mathrm{\nabla }·\mathbf{B}-\frac{1}{\tau }\frac{\psi }{c_h}-\frac{1}{2}\psi (\mathrm{\nabla }·\mathbf{v}).\end{array}$$(21)

where τ is the decay time and c_h is the wave cleaning speed:

$$\begin{array}{cc}& c_h=f_{\mathrm{clean}}{v}_{\mathrm{mhd}},\hfill \end{array}$$(22)

$$\begin{array}{cc}& v_{\mathrm{mhd}}=\sqrt{c_s^2+v_{\mathrm{A}}^2},\hfill \end{array}$$(23)

$$\begin{array}{cc}& v_{\mathrm{A}}=\sqrt{\frac{B^2}{\rho }}.\hfill \end{array}$$(24)

Here, c_s is the speed of sound, v_A the Alfvén velocity, and f_clean is an overcleaning factor. The f_clean factor can be used to increase the amount of divergence cleaning, however, this will reduce the timestep³ according to Δt → Δt/f_clean. Combining the cleaning equation with the induction equation produces a damped wave equation for the divergence (this form assumes constant c_h and τ):

$$\begin{array}{c}\hfill \frac{{\partial }^2(\mathrm{\nabla }·\mathbf{B})}{\partial t^2}-c_h^2{\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{B})+\frac{1}{\tau }\frac{\partial (\mathrm{\nabla }·\mathbf{B})}{\partial t}=0,\end{array}$$(25)

which effectively shows that the divergence is spread out and damped. The decay time is given by:

$$\begin{array}{c}\hfill {\tau }_a=\frac{h_a}{c_{h,a}{\sigma }_c}.\end{array}$$(26)

Here, σ_c is a dimensional constant, and was shown to be optimal with a value of 1 in 3D. Following Tricco & Price (2012), ∇ψ is discretized using the symmetric gradient operator (Eq. (7)) and ∇ ⋅ B using the anti-symmetric gradient operator (Eq. (8)). Within the GDSPH discretization, Eqs. (20) and (21) become:

$$\begin{array}{cc}& {\left(\frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}\right)}_{\psi ,a}=-{\displaystyle \sum _b}\frac{m_b}{{\rho }_b}({\psi }_a+{\psi }_b){\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}},\hfill \end{array}$$(27)

$$\begin{array}{cc}& \frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\psi }{c_h}\right)}_a=c_h^a{\displaystyle \sum _b}\frac{m_b}{{\rho }_b}{\mathbf{B}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}+\frac{{\psi }_a}{2c_h^a}{\displaystyle \sum _b}\frac{m_b}{{\rho }_b}{\mathbf{v}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}}-\frac{{\psi }_a}{c_h^a{\tau }_a}.\hfill \end{array}$$(28)

The divergence cleaning dissipates energy from the magnetic field. However, this term is so small compared to the other dissipation terms that it is not worth accounting for. We could, of course, add this energy to heat and conserve energy, however, as discussed by Tricco & Price (2012), the removal of magnetic energy and subsequent generation of thermal energy would be non-local due to the coupling of parabolic diffusion with hyperbolic transport. Due to this, we simply removed the energy.

To ensure that simulations are not affected by the divergence error, we monitor the normalized divergence error:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{div}B}=\frac{h|\mathrm{\nabla }·\mathbf{B}|}{|B|}.\end{array}$$(29)

The mean of this quantity should preferably remain below 10⁻². However, regions of locally high divergence error can occur, so careful inspection of the divergence error is required to ensure the quality of simulations.

2.5. Shock capturing

To correctly capture shocks in MHD, we need to modify the artificial viscosity term in the momentum equation (see Wadsley et al. 2017 for a detailed description of the artificial viscosity term in GASOLINE2). For MHD the sound speed is replaced by the fast magnetosonic speed (Eq. (23)). We also modify the gradient-based shock detector introduced in GASOLINE2, which determines the direction of the shock from the pressure gradient. For the MHD, we must include the Lorentz force to correctly determine the direction of the shock. A more general way to determine the direction of the shock is to estimate the acceleration of the MHD forces without the dissipation terms before the actual force calculation:

$$\begin{array}{c}\hfill \widehat{\mathbf{n}}=-{\left(\frac{\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}}{\left|\frac{\mathrm{d}v}{\mathrm{d}t}\right|}\right)}_{\mathrm{no}\mathrm{diss}}.\end{array}$$(30)

This addition improves the behavior of shock detection in convergent flows for MHD. In GASOLINE2, the diffusion of fluid scalar variables such as thermal energy, metals and so forth are modeled using subgrid turbulent mixing (Wadsley et al. 2008; Shen et al. 2010). However, we found that in strong shocks like the MHD blastwave, the thermal dissipation is not enough and can lead to incorrect velocity profiles. As such, we add a thermal shock dissipation similar to Eqs. (4.5) with (4.8) in Monaghan (1992), however, we use ${\overline{c}}_{\mathit{ab}}=0$ and with a larger constant (fitted parameter from the GASOLINE2 code):

$$\begin{array}{cc}& \frac{\mathrm{d}u}{\mathrm{d}{t}_{\mathrm{shock}}}=-{\displaystyle \sum _b}\frac{m_b}{{\overline{\rho }}_{\mathit{ab}}}{\mathrm{d}}_{\mathrm{shock}}{u}_{\mathit{ab}}\frac{({\widehat{\mathbf{r}}}_{\mathit{ab}}·{\mathrm{\nabla }}_a{\overline{W}}_{\mathit{ab}})}{\left|r_{\mathit{ab}}\right|},\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\mathrm{d}}_{\mathrm{shock}}=16{\overline{h}}_{\mathit{ab}}|{\mu }_{\mathit{ab}}|,\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\mu }_{\mathit{ab}}=\{\begin{array}{cc}\frac{{\overline{h}}_{\mathit{ab}}({\mathbf{v}}_{\mathit{ab}}·{\mathbf{r}}_{\mathit{ab}})}{r_{\mathit{ab}}^2+0.01{\overline{h}}_{\mathit{ab}}^2}\hfill & \text{for}{\mathbf{v}}_{\mathit{ab}}·{\mathbf{r}}_{\mathit{ab}}<0,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\hfill \end{array}$$(33)

3. Test problems

In this section, we present the results from our test cases. All the simulations were run with the MHD version of GASOLINE2. To remain consistent and show the production quality of the method, we decided to run all the tests in 3D and with a default set of code parameters (described below). While glass-like initial conditions should always be used to correctly capture the natural state of 3D SPH simulations, for the sake of comparison, we elected to follow the initial setups from other authors, which often use lattice-based initial conditions. Test cases that are originally 1D or 2D are made 3D by extending the non-active dimensions by a set number of particles. By default, GASOLINE2 sets the smoothing length based on a fixed number of neighbours. However, we found that in very uniform and precise tests, as in the circularized Alfvén wave test, this approach generates small force errors that generate perturbations in the traveling wave. As such, we made the smoothing length directly proportional to the density and simultaneously determined the density and smoothing length using an iterative summation (Springel & Hernquist 2002). We note, however, that in all other tests, no visible effect was seen. We ran simulations with both TSPH and GDSPH, and the only difference between them is the choice of ϕ in Eqs. (5) and (6); all the other numerical schemes and parameters remain the same. In many of these tests, we compare the results to the state-of-the-art SPMHD code PHANTOM (Price et al. 2018), and to the PSPH and the new meshless finite mass/volume (MFM/MFV) method of the GIZMO code (Hopkins 2015). The MFM/MFV method utilizes a Lagrangian Godunov type method that employs more complex gradient operators and calculates fluxes from Riemann solvers.

Default set of code parameters. For the smoothing kernel, we used a Wendland C4 kernel (Wendland 1995) with 200 neighbours⁴. Our artificial viscosity (AV) followed the prescription given in Wadsley et al. (2017), the AV parameters were set to α_max = 4, α_min = 0, $\tau =0.1\frac{h}{c}$ and β = 2. The artificial resistivity (AR) followed from the method outlined in Sect. 2.2, and the AR parameter was set to α_B = 0.5. The thermal diffusion followed the turbulent mixing model described in Wadsley et al. (2008) and Shen et al. (2010), with the turbulent diffusion coefficient set to C = 0.03.

3.1. Circularized polarized Alfvén wave

The circularized polarized Alfvén wave test was introduced by Tóth (2000) to serve as an analytical solution to the ideal MHD equations. Due to the waves being circularized, the gradient in magnetic pressure is zero and the wave should remain the same after each period. This proves to be a useful test for gauging the dissipation and dispersion of the MHD implementation. This test is sensitive to the tensile instability (Price & Monaghan 2005), so it also serves as a good test to see if the stability term (Eq. (13)) properly stabilizes the solution. The setup follows Gardiner & Stone (2008) and Price et al. (2018), in which the waves are traveling at an angle of θ = 30° with respect to the x axis, within a periodic box of length L = (l, l/2, l/2), where l = 3. The transverse velocities and magnetic fields are circularized:

$$\begin{array}{cc}& B_{\perp ,1}=v_{\perp ,1}=0.1sin(2\pi x_{\Vert }/\lambda ),\hfill \\ \hfill & B_{\perp ,2}=v_{\perp ,2}=0.1cos(2\pi x_{\Vert }/\lambda ),\hfill \end{array}$$

while the parallel components are set to:

$$\begin{array}{c}\hfill v_{\Vert }=0B_{\Vert }=1.\end{array}$$

Here, x_∥ is the direction of propagation, and λ = 1 is the wavelength. An adiabatic EOS (γ = 5/3) is used with uniform pressure P = 0.1 and density ρ = 1.0. The particles are set up on a close-packed lattice and the simulation is run for five periods (t = 5). As we have uniform density, there are no differences between GDSPH and TSPH in this test case. We plot the transverse component of the magnetic field in the direction of propagation, showing the results of different resolutions [n_x, n_y, n_z]=[128, 74, 78], [64, 36, 39], and [32, 18, 18] in the upper panel of Fig. 1. From the results, we can see that both the phase and amplitude converge towards the analytical solution as we increase resolution. For a more qualitative look at the convergence, we perform a convergence study for this test using different smoothing kernels and neighbour numbers. In the lower panel of Fig. 1, we show the L1 error norm for the transverse magnetic field at five different resolutions (n_x = 32, 48, 64, 96, 128), and as we can see all the kernels exhibits second-order convergence. The major outliner is at low resolution for the Wendland C4 kernel with more neighbours. This is caused by the larger smoothing length, which at low resolution becomes comparable to half the wave length of the Alfvén wave. This makes the MHD gradients more ill defined which causes force errors that shows itself predominately as a phase shift in the Alfvén wave as time goes by. The amplitude of the wave is only weakly affected by this. From the bottom panel in Fig. 1 we can also see that the Wendland kernel has slightly lower convergence speed then the quintic kernel at higher resolution. Despite this result, we chose to go with the Wendland kernel C4 with 200 neighbors as our code default for the forthcoming tests. This is for several reasons. First, the quintic kernel is susceptible to the pairing instability whereas the Wendland kernels are not (Dehnen & Aly 2012). In addition, the Wendland kernels tend to make the particle distribution remain well ordered in dynamical conditions, which improves the overall accuracy of the method (Rosswog 2015). While the computational cost increases roughly linearly with increased neighbour number, gravity is usually the more dominant cost in astrophysical simulations, which means that the increase in cost is usually not significant. In the end, the choice of kernel and neighbour number will depend on the application at hand. However, for simulations involving subsonic flows, a high neighbour number has been shown to be preferred (as showcased by the Gresho-Chan vortex test in Dehnen & Aly 2012 and Rosswog 2015).

Fig. 1. Results of the 3D circularly polarized Alfvén wave test. Top panel: transverse component of the magnetic field in the direction of propagation after five periods. The analytical/initial solution is plotted in black, and the simulation results with resolution [nx, ny, nz]=[128, 74, 78] in red, [nx, ny, nz]=[64, 36, 39] in green, and [nx, ny, nz]=[32, 18, 18] in blue. Both of these are with Wendland C4 kernel with 200 neighbours. Bottom panel: convergence study for the Alfvén wave test using different kernels and neighbour numbers. Shows how the L1 error scales with resolution (particles along the x-axis). The code default (Wendland C4 kernel with 200 neighbours) is shown in green, Wendland C4 kernel with 114 neighbours are shown in magenta, Quintic kernel with 114 neighbours are shown in blue and the dashed brown line shows the curve for second order convergence. Convergence towards the analytical solution for all kernels are close to second order. When the smoothing length becomes comparable to half the wave length of the Alfvén wave, the MHD gradients becomes more ill defined which causes the slower convergence speed for the Wendland kernel (${\overline{N}}_{\mathrm{neigh}}=200$) at low resolution.

3.2. Advection loop

The advecting current loop test was introduced by Gardiner & Stone (2005, 2008), in which a weak magnetic loop is advected by a constant velocity field. As the ratio between the thermal pressure and magnetic pressure is massive (β ≈ 10⁶), the magnetic field is dynamically unimportant and should simply be advected along the velocity field. This proves to be one of the more difficult tests for grid-based code due to intrinsic dissipation during advection. We followed the setup from Gardiner & Stone (2008), Hopkins & Raives (2016) and Price et al. (2018), and initialized a 3D thin periodic box with length $L=(2,1,2\frac{\sqrt{6}}{n_x})$, velocity $v=(2,1,\frac{0.1}{\sqrt{5}})$ and pressure P = 1. The magnetic field inside the loop was determined from the potential A_z = A₀(R₀ − r), where A₀ = 10⁻³, R₀ = 0.3, and r² = x² + y². The face-centered magnetic fields are then ${\mathit{B}}_{\mathbf{0}}=\mathrm{\nabla }\times A_z=\frac{A_0}{r}(y,-x,0)$ inside the loop and zero everywhere else. We set up two initial conditions, one with a uniform density ρ = 1 with resolution [n_x, n_y, n_z]=[128, 74, 12] and another with a density gradient ($\mathrm{\Delta }\equiv \frac{{\rho }_{\mathrm{in}}}{{\rho }_{\mathrm{out}}}=2$) between the inner loop (ρ_in = 2) and outer medium (ρ_out = 1), with resolution [n_x, n_y, n_z]=[256, 148, 12]. The particles are set up on a close-packed lattice and the loop is advected for twenty periods with all the default dissipation and divergence cleaning terms turned on. The results of the test can be seen in Figs 2 and 3.

Fig. 2. Results from the advection loop test in 3D, showing a rendering of |B| in units of the initial value |B0|. Top left panel: initial setup (same rendering for nx = 128 and nx = 256, and top right: uniform density case for GDSPH with resolution nx = 128 after twenty crossings t = 20. We can see that the current loop is conserved almost perfectly even with all the dissipation terms turned on. This shows an improvement comparing to the SPMHD results in PHANTOM, where dissipation is seen after five periods (Fig. 35 in Price et al. 2018), which plots the current density. The bottom two panels show the cases with a density gradient $\mathrm{\Delta }\equiv \frac{{\rho }_{\mathrm{in}}}{{\rho }_{\mathrm{out}}}=2$ between the inner loop and outer medium and with resolution of nx = 256, after twenty crossings t = 20. The bottom left panel shows the GDSPH case and the bottom right show the TSPH case (nx = 256). Both show similar dissipation of the magnetic field to the results from the MFM/MFV method in GIZMO (Hopkins & Raives 2016).

Fig. 3. Results from the advection loop test in 3D, showing the time evolution of the magnetic field energy in units of the initial value. After t = 20 our uniform case has dissipated about 0.3%, while the cases with the density gradient Δ = 2 have dropped around 5%, which is similar to the dissipation displayed by the MFM/MFV method in Hopkins & Raives (2016). There is a tiny difference between GDSPH and TSPH, but this owes itself to differences in the initial reordering.

From the figures, we can see that in the case with uniform density the field loop is closely conserved, resulting in only 0.3% reduction in magnetic energy after twenty periods. This is a significant improvement from the advection loop presented in Price et al. (2018), which starts to degrade after five periods. We have also tested this case using a quintic spline kernel, a smaller number of neighbours and an α_B = 1 for the AR similar to Price et al. (2018), while a little more degradation can be seen, the difference in magnetic energy is still small after twenty periods (1% instead of a 0.3% decrease in magnetic energy). Without any dissipation and divergence cleaning the advection loop with uniform density can be sustained for thousands of periods, as shown in Rosswog & Price (2007). This shows a significant advantage for Lagrangian codes compared to Eulerian codes, which suffer from resolution-dependent advection errors when the configuration is not aligned to the grid. In the case of the density gradient, we can see that there is now a faster dissipation of the magnetic energy. The sudden reduction in magnetic energy is largely due to the reordering of the initial particle lattice near the density contrast. Comparing to Hopkins & Raives (2016), we can see that we have a similar reduction in magnetic energy as in the results from the MFM/MFV method at t = 20. There is a tiny difference (<1%) between GDSPH and TSPH, however, this owes itself to the initial reordering, after that we can see that the rate of change in the two discretizations are practically the same. The averaged normalized divergence error, ⟨ϵ_divB⟩, is around 10⁻² for both the Δ = 1 case and Δ = 2 case.

3.3. Brio-Wu shocktube

The Brio-Wu shocktube (Brio & Wu 1988) is an MHD extension to the classic Sod shocktube test (the hydro setup is the same). It tests how well the implementation can handle different MHD shocks, rarefactions, and contact discontinuities. We followed the setup from Hopkins & Raives (2016) and Price et al. (2018) and initialized a 3D thin periodic box, with a total region of [n_x, n_y, n_z]=[1024, 24, 24] and an active region of [256, 24, 24] particles on the left side x_L = [ − 2, 0], and a total region of [n_x, n_y, n_z]=[512, 12, 12] and an active region of [128, 12, 12] on the right side x_R = [0, 2]. The left state was set to:

$$\begin{array}{c}\hfill ({\rho }_L,P_L,v_x,v_y,v_z,B_x,B_y,B_z)=(1,1,0,0,0,0.75,1,0),\end{array}$$

and the right state was set to:

$$\begin{array}{c}\hfill ({\rho }_R,P_R,v_x,v_y,v_z,B_x,B_y,B_z)=(0.125,0.1,0,0,0,0.75,-1,0).\end{array}$$

The adiabatic index is set to γ = 2. We ran the simulation up to t = 0.2 and the results can be seen in Fig. 4. The GDSPH and TSPH results are shown with black and red dots, respectively, and the analytical solution is shown in blue lines. Our results are very similar to those from the PHANTOM code default case with the same resolution (Fig. 30 in Price et al. 2018). However, there is noticeable wall heating in the internal energy u in our test, due to the conservative thermal dissipation term used in this work. A more aggressive thermal dissipation can be added to smooth out the wall heating, which improves the results in the density, thermal and pressure profiles. However, this often leads to over dissipation in cases involving gravitational fields and is thus not a preferable choice. Divergence errors are kept low with a maximum value of ∼10⁻³ at the shock, and B_x remains close to constant, which also indicates excellent divergence control. Varying artificial resistivity parameter α_B from α_B = 0.5 to α_B = 1 only shows minimal differences, and as we can see from the results, α_B = 0.5 is sufficient to capture the magnetic field structure. We also note that using constant artificial viscosity (AV) parameters decreases post-shock oscillations and improves the results in the velocity profile. From Fig. 4, we can see that there is very little difference between GDSPH and TSPH in this test.

Fig. 4. Results from the Brio-Wu shocktube in 3D, with an initial left state (ρL, PL, vx, vy, vz, Bx, By, Bz) = (1, 1, 0, 0, 0, 0.75, 1, 0) and right state (ρR, PR, vx, vy, vz, Bx, By, Bz) = (0.125, 0.1, 0, 0, 0, 0.75, −1, 0). The figure shows the active region of the shock after t = 0.2, which contains about nx ≈ 300−400 particles across the x-direction. The blue line shows the reference solution and the black dots show the result from the GDSPH simulation, while red dots show the result from the TSPH simulation. There are minimal differences between the GDSPH and TSPH result.

3.4. Orszag-Tang vortex

The Orszag-Tang vortex test was introduced by Orszag & Tang (1979) and is a standard test of MHD schemes, as it involves the development of super-sonic turbulence and the interaction of the different MHD shocks. We set up a 3D thin periodic box with $L=(1.0,1.0,\frac{2\sqrt{6}}{n_x})$ at varying resolutions ([n_x, n_y, n_z]=[128, 148, 12],[256, 296, 12] and [512, 590, 12]). The test consists out of a velocity vortex:

$$\begin{array}{c}\hfill [v_x,v_y,v_z]=v_0[-sin(2\pi (y-y_{min})),sin(2\pi (x-x_{min})),0],\end{array}$$

and a doubly periodic magnetic field:

$$\begin{array}{c}\hfill [B_x,B_y,B_z]=B_0[-sin(2\pi (y-y_{min})),sin(4\pi (x-x_{min})),0],\end{array}$$

where v₀ = 1, $B_0=1/\sqrt{4\pi },x_{min}=-0.5$ and y_min = −0.5. Setting the initial plasma beta β₀ = 10/3, the initial Mach number M = v₀/c_s = 1 and the adiabatic index γ = 5/3, we get the initial pressure $P_0=B_0^2{\beta }_0/2=0.133$ and density ρ₀ = γP₀M₀ = 0.221. We show the results of the different resolution runs in Fig. 5 after t = 0.5 (top row) and t = 1 (bottom row). The test was run with both GDSPH and TSPH, however, we found only very small differences between them, which is why we only show the result from GDSPH. This can additionally be seen in Fig. 6, which show the time evolution of the magnetic energy in all of our test cases. From the result at t = 0.5 we can see that we reproduce the shock structure well and capture the trapped dense filament in the centre of the domain for all resolutions. With increasing resolution, the shock structure and filament become more defined. At t = 1, a more turbulent flow has developed. Our simulations capture most of the key features and compare well with previous works (for example, Fig. 32 in Price et al. 2018). However, it appears that our method is unable to develop the central magnetic island, a feature that is supposed to form when the current sheet in the center becomes unstable and reconnects due to the tearing mode instability. This is the case for most previous implementations of SPMHD, however, in Wurster et al. (2017) the authors argue that with a less dissipative artificial resistivity switch the magnetic island can be reproduced, and this motivate them to use the signal speed v_sig, B given in Eq. (19).

Fig. 5. Results from the Orszag-Tang vortex in 3D done with GDSPH, which shows rendered density slices (z = 0) at t = 0.5 (top) and t = 1 (bottom) for varying resolution [nx, ny, nz]=[128, 148, 12],[256, 296, 12] and [512, 590, 12] (low to high from left to right). The simulations capture well most of the key features for all tested resolutions. With increasing resolution the flows are more defined and show increased complexity. There are no significant differences between GDSPH and TSPH in this case.

We use the same switch but the problem remains. Nevertheless, if we compare our evolution of the magnetic energy (the magenta curve in Fig. 6) to theirs (the grey curve) we can see that our simulations are actually less dissipative, likely because we use a smaller resistivity parameter α_B = 0.5. Increasing α_B to 1 leads to more dissipation and brings the final energy closer to the Phantom run. It is thus likely that the development of the magnetic island also depends on the other dissipation terms such as AV and artificial conductivity. The mean normalized divergence error in the simulations are of the order ⟨ϵ_divB⟩ = 10^−3.5–10^−2.5, decreasing with higher resolution.

Fig. 6. Evolution of the total magnetic energy in units of the initial value for the 3D Orszag-Tang vortex test. We plot the result from three different tested resolutions [nx, ny, nz]=[128, 148, 12],[256, 296, 12],[512, 590, 12] done in GDSPH. We also include a TSPH case with nx = 256 and the nx = 512 curve from Wurster et al. (2017) for comparison sake. We can see that there are no visible difference between the TSPH (purple curve) and GDSPH (brown curve) cases. From the figure, it is also clear that the GDSPH case of nx = 512 (magenta curve) is less dissipative than the simulation from Wurster et al. (2017) (grey curve). Significant differences between resolutions can be seen to occur at later times during the evolution.

3.5. MHD rotor

The MHD rotor test was introduced by Balsara & Spicer (1999), which tests the propagation of Alfvén waves generated by a magnetized rotor. We followed the setup from Tóth (2000) and Price et al. (2018), and initialized a 3D thin periodic box of $L=(1.0,1.0,\frac{2\sqrt{6}}{n_x})$ at two resolutions, [n_x, n_y, n_z]=[128, 148, 12] and [256, 296, 12]. A rotating dense cylinder (ρ = 10) was initiated with cylindrical radius R = 0.1, within a surrounding medium (ρ = 1). The cylinder was put into rotation with an initial velocity of

$$\begin{array}{c}\hfill v=\frac{v_0}{r_{\mathrm{cyl}}}[-(y-y_0),(x-x_0),0]r_{\mathrm{cyl}}<R,\end{array}$$

where $r_{\mathrm{cyl}}=\sqrt{x^2+y^2}$ and v₀ = 2. The initial pressure was set to P = 1, with an adiabatic index of γ = 1.4. The initial magnetic field was set to $B=[5/\sqrt{4},0,0]$. The particles were set up on a closed packed lattice and the simulation were run until t = 0.15. The density contrast was unsmoothed, which means that there will be some noise at the boundary initially, due to particle reordering. The results of the simulations done with GDSPH and TSPH can be seen in Fig. 7, which shows 30 contours and the rendering of the magnetic pressure, with limits taken to be the same as in Tóth (2000) and Price et al. (2018). From the results, we can see that the difference between GDSPH and TSPH is generally small. However, we do find that in GDSPH there are notable increases in magnetic pressure at the pressure maxima compared to TSPH, which also seen in Hopkins & Raives (2016) when the authors compared MFM with SPH (their Fig. 15). In general, the key features of the test are captured by both methods. The mean normalized divergence errors in the simulations are of the order ⟨ϵ_divB⟩ = 10⁻⁴–10⁻³.

Fig. 7. Result from the magnetic rotor in 3D, which shows rendered magnetic pressure slices (z = 0) at t = 0.15, for varying resolution [nx, ny, nz]=[128, 148, 12] and [256, 296, 12] and for both GDSPH and TSPH. The plot also shows 30 contours with limits taken to be the same as Tóth (2000) and Price et al. (2018) (Pmag = [0, 2.642]) for a more direct comparison. We can see that GDSPH develops a larger and broader magnetic pressure peak compared to TSPH.

3.6. Magnetized blastwave

The magnetized blastwave was introduced by Balsara & Spicer (1999) and Londrillo & Del Zanna (2000), in which a central over-pressurized region expands preferentially along the magnetic fields lines. We followed the setup from Stone et al. (2008) and Price et al. (2018), and initialized a 3D periodic box of L = [1.0, 1.0, 1.0] with uniform density at a resolution of N = 256³. An inner region of radius R = 0.125 was over-pressurized (P_in = 100) to 100 times the outer pressure, which was set to P_out = 1. The adiabatic index of the gas is set to γ = 5/3⁵. The initial magnetic field was set to:

$$\begin{array}{c}\hfill B=[10/\sqrt{2},0,10/\sqrt{2}].\end{array}$$

This sets the initial plasma beta to β_in = 2 in the inner region and β_out = 0.02 in the outer region. The simulation was run for t = 0.02 and the result can be seen in Fig. 8. The rendering and limits in the figure were set to the same as the results presented in Tóth (2000) and Price et al. (2018), so that they can be directly compared. We can see that our results agree well with the previous authors, capturing the inner and outer structure of the blast well. There are minimal differences between the GDSPH and TSPH results. The mean normalized divergence error in the simulations are of the order ⟨ϵ_divB⟩ = 10⁻⁵.

Fig. 8. Result from the magnetized blast in 3D in the 2563 GDSPH run, which shows rendered slices of different fluid quantities at t = 0.02. To the top left we can see the density rendering, top right the kinetic energy density, bottom left the thermal pressure and bottom right the magnetic pressure. The limits are taken to be the same as Stone et al. (2008) and Price et al. (2018) for a direct comparison: ρ = [0.19, 2.98], Ekin = [0, 33.1], P = [1, 42.4] and Pmag = [25.2, 65.9].

3.7. Kelvin-Helmholtz instability in MHD

Generation of the Kelvin-Helmholtz instability is crucial for efficient mixing in hydrodynamical simulations. Here, we look at the same problem but with a magnetic field applied parallel to the flow. This has a stabilizing effect on the shear flow due to the magnetic tension force. We followed the setup from McNally et al. (2012), but extend it to 3D, making a thin periodic box ($L=[1.0,1.0,\frac{2\sqrt{6}}{n_x}]$), with a resolution of [n_x, n_y, n_z]=[256, 296, 12]. We applied a uniform pressure of P = 5/2 with an adiabatic index of γ = 5/3. The hot outer stream has a density of ρ_out = 1 and velocity v_out = [ − 0.5, 0, 0]. The cold inner stream has a density of ρ_in = 2 and velocity v_in = [0.5, 0, 0]. A uniform magnetic field was added in the direction of the flow velocity B = [0.1, 0, 0].

The results for TSPH and GDSPH at t = 1.6 and 3.2 can be seen in Fig. 9. And in Fig. 10 we show the particle distribution of the surface boundary at t = 3.2. The TSPH result exhibits a very gloopy behaviour and shows a decreased growth of the KH mode. A strong artificial surface tension effect can clearly be seen between the hot and the cold phase in Fig. 10. With GDSPH this effect is largely eliminated and the growth rate improves significantly. This large improvement in GDSPH lends itself mainly to the reduction of the leading order errors, which we discussed in Sect. 2.2. Adding turbulent diffusion (that is, with the code default parameters) further improves the result, because it allows particles to effectively mix or reorder (as shown clearly in the particle distribution near the boundary regions in Fig. 10), but the growth rate remains similar to the GDSPH-only case. We note that the sharp contact discontinuities in the initial condition are not smoothed, unlike in McNally et al. (2012). We choose this because it represents an extreme situation for SPH where the initial particle ordering is not optimal (that is, zeroth order errors are relatively high), and we show that GDSPH performs well even in this extreme case. We also ran the setup using a smoothed contact discontinuity, and this is shown in the rightmost column in Fig. 9. The magnetic field effectively uncoils and stretches the vortex, which is in good agreement to the results shown in Hopkins & Raives (2016) with the same setup. Here TSPH and GDSPH develop indistinguishably until later time, where at the end only small differences can be seen. The mean normalized divergence error in the simulations are of the order of ⟨ϵ_divB⟩ = 10⁻³.

Fig. 9. Result from the magnetic Kelvin-Helmholtz instability in 3D, which shows rendered density slices (z = 0) at t = 1.6 (top) and t = 3.2 (bottom), for TSPH without diffusion (left), GDSPH without diffusion (middle left), GDSPH with diffusion (middle right) and GDSPH with an initial smoothed contact discontinuity (right). The TSPH result exhibits a very gloopy behaviour and shows decreased growth of the KH mode. This is mainly due to the artificial surface tension effect (see Fig. 10). GDSPH shows much better growth, the addition of diffusion only slightly improves the growth rate. The main effect from the magnetic field can be seen in all cases, which uncoils and stretches the vortex. The smoothed result closely resembles the MFM and grid result from Hopkins & Raives (2016).

Fig. 10. Surface boundary between the high and low density region in the magnetic Kelvin-Helmholtz instability at t = 3.2. The effect of the numerical surface tension can clearly be seen in the TSPH case, while GDSPH does not suffer from this issue. Adding thermal diffusion allows for local mixing between the cold and hot phase.

4. Collapse of a magnetized cloud

In this Section, we apply our method to an astrophysical problem and consider the collapse of a magnetized cloud. In this type of problem involving large dynamic scales, we see a substantial difference between GDSPH and TSPH. A rotating magnetized cloud is allowed to collapse under its own gravity. During the collapse, the cloud is compressed over several orders of magnitude, testing how the magnetic field develops and interacts with the gas during compression. The large-scale collapse is eventually halted by the formation of a pseudo-disk⁶, which then starts to slowly contract via magnetic braking. The collapse continues within the central region and as the first hydrostatic core starts to form, the magnetic field is twisted until it eventually launches a jet (Uchida & Shibata 1986; Lynden-Bell 1996; Ustyugova et al. 2000; Nakamura & Meier 2004). The formation and subsequent evolution of the first hydrostatic core stalls the collapse and a slow contraction phase begins. In this paper, we do not run the simulations far beyond the time of jet launching. The two main jet launching mechanisms are the magneto-centrifugal and the magnetic pressure driven mechanism. With a global poloidal magnetic field as in our model, both of these mechanisms play an important role. The resulting magnetic field structure of the jet consists of a poloidal dominated central core with a surrounding toroidal field which produces a strong current along the jet. We refer to this magnetic field structure as the magnetic tower throughout this paper. All these key aspects require the code to have excellent angular momentum conservation, small numerical dissipation and maintain low divergence errors (∇ ⋅ B).

We followed the setup outlined in Hennebelle & Fromang (2008) and Hopkins & Raives (2016) and set up a 3D periodic box L = [0.15 pc, 0.15 pc, 0.15 pc]. A cloud was initiated with a radius of $R_{\mathrm{c}}=0.015\mathrm{pc}$ and a mass of 1 solar mass (M_c = 1 M_⊙), within a surrounding medium that has 360 times lower density than the cloud (ρ_out = M_c/(360V_c)). The cloud was put in rotation with an orbital time of P = 4.7 × 10⁵ yr, which corresponds to a kinetic over potential energy ratio of about E_K/E_P ≈ 0.045. This is a higher ratio compared to the peak value of 0.02 from the observed distribution of rotation rates in molecular clouds (E_K/E_P ∈ (0.002, 1.4)) (Goodman et al. 1993). A constant magnetic field B₀ was initialized in the direction of the angular momentum vector ($\widehat{\mathit{z}}$), and we varied the strength in accordance to different mass-to-flux ratios. The mass-to-flux ratio μ is relative to the critical mass-to-flux ratio, (M_c/Φ)_crit, in which the cloud is fully supported by magnetic forces against gravity, that is,

$$\begin{array}{cc}& \mu =\left(\frac{M_{\mathrm{c}}}{\mathrm{\Phi }}\right){\left(\frac{M_{\mathrm{c}}}{\mathrm{\Phi }}\right)}_{\mathrm{crit}},\hfill \end{array}$$(34)

$$\begin{array}{cc}& \left(\frac{M_{\mathrm{c}}}{\mathrm{\Phi }}\right)=\frac{M_{\mathrm{c}}}{\pi R_{\mathrm{c}}^2{B}_0},{\left(\frac{M_{\mathrm{c}}}{\mathrm{\Phi }}\right)}_{\mathrm{crit}}=\frac{c_1}{3\pi }\sqrt{\frac{5}{G}}.\hfill \end{array}$$(35)

Here, c₁ = 0.53 is a parameter that can be determined numerically (Mouschovias & Spitzer 1976). We then get the corresponding initial magnetic field:

$$\begin{array}{c}\hfill B_0=\frac{610}{\mu }[\mu G].\end{array}$$(36)

The thermal pressure is determined by the following barotropic EOS,

$$\begin{array}{c}\hfill P=c_{s,0}^2\rho \sqrt{1+{(\rho /{\rho }_0)}^{4/3}},\end{array}$$(37)

with ρ₀ = 10⁻¹⁴ g cm⁻³ and c_s, 0 = 0.2 km s⁻¹. We looked at six different magnetic flux ratio values in our simulations, from weak to high (μ = ∞, 75, 20, 10, 5, 2). These were run with a moderate resolution of 50³ in the cloud, which corresponds to about 40³ particles in the low density medium, same as in the setup of Hopkins & Raives (2016). These six cases were run with both GDSPH and TSPH until the core has fully collapsed, close after the time of jet launching, which typically occurs when the maximum density hits a value in between ρ = 10⁻¹² ↔ 10⁻¹¹ g cm⁻³. This occurs near the free fall time $t_{\mathit{ff}}=\sqrt{\frac{3}{2\pi G\rho }}\approx 4\times {10}^4\mathrm{yr}$, at around t = 1.05t_ff ↔ 1.3t_ff depending on resolution/initial magnetic field strength. No sink particles were used in any of our simulations. The results of these simulations can be seen in Fig. 11 (GDSPH) and Fig. 12 (TSPH).

Fig. 11. Result of the magnetized cloud collapse for GDSPH at a resolution of 503 particles with varying magnetic flux ratio μ going left to right from high to low. We show figures at the time of jet formation (around t = tff), which occurs due to the winding of the magnetic field during the collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. The pure hydrodynamic run (μ = ∞) of GDSPH becomes gravitationally unstable and is very similar to that of TSPH in Fig. 12. We can see that a jet is launched in the cases of μ = 75, 20, 10, 5 while in the case of μ = 2 the interchange instability (see Fig. 16) disrupts the disk before jet launching.

Fig. 12. Result of the magnetized cloud collapse for TSPH at a resolution of 503 particles with varying magnetic flux ratio μ going left to right from high to low. We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. The pure hydrodynamic run (μ = ∞) of TSPH becomes gravitationally unstable and is very similar to that of GDSPH in Fig. 11. We can see that TSPH does not form a jet in any of the weak field cases (μ = 75, 20, 10) and there is no sign of a magnetic tower being formed. In the case of μ = 5, we can see a jet being launched, where a current dominated magnetic tower has formed, however, the central part of the tower has been completely quenched. The μ = 2 case also launches a jet, but collapses faster than in the high resolution case, which effectively leads to easier jet formation.

The pure hydrodynamic runs (μ = ∞) of both GDSPH and TSPH become gravitationally unstable and the resulting evolution is very similar (see Figs. 11 and 12). For GDSPH, we see in Fig. 11 that jet launching can be seen in the weak field regime (μ = 75, 20, 10), although it is very short-lived in the case of μ = 75. It is clear from the poloidal magnetic field and the current density (the fourth and fifth row in Fig. 11) that we have a developed magnetic tower in all these three cases. This is a remarkable achievement, especially for SPMHD, since the amplification of the magnetic field can easily be quenched by numerical dissipation. The jet strength, morphology and velocity structure resemble those in Hopkins & Raives (2016) with the same resolution using the MFM/MFV method with more complex gradient operators and Riemann solvers. In contrast, for TSPH (Fig. 14) we see neither jet launching nor a magnetic tower in the weak-field regime with the fiducial resolution of 50³. This is largely due to numerical dissipation which suppresses field amplification and hinders the formation of a jet. At the time of jet launching, fragmentation of the disc occurs in the two weakest cases μ = 75 and μ = 20 for TSPH, as the magnetic field is too weak to support the disc. For GDSPH, it does not occur until a later time in the simulation, however, the exact time of fragmentation is heavily dependent on other dissipation terms such as artificial viscosity.

For the more magnetized case with μ = 5 we can see that TSPH successfully develops a jet, however, closer inspection on |B_z| in Fig. 12 shows that the central portion of the magnetic tower is much less developed, with an order of magnitude smaller strength than the same run with GDSPH. Jet launching is also seen in the GDSPH case with μ = 5 with the magnetic tower intact, albeit weaker and less collimated than in the μ = 10 case. For μ = 2, the collapse proceeds differently in GDSPH and TSPH from an early stage. This can be seen in Fig. 15, which shows the density structure of the collapsing cloud in an early stage with different resolutions. As the magnetic field is very strong, accretion will occur primarily along the field lines, creating an elongated cloud structure. While at high resolution (250³) both GDSPH and TSPH runs converge to the correct structure, at the fiducial resolution of 50³, only GDSPH shows an elongated cloud. The cloud collapses faster in TSPH, likely due to excessive dissipation.

For μ = 2, we can see that GDSPH does not produce a coherent jet (Fig. 11). This is mainly due to the disk being disrupted by the magnetic interchange instability. This instability occurs due to an accumulation of magnetic flux near the accreting protostar, where the magnetic flux that would have been dragged into the protostellar core is redistributed to the surrounding medium by dissipative effects (AR). This builds a large magnetic pressure gradient which, together with the twisted magnetic field, eventually launches highly magnetized bubbles in the azimuthal direction. A density rendering together with a velocity map of the μ = 2 case after the launch of the magnetized bubbles is shown in Fig. 16. This is similar to the results seen in simulations using the AMR code ENZO (Krasnopolsky et al. 2012) and in SPH simulations with PHANTOM (Wurster et al. 2017). We would like to stress that, unlike the SPH runs in Hopkins & Raives (2016), the disk disruption is not due to divergence errors, but instead a consequence of the magnetic dissipation. At later times (for example, right panel of Fig. 16), we can see that the protostellar core remains centrally located, which indicates good divergence control and angular momentum conservation. As the formation of the interchange instability is driven by the redistribution of magnetic flux, it can depend heavily on the choice of AR prescription and the use of sink particles. Wurster et al. (2017) observed similar magnetic bubbles with the same AR prescription as ours, while other tested AR prescriptions did not launch magnetic bubbles. However, all other works that produce interchange instabilities use sink particles, which can artificially redistribute the flux as matter is accreted by the sink, while leaving the magnetic field close to the sink intact. The development of the interchange instability in our simulations without sink particles might indicate a more physical origin of the effect. Additional work will need to be done to determine if this is in fact a real effect or a consequence of the numerical scheme.

To investigate the convergence with resolution, we simulated the μ = 10 case across different resolutions (12³, 25³, 50³, 100³ and 250³) for both methods. The results are shown in Fig. 13 for GDSPH and Fig. 14 for TSPH. In GDSPH, resolved jet structures and fully developed magnetic towers are already evident in cases with 25³ resolution, and which increase in complexity as we increase the resolution. A weak outflow appears even in the lowest resolution of 12³. In contrast, the runs with TSPH shows slow convergence. The structure of the magnetic field is severely distorted, and magnetic tower and proper collimated jet are not developed in all cases except the highest resolution. Again, we note that our GDSPH results are very comparable to the MFM/MFV runs in Hopkins & Raives (2016), both in terms of jet properties and converging speed. The magnetic tower structure is also qualitatively similar to the cloud collapse simulation in the weak field regime from the moving-mesh code AREPO (Pakmor et al. 2011), although with a slightly different initial setup. We should note that the collapse of the 12³ is artificially suppressed and contract much slower then what is expected. This is because the local Jeans mass is not fully resolvable in these simulations. Bate & Burkert (1997) estimated that around 3 × 10⁴ particles where required to correctly resolve the Jeans mass in similar collapse cases⁷. The effect can partly be seen in the 25³ case as well, especially at later times. However, in this case, the cloud collapsing has a similar evolution up to the time of jet-launching as the higher resolution cases.

Fig. 13. Result of the resolution study of the magnetized cloud collapse for GDSPH with μ = 10. We vary the resolution from left to right, in the initial cloud (123, 253, 503, 1003, 2503) and medium (103, 203, 403, 803, 2003). We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2]; all quantities are shown in logarithmic scale. Jet formation and a proper magnetic tower can be seen to occur at very low resolution compared to TSPH. The jet structure and magnetic tower further increases in complexity as we increase the resolution.

Fig. 14. Result of the resolution study of the magnetized cloud collapse for TSPH with μ = 10. We vary the resolution from left to right, in the initial cloud (123, 253, 503, 1003, 2503) and medium (103, 203, 403, 803, 2003). We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. We can see that TSPH only forms a collimated jet and a proper magnetic tower at the highest resolution.

Fig. 15. Density rendered slice ([g cm−3]) through the rotation axis, showing the early cloud structure of the strong field case μ = 2 before the formation of the disk. As the magnetic field is very strong, accretion will occur primarily along the field lines, creating an elongated cloud structure. Both high resolution cases (2503) of TSPH and GDSPH creates an elongated cloud structure, while in the low-resolution case (503) only GDSPH forms the same cloud structure. TSPH instead forms a more compact central cloud, which as a consequence collapses faster than the GDSPH case and the high resolution cases.

Fig. 16. Magnetic interchange instability in the strong field case (μ = 2) for the GDSPH simulation with 503 resolution. Left panel: zoom-in of the disc structure seen face-on in Fig. 11 and right panel: same region at a later time. The white arrows show the direction of velocity and the colour scale indicates density. The instability launches magnetized bubbles in the azimuthal direction. At later times we can see that the central star starts to accrete again along the filamentary structure. These figures can be compared to the results from Krasnopolsky et al. (2012).

5. Discussion

In this paper, we present an SPMHD method, which utilizes the Geometric Density average force discretization (GDSPH) of the MHD equations. GDSPH has been shown in previous work to greatly improve the accuracy near density discontinuities and eliminate the surface tension problem. We show that MHD also benefits from this method. For a large part, the standard test problems (Sects. 3.1–3.6) both GDSPH and TSPH handle the problems very well, and the differences between the two methods are minimal. However, when the problem involves mixing such as in the case of Kevin-Helmholtz instabilities, GDSPH shows clear advantages. This is somewhat expected, and in agreement with earlier studies without magnetic fields (Wadsley et al. 2017). However, when we apply the method in the astrophysical test of a collapsing magnetized cloud, we see that GDSPH leads to a significant improvement. GDSPH not only realistically captures the development of a magnetic tower and jet launching in the weak field regime (μ ≳ 10), but also exhibits a fast improvement in the complexity and structure of the jet with increased resolution. In contrast, TSPH only manages to launch jets in the strong field regime with a resolution of 50³ and only develop a collimated jet in the highest resolution runs of the μ = 10 case. We also show that, in the strong field regime, GDSPH converges better than TSPH in accretion time and in the outer cloud structure. The results of TSPH is in agreement with previous studies using TSPH (Hopkins & Raives 2016; Bürzle et al. 2011) which used the same IC as this work, and other studies (Price & Bate 2007; Price et al. 2012) which employ a smaller initial cloud rotation (E_K/E_P = 0.005).

Overall, our new method shows improved or comparable results to previous SPMHD implementations such as in PHANTOM (Price et al. 2018) and in GIZMO (Hopkins & Raives 2016). In many test cases and particularly in the cloud collapse case, GDSPH produces qualitatively very similar result to that of the MFM/MFV method and achieves similar convergence speed. It is worth noting that all the simulations are run in 3D with the code default parameters listed in Sect. 3 without any adjustment by hand to specific problems. The success of GDSPH can most likely be attributed to the reduction of SPH “E₀ errors” (Eq. (10)) and linear errors by the geometric density average force formulation. As discussed in Sect. 2.2, the advantage of such discretization is more evident when the density gradient is large, such as in the cloud collapse case.

Using the constrained hyperbolic divergence cleaning scheme with variable cleaning speed from Tricco et al. (2016), we can keep the divergence error low in all cases. The mean normalized divergence error, ⟨ϵ_divB⟩=⟨h|∇ ⋅ B|/|B|⟩, is typically of order 10⁻⁵ − 10⁻³. In Fig. 17, we show the normalized divergence error maps for several test problems. Again we see that the divergence cleaning works extremely well here, the maximum error is generally around 10⁻². Comparing to Hopkins (2016, their Fig. 4), we find that the errors are smaller than their MFM simulations with the Dedner et al. (2002) cleaning in general, with an exception for the outskirt of the advection loop (where the magnetic field is essentially zero and thus not important for the result). The improvement is probably due to the more advanced constrained cleaning method (Tricco et al. 2016). The normalized divergence error for the μ = 10 cloud-collapse case at the jet launching time is shown in Fig. 18. Here, the divergence cleaning still performs very well in the disc, along the jets and for the majority of the regions where the outflow interacts with the ambient gas, especially when the divergence error is compared to the total gas pressure (right panel). The result is similar to the Dedner cleaning in Hopkins (2016), although our error is somewhat larger at the tip of the jets where the gas is shocked. However, we note that the comparison is not direct in this case as the jets may develop differently. Overall, the result from cleaning is still worse than the constrained transport or constrained gradient schemes (Hopkins 2016). For SPMHD, as shown in Tricco & Price (2012), divergence errors can be reduced to machine precision (or more practically to a certain tolerance value) using cleaning, with the help of a sub-cycling routine. However, local adjustments are required to determine the number of iterations for each particle to efficiently subcycle the cleaning in the simulation. This is because certain regions are more affected than others and because divergence is spread to nearby neighbours. Conceivably, if vector potentials (Stasyszyn & Elstner 2015) could work for a wider range of problems this could be an interesting avenue as well. However, the exploration of these methods in detail is beyond the scope of this work.

Fig. 17. Normalized divergence error (ϵdivB ≡ h|∇ ⋅ B|/|B|, Eq. (29)) of the different cases (all in 3D) with GDSPH. From left to right: first: advection of a current loop with a density ratio of Δ = 2 at t = 20 (Sect. 3.2), second: MHD rotor (nx = 256) at t = 0.15 (Sect. 3.5). Third: Orszag-Tang vortex (nx = 256) at t = 0.5 (Sect. 3.4). The MHD blastwave (n = 2563) at t = 0.02 (Sect. 3.5). The MHD Kelvin-Helmholtz instability (nx = 256) with an initial sharp contact discontinuity at t = 3.2 (Sect. 3.7).

Fig. 18. Normalized divergence error of the magnetized cloud collapse simulation with GDSPH at the jet launching time, t ≈ tff. The mass-to-flux ratio is μ = 10 and the resolution is 1003. Left: normalized divergence error as in (Eq. (29)), ϵdivB ≡ h|∇ ⋅ B|/|B|; right: divergence error normalized to total gas pressure, ${ϵ}_{\mathrm{div}B}/\sqrt{(1+\beta )}=h|\mathrm{\nabla }·\mathbf{B}|/\sqrt{({\mathbf{B}}^2+2P)}$, where β = 2P/B2 is the plasma beta.

From the tests, we can see that using a lower artificial resistivity coefficient (α_B = 0.5) than (α_B = 1) (Price et al. 2018) works well for all cases. However, the choice of artificial resistivity switch still remains somewhat ad-hoc as it is difficult to accurately detect the MHD discontinuity. A Godunov SPH scheme (Iwasaki & Inutsuka 2011) solves this problem by replacing artificial resistivity with Riemann solvers, which have been shown to produce minimal artificial diffusivity. However, this brings with it an increase in computational cost and it is unclear if the extra cost is worth it. Further improvements in the convergence of SPMHD may be found in the use of integral-based gradient estimates, which have been shown to be more accurate and less noisy than the standard SPH gradient estimate (García-Senz et al. 2012; Rosswog 2015; Cabezón et al. 2017). This could be especially beneficial for modeling subsonic turbulent flows (Valdarnini 2016). This gradient estimate can easily be implemented within the GDSPH framework and will be investigated in future work.

Meshless methods (SPH and MFM) were recently explored in local 3D simulations of the magnetorotational instability (MRI; Deng et al. 2019). The authors found that in the vertically stratified MRI simulations, SPH developed an unphysical state with strong toroidal field components and no sustained turbulence. In a forthcoming paper (Wissing et al., in prep.) we will show that GDSPH do not show this unphysical behavior and that they reproduce the characteristic periodic azimuthal magnetic field pattern (butterfly diagram) of the stratified MRI.

Acknowledgments

We thank the anonymous referee for the useful comments and suggestions which have improved the quality and clarity of the paper. We thank James Wadsley, Tom Quinn, Max Grönke and Ben Keller for very helpful and interesting discussions. The simulations were performed using the resources from the National Infrastructure for High Performance Computing and Data Storage in Norway, UNINETT Sigma2, allocated to Project NN9477K. We also acknowledge the support from the Research Council of Norway through NFR Young Research Talents Grant 276043.

References

All Figures

Fig. 1. Results of the 3D circularly polarized Alfvén wave test. Top panel: transverse component of the magnetic field in the direction of propagation after five periods. The analytical/initial solution is plotted in black, and the simulation results with resolution [nx, ny, nz]=[128, 74, 78] in red, [nx, ny, nz]=[64, 36, 39] in green, and [nx, ny, nz]=[32, 18, 18] in blue. Both of these are with Wendland C4 kernel with 200 neighbours. Bottom panel: convergence study for the Alfvén wave test using different kernels and neighbour numbers. Shows how the L1 error scales with resolution (particles along the x-axis). The code default (Wendland C4 kernel with 200 neighbours) is shown in green, Wendland C4 kernel with 114 neighbours are shown in magenta, Quintic kernel with 114 neighbours are shown in blue and the dashed brown line shows the curve for second order convergence. Convergence towards the analytical solution for all kernels are close to second order. When the smoothing length becomes comparable to half the wave length of the Alfvén wave, the MHD gradients becomes more ill defined which causes the slower convergence speed for the Wendland kernel (${\overline{N}}_{\mathrm{neigh}}=200$) at low resolution.

In the text

Fig. 2. Results from the advection loop test in 3D, showing a rendering of |B| in units of the initial value |B0|. Top left panel: initial setup (same rendering for nx = 128 and nx = 256, and top right: uniform density case for GDSPH with resolution nx = 128 after twenty crossings t = 20. We can see that the current loop is conserved almost perfectly even with all the dissipation terms turned on. This shows an improvement comparing to the SPMHD results in PHANTOM, where dissipation is seen after five periods (Fig. 35 in Price et al. 2018), which plots the current density. The bottom two panels show the cases with a density gradient $\mathrm{\Delta }\equiv \frac{{\rho }_{\mathrm{in}}}{{\rho }_{\mathrm{out}}}=2$ between the inner loop and outer medium and with resolution of nx = 256, after twenty crossings t = 20. The bottom left panel shows the GDSPH case and the bottom right show the TSPH case (nx = 256). Both show similar dissipation of the magnetic field to the results from the MFM/MFV method in GIZMO (Hopkins & Raives 2016).

In the text

Fig. 3. Results from the advection loop test in 3D, showing the time evolution of the magnetic field energy in units of the initial value. After t = 20 our uniform case has dissipated about 0.3%, while the cases with the density gradient Δ = 2 have dropped around 5%, which is similar to the dissipation displayed by the MFM/MFV method in Hopkins & Raives (2016). There is a tiny difference between GDSPH and TSPH, but this owes itself to differences in the initial reordering.

In the text

Fig. 4. Results from the Brio-Wu shocktube in 3D, with an initial left state (ρL, PL, vx, vy, vz, Bx, By, Bz) = (1, 1, 0, 0, 0, 0.75, 1, 0) and right state (ρR, PR, vx, vy, vz, Bx, By, Bz) = (0.125, 0.1, 0, 0, 0, 0.75, −1, 0). The figure shows the active region of the shock after t = 0.2, which contains about nx ≈ 300−400 particles across the x-direction. The blue line shows the reference solution and the black dots show the result from the GDSPH simulation, while red dots show the result from the TSPH simulation. There are minimal differences between the GDSPH and TSPH result.

In the text

Fig. 5. Results from the Orszag-Tang vortex in 3D done with GDSPH, which shows rendered density slices (z = 0) at t = 0.5 (top) and t = 1 (bottom) for varying resolution [nx, ny, nz]=[128, 148, 12],[256, 296, 12] and [512, 590, 12] (low to high from left to right). The simulations capture well most of the key features for all tested resolutions. With increasing resolution the flows are more defined and show increased complexity. There are no significant differences between GDSPH and TSPH in this case.

In the text

Fig. 6. Evolution of the total magnetic energy in units of the initial value for the 3D Orszag-Tang vortex test. We plot the result from three different tested resolutions [nx, ny, nz]=[128, 148, 12],[256, 296, 12],[512, 590, 12] done in GDSPH. We also include a TSPH case with nx = 256 and the nx = 512 curve from Wurster et al. (2017) for comparison sake. We can see that there are no visible difference between the TSPH (purple curve) and GDSPH (brown curve) cases. From the figure, it is also clear that the GDSPH case of nx = 512 (magenta curve) is less dissipative than the simulation from Wurster et al. (2017) (grey curve). Significant differences between resolutions can be seen to occur at later times during the evolution.

In the text

Fig. 7. Result from the magnetic rotor in 3D, which shows rendered magnetic pressure slices (z = 0) at t = 0.15, for varying resolution [nx, ny, nz]=[128, 148, 12] and [256, 296, 12] and for both GDSPH and TSPH. The plot also shows 30 contours with limits taken to be the same as Tóth (2000) and Price et al. (2018) (Pmag = [0, 2.642]) for a more direct comparison. We can see that GDSPH develops a larger and broader magnetic pressure peak compared to TSPH.

In the text

Fig. 8. Result from the magnetized blast in 3D in the 2563 GDSPH run, which shows rendered slices of different fluid quantities at t = 0.02. To the top left we can see the density rendering, top right the kinetic energy density, bottom left the thermal pressure and bottom right the magnetic pressure. The limits are taken to be the same as Stone et al. (2008) and Price et al. (2018) for a direct comparison: ρ = [0.19, 2.98], Ekin = [0, 33.1], P = [1, 42.4] and Pmag = [25.2, 65.9].

In the text

Fig. 9. Result from the magnetic Kelvin-Helmholtz instability in 3D, which shows rendered density slices (z = 0) at t = 1.6 (top) and t = 3.2 (bottom), for TSPH without diffusion (left), GDSPH without diffusion (middle left), GDSPH with diffusion (middle right) and GDSPH with an initial smoothed contact discontinuity (right). The TSPH result exhibits a very gloopy behaviour and shows decreased growth of the KH mode. This is mainly due to the artificial surface tension effect (see Fig. 10). GDSPH shows much better growth, the addition of diffusion only slightly improves the growth rate. The main effect from the magnetic field can be seen in all cases, which uncoils and stretches the vortex. The smoothed result closely resembles the MFM and grid result from Hopkins & Raives (2016).

In the text

	Fig. 10. Surface boundary between the high and low density region in the magnetic Kelvin-Helmholtz instability at t = 3.2. The effect of the numerical surface tension can clearly be seen in the TSPH case, while GDSPH does not suffer from this issue. Adding thermal diffusion allows for local mixing between the cold and hot phase.
In the text

Fig. 11. Result of the magnetized cloud collapse for GDSPH at a resolution of 503 particles with varying magnetic flux ratio μ going left to right from high to low. We show figures at the time of jet formation (around t = tff), which occurs due to the winding of the magnetic field during the collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. The pure hydrodynamic run (μ = ∞) of GDSPH becomes gravitationally unstable and is very similar to that of TSPH in Fig. 12. We can see that a jet is launched in the cases of μ = 75, 20, 10, 5 while in the case of μ = 2 the interchange instability (see Fig. 16) disrupts the disk before jet launching.

In the text

Fig. 12. Result of the magnetized cloud collapse for TSPH at a resolution of 503 particles with varying magnetic flux ratio μ going left to right from high to low. We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. The pure hydrodynamic run (μ = ∞) of TSPH becomes gravitationally unstable and is very similar to that of GDSPH in Fig. 11. We can see that TSPH does not form a jet in any of the weak field cases (μ = 75, 20, 10) and there is no sign of a magnetic tower being formed. In the case of μ = 5, we can see a jet being launched, where a current dominated magnetic tower has formed, however, the central part of the tower has been completely quenched. The μ = 2 case also launches a jet, but collapses faster than in the high resolution case, which effectively leads to easier jet formation.

In the text

Fig. 13. Result of the resolution study of the magnetized cloud collapse for GDSPH with μ = 10. We vary the resolution from left to right, in the initial cloud (123, 253, 503, 1003, 2503) and medium (103, 203, 403, 803, 2003). We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2]; all quantities are shown in logarithmic scale. Jet formation and a proper magnetic tower can be seen to occur at very low resolution compared to TSPH. The jet structure and magnetic tower further increases in complexity as we increase the resolution.

In the text

Fig. 14. Result of the resolution study of the magnetized cloud collapse for TSPH with μ = 10. We vary the resolution from left to right, in the initial cloud (123, 253, 503, 1003, 2503) and medium (103, 203, 403, 803, 2003). We show figures at the time of jet formation (around t = tff), which occur due to the winding of the magnetic field during collapse, which produces a magnetic tower structure. The top row shows a rendered face-on slice (Lxy = [2000 AU, 2000 AU]) of the density [g cm−3], the rest of the rows show rendered slices through the rotation axis (Lxz = [2000 AU, 2000 AU]), where the second shows density [g cm−3], the third show radial velocity [km s−1], the fourth show the absolute poloidal magnetic field [μG] and the fifth shows the current density [A m−2], all quantities are shown in logarithmic scale. We can see that TSPH only forms a collimated jet and a proper magnetic tower at the highest resolution.

In the text

Fig. 15. Density rendered slice ([g cm−3]) through the rotation axis, showing the early cloud structure of the strong field case μ = 2 before the formation of the disk. As the magnetic field is very strong, accretion will occur primarily along the field lines, creating an elongated cloud structure. Both high resolution cases (2503) of TSPH and GDSPH creates an elongated cloud structure, while in the low-resolution case (503) only GDSPH forms the same cloud structure. TSPH instead forms a more compact central cloud, which as a consequence collapses faster than the GDSPH case and the high resolution cases.

In the text

Fig. 16. Magnetic interchange instability in the strong field case (μ = 2) for the GDSPH simulation with 503 resolution. Left panel: zoom-in of the disc structure seen face-on in Fig. 11 and right panel: same region at a later time. The white arrows show the direction of velocity and the colour scale indicates density. The instability launches magnetized bubbles in the azimuthal direction. At later times we can see that the central star starts to accrete again along the filamentary structure. These figures can be compared to the results from Krasnopolsky et al. (2012).

In the text

Fig. 17. Normalized divergence error (ϵdivB ≡ h|∇ ⋅ B|/|B|, Eq. (29)) of the different cases (all in 3D) with GDSPH. From left to right: first: advection of a current loop with a density ratio of Δ = 2 at t = 20 (Sect. 3.2), second: MHD rotor (nx = 256) at t = 0.15 (Sect. 3.5). Third: Orszag-Tang vortex (nx = 256) at t = 0.5 (Sect. 3.4). The MHD blastwave (n = 2563) at t = 0.02 (Sect. 3.5). The MHD Kelvin-Helmholtz instability (nx = 256) with an initial sharp contact discontinuity at t = 3.2 (Sect. 3.7).

In the text

Fig. 18. Normalized divergence error of the magnetized cloud collapse simulation with GDSPH at the jet launching time, t ≈ tff. The mass-to-flux ratio is μ = 10 and the resolution is 1003. Left: normalized divergence error as in (Eq. (29)), ϵdivB ≡ h|∇ ⋅ B|/|B|; right: divergence error normalized to total gas pressure, ${ϵ}_{\mathrm{div}B}/\sqrt{(1+\beta )}=h|\mathrm{\nabla }·\mathbf{B}|/\sqrt{({\mathbf{B}}^2+2P)}$, where β = 2P/B2 is the plasma beta.

In the text