# REMIX SPH – improving mixing in smoothed particle hydrodynamics simulations
using a generalised, material-independent approach

T. D. Sandnes${}^{*,\,}$ V. R. Eke  J. A. Kegerreis  R. J. Massey  S. Ruiz-
Bonilla  M. Schaller  L. F. A. Teodoro

###### Abstract

We present REMIX, a smoothed particle hydrodynamics (SPH) scheme designed to
alleviate effects that typically suppress mixing and instability growth at
density discontinuities in SPH simulations. We approach this problem by
directly targeting sources of kernel smoothing error and discretisation error,
resulting in a generalised, material-independent formulation that improves the
treatment both of discontinuities within a single material, for example in an
ideal gas, and of interfaces between dissimilar materials. This approach also
leads to improvements in capturing hydrodynamic behaviour unrelated to mixing,
such as in shocks. We demonstrate marked improvements in three-dimensional
test scenarios, focusing on more challenging cases with particles of equal
mass across the simulation. This validates our methods for use-cases relevant
across applications spanning astrophysics and engineering, where particles are
free to evolve over a large range of density scales, or where emergent and
evolving density discontinuities cannot easily be corrected by choosing
bespoke particle masses in the initial conditions. We achieve these
improvements while maintaining sharp discontinuities; without introducing
additional equation of state dependence in, for example, particle volume
elements; and without contrived or targeted corrections. Our methods build
upon a fully compressible and thermodynamically consistent core-SPH
construction, retaining Galilean invariance as well as conservation of mass,
momentum, and energy. REMIX is integrated in the open-source, state-of-the-art
Swift code and is designed with computational efficiency in mind, which means
that its improved hydrodynamic treatment can be used for high-resolution
simulations without significant cost to run-speed.

###### keywords:

Smoothed particle hydrodynamics , Fluid dynamics , Mixing , Multi-material

††journal: Journal of Computational Physics

[Durham]organization=Institute for Computational Cosmology, Department of
Physics, Durham University, addressline=South Road, city=Durham, postcode=DH1
3LE, country=UK

[Ames]organization=NASA Ames Research Center, addressline=Moffett Field,
postcode=94035, state=CA, country=USA

[Lorentz]organization=Lorentz Institute for Theoretical Physics, Leiden
University, addressline=PO Box 9506, city=NL-2300 RA Leiden, country=the
Netherlands

[Leiden]organization=Leiden Observatory, Leiden University, addressline=PO Box
9513, city=NL-2300 RA Leiden, country=the Netherlands

[Oslo]organization=Centre for Space Sensors and Systems, Faculty of
Mathematics and Natural Sciences, University of Oslo, addressline=2007
Kjeller, country=Norway

[Glasgow]organization=School of Physics and Astronomy, University of Glasgow,
postcode=G12 8QQ, state=Scotland, country=UK

†† ∗Corresponding author<EMAIL_ADDRESS>
##  1 Introduction

Computational simulations are an invaluable tool for studying the inherently
complex behaviour of fluids. Smoothed particle hydrodynamics (SPH), first
developed by Lucy [1] and Gingold and Monaghan [2], is frequently utilised
across a range of applications spanning astrophysics [3, 4, 5] and engineering
[6, 7]. In astrophysics, it is used in particular for its geometry-independent
adaptive resolution, inherent conservation properties, and elegant coupling
with gravity solvers [8]. For engineering applications, it offers advantages
in the treatment of dynamic free surfaces, fluid–structure interactions, and
in simulating multiphase flow [9, 10]. This is in addition to its relatively
simple construction, numerical stability, and low computational cost. In this
work, our motivation and methods focus on conservative, fully-compressible,
gravity-coupled SPH schemes, where particles have unchanging masses and
material types throughout a simulation.

Two key concepts characterise SPH: the representation of a fluid as a discrete
set of interpolation points, or ‘particles’, that move with the fluid
velocity; and the use of a kernel function to estimate fluid fields and their
gradients at particle positions, by interpolation over neighbouring particles
[11]. However, specific errors are introduced with the assumptions that
underpin these core concepts. The discretisation of the continuous underlying
fluid results in leading-order error in the momentum equation, which is
sensitive to disorder in the local particle distribution [12]. Additionally,
the use of an extended kernel in the traditional, integral form of the SPH
density estimate leads to inadvertent smoothing of interpolated densities [11,
13]. In regions where variations in the underlying density field are not well
resolved by the instantaneous particle configuration, this can lead to the
calculation of spurious particle pressures, and subsequently to spurious
pressure gradients that are used in the equations of motion.

These errors combine particularly strongly at density discontinuities in
simulations where a fluid is represented by particles of fixed, equal mass –
which is the norm for many science applications. In such a case, a density
discontinuity constitutes a sharp change in particle spacing. Both
discretisation and kernel smoothing error combine to give rise to a spurious
surface tension-like effect that greatly suppresses both the mixing of fluid
across the interface and the growth of instabilities that would act to drive
turbulent mixing [14]. This is a well-established shortcoming of SPH, and a
range of approaches have been developed to address these sources of error.

First, we consider methods to reduce discretisation error. Using higher-order
kernel functions with more particle neighbours will generally reduce error
[15], and choices of free functions in the generalised form of the equations
of motion can be exploited to mitigate zeroth-order error [12, 16]. In
conjunction with these, improved gradient estimates from, for example,
reproducing kernels [17, 18, 19] or integral-based gradient estimates [20, 21,
22] have been demonstrated to improve the treatment of fluid mixing and
instability growth. These methods have no dependence on material or equation
of state (EoS) in their construction or underlying assumptions. We therefore
make use of some of these methods in this work.

Next, we consider ways to address kernel smoothing error at contact
discontinuities, many of which explicitly assume the use of a single, ideal
gas EoS. In particular, we note the use of artificial conduction for this
purpose, by which particle internal energies are smoothed over a similar
length scale to the inadvertent density smoothing [23]. This requires
thermodynamic behaviour such that smooth density and internal energy fields
result in a smooth pressure field. Therefore, this cannot reliably improve the
treatment of interfaces between dissimilar materials, represented by different
EoS. Alternatively, methods that use modified density estimates, weighted by a
simple thermodynamic quantity such as specific internal energy, also assume a
simple relationship between density and internal energy at constant pressure,
typically the inversely proportional relationship of an ideal gas [12, 24].

Dealing with kernel smoothing errors at interfaces between arbitrarily
different materials is more challenging since the simplicity of the ideal gas
equation cannot be exploited. A boundary between dissimilar materials in
thermal and pressure equilibrium will in general result in a density
discontinuity, so these problematic scenarios occur frequently in simulations
with multiple materials. Additionally, the surface tension-like effects caused
by the density smoothing are particularly strong for “stiff” EoS, for which
small changes in particle densities can result in large changes in calculated
pressures [25].

Methods to improve the treatment of material interfaces in fully-compressible
SPH formulations have been explored in the context of planetary impacts [26],
where density discontinuities between multiple, stiff materials are common and
can evolve across a range of thermodynamic phase space throughout the course
of a single simulation. The treatment of discontinuous free surface interfaces
is also important in this context [27]. Hosono et al. [28] present a “density
independent SPH” (DISPH [29]) scheme adapted for use with multiple materials.
Here, rather than being calculated from particle masses and densities, volume
elements are based on functions of pressure that are evolved in time and
recalculated to satisfy kernel normalisation in an additional iterative step
[30]. Although this approach prevents spurious pressures at density
discontinuities, specifically in regions with otherwise continuous pressures,
the extension of this method to arbitrary EoS leads to material-dependent
volume elements that intricately depend on fluid thermodynamics. Pearl et al.
[31] present an advanced scheme that, among other improvements, makes use of
Riemann solvers [32] and an optional slip condition at material interfaces.
Their choice of material-dependent density estimate effectively smooths volume
rather than density at material interfaces, in simulations where particles of
the same material have equal masses. If particles are deliberately set up with
equal volumes, then this density estimate will significantly reduce both
kernel-smoothing and discretisation errors. This improvement is evident in
Pearl et al. [31]’s mixing tests, where particles start on a single, ordered
grid. However, such a bespoke particle configuration is not representative of
many science applications. Therefore, these tests do not validate mixing with
their methods for general cases with emergent and evolving density
discontinuities, where particle configurations cannot easily be controlled in
this way throughout the simulation.

This approach of addressing smoothing error through the choice of particle
masses is also taken by Deng et al. [33], who demonstrate enhanced mixing in
their meshless finite-mass (MFM) [34, 35, 36, 37] simulations of planetary
giant impacts. Although MFM includes Riemann solvers and more advanced
gradient estimates that can improve on standard SPH formulations, densities
are still calculated with an interpolated estimate that in this case smooths
volume, and therefore is still subject to kernel smoothing error.
Additionally, a range of SPH modifications specific to material boundaries,
rather than arbitrary density discontinuities, have also been developed [38,
39, 40]. All the methods discussed here that address density smoothing
directly, rather than in the construction of initial conditions, rely on EoS-
or material-dependent treatments in, for example, the calculation of volume
elements and density estimates.

Here, we present the REMIX (Reduced Error MIXing) SPH scheme. REMIX is
constructed with the following goals in mind: (1) to improve the treatment of
density discontinuities and mixing in simulations with both one and multiple
EoS, by directly addressing sources of error in traditional SPH methods; (2)
to be able to achieve this for simulations with particles of equal mass, as
well as less challenging configurations; (3) to retain the key characteristics
of the SPH formalism; (4) to introduce no additional EoS dependence in, for
example, volume elements or density estimates; (5) for computational
efficiency, to require no more than three loops over particle neighbours and
no additional iterative steps compared with the traditional formulation. An
implementation of REMIX is publicly available as part of the open-source Swift
code111Swift is in open development including extensive documentation and
examples at swiftsim.com. [41].

This paper is structured as follows: in §2, we describe elements of the core
SPH formalism that we build on and key sources of error that we address in the
construction of REMIX, as well as the methods used in practice to run our
simulations; in §3, we present each component of the REMIX SPH scheme; in §4,
we validate REMIX in a range of hydrodynamic test simulations; and we
summarise our findings in §5.

##  2 Methods

###  2.1 Smoothed particle hydrodynamics

We first describe the key constituent components of SPH. Although REMIX
includes many improvements to traditional SPH222We use a tSPH formulation
based on that of Price [11], summarised in C, as a basis for our discussion
and for comparisons throughout. (tSPH), we do not deviate far from the core
SPH formalism. Additionally, we take this as an opportunity to describe both
the sources of error in SPH, whose reduction is central to the REMIX
formulation, and the nomenclature and notation that is used throughout. An
additional glossary of notation is included in A.

####  2.1.1 Kernel interpolation and the SPH density estimate

Kernel interpolation theory forms the framework for SPH estimates of fluid
fields and their gradients. In particular, the integral form of the density
estimate is a core component of many SPH schemes [11], by which a smoothed
density field at the position of particle $i$, $\langle\rho_{i}\rangle$, can
be reconstructed from the local spatial distribution of neighbouring particles
$j$, their masses, $m_{j}$, and a kernel function, $W_{ij}$ (described below),
via

$\langle\rho_{i}\rangle=\sum_{j}m_{j}W_{ij}\;.$ (1)

Throughout the governing equations of these SPH schemes, the interpolated
density, $\langle\rho_{i}\rangle$, is used as an estimate of the underlying
density field at the positions of particles, $\rho(\mathbf{r}_{i})$. This
density estimate is a specific application of kernel interpolation, which in
general can be used to reconstruct an arbitrary field, $F$, from its value
sampled at the positions of particle neighbours, via

$\langle F_{i}\rangle=\sum_{j}F_{j}W_{ij}V_{j}\;,$ (2)

where $V_{j}$ are volume elements of particle $j$. In Eqn. 1, volume elements
are taken to be $V_{j}=m_{j}/\rho_{j}$. Kernel interpolation can also be used
in estimates of the gradient333We make the choice of notation, here and
throughout, to express kernel gradients as total derivatives rather than with
“$\nabla$” which is often used to imply derivatives with fixed smoothing
length. In later sections, this allows us to more easily distinguish between
gradient estimates with and without grad-$h$ terms [42]. of $F$,

$\left\langle\frac{dF}{d\mathbf{r}}\bigg{|}_{i}\right\rangle=\sum_{j}F_{j}\frac{dW}{d\mathbf{r}}\bigg{\rvert}_{ij}V_{j}\;,$
(3)

such as in the calculation of pressure gradients and velocity divergences for
the SPH equations of motion.

The smoothing kernel, $W(\mathbf{r}-\mathbf{r^{\prime}},h(\mathbf{r}))$, is a
weighting function with radial extent characterised by the smoothing length
$h$. $W(\mathbf{r}-\mathbf{r^{\prime}},h(\mathbf{r}))$ approaches a delta
function in the limit $h(\mathbf{r})\rightarrow 0$. Traditionally, $W$ is a
positive function with approximately a truncated Gaussian-like shape; the
kernel is typically normalised, and spherically symmetric, as this ensures the
exact interpolation of linear fields in the continuum limit of kernel sampling
(number of neighbours, $N\rightarrow\infty$). For a particle pair $i$, $j$:
$W_{ij}\equiv W(\mathbf{r}_{ij},h_{i})\equiv
W(\mathbf{r}_{i}-\mathbf{r}_{j},h(\mathbf{r}_{i}))$, where
$\mathbf{r}_{ij}\equiv\mathbf{r}_{i}-\mathbf{r}_{j}$. Subscripts denote
quantities either sampled at the position of, or associated with, a particle.
Kernels with a compact support, $H\equiv H(h)$ such that $W(r>H)=0$, are used
to limit the number of neighbours to a finite number. We adopt the convention
of defining the smoothing length $h$ as twice the standard deviation of the
kernel444Despite the ubiquitous use of the nomenclature and notation of the
“smoothing length, $h$”, different definitions are frequently used for both
the relationship between $H$ and $h$, and the method used to calculate $h$ (in
our case, Eqn. 4). Although the differences are subtle, we draw attention to
this as an example of the difficulty of one-to-one comparisons between
simulation codes, especially as methods become increasingly complex. [15].
This relates the smoothing length to the compact support by a constant
multiplication factor $H/h$.

In the SPH construction presented here, a particle’s smoothing length is
evaluated iteratively to satisfy

$h_{i}=\eta_{\rm
kernel}\left(\frac{1}{\sum_{j}W(\mathbf{r}_{ij},\,h_{i})}\right)^{1/d}\;,$ (4)

where $d$ is the spatial dimensionality of the simulation and $\eta_{\rm
kernel}$ is a chosen constant. Eqn. 4 ensures that particles across the
simulation have an approximately constant number of neighbours, determined by
the form of the kernel function and the choice of $\eta_{\rm kernel}$.

We use the Wendland $C^{2}$ kernel [43] for validation of the REMIX SPH
scheme:

$W_{\text{WC2}}\left(|\mathbf{r}-\mathbf{r^{\prime}}|,\,H(h)\right)=\begin{cases}\dfrac{C}{H^{d}}\left(1-\dfrac{|\mathbf{r}-\mathbf{r^{\prime}}|}{H}\right)^{4}\left(1+4\dfrac{|\mathbf{r}-\mathbf{r^{\prime}}|}{H}\right)&\text{for
$|\mathbf{r}-\mathbf{r^{\prime}}|<H$}\;,\\\ 0&\text{otherwise}\;,\end{cases}$
(5)

