Astrophysical WPMHD

Astrophysical Weighted Particle Magnetohydrodynamics

Evghenii Gaburov, Keigo Nitadori
Leiden Observatory, Leiden University, the Netherlands
RIKEN, Tokyo, Japan

This paper presents applications of weighted meshless scheme for conservation laws to the Euler equations and the equations of ideal magnetohydrodynamics. The divergence constraint of the latter is maintained to the truncation error by a new meshless divergence cleaning procedure. The physics of the interaction between the particles is described by an one-dimensional Riemann problem in a moving frame. As a result, necessary diffusion which is required to treat dissipative processes is added automatically. As a result, our scheme has no free parameters that controls the physics of inter-particle interaction, with the exception of the number of the interacting neighbours which control the resolution and accuracy. The resulting equations have the form similar to SPH equations, and therefore existing SPH codes can be used to implement the weighed particle scheme. The scheme is validated in several hydrodynamic and MHD test cases. In particular, we demonstrate for the first time the ability of a meshless MHD scheme to model magneto-rotational instability in accretion disks.

1 Introduction

Computational magnetohydrodynamics (MHD) remains an important tool to understand complex behaviour of astrophysical plasmas. While many methods have been developed to solve equations of ideal MHD of Eulerian (cartesian) meshes, few successful Lagrangian meshless formulations exists. The latter, however, are desired for problems which lack particular symmetries, cover many length-scales or require adaptivity, for example stellar collisions star or star cluster formations.

Smoothed particle hydrodynamics (SPH) proved to be a successful Lagrangian meshless scheme to solve equations of fluid dynamics in variety research fields (Monaghan, 2005). Despite its limitations, it has been also successfully used in wide range of astrophysical problems, including evolution of gaseous disks around black holes or stars, star formation, stellar collisions, and cosmology. Because of its simplicity and versatility, several attempts, albeit with limited success, have been made to include magnetic fields into SPH, thereby formulating smoothed particle MHD, or SPMHD for short (Price & Monaghan, 2005; Børve et al., 2006; Rosswog & Price, 2007). It soon became clear that SPMHD equations where plagued with two main problems: tensile instability and maintenance of the divergence constraint, . The former is a general problem of SPH, namely the equations are unstable to tensile stresses (Swegle et al., 1995; Monaghan, 2000). In case of SPMHD, this instability manifest itself as particle clumping in the regions where magnetic pressure dominates gas pressure, and therefore rendering the simulation of strongly magnetised plasma unfeasible. It has been shown that the equations can be stabilised by either adding short-range repulsive forces between particles (e.g. Price & Monaghan, 2004), or sacrificing momentum conservation (e.g. Price & Monaghan, 2005). Alternatively, Børve et al. (2001) showed that the stability of SPMHD equations can also be achieved by adding a source term proportional to the divergence of magnetic field to the momentum equation. While these approaches appear to remove the tensile instability, the divergence constraint still remains an issue in these SPMHD formulations. Rosswog & Price (2007) were able to formulate manifestly divergence-free SPMHD equations which appears to work for a wide range of problems (e.g. Price & Rosswog, 2006, Price & Bate, 2008). However, they solve a limited form of the induction equation which only permits topologically trivial field configurations and is unable to model more complicated MHD phenomena, such as magneto-rotational instability. Dolag & Stasyszyn (2009) were also able to formulate stable SPMHD equations by using a combination of several techniques, such as the addition of the source term proportional to the magnetic divergence to the momentum equation, artificial dissipation and smoothing of the magnetic field, and modification of the induction equation. This approach introduces several free parameters, which the authors were able to constrain by fitting their results to the solutions of several shock tube problems computed with conservative Eulerian MHD schemes. Despite the progress in this field, the difficulties associated with formulating consistent SPMHD equations advocate the need of the alternative approaches to formulate meshless Lagrangian MHD schemes. In a gradient particle magnetohydrodynamics (GPMHD, Maron & Howes, 2003; Maron, 2005), the equations of ideal MHD were discretised on a set of particles by fitting a second or forth order polynomial into the data in order to obtain first and second order derivatives of the desired quantities. However, such approach still requires addition of artificial diffusion to the resulting equations to model dissipative and resistive processes across discontinuities, and this introduces additional free parameters which control dissipative processes. In addition, the local conservation, for example third Newton law, is satisfied to the truncation error only. This is a potential source of error for the problems with interacting strong shock waves, since the truncation error across a shock wave is always of order of unity. Nevertheless, GPMHD appears to be a viable, albeit noticeably more complex, alternative for SPMHD, at least for subsonic and weakly supersonic flows.

