rpSPH: a novel Smoothed Particle Hydrodynamics Algorithm
Abstract
We suggest a novel discretisation of the momentum equation for Smoothed Particle Hydrodynamics (SPH) and show that it significantly improves the accuracy of the obtained solutions. Our new formulation which we refer to as relative pressure SPH, rpSPH, evaluates the pressure force in respect to the local pressure. It respects Newtons first law of motion and applies forces to particles only when there is a net force acting upon them. This is in contrast to standard SPH which explicitly uses Newtons third law of motion continuously applying equal but opposite forces between particles. rpSPH does not show the unphysical particle noise, the clumping or banding instability, unphysical surface tension, and unphysical scattering of different mass particles found for standard SPH. At the same time it uses fewer computational operations. and only changes a single line in existing SPH codes. We demonstrate its performance on isobaric uniform density distributions, uniform density shearing flows, the Kelvin–Helmholtz and Rayleigh–Taylor instabilities, the Sod shock tube, the Sedov–Taylor blast wave and a cosmological integration of the Santa Barbara galaxy cluster formation test. rpSPH is an improvement these cases. The improvements come at the cost of giving up exact momentum conservation of the scheme. Consequently one can also obtain unphysical solutions particularly at low resolutions.
keywords:
Numerical Methods–Smoothed Particle Hydrodynamics–Hydrodynamics–Instabilities1 Motivation
The smoothed particle hydrodynamics method was invented by Lucy (1977) and Gingold and Monaghan (1977), both with interests in astrophysical applications. Besides an enormous literature of successful application also many shortcomings of it have been presented in the literature (Steinmetz and Mueller, 1993; Herant, 1994; Swegle et al., 1995; Imaeda and Inutsuka, 2002; Agertz et al., 2007; Read et al., 2009) often suggesting a fix to the reported problem (Monaghan, 1994; Cummins and Rudman, 1999; Rasio, 2000; Attwood et al., 2007; Hu and Adams, 2007; Graham and Hughes, 2008; Price, 2008; Børve et al., 2009; Rafiee and Thiagarajan, 2009; Xu et al., 2009, to name but a few) Similarly there have been troubling news of how seemingly small differences in the initial setup led to very unexpected results (e.g. Lombardi et al., 1999). This is all somewhat surprising given that in many of the cases where large inaccuracies have been found the only relevant equation (besides moving the particles ) stems from the pressure gradient accelerations
(1) 
where denotes the Lagrangian derivative, and the pressure. In what follows we describe a new discretization of the momentum equation that avoids essentially all of the previously known problems of SPH. We will refer to this new method as “relative pressure SPH” or abbreviated as rpSPH. We will first describe it, discuss implementation details and then present results for relevant tests highlighting the superior performance of this new approach.
All the simulations shown here are carried out with Gadget2 (Springel, 2005) (version 2.0.4) with only most minor changes explained in the text. The appendix describes how to convert Gadget to our rpSPH formalism.
2 Relative Pressure SPH
The equation of motion without viscous or gravitational forces in essentially all SPH codes and Gadget2 (Springel and Hernquist, 2002; Springel, 2005), in particular, is
(2) 
where the are defined by
(3) 
and the abbreviation has been used for the kernel function, . It employs variable smoothing lengths so that the number of neighbours for each particle with is maintained at a nearly fixed value . The compact cubic spline kernel is used which is summarized Monaghan (1992), extends to radii as large as the smoothing length and is zero outside. While many choices would exist to use a different discretisation here (Monaghan, 1992; Rosswog, 2009) most previous work we have found essentially retains a form very close to equation (2) or uses another symmetric form that sums over . The reason previously given for these choices are their symmetric form encapsulating Newtons third law of motion, the actionreaction law. By guaranteeing that particles give pairwise identical but reversed forces one ensures linear momentum conservation of the entire scheme. The key here is that particles are always pushing as soon as they have any pressure regardless of whether there is a pressure gradient. The ones with the highest pressure values are pushing the most. When one has a large number of particles in a perfectly symmetric configuration all the pushing will average out for an individual particle. This is to some extent what happens in real gases. The pressure itself is mediated by the collision of the molecules the gas is made of.
From a physical point of view of a Lagrangian fluid element, however, one should only be interested in the pressure forces of neighbouring fluid elements exerted on oneself, since the actual equation of motion is . This is the main idea of rpSPH, a particle is accelerated only if a force is acting upon it, i.e. Newtons first law of motion. rpSPH derives its equation of motion directly from equation (2) by subtracting the constant pressure of the particle under consideration from the pressures of all the particles being summed over. Since a gradient is computed, the subtraction of a constant does not change the mathematical meaning of the difference equation. However, as we will demonstrate it dramatically affects the error properties of the entire scheme. The resulting equation of motion reads
(4) 
One immediately notices that this formulation breaks the symmetry between the pairwise forces of particles. Particles that have a pressure difference are both accelerated into the same direction along the pressure gradient. Linear momentum conservation hence will only be achieved if the modeled pressure gradients are resolved. On the other hand if one does not resolved the relevant pressure gradients one cannot possibly get a correct solution to a hydrodynamic problem in any case.
After all, it is important to recall that when constructing conservative schemes one does not necessarily minimize the numerical errors but rather ensures that one is making symmetric errors so that the conserved quantity does not change. Consequently, in rpSPH monitoring the total angular and linear momentum is an indicator of whether one may have resolved the relevant length scales.
Many of the advantages of the entropic function based SPH formalism Springel and Hernquist (2002) stem from avoiding the term that generally is discretized analogous to equation (2). So in this formalism rpSPH is particularly trivial to implement. It involves setting the first term on the right hand side of equation (2) to zero and change the second by subtracting the pressure of the particle under consideration. This literally is achieved by modifying one line of code in Gadget2 Springel (2005) as shown in the last appendix. The resulting scheme saves two multiplies, one division and one addition for one additional subtraction in the main loop over neighbors. So there is no performance penalty in using rpSPH as compared to standard SPH.
rpSPH is seemingly close to equation (3.1) of Monaghan (1992) first discussed by Morris (1996a) which we will refer to as the Morris formulation. It reads,
(5) 
Monaghan dismissed his version for two reasons. The first is that “an isolated pair of particles with different pressures would bootstrap themselves to infinity” and the second is that it is difficult to construct a consistent energy equation. The latter is irrelevant in the formalism evolving an entropic function Lucy (1977); Springel and Hernquist (2002) in which the work does not enter. The first reason we find unappealing since it is actually the correct solution. The simulation having two particles estimates a pressure gradient. So over the model volume, i.e. the two particles and their smoothing volumes there exists a monotonic pressure gradient. Both particles hence should be accelerated along it. Interestingly, Monaghan did not discuss the equivalent case for the symmetric standard SPH. In this case both particles push each other to infinity no matter what. If they have the identical initial pressures their center of mass will not change if they vary their center of mass moves exactly as in rpSPH. In rpSPH they will move together while in SPH they will accelerate each other apart to infinity. We have tested this on a spherical blob of material in vacuum. We set the pressure of the particles after the densities have been computed from kernel smoothing. This way all particles have identical initial pressure. The configuration is completely stable in rpSPH yet blows itself apart in SPH in just a sound crossing time. The reason why we do not choose equation (3.1) of Monaghan (1992) is because we find it to be unstable at least with the leapfrog time integrator in Gadget2 (see Figure 5 below). Another formulation close to rpSPH discretization we could find in the literature is presented Morris et al. (1997a), who chose to subtract a background pressure. This still leaves an equation of motion in which the pressure of the particle under consideration remains part of the hydro force estimate.
An easy way to see why our discretization is valid (Wadsley, 2010, private communication) recognizes as the volume element and sees that equation 4 is equivalent to , which is the same as and is the term we want. The previous version discussed by Morris (1996a) in contrast is the discretized form rather than our . Note that our form is also not the general form suggested by equation of Monaghan (2005) and in this regard is a new formulation. The striking aspect is that in the actual difference form all that is new is that one index that used to be is now . So this literally is a one letter change to codes that implement the Morris formulation. How this can lead to a dramatic change in accuracy becomes obvious from the standard error analysis. Price (2004) gives the error of the summation interpolant in equation (3.9) and the error of the gradient operator in equation (3.11). In a formulation as by Morris (1996a) one discretises and sees that the errors in the interpolation of and the errors of the chosen discretisation multiply. What Morris (1996a) realized was that there are no error terms in his discretisation for constant functions in the gradient operator. However, the complete error terms still end up being the product of the density estimate and the pressure gradient estimate. The advantage of the rpSPH discretisation is that its error terms are the one of a gradient and are not further multiplied by errors of a density interpolation. It also retains the vanishing error for constant functions of the Morris discretisation.
Interestingly, a linear stability analysis reveals that rpSPH has the same dispersion relation as the form of (Morris, 1996a, his equation 10). He has shown that this form