with $\eta_{\rm kernel}=1.487$ (${\sim}100$ neighbours for $d=3$) [15]. Here
$C=21/(2\pi)$ is the normalisation constant for the Wendland $C^{2}$ kernel in
3D. Higher-order kernels can reduce error but require greater numbers of
neighbours, which can come at a significant cost to code speed. This kernel
offers a suitable compromise between improved accuracy and fast simulation
run-time, a relevant consideration for science applications. In D we
demonstrate the effect of the choice of kernel function in REMIX simulations,
finding that REMIX performs well even with the comparatively low-order cubic
spline kernel with fewer neighbours.

####  2.1.2 SPH equations of motion

The equations of motion govern the kinematic and thermodynamic evolution of
SPH particles. The Euler equations are used as the basis for the SPH equations
of motion for inviscid fluids. These consist of the continuity equation,
momentum equation, and energy equation, which are closed by an “equation of
state”, discussed in §2.2. The general, thermodynamically consistent SPH
equations of motion, where we additionally use the same kernel function across
the equations, take the form [12, 42]

$\displaystyle\frac{d\rho_{i}}{dt}$
$\displaystyle\,=\,\sum_{j}m_{j}\frac{\zeta_{i}}{\zeta_{j}}\mathbf{v}_{ij}\cdot\frac{dW}{d\mathbf{r}}\bigg{|}_{ij}\;,$
(6) $\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle\,=\,-\sum_{j}m_{j}\left(\frac{P_{i}}{\rho_{i}^{2}}\frac{\xi_{i}}{\xi_{j}}+\frac{P_{j}}{\rho_{j}^{2}}\frac{\xi_{j}}{\xi_{i}}\right)\frac{dW}{d\mathbf{r}}\bigg{|}_{ij}\;,$
(7) $\displaystyle\frac{du_{i}}{dt}$
$\displaystyle\,=\,\frac{P_{i}}{\rho_{i}^{2}}\sum_{j}m_{j}\frac{\zeta_{i}}{\zeta_{j}}\mathbf{v}_{ij}\cdot\frac{dW}{d\mathbf{r}}\bigg{|}_{ij}\;,$
(8)

where particle densities, $\rho_{i}$, velocities, $\mathbf{v}_{i}$, and
specific internal energies, $u_{i}$, are evolved in time based on gradients of
pressure, $P_{i}$, and velocity divergences calculated using the relative
velocity of particle pairs
$\mathbf{v}_{ij}\equiv\mathbf{v}_{i}-\mathbf{v}_{j}$. The free functions
$\zeta$ and $\xi$ are introduced in the process of discretisation. An SPH
scheme that explicitly conserves energy and momentum requires antisymmetric
kernel gradient terms in the exchange of particle pairs $i$ and $j$. The
integral form of the density estimate (Eqn. 1) is equivalent to the
differential form (Eqn. 6) in the continuum limit, for $\zeta_{i}=\zeta_{j}$
[12].

####  2.1.3 Kernel smoothing error

A fluid field reconstructed using an extended kernel with $h\neq 0$ will be
affected by smoothing error, even when sampled in the continuum limit [11]. In
the continuum limit, a reconstructed field, $\langle F\rangle$, is the
convolution of the underlying field, $F$, with a smoothing kernel $W$,

$\langle F(\mathbf{r},h)\rangle=\int
F(\mathbf{r^{\prime}})\,W(\mathbf{r}-\mathbf{r^{\prime}},h)\,dV^{\prime}\;.$
(9)

Eqn. 2 is the discretised form of this equation. Assuming a continuous,
infinitely differentiable field $F$, we can Taylor expand about the point
$\mathbf{r}$ to give

$\displaystyle\langle F(\mathbf{r},h)\rangle\,=\,$
$\displaystyle\;F(\mathbf{r})\,\cancelto{1}{\int\\!\\!W(\mathbf{r}-\mathbf{r^{\prime}},h)\,dV^{\prime}}~{}~{}~{}+~{}~{}~{}\frac{dF}{dr^{\,\alpha}}\bigg{|}_{\mathbf{r}}\,\,\cancelto{\mathbf{0}}{\int\\!(\mathbf{r^{\prime}}-\mathbf{r})^{\alpha}\,W(\mathbf{r}-\mathbf{r^{\prime}},h)\,dV^{\prime}}$
$\displaystyle~{}+~{}~{}\frac{1}{2}\frac{d^{2}\\!F}{dr^{\,\alpha}dr^{\,\beta}}\bigg{|}_{\mathbf{r}}\int\\!(\mathbf{r^{\prime}}-\mathbf{r})^{\alpha}\,(\mathbf{r^{\prime}}-\mathbf{r})^{\beta}\,W(\mathbf{r}-\mathbf{r^{\prime}},h)\,dV^{\prime}~{}~{}~{}+~{}~{}~{}\dots$
(10) $\displaystyle\,=\,$
$\displaystyle\;F(\mathbf{r})~{}~{}+~{}~{}\frac{1}{2}\frac{d^{2}\\!F}{dr^{\,\alpha}dr^{\,\beta}}\bigg{|}_{\mathbf{r}}\int\\!(\mathbf{r^{\prime}}-\mathbf{r})^{\alpha}\,(\mathbf{r^{\prime}}-\mathbf{r})^{\beta}\,W(\mathbf{r}-\mathbf{r^{\prime}},h)\,dV^{\prime}~{}~{}~{}+~{}~{}~{}\dots\;,$
(11)

where Greek letter superscripts correspond to spatial dimensions, and like
indices are summed over [11, 44]. We separate the first two terms in Eqn. 10
to demonstrate that, in the continuum limit, the choice of a normalised,
spherically symmetric kernel results in the zeroth- and first-order integrals
of the expansion taking values $1$ and $\mathbf{0}$ respectively.

In the continuum limit, kernel interpolation will only reproduce
$F(\mathbf{r})$ without error if the integrals of the second and higher order
terms are all equal to zero. Due to the assumed symmetry properties of the
kernel, integrals in odd terms of the expansion are trivially equal to 0,
while in general even terms will be non-zero. For a positive kernel these non-
zero terms act to smooth the reconstructed field. Although the integrals in
Eqn. 11 will be of order $h^{2}$ and higher powers of $h$, with exponents
corresponding to the term of the expansion [13], the errors become significant
in regions where second and higher order derivatives of the underlying field
are large over length scales of $h$ [45]. This is, in particular, the case for
an underlying field that approaches a discontinuity relative to $h$-length
scales. A discontinuity, where the field is not differentiable, will
inevitably be erroneously smoothed by kernel interpolation.

The integral SPH density estimate, Eqn. 1, is an example of the discrete form
of Eqn. 9. Through Eqn. 11, we see how a quantity calculated by kernel
interpolation in this way will experience smoothing error when the underlying
field varies sharply over $h$-length scales. At density discontinuities,
smoothing of the density field leads to spurious pressures that contribute to
surface tension-like effects that impede particle mixing across the interface.

####  2.1.4 Discretisation error

The kernel smoothing errors discussed above are in addition to, and separate
from, errors introduced by discretisation [11, 46]. Discretisation errors
manifest themselves both through the choice of free functions in the equations
of motion – affecting how closely Eqns. 6–8 approximate their continuous Euler
equation equivalents – and through the imperfect sampling of the kernel by a
finite number of particle neighbours, i.e. in the discretisation of integrals
like Eqn. 9.

The use of a normalised, spherically symmetric kernel leads to the exact
reconstruction of linear fields in the continuum limit by Eqn. 9, as the
higher-order derivatives in Eqn. 11 are zero by construction. However, in the
process of discretisation of the fluid into a finite set of particles, the
conditions

$\displaystyle\sum_{j}W_{ij}V_{j}$ $\displaystyle=1\;,$ (12)
$\displaystyle\sum_{j}\mathbf{r}_{ij}W_{ij}V_{j}$
$\displaystyle=\mathbf{0}\;,$ (13)

are no longer enforced. The exact reconstruction of fluid fields is therefore
lost, even to zeroth order. The amount of discretisation error is a function
of the disorder in the local particle distribution. This also applies to
gradient estimates, such as those used in the equations of motion.
Furthermore, in the equations of motion, gradient estimates are typically
modified to enforce conservation, so generally deviate further from exact
reproduction of underlying linear fields.

In SPH simulations where a fluid is represented by particles of equal mass, a
density discontinuity constitutes a sharp change in particle spacing and thus
large local anisotropies in particle distribution. This leads to
discretisation error also playing a considerable role in suppressing mixing at
density discontinuities [12].

###  2.2 Equations of state

The EoS characterises the thermodynamic behaviour of a material. In SPH
simulations, hydrodynamical evolution is tied directly to pressures and sound
speeds, calculated through the EoS. Many applications in astrophysics use
simulations with only a single, ideal gas EoS. However, in some cases,
multiple EoS are required to simulate dissimilar materials or phases, such as
for planetary impacts, where EoS are often highly complex [25]. The
improvements offered by the REMIX SPH scheme are EoS-independent, and so our
methods can be applied effectively to these simulations, as well as other
applications with multiphase fluids.

For the hydrodynamic test simulations presented in §4, we validate the REMIX
scheme using both ideal gases and more complex EoS. For ideal gas simulations,
the adiabatic index, $\gamma$, is problem-specific and chosen to draw
comparisons with past work. For simulations using more complex materials, we
use EoS typically used for planetary impact simulations. In most of these
tests, we consider iron and rock in conditions representative of the
core–mantle boundary in an Earth-like planet. We use the updated ANEOS
Fe85Si15 and forsterite EoS for these materials, respectively [47]. For
simplicity, we hereafter refer to these as “iron” and “rock”. In §4.9, we also
consider a Jupiter-like planet. For these simulations, we use the
hydrogen–helium EoS from Chabrier and Debras [48], with a helium mass fraction
of $Y=0.245$, and the AQUA EoS from Haldemann et al. [49] to represent heavy
elements or ice.

We note that in the simulations we present here, these materials are treated
as fluids without physical viscosities or strength properties.

###  2.3 The Swift code

Swift is a state-of-the-art, open-source hydrodynamics and gravity code that
specialises in SPH simulations for planetary applications as well as galaxy
formation and cosmology [41, 50]. By using task-based parallelism,
asynchronous communications, and graph-based decomposition of the work between
compute nodes, Swift can perform high-resolution simulations efficiently on
modern high-performance computing architectures [51]. REMIX is fully
integrated into and was developed using the Swift code, and is therefore
publicly available555Swift is available at www.swiftsim.com alongside
extensive documentation and a large suite of examples.. All simulations
presented here were carried out using the Swift code and all tests shown below
are shipped with the code package. Algorithms used for the neighbour-finding,
time-stepping, and gravity in our simulations are detailed in Schaller et al.
[41] and are used identically for simulations with both REMIX and traditional
SPH.

##  3 REMIX SPH

In this section, we detail the constitutive equations of the REMIX SPH
scheme666The full set of the final equations used in the REMIX scheme are
listed in B. The equations of the traditional SPH scheme that we use for
comparison simulations are listed in C.. We improve the treatment of mixing by
directly addressing the sources of SPH error discussed in §2.1. By targeting
both smoothing and discretisation error, we alleviate spurious surface
tension-like effects at density discontinuities, including in the more
challenging cases with equal-mass particles and at interfaces between
dissimilar, stiff materials. Note that we aim to address mixing at the
particle scale and not below. Therefore, we do not consider diffusion of
material type between particles, meaning that the material of each particle
remains fixed for the duration of the simulation.

We target error by exploiting three key freedoms in the SPH equations of
motion presented in §2.1.2: in the choice of density estimate (§3.1); in the
choice of free functions (§3.2); and in the form of the kernel function
(§3.3). Additionally, we develop a novel method that enables the appropriate
treatment of free surfaces when using these improved kernels (§3.4), and we
use improved artificial viscosity (§3.5) and artificial diffusion (§3.6)
formulations. These include new approaches both for the treatment of shocks
and to weakly smooth and mitigate accumulated noise on the particle scale. We
also include a term in the density evolution that re-ties densities to the
local particle distribution (§3.7). These components combine into the REMIX
equations of motion, given by

$\displaystyle\frac{d\rho_{i}}{dt}$
$\displaystyle=\sum\limits_{j}m_{j}\,\frac{\rho_{i}}{\rho_{j}}\,v_{ij}^{\alpha}\,\frac{1}{2}\left(\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ji}\right)\>+\>\left(\frac{d\rho_{i}}{dt}\right)_{\text{difn}}\>+\>\left(\frac{d\rho_{i}}{dt}\right)_{\text{norm}}\;,$
(14) $\displaystyle\frac{dv_{i}^{\alpha}}{dt}$
$\displaystyle=-\sum\limits_{j}m_{j}\,\frac{P_{i}+Q_{ij}+P_{j}+Q_{ji}}{\rho_{i}\,\rho_{j}}\,\frac{1}{2}\left(\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ji}\right)\;,$
(15) $\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=\sum\limits_{j}m_{j}\,\frac{P_{i}+Q_{ij}}{\rho_{i}\,\rho_{j}}\,v_{ij}^{\alpha}\,\frac{1}{2}\left(\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{dr^{\alpha}}\bigg{|}_{ji}\right)\>+\>\left(\frac{du_{i}}{dt}\right)_{\text{difn}}\;,$
(16)

where
$\left(d\widetilde{\mathcal{W}}\\!/d\mathbf{r}|_{ij}-d\widetilde{\mathcal{W}}\\!/d\mathbf{r}|_{ji}\right)/2$
are improved kernel gradient terms that are antisymmetric in the exchange of
$i$ and $j$ for explicit conservation of momentum and energy; $Q_{ij}$ and
$Q_{ji}$ are pairwise, artificial viscous pressures;
$\left(d\rho_{i}/dt\right)_{\text{difn}}$ and
$\left(du_{i}/dt\right)_{\text{difn}}$ are artificial diffusion of density and
internal energy; and $\left(d\rho_{i}/dt\right)_{\text{norm}}$ is the kernel
normalising term. Each of these are discussed in detail in their corresponding
sections below.

The equations of the REMIX scheme were developed to be implemented in just
three loops over particle neighbours, and without introducing any additional
iterative steps. In our test simulations, performed on the COSMA8 HPC
system777Simulations carried out on COSMA8 used 1 node with 128 cores and
those on COSMA7 used 1 node with 28 cores. These are both part of the DiRAC
cluster hosted by Durham University (https://dirac.ac.uk/memory-intensive-
durham/)., using REMIX led to a run-speed ${\sim}1.3$–$1.6$ times longer than
equivalent simulations performed with traditional SPH (and everything else
unchanged). The exact amount of slowdown is problem-dependent: this range
includes simulations both with and without gravity, and those using different
kernel functions888Simulations used to investigate the runtime were: 3D
Kelvin–Helmholtz instabilities (§4.3.2) and planets in hydrostatic equilibrium
(§4.9). These were tested with cubic spline and Wendland $C^{2}$ kernels.. On
the COSMA7 HPC system (which has fewer cores per node), simulations with the
overhead of gravity take ${\sim}1.6$–$1.8$ times longer, and simulations
without gravity take ${\sim}2$–$3$ times longer, depending on the test case.
We find that REMIX, in addition to dealing with density discontinuities that
are problematic in traditional SPH at all resolutions, is able to achieve an
improved treatment of non-discontinuous regions in simulations with over an
order of magnitude lower resolution compared with equivalent traditional SPH
results (§4.3.1). The effective slowdown from using REMIX is therefore much
smaller in practice than the ranges above suggest, since simulations with a
lower resolution (fewer SPH particles) could be used to obtain equivalent
results. As such, in many cases a science simulation with REMIX would run
faster than a traditional SPH simulation that would require a higher
resolution to achieve a comparable level of numerical convergence. For
example, the $2.9\times 10^{5}$ particle REMIX Kelvin–Helmholtz instability in
§4.3.1 runs over 20 times faster (on COSMA8) than the $4.7\times 10^{6}$
particle traditional SPH simulation, and is closer to the converged
solution999See REMIX, $N=128$ and tSPH, $N=512$ in Fig. 4..