Alternatively, one may attempt to utilise a Godunov approach. It has been shown that most of, if not all, SPH limitation can be removed by solving a hydrodynamic Riemann problem between each pair of interacting particles to obtain pressure forces, instead of computing separately pressure forces and artificial viscosity terms (Inutsuka, 2002; Cha & Whitworth, 2003; Cha et al., 2010). In such Godunov SPH (GSPH) formulations, the necessary mixing and dissipation is included into underlying Riemann problem which described the physics of the interaction. Borrowing these ideas, it is therefore conceivable that Godunov-like meshless MHD formulation will eliminate the problems that plague SPMHD formulation, thereby permitting formulation of consistent Lagrangian meshless MHD scheme.

In this paper we formulate a weighted particle MHD scheme. Our scheme is based on a meshless discretisation of the conservation law equations, which was pioneered by Vila (1999), and in Section 2 we present a heuristic derivation of the meshless conservative equations. We give the implementation details of our scheme in Section 3. In Section 4.1 we present applications to the equation of ideal hydrodynamics, and in Section 4.2 we show how our scheme can be applied to the equations of ideal MHD. In Section 5 we validate our meshless MHD scheme on several test problems, and finally, we present our conclusions in Section 6.

2 Methods

2.1 Meshless equations for conservation laws

In what follows, we present a heuristic derivation of meshless discretisation of a scalar conservation law. The readers interested in a rigorous mathematical formulation supplemented with convergence theorems are referred to the original papers by Vila (1999) and Lanson & Vila (2008a, b). Following these works, a weak solution to a scalar conservation law


is defined by


Here, is a scalar field, is its source, is its flux in a frame moving with velocity , and the integral is carried out over all space-time domain of a problem at hand. The function is an arbitrary differentiable function in space and time, is an advective derivative, and is an arbitrary smooth velocity field which will describe motion of particles. This integral is discretised on a set of particles with coordinates with the help of a partition of unity


where, is an estimate of the particle number density and the sum is carried out over all particles, is a smoothing length, and is a smoothing kernel111In this paper we use cubic-spline smoothing kernel, which is commonly used in SPH. with a compact support of size ; in what follows we assume that the kernel is normalised to unity. Inserting into an integral of an arbitrary function, we obtain


where is the effective volume of a particle , and in the third term we use first-order Taylor expansion of . In principle, a higher order discretisation is also possible, but for the purpose of this work such an one-point quadrature is sufficient; in fact, on a regular distribution of particles this discretisation is second order accurate (c.f. §3.2). Application of this discretisation to Eq. 2 gives


Here, the Einstein summation is assumed over Greek indexes, which refer to the components of a vector. Also, the gradient of a function at -particle location, is replaced by its discrete version , and this, for example, can be computed with SPH estimates of a gradient. However, a much better approach is to employ a more accurate meshless gradient estimate suggested by Lanson & Vila (2008a). They showed that a second order accurate meshless partial derivative is given by the following expression


where , and is a renormalisation matrix defined by its inverse


Finally, integrating the first term by parts with assumption that vanishes at boundaries, and applying the following rearrangement to the second term


permits separation of from the rest


The above is true for an arbitrary function if the expression in brackets vanishes, namely


The extension of this equation to a general vector field is straightforward: this equation is applied to each component of the field. These equations are similar to SPH equations, with the difference that the physics of the particle interaction is hidden in the fluxes and source terms. Indeed, if one uses the fluxes of ideal Lagrangian hydrodynamics, the equation of motions similar to SPH can be derived. As with SPH, such equations do not include dissipative processes, and therefore must be augmented with explicit diffusive terms, namely in the form of artificial viscosity, conductivity and resistivity.

However, the power of the new scheme becomes apparent with the realisation that one can utilise the fluxes produced by the solution of an appropriate Riemann problem between particles and , which automatically include necessary dissipation. Defining such a flux as , and setting gives


Finally, defining a vector and , the equations take the following form


It becomes clear, only the projection of the flux on the direction of the vector is required, and therefore for a wide range of problems the flux can be obtained by solving an appropriate 1D Riemann problem in a frame moving with mean velocity of the two particles, .

2.2 Linear monotonic reconstruction

The one dimensional flux in Eq. 12 is naturally computed at the midpoint between particles and , i.e. at . To achieve higher than the first order accuracy, it is necessary to linearly reconstruct left and right states of the Riemann problem to this location. The reconstruction step should in principle be done in characteristic variables, however for the second and third order schemes this can be done in primitive variables, , as well. In the case of MHD, the latter are density , pressure , velocity and magnetic field . An approximation of of an -paritcle state at is given by a first-order Taylor expansion from the point :


where is a gradient estimate of the primitive variables computed with Eq. 6, and is a vector of limiting functions which is required in order to assure non-oscillatory reconstruction (Balsara, 2004)


