# Dusty gas with SPH — II. Implicit timestepping and astrophysical drag
regimes

Guillaume Laibe, Daniel J. Price
Centre for Stellar and Planetary Astrophysics and School of Mathematical
Sciences, Monash University, Clayton, Vic 3800, Australia

###### Abstract

In a companion paper (Laibe & Price, 2011b), we have presented an algorithm
for simulating two-fluid gas and dust mixtures in Smoothed Particle
Hydrodynamics (SPH). In this paper, we develop an implicit timestepping method
that preserves the exact conservation of the both linear and angular momentum
in the underlying SPH algorithm, but unlike previous schemes, allows the
iterations to converge to arbitrary accuracy and is suited to the treatment of
non-linear drag regimes. The algorithm presented in Paper I is also extended
to deal with realistic astrophysical drag regimes, including both linear and
non-linear Epstein and Stokes drag. The scheme is benchmarked against the test
suite presented in Paper I, including i) the analytic solutions of the
dustybox problem and ii) solutions of the dustywave, dustyshock, dustysedov
and dustydisc obtained with explicit timestepping. We find that the implicit
method is 1–10 times faster than the explicit temporal integration when the
ratio $r$ between the the timestep and the drag stopping time is $1\lesssim
r\lesssim 1000$.

###### keywords:

hydrodynamics — methods: numerical — ISM: dust, extinction — protoplanetary
discs — planets and satellites: formation

††pagerange: Dusty gas with SPH — II. Implicit timestepping and astrophysical
drag regimes–References††pubyear: 2011

##  1 Introduction

Dust in cold astrophysical systems spans a huge range of sizes from sub-micron
sized grains in the interstellar medium to kilometre sized planetesimals
involved in planet formation. Moreover, the ratio of dust to gas as well as
the density and temperature of the gaseous environment in which dust is
embedded can also vary strongly. Handling this full range of physical
parameters presents a challenge to numerical schemes designed to simulate
dusty gas in astrophysics. The main challenges are i) that at high drag (e.g.
small grains), the small timestep required means that purely explicit
timestepping methods become prohibitive and ii) that a wide range of physical
drag prescriptions, including non-linear drag regimes, need to be handled by
the code.

Two main prescriptions for drag between gas and solid particles are applicable
to astrophysics: The Epstein regime — where the gas surrounding a grain can be
treated as a dilute medium — and the Stokes regime — where the grains can be
treated as solid bodies surrounded by a fluid —(see e.g. Baines et al. 1965;
Stepinski & Valageas 1996). The dependance of the drag term on local
parameters of the gas (density, temperature) and the dust (typical grain size,
mass) differ between the two regimes, in turn implying very different dynamics
for the dust grains. For example, in a protoplanetary disc, both of these
regimes may be applicable in different regions of the disc.

In a companion paper (Laibe & Price 2011b, hereafter Paper I), we have
developed a new algorithm for treating two-fluid gas-dust astrophysical
mixtures in Smoothed Particle Hydrodynamics (SPH). Benchmarking of the method
demonstrated that the algorithm gives accurate solutions on a range of test
problems relevant to astrophysics and substantially improve previous
algorithms (Monaghan & Kocharyan, 1995). However, in Paper I, we used only a
simple explicit time stepping and considered only linear drag regimes with a
constant drag coefficient. In this paper, we present an implicit timestepping
method that can be applied to both linear and non-linear drag regimes, which
is both more accurate and more general than the scheme proposed by Monaghan
(1997). We also discuss the SPH implementation of both the Epstein and the
Stokes regimes in their full generality.

The paper is organised as follows: The equations of motion and the
characteristics of the different astrophysical drag regimes for gas and dust
mixtures are given in Sec. 2. We summarise the SPH formalism used for
integrating these equations (derived in detail in Paper I) in Sec. 3 and
extend it to deal with the drag regimes encountered in astrophysics.
Particular attention is paid in Sec. 4 to improving the Monaghan (1997)
implicit timestepping scheme, including its generalisation for non-linear drag
regimes. Finally, the algorithm for non-linear drag regimes is tested against
the analytic solutions of the dustybox problem, as well as the dustywave,
dustysedov, dustyshock and dustydisc tests, in Sec. 5.

##  2 Gas and dust evolution in astrophysical systems

###  2.1 Evolution equations

The equations describing the evolution of astrophysical gas and dust mixtures,
where dust is treated as a pressureless, inviscid, continuous fluid have been
described in detail in Paper I. The equations in the continuum limit are given
by:

$\displaystyle\frac{\partial\hat{\rho}_{\mathrm{g}}}{\partial
t}+\nabla.\left(\hat{\rho}_{\mathrm{g}}\textbf{v}_{\mathrm{g}}\right)$
$\displaystyle=$ $\displaystyle 0,$ (1)
$\displaystyle\frac{\partial\hat{\rho}_{\mathrm{d}}}{\partial
t}+\nabla.\left(\hat{\rho}_{\mathrm{d}}\textbf{v}_{\mathrm{d}}\right)$
$\displaystyle=$ $\displaystyle 0,$ (2)
$\displaystyle\hat{\rho}_{\mathrm{g}}\left(\frac{\partial\textbf{v}_{\mathrm{g}}}{\partial
t}+\textbf{v}_{\mathrm{g}}.\nabla\textbf{v}_{\mathrm{g}}\right)$
$\displaystyle=$ $\displaystyle-\theta\phantom{.}\nabla P_{\rm
g}+\hat{\rho}_{\mathrm{g}}\textbf{f}-F_{\rm drag}^{\rm V},$ (3)
$\displaystyle\hat{\rho}_{\mathrm{d}}\left(\frac{\partial\textbf{v}_{\mathrm{d}}}{\partial
t}+\textbf{v}_{\mathrm{d}}.\nabla\textbf{v}_{\mathrm{d}}\right)$
$\displaystyle=$ $\displaystyle-\left(1-\theta\right)\nabla P_{\rm
g}+\hat{\rho}_{\mathrm{d}}\textbf{f}+F_{\rm drag}^{\rm V},$ (4)
$\displaystyle\frac{{\rm d}u_{{\rm g}}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle-\frac{P_{\mathrm{g}}}{\hat{\rho}_{\mathrm{g}}}\left[\theta\nabla\cdot{\bf
v}_{\rm g}+(1-\theta)\nabla\cdot{\bf v}_{\rm d}\right]+\Lambda_{\rm drag}.$
(5)

where the subscripts ${\rm g}$ and ${\rm d}$ refer to the gas and dust,
respectively such that $P_{\rm g}$ is the gas pressure,
$\textbf{v}_{\mathrm{g}}$ and $\textbf{v}_{\mathrm{d}}$ are the fluid
velocities and $u$ is the specific internal energy of gas. The volume
densities of gas and dust ($\hat{\rho}_{\mathrm{g}}$ and
$\hat{\rho}_{\mathrm{d}}$, respectively) are related to the corresponding
intrinsic densities ($\rho_{\rm{g}}$ and $\rho_{\rm{d}}$, respectively)
according to

$\displaystyle\hat{\rho}_{\mathrm{d}}$ $\displaystyle=$
$\displaystyle(1-\theta)\rho_{\mathrm{d}},$ (6)
$\displaystyle\hat{\rho}_{\mathrm{g}}$ $\displaystyle=$
$\displaystyle\theta\rho_{\mathrm{g}},$ (7)

where $\theta$ is the volume filling fraction of the dust. Finally, the drag
force and heating terms are given by:

$F_{\rm drag}^{\rm V}=K(\textbf{v}_{\mathrm{g}}-\textbf{v}_{\mathrm{d}}),$ (8)

and

$\Lambda_{\rm drag}=K(\textbf{v}_{\mathrm{g}}-\textbf{v}_{\mathrm{d}})^{2}.$
(9)

The drag coefficient K has dimensions of mass per unit volume per unit time
and is generally a function of the relative velocity between the two fluids
$\Delta v\equiv|\textbf{v}_{\mathrm{g}}-\textbf{v}_{\mathrm{d}}|$, implying a
non-linear drag regime with respect to the differential velocity between the
gas and the dust. In most of the astrophysical systems, the dust is diluted
enough into the gas so that the gas filling fraction is $\theta\simeq 1$ to a
very good level of approximation, such that the dust buoyancy term
$\left(1-\theta\right)\nabla P_{\rm g}$ is negligible.

###  2.2 Astrophysical drag regimes

Microscopic collisions of gas molecules on a single dust grain result in a net
exchange of momentum which is equivalent to a drag force $\bf{F}_{\rm{drag}}$
between the two phases. Two limiting cases occur when comparing the typical
geometrical size $s$ of a dust grain to the mean free path $\lambda_{\rm{g}}$
of the gas.

When the typical grain size is negligible compared to the collisional mean
free path of the gas particles ($s\ll\lambda_{\rm{g}}$), the grains are
surrounded by a dilute gas phase and may be treated using the Epstein drag
prescription. In this limit, the analytic expression of the resulting drag
force has been derived, assuming spherical, compact grains with homogeneous
composition, for both specular and diffuse reflections on the grain surface
(see Baines et al. 1965 for the complete derivation). These expressions have
been widely used in astrophysical studies (see e.g. Chiang & Youdin 2010 for
references), sometimes incorrectly where the grains are known to be porous and
have fractal structures (Blum & Wurm, 2008) (thus breaking the assumptions of
the Epstein prescription).

For grain sizes larger than the collisional mean free path
($s\gg\lambda_{\rm{g}}$), grains experience a local differential velocity with
respect to a uniform viscous flow and should be treated using the Stokes drag
prescription (Fan & Zhu, 1998). In this case, the momentum is diffused by
viscosity into the fluid, which implies that the drag expression strongly
depends on the local Reynolds number defined according to

$R_{\mathrm{d}}=\frac{2s\left|\textbf{v}_{\mathrm{d}}-\textbf{v}_{\mathrm{g}}\right|}{\nu},$
(10)

where $\nu$ is the kinematic viscosity of the gas. Analytic expressions for
the drag force can be derived at small Reynolds numbers. At higher Reynolds
number, the drag law is inferred from experiments. Rigorously, additional
contributions to the drag should arise from the grain acceleration (carried
mass and Basset contribution), grain rotation (Magnus) and in the presence of
strong local shear, pressure and temperature gradients (e.g. Fan & Zhu 1998).
These corrections are negligible in nearly all astrophysical contexts.

No current analytic theory describes how both the gas and the dust fluid
exchange momentum in the intermediate regime (i.e. $s\simeq\lambda_{\rm{g}}$).
Generally, an asymptotic continuous interpolation between the two limiting
Epstein and Stokes regimes is used. Stepinski & Valageas (1996) suggest
adopting $s=\frac{4}{9}\lambda_{\rm{g}}$ as a means of obtaining a smooth
transition. It should be noted that although this approach is convenient,
there is no clear measure of the physical accuracy of this assumption.

####  2.2.1 Epstein regime for dilute media

In a dilute medium ($\lambda_{\mathrm{g}}>{4s}/{9}$), grains are small enough
not to disturb the Maxwellian distribution of the gas velocity. Assuming
grains are spherical, that the mass of a gas molecule is negligible compared
to the mass of a dust grain and that the reflection of gas particles from
collisions with dust grains are specular, the expression of the drag force on
a single grain $\bf{F}_{\rm{drag}}$ (which differs from the volume force
$F_{\rm drag}^{\rm V}$ by a factor $\hat{\rho}_{\mathrm{d}}/m_{\rm d}$, see
Paper I) for the Epstein regime is given by

$\displaystyle\mathbf{F}_{\rm{drag}}=\displaystyle-2\pi
s^{2}\rho_{\mathrm{g}}\Delta v^{2}$
$\displaystyle\left[\frac{1}{2\sqrt{\pi}}\left\\{\left(\frac{1}{\psi}+\frac{1}{2\psi^{3}}\right)e^{-\psi^{2}}+\right.\right.$
(11)
$\displaystyle\left.\left.\left(1+\frac{1}{\psi^{2}}-\frac{1}{4\psi^{4}}\right)\sqrt{\pi}\,\mathrm{erf}\left(\psi\right)\right\\}\right]\textbf{x},$

where $s$ corresponds to the grain radius and $m$, $\rho_{\mathrm{g}}$, $T$
denote the mass of the gas molecules, the intrinsic gas density and the local
temperature of the mixture (the gas and the dust are supposed to have the same
temperature). The thermal sound speed of the gas is thus
$c_{\mathrm{s}}=\sqrt{\gamma k_{\mathrm{B}}T/m}$ and the mean thermal velocity
of the gas, $c_{\mathrm{s}}\sqrt{8/\pi\gamma}$. The dimensionless quantity
$\psi$ is defined according to

$\psi\equiv\displaystyle\sqrt{\frac{\gamma}{2}}\frac{\Delta
v}{c_{\mathrm{s}}},$ (12)

where $\Delta\textbf{v}=\textbf{v}_{\mathrm{d}}-\textbf{v}_{\mathrm{g}}=\Delta
v$ x is the differential velocity (x being a unit vector). However, depending
on the characteristics of the problem (i.e. low or high Mach numbers, or
both), simpler and computationally less expensive approximations may be used.
For $\psi\ll 1$, i.e. low Mach numbers, Eq. 11 can be expanded to third order
in $\psi$, giving

$\mathbf{F}_{\rm{drag}}=\displaystyle-\frac{4\pi}{3}\rho_{\mathrm{g}}s^{2}\sqrt{\frac{8}{\pi\gamma}}c_{\mathrm{s}}\Delta
v\left[1+\frac{\psi^{2}}{5}+\mathcal{O}\left(\psi^{4}\right)\right]\textbf{x},$
(13)

which is usually simplified to its linear term,

$\mathbf{F}_{\rm{drag}}=-\frac{4\pi}{3}\rho_{\mathrm{g}}s^{2}\sqrt{\frac{8}{\pi\gamma}}c_{\mathrm{s}}\Delta\mathbf{v}.$
(14)

For $\psi\gg 1$, i.e. high Mach numbers, the Taylor expansion in $1/\psi$ of
Eq. 11 gives

$\mathbf{F}_{\rm{drag}}=\displaystyle-\left[\pi\rho_{\mathrm{g}}s^{2}\Delta
v^{2}\left(1+\frac{1}{\psi^{2}}-\frac{1}{4\psi^{4}}\right)+\mathcal{O}\left(e^{-\psi^{2}}\right)\right]\textbf{x},$
(15)

which is usually reduced to its quadratic term,

$\mathbf{F}_{\rm{drag}}=-\pi\rho_{\mathrm{g}}s^{2}\Delta v\Delta\mathbf{v}.$
(16)

A convenient way to handle Epstein drag at both low and high Mach numbers is
to use an interpolation between the two asymptotic regimes given by Eqs. 14
and 16 as derived in Kwok (1975) (cf. Paardekooper & Mellema 2006), giving

$\mathbf{F}_{\rm{drag}}=-\frac{4\pi}{3}\rho_{\mathrm{g}}s^{2}\sqrt{\frac{8}{\pi\gamma}}c_{\mathrm{s}}\sqrt{1+\frac{9\pi}{128}\frac{\Delta
v^{2}}{c_{\rm{s}}^{2}}}\Delta\mathbf{v}.$ (17)

The deviation of Eq. 17 from the full expression (Eq. 11) is $\lesssim 1\%$
(Kwok, 1975). Thus, in general, we adopt Eq. 17 for the Epstein regime. We
compare the differences between the various Epstein expressions in Sec. 5.

####  2.2.2 Stokes regime for dense media

In a dense medium ($\lambda_{\mathrm{g}}>{4s}/{9}$), grains should be treated
with the Stokes drag regime, for which the expression of the drag force
$\bf{F}_{\rm{drag}}$ is:

$\mathbf{F}_{\rm{drag}}=-\frac{1}{2}C_{\mathrm{D}}\pi
s^{2}\rho_{\mathrm{g}}\Delta v\Delta\textbf{v},$ (18)

where the coefficient $C_{\mathrm{D}}$ is a piecewise function of the local
Reynolds number:

$C_{\mathrm{D}}=\begin{cases}24R_{\mathrm{d}}^{-1},&R_{\mathrm{d}}<1;\\\
24R_{\mathrm{d}}^{-0.6},&1<R_{\mathrm{d}}<800;\\\
0.44,&800<R_{\mathrm{d}},\end{cases}$ (19)

where $R_{\mathrm{d}}$ is defined in Eq. 10. Equation 19 indicates that at
small Reynolds numbers ($R_{\mathrm{d}}<1$), the drag force is linear with
respect to the local differential velocity between the grain and the gas. The
relation transitions to a power-law regime
($\mathbf{F}_{\rm{drag}}\propto\Delta v^{0.4}\Delta\textbf{v}$) at
intermediate Reynolds numbers ($1<R_{\mathrm{d}}<800$) and becomes quadratic
at large Reynolds numbers ($R_{\mathrm{d}}>800$). When the local concentration
of dust grains becomes very large (i.e., average distance between the
particles comparable to the grain size), the coefficient $C_{\mathrm{D}}$
should also depend on the local concentration of particles. However, this
extreme situation is not encountered in astrophysical situations.

Assuming gas molecules interact as hard spheres, the dynamic viscosity of the
gas can be computed according to (Chapman & Cowling, 1970):

$\mu=\frac{5m}{64\sigma_{\mathrm{s}}}\sqrt{\frac{\pi}{\gamma}}c_{\mathrm{s}},$
(20)

where $m=2m_{\mathrm{H}}$ and $\sigma_{\mathrm{s}}$ is the geometric cross
section of the molecule ($\sigma_{\mathrm{s}}=2.367\times 10^{-15}$ cm2 for
H2). The gas mean free path $\lambda_{\mathrm{g}}$ and the kinematic viscosity
$\nu$ of the gas are deduced from $\mu$ using

$\lambda_{\mathrm{g}}=\displaystyle\sqrt{\frac{\pi\gamma}{2}}\frac{\mu}{\rho_{\mathrm{g}}c_{\mathrm{s}}},$
(21)

and

$\nu=\displaystyle\frac{\mu}{\rho_{\mathrm{g}}}.$ (22)

##  3 Asytrophysical dust and gas mixtures in SPH

###  3.1 SPH evolution equations

The SPH version of the continuity equations Eqs. 1 – 2 are given by the
density summations for both the gas and the dust phase, computed according to:

$\displaystyle\hat{\rho}_{a}=\sum_{b}m_{b}W_{ab}(h_{a});$ $\displaystyle
h_{a}=\eta\left(\frac{m_{a}}{\hat{\rho}_{a}}\right)^{1/\nu},$ (23)
$\displaystyle\hat{\rho}_{i}=\sum_{j}m_{j}W_{ij}(h_{i});$ $\displaystyle
h_{i}=\eta\left(\frac{m_{j}}{\hat{\rho}_{i}}\right)^{1/\nu},$ (24)

where as in Paper I, the indices $a,b,c$ refer to quantities computed on gas
particles and $i,j,k$ refer to quantities computed on dust particles. The
volume filling fraction $\theta$, is defined on a _gas_ particle, $a$,
according to

$\theta_{a}=1-\frac{\hat{\rho}_{{\rm d},a}}{\rho_{\rm d}},$ (25)

where $\hat{\rho}_{{\rm d},a}$ is the density of _dust_ at the _gas_ particle
location, calculated using

$\hat{\rho}_{{\rm d},a}=\sum^{N_{neigh,dust}}_{j=1}m_{j}W_{aj}(h_{a}),$ (26)

where $h_{a}$ is the smoothing length of the _gas_ particle computed using gas
neighbours. The local density of dust at the gas location can thus be zero
(giving $\theta=1$) if no dust particles are found within the kernel radius
computed with the gas smoothing length. Importantly, as $\hat{\rho}$ and $h$
are mutually dependent, they have to be simultaneously calculated for each
type of particle, e.g. by the iterative procedure described in Price &
Monaghan (2007).

The SPH equations of motion for the gas and the dust particles, corresponding
to the SPH translation of Eqs. 3 and 4, are given by

$\displaystyle\frac{\mathrm{d}\textbf{v}_{a}}{\mathrm{d}t}=$
$\displaystyle-\sum_{b}m_{b}\left[\frac{P_{a}\tilde{\theta}_{a}}{\Omega_{a}\hat{\rho}_{a}^{2}}\nabla_{\\!\\!\
a}W_{\\!\\!\ a\\!\\!\
b}\left(h_{a}\right)+\frac{P_{b}\tilde{\theta}_{b}}{\Omega_{b}\hat{\rho}_{b}^{2}}\nabla_{\\!\\!\
a}W_{\\!\\!\ a\\!\\!\ b}\left(h_{b}\right)\right]$
$\displaystyle-\sum_{j}m_{j}\frac{P_{a}\left(1-\theta_{a}\right)}{\hat{\rho}_{a}\hat{\rho}_{{\rm
d},a}}\nabla_{\\!\\!\ a}W_{\\!\\!\ a\\!\\!\ j}\left(h_{a}\right)$
$\displaystyle+\nu\sum_{j}m_{j}\frac{K_{aj}}{\hat{\rho}_{a}\hat{\rho}_{j}}\left({\bf
v}_{aj}\cdot\hat{\textbf{r}}_{aj}\right)\hat{\textbf{r}}_{aj}D_{aj}(h_{a}),$
(27)

for an SPH gas particle and

$\displaystyle\frac{\mathrm{d}{\bf v}_{i}}{\mathrm{d}t}=$
$\displaystyle\sum_{b}m_{b}\frac{P_{b}\left(1-\theta_{b}\right)}{\hat{\rho}_{b}\hat{\rho}_{{\mathrm{d},b}}}\nabla_{\\!\\!\
i}W_{\\!\\!\ b\\!\\!\ i}\left(h_{b}\right)$ (28)
$\displaystyle-\nu\sum_{b}m_{b}\frac{K_{bi}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left({\bf
v}_{bi}\cdot\hat{\textbf{r}}_{bi}\right)\hat{\textbf{r}}_{bi}D_{ib}(h_{i}),$

for an SPH dust particle. $\Omega$ is the usual variable smoothing length term

$\Omega_{b}\equiv 1-\frac{\partial
h_{b}}{\partial\hat{\rho}_{b}}\sum_{c}m_{c}\frac{\partial W_{\\!\\!\ b\\!\\!\
c}\left(h_{b}\right)}{\partial h_{b}}.$ (29)

It should be noted that $\Omega_{\rm{d}}$ is computed only using dust particle
neighbours according to:

$\Omega_{\mathrm{d},b}=1-\frac{\partial
h_{b}}{\partial\hat{\rho}_{{\mathrm{d},b}}}\sum_{j}m_{j}\frac{\partial
W_{\\!\\!\ b\\!\\!\ j}\left(h_{b}\right)}{\partial h_{b}}.$ (30)

$\tilde{\theta}$ is defined according to

$\tilde{\theta}\equiv\theta+\frac{\hat{\rho}_{\rm
g}}{\hat{\rho}_{\mathrm{d}}}(1-\theta)(1-\Omega_{\rm d}).$ (31)

At this stage, no assumptions are made with respect to the functional form of
the drag coefficient $K$. The evolution of the internal energy for an SPH gas
particle is given by

$\displaystyle\frac{{\rm d}u_{a}}{{\rm d}t}$
$\displaystyle=\frac{\tilde{\theta}_{a}P_{a}}{\Omega_{a}\hat{\rho}_{a}^{2}}\sum_{b}m_{b}\left(\textbf{v}_{a}-\textbf{v}_{b}\right)\cdot\nabla_{a}W_{ab}(h_{a})$
(32) $\displaystyle+\frac{(1-\theta_{a})P_{a}}{\hat{\rho}_{a}\hat{\rho}_{{\rm
d},a}}\sum^{N_{neigh,dust}}_{j=1}m_{j}\left({\bf v}_{a}-{\bf
v}_{j}\right)\cdot\nabla_{a}W_{aj}(h_{a})$
$\displaystyle+\nu\sum_{k}m_{k}\frac{K_{ak}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\textbf{v}_{ak}\cdot\hat{\textbf{r}}_{ak}\right)^{2}D_{ak}(h_{a}).$

In Paper I, we showed that the total linear and angular momentum as well as
the total energy are exactly conserved. Thermal coupling terms have been
neglected in this paper.

###  3.2 Kernel functions

Two different kernels are employed to perform the SPH interpolations. First, a
standard bell-shaped kernel $W$:

$W\left(r,h\right)=\frac{\sigma}{h^{\nu}}f\left(q\right),$ (33)

where $h$ denotes the smoothing lengths of each phases, $\nu$ the number of
spatial dimensions and $q\equiv|{\bf r}-{\bf r}^{\prime}|/h$ is the
dimensionless variable used to calculate the densities and the buoyancy terms.
The function $f$ is usually the $M_{4}$ cubic spline kernel (Monaghan, 2005).
The drag interpolation is performed using a second kernel $D$. As shown in
Paper I, double-hump shaped kernels given by

$D\left(r,h\right)=\frac{\tilde{\sigma}}{h^{\nu}}q^{2}f(q),$ (34)

significantly improve the accuracy of the drag interpolation — for the same
computational cost — compared to bell-shaped kernels. The normalisation
constants $\tilde{\sigma}$ for various double hump kernels are given in Paper
I. We adopt the double hump cubic for the drag terms in this paper.

###  3.3 Astrophysical drag regimes in SPH

####  3.3.1 Gas viscosity and mean free path

The drag coefficients $K_{ak}$ involved in Eqs. 3.1 – 28 and 32 are computed
independently for each pair of any gas particle $a$ and dust particle $k$. We
first use the sound speed $c_{\rm{s},a}$ to estimate the viscosity $\mu_{a}$
on the gas particle $a$ using (see Eq. 20)

$\mu_{a}=\frac{5m}{64\sigma_{\mathrm{s}}}\sqrt{\frac{\pi}{\gamma}}c_{\rm{s},a}.$
(35)

The mean free path is then computed according to Eq. 21, giving

$\lambda_{\mathrm{g},a}=\displaystyle\sqrt{\frac{\pi\gamma}{2}}\frac{\mu_{\rm{a}}}{\hat{\rho}_{a}c_{\rm{s},a}}.$
(36)

Finally, $\lambda_{\mathrm{g},a}$ is compared to the quantity $4s_{k}/9$ —
$s_{k}$ being the grain size of the dust particle — to determine whether the
drag coefficient of the SPH pair $K_{ak}$ is calculated using the Epstein or
the Stokes drag regimes.

####  3.3.2 Epstein regime

If $4s_{k}/9\leq\lambda_{\mathrm{g},a}$, the drag coefficient $K_{ak}$ is
calculated using the Epstein prescription. Introducing the SPH quantity
$\psi_{ak}$ calculated on a pair of gas and dust SPH particles and defined by

$\psi_{ak}\equiv\sqrt{\frac{\gamma}{2}}\frac{|{\bf
v}_{ak}|}{c_{\mathrm{s},a}},$ (37)

Eq. 11 can be straightforwardly translated to get the drag coefficient
$K_{ak}$ involved in the SPH drag force

$\displaystyle
K_{ak}=\displaystyle-\sqrt{\pi}s^{2}\rho_{\mathrm{g}}\frac{\hat{\rho}_{\mathrm{d}}}{m_{\rm{d}}}|\textbf{v}_{ak}|$
$\displaystyle\left[\frac{1}{2\sqrt{\pi}}\left\\{\left(\frac{1}{\psi_{ak}}+\frac{1}{2\psi_{ak}^{3}}\right)e^{-\psi_{ak}^{2}}+\right.\right.$
(38)
$\displaystyle\left.\left.\left(1+\frac{1}{\psi_{ak}^{2}}-\frac{1}{4\psi_{ak}^{4}}\right)\sqrt{\pi}\,\mathrm{erf}\left(\psi_{ak}\right)\right\\}\right]\textbf{x},$

where $s$ is the grain radius, $m_{\rm{d}}$ is the grain mass and $\gamma$ is
the adiabatic index. Eq. 38 is computationally expensive as it involves
exponential and error functions. The SPH equivalent of Eq.17 is given by

$K_{ak}=\frac{4}{3}\pi\sqrt{\frac{8}{\pi\gamma}}\frac{\hat{\rho}_{k}}{m_{\rm{d}}}\frac{\hat{\rho}_{a}}{\theta_{a}}s^{2}c_{\mathrm{s},a}\sqrt{1+\frac{9\pi}{128}\frac{v_{ak}^{2}}{c_{\mathrm{s},a}^{2}}}.$
(39)

Both Eqs. 38 and 39 reduce to the linear Epstein regime at low Mach numbers
(equivalent of Eq. 14) for which the coefficient $K_{ak}$ is

$K_{ak}=\frac{4}{3}\pi\sqrt{\frac{8}{\pi\gamma}}\frac{\hat{\rho}_{k}}{m_{\rm{d}}}\frac{\hat{\rho}_{a}}{\theta_{a}}s^{2}c_{\mathrm{s},a},$
(40)

and to the quadratic drag regime at high Mach numbers (equivalent of Eq. 16),
for which the coefficient $K_{ak}$ is

$K_{ak}=\pi\rho_{\mathrm{g}}s^{2}\frac{\hat{\rho}_{\mathrm{d}}}{m_{\rm{d}}}|{\bf
v}_{ak}|.$ (41)

####  3.3.3 Stokes regime

If $4s_{k}/9>\lambda_{\mathrm{g},a}$, the drag coefficient $K_{ak}$ is
calculated using the Stokes prescription (see Eqs. 18–19). The local Reynolds
number $R_{\mathrm{d},ak}$ is computed for each pair of gas and dust particles
using

$R_{\mathrm{d},ak}\equiv\frac{2s\hat{\rho}_{a}\left|\textbf{v}_{ak}\right|}{\mu_{a}\theta_{a}},$
(42)

such that the drag coefficient $K_{ak}$ can be computed according to

$K_{ak}=\begin{cases}\displaystyle
6\pi\frac{\hat{\rho}_{k}}{m_{\rm{d}}}\mu_{a}s&R_{\mathrm{d},ak}<1,\\\
\displaystyle\frac{12\pi}{2^{0.6}}\frac{\hat{\rho}_{k}}{m_{\rm{d}}}\frac{\mu_{a}^{0.6}}{\theta_{a}^{0.4}\hat{\rho}_{a}^{0.6}}s^{1.4}\left|\textbf{v}_{ak}\right|^{0.4}&1<R_{\mathrm{d},ak}<800,\\\
\displaystyle
0.22\pi\frac{\hat{\rho}_{k}}{m_{\rm{d}}}\frac{\hat{\rho}_{a}}{\theta_{a}}s^{2}\left|\textbf{v}_{ak}\right|&R_{\mathrm{d},ak}>800.\\\
\end{cases}$ (43)

These expressions have been used by Ayliffe et al. (2011) to compute the drag
on planetesimals in a protoplanetary disc.

##  4 Timesteping

###  4.1 Explicit timesteping

The simplest method to evolve the evolution equations for the SPH particles is
to use an explicit integrator (e.g. the standard Leapfrog). The stability of
the system is guaranteed provided the timestep remains smaller than a critical
value $\Delta t_{\rm{c}}$. In Paper I, we performed a Von Neumann analysis of
the continuous equations, deriving the explicit timestepping criterion

$\Delta
t_{\rm{c},a}=\min_{k}\left[\frac{\hat{\rho}_{a}\hat{\rho}_{k}}{K_{ak}(\hat{\rho}_{a}+\hat{\rho}_{k})}\right];\hskip
14.22636pt\Delta
t_{\rm{c},i}=\min_{b}\left[\frac{\hat{\rho}_{b}\hat{\rho}_{i}}{K_{bi}(\hat{\rho}_{b}+\hat{\rho}_{i})}\right];$
(44)

for gas and dust particles, respectively, with the minimum being taken over
all the particle’s neighbours. Although this criterion was derived in Paper I
for linear drag regimes only, it remains valid even for non-linear drag
regimes where the drag coefficients depend on the differential velocity
between the particles, i.e. $K_{ak}=K_{ak}\left(|\mathbf{v}_{ak}|\right)$.

###  4.2 Implicit timestepping

When the drag timescale becomes smaller than other time scales in the system
(e.g. the Courant condition or the orbital timescale), the timestep
restriction of the explicit methods may become prohibitive and implicit
methods are required. Monaghan (1997) considered the application of two
implicit schemes (the first-order Backward-Euler and second-order Tischer
scheme) to SPH dust-gas mixtures. Both schemes are unconditionally stable, but
a higher accuracy is achieved with second-order schemes.

####  4.2.1 Backward-Euler method

The Backward-Euler scheme applied to the drag interaction between SPH dust and
gas particles is given by

$\displaystyle\frac{\mathbf{v}^{n+1}_{a}-\mathbf{v}^{n}_{a}}{\Delta t}$
$\displaystyle=$
$\displaystyle-\nu\sum_{k}m_{k}\frac{K_{ak}^{n+1}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\mathbf{v}^{n+1}_{ak}\cdot\hat{\mathbf{r}}_{ak}\right)\hat{\mathbf{r}}_{ak}D_{ak},$
(45) $\displaystyle\frac{\mathbf{v}^{n+1}_{i}-\mathbf{v}^{n}_{i}}{\Delta t}$
$\displaystyle=$
$\displaystyle+\nu\sum_{b}m_{b}\frac{K_{bi}^{n+1}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left(\mathbf{v}^{n+1}_{bi}\cdot\hat{\mathbf{r}}_{bi}\right)\hat{\mathbf{r}}_{bi}D_{bi}.$
(46)

Although the scheme is unconditionally stable, the implicit equation with
respect to the velocities $\textbf{v}^{n+1}$ must be solved at each time step.
Direct numerical inversion of this linear system would be prohibitive given
the typical number of neighbour interactions for each SPH particle. Thus,
approximate or iterative solutions to Eqs. 45 – 46 are required.

####  4.2.2 Monaghan (1997) scheme

Monaghan (1997) suggested approximating the velocities $\textbf{v}^{n+1}$ of
Eqs. 45 – 46 using a pairwise treatment in order to preserve the exact
conservation of linear and angular momentum in the SPH formalism. Considering
the interaction between the SPH gas particle $a$ the dust particle $i$,
Monaghan (1997) introduced pairwise auxiliary velocities $\tilde{\bf{v}}$
defined by:

$\displaystyle\tilde{\textbf{v}}_{a}$ $\displaystyle=$ $\displaystyle{\bf
v}_{a}^{n}-m_{i}~{}\Delta t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left(\tilde{\textbf{v}}_{ai}\cdot\hat{\bf{r}}_{ai}\right)\hat{\bf{r}}_{ai},$
(47) $\displaystyle\tilde{\textbf{v}}_{i}$ $\displaystyle=$ $\displaystyle{\bf
v}_{i}^{n}+m_{a}~{}\Delta t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left(\tilde{\textbf{v}}_{ai}\cdot\hat{\bf{r}}_{ai}\right)\hat{\bf{r}}_{ai},$
(48)

Eqs. 47 and 48 are solved, for a given pair of particles, by taking the scalar
product by $\hat{\bf{r}}_{ai}$ of the difference of the two equations, giving

$\tilde{\textbf{v}}_{ai}\cdot\hat{\bf{r}}_{ai}=\frac{\textbf{v}_{ai}^{n}\cdot\hat{\bf{r}}_{ai}}{1+\Delta
t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left(m_{a}+m_{i}\right)}.$ (49)

Substituting this expression into Eq. 47 and 48 gives expressions for
$\tilde{\textbf{v}}_{a}$ and $\tilde{\textbf{v}}_{i}$. Iterating this pairwise
process by looping over all the SPH particles provides an approximate solution
for the velocities $\bf{v}^{n+1}$, i.e.

$\displaystyle\frac{\mathbf{v}^{n+1}_{a}-\mathbf{v}^{n}_{a}}{\Delta t}$
$\displaystyle\simeq$
$\displaystyle-\nu\sum_{k}m_{k}\frac{K_{ak}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\tilde{\mathbf{v}}_{ak}\cdot\hat{\mathbf{r}}_{ak}\right)\hat{\mathbf{r}}_{ak}D_{ak},$
(50) $\displaystyle\frac{\mathbf{v}^{n+1}_{i}-\mathbf{v}^{n}_{i}}{\Delta t}$
$\displaystyle\simeq$
$\displaystyle+\nu\sum_{b}m_{b}\frac{K_{bi}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left(\tilde{\mathbf{v}}_{bi}\cdot\hat{\mathbf{r}}_{bi}\right)\hat{\mathbf{r}}_{bi}D_{bi}.$
(51)

The main drawback of this method is that the approximation given by Eqs. 50
and 51 is inexact – that is, it provides only an approximate solution to Eqs.
45 and 46. Furthermore the accuracy of the approximation is not known _a
priori_ and there is no possibility of performing repeated sweeps in order to
converge to a more accurate solution. In practice, we find that the velocities
obtained by this scheme (for example on the dustybox test) can be
significantly in error, with no possibility of improving the convergence (for
example, by doing several iterations/sweeps).

####  4.2.3 Alternative pairwise treatment for linear drag regimes

We propose a more consistent method for solving Eqs. 45–46 on a given gas or
dust particle ($a$ and $i$, respectively) by sweeping over all particle pairs
and updating the velocities iteratively according to

$\displaystyle{\textbf{v}}_{a}^{**}$ $\displaystyle=$ $\displaystyle{\bf
v}_{a}^{n}+\Delta t\textbf{F}_{a,\mathrm{drag}}^{*}$ (52) $\displaystyle-$
$\displaystyle m_{i}~{}\Delta t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left[\left({\textbf{v}}_{ai}^{**}-\textbf{v}_{ai}^{*}\right)\cdot\hat{\bf{r}}_{ai}\right]\hat{\bf{r}}_{ai},$
$\displaystyle{\textbf{v}}_{i}^{**}$ $\displaystyle=$ $\displaystyle{\bf
v}_{i}^{n}+\Delta t\textbf{F}_{i,\mathrm{drag}}^{*}$ (53) $\displaystyle+$
$\displaystyle m_{a}~{}\Delta t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left[\left({\textbf{v}}_{ai}^{**}-\textbf{v}_{ai}^{*}\right)\cdot\hat{\bf{r}}_{ai}\right]\hat{\bf{r}}_{ai},$

where ${\bf v}^{**}$ refers to the improved approximation to ${\bf v}^{n+1}$
obtained after updating each pair and ${\bf v}^{*}$ to the previous iteration
value of ${\bf v}^{**}$. Eqs. 52 and 53 are solved for each pair of particles
by taking the dot product of $\hat{\bf{r}}_{ai}$ with the difference of the
two equations, giving

$\left({\textbf{v}}_{ai}\cdot\hat{\bf{r}}_{ai}\right)^{**}=\frac{\left(\textbf{v}_{ai}^{n}+\Delta
t\textbf{F}_{ai,\mathrm{drag}}^{*}+\Delta t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left(m_{a}+m_{i}\right)\textbf{v}_{ai}^{*}\right)\cdot\hat{\bf{r}}_{ai}}{1+\Delta
t\frac{\nu
K_{ai}D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}\left(m_{a}+m_{i}\right)}.$ (54)

Substituting Eq. 54 in Eqs. 52 and 53 gives the updated velocities for the
pair. Note that during the global sweep over particle pairs ${\bf v}^{*}$
begins as ${\bf v}^{n}$ at the first iteration but is updated as soon as new
values become available.

This pairwise correction ensures that i) both the linear and the angular
momentum are exactly conserved and ii) the velocities converge to the correct
solution of the implicit scheme given by Eqs. 45 and 46, since the last term
of Eqs. 52 and 53 tends to zero as the number of iterations increases. We thus
refine our approximation to the solution by performing as many successive
iterations as are required to reach a suitable convergence criterion.

####  4.2.4 Convergence criterion

We consider that the approximation we obtain from the implicit scheme
described above is accurate enough when

$\frac{|\bf{v}^{k+1}-\bf{v}^{k}|}{\min c_{\rm{s}}}<\varepsilon,$ (55)

is satisfied for each particle. Typically, we adopt $\varepsilon=10^{-4}$,
which ensures that the approximation we make on the time stepping is
negligible compared to the $\mathcal{O}(h^{2})$ truncation error of the
underlying SPH scheme.

####  4.2.5 Implementation into Leapfrog

The Leapfrog scheme is well suited to the evolution of particle methods
because, for position-dependant forces, it preserves geometric properties of
particle orbits and requires only one evaluation per timestep to give second
order accuracy. In the standard formulation, the evolution is computed
according to

$\begin{array}[]{lrrcll}{\rm Kick}&\left[\right.&{\bf v}^{1/2}&=&{\bf
v}^{0}+\frac{\Delta t}{2}{\bf f}^{0}\left({\bf x}^{0},{\bf
v}^{0}\right),&\left.\right]\\\\[5.0pt] {\rm Drift}&\left[\right.&{\bf
x}^{1}&=&{\bf x}^{0}+\Delta t{\bf v}^{1/2},&\left.\right]\\\\[5.0pt] {\rm
Kick}&\left[\right.&{\bf v}^{1}&=&{\bf v}^{1/2}+\frac{\Delta t}{2}{\bf
f}^{1}\left({\bf x}^{1},{\bf v}^{1}\right),&\left.\right]\end{array}$ (56)

corresponding to Kick, Drift and Kick steps respectively. Adapting Leapfrog to
deal with velocity dependent forces (e.g. drag) is a priori more difficult
since for velocity-dependent forces, the last Kick is implicit in ${\bf
v}^{1}$. For our present purposes, this does not present a major problem since
the drag is already computed implicitly. Splitting the forces into position-
dependent (${\bf f}_{\rm SPH}$) and drag (${\bf f}_{\rm drag}$) contributions,
the scheme becomes

$\begin{array}[]{lrrcll}{\rm Kick}&\left[\right.&{\bf\tilde{v}}^{1/2}&=&{\bf
v}^{0}+\frac{\Delta t}{2}{\bf f}^{0}_{\rm SPH}\left({\bf
x}^{0}\right),&\left.\right]\\\\[5.0pt] {\rm Drift}&\left[\right.&{\bf
x}^{1/2}&=&{\bf x}^{0}+\frac{\Delta t}{2}\tilde{{\bf
v}}^{1/2},&\left.\right]\\\\[5.0pt] {\rm Drag}&\left[\right.&{\bf
v}^{1/2}&=&{\bf\tilde{v}}^{1/2}+\frac{\Delta t}{2}{\bf f}^{1/2}_{\rm
drag}\left({\bf x}^{1/2},{\bf v}^{1/2}\right)&\left.\right]\\\\[5.0pt] {\rm
Drift}&\left[\right.&{\bf x}^{1}&=&{\bf x}^{0}+\Delta t{\bf
v}^{1/2},&\left.\right]\\\\[5.0pt] {\rm
Kick}&\left[\right.&{\bf\tilde{v}}^{1}&=&{\bf v}^{1/2}+\frac{\Delta t}{2}{\bf
f}^{1}_{\rm SPH}\left({\bf x}^{1}\right),&\left.\right]\\\\[5.0pt] {\rm
Drag}&\left[\right.&{\bf v}^{1}&=&{\bf\tilde{v}}^{1}+\frac{\Delta t}{2}{\bf
f}^{1}_{\rm drag}\left({\bf x}^{1},{\bf
v}^{1}\right)&\left.\right]\end{array}$ (57)

where the Drag steps represent the implicit updates computed as described in
Sec. 4.2.3. The disadvantage of Eq. 57 is that two drag force evaluations are
required, removing one of the advantages of the Leapfrog integrator.
Inspection of 57 reveals that an alternative version that requires only one
Drag step can be constructed according to

$\begin{array}[]{lrlr}{\rm Kick}&\bigl{[}&\begin{array}[]{rcl}{\bf
v}^{1/2}&=&{\bf v}^{0}+\frac{\Delta
t_{0}}{2}{\bf\tilde{f}},\end{array}&\bigr{]}\\\\[15.00002pt] {\rm
Drift}&\bigl{[}&\begin{array}[]{rcl}{\bf x}^{1}&=&{\bf x}^{0}+\Delta t_{0}{\bf
v}^{1/2},\end{array}&\bigr{]}\\\\[10.00002pt] {\rm
Drag}&&\left\\{\begin{array}[]{rcl}{\tilde{\bf v}}^{3/2}&=&{\bf
v}^{1/2}+\frac{\Delta t_{0}+\Delta t_{1}}{2}{\bf f}^{1}_{\rm SPH}\left({\bf
x}^{1}\right),\\\\[5.0pt] {\bf v}^{3/2}&=&{\tilde{\bf v}}^{3/2}+\frac{\Delta
t_{0}+\Delta t_{1}}{2}{\bf f}^{1}_{\rm drag}\left({\bf x}^{1},{\bf
v}^{3/2}\right),\\\\[5.0pt] {\tilde{\bf f}}&=&2\left({\tilde{\bf
v}}^{3/2}-{\tilde{\bf v}}^{1/2}\right)/\left(\Delta t_{0}+\Delta
t_{1}\right),\end{array}\right.&\\\\[30.00005pt] {\rm
Kick}&\bigl{[}&\begin{array}[]{rcl}{\bf v}^{1}&=&{\bf v}^{1/2}+\frac{\Delta
t_{0}}{2}{\bf\tilde{f}}.\end{array}&\bigr{]}\end{array}$ (58)

where we have combined the drag steps by predicting the velocity ${\bf
v}^{3/2}$. Note that strictly, the Drag step in this method is semi-implicit
since the force is evaluated using ${\bf x}^{1}$ rather than ${\bf x}^{3/2}$.
However, we expect this approximation to be reasonable as at high drag (for
which the implicit method is designed), the drag mainly changes the
differential velocity between the fluids and has less of an effect on the
positions. Care is also required when the timestep changes between the steps.
We have indicated the correct procedure by specifying $\Delta t_{0}$ and
$\Delta t_{1}$ where $\Delta t_{1}$ is the timestep computed based on ${\bf
x}^{1}$. Finally, Eq. 58 requires that ${\tilde{\bf f}}$ is known at the
beginning of the integration. This can be easily achieved by performing the
Drag step in Eq. 58 with ${\bf v}^{1/2}={\bf v}^{0}$, $\Delta t_{0}=0$ and
$\Delta t_{1}$ equal to the timestep calculated using the initial particle
positions.

####  4.2.6 Generalisation to non-linear drag regimes

To extend this alternative pairwise treatment to any non-linear drag regime,
two additional points have to be considered. Firstly, although in principle
six quantities ($v_{x,y,z}$ for each particle) have to be determined for each
pair, this can be reduced to a single unknown quantity since the drag
coefficient depends only on the modulus of the differential velocity and the
exchange of momentum is directed along the line of sight joining the
particles. The system of equations for a single pair thus reduces to

$\displaystyle{\textbf{v}}_{a}^{**}$ $\displaystyle=$ $\displaystyle{\bf
v}_{a}^{n}+\Delta t\textbf{F}_{a,\mathrm{drag}}^{*}$ (59) $\displaystyle-$
$\displaystyle m_{i}~{}\Delta t\frac{\nu
D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}K_{ai}\sqrt{\left[\left({\textbf{v}}_{ai}^{**}-\textbf{v}_{ai}^{*}\right)\cdot\hat{\bf{r}}_{ai}\right]^{2}+V^{2,n}_{\rm
orth}},$ $\displaystyle{\textbf{v}}_{i}^{**}$ $\displaystyle=$
$\displaystyle{\bf v}_{i}^{n}+\Delta t\textbf{F}_{i,\mathrm{drag}}^{*}$ (60)
$\displaystyle+$ $\displaystyle m_{a}~{}\Delta t\frac{\nu
D_{ai}}{\hat{\rho}_{a}\hat{\rho}_{i}}K_{ai}\sqrt{\left[\left({\textbf{v}}_{ai}^{**}-\textbf{v}_{ai}^{*}\right)\cdot\hat{\bf{r}}_{ai}\right]^{2}+V^{2,n}_{\rm
orth}},$

where

$V^{2,n}_{\rm
orth}={\textbf{v}}_{ai}^{n}\cdot{\textbf{v}}_{ai}^{n}-\left({\textbf{v}}_{ai}^{n}\cdot\hat{\bf{r}}_{ai}\right)^{2}.$
(61)

Secondly, taking the dot product of $\hat{\bf{r}}_{ai}$ with the difference of
the two equations Eqns. 59 and 60 does not lead in general to an equation
which can be solved analytically. The values of
${\textbf{v}}_{ai}\cdot\hat{\bf{r}}_{ai}$ must therefore be determined using a
numerical rootfinding procedure (we use a Newton-Raphson scheme) before being
substituted in Eqns. 59 and 60 to determine the velocities for both the gas
and the dust particles.

####  4.2.7 Performance of the implicit scheme

The computational cost of a timestep with the implicit pairwise treatment is
more expensive than an explicit timestep since at least two iterations have to
be performed to ensure that the scheme is converged. However, the implicit
pairwise treatment will be more efficient provided that the number of
iterations is much smaller than the number of explicit timesteps that would
otherwise be required.

We find in practice that the efficiency of the pairwise treatment is mainly
determined by the number of iterations required to satisfy Eq. 55 (this aspect
was not addressed in the Monaghan (1997) scheme where only one iteration is
ever taken in the hope that the approximation is sufficiently accurate). The
rapidity of the convergence depends primarily on the ratio $r=\Delta t/t_{\rm
s}$ of the timestep over the drag stopping time (defined in Eq. (96) of Paper
I) and on the value of $\varepsilon$. Empirically, we have found that, for
$\varepsilon=10^{-4}$ and $1\lesssim r\lesssim 10$, the implicit pairwise
treatment converges efficiently, the ratio $|\bf{v}^{k+1}-\bf{v}^{k}|/(\min
c_{\rm{s}})$ decreasing by $\sim$ two orders of magnitude at each iterations.
Thus, the implicit pairwise treatment improves the computational time by a
factor of $\sim 1$–$10$. However, this rapidity of convergence decreases as
$r$ increases. At very high drag ($r\gtrsim 1000$), we find that the implicit
scheme becomes less efficient than explicit timestepping due to the large
number of iterations required. A similar behaviour has been found using the
Gauss-Seidel iterative scheme developed by Whitehouse et al. (2005) to treat
SPH radiative transfer in the flux-limited diffusion approximation (Bate 2011,
private communication), so this issue is not specific to the pairwise
treatment.

It is important to note that the computational gain obtained with the pairwise
scheme does not solve the resolution issue at high drag extensively discussed
in Paper I. Both of these problems suggest that a more efficient method for
handling high drag regimes is required. Such a method is beyond the scope of
the present paper.

####  4.2.8 Higher order implicit schemes

Higher temporal accuracy may be achieved by using second instead of first
order implicit schemes. The gain in accuracy is obtained by dividing the drag
timestep $\Delta t$ into two half timesteps. Monaghan (1997) suggested the
‘Tischer’ scheme, where the two half timesteps are given by

$\displaystyle\frac{\mathbf{v}^{n+\frac{1}{2}}_{a}-\mathbf{v}^{n}_{a}}{\Delta
t/2}$ $\displaystyle=$
$\displaystyle-0.6\nu\sum_{k}m_{k}\frac{K_{ak}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\mathbf{v}^{n+\frac{1}{2}}_{ak}\cdot\hat{\mathbf{r}}_{ak}\right)\hat{\mathbf{r}}_{ak}D_{ak},$
(62)
$\displaystyle-0.4\nu\sum_{k}m_{k}\frac{K_{ak}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\mathbf{v}^{n}_{ak}\cdot\hat{\mathbf{r}}_{ak}\right)\hat{\mathbf{r}}_{ak}D_{ak},$
$\displaystyle\frac{\mathbf{v}^{n+\frac{1}{2}}_{i}-\mathbf{v}^{n}_{i}}{\Delta
t/2}$ $\displaystyle=$
$\displaystyle+0.6\nu\sum_{b}m_{b}\frac{K_{bi}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left(\mathbf{v}^{n+\frac{1}{2}}_{bi}\cdot\hat{\mathbf{r}}_{bi}\right)\hat{\mathbf{r}}_{bi}D_{bi},$
(63)
$\displaystyle+0.4\nu\sum_{b}m_{b}\frac{K_{bi}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left(\mathbf{v}^{n}_{bi}\cdot\hat{\mathbf{r}}_{bi}\right)\hat{\mathbf{r}}_{bi}D_{bi},$

and then

$\displaystyle\frac{\mathbf{v}^{n+1}_{a}-\left[1.4\mathbf{v}^{n+\frac{1}{2}}_{a}-0.4\mathbf{v}^{n}_{a}\right]}{\Delta
t/2}$ $\displaystyle=$
$\displaystyle-0.6\nu\sum_{k}m_{k}\frac{K_{ak}}{\hat{\rho}_{a}\hat{\rho}_{k}}\left(\mathbf{v}^{n+1}_{ak}\cdot\hat{\mathbf{r}}_{ak}\right)\hat{\mathbf{r}}_{ak}D_{ak},$
$\displaystyle\frac{\mathbf{v}^{n+1}_{i}-\left[1.4\mathbf{v}^{n+\frac{1}{2}}_{i}-0.4\mathbf{v}^{n}_{i}\right]}{\Delta
t/2}$ $\displaystyle=$
$\displaystyle+0.6\nu\sum_{b}m_{b}\frac{K_{bi}}{\hat{\rho}_{b}\hat{\rho}_{i}}\left(\mathbf{v}^{n+1}_{bi}\cdot\hat{\mathbf{r}}_{bi}\right)\hat{\mathbf{r}}_{bi}D_{bi}.$

The last terms of Eqs. 62–63 correspond to the explicit drag force involved in
the Forward-Euler scheme. These quantities are computed form the velocities at
the timestep $n$ at the beginning of the scheme, as described in Sec. 4.1.
Then, the successive determination of $\bf{v}^{n+\frac{1}{2}}$ and
$\bf{v}^{n}$ as given by Eqs.
LABEL:eq:Tischer2_scheme1–LABEL:eq:Tischer2_scheme2 consists of two Backward-
Euler steps with step size $\frac{\Delta t}{2}$. They are therefore computed
using our alternative pairwise scheme, until the iterations for each half time
step have converged. Eqs. 62 and 63 concerns the specific case of a linear
drag regime, but this scheme can easily be extended to non-linear drag regimes
as in Sec. 4.2.6.

##  5 Numerical tests

###  5.1 dustybox: Two fluid drag in a periodic box

The dustybox problem presented by Laibe & Price (2011a) and described in
detail in Paper I consists of two fluids in a periodic box moving with a
differential velocity ($\Delta v_{0}=v_{d,0}-v_{g,0}$). This is the only test
where analytic solutions are known for several functional forms corresponding
to non-linear drag regimes (see Laibe & Price 2011a). These represent the
functional forms of the Epstein and Stokes prescription. We thus use the
dustybox problem to benchmark the accuracy of our algorithm for non-linear
drag regimes using both explicit and implicit timestepping.

####  5.1.1 dustybox: setup

We set up constant densities $\hat{\rho}_{\mathrm{g}}$ and
$\hat{\rho}_{\mathrm{d}}$ and gas pressure $P_{\mathrm{g}}$ in a 3D periodic
domain $x,y,z\in[0,1]$ filled by $20^{3}$ gas particles set up on a regular
cubic lattice and $20^{3}$ dust particles shifted by half of the lattice
spacing in each direction (as in Paper I, we verified that the results are
independent of the offset of the dust lattice). The gas sound speed, the gas
and the dust densities are set to unity in code units and no artificial
viscosity is applied. The intrinsic dust volume is neglected by assuming
$\theta=1$.

Simulations have been performed using both the explicit timestepping presented
in Paper I and the implicit pairwise timestepping described in Sec. 4. For the
latter, we verified that both the total linear and angular momentum are
exactly conserved as expected.

####  5.1.2 dustybox: different drag regimes

Figure 1: Dust velocity (solid lines) as a function of time in the dustybox
test, using $2\times 20^{3}$ particles, a dust-to-gas ratio of unity and five
different linear and non-linear drag regimes — quadratic, power-law, linear,
third order expansion and mixed, from top to bottom — compared to the exact
solution for each case (long dashed/red lines). The initial velocities are set
to $v_{\mathrm{d},i}=1$, $v_{\mathrm{g},i}=0$ and the time integration is
performed using using the pairwise implicit treatment described in Sec. 4. The
accuracy ($\lesssim 0.1\%$) of the SPH treatment for dust-gas mixtures is
obtained by using the double-hump cubic kernel.

Fig. 1 shows the results of the dustybox test for the five different regimes
given in Table 1 of Laibe & Price (2011a): linear ($K=K_{0}$), quadratic
($K=K_{0}|\Delta v|$), power-law ($K=K_{0}|\Delta v|^{a}$, with $a=0.4$),
third order expansion ($K=K_{0}[1+a_{3}|\Delta v|^{2}]$, with $a_{3}=0.5$) and
mixed ($K=K_{0}\sqrt{1+a_{2}|\Delta v|^{2}}$, with $a_{2}=5$) where we have
used $K_{0}=1$ in each case. The analytic solutions are reproduced within an
accuracy comprised between $0.1\%$ and $1\%$ in every case — both linear and
non-linear. The implicit scheme was find to converge quickly for this problem,
requiring no more than two iterations at every stage of the evolution in each
case.

The efficiency of the damping for the dustybox problem decreases when the
exponent of the drag regime increases, since $||\Delta v||<1$. On the
contrary, additional non-linear terms give an additional contribution to the
drag for the mixed and third order drag regimes, leading to a differential
velocity that is more efficiently damped compared to the linear case.

###  5.2 dustywave: Sound waves in a dust-gas mixture

The exact solution for linear waves propagating in a dust-gas mixture
(dustywave) was derived in Laibe & Price (2011a) assuming a linear drag
regime. Unfortunately, exact solutions prove difficult to obtain for the case
of non-linear drag. Instead, we have verified simply that our simulations of
the dustywave problem for non-linear drag regimes are converged in both space
and time with an explicit timestepping scheme. We have then used these results
to benchmark our simulations using implicit timestepping.

We run the dustyywave problem for the same five non-linear drag regimes used
above. Strictly speaking, It should be noted that assuming $K_{0}$ is constant
(which we assume for this test problem) corresponds to Epstein and Stokes drag
only to first order for the dustywave problem.

####  5.2.1 dustywave: Setup

The dustywave test is performed in a 1D periodic box, placing equally spaced
particles in the periodic domain $x\in[0,1]$ such that the gas and dust
densities are unity in code units. We do not apply any form of viscosity and
the gas sound speed is set to unity. To remain in the linear acoustic regime,
the relative amplitude of the perturbation of both velocity and density are
set to $10^{-4}$.

####  5.2.2 dustywave: Different drag regimes

Figure 2: Solution of the 1D dustywave problem showing the SPH gas (solid
black lines) and dust (dashed black) velocities using three different drag
regimes: linear (top panel), power-law with an exponent of $0.4$ (center
panel) and quadratic (bottom panel). The damping is strongly reduced for non-
linear drag regimes compared to the linear case. Results are shown after 5
wave periods and the linear case (top panel) may be compared to the exact
analytic solution (red lines). The dust to gas ratio and $K_{0}$ are set to 1
in code units.

Fig. 2 shows the velocity profiles after 5 periods for three drag regimes —
linear, power law and quadratic — in the 1D dustywave problem using $K_{0}=1$
and a dust-to-gas ratio of unity. The solution obtained for the linear drag
regime shows an efficiently damped perturbation at $t=5$, consistent with a
stopping time of order unity. By comparison, the perturbation is only weakly
damped for the case of the power-law drag regime at the same time. The drag is
weaker still in the quadratic regime, for which both the dust and the gas are
mostly decoupled. Indeed, since the drag stopping time is a decreasing
function of the differential velocity for non-linear drag, and the
differential velocity is small with respect to the sound speed, the damping is
inefficient.

We have also performed dustywave simulations using the third order expansion
and mixed drag regimes. However, for these cases, non-linear terms represent
only negligible corrections compared to the linear term, thus giving the same
results as for the linear case.

###  5.3 dustyshock: Two fluid dust-gas shocks

The dustyshock problem (Paper I) is a two fluid version of the standard Sod
(1978) shock tube problem. Results were presented in Paper I using a constant
drag coefficient $K$ and no heat transfer between the gas and the dust phase.
Here we extend the test to non-linear drag regimes. While the evolution during
the transient stage is dependent on the drag regime, the solution during the
stationary stage remains unchanged, being a fixed function of the gas sound
speed and the dust-to-gas ratio. To facilitate the comparison between a
physical Epstein drag and the case of a constant drag coefficient (Paper I),
we fix the ratio $s^{2}/m_{\rm{d}}$ — $s$ being the grain size and
$m_{\rm{d}}$ the grain mass— to unity in code units, such that the Epstein
drag coefficient is unity in regions where $\hat{\rho}_{\mathrm{g}}=1$.

####  5.3.1 dustyshock: Setup

Equal mass particles are placed in the 1D domain $x\in[-0.5,0.5]$, where for
$x<0$ we use $\rho_{\mathrm{g}}=\rho_{\mathrm{d}}=1$,
$v_{\rm{g}}=v_{\rm{d}}=0$ and $P_{\rm{g}}=1$, while for $x>0$
$\rho_{\mathrm{g}}=\rho_{\mathrm{d}}=0.125$, $v_{\rm{g}}=v_{\rm{d}}=0$ and
$P_{\rm{g}}=0.1$. We use an ideal gas equation of state $P=(\gamma-1)\rho u$
with $\gamma=5/3$. Initial particle spacing to the left of the shock in both
fluids is $\Delta x=0.001$ while to the right it is $\Delta x=0.008$, giving
$569$ equal mass particles in each phase. Standard SPH artificial viscosity
and conductivity terms are applied as in Paper I.

####  5.3.2 dustyshock: Different drag regimes

Figure 3: Results of the dustyshock problem with a linear constant drag
coefficient ($K=1$) (left panel) and a non-linear Epstein drag (right panel)
where the drag coefficients are initially the same ahead of the shock. The
dust-to-gas ratio is set to unity. At $t=0.2$, the solutions are in the
transient stage where the analytic solution is not known. As an indication,
the solution for the later stationary stage is shown by the dotted/red lines.
The profiles differ as the damping is less efficient using the non-linear
Epstein regime.

Fig. 3 illustrates how the transient regime of the dustyshock is affected when
treating the drag with an astrophysical prescription where the drag
coefficient depends on the local density and the gas sound speed (right panel)
rather than a constant coefficient (left panel). Specifically, we use the non-
linear Epstein drag regime given by Eq. 39. In this dustyshock test, the non-
linear terms constitute a small correction to the linear Epstein drag regime
given by Eq. 40. Fig. 3 shows that in the case of the Epstein regime, the drag
is less efficient than for the constant coefficient case, leading to a larger
($\sim$ by a factor $7$) differential velocity between the gas and the dust
after $t=0.2$ in code units. The dust velocity profile is also smoother than
in the constant coefficient case. As a result, the kinetic energy is less
efficiently dissipated by the drag, leading to a less sharp peak in the
internal energy of the gas. The density profile of the dust is also closer to
its initial profile behind and ahead of the shock.

###  5.4 dustysedov: Two fluid dust-gas blast wave

The dustysedov problem (Paper I) involves the propagation of a blast wave in
an astrophysical mixture of dust and gas. We adopt physical units for this
problem, assuming a box size of $1$ pc, an ambient sound speed of $2\times
10^{4}$ cm/s and a gas density of $\rho_{0}=6\times 10^{-23}$ g/cm3 the energy
of the blast is $2\times 10^{51}$ erg and time is measured in units of $100$
years, roughly corresponding to a supernova blast wave propagating into the
interstellar medium. We these units, we choose the grain size, $0.1\mu$m and
the dust-to-gas ratio, $0.01$, to be typical of the interstellar medium. In
code units, this corresponds to an initial drag coefficient of $K=1$ outside
the blast radius. As for the dustyshock, we compare the results using a non-
linear Epstein drag prescription with with the constant coefficient case
described in Paper I.

####  5.4.1 dustysedov: Setup

The problem is set up in a 3D periodic box ($x,y,z\in[-0.5,0.5]$), filled by
$50^{3}$ particles for both the gas and the dust. Gas particles are set up on
a regular cubic lattice, with the dust particles also on a cubic lattice but
shifted by half of the lattice step in each direction. For shock-capturing, we
set $\alpha_{\rm SPH}=1$ and $\beta_{\rm SPH}=2$ for the artificial viscosity
terms and $\alpha_{u}=1$ for the artificial conductivity term. An ideal gas
equation of state $P=(\gamma-1)\rho u$ is adopted with $\gamma=5/3$.

The internal energy is distributed of the gas over the particles located
inside a radius $r<r_{\rm{b}}$ where $r_{\rm b}$ is set to 2h (i.e., the
radius of the smoothing kernel which for $50^{3}$ particles and $\eta=1.2$ is
$0.048$). In code units the total blast energy is $E=1$, with
$\hat{\rho}_{\mathrm{g}}=1$ and $\hat{\rho}_{\mathrm{d}}=0.01$. For
$r>r_{\rm{b}}$, the gas sound speed is set to be $2\times 10^{-5}$ in code
units.

####  5.4.2 dustysedov: Different drag regimes

Figure 4: Results of the 3D dustysedov test, showing the density in the gas
(left figure) and dust (right figure) from a Sedov blast wave propagating in
an astrophysical ($1\%$ dust-to-gas ratio) mixture of gas and $0.1\mu$m dust
grains in a $1$ pc box. The drag coefficient is constant ($K=1$, top panels)
or given by the Epstein regime (bottom panels). The low dust-to-gas ratio
means that the gas is only weakly affected by the drag from the dust, and is
thus close to the self-similar Sedov solution (dotted/red line). In the
Epstein case, the drag is much higher inside the blast radius and the dust
particles are efficiently piled up by the passage of the gas over-density.
Figure 5: Cross-section slice showing density in the midplane in the 3D
dustysedov problem, for both the gas (left panel) and the dust (right panel)
at $t=0.1$. Initially, the dust-to-gas ratio is $0.01$ and the drag
coefficient is given by the Epstein regime for grains of $0.1\mu$m in size.
$50^{3}$ SPH particles have been used in each phase.

Figs. 4 and 5 show the evolution of the gas and dust mixture where a constant
drag coefficient is used (top panels of Fig. 4) compared to a drag prescribed
by the Epstein regime (bottom panels of Fig. 4, Fig. 5). The gas profiles are
similar in both cases since the gas is poorly affected by the dust given the
low dust-to-gas ratio. However, the dust density profiles differ, essentially
due to the fact that the drag coefficient scales with the sound speed and is
thus higher in the inner blast region for the Epstein case. Thus, the dust is
efficiently piled up and accumulates in the gas over-density. As a result, the
dust is cleaned up by the gas in the inner regions of the blast, but is more
concentrated (by $\sim 10\%$) close to the gas over-density than for the
constant drag coefficient case.

The results using either explicit or implicit timestepping were found to be
indistinguishable. For the Epstein case, we found that roughly ten iterations
were required for the implicit scheme to converge on this problem.

###  5.5 dustydisc

The dustydisc problem concerns the evolution of a dusty gas mixture in a
protoplanetary disc (see Paper I for details). For our test case, we study how
the dust distribution is affected when considering a general non-linear
Epstein drag instead of the standard linear regime. The results obtained when
implicitly integrating the non-linear drag regime have been found to be
similar to benchmark tests performed with explicit integration.

####  5.5.1 dustydisc: Setup

We setup $10^{5}$ gas particles and $10^{5}$ dust particles in a
$0.01M_{\odot}$ gas disc (with $0.0001M_{\odot}$ of dust) surrounding a
$1M_{\odot}$ star. The disc extends from 10 to 400 AU. Both gas and dust
particles are placed using a Monte-Carlo setup such that the surface density
profiles of both phases are $\Sigma\left(r\right)\propto r^{-1}$. The radial
profile of the gas temperature is taken to be $T\left(r\right)\propto
r^{-0.6}$ with a flaring $H/r=0.05$ at 100 AU. One code unit of time
corresponds to $10^{3}$ yrs.

####  5.5.2 dustydisc: Evolution of the particles

Figure 6: Rendering of the density for the dust of a typical T-Tauri Star
protoplanerary disc using $2\times 10^{5}$ SPH particles, using an explicit
time integration in the linear Epstein regime (left panel) and an implicit
integrator in the full non-linear Epstein drag regimes (right panel).

Fig. 6 shows a face-on view of a protoplanetary disc, integrating the linear
Epstein regime (left panel) and the full non-linear Epstein drag (right
panel). The dust distributions are not found to exhibit significant
discrepancies. In the non-linear drag regime case however, the dust
distribution is slightly smoother since the drag (and thus, the coupling with
the gas phase) is more efficient.

Figure 7: Vertical settling of a dust grain ($1$cm in size) initially located
at $r_{0}=100$ AU and $z_{0}=2$ AU (solid/black), integrating implicitly the
non-linear Epstein regime. SPH results are compared to the explicit
integration of the linear Epstein regime (dashed/red) and the estimation given
by the damped harmonic oscillator approximation (pointed/red). In the full
non-linear drag regime, the settling is more efficient than for the linear
case since the vertical oscillations in the dust motion reaches a fraction
$z_{0}/H$ of the sound speed.

Fig. 7 compares the vertical motion of a dust grain initially located at
$z=z_{0}$ using the linear (explicit integration) and the full non linear
(implicit integration) Epstein regimes. In the full non-linear case, the
settling is more efficient since the vertical differential velocity between
the dust grains and the gas in the mid plane of the disc reaches a fraction
$z_{0}/H$ of the sound speed, meaning that the non-linear terms are no longer
negligible.

##  6 Conclusions

We have extended the SPH formalism for two-fluid dust and gas mixtures
developed in Paper I to handle the drag regimes usually encountered in a large
range of astrophysical contexts. Specifically, our algorithm is now designed
to treat the dynamics of grains surrounded by a dilute medium (Epstein regime)
or dense fluid (Stokes regime), for which the drag force can be either linear
or non-linear with respect to the differential velocity between the gas and
the dust.

Particular attention has been paid to developing an implicit timestepping
scheme to efficiently simulate the case of high drag, extending the scheme
proposed by Monaghan (1997) which we found to be unsatisfactory. We have
presented a new pairwise implicit scheme that, like the Monaghan (1997)
scheme, preserves the exact conservation of linear and angular momentum but,
unlike the Monaghan (1997) scheme, i) provides control over the accuracy of
the iterative procedure and ii) can incorporate non-linear terms for both
Epstein and Stokes drag. We found that when the ratio $r$ between the the
timestep and the drag stopping time is $1\lesssim r\lesssim 1000$, the
implicit timestepping is faster than a standard explicit integration. However,
at higher values of $r$, the algorithm is less efficient.

The accuracy of the generalised algorithm is benchmarked against the suite of
test problems presented in Paper I. In particular, the solutions obtained for
the dustybox problem are compared to their known analytic solutions for a
large range of non-linear drag regimes and the solutions of the dustywave,
dustyshock, dustysedov and dustydisc problems are benchmarked against
converged results obtained with explicit timestepping.

The two key issues addressed in this paper complete the study of our algorithm
developed in Paper I. Our intention is to apply it to various astrophysical
problems involving gas and dust mixtures in star and planet formation. A first
application is given in Ayliffe et al. (2011).

## Acknowledgments

We thank Ben Ayliffe and Matthew Bate and Joe Monaghan for useful discussions
and comments. Figures have been produced using splash (Price, 2007) with the
new giza backend by DP and James Wetter. We are grateful to the Australian
Research Council for funding via Discovery project grant DP1094585.

## References

  * Ayliffe et al. (2011) Ayliffe B., Laibe G., Price D. J., Bate M. R., 2011, MNRAS, submitted
  * Baines et al. (1965) Baines M. J., Williams I. P., Asebiomo A. S., 1965, MNRAS, 130, 63
  * Blum & Wurm (2008) Blum J., Wurm G., 2008, ARA&A, 46, 21
  * Chapman & Cowling (1970) Chapman C., Cowling T., 1970, The mathematical theory of non-uniform gases. Cambridge at the university press
  * Chiang & Youdin (2010) Chiang E., Youdin A. N., 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Fan & Zhu (1998) Fan L.-S., Zhu C., 1998, Principles of Gas-Solid Flows. Cambridge University Press
  * Kwok (1975) Kwok S., 1975, ApJ, 198, 583
  * Laibe & Price (2011b) Laibe G., Price D. J., 2011b, MNRAS, submitted
  * Laibe & Price (2011a) Laibe G., Price D. J., 2011a, MNRAS, in press
  * Monaghan (1997) Monaghan J. J., 1997, Journal of Computational Physics, 138, 801
  * Monaghan (2005) Monaghan J. J., 2005, Reports on Progress in Physics, 68, 1703
  * Monaghan & Kocharyan (1995) Monaghan J. J., Kocharyan A., 1995, Computer Physics Communications, 87, 225
  * Paardekooper & Mellema (2006) Paardekooper S.-J., Mellema G., 2006, A&A, 453, 1129
  * Price (2007) Price D. J., 2007, Publications of the Astronomical Society of Australia, 24, 159
  * Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  * Sod (1978) Sod G. A., 1978, Journal of Computational Physics, 27, 1
  * Stepinski & Valageas (1996) Stepinski T. F., Valageas P., 1996, A&A, 309, 301
  * Whitehouse et al. (2005) Whitehouse S. C., Bate M. R., Monaghan J. J., 2005, MNRAS, 364, 1367