###  3.1 Density estimate

In the REMIX SPH scheme we use a differential form of the density estimate: we
evolve the density in time with Eqn. 14 rather than recalculating it each
timestep (e.g. Eqn. 1), similarly to internal energy in traditional SPH
schemes. There are three key benefits of this treatment: (1) we directly
address systematic smoothing error in particle densities, which is
particularly significant at density discontinuities, including those at free
surfaces; (2) it allows us to constrain zeroth-order error in the equations of
motion while starting from a basis of thermodynamic consistency (§3.2); (3) we
do not require an additional loop over particle neighbours to calculate a new
density each timestep. We note that particle mass is fixed throughout the
simulation, so the evolution of densities is equivalent to an evolution of
volumes. In §4.4, we show the differences in Kelvin–Helmholtz instability
simulations when using the full REMIX scheme, and the REMIX scheme modified to
use a traditional integral density estimate. Using our evolved density
estimate, both to calculate thermodynamic quantities and in volume elements,
leads to a considerable improvement in addressing spurious surface tension-
like effects that suppress instability growth and mixing on the particle
scale.

In practice, we set a density floor $\rho_{\rm
min,\,i}\equiv\langle\rho_{i}\rangle_{\rm min}=m_{i}W(\mathbf{0},h_{i})$ such
that $\rho_{i}=\rho_{\rm min,\,i}$ if the density would evolve below the
minimum value. This prevents EoS extrapolation issues that arise for tiny
densities in simulations involving a vacuum region.

Evolved density estimates are used frequently in SPH schemes developed for
engineering applications [52] as well as in some astrophysical SPH schemes, in
particular those that include material strength models [53]. However, in most
astrophysical SPH schemes, an integral density estimate is preferred for its
robustness: the accumulation of error in an evolved density estimate is less
predictable than the relatively controlled errors in a density estimate
calculated each timestep from the instantaneous local particle distribution.
For instance, if left to evolve freely over many timesteps, densities could in
principle take values such that volume elements $m_{j}/\rho_{j}$ are far from
normalising the kernel $W_{ij}$, despite the kernel being a normalised
function101010Volume elements that use the interpolated density,
$V_{j}=m_{j}/\langle\rho_{j}\rangle$, are inherently tied to kernel
normalisation. The equations for kernel normalisation, Eqn. 12, and the
integral density estimate, Eqn. 1, are equivalent to each other in the limit
of constant density on the kernel length scale,
$\langle\rho_{j}\rangle\rightarrow\langle\rho_{i}\rangle$ for all $j$.. We
address these concerns with four approaches: (1) by introducing a novel term
in the density evolution that re-ties densities to the local particle
distribution (§3.7); (2) by using kernels that are normalised to the evolved
densities (§3.3); (3) by including a weak density diffusion to smooth out
accumulated noise (§3.6); (4) and by taking preventative measures in reducing
error that could accumulate with time, reflected in the choices of our
equations of motion (§3.2), the use of kernel functions constructed to reduce
discretisation error (§3.3), and our improved viscosity formulation (§3.5).

Evolved densities are used wherever density appears in the equations of the
REMIX scheme. This includes for calculating thermodynamic quantities, using
the equation of state, and in all volume elements.

###  3.2 Free functions in the equations of motion

In traditional SPH formulations, the free functions, $\zeta$ and $\xi$, in the
equations of motions (Eqns. 6–8) typically take equal values for all particles
and cancel. An alternate formulation with $\zeta=\xi=\rho$, such that the
equations of motion include ratios of the densities of particles $i$ and $j$,
helps to constrain error in the equations of motion at density discontinuities
and for irregular particle distributions on the kernel scale [12]. This choice
avoids the use of gradients of density in the derivation of the momentum
equation, by using the identity

$\frac{\nabla P}{\rho}=\frac{\nabla P}{\rho}+\frac{P}{\rho}\nabla 1\;,$ (17)

rather than

$\frac{\nabla
P}{\rho}=\nabla\left(\frac{P}{\rho}\right)+\frac{P}{\rho^{2}}\nabla\rho\;.$
(18)

SPH formulations using the density as the free functions have been shown to
improve the treatment of mixing [16]. For simulations using only a single
ideal gas, the choice of $1/u$ as a free function is equivalent to this, with
the additional assumption of constant pressure on the kernel scale [24].

Using density as the free function in the integral form of the density
estimate (Eqn. 1) for simulations with arbitrary EoS is not possible without
iteration, since the density would be needed in the density calculation.
However, using the differential form to evolve the density (Eqn. 14) enables
us to develop the REMIX SPH scheme from a basis of full thermodynamic
consistency with $\zeta_{i}=\rho_{i}$. We also use $\xi_{i}=\rho_{i}$ to
reduce zeroth-order error in the momentum equation [12, 16]. All densities
used are the evolved densities of particles.

In §4.4, we demonstrate the improvements in REMIX simulations of the
Kelvin–Helmholtz instability from using $\zeta_{i}=\xi_{i}=\rho_{i}$, compared
with the REMIX scheme modified to use traditional equal-valued free functions.

###  3.3 Linear-order reproducing kernels

To reduce discretisation error, we construct kernels that explicitly satisfy
the conditions given by Eqns. 12 and 13. Therefore, these kernels reproduce
exact values for fields that are spatially constant or that vary linearly with
position. This methodology is largely based on that of Frontiere et al. [19].
To account for spatial variations of the smoothing length, we include grad-$h$
terms that were previously neglected. These grad-$h$ terms take a non-standard
form, compared with Hopkins [42], since our evolved density is not tied
directly to smoothing lengths through the instantaneous distribution of
particles. We also modify our kernels to include a free-surface treatment
(§3.4) to allow them to appropriately handle vacuum boundaries.

The modified kernel, $\mathcal{W}_{ij}$, is constructed so that the sum over
neighbours always satisfies

$\displaystyle\sum_{j}\mathcal{W}_{ij}V_{j}$ $\displaystyle=1\;,$ (19)
$\displaystyle\sum_{j}\mathbf{r}_{ij}\mathcal{W}_{ij}V_{j}$
$\displaystyle=\mathbf{0}\;.$ (20)

We use volume elements $V_{j}=m_{j}/\rho_{j}$, where $\rho_{j}$ are the
evolved densities. We stress that for use in the equations of motion we must
undergo a necessary step to make the kernel gradient terms antisymmetric in
exchanges of particle pairs, to enforce the conservation of energy and
momentum, as is also done by Frontiere et al. [19]. Therefore, the gradient
estimates used in the equations of motion end up being not exactly first-order
reproducing. Despite this, these gradient estimates show significant
improvements when compared with unmodified kernels (as seen directly in §4.4).

To construct $\mathcal{W}_{ij}$, an unmodified SPH kernel is multiplied by a
linear polynomial

$\mathcal{W}_{ij}\equiv
A_{i}\left(1+B_{i}^{\alpha}r_{ij}^{\alpha}\right)\overline{W}_{ij}\;,$ (21)

where
$\overline{W}_{ij}\equiv[W(\mathbf{r}_{ij},h_{i})+W(\mathbf{r}_{ji},h_{j})]/2$
is a symmetrised kernel111111We find this to be beneficial when we enforce the
antisymmetrisation required for use in the equations of motion, as
demonstrated in E. Note that for certain computational steps, this choice
extends the definition of particle $i$’s “neighbours”, $j$, to be those that
satisfy either $|\mathbf{r}_{ij}|<H_{i}$ or $|\mathbf{r}_{ij}|<H_{j}$ rather
than just the first condition., and $A_{i}$ and $\mathbf{B}_{i}$ are
coefficients that satisfy Eqns. 19 and 20, as shown in Appendix A of Frontiere
et al. [19]:

$\displaystyle A_{i}$
$\displaystyle=\left(\overline{m}_{0,\,i}-\left(\overline{m}_{2,\,i}^{\,-1}\right)^{\alpha\beta}\overline{m}_{1,\,i}^{\,\alpha}\,\overline{m}_{1,\,i}^{\,\beta}\right)^{-1}\;,$
(22) $\displaystyle B_{i}^{\alpha}$
$\displaystyle=-\left(\overline{m}_{2,\,i}^{\,-1}\right)^{\alpha\beta}\overline{m}_{1,\,i}^{\,\beta}\;,$
(23)

where the geometric moments are defined as

$\displaystyle\overline{m}_{0,\,i}$
$\displaystyle=\sum_{j}\overline{W}_{ij}V_{j}\;,$ (24)
$\displaystyle\overline{m}_{1,\,i}^{\,\alpha}$
$\displaystyle=\sum_{j}r_{ij}^{\alpha}\overline{W}_{ij}V_{j}\;,$ (25)
$\displaystyle\overline{m}_{2,\,i}^{\,\alpha\beta}$
$\displaystyle=\sum_{j}r_{ij}^{\alpha}r_{ij}^{\beta}\overline{W}_{ij}V_{j}\;.$
(26)

Greek letter indices correspond to spatial dimensions and like indices are
summed over. Bars indicate the use of the symmetrised kernel in the kernel
interpolation. This distinction becomes important since we use $m_{0,\,i}$,
calculated similarly but using an unsymmetrised kernel, for alternative
gradient estimates used later in this section and in §3.7.

To calculate gradient terms for the equations of motion, we require the
spatial derivative of $\mathcal{W}$. We include terms that depend on the
gradient of smoothing lengths, unlike Frontiere et al. [19]. We find the
effects of these to be small in practice, but include them for completeness of
the method – without assuming these to be negligible.

The smoothing length dependence of Eqns. 21–26 is contained within
$\overline{W}_{ij}$. We therefore express the derivatives with the
parameterisation
$\mathcal{W}(\mathbf{r},\,\mathbf{r}_{j})\equiv\mathcal{W}\left(\mathbf{r}-\mathbf{r}_{j},\;\overline{W}(\mathbf{r}-\mathbf{r}_{j},h(\mathbf{r}),h(\mathbf{r}_{j})),\;A(\mathbf{r}),\;\mathbf{B}(\mathbf{r})\right)$,
giving

$\displaystyle\frac{d\mathcal{W}}{dr^{\gamma}}$
$\displaystyle=AB^{\gamma}\overline{W}+\frac{\partial\mathcal{W}}{\partial\overline{W}}\frac{d\overline{W}}{dr^{\gamma}}+\frac{\partial\mathcal{W}}{\partial
A}\frac{dA}{dr^{\gamma}}+\frac{\partial\mathcal{W}}{\partial
B^{\alpha}}\frac{dB^{\alpha}}{dr^{\gamma}}\;.$ (27)

When evaluated for a particle pair $i,j$ this becomes121212We use the notation
$\frac{d\mathcal{W}}{dr^{\gamma}}\big{|}_{ij}\equiv\frac{d\mathcal{W}}{dr^{\gamma}}\left(\mathbf{r}_{i}-\mathbf{r}_{j},\;\overline{W}(\mathbf{r}_{i}-\mathbf{r}_{j},h(\mathbf{r}_{i}),h(\mathbf{r}_{j})),\;A(\mathbf{r}_{i}),\;\mathbf{B}(\mathbf{r}_{i})\right)\equiv\frac{d\mathcal{W}}{dr^{\gamma}}\left(\mathbf{r}_{ij},\;\overline{W}_{ij},\;A_{i},\;\mathbf{B}_{i}\right)$.

$\displaystyle\frac{d\mathcal{W}}{dr^{\gamma}}\bigg{|}_{ij}$
$\displaystyle=A_{i}B_{i}^{\alpha}\overline{W}_{ij}+A_{i}\left(1+B_{i}^{\alpha}r_{ij}^{\alpha}\right)\frac{d\overline{W}}{dr^{\gamma}}\bigg{|}_{ij}+\left(1+B_{i}^{\alpha}r_{ij}^{\alpha}\right)\overline{W}_{ij}\frac{dA}{dr^{\gamma}}\bigg{|}_{i}+A_{i}r_{ij}^{\alpha}\overline{W}_{ij}\frac{dB^{\alpha}}{dr^{\gamma}}\bigg{|}_{i}\;.$
(28)

Equations to calculate the gradients of $A$, $B$, and the geometric moments
are included in B. The derivative of the symmetrised kernel is given
by131313Note that we are taking the derivative of the continuous function
$\overline{W}(\mathbf{r},\mathbf{r}_{j})=[W(\mathbf{r}-\mathbf{r}_{j},h(\mathbf{r}))+W(\mathbf{r}_{j}-\mathbf{r},h(\mathbf{r}_{j}))]/2$
with respect to $\mathbf{r}$, with fixed neighbour positions $\mathbf{r}_{j}$,
and evaluating it at $\mathbf{r}_{i}$. Therefore, there is only a grad-$h$
term associated with the first term in the brackets.

$\frac{d\overline{W}}{dr^{\gamma}}\bigg{|}_{ij}=\frac{1}{2}\left(\frac{\partial
W}{\partial r^{\gamma}}\bigg{|}_{ij}+\frac{\partial W}{\partial
h}\bigg{|}_{ij}\frac{dh}{dr^{\gamma}}\bigg{|}_{i}-\frac{\partial W}{\partial
r^{\gamma}}\bigg{|}_{ji}\right)\;,$ (29)

and so the inclusion of grad-$h$ terms in the gradient calculations in
practice only takes the form of the additional term in Eqn. 29. Both $\partial
W/\partial r^{\gamma}$ and $\partial W/\partial h$ can be calculated directly
from the kernel function [54].

Finally, we require $dh/dr^{\gamma}$. In SPH schemes that use the traditional
density estimate, $dh/dr^{\gamma}$ do not need to be calculated explicitly
[42], since smoothing lengths and densities are inherently linked. However,
for the scheme presented here, where we use an evolved density estimate, we
must calculate this explicitly. One approach is to directly differentiate Eqn.
4. However, we find that zeroth-order error from calculating grad-$h$ terms in
this way leads to spurious behaviour in simulations. We therefore calculate
these by kernel interpolation. Since, in practice, $d\overline{W}/dr^{\gamma}$
has not been constructed yet due to the order of these operations in the loops
over particle neighbours, we are unable to use these improved gradient terms
for $dh/dr^{\gamma}$ if we want to avoid introducing a 4th loop. This also
applies for gradient estimates in our viscosity (§3.5) and diffusion (§3.6)
schemes, discussed later. We therefore require an alternative gradient
estimate for these calculations. However, we must be mindful of kernel
normalisation in these alternative gradient estimates, since we use evolved
densities for volume elements throughout. We therefore use the kernel gradient
term