Here, and are the minimal and maximal respectively over all neighbours that particles interacts with, and and are the minimal and maximal resulted from the reconstruction in Eq. 13 for each of these neighbours. The scalar constant vector should have values between and in order to achieve second order of accuracy. Following suggestions of Balsara (2004), value should be used for both pressure and velocity, and for both density and magnetic field; however, we find no problems while using for all fluid quantities. The frame velocity at is approximated as . Finally, the reconstructed left and right states, the frame velocity , and the unit vector are used to obtain the flux from the 1D Riemann problem.

Instead of a linear, a piecewise parabolic reconstruction can also be used to achieve third-order spatial accuracy. In Appendix A, we describe parabolic reconstruction of a scalar field . Due to large operation count, this reconstruction is presented only for purpose of completeness, and in the test that will follow later, only linear reconstruction is used.

3 Implementation

3.1 Smoothing length

The smoothing length in our scheme is a property of the particle distribution and, in contrast to SPH, does not depend on the fluid state. In principle, a constant smoothing length can be used throughout the whole space and time domain. In practice however this cause difficulties due to the possible development of wide range in particle number densities as the simulation progresses. Similarly to SPH, this can lead to under- or oversampling in low and high particle density regions respectively. The approach used here is inspired by conservative SPH formulations (Monaghan, 2002; Springel & Hernquist, 2002). The idea is to constrain the smoothing length of a particle , , to the particle number density at this location, , i.e.


This tends to maintain approximately number of neighbours for each particle; here , and for , and dimensions respectively, and , where is defined in Eq. 3. As in conservative SPH equations, is obtained by iteratively solving Eq. 15, for example via Newton-Raphson method (e.g. Press et al., 1992). One might be also tempted to use the continuity equation


to compute time evolution of the smoothing length from its initial value. This, however, is undesirable for two main reasons: a) the result depends on the functional form of the divergence operator, and b) in discontinuous flows the may be undefined at some points, which can result in unexpected behaviour. As a result, in our tests we chose to iteratively solve Eq. 15, but we use the differential form to predict as a first guess to an iterative solver.

Finally, knowledge of smoothing length permits calculation of the rest of geometric quantities, such as effective volume of a particle, . It is possible to use numerical quadrature to evaluate with a desired accuracy, however we find that defining works fine for our purpose, and therefore we decided not to perform more accurate volume estimates.

3.2 Particle regularity

Particle regularity is an important aspect of the scheme. If particles are randomly sampled within a domain, there is non-zero probability that particle’s smoothing length, , will differ significantly from the average in its neighbourhood. Furthermore, the resulting -distribution will not be a smooth function of position, and therefore will not be differentiable. This will break the approximation which lead to Eq. 10. Namely, the variation of within the neighbour sphere will be large enough that the estimate in Eq. 4 will result in intolerable errors which produces unexpected behaviour, such as negative values of density or pressure. To avoid these situations, the particle distribution must be first regularised. If the initial particle distribution is regular, it will maintain its regularity during the simulation except in the regions where particle velocity field, , is discontinuous, e.g. across shock waves (Vila, 1999; Lanson & Vila, 2008a). The criteria which determines regularity of the particle distribution depends on the approximations of Eq. 4. Expanding to the first order, gives


The first term is , and we can rewrite the integral in the second term in the following form


where in the right hand side we changed variables from to . If we require that is approximately constant in the neighbourhood of an -particle, we can discretise the integral on the right hand side to obtain . Hence, if the particles are relaxed such that this sum is minimised, the start up noise becomes negligible. Empirically, we found that the particle distribution is regularised if the following quantity is minimised




Ideally, should be equal to zero, but this is appears to be only possible if particles are arranged on a lattice, for example a cubic or hexagonal close-packed lattice. If particles are sampled randomly, which is more desirable in many problems, their positions must be adjusted until Eq. 19 reaches its (approximate) minimum before assigning fluid state to the particles; this will reduce the start-up noise in a simulation. Afterwards, this particle distribution can be used to assign initial conditions for a problem at hand. To regularise, or to relax, particle distribution, we use the following iterative procedure. First, we compute Eq. 20 for all particles. Afterwards, -particle position is updated: , where ; such update reduces . This operation is repeated until is reached its minimum, or a desired minimal value.

3.3 Time marching

Due to Godunov’s nature of the particle conservation laws, it is tempting to employ Hancock scheme (e.g. van Leer, 2006) to achieve a single-stage second-order accurate time integration. However, such scheme does not include a predictor for transverse waves in multi-dimensional Riemann problem, and therefore these remain only first-order accurate. In multiple dimensions, unsplit numerical schemes usually adopt Corner-Transport-Upwind method (CTU, Colella, 1990; Stone et al., 2008) for one-step second-order time integration, but applicability of CTU to meshless schemes is not clear. Nevertheless, higher than the first order accurate time integration can be achieved with multi-stage total-variation diminishing (TVD) Runge-Kutta methods (Gottlieb & Shu, 1998). Here, we use a two-stage second order TVD Runge-Kutta time marching scheme