is always stable, independent of the background pressure,

has a numerical sound speed that depends less on the particle spacing as compared to standard SPH,

and does not have unphysical unstable transverse waves in two nor three dimensions when using kernels with compact support.
rpSPH retains all of these advantages while at the same time reducing the discretisation error.
In the standard SPH approach there are in fact infinitely many possible choices for the discretisation of the pressure equation (equation 3.5 in Monaghan, 1992). This is also true for rpSPH. E.g. which suggest the discretisation
(6) 
for any different from zero. We have verified that many choices of work for a variety of test problems. In the following, however, we restrict our attention to the case of .
Whether these theoretically advantageous properties of rpSPH hold up in practice is assessed in a range of test problems in the following section.
3 Tests of rpSPH
We will employ a Courant factor of (i.e. in Gadget where the kernel has a maximal radius of ).
3.1 Reduced Velocity Noise
We start with particles on a periodic regular lattice with , a sound speed and uniform density of unity and zero initial velocities. The particles should stay at rest. However, as we can see in Figure 1 the total kinetic energy in the volume grows rapidly. The total energy in the system is, however, conserved to better than of the initial value for these tests at a value . So the kinetic energy growth in the particle distribution only corresponds to about less than one in one thousand of the total. I.e. the kinetic energy the particles obtain is taken from a slightly decreasing internal energy allowing the total to be conserved to high precision. The lower the neighbour number the faster that growth. The maximum noise reached is controlled by the artificial viscosity. The noise also decreases only very slowly over time after reaching the maximum. This is one of the main reasons why particle settling is so important in SPH simulations. The slow decline also shows why in general settling procedures can be computationally quite intensive. In the same figure we also plot results using rpSPH which dramatically reduces this spurious kinetic energy creation keeping it at zero to machine precision. Rasio (2000) caution that it makes no sense in standard SPH to increase particle numbers while keeping the number of neighbors fixed. Once a neighbour number is reached that keeps noise in the force calculation to a minimum we find rpSPH to be stable while only increasing the particle number. Note that we also have ran these tests dramatically reducing the Courant factor without any improvement in the case of standard SPH.
The thick solid line in Figure 1 uses 20 neighbors which seems optimal for this 2D calculation with the cubic spline kernel. Here one has enough neighbors to estimate the gradients more accurately while still having too few neighbours to show its pairing instability. So one may be tempted to dismiss the finding that one has the large velocity noise as long as one uses the “correct” number of neighbours in one simulation. Unfortunately, this best choice, however, is only applicable at the uniform density. To show this we perturb the positions by a small amount so that the initially uniform positions are changed by adding to them which gives central densities that are about 30% above the mean. We keep again the pressure to be exactly constant by setting the entropy of the gas only once the density has been estimated from kernel smoothing. The thick long dashed line in Figure 1 gives the associated velocity noise. It again is of order one percent of the sound speed and grew very rapidly.
One may also be tempted to to dismiss this particle noise as irrelevant as it only contains less than a tenth of a percent of the total energy of the system. However, we will see in the following it is what leads to unphysical shear viscosity once one considers shear flows further below.
3.2 Absence of the Pairing Instability
The velocity noise we just discussed is unfortunately not isotropic nor is it random. It has a dominant component for velocities towards directions of other particles and is an effect that aids the pairing instability. Our reasoning here is somewhat contrary to the explanations in the literature as e.g. in Schuessler and Schmitt (1981), Vaughan et al. (2008) or Read et al. (2009). These studies suggest that it is the shape of the kernel function that causes clumping instability.
Given that rpSPH does not show any spurious velocities in the uniform density test given above while standard SPH develops clumps within a few sound crossing times already shows that it cannot be the shape of the kernel alone that is relevant here. In the following tests we have looked for any sign of the clumping instability but have not found any evidence for it independent of the number of neighbors we used. This fact is certainly one of the factors that makes our formulation have radically reduced errors in general.
The clumping instability stems from particles pushing each other closer to other particles. With a smaller distance to the other particles the gradient of the kernel becomes smaller and in the next time step the particle gets pushed further away from its initial position. This way particles can pile up in the flat central part of the smoothing kernel. That rpSPH does not show the clumping instability makes one hopeful that even higher order kernels could now be used to further improve the accuracy of the obtained solution.
3.3 Dramatically Reduced Numerical Shear Viscosity
An easy two dimensional setup uses an adiabatic index of in a unit square domain , with periodic boundary conditions. Particles are initialised on an exactly square lattice with a density of so that the initial density estimate from the SPH kernel gives a density estimate of unity to better than four parts in one thousand. We then add a sinusoidal velocity perturbations to this uniform distribution. We set the pressure to to have a sound speed of unity. For the first tests here we only use particles as there are no features to resolve. In all cases we evolve to time .
We choose with . This shear flow setup gives an average kinetic energy of . A detailed discussion of how SPH behaves on this test for different neighbor numbers, particle numbers and viscosity prescriptions is given in the Appendix. In summary, it does very poorly and transfers kinetic energy into heat very rapidly loosing tens of percent in only four sound crossing times (two crossing times of the fastest particles). It also gives more dissipation when using more particles which makes a convergence study at fixed neighbour number impossible. Below, when we discuss rpSPH for viscous flow, we show that the effective numerical viscosity of standard SPH is nonNewtonian and very large which explains why standard SPH is inadequate to modeling fluids in general.
The results for rpSPH are summarized in Figure 2 which is to be compared to the bottom panel for SPH, plotting the kinetic energy in the box as a function of sound crossing times. Note that the axis in the two panels of Figure 2 is different by a factor of 30.
We should note that this test problem when used in typical Cartesian grid codes will give zero numerical dissipation to machine precision because the uniform nature of the flow along the grid axes.
Price (2004) considered seemingly similar tests in his thesis. However, note that there an isothermal equation of state was used. He shows one case with , i.e. a pressure less fluid and another case with . His initial shear profile is like the one we choose here but with twice the amplitude so the fastest particles move with unit speed. The pressureless case is irrelevant here since without pressure forces the discretization of the momentum equation cannot make a difference. For the second case with he gives the result only after one sound crossing time at . The density fluctuations have grown to order of a percent of the initial value after that single sound crossing time. As Figure 2 shows for most reasonable choices of neighbor numbers not much kinetic energy is dissipated over this time also in our simulations with an adiabatic index of . We also repeated this isothermal shear test and find results consistent with Price (2004). This again emphasises that as long one is interested in very few sound crossing times, SPH can give correct results and gives an indication that this is not specific to Gadget.
This problem also allows us to measure the effective Reynolds numbers one can hope to model with standard SPH.
Solving analytically the incompressible NavierStokes equations for our initial conditions because to good approximation only the viscous term is relevant we have
Since the variables are separated we can easily find that
(7) 
showing that the initial shape of the profiles should not change. We can now use equation 7 to get an estimate of the kinematic viscosity from the fraction of kinetic energy remaining
(8) 
where denotes the fraction of the kinetic energy remaining up to time, . This is at four crossing times.
The Reynolds number measures the ratio of inertial forces, , to viscous ones, , where is the characteristic length scale, the mean velocity and is the dynamic viscosity. So with the kinematic viscosity. Despite the ambiguities we may take the quarter wavelength of the velocity perturbation, the root mean squared velocity.
So for the typical values of we found for SPH say 70 (97) per cent remaining equation (8) gives a kinematic viscosity of (). Consequently the numerical Reynolds number . This is very low and lower than most observed transitions between laminar and turbulent flows in the laboratory or terrestrial applications.
The analytic solution in equation (7) only holds if the fluid is Newtonian so that the shear stress can be described by the single number of the kinematic viscosity . Because, Figure 13, shows that the initial velocity profile changed strongly one also concludes that the fluid flow as modeled by SPH is that of a nonNewtonian fluid. So while it would have been interesting to think of SPH as solving effectively a the Navier Stokes rather than the Euler equations we see that this is not the case. The effective numerical viscosity is nonNewtonian and does not in general decrease with numerical resolution.
3.4 Kelvin–Helmholtz Instability
The Kelvin–Helmholtz instability occurs at the interface of two shearing fluids of different densities when velocity perturbations perpendicular to the interface grow to eventually mix the layers. In the inviscid case this is understood analytically (Chandrasekhar, 1961) and the growth time scale for a sinusoidal mode of wavelength between two fluids of density and with a shear velocity between them is
(9) 
This problem is typically setup with an adiabatic index of in a unit square domain , with (e.g. Read et al., 2009)
(10) 
We choose and and the uniform pressure which gives a sound speed of in the low density surrounding medium. This standard setup then perturbs the interface with
(11)  
where we choose and vary . So for our density contrasts the growth times are , 2.68 and 3.49 for the density contrasts 1:2, 1:5, and 1:10, respectively.
It would seem prudent to compare these results to the many investigations that recently have addressed the KH instability using SPH (Agertz et al., 2007; Wadsley et al., 2008; Price, 2008; Read et al., 2009; Hess and Springel, 2009). However, all of them chose a setup that Robertson et al. (2009) showed to be ill defined. While all these studies where concerned with the issue of whether SPH can handle KH instabilities at all, they do not ask whether it actually converges to a correct solution. This can be seen e.g. in Read et al. (2009) where their modified SPH solution compares already visually poorly to the corresponding Eulerian result.
Following Robertson et al. (2009) we opt for a more well defined setup for which they explicitly showed convergence. We modify the initial density, and velocity profile using the “ramp” function
(12) 
choosing so that and the velocity shear is . For the initial velocity perturbation we take setting and .
Figure 3 compares results obtained with Enzo (Bryan and Norman, 1997; Bryan et al., 2001; O’Shea et al., 2004), standard SPH and rpSPH for the two different initial velocity perturbations. The improvement of rpSPH over SPH is dramatic. The billows grow at the correct rate and show a minimum of artificial small scale structure. The figure shows the highest resolution we have calculated. However, even with 120 particles neighbours one can obtain correct results using rpSPH. Also our choice of a high is inconsequential in rpSPH in this incompressible setup since the Balsara switch reduces it dramatically in practice. We left this high value just to show that one can get an accurate solution without having to adjust his viscosity parameter.
3.5 Rayleigh Taylor instability instabilities
Another classic test of a code’s ability to handle subsonic perturbations is the RT instability (e.g Fryxell et al., 2000; Stone and Gardiner, 2007; Stone et al., 2008). SPH has been used to model supernova explosion previously and RT instabilities have been seen (Herant, 1994) as well as global convective instabilities (Fryer, 2004). The physical situation (Chandrasekhar, 1961) here consists of a heavy fluid being supported by a lighter fluid against which initially are in pressure equilibrium with a constant acceleration (e.g. gravity). Remarkably, all idealised test cases that we are aware of use an initially unresolved contact discontinuity and consequently no converged results independent of the method of solution have so far been presented. Instead the differences between different reconstructions schemes, Riemann solvers, meshing, etc. all contribute to the final structures produced in these simulations.
Similarly as the Kelvin–Helmholtz problem above, we chose here initial conditions that try to minimize the computational requirements while yielding converged results. The two dimensional domain is setup with periodic boundary conditions along the direction with and reflecting boundary conditions along with . To achieve the reflecting boundary conditions in Gadget2 we set up the density distribution in and make all particles with and stationary (DSPH_BND_PARTICLES compile option). These particles are not allowed to change their entropy or positions consequently they retain their initial pressure and density. Particles that are trying to penetrate through these walls have their positions changed to be at the wall and velocity vectors reversed. This is only a crude way of modeling reflecting boundaries with SPH but will suffice here to compare between SPH and rpSPH.
In keeping with past literature we use an adiabatic index of , and set up the density at the top to be and at the bottom. So the density profile is with in the cases presented here. The velocity perturbation is applied in direction with and the velocities are set to zero for positions above 0.7 and below 0.3. The pressure is set to to give a sound speed of one at the interface and is set into hydrostatic equilibrium with the constant acceleration in the negative direction with . This gives a pressure difference of 60% between the top and the bottom of the domain. The smaller this pressure difference the more difficult it becomes for SPH to model it. Similarly, the initial velocity perturbation again should not be very much smaller than the sound speed in order to survive viscous damping before a growth time of the instability.
We present the results for this test for velocity perturbations of and in Figure 4.
The differences are dramatic. Where SPH fails completely to see growth of the instability rpSPH gives the expected behaviour for both perturbation strengths.
That rpSPH is dramatic improvement over Morris’ formulation despite only differing in one index is seen in Figure 5. There we give a Rayleigh Taylor problem at low resolution of particles and a density ratio of as further discussed in section 4.1. All parameters were the same. A courant factor of is used, a neighbour number of 40, , Balsara switch is on, and the initial velocity perturbation amplitude is . Clearly our formulation is more stable lending support to our discussion on the different error properties of the two discretisations given above.
3.6 Shock tubes
3.6.1 Sod shock tube
So far we have tested our new formalism only in very weakly compressible situations. We will use the classic Sod shock tube (Sod, 1978) to compare rpSPH to standard SPH here. Rosswog (2009) recently, gave the results for varying viscosity prescriptions and including artificial conduction terms. We change the setup only slightly. The left state has a density and pressure of unity while the right state has a quarter of the density and a pressure of . This test is evolved with an adiabatic index of and we set it up as a two dimensional problem with equal mass particles in a box that extends from zero to ten in and zero to one in the direction. We choose 40 rows of particles in the direction and vary the spacing along to achieve the given densities using a total of particles which are initially at rest. Additionally we set the interface to be at and smooth it with an exponential ramp such that all hydrodynamic variables are given by where we take and and denote the left and right states. We employ periodic boundary conditions which gives us another interface at which has the reversed left and right states but has an initially discontinuous state. This will give us the opportunity to show the difference between smoothed interfaces and artificially sharp to be visible in one figure. For the first results we present we have used 80 neighbors and an artificial viscosity parameter of for both the SPH and the rpSPH calculation. Both employ the Balsara switch to limit the viscosity which will play no role here because there is no rotational component to the flow. Figure 6 summarizes the hydrodynamic state variables at time . To first approximation one gets identical results with the new formalism as compared to the standard approach. Linear momentum is not conserved in rpSPH and we find a linear excess velocity of so per particle an error on the velocity of in the direction and a completely negligible component along the direction. This is at a time when the r.m.s. velocity is so just slightly above one half of a per cent error in the dominant velocity. SPH has poor behaviour at the contact discontinuities. For both the one originating from the initially smoothed and the the discontinuous interface at the right boundary. Both contacts at and are better captured by rpSPH. The Sod shock tube has few features and it is reassuring that using as many as particles can give an excellent answer.
There are only slight differences in how rpSPH handles one dimensional shock tubes. We will discuss one very popular application taken from a cosmological context after testing a very strong shock next.
3.6.2 Strong Shock
Here we give another test of a much stronger shock than the one by sod. This one has a Mach number close to one hundred. We also use the chance to compare this to the difference formulation studied by Morris (1996a). The density and pressure are on the left and on the right. This is very similar to the one studied by Pfrommer et al. (2006) and is well known to work well with standard SPH. Here we use neighbours, and 5000 particles. This is a good example where one can make rpSPH and the Morris formulation give unphysical results. These methods require the pressure gradient to be resolved. So if you start with completely discontinuous left right states one will get unphysical waves giving unexpected results. However, this is not a shortcoming of the method but simply are errors that come from not resolving the initial conditions. We again use the ramp function from above with a width of 4 in this very long domain ranging from 0 to 500 in and to in .
We cannot confirm Morris’ claim that his formulation gives large post shock oscillations in this method and suspect that he may have set up discontinuous initial conditions. We can see that our new formulation performs somewhat better than the Morris formulation as it does not overshoot the analytical density jump of 4 raising the density from to in Figure 7. Otherwise both approaches work fine and have no problem in modeling strong shocks and evolving it for large distances.
3.7 Sedov–Taylor Blast Wave
Another particularly strong shock is formed in the SedovTaylor blast wave (Landau and Lifshitz, 1959) presenting a difficult test problem for incompressible hydrodynamics codes. One the one hand it is a self similar solution which makes it insensitive to how exactly one sets it up as long as one evolves the system for a very long time. On the other hand it is the solution for a point explosion. For a given finite resolution, however, there is no unique way of specifying the initial conditions. Here is where exact momentum and energy conservation is very helpful as one can set up the initial conditions at will and even if one were to make very large errors in the time evolution the method will still arrive at the self similar solution. In conservative grid codes this still can lead to aspherical solutions if one did not resolve the spherical central hot region. Since SPH however uses spherical kernels one can get away sometimes even by just heating one single particle (Springel and Hernquist, 2002). This is very useful in applications such as galaxy formation simulations where one is always far from resolving the relevant length scales of an explosion. On the other hand any physics that were to occur at a scale of the shell thickness would be impossible to resolve in such a single particle energy ejection. For rpSPH and the Morris formulation we need to resolve the pressure gradients in the initial conditions as we saw above in the strong shock setup.
We set up a square lattice of particles with 300 particles on a side in the unit domain. For resolved initial conditions we set a spherical region in the center of radius with the same ramp function as above using a width of , to have a sound speed of one for an adiabatic index of . For both simulations we used a Courant number of 0.2 ( in Gadget), 80 neighbors, artificial viscosity and had the Balsara switch off. Figure 8 shows that there potentially is also an advantage to rpSPH simulations when modeling shocks. We have failed to get standard SPH to give a stable correct density jump of for our setup. Also the three dimensional versions shown by Springel and Hernquist (2002) always seem to be too low by as much as a factor of two. However, standard SPH is much less sensitive to how one sets up the initial conditions and performs much better at low resolutions.
These strong shock problems can work but clearly are not the biggest strengths of rpSPH. However, at this point we simply reused the artificial viscosity prescription which was designed for standard SPH. We believe it is likely that one can find an alternative formulation for the artificial viscosity that fits better into the rpSPH discretisation which may improve its behaviour for highly supersonic conditions. Until a better artificial viscosity prescription is designed one may opt to switch between standard SPH and rpSPH based on the local divergence. We have successfully applied this strategy by using a switch that evaluates the standard SPH sum if and the rpSPH sum for less strongly convergent flow. Here and denote the smoothing length and the current sound speed of the particle. This formulation is robust in all our tests.
3.8 Cosmological Integration of the Santa Barbara Galaxy Cluster
In 1995 a comparison project was initiated that aimed to compare all numerical cosmology codes at that time for relevant realistic initial conditions. The study focused on three dimensional calculations of the formation of a galaxy cluster in the standard CDM scenario of structure formation. The choice was a setup which does not include any other physics than cosmological hydrodynamics with an ideal gas equation of state (often referred to as adiabatic simulations despite the entropy generation in shocks). The study produced a detailed report in (Frenk et al., 1999, F99, herafter). One of the most surprising findings of the study was that while there was very good agreement between the six different SPH codes used in the study they did not agree with the solutions of the grid based codes. The central entropy of the simulated galaxy cluster differed markedly between the grid and SPH codes. In particular, the only AMR code in the study by Bryan and Norman (1997) which is now called Enzo (Bryan et al., 2001; O’Shea et al., 2004) found a flat entropy core while the particle codes found the central entropy to continue to rise towards smaller radii. Note that there has been a significant debate on what the correct solution may be and potential sources between the differences between grid and particle based methods (Springel, 2005; Dolag et al., 2005; Kawata et al., 2009; Wadsley et al., 2008; Mitchell et al., 2009; Agertz et al., 2007; Springel, 2010) is the real reason for the difference.
We do not attempt a resolution study here but simply show how an rpSPH solution compares to an SPH run with otherwise identical parameters and the solution derived with a cosmological AMR code. We use gas as well as dark matter particles for the particles based approach. For the AMR code we again use Enzo already used in F99 using dark matter particles a root grid of cells and seven additional refinement levels. Refinement is based on density thresholds in the baryons and dark matter component. The viscosity parameter in the particle based runs was , and a neighbor number of 300 with tolerance of one was used. The initial redshift is and the initial temperature Kelvin.
A two dimensional slice through the temperature field is given in Figure 9 comparing SPH to rpSPH. Only relatively subtle differences are found. There are perhaps slightly more small scale features visible in the rpSPH calculation and some slight differences in the post shock gas in the main cluster are visible. Slightly more shocking occurs at larger radii towards low density voids in the rpSPH vs. SPH calculation.
Figure 10 compares the solutions using spherically averaged radial profiles as described in F99. The entropy profiles of rpSPH agree better with the AMR than the SPH results. The rpSPH solutions shows the lowest central densities of the three methods and agrees better with the slightly shallower density profile of the grid code. Clearly the differences between all three, however, are rather subtle given that one evolved this system for 13 billion years which are many tens of sound crossing times for the central part of the resulting cluster.
Both SPH and rpSPH simulations have a final linear momentum corresponding to 4.1 and 3.2 km per second per gas particle, respectively. The difference vector between the final gas momenta of the simulation has a magnitude of 2.6 km/s per gas particle. This is a difference of order one half of a percent of the mass weighted mean r.m.s. velocity of km/s. Obviously, giving up the strict linear momentum conservation in our equation of motion has not lead to any noticeable difference in this measure but has improved the comparison with results from adaptive mesh refinement codes.
However, this particular application is relatively easy as dark matter dominates the gravitational potential. As we will show further below rpSPH is quite easy to break with selfgravitating fluids.
4 Further Benefits of rpSPH
4.1 Variable masses
Next we demonstrate that using our pressure force discretisation give another very important advantage. Simulations with drastically varying particle masses continue to give correct results. This is markedly different compared to previous SPH simulations employing particle splitting. The latter only worked reasonably as long as different particle masses were very well separated spatially. As an explicit example we revisit the Rayleigh–Taylor problem from above. This time we initialize particles on a uniform lattice and model the density contrast by changing the particle masses according to the density profile. We employ a density at the top ten times the one of the one at the bottom fluid to demonstrate that this is not just a marginally better aspect of rpSPH. Figure 11 compares the SPH and the rpSPH solution again for . Instead of showing the density field we show the particles painted by circles and colored by their density. This gives us an opportunity to see that rpSPH does not show any clumping instability while it is severe for SPH. In this run we only use 100 by 200 particles demonstrating that rpSPH handles the Rayleigh–Taylor problem very well at small initial perturbations and low particle numbers.
The SedovTaylor blast wave above was also carried out with a staggered grid of particles of varying mass and retains a nice spherical shape despite the multiple squares introduced in the staggered “mesh” of the initial particle distribution.
It is of great interest for a method to be stable under largely varying particle masses. If it is one can use particle splitting likely without worrying too much of how to place the new particles and keep them separate from particles with different masses (e.g. Kitsionas and Whitworth, 2002).
4.2 Formulation for Magnetic Forces
There have been many attempts to implement ideal MHD into the standard SPH formalism (e.g. Gingold and Monaghan, 1977; Phillips and Monaghan, 1985; Price and Monaghan, 2004; Dolag and Stasyszyn, 2009) In our experience rpSPH performs in the hydro part better than the Morris formulation. The latter has been used in implementing ideal MHD into SPH (Morris, 1996b; Price and Monaghan, 2004; Dolag and Stasyszyn, 2009). Therefore, we expect that our new discretization may be of use in this case as well.
A symmetric conservative form of the Lorenz force is generally implemented using the magnetic stress tensor (Phillips and Monaghan, 1985),
(13) 
So the acceleration from the magnetic fields on the th particle is then written as
In the limit where we only have forces from the magnetic pressure gradient only the diagonal of the tensor has terms which are and we recognize this discretisation as exactly the form of equation (2) above. So also the Lorentz force lends itself to be discretised following our new approach. We split the force into the tension component and the magnetic pressure component without loss of generality into
(14) 
In both terms one takes spatial derivatives and hence is allowed to subtract a constant. Choosing for the first one and for the second one avoids finite pairwise forces between particles in regions of constant field. So the isotropic part becomes
(15) 
to which the tension force is added
(16) 
Both the terms simply add to the accelerations in the force calculation. All that is left to do is to replace the fastest signal velocity with the one given in equation (46) of Price and Monaghan (2004) and one has a MHD implementation of SPH. On some simple initial tests this formulation seems to work quite well even without any regularization technique (e.g. Dolag and Stasyszyn, 2009, and references therein) or artificial B field dissipation. A full exploration of the performance of this discretization though is beyond the scope of this contribution.
5 How to break rpSPH
From the tests above it is clear that at least in the weakly compressible limit rpSPH is a very useful improvement over the standard formulation. However, giving up momentum and energy conservation is a big price to pay for those advances and rpSPH cannot possibly replace standard SPH in all problems of interest. rpSPH should be easy to break at low resolutions when one does not resolve the pressure gradients adequately. We now give an illustrative example that makes rpSPH give very bad results which shall serve as a cautionary note and things to look for when applying rpSPH.
The Evrard collapse of a cold gas sphere (Evrard, 1988) has been extensively used for verification of astrophysical SPH codes (e.g. Wadsley et al., 2004; Springel, 2005). In the version that is part of the Gadget distribution it is realized with only 1472 equal mass particles. A centrally concentrated cloud of cold gas collapses, bounces and eventually virializes. Vacuum boundaries are assumed. Very clearly this is only meant to investigate energy conservation of the code and is a simple test running in seconds on ones laptop.
Figure 12 shows the standard SPH solution together with a completely failed solution of rpSPH. We should not that as we increase the resolution rpSPH does converge to the same solution as standard SPH. However, this is a clear example where too little resolution coupled with a solver that does not conserve momentum or energy will fail completely.
Figure 12 also gives the energies in the Evrard collapse for which we formally added a zero and label it rpSPH+0. The term added to the momentum equations right hand side is , which in SPH reads (e.g. Ritchie and Thomas, 2001)
This term which tends to zero when the density and pressure gradients are well resolved. For our poor resolution setup, however, we see that it gives better energy conservation. We do not recommend this formalism over rpSPH, however, since it still has the problems with contact discontinuities and weakly compressible flows. However, it highlights that as for any numerical study resolution studies are crucial and that rpSPH will likely be of little use when one cannot afford sufficient resolution for the particular problem one is interested in.
6 Conclusions
We have presented a novel discretization of the pressure equation for the smoothed particle hydrodynamics, which we call rpSPH that removes the local pressure from the scheme and only considers pressure gradients. This methodology avoids the clumping and banding instability, artificial surface tension, unphysical particle noise, dramatically reduces inherent shear viscosity and numerical dissipation, and allows to realistically evolve density distributions sampled even with disparate particle masses. We have discussed a large number of test in all of which our new discretization outperforms the traditional SPH results. While our approach is not manifestly momentum conserving and easy to break wiht selfgravity and or low resolutions it clearly seems much more accurate than previous approaches to Lagrangian hydrodynamics using SPH. Since our formulation is more accurate and requires as little as one line of code to be changed in previous implementations we do believe it to likely to be useful. The caveat of rpSPH remains that if one knows that one cannot afford to resolve the pressure gradient in ones initial conditions that because it is not momentum conserving it can give very wrong results. Fortunately, a resolution and convergence study can reveal whether one is in this limit.
In summary, some of the biggest shortcomings of SPH can in some circumstances can be overcome if one gives up the idea of applying equal but opposite forces to particle pairs. While the latter is what happens physically to the atoms or molecules making up the fluid it is simply incorrect for Lagrangian fluid elements the particles are meant to represent. Physically it also does not make sense to introduce repulsive forces for two spatially separated points even when there is no pressure gradient between them. To require such symmetry between particles neglects that they are spatially separated and that the gradient of the pressure field is different at the two locations in general. Our new discretization avoids these unphysical forces and allows the SPH particles to behave as Lagrangian volume elements recovering fluid behaviour in a large number of tests.
We have successfully used a fifth order spline kernel giving smaller errors on the uniform shear problem and the RayleighTaylor problems discussed above. Consequently, we believe that further improvements to rpPSPH should be possible in the future.
We have also studied multiple forms of discretising the specific internal energy equation
(17) 
The simplest version that we successfully applied to some of our test problems is given by
(18) 
While we prefer the entropy formulation this form here may be useful for codes that start from an internal energy formulation.
While standard SPH is conservative it fails to correctly capture fluid instabilities and shows large nonNewtonian viscosity. rpSPH, on the other hand is more accurate, but is not inherently momentum or energy conserving. Consequently it is a useful modification to the SPH algorithm when one is studying problems where one can afford to resolve the relevant pressure gradients and the density field.
Acknowledgements
I would like to acknowledge the hospitality of the Zentrum für Astronomie Heidelberg and the Institut für Theoretische Astrophysik in particular where this work was started. I am particularly indebted to Volker Springel for making Gadget public making this work possible in the first place, and thank Paul Clarke, James Wadsley and Ralf Klessen and my group and colleagues at KIPAC for useful discussions. Furthermore, I am indebted to Landesstiftung Baden Württemberg for travel support. This work was partially supported by NASA ATFP grant NNX08AH26G and NSF AST0807312. This research has made use of NASA’s Astrophysics Data System Bibliographic Services and splash developed and described by Price (2007).
Appendix A Shearing flows of Uniform Density with standard SPH
While rpSPH overcomes the problems found for shearing flows this appendix is more of historical interest. However, some readers may find it worthwhile since to our knowledge the surprisingly large shear viscosity and its erratic behaviour with increasing particle numbers has not been documented previously.
A simple two dimensional setup uses an adiabatic index of in a unit square domain , with periodic boundary conditions. All particles are set up on an exactly square lattice with a density of so that the initial density estimate from the SPH kernel in fact gives a density estimate of unity to better than four parts in one thousand. We then add different velocity perturbations to this uniform distribution. We set the pressure to to have a sound speed of unity. For the first tests here we only use particles as there are no features to resolve. In all cases we evolve to time .
First we start with no velocity perturbation. I.e. a completely static uniform density distribution evolved over four sound crossing times. Using 30 neighbors the density estimate by all particles is . For different neighbor numbers this fluctuates around 1 and is close enough. We will use 30 neighbors for most of the rest of this section. Initially all velocities are zero yet after four crossing times we have an r.m.s. velocity with no obvious preferred direction. This results is obtained for the typical viscosity value of . For lower values this random noise increases to for . So clearly even under the most quiet conditions imaginable, a uniform density in pressure equilibrium, we could not represent velocities of order a few percent of the sound speed. Now let us perturb the velocity along the direction and set it to a uniform value of . This should be exactly identical to the previous setup given that SPH is formulated to be Galilean invariant. Now the random noise for and again for .
Now for our next experiment with this uniform density setup we use with . This shear flow setup gives an average kinetic energy of . After only 4 sound crossing times (or two crossing times of the fastest particles) the mean kinetic energy of particles has decreased by 15% and the r.m.s. velocity in direction is already when using a viscosity parameter of . Using the standard value we have a lower r.m.s. velocity in the direction of yet at the same time the total kinetic energy has decreased by a as much as 27 percent suggesting that the standard value does convert unacceptable levels of the shear into heat. This hardly is inviscid flow! Again for but particles which allow the shear to be better resolved one would hope for less dissipation. Yet we find that the total kinetic energy still decreases by 30% and the r.m.s. velocity in the y direction becomes . We have also run this test with particles and find that the kinetic energy dissipation is approximately independent of resolution up to this particle number. The kinetic energy lost after two crossing times was 30.3% and the final r.m.s. velocity fluctuations in the ydirection was i.e. a fiftieth of the sounds speed. The latter velocity dispersion became as high as a thirtieth during the first time interval , decreased and then stayed stable afterwards at . The maximal vertical velocities are as much as one tenth of the sound speed for a problem which should not develop any perpendicular velocities.
Perhaps using a lower viscosity parameter could help? So with and particles this indeed gives lower overall dissipation of the kinetic energy of “only” 16% but then gives random motions perpendicular to the shear of . Figure 2 summarizes the loss of kinetic energies for simulations with and particles for 30 and 48 neighbors and the artificial viscosity parameters of and .
From the figure it is very clear that using more particles actually leads to more dissipation. How is one supposed to carry out a resolution study when the numerical dissipation can keep increasing when using more and more computational resources?
A histogram over all particles showing their velocity as a function of their coordinate in Figure 3 reveals how strongly the shear viscosity turned the initial sinusoidal perturbation into flattened extrema with linear profiles between them. This graph looks similar for different viscosity values. The bottom panel of Figure 3 visualizes the particles making the ones that received the entropy clearly visible. The lowest (initial) value one should expect is white and in fact below the minimum on that image . However, one can see clear bands in the places where one finds the largest gradients in the shear velocities. The clumping instability is clearly visible through the bunching of entropy values in the plot.
It is worth noting that the Balsara switch (Balsara, 1995) which is designed to limit this shear dissipation indeed helps. Without it we find that after two crossing times 48% of the kinetic energy are already artificially dissipated in the test with the sinusoidal shear at Mach one half and a uniform density. These 48% are to be compared to the 30% which were dissipated using the Balsara switch.
The viscosity limiter implemented in Gadget2 does not influence the results here. Also decreasing the Courant factor by an order of magnitude can change the exact amount of dissipation but does in general not decrease it appreciably.
That the effective shear viscosity changed little with increasing particle number is very unfortunate. However, if we increase the neighbor number employed to 60 neighbors the effective shear viscosity drops dramatically and only 1% of the total kinetic energy is artificially dissipated over the same time interval. However this comes at the price of particles clumping into bands through the well known tensile instability. As discussed in some detail by (Read et al., 2009) the amount of clumping is specific to the kernel choice.
One compromise for this uniform shear problem is a neighbor number of 48 which leads to banding and only about 3% kinetic energy dissipation in the two crossing times. Simple scaling implies then a choice of neighbors for three dimensional calculations. A number much larger than typically employed.
We have evolved this same test to many more crossing times and find that at the larger neighbor numbers dissipation simply occurs later but in fact looks qualitatively just like in the low neighbor number case. This enables now a further discussion of the origin of this artificial shear viscosity. In the top panel of Figure 13 we can see how the extrema in the velocity are clipped by the viscosity and that the particles positions which were one a regular lattice in the direction spread out and led to banding. This in large part comes from the nonzero velocity the particles obtain leading them to artificially mix into regions perpendicular to the velocities they have. So small fluctuations in the pressure forces enable a coupling taking energy from the velocity to stimulate motions in the direction. Particles then artificially mix into regions of the flow where they start to interact with fluid parcels of different shear velocities triggering the artificial viscosity which then tries to damp this noise. This is why smaller artificial viscosity parameters will lead to larger r.m.s. velocities. This also explains why increasing the neighbor number delays this artificial shear viscosity in that it decreases the amplitude of the forces leading to the velocities.
Interestingly the uniform density tests presented here are similar as the one presented by Monaghan (2006) in terms of the velocity profile and in the sense that it is a low mach number flow. Monaghan (2006), however, chose to pick making the EOS very stiff. While this may be a useful trick to model incompressible flow with SPH it is not something we repeat here since for most applications in astrophysics we have . However, even with the stiff equation of state his Figure 2 shows the same clipping of the maximal velocity amplitudes as our Figure 13 after only one crossing time. This agrees with our findings that only for short time scales (as compared to the crossing time) the shear viscosity may be negligible. We just differ in the interpretation of whether this is an acceptable level of dissipation or not.
Appendix B Modifying Gadget2.0.4 to rpSPH
For the convenience of other researchers we give the details of what to do to convert Gadget2.0.4^{1}^{1}1http://www.mpagarching.mpg.de/gadget/ to take advantage of the rpSPH discretization. In
ydra.c find te line that reads
and change it to
and the conversion is complete.
Another form that fits more closely to the artificial viscosity prescription useful for problems with large density gradients is
In order to keep standard SPH for very strong shocks matching the standard viscosity implementation one may choose to keep both lines but preface the latter with
, or other criteria that trigger at strong shocks.
References
 Agertz et al. (2007) Agertz, O., Moore, B., Stadel, J., Potter, D., Miniati, F., Read, J., Mayer, L., Gawryszczak, A., Kravtsov, A., Nordlund, Å., Pearce, F., Quilis, V., Rudd, D., Springel, V., Stone, J., Tasker, E., Teyssier, R., Wadsley, J., and Walder, R.: 2007, MNRAS 380, 963
 Artymowicz and Lubow (1994) Artymowicz, P. and Lubow, S. H.: 1994, ApJ 421, 651
 Attwood et al. (2007) Attwood, R. E., Goodwin, S. P., and Whitworth, A. P.: 2007, A&A 464, 447
 Balsara (1995) Balsara, D. S.: 1995, Journal of Computational Physics 121, 357
 Basa et al. (2009) Basa, M., Quinlan, N. J., and Lastiwka, M.: 2009, International Journal for Numerical Methods in Fluids 60(10), 1127
 Børve et al. (2009) Børve, S., Speith, R., and Trulsen, J.: 2009, ApJ 701, 1269
 Bryan et al. (2001) Bryan, G. L., Abel, T., and Norman, M. L.: 2001, in Proceedings of Supercomputing 2001, Vol. arXiv:astroph/0112089
 Bryan and Norman (1997) Bryan, G. L. and Norman, M. L.: 1997, in D. A. Clarke and M. J. West (eds.), Computational Astrophysics; 12th Kingston Meeting on Theoretical Astrophysics, Vol. 123 of Astronomical Society of the Pacific Conference Series, pp 363–+
 Cartwright et al. (2009) Cartwright, A., Stamatellos, D., and Whitworth, A. P.: 2009, MNRAS 395, 2373
 Chandrasekhar (1961) Chandrasekhar, S.: 1961, Hydrodynamic and hydromagnetic stability
 Cummins and Rudman (1999) Cummins, S. J. and Rudman, M.: 1999, Journal of Computational Physics 152(2), 584
 Dolag and Stasyszyn (2009) Dolag, K. and Stasyszyn, F.: 2009, MNRAS 398, 1678
 Dolag et al. (2005) Dolag, K., Vazza, F., Brunetti, G., and Tormen, G.: 2005, MNRAS 364, 753
 Evrard (1988) Evrard, A. E.: 1988, MNRAS 235, 911
 Frenk et al. (1999) Frenk, C. S., White, S. D. M., Bode, P., Bond, J. R., Bryan, G. L., Cen, R., Couchman, H. M. P., Evrard, A. E., Gnedin, N., Jenkins, A., Khokhlov, A. M., Klypin, A., Navarro, J. F., Norman, M. L., Ostriker, J. P., Owen, J. M., Pearce, F. R., Pen, U., Steinmetz, M., Thomas, P. A., Villumsen, J. V., Wadsley, J. W., Warren, M. S., Xu, G., and Yepes, G.: 1999, ApJ 525, 554
 Fryer (2004) Fryer, C. L.: 2004, ApJL 601, L175
 Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W., and Tufo, H.: 2000, ApJS 131, 273
 Gingold and Monaghan (1977) Gingold, R. A. and Monaghan, J. J.: 1977, MNRAS 181, 375
 Graham and Hughes (2008) Graham, D. I. and Hughes, J. P.: 2008, International Journal for Numerical Methods in Fluids 56, 1261
 Herant (1994) Herant, M.: 1994, Memorie della Societa Astronomica Italiana 65, 1013
 Hess and Springel (2009) Hess, S. and Springel, V.: 2009, ArXiv eprints
 Hu and Adams (2007) Hu, X. and Adams, N.: 2007, Journal of Computational Physics 227(1), 264
 Imaeda and Inutsuka (2002) Imaeda, Y. and Inutsuka, S.: 2002, ApJ 569, 501
 Kawata et al. (2009) Kawata, D., Okamoto, T., Cen, R., and Gibson, B. K.: 2009, ArXiv eprints
 Kitsionas and Whitworth (2002) Kitsionas, S. and Whitworth, A. P.: 2002, MNRAS 330, 129
 Landau and Lifshitz (1959) Landau, L. D. and Lifshitz, E. M.: 1959, Fluid mechanics
 Lodato and Price (2010) Lodato, G. and Price, D. J.: 2010, MNRAS 405, 1212
 Lombardi et al. (1999) Lombardi, J. C., Sills, A., Rasio, F. A., and Shapiro, S. L.: 1999, Journal of Computational Physics 152(2), 687
 Lucy (1977) Lucy, L. B.: 1977, AJ 82, 1013
 Mitchell et al. (2009) Mitchell, N. L., McCarthy, I. G., Bower, R. G., Theuns, T., and Crain, R. A.: 2009, MNRAS 395, 180
 Monaghan (1992) Monaghan, J. J.: 1992, ARA&A 30, 543
 Monaghan (1994) Monaghan, J. J.: 1994, Journal of Computational Physics 110(2), 399
 Monaghan (2005) Monaghan, J. J.: 2005, Reports on Progress in Physics 68, 1703
 Monaghan (2006) Monaghan, J. J.: 2006, MNRAS 365, 199
 Morris (1996a) Morris, J. P.: 1996a, Publications of the Astronomical Society of Australia 13, 97
 Morris (1996b) Morris, J. P.: 1996b, Ph.D. thesis, Monash University
 Morris et al. (1997a) Morris, J. P., Fox, P. J., and Zhu, Y.: 1997a, J. Comput. Phys. 136(1), 214
 Morris et al. (1997b) Morris, J. P., Fox, P. J., and Zhu, Y.: 1997b, Journal of Computational Physics 136(1), 214
 Murray (1996) Murray, J. R.: 1996, MNRAS 279, 402
 O’Shea et al. (2004) O’Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., Harkness, R., and Kritsuk, A.: 2004, ArXiv Astrophysics eprints
 Pfrommer et al. (2006) Pfrommer, C., Springel, V., Enßlin, T. A., and Jubelgas, M.: 2006, MNRAS 367, 113
 Phillips and Monaghan (1985) Phillips, G. J. and Monaghan, J. J.: 1985, MNRAS 216, 883
 Price (2004) Price, D.: 2004, Ph.D. thesis, University of Cambridge
 Price (2007) Price, D. J.: 2007, Publications of the Astronomical Society of Australia 24, 159
 Price (2008) Price, D. J.: 2008, Journal of Computational Physics 227, 10040
 Price and Monaghan (2004) Price, D. J. and Monaghan, J. J.: 2004, MNRAS 348, 123
 Rafiee and Thiagarajan (2009) Rafiee, A. and Thiagarajan, K. P.: 2009, Computer Methods in Applied Mechanics and Engineering 198(3336), 2785
 Rasio (2000) Rasio, F. A.: 2000, Progress of Theoretical Physics Supplement 138, 609
 Read et al. (2009) Read, J. I., Hayfield, T., and Agertz, O.: 2009, ArXiv eprints
 Ritchie and Thomas (2001) Ritchie, B. W. and Thomas, P. A.: 2001, MNRAS 323, 743
 Robertson et al. (2009) Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., and Rudd, D. H.: 2009, ArXiv eprints
 Rosswog (2009) Rosswog, S.: 2009, New Astronomy Review 53, 78
 Schuessler and Schmitt (1981) Schuessler, I. and Schmitt, D.: 1981, A&A 97, 373
 Sijacki and Springel (2006) Sijacki, D. and Springel, V.: 2006, MNRAS 371, 1025
 Sod (1978) Sod, G. A.: 1978, Journal of Computational Physics 27, 1
 Springel (2005) Springel, V.: 2005, MNRAS 364, 1105
 Springel (2010) Springel, V.: 2010, MNRAS 401, 791
 Springel and Hernquist (2002) Springel, V. and Hernquist, L.: 2002, MNRAS 333, 649
 Steinmetz and Mueller (1993) Steinmetz, M. and Mueller, E.: 1993, A&A 268, 391
 Stone and Gardiner (2007) Stone, J. M. and Gardiner, T.: 2007, ApJ 671, 1726
 Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., and Simon, J. B.: 2008, ApJS 178, 137
 Swegle et al. (1995) Swegle, J. W., Hicks, D. L., and Attaway, S. W.: 1995, Journal of Computational Physics 116(1), 123
 Vaughan et al. (2008) Vaughan, G. L., Healy, T. R., Bryan, K. R., Sneyd, A. D., and Gorman, R. M.: 2008, International Journal for Numerical Methods in Fluids 56, 37
 Wadsley et al. (2004) Wadsley, J. W., Stadel, J., and Quinn, T.: 2004, New Astronomy 9, 137
 Wadsley et al. (2008) Wadsley, J. W., Veeravalli, G., and Couchman, H. M. P.: 2008, MNRAS 387, 427
 Xu et al. (2009) Xu, R., Stansby, P., and Laurence, D.: 2009, Journal of Computational Physics 228(18), 6703