$\partial_{i}^{\gamma}\hat{W}_{ij}\equiv\frac{\partial_{i}^{\gamma}W_{ij}}{m_{0,\,i}}-\frac{W_{ij}}{m_{0,\,i}^{2}}\partial_{i}^{\gamma}m_{0}\;,$
(30)

where we note that the lack of bars throughout indicates the use of a standard
(e.g. Wendland $C^{2}$) kernel, rather than one symmetrised by averaging with
neighbouring kernels, and ‘$\partial$’, rather than total derivatives,
indicates a lack of grad-$h$ terms. These choices allow us to calculate these
kernel gradients in two loops over particle neighbours, so they can be used
here and in the artificial viscosity and diffusion schemes. Circumflexes, here
and throughout, indicate the use of the normalised kernel $\hat{W}_{ij}\equiv
W_{ij}\,/\,m_{0,\,i}$.

We then calculate

$\partial_{i}^{\gamma}\hat{h}=\sum_{j}(h_{j}-h_{i})\,\partial_{i}^{\gamma}\hat{W}_{ij}\frac{m_{j}}{\rho_{j}}\;,$
(31)

and use this in place of $dh/dr^{\gamma}$.

All these equations combine in Eqn. 28 to give the gradients of the linear-
order reproducing kernels. The use of these kernels reduces discretisation
error in the equations of motion. In §4.4, we show the effect of these kernels
on simulations of the Kelvin–Helmholtz instability by using either the full
REMIX scheme or the REMIX scheme with unmodified, Wendland $C^{2}$ kernels.

###  3.4 Vacuum boundary treatment

We develop a method to switch the kernel gradients constructed in §3.3 to the
unmodified spherically symmetric kernel gradients in regions identified as
vacuum boundaries. We stress that this method is not applying a targeted
correction to vacuum boundaries as done by, for example, Reinhardt and Stadel
[27]. In fact, our evolved density estimate corrects density smoothing at
discontinuous free surfaces without any need for a targeted approach. Instead,
the vacuum treatment we present here is just an expansion of the form of the
linear-order reproducing kernels (§3.3) to allow them to capture free surfaces
as vacuum boundaries, a case not considered – rather than handled poorly – in
their general construction.

A region with no SPH particles is not trivially equivalent to the
representation of a vacuum. Since SPH particles are moving interpolation
points, a region not sampled by SPH particles can be seen as analogous to a
region in a grid-based code where the grid points have been removed. There is
therefore no inherent information associated with these regions that would
make them equivalent to a region with zero pressure, rather than a region to
extrapolate into. However, if a spherically symmetric kernel, normalised to
the continuum, is used to calculate pressure gradients in the equation of
motion, vacuum-like behaviour is achieved. At a free surface, a particle with
a spherically symmetric kernel will calculate pressure gradients equivalent to
those calculated if the vacuum region were built up of particles with
appropriate volumes but zero pressure141414These gradients may not be fully
equivalent in the equations of motion where the additional condition of
antisymmetry in exchange of neighbours is imposed, however, they remain
closely related..

This is not the case for the linear-order reproducing kernels described in
§3.3. Since kernels are constructed to satisfy Eqns. 19 and 20 for volumes
built up by particles only, the vacuum region is treated as a region to
extrapolate into. SPH applications typically require the treatment of a region
without SPH particles as a vacuum, or a region with negligible pressure. We
therefore switch our kernel gradient terms to gradients of unmodified kernels
at free surfaces:

$\frac{d\widetilde{\mathcal{W}}}{dr^{\gamma}}\bigg{|}_{ij}=s_{i}\frac{d\mathcal{W}}{dr^{\gamma}}\bigg{|}_{ij}+\left(1-s_{i}\right)\frac{dW}{dr^{\gamma}}\bigg{|}_{ij}\;,$
(32)

where $s$ is a function that switches from 1 in regions where no vacuum
boundary is detected, to 0 in regions near a vacuum boundary. Note that we
smoothly switch between kernel gradients rather than the kernels themselves.
This is to avoid terms with gradients of $s$. A switch that is accurate in
identifying vacuum boundary particles only will inevitably have sharp spatial
gradients, which could significantly influence the evolution of particles.
Since we do not calculate densities by Eqn. 1, we do not require the direct
calculation of the function whose derivative is given by Eqn. 32 to maintain
thermodynamic consistency.

We modify the kernel gradient terms based only on parameters of the kernel
function itself. Therefore, conceptually, we adapt the kernel function rather
than making the kernel respond to the physical system simulated. For $s$, we
use a Gaussian switch,

$s\left(h_{i}|\mathbf{B}_{i}|\right)=\begin{cases}\exp{\left[-\,\dfrac{\left(0.8-h_{i}|\mathbf{B}_{i}|\right)^{2}}{0.08}\right]}&\;\;\text{for
$h_{i}|\mathbf{B}_{i}|\geq 0.8$}\;,\\\ 1&\;\;\text{otherwise}\;,\\\
\end{cases}$ (33)

where the offset, $0.8$, and denominator, $0.08$, of the switch are chosen
empirically to identify boundary particles as those with a large
$|\mathbf{B}_{i}|$ (Eqn. 23) greater than ${\sim}1/h_{i}$. These are particles
whose kernels would have to drastically change shape to deal with large
anisotropies in the volume elements of particle neighbours. We find that using
these values allows the switch to identify particles near free surfaces
reliably without misidentifying particles in non-vacuum regions, as we show in
§4.9, where we also demonstrate the need for this vacuum boundary treatment.
In the example presented, the free surface of a Jupiter-like planet in
hydrostatic equilibrium is unstable when the vacuum boundary treatment is not
included. As well as its use in switching the kernel function, $s$ is also
used in the kernel normalisation term in the density evolution, as detailed in
§3.7.

###  3.5 Artificial viscosity

Artificial viscosity is required to capture shocks in SPH simulations, whose
constituent equations otherwise model adiabatic and dissipationless evolution
[55]. A difficulty faced by artificial viscosity constructions is over-
dissipation in regions not affected by a shock. Artificial viscosity switches,
like the Balsara switch [56],

$\mathcal{B}_{i}=\frac{\left|\nabla\cdot\mathbf{v}_{i}\right|}{\left|\nabla\cdot\mathbf{v}_{i}\right|\,+\,\left|\nabla\times\mathbf{v}_{i}\right|\,+\,0.0001\,c_{i}/h_{i}}\;,$
(34)

where $c$ is the sound speed, or higher-order switches like that of Read and
Hayfield [57] are used to switch artificial viscosity off in shearing regions.
Time-dependent viscosity parameters have also been developed [58, 59, 60] to
reduce over-dissipation.

Recently, the limiting of artificial viscosity by the use of reconstructed
velocities at particle-pair midpoints has been demonstrated to be an effective
alternative approach [19, 22, 31]. For each particle pair, two velocities are
estimated at the midpoint of the pair based on Taylor expansions from each
particle using their individual velocities and estimated velocity gradients.
The difference between these velocities is then used in the viscosity scheme
instead of the relative velocity of the particles themselves. This is the
approach taken in REMIX. We use linear reconstruction as we find further
improvements due to quadratic reconstruction to be small, as also noted by
Rosswog [22]. If the velocity field is locally linear, artificial viscosity
would effectively be switched off with linear reconstruction. For schemes that
use linear reconstruction, the viscosity in shearing regions where the
velocity field is not exactly locally linear is not negligible and will still
influence the fluid behaviour. However, this results in a helpful effect,
acting as a weak artificial diffusion of momentum that smooths particle noise
in the velocity field by guiding it towards being locally linear on the
particle scale.

Our artificial viscosity treatment is largely based on those of Frontiere et
al. [19] and Rosswog [22], with some additional, novel approaches. As detailed
below, a slope limiter is used to prevent reconstruction at discontinuities,
thereby increasing artificial viscosity where it is required for shock
capturing. However, we find that a slope limiter alone does not effectively
switch off reconstruction, because the velocity gradients used to construct it
are inherently smoothed by their calculation using a smoothing kernel.
Therefore, they do not identify sharp discontinuities well. We introduce a
Balsara switch (Eqn. 34) into the slope limiter term to switch off
reconstruction at shocks more effectively. Here we calculate
$\left|\nabla\cdot\mathbf{v}_{i}\right|$ and
$\left|\nabla\times\mathbf{v}_{i}\right|$ in the Balsara switch using the
kernel gradient term given by Eqn. 30, and also use these same gradient
estimates for the velocity gradients used in the linear reconstruction,

$\partial_{i}^{\gamma}\hat{v}^{\alpha}=\sum_{j}(v_{j}^{\alpha}-v_{i}^{\alpha})\,\partial_{i}^{\gamma}\hat{W}_{ij}\frac{m_{j}}{\rho_{j}}\;.$
(35)

The velocity reconstructed to the midpoint of a particle pair is given by

$\tilde{v}_{ij}^{\,\alpha}=v_{i}^{\alpha}+\frac{1}{2}\left(1-\mathcal{B}_{i}^{\text{SL}}\right)\Phi_{v,\;ij}\left(r_{j}^{\gamma}-r_{i}^{\gamma}\right)\partial_{i}^{\gamma}\hat{v}^{\alpha}\;,$
(36)

where $\mathcal{B}_{i}^{\text{SL}}$ is the standard Balsara switch (Eqn. 34),
and the SL (slope limiter) superscript just indicates its use in conjunction
with the slope limiter. $\Phi_{ij}$ is the van Leer slope limiter [61], given
by

$\Phi_{ij}=\begin{cases}0&\;\;\text{for $A_{ij}<0$}\;,\\\
\dfrac{4A_{ij}}{(1+A_{ij})^{2}}\exp{\left[-\left(\dfrac{\eta_{ij}^{\rm
min}-\eta_{\text{crit}}}{0.2}\right)^{2}\right]}&\;\;\text{for
$\eta_{ij}<\eta_{\text{crit}}$}\;,\\\
\dfrac{4A_{ij}}{(1+A_{ij})^{2}}&\;\;\text{otherwise}\;,\end{cases}$ (37)

where the additional Gaussian term in Eqn. 37 switches the slope limiter to 0
for particle pairs with a small separation. $\eta_{ij}^{\rm min}$ is the
smaller value of $|\boldsymbol{\eta}_{ij}|$ and $|\boldsymbol{\eta}_{ji}|$,
where $\boldsymbol{\eta}_{ij}=(\mathbf{r}_{i}-\mathbf{r}_{j})/h_{i}$ and
similarly for the exchanged particle indices. $\eta_{\rm crit}$ represents a
separation closer than one would expect from the distribution of the rest of
the particle’s neighbours. For viscosity calculations, we use the ratio of
projected velocity gradients $A_{ij}\equiv A_{v,\;ij}$ given by

$A_{v,\;ij}=\frac{\partial_{i}^{\alpha}\hat{v}^{\,\beta}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\beta}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\alpha}}{\partial_{j}^{\gamma}\hat{v}^{\,\phi}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\gamma}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\phi}}\;.$
(38)

For $\eta_{\rm crit}$ we use

$\eta_{\rm
crit}=\frac{1}{h_{i}}\left(\frac{1}{\sum_{j}W_{ij}}\right)^{1/d}\equiv\frac{1}{\eta_{\rm
kernel}}\;,$ (39)

where the equivalency is due to the definition of the smoothing length in Eqn.
4. Note that the term in brackets is an approximation of the particle volume
assuming neighbours with equal volumes.

The reconstructed velocities appear in the artificial viscosity formulation
through

$\mu_{ij}=\begin{cases}\dfrac{\tilde{\mathbf{v}}_{ij}\cdot\boldsymbol{\eta}_{ij}}{\boldsymbol{\eta}_{ij}\cdot\boldsymbol{\eta}_{ij}+\epsilon^{2}}&\;\text{for
$\tilde{\mathbf{v}}_{ij}\cdot\boldsymbol{\eta}_{ij}<0$}\;,\\\
0&\;\text{otherwise}\;,\end{cases}$ (40)

and similarly for $\mu_{ji}$ with all particle indices exchanged throughout
the calculations. $\epsilon=0.1$ is a small constant. Similarly to the
artificial viscosity of Monaghan and Gingold [55], each pressure term in the
equations of motion is modified with the addition of a pairwise viscous
pressure151515$P_{i}$ becomes $P_{i}+Q_{ij}$ and $P_{j}$ becomes
$P_{j}+Q_{ji}$.. The viscous pressure terms $Q_{ij}$ combine a linear bulk
viscosity term and a quadratic Von Neumann–Richtmyer viscosity term [62],

$Q_{ij}=\frac{1}{2}\left(a_{\mathrm{visc}}+b_{\mathrm{visc}}\mathcal{B}_{i}^{\text{visc}}\right)\rho_{i}\left(-\alpha
c_{i}\mu_{ij}+\beta\mu_{ij}^{2}\right)\;,$ (41)

and similarly for $Q_{ji}$ with all particle indices exchanged throughout the
calculations. The constants $\alpha$ and $\beta$ set the strengths of the bulk
and Von Neumann–Richtmyer terms. The constants $a_{\mathrm{visc}}$ and
$b_{\mathrm{visc}}$ set the strength of the viscosity in regions of different
flow, based on the Balsara switch, $\mathcal{B}_{i}^{\text{visc}}$.

The REMIX artificial viscosity scheme differs from those of Frontiere et al.
[19] and Rosswog [22] in some notable aspects: firstly, the Balsara switch,
$\mathcal{B}_{i}^{\text{SL}}$, is included in the slope limiter term (in Eqn.
36). This avoids reducing the artificial viscosity where it is needed, leading
to a more effective targeting of shocks. This allows us to introduce a factor
of $1/2$ in Eqn. 41 to recover equations more closely equivalent to those in
Price [11]. Otherwise, the contributions from both $Q_{ij}$ and $Q_{ji}$ would
effectively lead to this being a factor of $2$ stronger, which is to some
extent mitigated by those schemes being ineffective at switching off velocity
reconstruction in shocks. Secondly, we use $\alpha=1.5$ and $\beta=3$ as we
find that these slightly larger constants, compared with $\alpha=1$ and
$\beta=2$ as used by Frontiere et al. [19] and Rosswog [22], help to dissipate
spurious oscillations in shocks in 3D. This is consistent with typical values
used in planetary impact simulations [e.g. 27, 63]. Thirdly, we use an
additional Balsara switch directly in Eqn. 41, which, combined with the values
we use for $a_{\mathrm{visc}}=2/3$ and $b_{\mathrm{visc}}=1/3$, acts to switch
between $\alpha=1.5$ and $\beta=3$ in shocks and $\alpha=1$ and $\beta=2$ in
shearing regions. Here we make relatively conservative choices to limit the
effect of artificial viscosity in smoothing particle noise in shearing
regions, despite finding it to be a useful effect, owing to the velocity
reconstruction to particle midpoints. Our artificial viscosity scheme is
constructed to be less dissipative in shearing regions and to target shocks
more effectively than similar schemes. These choices are all discussed in more
detail in F.

###  3.6 Artificial diffusion

Artificial diffusion of internal energy, or “artificial conduction161616In
later sections, we use “artificial diffusion” to refer to cases that include
the diffusion of both density and internal energy and “artificial conduction”
where there is only diffusion of internal energy.”, is frequently used to
smooth accumulated noise in particle internal energies [64] or entropies [65],
and to improve the treatment of density discontinuities in ideal gas-only
simulations [23]. As with artificial viscosity, a targeted approach is
desirable to avoid artificial conduction playing a dominant role in the
thermodynamic evolution, instead of acting as a correction on the particle
scale [57].