where and are time derivatives of an -particle computed from and respectively via Eq. 12, is the effective volume of -particle, and , where , and for 1D, 2D and 3D respectively is a measure of particle linear size, is particle’s signal speed (speed of sound for HD and of the fast magnetosonic wave for MHD), and is a usual Courant-Fridrisch-Levy number.

Particle positions obey the following equation of motion


where is particle’s velocity. In practice we set it equal to the fluid velocity, , and the time integration is carried out in drift-kick-drift approach. Namely, the particles are first drifted from their current positions, , to the position at half time-step


This particle distribution is used to compute smoothing lengths via Eq. 15 and other geometric quantities. Afterwards, we apply a two stage TVD Runge-Kutta method to perform an update from to while keeping particles fixed in space and setting , where is -particle volume at half time-step. Finally, the particles are drifted for another half time-step with the updated fluid velocity, ,


3.4 Non-conservative formulation

In the case of an HD or MHD system, the conservative formulation updates total energy instead of thermal; thermal energy is obtained by subtracting magnetic and kinetic energies from the total energy. When the supersonic advection is present, the sum of thermal and magnetic energies, , is obtained by subtracting two large numbers, namely total energy and kinetic energy. To avoid this, we suggest an alternative non-conservative formulation, which evolves instead of the total energy. Writing , where and are momentum and mass respectively, gives


The latter term can be rewritten as resulting in the following equations for


Here, and are time derivatives computed from conservative meshless equations. In this form, the total energy will not be conserved to the machine accuracy, but rather to the truncation error of the time-marching scheme. In other words, the integration error is now lost from the system instead of appearing in the thermal energy. Furthermore, the total energy can now be used as a quality control indicator of a simulation.

3.5 Modification of existing SPH codes

The existing SPH codes which use conservative SPH formulation (Monaghan, 2002; Springel & Hernquist, 2002) can be straightforwardly modified to implement our weighted particle scheme due to similarity of SPH equations of motions and our equations of meshless conservation laws, Eq. 12. In particular, the neighbour search should be modified such that Eq. 15 is solved, which constrains the number of particles, instead of the enclosed mass, in the neighbour sphere. At the end of this step, the for each particle will be known that permits to compute volume of a particle, . This volume is required to convert conservative fluid variables, , to primitive ones . Afterwards, the first loop is carried out over gather neighbours which computes the renormalisation matrix, Eq. 7, and gradients of primitive fluid variables. The second neighbour loop is required to compute the limiting functions, Eq. 14, and this loop must be carried out over both gather and scatter neighbours due to need to limit reconstruction to each of the interacting particles. Finally, in the third neighbour loop the interactions between the particles are computed. This is done in exactly the same way as in SPH, except that for every -neighbour of an -particle, the fluid states are reconstructed at the midpoint, through Eq. 13. Finally, these together with the vector in Eq. 12 and the interface velocity , are used in the Riemann solver to compute the interface fluxes, , which we describe in the following sections. At the end of this final neighbour loop, one will have time derivative that should be used to update the conservative fluid variables. The global conservation of mass and other conservative quantities, in the absence of source terms, is maintained to the machine precision independently of the particle distribution. This also includes total energy, unless Eq. 27 is used, in which case the total energy is conserved to the truncation error of a time integration scheme.

4 Fluid dynamics

Our meshless conservative equations can be applied to a system which can be written in the form of conservation laws. The physics in this case is completely described by the source terms, , and fluxes, . The latter can be obtained by solving an appropriate one-dimensional Riemann problem between a particle and its neighbours. Hence, the Eq. 12 can be applied to a variety of problems, such as hydrodynamics, magnetohydrodynamics, and radiative transfer in the flux-limited diffusion approximation. In what follows, we present application of our scheme to both ideal hydrodynamics and magnetohydrodynamics.

4.1 Ideal hydrodynamics

The Euler equations of ideal hydrodynamics in a frame moving with the velocity read




where is a unit tensor, and other symbols have their usual meaning.

To obtain fluxes, an 1D Riemann problem is solved between two particles in the direction. This is accomplished by defining the rotation matrix, , such that in the new coordinate system coincides with the -axis, i.e. . This transformation is applied to all vector quantities from both the left and the right states, the scalar quantities are left untouched. For example, the velocity transformation results in , where for the left and for the right state. These transformed states are the initial conditions of the Riemann problem for 1D Euler equations


where, and


The space discretisation of this equation has the following form


where is size of the grid cell in 1D, and is an interface flux between cells and . This is the flux required to substitute into Eq. 12, and it can be obtained from the solution of an appropriate 1D Riemann problem.

