Dealing with density discontinuities in planetary SPH simulations
Density discontinuities cannot be precisely modelled in standard formulations of smoothed particles hydrodynamics (SPH) because the density field is defined smoothly as a kernel-weighted sum of neighbouring particle masses. This is a problem when performing simulations of giant impacts between proto-planets, for example, because planets typically do have density discontinuities both at their surfaces and at any internal boundaries between different materials. The inappropriate densities in these regions create artificial forces that effectively suppress mixing between particles of different material and, as a consequence, this problem introduces a key unknown systematic error into studies that rely on SPH simulations. In this work we present a novel, computationally cheap method that deals simultaneously with both of these types of density discontinuity in SPH simulations. We perform standard hydrodynamical tests and several example giant impact simulations, and compare the results with standard SPH. In a simulated Moon-forming impact using $10^7$ particles, the improved treatment at boundaries affects at least 30% of the particles at some point during the simulation.