In some SPH schemes, relatively strong artificial conduction is used to
address kernel smoothing at density discontinuities by smoothing particle
internal energies over kernel length scales [54]. For a single equation of
state, with no phase transitions, this leads to a smooth pressure field in the
continuous limit. However, this is not an appropriate treatment in simulations
with multiple and/or complex materials, where smooth density and internal
energy fields do not necessarily lead to smooth pressures. Additionally, even
in ideal gas-only simulations, this does not completely solve the issue, since
(1) artificial conduction becomes a less effective correction at large density
discontinuities; (2) in simulations with gravity, strong diffusion will
disturb a system’s hydrostatic equilibrium; (3) artificial conduction does not
attempt to address the source of kernel smoothing error directly, instead it
alters the physical system itself to one without discontinuities.

In simulations that use an evolved density estimate (Eqn. 6), a similar
artificial diffusion term can be used in the evolution of densities, for
example, in the $\delta$-SPH formulation, used predominantly for engineering
applications [52, 66, 67].

In REMIX, we include artificial diffusion of specific internal energy and of
density, both to improve the treatment of shocks and to smooth accumulated
noise on the particle scale, using reconstruction to particle midpoints [22,
52]. Similarly to the phase dependence in the diffusion schemes of Sun et al.
[67] and Pearl et al. [31], we only allow diffusion between particles of the
same material type. Without this distinction, artificial diffusion of internal
energy between different materials would cause unphysical evolution, since
smoothing would be based on internal energy and not temperature. Diffusing
density between different materials would lead to density discontinuities at
material interfaces returning to a similar, smoothed state as in simulations
with smoothing error in the density estimate.

The diffusion terms in the equations of motion take the form

$\displaystyle\left(\frac{du_{i}}{dt}\right)_{\text{difn}}$
$\displaystyle=\sum\limits_{j}\kappa_{ij}\left(a_{u}+b_{u}\mathcal{B}_{ij}^{\text{difn}}\right)v_{\text{sig},ij}\>(\tilde{u}_{j}-\tilde{u}_{i})\frac{m_{j}}{\rho_{ij}}\frac{1}{2}\left|\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ji}\right|\;,$
(42) $\displaystyle\left(\frac{d\rho_{i}}{dt}\right)_{\text{difn}}$
$\displaystyle=\sum\limits_{j}\kappa_{ij}\left(a_{\rho}+b_{\rho}\mathcal{B}_{ij}^{\text{difn}}\right)v_{\text{sig},ij}\>(\tilde{\rho}_{j}-\tilde{\rho}_{i})\frac{\rho_{i}}{\rho_{j}}\frac{m_{j}}{\rho_{ij}}\frac{1}{2}\left|\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ji}\right|\;,$
(43)

where $\kappa_{ij}=1$ for particles of the same material and $\kappa_{ij}=0$
otherwise. The average Balsara switch for each particle pair is used,
$\mathcal{B}_{ij}^{\text{difn}}=(\mathcal{B}_{i}+\mathcal{B}_{j})/2$, for
conservation. We take the signal velocity to be
$v_{\text{sig},\;ij}=\left|\tilde{\mathbf{v}}_{i}-\tilde{\mathbf{v}}_{j}\right|$
and do not draw any distinctions between simulations with and without gravity
(unlike some previous works [54, 22]), since we aim to validate the full REMIX
formulation independently of specific simulation properties. The parameters
$a_{u}$ and $a_{\rho}$ set the strength of the artificial conduction in
shearing regions (where $\mathcal{B}_{ij}^{\text{difn}}\rightarrow 0$) and are
increased to $a_{u}+b_{u}$ and $a_{\rho}+b_{\rho}$ in shocks. In shearing
regions we choose to have low amounts of diffusion to avoid this strongly
influencing thermodynamic evolution, and to allow for persisting and emergent
discontinuities. We therefore use $a_{u}=a_{\rho}=0.05$, similarly to Rosswog
[22]. In the presence of shocks we find that we need a much larger amount of
diffusion to prevent spikes in density and internal energy, and so we use
$b_{u}=b_{\rho}=0.95$. We motivate and test the sensitivity of these choices
in F. The volume elements in Eqn. 42 are chosen to conserve energy. In Eqn.
43, they include an additional ratio of densities, to conserve volume in each
pairwise interaction171717Substituting
$\frac{d\rho_{i}}{dt}=-\frac{\rho_{i}^{2}}{m_{i}}\frac{dV_{i}}{dt}$ and
solving for $\frac{dV_{i}}{dt}$ gives an equation antisymmetric in the
exchange of particles.. Although conserving volume in a pairwise interaction
between particles is not strictly necessary, we find that it improves the
treatment of the density diffusion in shocks.

When calculating the artificial diffusion terms, internal energies and
densities are reconstructed to particle midpoints similarly to the velocities
in the artificial viscosity scheme via

$\tilde{u}_{i}=u_{i}+\frac{1}{2}\Phi_{u,\;ij}\left(r_{j}^{\gamma}-r_{i}^{\gamma}\right)\partial_{\kappa,\,i}^{\gamma}\,\hat{u}\;,$
(44)

$\tilde{\rho}_{i}=\rho_{i}+\frac{1}{2}\Phi_{\rho,\;ij}\left(r_{j}^{\gamma}-r_{i}^{\gamma}\right)\partial_{\kappa,\,i}^{\gamma}\,\hat{\rho}\;.$
(45)

The derivatives are calculated using only particles of the same material
species as

$\partial_{\kappa,\,i}^{\gamma}\,\hat{u}=\sum_{j}\kappa_{ij}\,(u_{j}-u_{i})\,\partial_{i}^{\gamma}\hat{W}_{ij}\frac{m_{j}}{\rho_{j}}\;,$
(46)

$\partial_{\kappa,\,i}^{\gamma}\,\hat{\rho}=\sum_{j}\kappa_{ij}\,(\rho_{j}-\rho_{i})\,\partial_{i}^{\gamma}\hat{W}_{ij}\frac{m_{j}}{\rho_{j}}\;.$
(47)

The material dependence of these gradients helps to preserve real
discontinuities at material boundaries.

The slope limiter is calculated in the same way as for the viscosity, Eqn. 37,
but with $A_{ij}=A_{u,\;ij}$ and $A_{ij}=A_{\rho,\;ij}$ given by

$A_{u,\;ij}=\frac{\partial_{\kappa,\,i}^{\alpha}\,\hat{u}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\alpha}}{\partial_{\kappa,\,j}^{\beta}\,\hat{u}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\beta}}\;,$
(48)

$A_{\rho,\;ij}=\frac{\partial_{\kappa,\,i}^{\alpha}\,\hat{\rho}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\alpha}}{\partial_{\kappa,\,j}^{\beta}\,\hat{\rho}(\mathbf{r}_{j}-\mathbf{r}_{i})^{\beta}}\;.$
(49)

Although our diffusion scheme technically includes material dependence, this
is not a correction targeted at material boundaries, nor with any dependence
on the actual EoS. Rather, we actively turn off these parts of our method for
particles of different species. Our artificial diffusion scheme is used to
improve the treatment of shocks, and to weakly smooth accumulated noise. It is
not used to address surface tension-like effects that prevent mixing and
instability growth at density discontinuities, even in our ideal gas-only
simulations.

###  3.7 Normalising term

We add a normalising term to the density evolution equation. This aims to
evolve densities to reflect the distribution of mass in nearby particles,
particularly in regions where particle volume elements systematically fail to
satisfy the normalisation of the kernel. Since error accumulates in the
evolution of densities based on timescales set by the divergence operator used
in the equations of motion, we set the normalising term to act over timescales
determined by the motion of particles. This also allows particles to move in
response to changes in density caused by the normalising term.

Particle volume elements should approximately satisfy $\sum_{j}W_{ij}V_{j}=1$
(Eqn. 12) for a normalised kernel function. However, this condition will not
be satisfied either if particle densities are poor estimates of the underlying
field or if particle masses do not appropriately represent the mass
distribution of the fluid. Our methods inherently conserve mass, as particle
masses do not evolve during the simulation, and are fully Lagrangian.
Therefore, we choose to maintain the simplicity and computational stability of
this construction, and address discrepancies in volume elements through
particle densities rather than through particle masses or their distribution.
We do this by including an additional term in the density evolution, which we
refer to as the “normalising term”, that evolves densities towards a set of
volume elements that aim to appropriately build up the continuous simulation
volume. We note that the role of this term is not to obtain volume elements
that exactly satisfy normalisation for all particle kernels at any given time,
but rather to keep volume elements loosely tied to kernel normalisation and to
address regions with systematic discrepancies.

To construct our normalising term, we consider the zeroth geometric moment of
the unmodified kernel,

$m_{0,\,i}=\sum_{j}W_{ij}V_{j}\;,$ (50)

where $m_{0,\,i}=1$ if the kernel $W_{ij}$ is normalised over the volume
elements $V_{j}=m_{j}/\rho_{j}$. For a single particle $i$, we could trivially
satisfy this condition by modifying the density of the particle and all its
neighbours, $j$, by replacing $\rho_{j}$ with $m_{0,\,i}\,\rho_{j}$. However,
this does not imply that $m_{0,\,j}=1$ for all $j$, which will all have
different $m_{0,\,j}$ and different sets of neighbours. But if there are
systematic discrepancies in $m_{0}$ for many neighbouring particles, then
modifying densities in a similar way for all these particles will move them
closer to $m_{0,\,j}=1$. For instance, consider a region where particles have
systematically too low density, leading to a local trend of $m_{0,\,j}>1$.
Here, increasing the densities will evolve these particles towards
$m_{0,\,j}=1$ and towards a density field that better represents the local
mass distribution. In practice, we capture this behaviour with a smooth
evolution in time. Unlike in the initial naïve example of modifying the
densities of all $j$ to satisfy $m_{0,\,i}=1$ for $i$ only, we evolve the
density of $i$ only, based on its own $m_{0,\,i}$. This reduces the risk of
emergent chaotic behaviour and still captures the desired behaviour in regions
of systematic trends away from kernel normalisation. The normalising term in
the density evolution equation takes the form

$\left(\frac{d\rho_{i}}{dt}\right)_{\text{norm}}=\alpha_{\text{norm}}\,s_{i}\,(m_{0,\,i}-1)\,\rho_{i}\sum\limits_{j}v_{\text{norm},\;ij}\>\frac{m_{j}}{\rho_{ij}}\frac{1}{2}\left|\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ij}-\frac{d\widetilde{\mathcal{W}}}{d\mathbf{r}}\bigg{|}_{ji}\right|\;.$
(51)

where $\alpha_{\text{norm}}=1$ is a constant and
$v_{\text{norm},\;ij}=\left|\mathbf{v}_{i}-\mathbf{v}_{j}\right|$ is the
effective signal velocity. Eqn. 51 aims to contribute to a weak evolution of
$\rho_{i}$ towards $m_{0,\,i}\,\rho_{i}$. We include the vacuum switch, $s$,
described in §3.4, since the kernel should not be normalised by particle
volume elements at vacuum boundaries181818At vacuum boundaries, one would
instead expect $m_{0,\,i}\approx 1/2$.. Here, we use the same volume elements
and kernel gradient terms as are used in the diffusion of internal energy
(Eqn. 42), despite not being motivated by conservation in this term, since it
does not represent the exchange of a quantity between particles. We use these
so that the timescale of the normalising evolution is based on terms in the
sum that are equal for both particles in each pairwise interaction. This
prevents individual particles dominating in the corrective evolution. Using a
timescale that depends on particle motion $v_{\text{norm},\;ij}$ rather than,
for example the sound speed, allows particles to react and move in response to
changes in density caused by the normalising term. We find that using an
effective signal velocity that depends on the sound speed, even with a small
multiplication factor, can lead to spatial oscillations in density, because
densities change to attempt to satisfy normalisation faster than particles can
respond to these changes.

We show the effect of this term in simulations in §4.3.2 and §4.9. In
particular, we show that without this term, an example Jupiter-like planet in
hydrostatic equilibrium will develop numerical instabilities as particles with
low evolved densities, but are in regions of high particle number density,
move from the planet’s surface towards its core (§4.9). In less extreme cases,
the normalising term does not have a significant effect on hydrodynamics,
although it does lead to particle densities that are generally closer to
satisfying kernel normalisation, $m_{0,\,i}=1$.

##  4 Hydrodynamic Tests

In this section, we validate REMIX in simulations to test its ability to
capture physically realistic fluid behaviour. The primary tests are performed
with particles of equal mass across the simulation, though we also include a
subset of additional simulations for direct comparisons with past work, where
particles are placed onto a regular grid but have different masses. We refer
to these two cases as “equal mass” and “equal spacing” throughout the
following sections. The choice to focus on simulations with equal mass is made
to validate our methods for science applications where particle densities and
configurations can evolve significantly from their initial states, so particle
masses cannot be easily chosen in the initial configuration to address errors.
All simulations are performed in 3D, to account for effects that do not change
predictably when increasing the number of dimensions, such as due to more
freedom in particle configurations, or the change in scaling between neighbour
number and length scale of particle interactions191919A 2D simulation will
have a lower number of neighbours for a given smoothing length than the
equivalent 3D simulation. Increasing $\eta_{\rm kernel}$ to compensate for
this would lead to kernel smoothing over a larger length scale.. Additionally,
in figures showing simulation snapshots, we deliberately plot individual
particles rather than the smooth, reconstructed fields shown in some works. It
is particularly important to visualise small-scale behaviour of simulations
that aim to improve the treatment of density discontinuities where the effects
that suppress mixing act on the particle scale.

We present results for the following hydrodynamic test scenarios:

  * 1.

the square test (§4.1), where we investigate the treatment of density
discontinuities in static equilibrium;

  * 2.

the Sod shock tube (§4.2), where we investigate the treatment of shocks;

  * 3.

the Kelvin–Helmholtz instability both with an ideal gas EoS (§4.3) and between
different, stiff materials set up to be representative of iron & rock material
boundaries in an Earth-like planet (§4.4);

  * 4.

the Rayleigh–Taylor instability, also both with an ideal gas EoS (§4.5) and
with iron & rock (§4.6);

  * 5.

the blob test (§4.7), with which we investigate the onset of turbulence due to
unseeded instabilities in both subsonic and supersonic regimes;

  * 6.

the Evrard collapse (§4.8), which is used to test the interaction of our
hydrodynamic treatment with gravity and shocks;

  * 7.

and finally, planets in hydrostatic equilibrium (§4.9), which we consider as a
test scenario that combines gravity, complex-material boundaries, and a vacuum
boundary.

The initial conditions needed to perform these tests are included as examples
in the open-source Swift code.

We include comparisons with simulations carried out both using a traditional
SPH formulation (“tSPH”) and a traditional formulation that includes
artificial conduction of internal energy (“tSPH $+$ cond.”), with full details
in C. These are used to demonstrate the motivation and need for many of the
improvements in REMIX. We note that in most ideal gas tests, we follow the
convention of past work and leave quantities unitless.

###  4.1 Square test