The solution of this 1D Riemann problem gives the fluxes, , for each of the component of in direction . While the fluxes of the scalar components can be directly used in Eq. 12, those of a spatial vector, however, must be rotated back to the original coordinate system because the flux of an -component of the vector, rather than , is required in the direction . In the case of the velocity vectory, one first combines -flux into a vector . Then an inverse transformation is applied to compute -flux, namely , and these fluxes are then substituted into Eq. 12.

This approach permits the use of Riemann solvers that directly approximate interface flux, rather than fluid states. In what follows, the HLLC Riemann solver is used to compute the flux. The HLLC solver requires only velocity estimates for certain characteristic waves, and is oblivious to the exact form of the equation of state. As a result it can be easily extended to equations of states other than that of ideal gas. Here, we present formulae of the HLLC Riemann solver in moving frame which can be directly used in Eq. 12; for the detailed derivation we refer the interested reader to published literature (e.g. Toro, 1999 (§10), Miyoshi & Kusano, 2005).

Figure 1: This figure shows the wave-structure of HLLC Riemann solver. Namely, HLLC solver resolves all three HD waves: left, right, and middle wave also known as contact discontinuity waves. The left and right wave can be either shock or rarefaction wave, or both. This solver approximates the structure of the two intermediate states, between and , and and , by a constant state.

A typical wave-structure of HLLC solver is shown in Fig. 1. The velocity of left, right and middle waves are , and respectively, and the constant states sandwiched between these waves are and respectively. The HLLC-flux at the interface moving with velocity is


where and . The intermediate states are defined by


for and . Finally, estimates of wave speeds are and , where . Other wave-speed estimates, such as those based on Roe-averages, can be used as well. Finally, the speed of the middle wave and the pressure in the -states is given by the following formulae


with and .

4.2 Ideal MHD

Our meshless conservative equations can also be applied to the equations of ideal MHD. In particular, we can rewrite MHD equations in the conservative form with the following conservative variables, fluxes and source terms


here is sum of the thermal and magnetic pressures. The biggest difficulty in solving these equations is to maintain constraint. While there are several ways to satisfy this constraint in finite-difference methods, it is not clear how this can be done in meshless schemes. In SPH, Rosswog & Price (2007) were able to circumvent this problem via use of the Euler potentials, and , which in smooth flows are advected with the flow, and are used to compute magnetic field . While mathematically this guarantees zero-divergence, in practice, however, the divergence is non-zero because of non-commuting nature of SPH cross-derivatives. This is also holds for the vector potential, , which may explain unstable behaviour of SPMHD equations with vector potential (Price, 2010). While SPH equations appear to be stable in the Euler potential formulation, they are unable to model topologically non-trivial field configurations (Brandenburg, 2010), and therefore are not suited to simulate complex magnetic field evolution, such as winding of magnetic field lines and magneto-rotational instability. Motivated by this, and by the fact that the existing MHD Riemann solvers use magnetic field as a primary quantity, we chose to evolve instead, and apply one of the divergence cleaning methods known to work in Godunov finite-difference MHD schemes.

In their paper, Powell et al. (1999) showed that a self-consistent MHD equations must include source terms proportional to to insure stability and Galilean invariance


While this formulation both removes the instabilities associated with non-zero divergence and tends to keep the divergence at the truncation level, the divergence can still grow in certain situations. To avoid this, we also include a Galilean invariant form of the hyperbolic-parabolic divergence cleaning method due to Dedner et al. (2002), which results in the following system


Here, and is the speed and the damping timescale of the divergence wave, also known as 8-wave, respectively. Usually, is set to be highest characteristic speed, i.e. that of the fast magnetosonic wave, , where is an effective size of the particle defined in §3.3, and is a constant, which, following Mignone & Tzeferacos (2010), we set .

In their paper, Dedner et al. (2002) also presented a method to generalise an arbitrary 1D MHD Riemann solver to include an additional scalar field which transports the divergence away form the source and damps it. In general, the left and right reconstructed states in coordinate system (see §4.1) have discontinuous –the magnetic field component normal to the interface. The 1D Riemann solvers, however, require this field to be continuous, and this can be computed with the following equations


Here, are the left () and right () reconstructed states for and , and . These interface values of and are used to compute and respectively, which are required for the source terms in Eq. 39,


In Appendix B, we provide standard formulae to compute HLL and HLLD fluxes for the 1D MHD Riemann problem, which take as the continuous normal component of the magnetic field. Finally, the advection flux for the scalar is


where, is a mass flux given by the Riemann solver. This completes the description of our meshless MHD scheme.

5 Scheme validation