The “square test” is used to investigate spurious surface tension-like effects
from sharp discontinuities in a system that should be in static equilibrium
[29]. Here we test both an equal spacing scenario, i.e., with different
particle masses in the two regions, and an equal mass scenario. The
significant contributions from both smoothing and discretisation error
(§2.1.3, §2.1.4) at the density discontinuity makes the equal mass test
particularly challenging for SPH.

A square (or cube) of fluid of higher density is initiated in pressure
equilibrium with the surrounding region of low density fluid. Since the fluid
experiences no gradients in pressure, other than those created by numerical
errors, the shape of the square should not distort with time. In tSPH
simulations, spurious surface tension-like effects at the density
discontinuity leads to non-zero accelerations and a deformation of the square
[5]. Typically, this test is carried out in 2D however, here we simulate a
more challenging 3D cube with its effectively “sharper” higher-dimension
corners, similarly to Rosswog [22].

First, for the equal spacing scenario, we use initial conditions set up to
match those of Rosswog [22]. $40^{3}$ particles are placed in a simple cubic
lattice with spacing $1/40$ between adjacent particles. The simulation box is
periodic and has length 1 in each of the $x$, $y$, $z$ directions. Masses are
chosen such that $m_{i}=\rho(\mathbf{r}_{i})/40^{3}$, with densities $\rho=4$
in the region $-0.25<x,y,z<0.25$ and $\rho=1$ otherwise. An ideal gas EoS with
$\gamma=5/3$ is used for all particles. Initial internal energies are set to
give a uniform pressure202020We note the use of the unsmoothed density $\rho$
rather than the smoothed $\langle\rho\rangle$ used to set the internal
energies of the initial conditions. Therefore tSPH simulations are not
initialised in pressure equilibrium, due to smoothing error in the density
estimate. of $P(\rho,u)=2.5$.

Figure 1: Central cross-sections from 3D “square test” simulations. Snapshots
of the cube are shown at three times from simulations with equal initial
particle spacing (top full-row) and equal particle mass across the simulation
(bottom full-row). Simulations were carried out using three SPH formulations:
tSPH, tSPH with artificial conduction, and REMIX. Individual particles are
plotted at their positions in the $x$–$y$ plane on a grey background, and
coloured by their density. Particles in a slice of thickness 0.1 are plotted,
so the grey background is visible in regions that have maintained their grid
alignment in $z$ from the initial conditions.

The evolution of the equal spacing square test, carried out using each of
tSPH; tSPH with artificial conduction; and REMIX, is shown in the top panels
of Fig. 1. In the equal spacing scenario, the major contribution to spurious
surface tension is due to the smoothing of the density field. The contribution
of discretisation error is small, due to the well-ordered particle
distribution and use of a relatively high-order kernel. With tSPH, the cube
quickly deforms to a more stable, spherical shape, as illustrated by the
upper, top-left panels of Fig. 1. Artificial conduction reduces the effects of
smoothing error and so a square shape persists for longer, although the
sharpness of the discontinuity is not maintained (Fig. 1 lower, top-left
panels). With REMIX, particle motion is negligible, relative to the particle
separation, and the cube retains its shape (Fig. 1 top-right panels). This is
in large part due to the the use of the evolved density estimate, which
prevents density smoothing – and therefore spurious pressures – at the
discontinuities. Our choice of the free functions in the equations of motion
and kernel construction also help in reducing discretisation error to achieve
these results.

Next, we consider the more challenging case for SPH: the use of equal mass
particles, which leads to particles set up in considerably different grid-
spacings interacting at the density discontinuity. Particles in the low
density region are placed in the same configuration as in the equal spacing
scenario. Then, instead of increasing particle masses in the high density
region, the particle spacing is decreased and masses are kept the same as in
the low density region. To satisfy these conditions while closely matching the
density ratio in the equal spacing test, the high density region is given a
grid-spacing of a factor 0.625 finer than the grid-spacing of the low density
region. This corresponds to a density of 4.096. The new spacing of high
density particles is chosen such that the layers of particles on either side
of discontinuities are separated by the mean of the two grid-spacings, for all
cube faces.

The evolution of this square test with equal mass particles is shown in the
bottom panels of Fig. 1. There is now a large contribution of both smoothing
and discretisation error in both of the traditional SPH formalisms. As such,
the cube quickly deforms, even with conduction acting to reduce smoothing
error. In the REMIX formulation, some minor deformation can be observed over
these timescales. However, the general shape is maintained (Fig. 1 bottom-
right panels). We note that although past work typically shows 2D square test
evolution over longer timescales than those of our plotted snapshots, our
plots show times later than the comparable 3D tests in Rosswog [22], beyond
the time at which their equal spacing cubes have deformed. Reducing the
effects of artificial surface tension requires all of (1) a density estimate
that does not smooth density discontinuities, (2) our choice of equations of
motion, and (3) improved gradient estimates. In the REMIX simulation,
artificial diffusion is not the dominant source of correction, as
discontinuities in both density and internal energy remain sharp.

If the linear-order reproducing kernels are used in the equations of motion
without the antisymmetrisation, which is needed to enforce conservation, the
square will remain undisturbed over much longer timescales, even in the equal
mass case. The difference in outcome between using the conservative,
antisymmetric construction and the exactly linear reproducing construction is
sensitive to the kernel function used to construct the linear reproducing
kernel. Therefore, reducing the additional error introduced in
antisymmetrisation becomes an important consideration when choosing the form
of the kernel from which the linear-order reproducing kernels are constructed.
This can be seen in E, where we present sensitivities in these results to
different elements of the REMIX construction.

###  4.2 Sod shock tube

Figure 2: 3D Sod shock tube at time $t=0.2$ simulated using tSPH and REMIX.
Plots show velocity in the $x$-direction, $v_{x}$, density, $\rho$, specific
internal energy, $u$, and pressure, $P$, of individual particles plotted
against their $x$-position for tSPH (a–d) and REMIX (e–h) respectively. The
red, dashed line shows a reference solution, solved for directly by using a
Riemann solver. All particles are plotted.

The “Sod shock tube” [68] is used to assess the shock capturing capabilities
of our hydrodynamic scheme. This is a classic Riemann problem with a known
analytic solution. Since the inclusion of artificial viscosity and diffusion
are necessary to deal with shocks in the REMIX scheme, we also use this test
to motivate choices made in the artificial viscosity and diffusion
formulations, as detailed in F. The choices made in the viscosity scheme
relating to this test focus on reducing ringing oscillations behind the shock.
The diffusion scheme focuses on reducing the size of spikes in density and
internal energy at the discontinuity.

Ideal gas, $\gamma=5/3$, particles of equal mass are placed in a periodic 3D
domain with size 2 in each of $x,y,z$ directions, centred at (0, 0, 0). We use
two glass configurations, scaled appropriately for the two regions of
different initial density: $\rho_{1}=1$ in the region $x<0$ and $\rho_{2}=1/8$
in the region $x>0$. Initial internal energies are set such that $P_{1}=1$ and
$P_{2}=0.1$. Simulations have a total of 589,824 equal mass particles.

The particle velocities in the $x$-direction, densities, internal energies,
and pressures at a time $t=0.2$ are shown in Fig. 2. The shock is captured
well with REMIX, and the particles align with the reference solution. Noise in
particle velocities is reduced compared with tSPH. The size of the spike in
internal energy is also reduced. The pressure blip could be further smoothed
by increasing the strength of our artificial diffusion scheme, through choices
of the $a$ and $b$ factors. However, we choose to take a conservative approach
to artificial diffusion to avoid deviating far from the thermodynamically
consistent core equations of motion.

###  4.3 Kelvin–Helmholtz instability – ideal gas

The Kelvin–Helmholtz instability (KHI) is the first test we use to investigate
the treatment of mixing and dynamic instability growth in our simulations. The
KHI arises at shearing interfaces in fluids [69]. Perturbations at the
interface grow to form vortices that act to cascade energy to shorter length
scales. As such, the KHI plays a significant role in the onset of turbulence
in physical systems. Capturing the growth of the KHI has therefore been widely
adopted as a benchmarking test to assess a numerical method’s ability to
simulate turbulence-driven mixing, as well as mixing on the particle scale.
However, unlike the other tests above, an analytical solution does not exist
for the KHI.

Here we first consider the growth of these instabilities at shearing density
contrasts in an ideal gas. All simulations presented are carried out in 3D,
with a thin $z$ direction depth relative to the other dimensions, similarly to
Hopkins [42], Read et al. [12], and Rosswog [22]. We focus primarily on cases
with a sharp density discontinuity and equal mass particles. This is in
contrast with an alternative setup with which we directly compare our results
with a reference solution [70], where we consider an initially smoothed
discontinuity and equal particle spacing. Although the use of this second form
of initial conditions with smooth initial densities and velocities is
motivated by the existence of a converged solution, these choices change the
physical system to one with inherently less smoothing and discretisation
error, which are the main effects of interest that normally suppress
instability growth in SPH simulations. These smooth initial conditions
therefore do not give the full picture of an SPH scheme’s ability to capture
KHI growth at sharp density discontinuities, where these sources of error can
play a dominant role. This is particularly important at material boundaries,
where smoothing the density discontinuity between different materials may lead
to particles of both materials occupying extreme regions of their EoS phase
space, so considering deliberately smoothed, equilibrium initial conditions
would not be representative of a physical system.

Traditional formulations of SPH struggle to capture the KHI [14], with the
growth of the instability being strongly suppressed. In particular, for
shearing density discontinuities, smoothing in the density estimate leads to
surface tension-like effects that act to artificially stabilise the interface.
Additionally, for simulations where SPH particles in both density regions have
equal mass, or configurations that give similarly anisotropic local particle
distributions at the interface, leading-order error in the momentum equation
will also contribute significantly to this spurious surface tension-like
effect. Not only do these effects act to suppress mixing by hampering the
large-scale evolution of naturally arising instabilities that should act to
drive mixing, but they will also impede particles crossing interfaces, thereby
suppressing mixing both indirectly and directly.

The growth of a mode of wavelength $\lambda$ is characterised by the timescale
[22, 23]

$\tau_{\rm
KH}=\frac{(\rho_{1}+\rho_{2})~{}\lambda}{\sqrt{\rho_{1}\rho_{2}}\,|v_{1}-v_{2}|}\;,$
(52)

where $\rho_{1}$ and $\rho_{2}$ are the densities in regions separated by the
shearing interface and $|v_{1}-v_{2}|$ is their relative speed. We use this
parameterisation so that comparisons can be drawn at the same $\tau_{\rm KH}$
between simulations with different initial conditions, since we consider KHIs
with both smoothed and sharp interfaces, for different density ratios, and
between different materials. We note that initial conditions with and without
initial smoothing of fields at the interface are physically different systems,
so we do not expect converged results between the two.

In the absence of stabilising influences such as physical surface tension or
gravity, a shearing discontinuity is unstable to perturbation modes of all
wavenumbers [69]. In a realistic system satisfying these conditions,
instability will always be triggered, as even the smallest local
inhomogeneities will seed mode growth. Similarly, in a simulation, numerical
error will inevitably trigger instability at shearing discontinuities. The
wavenumbers of error-seeded modes are sensitive not only to the numerical
methods used and the construction of initial conditions, but also to the
resolution of the simulation: a higher resolution simulation will be able to
resolve the excitation of a wider range of mode wavelengths [71]. The growth
of KHIs at sharp discontinuities can therefore not be used reliably for
convergence studies.

McNally et al. [70] and Robertson et al. [71] construct KHI initial conditions
with smooth initial velocities and densities across the shearing interface.
They show that the inclusion of a well-resolved transition region acts to
stabilise the system, suppressing modes other than those deliberately seeded
in the initial conditions. They demonstrate convergence and present a well-
posed method to benchmark the early evolution of KHI simulations. In §4.3.1
below, we present REMIX simulations using the initial conditions of McNally et
al. [70], including quantitative comparisons of mode growth with their
converged reference solution. In §4.3.2, we present KHI simulations with sharp
discontinuities in density and velocity across the interface. Although we
cannot make quantitative comparisons of this more challenging case with
converged reference solutions, useful comparisons can still be drawn between
simulations and the expected qualitative behaviour of the instability, with a
motivation of reducing the clear suppression of the KHI observed when using
traditional SPH. We additionally use equal mass particles across the
simulation, making this setup particularly challenging for SPH schemes, but
more applicable to most science applications. In §4.3.3 we present KHI
simulations with a larger density ratio, a discontinuous interface, and equal
mass particles. This system is even more challenging again for SPH schemes:
both smoothing and discretisation errors are increased here due to the larger
density-smoothing effects and the even more extreme local anisotropy in
particle distribution at the interface. After considering these ideal gas
scenarios, we present KHI simulations at interfaces between dissimilar, stiff
materials in §4.4.

####  4.3.1 KHI with smooth initial conditions

Figure 3: Ideal gas Kelvin–Helmholtz instabilities with smoothed initial
density and velocity profiles. Columns correspond to simulations of different
resolutions, with the top row showing results from simulations using tSPH and
the bottom row from simulations using REMIX. We plot snapshots at
$t=2~{}\tau_{\rm KH}$ from 3D, ideal gas KHI simulations. The density ratio
between the two regions is 1:2. Here particles of different mass are used to
match consistent initial particle spacing and volume across the simulation.
Individual particles are plotted on a grey background and coloured by their
density. Particles at all $z$ are plotted, so the grey background is visible
in regions that have maintained their grid alignment in $z$ from the initial
conditions.
Figure 4: Evolution of mode amplitude, $M$, in Kelvin–Helmholtz simulations
with smoothed initial density and velocity profiles. We plot mode growth for
simulations with three different resolutions ($N=128$, $256$, $512$) using
both tSPH (grey, dashed) and REMIX (blue, dashed). Mode amplitude is
normalised to the initial amplitude of the excited mode, $M_{0}$, and time is
normalised to the characteristic timescale of KHI growth, $\tau_{\rm KH}$. The
reference solution (red, solid) corresponds to the 4096 $\times$ 4096
simulation of McNally et al. [70] using the Pencil code.

McNally et al. [70] present converged, high-resolution simulations of the
early linear growth of the KHI. They use initial conditions with smooth
initial velocity and density fields across the shearing interface. Similarly
to Rosswog [22], here we use these smooth initial conditions, adapted to 3D,
and use the mode growth of the reference solution of McNally et al. [70] to
quantitatively assess the accuracy of our numerical methods.

Particles are initialised in a 3D cubic lattice in a periodic box with
$N\times N\times 18$ particles in $x,y,z$ directions (i.e. a thin slice in the
$z$ direction relative to $x$ and $y$). We run simulations with resolutions
$N=128$, $256$, $512$. Spatial dimensions are normalised to the size of the
simulation box length in the $x$ and $y$ directions. A low density region of
$\rho_{1}=1$ shears against a high density region of $\rho_{2}=2$ with speeds
in the $x$ direction of $v_{1}=-0.5$ and $v_{2}=0.5$ such that the relative
velocity is $|v_{1}-v_{2}|=1$. The regions are layered in $y$ and have
relative velocities in $x$. However, both density and shearing velocity are
smoothed at the shearing interface such that initial densities are given by

$\rho(y)=\begin{cases}\rho_{1}-\rho_{m}e^{(y-0.25)/\Delta}&\text{for $0.00\leq
y<0.25$}\,,\\\ \rho_{2}+\rho_{m}e^{(0.25-y)/\Delta}&\text{for $0.25\leq
y<0.50$}\,,\\\ \rho_{2}+\rho_{m}e^{(y-0.75)/\Delta}&\text{for $0.50\leq
y<0.75$}\,,\\\ \rho_{1}-\rho_{m}e^{(0.75-y)/\Delta}&\text{for $0.75\leq
y<1.00$}\;,\end{cases}$ (53)

and initial velocities in the $x$ direction are given by

$v_{x}(y)=\begin{cases}v_{1}-v_{m}e^{(y-0.25)/\Delta}&\text{for $0.00\leq
y<0.25$}\,,\\\ v_{2}+v_{m}e^{(0.25-y)/\Delta}&\text{for $0.25\leq
y<0.50$}\,,\\\ v_{2}+v_{m}e^{(y-0.75)/\Delta}&\text{for $0.50\leq
y<0.75$}\,,\\\ v_{1}-v_{m}e^{(0.75-y)/\Delta}&\text{for $0.75\leq
y<1.00$}\;.\end{cases}$ (54)

Here $\rho_{m}=(\rho_{1}-\rho_{2})/2$, $v_{m}=(v_{1}-v_{2})/2$, and
$\Delta=0.025$. Since particle positions are initialised in a single cubic
lattice, particle masses are set by $m_{i}=\rho(y_{i})/N^{3}$. Particle
internal energies are set to give a pressure of $P(\rho,u)=2.5$ across the
simulation for an ideal gas with $\gamma=5/3$. A small velocity perturbation,
$v_{y}=0.01\sin{(2\pi x/\lambda)}$, is added in the $y$ direction with
wavelength $\lambda=0.5$, to seed the primary instability.

The simulated KHI with these initial conditions is shown in Fig. 3. We plot
particle densities at particle positions for simulations of resolution
$N=128$, $256$, $512$. Top row plots correspond to tSPH and bottom to REMIX.
All snapshots are shown at simulation time $t=2~{}\tau_{\rm KH}$. Traditional
SPH struggles to capture this instability at low resolutions. In REMIX
simulations the seeded mode is not suppressed and grows at a close to
resolution-independent rate. We find, however, that at later times secondary
modes will eventually grow and disturb the evolution of the primary mode, so
we do not observe strict convergence over long timescales. For an SPH scheme
aiming to model an inviscid fluid with realistic turbulence-driven mixing, a
compromise on this is difficult to avoid.

The evolution of the amplitude of the seeded mode is shown in Fig. 4, for
these simulations. This quantity, $M$, is calculated from Eqns. 10–13 of
McNally et al. [70]. We normalise the mode amplitude to $M_{0}\equiv M(t=0)$
to allow for more direct comparisons between simulations with different
initial conditions, presented later. The reference solution is from the high-
resolution $4096^{2}$ cell KHI simulation performed by McNally et al. [70]
using the Eulerian mesh, finite-difference code Pencil. The mode growth of the
tSPH simulations fall significantly short of the reference solution. This
result is consistent with the SPH simulations used for comparisons by McNally
et al. [70]. In contrast, the mid- and high-resolution REMIX simulations
closely match the reference solution, and even the lowest resolution REMIX
simulation is considerably closer to the reference solution than the highest
resolution tSPH simulation.

####  4.3.2 KHI with discontinuous initial conditions

Figure 5: Growth of 3D Kelvin–Helmholtz instabilities in the more challenging
case of discontinuous initial density and velocity profiles and equal mass
particles. Columns show snapshots at different times, with the top rows
showing results from simulations using tSPH – without and with artificial
conduction – and the bottom row from simulations using REMIX. These
simulations are both relatively low resolution, with $N_{1}=128$. The density
ratio between the two regions is 1:1.91. Figure 6: Evolution of
Kelvin–Helmholtz instability mode amplitude, $M$, in sharp-discontinuity,
equal mass KHI simulations with (a) an ideal gas EoS (§4.3.2) and (b) between
dissimilar, stiff EoS (§4.4). Mode growth in simulations using both tSPH
(grey, dashed) and REMIX (blue, dashed) are shown for three different
resolutions ($N_{1}=128$, $256$, $512$). The mode amplitude is normalised to
the initial amplitude of the excited mode, $M_{0}$, and time is normalised to
the characteristic timescale of KHI growth, $\tau_{\rm KH}$, which are both
different for each case.

Next we consider KHI growth at an interface that is discontinuous in density
and velocity. A shearing discontinuous interface is unstable to modes of all
wavelengths, so noise- or error-seeded secondary modes will inevitably lead to
a turbulent system and preclude numerical convergence. Although no converged
reference solution exists for this problem, we use this system to
qualitatively demonstrate the suppression in tSPH of both instability growth
at, and mixing across, density discontinuities, and the effectiveness of REMIX
in alleviating these issues. We deliberately constrain our analysis to low-
resolution simulations, where the primary, intentionally seeded mode remains
relatively undisturbed by secondary modes during the early growth of the
instability (as discussed further in G).

Similarly to in §4.3.1, we consider shearing between a low density region of
$\rho_{1}=1$ and a high density region of $\rho_{2}\approx 2$ with relative
speeds of $v_{1}=-0.5$ and $v_{2}=0.5$. Here we initialise particles with a
sharp discontinuity in both density and shearing velocity. The low density
region is set up in the same cubic lattice as in the smoothed simulations of
the previous section. We use particles of equal mass across the simulation. We
refer to the resolution of these simulations by the effective resolution of
the low density region: $N_{1}=128$, $256$, $512$; if the box were filled with
a cubic lattice of the low density material only, then this lattice would
consist of $N_{1}\times N_{1}\times 18$ particles in $x,y,z$ directions.
Particles in the high density region are arranged in a cubic lattice of
shorter spacing. A cubic lattice corresponding to a density $\rho_{2}=2$ is
adjusted to allow a continuous grid in the $x$ dimension of the periodic
box212121We also enforce that the effective resolution in this region,
$N_{2}$, is divisible by $4$ (the number of vortices formed by the evolution
of the seeded mode) to avoid the possibility of asymmetric evolution of
individual vortices triggering an early onset of secondary modes. In practice
this has no noticeable effect here, but similar considerations do matter for
the Rayleigh–Taylor instability simulations we examine in §4.5 and §4.6.. The
spacing of particles in $z$ is slightly adjusted away from a perfectly cubic
lattice such that particle spacing in this dimension is also continuous across
the boundary of the box. The regions are shifted in the $y$ direction such
that the layers of particles across the interface from each other, directly
adjacent to, and parallel with, the discontinuity are separated by the mean of
the two grid-spacings. The size of the simulation box is adjusted in the $y$
direction to compensate for this and to maintain continuity across boundaries
of the periodic box. The density $\rho_{2}$ is recalculated based on these
grid modifications and the use of equal mass particles. To satisfy these
conditions, in the high density region we use $\rho_{2}=1.91$ in a lattice
with, for example, $N_{2}=160$ and $22$ particles in the $z$ direction for
$N_{1}=128$. Initial internal energies are calculated such that particles have
a uniform initial pressure222222We note here that the density used in these
initial conditions are unsmoothed, so the tSPH simulations will not be in
pressure equilibrium due to their smoothing of the densities at the
discontinuities. of $P(\rho,u)=2.5$ by the ideal gas equation with
$\gamma=5/3$. We seed a small velocity perturbation, $v_{y}=0.01\sin{(2\pi
x/\lambda)}$, in the $y$ direction with $\lambda=0.5$.

In Fig. 5 we show the growth of these KHIs simulated using tSPH, tSPH with
artificial conduction, and REMIX. We plot individual particles, coloured by
their densities, at three times through the evolution of the instability. In
the tSPH simulations, surface tension-like effects act both to suppress the
growth of the instability and to prevent mixing of particles across the
interface. As noted by Agertz et al. [14], particles form ordered bands with
large gaps at the interface, which act as barriers to mixing. Artificial
conduction helps to enable some mixing on the particle scale, allowing the
boundary to become diffuse with time. However, the evolution of the
instability is slow, as can be seen when comparing with the similar sharp-
interface KHI simulations of Hopkins [37], performed with their improved
methods, at comparable scaled times (their Fig. 21). While we note that
differences in initial conditions construction means that we cannot make
direct comparisons, the growth of the instability in both traditional cases is
clearly too slow. The REMIX simulation shows a clear improvement: not only do
the characteristic vortices of the KHI form without impedance by surface
tension-like effects, but interfaces are maintained as sharp discontinuities
as the system evolves. Particles do not align themselves in bands separated by
gaps that would prevent mixing across the discontinuity.

The mode amplitude growth of these KHIs and equivalent higher resolution
simulations are plotted in Fig. 6(a). Since this system is constructed
differently from that in §4.3.1, we cannot draw direct comparisons between
these results and the converged reference solution for a smoothed interface.
For example, the instability grows more quickly in this case where the
shearing velocity is discontinuous. However, we do observe qualitatively
similar behaviour when comparing Fig. 6(a) with Fig. 4. The seeded mode grows
more quickly in REMIX simulations than in those using tSPH. The early growth
rate of modes in REMIX simulations is slightly steeper as resolution is
increased, mirroring the behaviour of the analogous simulations in Fig. 4. The
approach of the mode evolution of tSPH simulations towards the REMIX
simulations is also similar here, and again, the lowest resolution REMIX
simulation grows more quickly than the highest resolution tSPH simulation.
Despite this behaviour with increased resolution, high-resolution tSPH
simulations still fail to form spiralling plumes, as surface tension-like
effects continue to dominate, as shown in G.

Figure 7: The evolution of the distribution of the $m_{0}$ kernel geometric
moment in ideal gas KHI simulations with sharp discontinuities. $m_{0}\approx
1$ corresponds to a local distribution of densities that reflects the particle
configuration. Plots show results from simulations with resolution $N_{1}=128$
using the REMIX scheme both (a) without and (b) with the normalising term in
the density evolution.

The effect of the normalising term (§3.7) in a KHI simulation with sharp
discontinuities is demonstrated in Fig. 7. Without it, as the simulation
evolves, $m_{0}$ of some particles drifts away from 1, the value corresponding
to normalisation of the unmodified kernel (see Eqn. 50). The normalising term
ties the density evolution to kernel normalisation, so as the system evolves,
volume elements continue to accurately build up the continuum over which the
kernel function is normalised. In these simulations, the drift in $m_{0}$ does
not noticeably affect the simulation outcome, however, in §4.9, we show an
example where the inclusion of the normalising term is necessary to simulate a
system in hydrostatic equilibrium.

Figure 8: Growth of Kelvin–Helmholtz instabilities between layers with a large
density contrast. We plot snapshots from 3D, ideal gas KHI simulations with
discontinuous initial density and velocity profiles. We show results from
simulations carried out using tSPH, tSPH with artificial conduction, and
REMIX. Snapshots are plotted at two times, from simulations with two
resolutions. The density ratio between the two regions is 1:10.4 and particles
have the same mass across the simulation. Insets show a magnified view of a
KHI plume.

####  4.3.3 KHI with a large density ratio

Capturing the KHI at interfaces in fluids with a large density jump is
additionally challenging for SPH. More smoothing in the density estimate and
larger discretisation error, in equal mass particle simulations, will make
surface tension-like effects stronger at larger density contrasts.
Additionally, using artificial conduction as a method for correcting density
discontinuity is not as effective at larger jumps in density [23]. Our initial
conditions aim to follow those of Price [23] with a density ratio of 1:10,
however, we continue to use 3D simulations to validate our methods for more
typical applications.

Here we construct initial conditions similarly to §4.3.2: sharply
discontinuous in both density and shearing velocity. The low density region is
constructed exactly equivalently with $\rho_{1}=1$, while resolution is
increased in the high density region, following the same method as outlined
previously, such that this region has a density of $\rho_{2}=10.4$. Speeds in
the $x$ direction are again set to $v_{1}=-0.5$ and $v_{2}=0.5$, however the
wavelength of the initial perturbation in the $y$ direction is decreased to
$\lambda=0.128$, although with the same amplitude of $0.01$ [23]. Comparisons
of resolution can not be directly drawn to the previous section, as here fewer
particles will make up individual vortices at a given time due to the
decreased perturbation wavelength.

In Fig. 8, we plot snapshots showing the evolution of these initial conditions
in tSPH, tSPH with conduction, and REMIX simulations for two resolutions. Due
to the lower wavelength of the seeded mode compared with that in previous
sections, we consider simulations with overall higher resolutions, although
this does not necessarily correspond to higher resolution in each individual
vortex, which now occupies a smaller region in the simulation box. The
instability fails to grow with tSPH and grows only slowly in the higher
resolution simulation with conduction. However, the instability is captured
successfully with REMIX, in particular at the higher-resolution, where we
capture spiralling within the plume.

###  4.4 Kelvin–Helmholtz instability – Earth-like iron & rock

Figure 9: Kelvin–Helmholtz instability growth between dissimilar, stiff
materials. We plot snapshots from 3D KHI simulations with multiple, complex
equations of state at densities and pressures representing those at material
boundaries within the Earth. Columns show snapshots at different times with
the top row showing results from simulations using tSPH and the bottom row
using REMIX. The initial density and velocity profiles are discontinuous and
particles have equal mass. These simulations are both relatively low
resolution, with $N_{1}=128$. Individual particles are plotted on a grey
background and coloured by their material type and density. Particles at all
$z$ are plotted, so the grey background is visible in regions that have
maintained their grid alignment in $z$ from the initial conditions.
Figure 10: Kelvin–Helmholtz instability simulations illustrating the interplay
between multiple component methods of the REMIX scheme. Here we show results
from REMIX simulations when using: (a) a more traditional form of the
equations of motion; (b) the traditional, integral-form of the density
estimate, instead of the evolved density; (c) an unmodified Wendland $C^{2}$
kernel, instead of the linear-order reproducing kernels; (d) the full REMIX
SPH scheme. Removing any one of these affects the growth of the instability
significantly. We plot snapshots at $t=2~{}\tau_{\rm KH}$ from 3D simulations
with multiple, stiff equations of state at densities and pressures
representative of those at the core-mantle boundary within the Earth.

Since the evolution of the KHI is predominantly inertial, we expect
instabilities to grow similarly between shearing fluids of different
materials, represented in our simulations as inviscid fluids only differing in
the calculation of pressures and sound speeds through the EoS (§2.2). We
construct similar initial conditions to those used in §4.3.2, but using the
ANEOS Fe85Si15 (iron) and forsterite (rock) EoS with densities and pressures
comparable with those of the Earth’s core-mantle interface [47].

We simulate the KHI at a discontinuity between low-density rock at
$\rho_{1}=5000$ kg m-3 and high-density iron at $\rho_{2}=9550$ kg m-3.
Particles are placed in a periodic box in a configuration exactly matching
that of §4.3.2. These simulations use particles of equal mass. Spatial
dimensions are scaled such that the box spans a length of
$1~{}R_{\oplus}=6371$ km in the $x$ and $y$ dimensions. The velocities in $x$
of the two layers are initialised to $v_{1}=-10^{-4}~{}R_{\oplus}$ s-1,
$v_{2}=10^{-4}~{}R_{\oplus}$ s-1 and the seeded mode has the form
$v_{y}=0.01|v_{1}-v_{2}|\sin{(2\pi x/\lambda)}$ with
$\lambda=0.5~{}R_{\oplus}$. Initial internal energies are calculated through
each material’s EoS such that the regions are in pressure equilibrium with
$P(\rho,u)=1.2\times 10^{11}$ Pa.