The weighted particle scheme is validated on several standard problems for ideal hydrodynamics and MHD problems, which are usually used to test hydrodynamic and MHD schemes (e.g. Tóth, 2000, Stone et al., 2008. In all simulations what follows, we use and in Eq. 15 for 2D and 3D simulations respectively. This choice is motivated by the analogy with finite-difference schemes. In one dimension, the number of the neighbouring cells that a given cell interacts with depends on the order of the numerical scheme, and is usually two for the second order scheme. In SPH, two to four neighbouring particles are usually used in one-dimensional simulations. If this number is scaled to three dimensions, and taking into account that a kernel has spherical shape, the estimated number of neighbours are 16 and 32 for 2D and 3D respectively. The MHD simulations in 2D appear noisy with , which motivated us to use larger . In principle, large can be used, should this be necessary, but this will decrease the resolution due to larger smoothing length.

5.1 Shock tubes

5.1.1 Brio-Wu shock tube

This problem was introduced by Brio & Wu (1988) to test the ability of an MHD scheme to accurately model shock waves, contact discontinuities and compound structures of MHD. Here, the problem is solved in a periodic 2D domain of size with randomly sampled particles; this results in an effective resolution of . The particles are initially relaxed before the initial conditions are set. The left state, , is set with the following values: , , . The right state has , and . Both states have zero initial velocities, and . The problem is solved with an ideal gas equation of states and . Various profiles, e.g. density and pressure, at time are show in Fig. 2.

Figure 2: Solution to the Brio-Wu MHD shock tube problem at time . The panels display (from left to right, top to bottom) density, gas pressure, - and - velocity profiles, , thermal energy , and the logarithm of , where is an effective size of a particle (§3.3).

The results are overall consistent with 1D Godunov-MHD scheme, namely jumps across the discontinuities and location of discontinuities. The two bottom right panels show the value of field and the as a function of particle -coordinate. The parallel magnetic field slightly deviates from its constant value, except near discontinuities where it exhibits jumps. The divergence, however, remains small, even across discontinuities. The existence of blip in pressure at location of contact discontinuity, , and shock waves, has the same origin as in SPH: the particle distribution across a discontinuity is less regular in a sense that the approximation in Eq. (4 is not sufficient to provide accurate results. This is a known issue in Godunov SPH, and higher order approximations are able to reduce the amplitude of the blip (Inutsuka, 2002). Overall, the solution obtained by the meshless scheme is in a good agreement with high-resolution 1D Euleriean schemes.

5.1.2 Tóth shock tube

Another challenging shock tube problem was introduced by Tóth (2000). In this problem, two streams of magnetised gas supersonically collide with each other. Tóth (2000) showed that some MHD schemes with source terms proportional to produce wrong jump conditions. This was challenged by Mignone & Tzeferacos (2010), who showed that if divergence is cleaned in hyperbolic-parabolic manner (Dedner et al., 2002), the jump conditions are correct even with the presence of the source terms. The problem is set in 2D periodic domain with size with randomly distributed particles. After the particle distribution is relaxed, the following initial conditions are set. The left states has , , , and the right state has , , . Both states have . This problem is solved with an ideal gas equation of state and . The solution to this problem is shown in Fig. 3, which can be compared to solutions obtained by (Tóth, 2000) and Mignone & Tzeferacos (2010). Our meshless scheme is able to recover correct jump conditions and to maintain constant field within few percent accuracy, except across discontinuities. Furthermore, the divergence in this problem remains small.

Figure 3: Solution to the Tóth MHD shock tube problem at time . The panels display (from left to right, top to bottom) density, gas pressure, - and - velocity profiles, , thermal energy , and , where is an effective size of a particle.

5.2 Advection of a magnetic field loop

This problem tests the ability of the scheme to transport magnetic loop across computational domain. This problem is proven to be a stringent test for finite-difference schemes. The computational domain in this test is a cuboid with dimensions which is initially filled with particles on cubic lattice to assure zero noise. The fluid state is set to inside the loop and outside the loop, , and , where and


Here, is a radius of the loop, is the initial magnetic field strength that results in . With such high magnetic field does not play dynamical role and should be transported as a passive scalar. Periodic boundary conditions are used in this test, and the ideal gas equation of state with .

In Fig. 4 we shows magnetic field structure at the and , which corresponds to ten crossings of the computational domain, and in Fig. 5 we show the magnetic energy as a function of time. This figures demonstrate the ability of the meshless scheme to advect magnetic loop quite well. Furthermore, the decay of the magnetic energy is at least as slow as resulted from high-order Godunov MHD schemes. This is certainly expected in light of semi-Lagrangian nature of the scheme. One may however expect that the energy should not decay at all since the magnetic field is transported as a passive scalar. Indeed, Fig. 6 demonstrates that the scheme is able to transport mass, and therefore passive scalar, without any diffusion. The difference with magnetic field stems from the different nature of equation that transport mass, or advect passive scalar, and the induction equation as implemented in the scheme. In fact, the decay is caused solely by the diffusion of magnetic field due to divergence cleaning equations, i.e. Eq. (40 and Eq. 41. The fluxes resulting from HLLD Riemann solver are in fact zero since both particles and the frame move with the same velocity. The divergence of the magnetic field, even though is zero analytically, is not necessary zero in the discretisation set by Eq. 40 andEq. 42, and this produces evolution of magnetic field due to non zero value of in the source terms of Eq. 39. Among all discretisations studied by Tóth (2000), this is the special one because of its use in the discretisation of Maxwell stress term to obtain Lorentz force. If the in this discretisation, no force parallel to magnetic field exists, and therefore the source terms proportional to this divergence vanishes. Incidentally, this is the discretisation of divergence that is enforced to zero in constrained transport (CT) formalism (Evans & Hawley, 1988). More importantly, the maintenance of zero divergence in other discretisations, such as cell-centred, does not guarantee vanishing divergence in CT-discretisation, but this will probably be small thus giving minimal damage to the solution.

Figure 4: Amplitude (colours) and direction (vectors) of magnetic field at the beginning (left panel, ) and at the end (right panel, ) of the simulation. The magnetic field distribution is show in the -plane passing through . In the right panel, the strength of magnetic field in the centre is close to zero due to numerical resistivity, which is consistent with finite-difference calculations.
Figure 5: This figure demonstrates decay of magnetic energy during advection of magnetic field loop. The magnetic energy decay is due to numerical resistivity introduced by divergence cleaning procedure, and not through the solution of an MHD Riemann problem.
Figure 6: This figure shows that our scheme is able to transport mass, and therefore passive scalars, without any numerical diffusion in a constant velocity field. The right and left panels show plots of the density field in the -plan passing through , at the beginning and the end of the simulation respectively.

5.3 Blob test

An interesting and challenging problem to test particle hydrodynamic schemes has been proposed by Agertz et al. (2007). This problem demonstrates the destructive property of the Kelvin-Helmholtz instability (KHI) in three-dimensional simulations. The setup consist of a dense cloud with moving supersonically through a less dense ambient fluid with . Initially, the cloud is at rest and in the pressure equilibrium with the ambient fluid, at . The velocity of the ambient medium is , where is its sound speed. The cloud radius is , and an ideal gas equation of state is used with . The initial magnetic field is set to zero, and it will stay so throughout the simulation. This problem is set in a periodic domain with dimensions , in which particles are sampled in the strip and outside. This setup permits to resolve cloud and the impacting ambient fluid with high resolution () for a total of particles. The particles are sampled from a three-dimensional Sobol quasi-random sequence (Press et al., 1992). To remove the initial noise, this initial particle distribution is relaxed before the initial values for the density field, pressure and velocity are set.

Following Agertz et al. (2007), the Kelvin-Helmholtz timescale is defined , where . For parameters used in this simulations, this gives . The snapshots of the density distribution in plane at and are shown in Fig. 7. These are in an excellent agreement with those presented in Agertz et al. (2007). The time dependence of cloud mass is shown in Fig. 8, where a particle is considered to be part of a cloud if and . In agreement with finite-difference methods, the cloud lost nearly 90% of its mass within .

Figure 7: Density field plot in -plane passing through . The panels show density profile at time , , and , where . At (top panel), the cloud shape is distorted by RAM pressure. At (second panel from the top), the KHI instability deforms the shape of the cloud further. The cloud destruction begins at , and by the cloud lost more than 90% of its mass.
Figure 8: This plot show cloud mass as a function of time. In agreement with finite-difference scheme the cloud looses most of its mass at , and is completely destroyed at .

5.4 Spherical blast-wave

The problem is initiated with an overpressured central region in a uniform density and magnetic field. The computational domain is a unit square filled with fluid with . Within , the pressure is set to , whereas outside . In magnetised case, there is also a uniform magnetic field . The equation of state is that of an ideal gas with . The particles are sampled in a periodic box , such that particles are randomly sampled in three nested rectangles: , and . Before the initial conditions are set, the particles are relaxed into a regular distribution to reduce start-up noise. In figures 9 and 10 we show density plots for non-magnetised and magnetised cases respectively. Of particular interest here is the ability of the particle weighted method to resolve Richtmyer-Meshkov instability, shown in right panel of Fig. 9. In the magnetised case, however, the presence of strong magnetic field inhibits development of this instability (the right panel of Fig. 10).

Figure 9: This figure show the density at and for non-magnetised spherical blast wave problem. The location of the shock front is in agreement with high-resolution conservative numerical schemes (Stone et al., 2008). Right panel show dense fingers in rarefied media which are formed by Richtmyer-Meshkov instability. This demonstrates that our scheme is able to capture important fluid instabilities without any fine-tuning.
Figure 10: This figure show the density at and for magnetised spherical blast wave problem. The initial magnetic field in this problem has angle of with -axis. Since the shock moves more easily along the magnetic field lines, this explains the elongated shape of the shock-front, which contrasts with hydrodynamical case. As in hydrodynamical test, here the location of shock from is in excellent agreements with conservative Eulerian MHD schemes (Stone et al., 2008). The right panel, shows further evolution of the black wave, after shock reached the rarefied medium. In contrast to the hydrodynamical case, magnetic field inhibits development of Richtmyer-Meshkov instability.

5.5 Orszag-Tang vortex

The Orszag-Tang vortex (Orszag & Tang, 1979) is a standard test problem that is used to validate many numerical MHD schemes. The setup involves periodic domain of size with an adiabatic equation of state with . The initial density and pressure are set in all computational domain to and respectively. The velocity and magnetic field , where . The second simulation involves the same initial conditions, with the exception that the problem is solved a boosted frame, with the initial velocity , where . The simulation is solved with particles, where first where randomly sampled within the square and the second in the square of half size, . This permitted to achieve high resolution in the central region of the problem. The initial conditions were set, as soon as distribution was relaxed.

In the Fig. 11 we show density at for both rest-frame and boost-frame initial conditions, and in Fig. 12 we show particle distribution. As expected, there are no discernible differences, and both simulations resolve discontinuities well. Furthermore, due to Lagrangian nature of the method, the particle number density correlates with the mass density. It is worth pointing out that in this particular simulation, the particle distribution appears to remain regular even in the vicinity of the shocks. The Fig. 13 show 1D pressure profile at (top) and (bottom) for both boosted and rest-frame simulations. The agreement with finite-difference scheme is excellent, ac can be compared to published results (e.g. Tóth, 2000, Rosswog & Price, 2007, Stone et al., 2008). In contrast to shock-tube problems, no pressure blips are visible here.

Figure 11: Density distribution of Orszag-Tang vortex at . The left and right panels show the result of the simulation in the rest and boosted frame respectively. We point out that even at low resolution our scheme is able to resolve the fine structure in the bottom left quadrant; in this quadrant, the effective resolution is , compared to in the centre.
Figure 12: Particle distribution of Orszag-Tang vortex at . The left and right panels show the particle distribution in all domain and in the central region respectively. Due to Lagrangian nature of the simulation, the particle distribution correlated with that of the density. We would like to note that even across the shocks, the particle distribution remains remarkably regular without any observable signs of clumping.
Figure 13: Pressure distribution of Orszag-Tang vortex at and (top) and (bottom). This can be compare to the results of other finite-difference or meshless MHD schemes (e.g. Rosswog & Price, 2007; Stone et al., 2008). In particular, our scheme has little noise in pressure, and the shock waves are resolve within few particles, without observable oscillatory behaviour in post-shock regions.

5.6 MHD rotor

The rotor problem, introduced by Balsara & Spicer (1999) to test propagation of strong torsional Alfvén waves, is also considered as one of the standard candles to validate numerical MHD schemes. Here, the computational domain is a unit square, . The initial pressure and magnetic field are uniform with values and . Inside there is a dense uniformly rotating disk with and , where . The ring is occupied by a transition region with and , where


Outside , the velocity is set to zero and . This problem uses an ideal gas equation of state with . The particles are distributed in the same way as in strong blast wave problem (§5.4). The problem is solved in both rest and boosted frames, with . The plots of density, gas and magnetic pressure, and mach number are shown in Fig. 14. The 1D slices of magnetic field at and are show in Fig. 15 to facilitate comparison to Eulerian MHD schemes.

Figure 14: Density (top left), gas pressure (bottom left), magnetic pressure (top right) and Mach number (bottom right) plots for MHD rotor problem at . The Mach number panel plots Mach number using rest-frame velocities, that is first subtracting the boost velocity and afterwards computing mach number. The result of the rest frame simulation have little difference, and therefore not presented here. The colour coding is chosen to facilitate the comparison with published results produced by FLASH3 code for this problem (Fryxell et al., 2000, Sect 21.2.3 in FLASH 3 user guide manual
Figure 15: Distribution of in slice (top pannel) and in slice (bottom) for MHD rotor problem at . This plots can be readily compared to Fig.26 of (Stone et al., 2008). The agreement is rather good everywhere, except for near and , where the meshless scheme appears to diffuse the discontinuity in magnetic field.

5.7 Magneto-rotational instability

Magneto-rotational instability, or MRI for short, is a powerful local instability in weakly magnetised disks that plays an important role in astrophysics (Balbus & Hawley, 1991, 1998). The ability of the scheme to model this instability is crucial for simulation of magnetised accretion disks, or any other simulation where non-trivial MHD effects are expected to appear. Below we conduct two simulations that test the ability of the meshless MHD scheme to faithfully model MRI.

5.7.1 2D axisymmetric shearing box

The first simulation is solved in a local 2D axisymmetric shearing sheet with the initial conditions exactly the same as in the fiducial model of Guan & Gammie (2008). To repeat, the unit box has a size with an initially uniform density fluid with unit density, , embedded in a vertical magnetic field