In Fig. 9 we show the evolution of a KHI with these initial conditions using
tSPH and REMIX. In the tSPH simulation, surface tension-like effects are
strong. Undesired smoothing of the discontinuity in the SPH density estimate
combined with the stiff equations of state leads to strong artificial forces
at the interface, which both prevent mixing of particles of different
materials and strongly suppress the growth of the instability. These effects
as well as their contributions from zeroth-order error in the momentum
equation are addressed in the construction of the REMIX SPH scheme, so the
instability is allowed to grow and particles of different materials are able
to intermix in a qualitatively similar way to the ideal gas cases.

The mode amplitude growth of these simulations is plotted in Fig. 6(b). We
find strong quantitative similarities between these and the mode growth of the
ideal gas simulations potted in Fig. 6(a). Although we have no experimental or
analytical predictions for the growth of the KHI in these conditions and with
these materials, we find that: (1) spurious surface tension analogous to that
in tSPH KHI simulations with ideal gas is also clearly visible and strong in
tSPH simulations with multiple materials; (2) the construction of the REMIX
scheme is general in, and shown to be effective in, its reduction of
established sources of error in the SPH formalism; (3) without any tuning of
the method to material-specific boundaries, improvements that alleviate
surface tension-like effects in ideal gas KHI simulations also allow the KHI
to form in a qualitatively similar manner in the multi-material case.

To achieve these improved results of the REMIX scheme demonstrated in Fig. 9,
we require interplay between a combination of its constitutive methods (§3).
We use this KHI with iron & rock to highlight the importance of individual
methods included in the REMIX SPH scheme as, while their effects are visible
in all simulations, they present particularly clearly in this case. Fig. 10
shows Earth-like KHI simulations that use the REMIX SPH scheme with different
ones of its constituent methods removed from the construction and reverted to
its traditional SPH analogue in each panel. We show simulations that: (a) use
a more standard form of the equations of motion with equal-valued free
functions (§3.2); (b) use the integral rather than differential form of the
density estimate (§3.1); (c) use an unmodified Wendland $C^{2}$ kernel rather
than linear-order reproducing kernels (§3.3); (d) the full REMIX SPH scheme.
Taking a more traditional approach in any one of these methods leads to much
stronger surface tension-like effects, such that only the full scheme enables
the expected spirals to form. The improvements of the REMIX scheme are in many
cases due to interplay between its constitutive methods all together, rather
than individual components solving separate issues.

###  4.5 Rayleigh–Taylor instability – ideal gas

Figure 11: Rayleigh–Taylor instabilities in an ideal gas. The RTIs are shown
at two times, in simulations using tSPH, tSPH with artificial condition, and
REMIX. These simulations have a resolution of $N_{1}=256$, equal mass
particles, and are performed in 3D. The regions of fixed particles at the top
and bottom of the simulations have been cropped from the figure; their
positions and densities do not change.

We next consider the Rayleigh–Taylor instability (RTI) as an additional
scenario to test the treatment of instability growth and mixing, which, unlike
the previous tests, also includes gravity.

The RTI arises at the interface between a high density fluid being displaced
by a low density fluid [69]. We simulate the gravity-driven growth of the RTI,
where a layer of high density fluid is initially positioned above a layer of
low density fluid (relative to the downward direction of gravity). Hydrostatic
equilibrium is disturbed by a small velocity perturbation. Similarly to the
KHI, surface tension-like effects in traditional SPH formulations strongly
suppress the growth of the RTI.

Our initial conditions are based on those of Frontiere et al. [19]. However,
as with the KHI tests, these simulations are carried out in 3D, with particles
of equal mass, and without deliberate smoothing of the initial density
discontinuity. Particles are placed in a periodic simulation domain in two
cubic lattices, each a square in the $x,y$ dimensions and thin in $z$. The box
has dimensions of 0.5, 1 in the $x$ and $y$ directions, with a thin and
resolution-dependent $z$ box size. The low density region has $N_{1}\times
N_{1}\times 18$ particles with density $\rho_{1}=1$ and occupies the bottom
half of the domain. The high density region is constructed similarly to that
in §4.3.2, giving a density of $\rho_{2}=1.91$ for the upper region while also
ensuring a lattice that is consistent across the periodic simulation box
edges. Particles in the top and bottom 0.05 of the box are fixed in place
throughout the course of the simulation. Initial internal energies are set to
satisfy hydrostatic equilibrium using an ideal gas EoS with $\gamma=7/5$,
constant gravitational acceleration $g=-0.5$, and a pressure at the interface
of $P_{0}=\rho_{2}/\gamma$. Particles are initially at rest, other than an
initial velocity perturbation that seeds the instability,

$v_{y}(x,y)=\begin{cases}\delta_{y}\left[1+\cos\left(8\pi\left(x+0.25\right)\right)\right]\left[1+\cos\left(5\pi\left(y-0.5\right)\right)\right]&\text{for
$0.3<y<0.7$}\;,\\\ 0&\text{otherwise.}\end{cases}$ (55)

We use a perturbation amplitude of $\delta_{y}=0.025$.

In Fig. 11, we show snapshots from RTI simulations with resolution
$N_{1}=256$, simulated using tSPH, tSPH with artificial conduction, and REMIX.
The growth of this instability is strongly suppressed, even with artificial
conduction. REMIX is able to capture the growth of the RTI well. Additionally,
we are able to maintain discontinuities as the simulation evolves. As in the
KHI, these discontinuities are inherently unstable to perturbation modes of
all wavelengths and so we see the growth of secondary, unseeded KHIs and RTIs
that contribute to an onset of turbulent mixing. As the simulation progress,
we observe turbulence driving mixing on the particle scale, the scale of the
primary instability, and in between.

###  4.6 Rayleigh–Taylor instability – Earth-like iron & rock

Figure 12: Rayleigh–Taylor instabilities between dissimilar, stiff materials.
The RTIs are shown at three times, in simulations using tSPH and REMIX. These
simulations have a resolution of $N_{1}=256$, equal mass particles, and are
performed in 3D.

We now consider the treatment of the RTI at an interface between different
materials. The stiff iron & rock EoS that we use makes this an even more
challenging scenario for traditional SPH.

The high density iron layer is placed above the low density rock layer. The
particle configuration is constructed just as in the ideal gas case above,
since the density ratio is taken to be the same. However the box is scaled to
have dimensions $0.05~{}R_{\oplus}$ and $0.1~{}R_{\oplus}$ in the $x$ and $y$
dimensions. The velocity perturbation is similar, although scaled to the box
size and with an amplitude $\delta_{y}=2.5\times 10^{-5}~{}R_{\oplus}$ s-1.
Again, particles are initially in hydrostatic equilibrium, other than due to
the seeded perturbation. Internal energies are chosen to satisfy this for the
constant gravitational acceleration $g=-9.9$ m s-2 and interface pressure
$P_{0}=120$ GPa, representative of the gravitational acceleration and pressure
at the Earth’s core-mantle boundary.

In Fig. 12 we show snapshots from RTI simulations with Earth-like materials
with resolution $N_{1}=256$, with tSPH and REMIX. The RTI does not grow in the
tSPH simulation. In contrast, the behaviour of the REMIX simulation is similar
to the equivalent ideal gas case: unimpeded evolution of the instability,
mixing at different length scales, onset of turbulence, and growth of unseeded
secondary modes.

###  4.7 Blob test

In a physical system, mixing due to fluid instabilities is typically much less
controlled and isolated than in the deliberately seeded scenarios of the
previous sections. The “blob test” [14] is used to investigate the treatment
of turbulent mixing at density discontinuities due to unseeded instabilities.

A spherical cloud of high-density fluid, initially at rest, is placed in an
uniform flow of low-density fluid. Emergent Kelvin–Helmholtz and
Rayleigh–Taylor instabilities at the interface, as well as ram-pressure
stripping, should act to break up the cloud, driving its evolution to a well-
mixed state. As with the instability tests presented in previous sections,
traditional SPH schemes struggle to capture instability growth at density
discontinuities and so the mixing of the cloud into the surrounding fluid is
strongly suppressed. Typically, blob tests are carried out in a supersonic
regime, where interactions between shock waves and the cloud can also be
assessed, applicable to a range of astrophysical scenarios. However, here we
additionally simulate blob tests in a subsonic regime to demonstrate the
ability of the REMIX SPH scheme in capturing subsonic turbulent mixing, which
is even more strongly suppressed in the tSPH formalism [72].

Braspenning et al. [73] compare blob test simulations using seven
hydrodynamical solvers, including SPH schemes and mesh-based methods. We
reproduce their initial conditions to allow direct comparisons with their
simulations. Particles are placed in a 3D periodic box with axes aligned such
that the initially uniform wind flows in the $x$ direction. The length of the
box in $y$ and $z$ is chosen to be 1 pc, and the length in the $x$ direction
is 4 pc. Particles in the ambient wind are initially placed in a cubic lattice
with $4N\times N\times N$ in the $x,y,z$ directions, where we parameterise the
simulation resolution by $N$. We carry out simulations with
$N=16,\,32,\,64,\,128,\,256$. Particles in both the cloud and surrounding wind
have equal masses and so particles in the cloud are placed in a cubic lattice
of higher number density corresponding to the chosen density contrast. We
simulate blob tests with initial density contrasts $\chi=10,\,100$ and the
initial density of the surrounding medium is $10^{-4}$ mp cm-3, where mp is
the proton mass. Clouds are spherical and have a radius of $R_{\rm cloud}=0.1$
pc. Both the cloud and surrounding medium are an ideal gas with $\gamma=5/3$
and internal energies are chosen so that the cloud and surrounding medium are
in pressure equilibrium with each other and the cloud has an initial
temperature of $10^{4}$ K. We carry out simulations with three wind speeds,
characterised by the Mach number $\mathcal{M}\equiv v_{\rm wind}/c_{\rm
wind}$: $\mathcal{M}=1.5$ for a direct comparison to the simulations of
Braspenning et al. [73], $\mathcal{M}=2.7$ the value used most frequently in
validating hydrodynamic methods [14, 19], and $\mathcal{M}=0.5$ to test mixing
in the subsonic regime. We use units of the cloud crushing timescale

$t_{\rm cc}=\frac{\sqrt{\chi}\,R_{\rm cloud}}{v_{\rm wind}}\;,$ (56)

to compare simulations with different initial density contrasts and wind
speeds and for direct comparisons with the results of Braspenning et al. [73].

Figure 13: REMIX simulations of a high-resolution ($N=256$) 3D blob test, at
time $t\approx 5~{}t_{\rm cc}$. The disruption of the cloud is shown in
scenarios with two different initial density contrasts (columns) and three
Mach numbers (rows). Individual particles in a central cross-section in $z$
are plotted, and coloured by their density, $\rho$, relative to the wind
density, $\rho_{\rm wind}$.

In Fig. 13 we plot particle densities from a central cross-section of REMIX
blob test simulations with $N=256$. Results from simulations with three
initial wind speeds and two initial density contrasts are shown for a time
$t\approx 5~{}t_{\rm cc}$. The middle row therefore corresponds directly to
results from simulations plotted in Fig. 1 of Braspenning et al. [73]. REMIX
captures disruption of the cloud in both a subsonic and supersonic regime.
With time, the cloud reaches a well-mixed state with the surrounding medium.
This contrasts with the Braspenning et al. [73] SPH simulations, in which
clouds with $\chi=100$ do not fully mix (their Fig. A1). Additionally, the
onset of turbulence does not produce a highly symmetric structure like that
provoked by the use of a regular grid in simulations using an adaptive mesh
refinement (AMR) method (seen most clearly in Figs. 5 and 6 Braspenning et al.
[73]).

To facilitate direct quantitative comparisons to the simulations of
Braspenning et al. [73], we consider the evolution of the mass of dense gas
and the mass of intermediate-temperature gas. The mass of dense gas, $M_{\rm
dense}$, is defined as the sum of masses of particles with density above a
threshold of $\rho_{\rm cloud}/\,3$, where $\rho_{\rm cloud}$ is the initial
cloud density. The mass of intermediate-temperature gas, $M_{\rm mix}$, is
defined as the sum of particle masses, $m_{i}$, of particles whose
temperature, $T_{i}$, lies within half the logarithmic temperature range
between the cold cloud and the hot wind, centred on the geometric mean
temperature, i.e.

$M_{\rm mix}=\sum_{i}m_{\text{mix},\;i}\quad\text{where}\quad
m_{\text{mix},\;i}=\begin{cases}m_{i}&\text{for }\log(T_{\rm
mix})-\frac{1}{4}\log(\chi)<\log(T_{i})<\log(T_{\rm
mix})+\frac{1}{4}\log(\chi)\;,\\\ 0&\text{otherwise,}\end{cases}$ (57)

where $T_{\rm mix}$ is the geometric mean of the cloud and wind temperatures,
$T_{\rm mix}=\sqrt{T_{\rm cloud}T_{\rm wind}}$. We normalise both $M_{\rm
dense}$ and $M_{\rm mix}$ to the initial cloud mass.

The evolution of $M_{\rm dense}$ and $M_{\rm mix}$ for simulations varying
initial wind speed, initial density contrast, and resolution is plotted in
Fig. 14. We find that REMIX is able to capture the disruption of the cloud in
all these simulations, as shown in Fig. 14(a–f). The middle row corresponds
directly to Figs. 2 and 3 of Braspenning et al. [73]. We see strong
similarities between the behaviour of our REMIX simulations and the
simulations of Braspenning et al. [73] with hydrodynamic solvers that they
find demonstrate good mixing. The evolution of these quantities is well
parameterised by the cloud crushing timescale for this range of Mach numbers,
with features appearing at approximately the same scaled time for all rows.
Increasing resolution results in behaviour that indicates an approach towards
numerical convergence for both $M_{\rm dense}$ and $M_{\rm mix}$, despite the
scenario itself being highly turbulent with no true converged solution.

Figure 14: Time evolution of the mass of dense gas, $M_{\rm dense}$ (a–f), and
the mass of intermediate-temperature gas, $M_{\rm mix}$ (g–l), in blob test
simulations. These quantities are plotted for simulations for two different
initial density contrasts ($\chi$; columns) and three Mach numbers
($\mathcal{M}$; rows). Line colour corresponds to simulation resolution.

###  4.8 Evrard collapse

The Evrard collapse [74] considers the collapse of an isothermal, spherical
cloud of gas under its self-gravity. A shock is formed and moves outwards as
the cloud collapses. We use this test to investigate the coupling of gravity
and hydrodynamics, with large transformations of energy between gravitational,
kinetic and thermal forms.

Figure 15: Evrard collapse with a resolution of ${\sim}10^{7}$ SPH particles
at time $t=0.8$, simulated using REMIX. Plots show (a) radial velocity,
$v_{r}$, (b) density, $\rho$, and (c) specific internal energy, $u$, plotted
against radial distance from the cloud centre, $r$. Individual particles are
plotted in blue, and the dashed red line shows a reference solution from a
high-resolution grid code simulation [60]. Figure 16: Energy evolution in the
Evrard collapse. (a) Different forms of energy, and (b) fractional deviation