Hostname: page-component-7d684dbfc8-dh8xm Total loading time: 0 Render date: 2023-09-30T02:10:34.821Z Has data issue: false Feature Flags: { "corePageComponentGetUserInfoFromSharedSession": true, "coreDisableEcommerce": false, "coreDisableSocialShare": false, "coreDisableEcommerceForArticlePurchase": false, "coreDisableEcommerceForBookPurchase": false, "coreDisableEcommerceForElementPurchase": false, "coreUseNewShare": true, "useRatesEcommerce": true } hasContentIssue false

Phantom: A Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code for Astrophysics

Published online by Cambridge University Press:  25 September 2018

Daniel J. Price*
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia
James Wurster
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia School of Physics, University of Exeter, Exeter EX4 4QL, UK
Terrence S. Tricco
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia Canadian Institute for Theoretical Astrophysics (CITA), University of Toronto, Toronto, ON M5S 3H8, Canada
Chris Nixon
Affiliation:
Theoretical Astrophysics Group, Department of Physics & Astronomy, University of Leicester, Leicester LE1 7RH, UK
Stéven Toupin
Affiliation:
Institut d’Astronomie et d’Astrophysique (IAA), Université Libre de Bruxelles (ULB), CP226, Boulevard du Triomphe B1050 Brussels, Belgium
Alex Pettitt
Affiliation:
Department of Cosmosciences, Hokkaido University, Sapporo 060-0810, Japan
Conrad Chan
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia
Daniel Mentiplay
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia
Guillaume Laibe
Affiliation:
Centre de Recherche Astrophysique de Lyon, Univ Lyon, ENS de Lyon, CNRS, Saint-Genis-Laval F-69230, France
Simon Glover
Affiliation:
Institut für Theoretische Astrophysik, Zentrum für Astronomie der Universität Heidelberg, D-69120 Heidelberg, Germany
Clare Dobbs
Affiliation:
School of Physics, University of Exeter, Exeter EX4 4QL, UK
Rebecca Nealon
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia Theoretical Astrophysics Group, Department of Physics & Astronomy, University of Leicester, Leicester LE1 7RH, UK
David Liptai
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia
Hauke Worpel
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia AIP Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany
Clément Bonnerot
Affiliation:
Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands
Giovanni Dipierro
Affiliation:
Theoretical Astrophysics Group, Department of Physics & Astronomy, University of Leicester, Leicester LE1 7RH, UK
Giulia Ballabio
Affiliation:
Theoretical Astrophysics Group, Department of Physics & Astronomy, University of Leicester, Leicester LE1 7RH, UK
Enrico Ragusa
Affiliation:
Dipartimento di Fisica, Università Degli Studi di Milano, Via Celoria 16, Milano 20133, Italy
Christoph Federrath
Affiliation:
Research School of Astronomy and Astrophysics, Australian National University, Canberra ACT 2611, Australia
Roberto Iaconi
Affiliation:
Department of Physics and Astronomy, Macquarie University, 2109 Sydney, Australia
Thomas Reichardt
Affiliation:
Department of Physics and Astronomy, Macquarie University, 2109 Sydney, Australia
Duncan Forgan
Affiliation:
St Andrews Centre for Exoplanet Science and School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK
Mark Hutchison
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia Physikalisches Institut, Universität Bern, Gesellschaftstrasse 6, 3012 Bern, Switzerland Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland
Thomas Constantino
Affiliation:
School of Physics, University of Exeter, Exeter EX4 4QL, UK
Ben Ayliffe
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia Met Office, FitzRoy Road, Exeter EX1 3PB, UK
Kieran Hirsh
Affiliation:
Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Vic. 3800, Australia
Giuseppe Lodato
Affiliation:
Dipartimento di Fisica, Università Degli Studi di Milano, Via Celoria 16, Milano 20133, Italy
Rights & Permissions [Opens in a new window]

Abstract

We present Phantom, a fast, parallel, modular, and low-memory smoothed particle hydrodynamics and magnetohydrodynamics code developed over the last decade for astrophysical applications in three dimensions. The code has been developed with a focus on stellar, galactic, planetary, and high energy astrophysics, and has already been used widely for studies of accretion discs and turbulence, from the birth of planets to how black holes accrete. Here we describe and test the core algorithms as well as modules for magnetohydrodynamics, self-gravity, sink particles, dust–gas mixtures, H2 chemistry, physical viscosity, external forces including numerous galactic potentials, Lense–Thirring precession, Poynting–Robertson drag, and stochastic turbulent driving. Phantom is hereby made publicly available.

Type
Research Article
Copyright
Copyright © Astronomical Society of Australia 2018 

1 INTRODUCTION

Numerical simulations are the ‘third pillar’ of astrophysics, standing alongside observations and analytic theory. Since it is difficult to perform laboratory experiments in the relevant physical regimes and over the correct range of length and timescales involved in most astrophysical problems, we turn instead to ‘numerical experiments’ in the computer for understanding and insight. As algorithms and simulation codes become ever more sophisticated, the public availability of simulation codes has become crucial to ensure that these experiments can be both verified and reproduced.

Phantom is a smoothed particle hydrodynamics (SPH) code developed over the last decade. It has been used widely for studies of turbulence (e.g. Kitsionas et al. Reference Kitsionas2009; Price & Federrath Reference Price and Federrath2010; Price, Federrath, & Brunt Reference Price, Federrath and Brunt2011), accretion (e.g. Lodato & Price Reference Lodato and Price2010; Nixon, King, & Price Reference Nixon, King and Price2012a; Rosotti, Lodato, & Price Reference Rosotti, Lodato and Price2012), star formation including non-ideal magnetohydrodynamics (MHD) (e.g. Wurster et al. Reference Wurster, Price and Bate2016, Wurster, Price, & Bate Reference Wurster, Price and Bate2017), star cluster formation (Liptai et al. Reference Liptai, Price, Wurster and Bate2017), and for studies of the Galaxy (Pettitt et al. Reference Pettitt, Dobbs, Acreman and Price2014; Dobbs et al. Reference Dobbs, Price, Pettitt, Bate and Tricco2016), as well as for simulating dust–gas mixtures (e.g. Dipierro et al. Reference Dipierro, Price, Laibe, Hirsh and Cerioli2015; Ragusa et al. Reference Ragusa, Dipierro, Lodato, Laibe and Price2017; Tricco, Price, & Laibe Reference Tricco, Price and Laibe2017). Although the initial applications and some details of the basic algorithm were described in Price & Federrath (Reference Price and Federrath2010), Lodato & Price (Reference Lodato and Price2010), and Price (Reference Price2012a), the code itself has never been described in detail and, until now, has remained closed source.

One of the initial design goals of Phantom was to have a low memory footprint. A secondary motivation was the need for a public SPH code that is not primarily focused on cosmology, as in the highly successful Gadget code (Springel, Yoshida, & White Reference Springel, Yoshida and White2001; Springel Reference Springel2005). The needs of different communities produce rather different outcomes in the code. For cosmology, the main focus is on simulating the gravitational collapse of dark matter in large volumes of the universe, with gas having only a secondary effect. This is reflected in the ability of the public Gadget-2 code to scale to exceedingly large numbers of dark matter particles, yet neglecting elements of the core SPH algorithm that are essential for stellar and planetary problems, such as the Morris & Monaghan (Reference Morris and Monaghan1997) artificial viscosity switch [c.f. the debate between Bauer & Springel (Reference Bauer and Springel2012) and Price (Reference Price2012b)], the ability to use a spatially variable gravitational force softening (Bate & Burkert Reference Bate and Burkert1997; Price & Monaghan Reference Price and Monaghan2007) or any kind of artificial conductivity, necessary for the correct treatment of contact discontinuities (Chow & Monaghan Reference Chow and Monaghan1997; Price & Monaghan Reference Price and Monaghan2005; Rosswog & Price Reference Rosswog and Price2007; Price Reference Price2008). Almost all of these have since been implemented in development versions of Gadget-3 (e.g. Iannuzzi & Dolag Reference Iannuzzi and Dolag2011; Beck et al. Reference Beck2016; see recent comparison project by Sembolini et al. Reference Sembolini2016) but remain unavailable or unused in the public version. Likewise, the implementation of dust, non-ideal MHD, and other physics relevant to star and planet formation is unlikely to be high priority in a code designed for studying cosmology or galaxy formation.

Similarly, the sphng code (Benz et al. Reference Benz, Cameron, Press and Bowers1990; Bate Reference Bate1995) has been a workhorse for our group for simulations of star formation (e.g. Price & Bate Reference Price and Bate2007, Reference Price and Bate2009; Price, Tricco, & Bate Reference Price, Tricco and Bate2012; Lewis, Bate, & Price Reference Lewis, Bate and Price2015) and accretion discs (e.g. Lodato & Rice Reference Lodato and Rice2004; Cossins, Lodato, & Clarke Reference Cossins, Lodato and Clarke2009), contains a rich array of physics necessary for star and planet formation including all of the above algorithms, but the legacy nature of the code makes it difficult to modify or debug and there are no plans to make it public.

Gasoline (Wadsley, Stadel, & Quinn Reference Wadsley, Stadel and Quinn2004) is another code that has been widely and successfully used for galaxy formation simulations, with its successor, Gasoline 2 (Wadsley et al. Reference Wadsley, Keller and Quinn2017), recently publicly released. Hubber et al. (Reference Hubber, Batty, McLeod and Whitworth2011) have developed Seren with similar goals to Phantom, focused on star cluster simulations. Seren thus presents more advanced N-body algorithms compared to what is in Phantom but does not yet include magnetic fields, dust, or H2 chemistry.

A third motivation was the need to distinguish between the ‘high performance’ code used for 3D simulations from simpler codes used to develop and test algorithms, such as our already-public ndspmhd code (Price Reference Price2012a). Phantom is designed to ‘take what works and make it fast’, rather than containing options for every possible variation on the SPH algorithm. Obsolete options are actively deleted.

The initial release of Phantom has been developed with a focus on stellar, planetary, and Galactic astrophysics, as well as accretion discs. In this first paper, coinciding with the first stable public release, we describe and validate the core algorithms as well as some example applications. Various novel aspects and optimisation strategies are also presented. This paper is an attempt to document in detail what is currently available in the code, which include modules for MHD, dust–gas mixtures, self-gravity, and a range of other physics. The paper is also designed to serve as guide to the correct use of the various algorithms. Stable releases of Phantom are posted on the webFootnote 1, while the development version and wiki documentation are available on the Bitbucket platformFootnote 2.

The paper is organised as follows: We describe the numerical methods in Section 2 with corresponding numerical tests in Section 5. We cover SPH basics (Section 2.1), our implementation of hydrodynamics (Sections 2.2 and 5.1), the timestepping algorithm (Section 2.3), external forces (Sections 2.4 and 5.2), turbulent forcing (Sections 2.5 and 6.1), accretion disc viscosity (Sections 2.6 and 5.3), Navier–Stokes viscosity (Sections 2.7 and 5.4), sink particles (Sections 2.8 and 5.5), stellar physics (Section 2.9), MHD (Sections 2.10 and 5.6), non-ideal MHD (Sections 2.11 and 5.7), self-gravity (Sections 2.12 and 5.8), dust–gas mixtures (Sections 2.13 and 5.9), ISM chemistry and cooling (Sections 2.14 and 5.10), and particle injection (Section 2.15). We present the algorithms for generating initial conditions in Section 3. Our approach to software engineering is described in Section 4. We give five examples of recent applications highlighting different aspects of Phantom in Section 6. We summarise in Section 7.

2 NUMERICAL METHOD

Phantom is based on the SPH technique, invented by Lucy (Reference Lucy1977) and Gingold & Monaghan (Reference Gingold and Monaghan1977) and the subject of numerous reviews (Benz Reference Benz and Buchler1990; Monaghan Reference Monaghan1992, Reference Monaghan2005, Reference Monaghan2012; Rosswog Reference Rosswog2009; Springel Reference Springel2010; Price Reference Price2012a).

In the following, we adopt the convention that a, b, and c refer to particle indices; i, j, and k refer to vector or tensor indices; and n and m refer to indexing of nodes in the treecode.

2.1. Fundamentals

2.1.1. Lagrangian hydrodynamics

SPH solves the equations of hydrodynamics in Lagrangian form. The fluid is discretised onto a set of ‘particles’ of mass m that are moved with the local fluid velocity ${\bm v}$. Hence, the two basic equations common to all physics in Phantom are

(1) \begin{align} \frac{{\rm d}{\bm r}}{{\rm d}t} = {\bm v}, \end{align}
(2) \begin{align} \frac{{\rm d}\rho }{{\rm d}t} = -\rho (\nabla \cdot {\bm v}), \end{align}

where ${\bm r}$ is the particle position and ρ is the density. These equations use the Lagrangian time derivative, ${\rm d}/{\rm d}t \equiv \partial / \partial t + {\bm v} \cdot \nabla$, and are the Lagrangian update of the particle position and the continuity equation (expressing the conservation of mass), respectively.

2.1.2. Conservation of mass in SPH

The density is computed in Phantom using the usual SPH density sum

(3) $$\begin{equation} \rho _{a} = \sum _{b} m_{b} W(\vert {\bm r}_{a} - {\bm r}_{b}\vert , h_{a}), \end{equation}$$

where a and b are particle labels, m is the mass of the particle, W is the smoothing kernel, h is the smoothing length, and the sum is over neighbouring particles (i.e. those within R kernh, where R kern is the dimensionless cut-off radius of the smoothing kernel). Taking the Lagrangian time derivative of (3), one obtains the discrete form of (2) in SPH

(4) $$\begin{equation} \frac{{\rm d}\rho _{a}}{{\rm d}t} = \frac{1}{\Omega _{a}} \sum _{b} m_{b} ({\bm v}_{a} - {\bm v}_{b}) \cdot \nabla _{a} W_{ab} (h_{a}), \end{equation}$$

where $W_{ab}(h_{a})\equiv W(\vert {\bm r}_{a} - {\bm r}_{b}\vert , h_{a})$ and Ωa is a term related to the gradient of the smoothing length (Springel & Hernquist Reference Springel and Hernquist2002; Monaghan Reference Monaghan2002) given by

(5) $$\begin{equation} \Omega _{a} \equiv 1 - \frac{\partial h_{a}}{\partial \rho _{a}}\sum _{b} m_{b} \frac{\partial W_{ab} (h_{a})}{\partial h_{a}}. \end{equation}$$

Equation (4) is not used directly to compute the density in Phantom, since evaluating (3) provides a time-independent solution to (2) (see e.g. Monaghan Reference Monaghan1992; Price Reference Price2012a for details). The time-dependent version (4) is equivalent to (3) up to a boundary term (see Price Reference Price2008) but is only used in Phantom to predict the smoothing length at the next timestep in order to reduce the number of iterations required to evaluate the density (see below).

Since (3)–(5) all depend on the kernel evaluated on neighbours within R kern times h a, all three of these summations may be computed simultaneously using a single loop over the same set of neighbours. Details of the neighbour finding procedure are given in Section 2.1.7.

2.1.3. Setting the smoothing length

The smoothing length itself is specified as a function of the particle number density, n, via

(6) $$\begin{equation} h_{a} = h_{\rm fact} n_{a}^{-1/3} = h_{\rm fact} \left( \frac{m_{a}}{\rho _{a}} \right)^{1/3}, \end{equation}$$

where h fact is the proportionality factor specifying the smoothing length in terms of the mean local particle spacing and the second equality holds only for equal mass particles, which are enforced in Phantom. The restriction to equal mass particles means that the resolution strictly follows mass, which may be restrictive for problems involving large density contrasts (e.g. Hutchison et al. Reference Hutchison, Price, Laibe and Maddison2016). However, our view is that the potential pitfalls of unequal mass particles (see e.g. Monaghan & Price Reference Monaghan and Price2006) are currently too great to allow for a robust implementation in a public code.

As described in Price (Reference Price2012a), the proportionality constant h fact can be related to the mean neighbour number according to

(7) $$\begin{equation} \overline{N}_{\rm neigh} = \frac{4}{3} \pi (R_{\rm kern} h_{\rm fact})^{3}, \end{equation}$$

however, this is only equal to the actual neighbour number for particles in a uniform density distribution (more specifically, for a density distribution with no second derivative), meaning that the actual neighbour number varies. The default setting for h fact is 1.2, corresponding to an average of 57.9 neighbours for a kernel truncated at 2h (i.e. for R kern = 2) in three dimensions. Table 1 lists the settings recommended for different choices of kernel. The derivative required in (5) is given by

(8) $$\begin{equation} \frac{\partial h_{a}}{\partial \rho _{a}} = -\frac{3 h_{a}}{\rho _{a}}. \end{equation}$$

Table 1. Compact support radii, variance, standard deviation, recommended ranges of h fact, and recommended default h fact settings (hd fact) for the kernel functions available in Phantom.

2.1.4. Iterations for h and ρ

The mutual dependence of ρ and h means that a rootfinding procedure is necessary to solve both (3) and (6) simultaneously. The procedure implemented in Phantom follows Price & Monaghan (Reference Price and Monaghan2004b, Reference Price and Monaghan2007), solving, for each particle, the equation

(9) $$\begin{equation} f(h_{a}) = \rho _{\rm sum} (h_{a}) - \rho (h_{a}) = 0, \end{equation}$$

where ρsum is the density computed from (3) and

(10) $$\begin{equation} \rho (h_{a}) = m_{a} (h_{\rm fact}/h_{a})^{3}, \end{equation}$$

from (6). Equation (9) is solved with Newton–Raphson iterations:

(11) $$\begin{equation} h_{a, {\rm new}} = h_{a} - \frac{f(h_{a})}{f^{\prime }(h_{a})}, \end{equation}$$

where the derivative is given by

(12) \begin{align} f^{\prime }(h_{a}) = \sum _{b} m_{b} \frac{\partial W_{ab} (h_{a})}{\partial h_{a}} - \frac{\partial \rho _{a}}{\partial h_{a}} = -\frac{3\rho _{a}}{h_{a}} \Omega _{a}. \end{align}

The iterations proceed until |h a, newh a|/h a, 0 < εh, where h a, 0 is the smoothing length of particle a at the start of the iteration procedure and εh is the tolerance. The convergence with Newton–Raphson is fast, with a quadratic reduction in the error at each iteration, meaning that no more than 2–3 iterations are required even with a rapidly changing density field. We avoid further iterations by predicting the smoothing length from the previous timestep according to

(13) $$\begin{equation} h^{0}_{a} = h_{a} + \Delta t \frac{{\rm d}h_{a}}{{\rm d}t} = h_{a} + \Delta t \frac{\partial h_{a}}{\partial \rho _{a}} \frac{{\rm d}\rho _{a}}{{\rm d} t}, \end{equation}$$

where dρa/dt is evaluated from (4).

Since h and ρ are mutually dependent, we store only the smoothing length, from which the density can be obtained at any time via a function call evaluating ρ(h). The default value of εh is 10−4 so that h and ρ can be used interchangeably. Setting a small tolerance does not significantly change the computational cost, as the iterations quickly fall below a tolerance of ‘one neighbour’ according to (7), so any iterations beyond this refer to loops over the same set of neighbours which can be efficiently cached. However, it is important that the tolerance may be enforced to arbitrary precision rather than being an integer as implemented in the public version of Gadget, since (9) expresses a mathematical relationship between h and ρ that is assumed throughout the derivation of the SPH algorithm (see discussion in Price Reference Price2012a). The precision to which this is enforced places a lower limit on the total energy conservation. Fortunately, floating point neighbour numbers are now default in most Gadget-3 variants also.

2.1.5. Kernel functions

We write the kernel function in the form

(14) $$\begin{equation} W_{ab}(r,h) \equiv \frac{C_{\rm norm}}{h^{3}} f(q), \end{equation}$$

where C norm is a normalisation constant, the factor of h 3 gives the dimensions of inverse volume, and f(q) is a dimensionless function of $q \equiv \vert {\bm r}_{a} - {\bm r}_{b} \vert / h$. Various relations for kernels in this form are given in Morris (Reference Morris1996a) and in Appendix B of Price (Reference Price2010). Those used in Phantom are the kernel gradient

(15) $$\begin{equation} \nabla _{a} W_{ab} = \hat{\bm r}_{ab} F_{ab},\textrm { where } F_{ab} \equiv \frac{C_{\rm norm}}{h^{4}} f^{\prime }(q), \end{equation}$$

and the derivative of the kernel with respect to h,

(16) $$\begin{equation} \frac{\partial W_{ab}(r,h)}{\partial h} = - \frac{C_{\rm norm}}{h^{4}} \left[ 3 f(q) + qf^{\prime }(q) \right]. \end{equation}$$

Notice that the ∂W/∂h term in particular can be evaluated simply from the functions needed to compute the density and kernel gradient and hence does not need to be derived separately if a different kernel is used.

2.1.6. Choice of smoothing kernel

The default kernel function in SPH for the last 30 yr (since Monaghan & Lattanzio Reference Monaghan and Lattanzio1985) has been the M 4 cubic spline from the Schoenberg (Reference Schoenberg1946) B-spline family, given by

(17) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{ll}1 - \displaystyle\frac{3}{2} q^{2} + \displaystyle\frac{3}{4} q^{3}, & 0 \le q < 1; \\[3pt] \displaystyle\frac{1}{4}(2-q)^3, & 1 \le q < 2; \\[3pt] 0. & q \ge 2, \end{array} \right. \end{equation}$$

where the normalisation constant C norm = 1/π in 3D and the compact support of the function implies that R kern = 2. While the cubic spline kernel is satisfactory for many applications, it is not always the best choice. Most SPH kernels are based on approximating the Gaussian, but with compact support to avoid the $\mathcal {O}(N^{2})$ computational cost. Convergence in SPH is guaranteed to be second order (∝h 2) to the degree that the finite summations over neighbouring particles approximate integrals (e.g. Monaghan Reference Monaghan1992, Reference Monaghan2005; Price Reference Price2012a). Hence, the choice of kernel and the effect that a given kernel has on the particle distribution are important considerations.

In general, more accurate results will be obtained with a kernel with a larger compact support radius, since it will better approximate the Gaussian which has excellent convergence and stability properties (Morris Reference Morris1996a; Price Reference Price2012a; Dehnen & Aly Reference Dehnen and Aly2012). However, care is required. One should not simply increase h fact for the cubic spline kernel because even though this implies more neighbours [via (7)], it increases the resolution length. For the B-splines, it also leads to the onset of the ‘pairing instability’ where the particle distribution becomes unstable to transverse modes, leading to particles forming close pairs (Thomas & Couchman Reference Thomas and Couchman1992; Morris Reference Morris1996a, Reference Morris1996b; Børve, Omang, & Trulsen Reference Børve, Omang and Trulsen2004; Price Reference Price2012a; Dehnen & Aly Reference Dehnen and Aly2012). This is the motivation of our default choice of h fact = 1.2 for the cubic spline kernel, since it is just short of the maximum neighbour number that can be used while remaining stable to the pairing instability.

A better approach to reducing kernel bias is to keep the same resolution lengthFootnote 3 but to use a kernel that has a larger compact support radius. The traditional approach (e.g. Morris Reference Morris1996a, Reference Morris1996b; Børve et al. Reference Børve, Omang and Trulsen2004; Price Reference Price2012a) has been to use the higher kernels in the B-spline series, i.e. the M 5 quartic which extends to 2.5h

(18) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{ll} \left(\displaystyle\frac{5}{2} -q\right)^4 - 5\left(\displaystyle\frac{3}{2} -q\right)^4 & 0 \le q < \displaystyle\frac{1}{2}, \\[6pt] \quad + 10\left(\displaystyle\frac{1}{2}-q\right)^4,\\[6pt] \left(\displaystyle\frac{5}{2} -q\right)^4 - 5\left(\displaystyle\frac{3}{2} -q\right)^4, & \displaystyle\frac{1}{2} \le q < \displaystyle\frac{3}{2}, \\[6pt] \left(\displaystyle\frac{5}{2} -q\right)^4, & \displaystyle\frac{3}{2} \le q < \displaystyle\frac{5}{2}, \\[3pt] 0, & q \ge \frac{5}{2}, \end{array} \right. \end{equation}$$

where C norm = 1/(20π), and the M 6 quintic extending to 3h,

(19) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{ll}(3-q)^5 - 6(2-q)^5 + 15(1-q)^5, & 0 \le q < 1, \\ (3-q)^5 - 6(2-q)^5, & 1 \le q < 2, \\[3pt] (3-q)^5, & 2 \le q < 3, \\[3pt] 0, & q \ge 3, \end{array} \right. \end{equation}$$

where C norm = 1/(120π) in 3D. The quintic in particular gives results virtually indistinguishable from the Gaussian for most problems.

Recently, there has been tremendous interest in the use of the Wendland (Reference Wendland1995) kernels, particularly since Dehnen & Aly (Reference Dehnen and Aly2012) showed that they are stable to the pairing instability at all neighbour numbers despite having a Gaussian-like shape and compact support. These functions are constructed as the unique polynomial functions with compact support but with a positive Fourier transform, which turns out to be a necessary condition for stability against the pairing instability (Dehnen & Aly Reference Dehnen and Aly2012). The 3D Wendland kernels scaled to a radius of 2h are given by C 2:

(20) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}\left(1 - \displaystyle\frac{q}{2} \right)^{4} \left(2 q + 1\right), & q < 2, \\[3pt] 0, & q \ge 2, \end{array}\right. \end{equation}$$

where C norm = 21/(16π), the C 4 kernel

(21) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}\left(1 - \displaystyle\frac{q}{2}\right)^{6} \left(\displaystyle\frac{35 q^{2}}{12} + 3 q + 1\right), & q < 2, \\[3pt] 0, & q\ge 2, \end{array}\right. \end{equation}$$

where C norm = 495/(256π), and the C 6 kernel

(22) $$\begin{equation} f(q) = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}\left(1 - \displaystyle\frac{q}{2}\right)^{8} \left(4 q^{3} + \displaystyle\frac{25 q^{2}}{4} + 4 q + 1\right), & q < 2, \\[3pt] 0, & q \ge 2, \end{array}\right. \end{equation}$$

where C norm = 1365/(512π). Figure 1 graphs f(q) and its first and second derivative for each of the kernels available in Phantom.

Figure 1. Smoothing kernels available in Phantom (solid lines) together with their first (dashed lines) and second (dotted lines) derivatives. Wendland kernels in Phantom (bottom row) are given compact support radii of 2, whereas the B-spline kernels (top row) adopt the traditional practice where the support radius increases by 0.5. Thus, use of alternative kernels requires adjustment of h fact, the ratio of smoothing length to particle spacing (see Table 1).

Several authors have argued for use of the Wendland kernels by default. For example, Rosswog (Reference Rosswog2015) found best results on simple test problems using the C 6 Wendland kernel. However, ‘best’ in that case implied using an average of 356 neighbours in 3D (i.e. h fact = 2.2 with R kern = 2.0) which is a factor of 6 more expensive than the standard approach. Similarly, Hu et al. (Reference Hu, Naab, Walch, Moster and Oser2014) recommend the C 4 kernel with 200 neighbours which is 3.5 times more expensive. The large number of neighbours are needed because the Wendland kernels are always worse than the B-splines for a given number of neighbours due to the positive Fourier transform, meaning that the kernel bias (related to the Fourier transform) is always positive where the B-spline errors oscillate around zero (Dehnen & Aly Reference Dehnen and Aly2012). Hence, whether or not this additional cost is worthwhile depends on the application. A more comprehensive analysis would be valuable here, as the ‘best’ choice of kernel remains an open question (see also the kernels proposed by Cabezón, García-Senz, & Relaño Reference Cabezón, García-Senz and Relaño2008; García-Senz et al. Reference García-Senz, Cabezón, Escartín and Ebinger2014). An even broader question regards the kernel used for dissipation terms, for gravitational force softening and for drag in two-fluid applications (discussed further in Section 2.13). For example, Laibe & Price (Reference Laibe and Price2012a) found that double-hump-shaped kernels led to more than an order of magnitude improvement in accuracy when used for drag terms.

A simple and practical approach to checking that kernel bias does not affect the solution that we have used and advocate when using Phantom is to first attempt a simulation with the cubic spline, but then to check the results with a low resolution calculation using the quintic kernel. If the results are identical, then it indicates that the kernel bias is not important, but if not, then use of smoother but costlier kernels such as M 6 or C 6 may be warranted. Wendland kernels are mainly useful for preventing the pairing instability and are necessary if one desires to employ a large number of neighbours.

2.1.7. Neighbour finding

Finding neighbours is the main computational expense to any SPH code. Earlier versions of Phantom contained three different options for neighbour-finding: a Cartesian grid, a cylindrical grid, and a kd-tree. This was because we wrote the code originally with non-self-gravitating problems in mind, for which the overhead associated with a treecode is unnecessary. Since the implementation of self-gravity in Phantom the kd-tree has become the default, and is now sufficiently well optimised that the fixed-grid modules are more efficient only for simulations that do not employ either self-gravity or individual particle timesteps, which are rare in astrophysics.

A key optimisation strategy employed in Phantom is to perform the neighbour search for groups of particles. The results of this search (i.e. positions of all trial neighbours) are then cached and used to check for neighbours for individual particles in the group. Our kd-tree algorithm closely follows Gafton & Rosswog (Reference Gafton and Rosswog2011), splitting the particles recursively based on the centre of mass and bisecting the longest axis at each level (Figure 2). The tree build is refined until a cell contains less than N min particles, which is then referred to as a ‘leaf node’. By default, N min = 10. The neighbour search is then performed once for each leaf node. Further details are given in Appendix A.3.1.

Figure 2. Example of the kd-tree build. For illustrative purposes only, we have constructed a 2D version of the tree on the projected particle distribution in the xy plane of the particle distribution from a polytrope test with 13 115 particles. Each level of the tree recursively splits the particle distribution in half, bisecting the longest axis at the centre of mass until the number of particles in a given cell is <N min. For clarity, we have used N min = 100 in the above example, while N min = 10 by default.

2.2. Hydrodynamics

2.2.1. Compressible hydrodynamics

The equations of compressible hydrodynamics are solved in the form

(23) $$\begin{eqnarray} \frac{{\rm d}{\bm v}}{{\rm d}t} &=& -\frac{\nabla P}{\rho } + \Pi _{\rm shock} + {\bm a}_{\rm ext}({\bm r}, t) \nonumber \\ && +\, {\bm a}_{\rm sink-gas} + {\bm a}_{\rm selfgrav}, \end{eqnarray}$$
(24) \begin{align} \frac{{\rm d}u}{{\rm d}t} = -\frac{P}{\rho } \left(\nabla \cdot {\bm v}\right) + \Lambda _{\rm shock} - \frac{\Lambda _{\rm cool}}{\rho }, \end{align}

where P is the pressure, u is the specific internal energy, ${\bm a}_{\rm ext}$, ${\bm a}_{\rm sink-gas}$ and ${\bm a}_{\rm selfgrav}$ refer to (optional) accelerations from ‘external’ or ‘body’ forces (Section 2.4), sink particles (Section 2.8), and self-gravity (Section 2.12), respectively. Πshock and Λshock are dissipation terms required to give the correct entropy increase at a shock front, and Λcool is a cooling term.

2.2.2. Equation of state

The equation set is closed by an equation of state (EOS) relating the pressure to the density and/or internal energy. For an ideal gas, this is given by

(25) $$\begin{equation} P = (\gamma - 1)\rho u, \end{equation}$$

where γ is the adiabatic index and the sound speed c s is given by

(26) $$\begin{equation} c_{\rm s} = \sqrt{\frac{\gamma P}{\rho }}. \end{equation}$$

The internal energy, u, can be related to the gas temperature, T, using

(27) $$\begin{equation} P = \frac{\rho k_{\rm B} T}{\mu m_{\rm H}}, \end{equation}$$

giving

(28) $$\begin{equation} T = \frac{\mu m_{\rm H}}{k_{\rm B}} (\gamma -1) u, \end{equation}$$

where k B is Boltzmann’s constant, μ is the mean molecular weight, and m H is the mass of a hydrogen atom. Thus, to infer the temperature, one needs to specify a composition, but only the internal energy affects the gas dynamics. Equation (25) with γ = 5/3 is the default EOS in Phantom.

In the case where shocks are assumed to radiate away all of the heat generated at the shock front (i.e. Λshock = 0) and there is no cooling (Λcool = 0), (24) becomes simply, using (2)

(29) $$\begin{equation} \frac{{\rm d} u}{{\rm d} t} = \frac{P}{\rho ^{2}} \frac{{\rm d}\rho }{{\rm d}t}, \end{equation}$$

which, using (25) can be integrated to give

(30) $$\begin{equation} P = K \rho ^{\gamma }, \end{equation}$$

where K is the polytropic constant. Even more simply, in the case where the temperature is assumed constant, or prescribed as a function of position, the EOS is simply

(31) $$\begin{equation} P = c_{\rm s}^{2} \rho . \end{equation}$$

In both of these cases, (30) and (31), the internal energy does not need to be stored. In this case, the temperature is effectively set by the value of K (and the density if γ ≠ 1). Specifically,

(32) $$\begin{equation} T = \frac{\mu m_{\rm H}}{k_{\rm B}} K \rho ^{\gamma - 1}. \end{equation}$$

2.2.3. Code units

For pure hydrodynamics, physical units are irrelevant to the numerical results since (1)–(2) and (23)–(24) are scale free to all but the Mach number. Hence, setting physical units is only useful when comparing simulations with Nature, when physical heating or cooling rates are applied via (24), or when one wishes to interpret the results in terms of temperature using (28) or (32).

In the case where gravitational forces are applied, either using an external force (Section 2.4) or using self-gravity (Section 2.12), we adopt the standard procedure of transforming units such that G = 1 in code units, i.e.

(33) $$\begin{equation} u_{\rm time} = \sqrt{\frac{u_{\rm dist}^{3}}{{\rm G}u_{\rm mass}}}, \end{equation}$$

where u time, u dist, and u mass are the units of time, length, and mass, respectively. Additional constraints apply when using relativistic terms (Section 2.4.5) or magnetic fields (Section 2.10.3).

2.2.4. Equation of motion in SPH

We follow the variable smoothing length formulation described by Price (Reference Price2012a), Price & Federrath (Reference Price and Federrath2010), and Lodato & Price (Reference Lodato and Price2010). We discretise (23) using

(34) $$\begin{eqnarray} \frac{{\rm d}{\bm v}_{a}}{{\rm d}t} &= & - \sum _{b} m_{b} \left[ \frac{P_{a} + q^{a}_{ab}}{\rho _{a}^{2}\Omega _{a}}\nabla _a W_{ab}(h_{a}) + \frac{P_{b}+q^{b}_{ab}}{\rho _{b}^{2}\Omega _{b}}\nabla _a W_{ab}(h_{b}) \right] \nonumber \\ && +\, {\bm a}_{\rm ext} ({\bm r}_{a}, t) + {\bm a}^{a}_{\rm sink-gas} + {\bm a}^{a}_{\rm selfgrav}, \end{eqnarray}$$

where the qaab and qbab terms represent the artificial viscosity (discussed in Section 2.2.7).

2.2.5. Internal energy equation

The internal energy equation (24) is discretised using the time derivative of the density sum (c.f. 29), which from (4) gives

(35) $$\begin{equation} \frac{{\rm d} u_{a}}{{\rm d} t} = \frac{P_{a}}{\rho _{a}^{2}\Omega _{a}} \sum _{b} m_{b} {\bm v}_{ab} \cdot \nabla _{a} W_{ab} (h_{a}) + \Lambda _{\rm shock} - \frac{\Lambda _{\rm cool}}{\rho }, \end{equation}$$

where ${\bm v}_{ab} \equiv {\bm v}_a - {\bm v}_b$. In the variational formulation of SPH (e.g. Price Reference Price2012a), this expression is used as a constraint to derive (34), which guarantees both the conservation of energy and entropy (the latter in the absence of dissipation terms). The shock capturing terms in the internal energy equation are discussed below.

By default, we assume an adiabatic gas, meaning that PdV work and shock heating terms contribute to the thermal energy of the gas, no energy is radiated to the environment, and total energy is conserved. To approximate a radiative gas, one may set one or both of these terms to zero. Neglecting the shock heating term, Λshock, gives an approximation equivalent to a polytropic EOS (30), as described in Section 2.2.2. Setting both shock and work contributions to zero implies that du/dt = 0, meaning that each particle will simply retain its initial temperature.

2.2.6. Conservation of energy in SPH

Does evolving the internal energy equation imply that total energy is not conserved? Wrong! Total energy in SPH, for the case of hydrodynamics, is given by

(36) $$\begin{equation} E = \sum _a m_a \left(\frac{1}{2} v_a^2 + u_a\right). \end{equation}$$

Taking the (Lagrangian) time derivative, we find that conservation of energy corresponds to

(37) $$\begin{equation} \frac{{\rm d}E}{{\rm d}t} = \sum _{a} m_{a} \left( {\bm v}_{a} \cdot \frac{{\rm d} {\bm v}_{a}}{{\rm d} t} + \frac{{\rm d} u_{a}}{{\rm d} t} \right) = 0. \end{equation}$$

Inserting our expressions (34) and (35), and neglecting for the moment dissipative terms and external forces, we find

(38) $$\begin{eqnarray} \frac{{\rm d}E}{{\rm d}t} &=& -\sum _{a} \sum _{b} m_{a} m_{b} \left[ \frac{P_{a}\bm {v}_b}{\rho _{a}^{2}\Omega _{a}}\cdot \nabla _a W_{ab}(h_{a}) \right. \nonumber \\ && \left. + \frac{P_{b}\bm {v}_a}{\rho _{b}^{2}\Omega _{b}}\cdot \nabla _a W_{ab}(h_{b}) \right] = 0. \end{eqnarray}$$

The double summation on the right-hand side equals zero because the kernel gradient, and hence the overall sum, is antisymmetric. That is, ∇aW ab = −∇bW ba. This means one can relabel the summation indices arbitrarily in one half of the sum, and add it to one half of the original sum to give zero. One may straightforwardly verify that this remains true when one includes the dissipative terms (see below).

This means that even though we employ the internal energy equation, total energy remains conserved to machine precision in the spatial discretisation. That is, energy is conserved irrespective of the number of particles, the number of neighbours or the choice of smoothing kernel. The only non-conservation of energy arises from the ordinary differential equation solver one employs to solve the left-hand side of the equations. We thus employ a symplectic time integration scheme in order to preserve the conservation properties as accurately as possible (Section 2.3.1).

2.2.7. Shock capturing: momentum equation

The shock capturing dissipation terms are implemented following Monaghan (Reference Monaghan1997), derived by analogy with Riemann solvers from the special relativistic dissipation terms proposed by Chow & Monaghan (Reference Chow and Monaghan1997). These were extended by Price & Monaghan (Reference Price and Monaghan2004b, Reference Price and Monaghan2005) to MHD and recently to dust–gas mixtures by Laibe & Price (Reference Laibe and Price2014b). In a recent paper, Puri & Ramachandran (Reference Puri and Ramachandran2014) found this approach, along with the Morris & Monaghan (Reference Morris and Monaghan1997) switch (which they referred to as the ‘Monaghan–Price–Morris’ formulation) to be the most accurate and robust method for shock capturing in SPH when compared to several other approaches, including Godunov SPH (e.g. Inutsuka Reference Inutsuka2002; Cha & Whitworth Reference Cha and Whitworth2003).

The formulation in Phantom differs from that given in Price (Reference Price2012a) only by the way that the density and signal speed in the q terms are averaged, as described in Price & Federrath (Reference Price and Federrath2010) and Lodato & Price (Reference Lodato and Price2010). That is, we use

(39) $$\begin{equation} \Pi _{\rm shock}^{a} \equiv -\sum _{b} m_{b} \left[ \frac{q^{a}_{ab}}{\rho _{a}^{2}\Omega _{a}}\nabla _a W_{ab}(h_{a}) + \frac{q^{b}_{ab}}{\rho _{b}^{2}\Omega _{b}}\nabla _a W_{ab}(h_{b}) \right], \end{equation}$$

where

(40) $$\begin{equation} q^{a}_{ab} = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}-\displaystyle\frac{1}{2} \rho _{a} v_{{\rm sig}, a} {\bm v}_{ab} \cdot \hat{\bm r}_{ab}, & {\bm v}_{ab} \cdot \hat{\bm r}_{ab} < 0 \\[6pt] 0 & \text{otherwise,} \end{array}\right. \end{equation}$$

where ${\bm v}_{ab} \equiv {\bm v}_{a} - {\bm v}_{b}$, $\hat{\bm r}_{ab} \equiv ({\bm r}_{a} - {\bm r}_{b})/\vert {\bm r}_{a} - {\bm r}_{b} \vert$ is the unit vector along the line of sight between the particles, and v sig is the maximum signal speed, which depends on the physics implemented. For hydrodynamics, this is given by

(41) $$\begin{equation} v_{{\rm sig},a} = \alpha ^{\rm AV}_{a} c_{{\rm s},a} + \beta ^{\rm AV} \vert {\bm v}_{ab} \cdot \hat{\bm r}_{ab} \vert , \end{equation}$$

where in general αAVa ∈ [0, 1] is controlled by a switch (see Section 2.2.9), while βAV = 2 by default.

Importantly, α does not multiply the βAV term. The βAV term provides a second-order Von Neumann & Richtmyer-like term that prevents particle interpenetration (e.g. Lattanzio et al. Reference Lattanzio, Monaghan, Pongracic and Schwarz1986; Monaghan Reference Monaghan1989), and thus βAV ⩾ 2 is needed wherever particle penetration may occur. This is important in accretion disc simulations where use of a low α may be acceptable in the absence of strong shocks, but a low β will lead to particle penetration of the disc midplane, which is the cause of a number of convergence issues (Meru & Bate Reference Meru and Bate2011, Reference Meru and Bate2012). Price & Federrath (Reference Price and Federrath2010) found that βAV = 4 was necessary at high Mach number (M ≳ 5) to maintain a sharp shock structure, which despite nominally increasing the viscosity was found to give less dissipation overall because particle penetration no longer occurred at shock fronts.

2.2.8. Shock capturing: internal energy equation

The key insight from Chow & Monaghan (Reference Chow and Monaghan1997) was that shock capturing not only involves a viscosity term but involves dissipating the jump in each component of the energy, implying a conductivity term in hydrodynamics and resistive dissipation in MHD (see Section 2.10.5). The resulting contribution to the internal energy equation is given by (e.g. Price Reference Price2012a)

(42) $$\begin{eqnarray} \Lambda _{\rm shock} & \equiv& - \frac{1}{\Omega _{a}\rho _a} \sum _{b} m_{b} v_{{\rm sig}, a} \frac{1}{2} ({\bm v}_{ab} \cdot \hat{\bm r}_{ab})^{2} F_{ab}(h_{a}) \nonumber \\ && +\, \sum _{b} m_{b} \alpha _{u} v^{u}_{\rm sig} (u_{a} - u_{b}) \frac{1}{2} \left[ \frac{F_{ab} (h_{a})}{\Omega _{a} \rho _{a}} + \frac{F_{ab} (h_{b})}{\Omega _{b} \rho _{b}} \right] \nonumber \\ && +\, \Lambda _{\rm artres}, \end{eqnarray}$$

where the first term provides the viscous shock heating, the second term provides an artificial thermal conductivity, F ab is defined as in (15), and Λartres is the heating due to artificial resistivity [Equation (182)]. The signal speed we use for conductivity term differs from the one used for viscosity, as discussed by Price (Reference Price2008, Reference Price2012a). In Phantom, we use

(43) $$\begin{equation} v^{u}_{\rm sig} = \sqrt{\frac{\vert P_{a} - P_{b} \vert }{\overline{\rho }_{ab}}} \end{equation}$$

for simulations that do not involve self-gravity or external body forces (Price Reference Price2008), and

(44) $$\begin{equation} v^{u}_{\rm sig} = \vert {\bm v}_{ab} \cdot \hat{\bm r}_{ab} \vert \end{equation}$$

for simulations that do (Wadsley, Veeravalli, & Couchman Reference Wadsley, Veeravalli and Couchman2008). The importance of the conductivity term for treating contact discontinuities was highlighted by Price (Reference Price2008), explaining the poor results found by Agertz et al. (Reference Agertz2007) in SPH simulations of Kelvin–Helmholtz instabilities run across contact discontinuities (discussed further in Section 5.1.4). With (44), we have found there is no need for further switches to reduce conductivity (e.g. Price Reference Price2004; Price & Monaghan Reference Price and Monaghan2005; Valdarnini Reference Valdarnini2016), since the effective thermal conductivity κ is second order in the smoothing length (∝h 2). Phantom therefore uses αu = 1 by default in (42) and we have not yet found a situation where this leads to excess smoothing.

It may be readily shown that the total energy remains conserved in the presence of dissipation by combining (42) with the corresponding dissipative terms in (34). The contribution to the entropy from both viscosity and conductivity is also positive definite (see the appendix in Price & Monaghan Reference Price and Monaghan2004b for the mathematical proof in the case of conductivity).

2.2.9. Shock detection

The standard approach to reducing dissipation in SPH away from shocks for the last 15 yr has been the switch proposed by Morris & Monaghan (Reference Morris and Monaghan1997), where the dimensionless viscosity parameter α is evolved for each particle a according to

(45) $$\begin{equation} \frac{{\rm d}\alpha _{a}}{{\rm d}t} = \max (-\nabla \cdot {\bm v}_{a}, 0) - \frac{(\alpha _{a} - \alpha _{\rm min})}{\tau _{a}}, \end{equation}$$

where τ ≡ h/(σdecayv sig) and σdecay = 0.1 by default. We set v sig in the decay time equal to the sound speed to avoid the need to store dα/dt, since $\nabla \cdot {\bm v}$ is already stored in order to compute (4). This is the switch used for numerous turbulence applications with Phantom (e.g. Price & Federrath Reference Price and Federrath2010; Price et al. Reference Price, Federrath and Brunt2011; Tricco et al. Reference Tricco, Price and Federrath2016b) where it is important to minimise numerical dissipation in order to maximise the Reynolds number (e.g. Valdarnini Reference Valdarnini2011; Price Reference Price2012b).

More recently, Cullen & Dehnen (Reference Cullen and Dehnen2010) proposed a more advanced switch using the time derivative of the velocity divergence. A modified version based on the gradient of the velocity divergence was also proposed by Read & Hayfield (Reference Read and Hayfield2012). We implement a variation on the Cullen & Dehnen (Reference Cullen and Dehnen2010) switch, using a shock indicator of the form

(46) $$\begin{equation} A_{a} = \xi _{a} \max \left[-\frac{{\rm d}}{{\rm d}t}(\nabla \cdot {\bm v}_{a}), 0 \right] , \end{equation}$$

where

(47) $$\begin{equation} \xi = \frac{\vert \nabla \cdot {\bm v} \vert ^{2}}{\vert \nabla \cdot {\bm v} \vert ^{2} + \vert \nabla \times {\bm v} \vert ^{2}} \end{equation}$$

is a modification of the Balsara (Reference Balsara1995) viscosity limiter for shear flows. We use this to set α according to

(48) $$\begin{equation} \alpha _{{\rm loc}, a} =\min \left( \frac{10 h_{a}^{2} A_{a}}{c_{{\rm s}, a}^{2}}, \alpha _{\rm max} \right), \end{equation}$$

where c s is the sound speed and αmax = 1. We use c s in the expression for αloc also for MHD (Section 2.10) since we found using the magnetosonic speed led to a poor treatment of MHD shocks. If αloc, a > αa, we set αa = αloc, a, otherwise we evolve αa according to

(49) $$\begin{equation} \frac{{\rm d}\alpha _{a}}{{\rm d}t} = - \frac{(\alpha _{a} - \alpha _{{\rm loc}, a})}{\tau _{a}}, \end{equation}$$

where τ is defined as in the Morris & Monaghan (Reference Morris and Monaghan1997) version, above. We evolve α in the predictor part of the integrator only, i.e. with a first-order time integration, to avoid complications in the corrector step. However, we perform the predictor step implicitly using a backward Euler method, i.e.

(50) $$\begin{equation} \alpha ^{n+1}_{a} = \frac{\alpha _{a}^{n} + \alpha _{{\rm loc}, a} \Delta t / \tau _{a}}{1 + \Delta t / \tau _{a}}, \end{equation}$$

which ensures that the decay is stable regardless of the timestep (we do this for the Morris & Monaghan method also).

We use the method outlined in Appendix B3 of Cullen & Dehnen (Reference Cullen and Dehnen2010) to compute ${\rm d}(\nabla \cdot {\bm v}_{a})/{{\rm d}t}$. That is, we first compute the gradient tensors of the velocity, ${\bm v}$, and acceleration, ${\bm a}$ (used from the previous timestep), during the density loop using an SPH derivative operator that is exact to linear order, that is, with the matrix correction outlined in Price (Reference Price2004, Reference Price2012a), namely

(51) $$\begin{equation} R^{ij}_{a} \frac{\partial v^{k}_{a}}{\partial x^{j}_{a}} = \sum _{b} m_{b} \left(v^{k}_{b} - v^{k}_{a}\right) \frac{\partial W_{ab}(h_{a})}{\partial x^{i}} , \end{equation}$$

where

(52) $$\begin{equation} R^{ij}_{a} = \sum _{b} m_{b} \left(x^{i}_{b} - x^{i}_{a}\right) \frac{\partial W_{ab}(h_{a})}{\partial x^{j}} \approx \delta ^{ij}, \end{equation}$$

and repeated tensor indices imply summation. Finally, we construct the time derivative of the velocity divergence according to

(53) $$\begin{equation} \frac{{\rm d}}{{\rm d}t}\left(\frac{\partial v_{a}^{i}}{\partial x_{a}^{i}}\right) = \frac{\partial a_{a}^{i}}{\partial x_{a}^{i}} - \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}} \frac{\partial v_{a}^{j}}{\partial x_{a}^{i}}, \end{equation}$$

where, as previously, repeated i and j indices imply summation. In Cartesian coordinates, the last term in (53) can be written out explicitly using

(54) $$\begin{eqnarray} \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}} \frac{\partial v_{a}^{j}}{\partial x_{a}^{i}} & =& \left(\frac{\partial v^{x}}{\partial x} \right)^{2} + \left(\frac{\partial v^{y}}{\partial y} \right)^{2} + \left(\frac{\partial v^{z}}{\partial z} \right)^{2} \nonumber \\ & &+\, 2 \left[\frac{\partial v^{x}}{\partial y} \frac{\partial v^{y}}{\partial x} + \frac{\partial v^{x}}{\partial z} \frac{\partial v^{z}}{\partial x} + \frac{\partial v^{z}}{\partial y} \frac{\partial v^{y}}{\partial z}\right]. \end{eqnarray}$$

2.2.10. Cooling

The cooling term Λcool can be set either from detailed chemical calculations (Section 2.14.1) or, for discs, by the simple ‘β-cooling’ prescription of Gammie (Reference Gammie2001), namely

(55) $$\begin{equation} \Lambda _{\rm cool} = \frac{\rho u}{t_{\rm cool}}, \end{equation}$$

where

(56) $$\begin{equation} t_{\rm cool} \equiv \frac{\Omega (R)}{\beta _{\rm cool}}, \end{equation}$$

with βcool an input parameter to the code specifying the cooling timescale in terms of the local orbital time. We compute Ω in (56) using Ω ≡ 1/(x 2 + y 2 + z 2)3/2, i.e. assuming Keplerian rotation around a central object with mass equal to unity, with G = 1 in code units.

2.2.11. Conservation of linear and angular momentum

The total linear momentum is given by

(57) $$\begin{equation} \bm {P} = \sum _{a} m_{a} {\bm v}_{a}, \end{equation}$$

such that conservation of momentum corresponds to

(58) $$\begin{equation} \frac{{\rm d}\bm {P}}{{\rm d}t} = \sum _{a} m_{a} \frac{{\rm d}{\bm v}_{a}}{{\rm d}t} = 0. \end{equation}$$

Inserting our discrete equation (34), we find

(59) $$\begin{eqnarray} \frac{{\rm d}\bm {P}}{{\rm d}t} &=& \sum _{a} \sum _{b} m_a m_{b} \left[ \frac{P_{a} + q^{a}_{ab}}{\rho _{a}^{2}\Omega _{a}}\nabla _a W_{ab}(h_{a}) \right. \nonumber \\ && \left. +\, \frac{P_{b}+q^{b}_{ab}}{\rho _{b}^{2}\Omega _{b}}\nabla _a W_{ab}(h_{b}) \right] = 0, \end{eqnarray}$$

where, as for the total energy (Section 2.2.6), the double summation is zero because of the antisymmetry of the kernel gradient. The same argument applies to the conservation of angular momentum,

(60) $$\begin{equation} \sum _{a} m_{a} {\bm r}_{a} \times {\bm v}_{a} \end{equation}$$

(see e.g. Price Reference Price2012a for a detailed proof). As with total energy, this means linear and angular momentum are exactly conserved by our SPH scheme to the accuracy with which they are conserved by the timestepping scheme.

In Phantom, linear and angular momentum are both conserved to round-off error (typically ~10−16 in double precision) with global timestepping, but exact conservation is violated when using individual particle timesteps or when using the kd-tree to compute gravitational forces. The magnitude of these quantities, as well as the total energy and the individual components of energy (kinetic, internal, potential, and magnetic), should thus be monitored by the user at runtime. Typically with individual timesteps, one should expect energy conservation to ΔE/E ~ 10−3 and linear and angular momentum conservation to ~10−6 with default code settings. The code execution is aborted if conservation errors exceed 10%.

2.3. Time integration

2.3.1. Timestepping algorithm

We integrate the equations of motion using a generalisation of the Leapfrog integrator which is reversible in the case of both velocity dependent forces and derivatives which depend on the velocity field. The basic integrator is the Leapfrog method in ‘Kick–Drift–Kick’ or ‘Velocity Verlet’ form (Verlet Reference Verlet1967), where the positions and velocities of particles are updated from time t n to t n + 1 according to

(61) $$\begin{eqnarray} {\bm v}^{n+\frac{1}{2}} = {\bm v}^{n} + \frac{1}{2} \Delta t {\bm a}^{n}, \end{eqnarray}$$
(62) $$\begin{eqnarray} {\bm r}^{n+1} = {\bm r}^{n} + \Delta t {\bm v}^{n+\frac{1}{2}}, \end{eqnarray}$$
(63) $$\begin{eqnarray} {\bm a}^{n + 1} = {\bm a}({\bm r}^{n+1}), \end{eqnarray}$$
(64) $$\begin{eqnarray} {\bm v}^{n + 1} = {\bm v}^{n+\frac{1}{2}} + \frac{1}{2} \Delta t {\bm a}^{n+1}, \end{eqnarray}$$

where Δtt n + 1t n. This is identical to the formulation of Leapfrog used in other astrophysical SPH codes (e.g. Springel Reference Springel2005; Wadsley et al. Reference Wadsley, Stadel and Quinn2004). The Verlet scheme, being both reversible and symplectic (e.g. Hairer, Lubich, & Wanner Reference Hairer, Lubich and Wanner2003), preserves the Hamiltonian nature of the SPH algorithm (e.g. Gingold & Monaghan Reference Gingold and Monaghan1982b; Monaghan & Price Reference Monaghan and Price2001). In particular, both linear and angular momentum are exactly conserved, there is no long-term energy drift, and phase space volume is conserved (e.g. for orbital dynamics). In SPH, this is complicated by velocity-dependent terms in the acceleration from the shock-capturing dissipation terms. In this case, the corrector step, (64), becomes implicit. The approach we take is to notice that these terms are not usually dominant over the position-dependent terms. Hence, we use a first-order prediction of the velocity, as follows:

(65) $$\begin{eqnarray} {\bm v}^{n+\frac{1}{2}} = {\bm v}^{n} + \frac{1}{2} \Delta t {\bm a}^{n}, \end{eqnarray}$$
(66) $$\begin{eqnarray} {\bm r}^{n+1} = {\bm r}^{n} + \Delta t {\bm v}^{n+\frac{1}{2}}, \end{eqnarray}$$
(67) $$\begin{eqnarray} {\bm v}^{*} = {\bm v}^{n+\frac{1}{2}} + \frac{1}{2} \Delta t {\bm a}^{n}, \end{eqnarray}$$
(68) $$\begin{eqnarray} {\bm a}^{n + 1} = {\bm a}({\bm r}^{n+1}, {\bm v}^{*}), \end{eqnarray}$$
(69) $$\begin{eqnarray} {\bm v}^{n + 1} = {\bm v}^{*} + \frac{1}{2} \Delta t \left[ {\bm a}^{n+1} - {\bm a}^{n} \right]. \end{eqnarray}$$

At the end of the step, we then check if the error in the first-order prediction is less than some tolerance ε according to

(70) $$\begin{equation} e = \frac{ \vert {\bm v}^{n+1} - {\bm v}^{*} \vert }{ \vert {\bm v}^{\rm mag} \vert } < \epsilon _{\rm v}, \end{equation}$$

where ${\bm v}^{\rm mag}$ is the mean velocity on all SPH particles (we set the error to zero if $\vert {\bm v}^{\rm mag}\vert = 0$) and by default εv = 10−2. If this criterion is violated, then we recompute the accelerations by replacing ${\bm v}^{*}$ with ${\bm v}^{n+1}$ and iterating (68) and (69) until the criterion in (70) is satisfied. In practice, this happens rarely, but occurs for example in the first few steps of the Sedov problem where the initial conditions are discontinuous (Section 5.1.3). As each iteration is as expensive as halving the timestep, we also constrain the subsequent timestep such that iterations should not occur, i.e.

(71) $$\begin{equation} \Delta t = \min \left( \Delta t, \frac{\Delta t}{\sqrt{e_{\rm max} / \epsilon }} \right), \end{equation}$$

where e max = max(e) over all particles. A caveat to the above is that velocity iterations are not currently implemented when using individual particle timesteps.

Additional variables such as the internal energy, u, and the magnetic field, ${\bm B}$, are timestepped with a predictor and trapezoidal corrector step in the same manner as the velocity, following (65), (67) and (69).

Velocity-dependent external forces are treated separately, as described in Section 2.4.

2.3.2. Timestep constraints

The timestep itself is determined at the end of each step, and is constrained to be less than the maximum stable timestep. For a given particle, a, this is given by (e.g. Lattanzio et al. Reference Lattanzio, Monaghan, Pongracic and Schwarz1986; Monaghan Reference Monaghan1997)

(72) $$\begin{equation} \Delta t_{{\rm C},a} \equiv C_{\rm cour} \frac{h_a}{v^{\rm dt}_{{\rm sig},a}}, \end{equation}$$

where C cour = 0.3 by default (Lattanzio et al. Reference Lattanzio, Monaghan, Pongracic and Schwarz1986) and v dtsig is taken as the maximum of (41) over the particle’s neighbours assuming αAV = max(αAV, 1). The criterion above differs from the usual Courant–Friedrichs–Lewy condition used in Eulerian codes (Courant, Friedrichs, & Lewy Reference Courant, Friedrichs and Lewy1928) because it depends only on the difference in velocity between neighbouring particles, not the absolute value.

An additional constraint is applied from the accelerations (the ‘force condition’), where

(73) $$\begin{equation} \Delta t_{{\rm f}, a} \equiv C_{\rm force} \sqrt{\frac{h_{a}}{\vert {\bm a}_{a} \vert }}, \end{equation}$$

where C force = 0.25 by default. A separate timestep constraint is applied for external forces

(74) $$\begin{equation} \Delta t_{{\rm ext},a} \equiv C_{\rm force} \sqrt{\frac{h}{\vert {\bm a}_{{\rm ext},a} \vert }}, \end{equation}$$

and for accelerations to SPH particles to/from sink particles (Section 2.8)

(75) $$\begin{equation} \Delta t_{{\rm sink-gas},a} \equiv C_{\rm force} \sqrt{\frac{h_{a}}{\vert {\bm a}_{{\rm sink-gas},a} \vert }}. \end{equation}$$

For external forces with potentials defined such that Φ → 0 as r → ∞, an additional constraint is applied using (Dehnen & Read Reference Dehnen and Read2011)

(76) $$\begin{equation} \Delta t_{{\Phi },a} \equiv C_{\rm force} \eta _\Phi \sqrt{\frac{\vert \Phi _{a}\vert }{\vert \nabla \Phi \vert _{a}^{2}}}, \end{equation}$$

where $\eta _\Phi = 0.05$ (see Section 2.8.5).

The timestep for particle a is then taken to be the minimum of all of the above constraints, i.e.

(77) $$\begin{eqnarray} \Delta t_{a} = \min & \left(\Delta t_{\rm C}, \Delta t_{\rm f}, \Delta t_{\rm ext}, \Delta t_{\rm sink-gas}, \Delta t_{\Phi } \right)_a, \end{eqnarray}$$

with possible other constraints arising from additional physics as described in their respective sections. With global timestepping, the resulting timestep is the minimum over all particles

(78) $$\begin{equation} \Delta t = \min _{a} ( \Delta t_a). \end{equation}$$

2.3.3. Substepping of external forces

In the case where the timestep is dominated by any of the external force timesteps, i.e. (74)–(76), we implement an operator splitting approach implemented according to the reversible reference system propagator algorithm (RESPA) derived by Tuckerman, Berne, & Martyna (Reference Tuckerman, Berne and Martyna1992) for molecular dynamics. RESPA splits the acceleration into ‘long-range’ and ‘short-range’ contributions, which in Phantom are defined to be the SPH and external/point-mass accelerations, respectively.

Our implementation follows Tuckerman et al. (Reference Tuckerman, Berne and Martyna1992) (see their Appendix B), where the velocity is first predicted to the half step using the ‘long-range’ forces, followed by an inner loop where the positions are updated with the current velocity and the velocities are updated with the ‘short-range’ accelerations. Thus, the timestepping proceeds according to

(79) $$\begin{equation} {\bm v} = {\bm v} + \frac{\Delta t_{\rm sph}}{2} {\bm a}_{\rm sph}^{n},\vspace*{-6pt} \end{equation}$$
\begin{eqnarray*} {\rm over\ substeps}\left\{\begin{array}{rcl} \bm v &=& {\bm v} + \displaystyle\frac{\Delta t_{\rm ext}}{2} {\bm a}_{\rm ext}^{m}, \\ {\bm r} &=& {\bm r} + \Delta t_{\rm ext} {\bm v}, \\ &&\qquad\textrm {get } {\bm a}_{\rm ext}({\bm r}),\\ {\bm v} &=& {\bm v} + \frac{\Delta t_{\rm ext}}{2} {\bm a}_{\rm ext}^{m+1}, \end{array}\right. \end{eqnarray*}
(84) $$\begin{equation} \hspace*{-15.53pt}\textrm {get } {\bm a}_{\rm sph}({\bm r}),\hspace*{15.53pt} \vspace*{-6pt} \end{equation}$$
(85) $$\begin{equation} {\bm v} = {\bm v} + \frac{\Delta t_{\rm sph}}{2} {\bm a}_{\rm sph}^{n}, \end{equation}$$

where ${\bm a}_{\rm sph}$ indicates the SPH acceleration evaluated from (34) and ${\bm a}_{\rm ext}$ indicates the external forces. The SPH and external accelerations are stored separately to enable this. Δt ext is the minimum of all timesteps relating to sink–gas and external forces [equations (74)–(76)], while Δt sph is the timestep relating to the SPH forces [equations (72), (73), and (288)]. Δt ext is allowed to vary on each substep, so we take as many steps as required such that ∑m − 1jΔt ext, j + Δt ext, f = Δt sph, where Δt ext, f < Δt ext, j is chosen so that the sum will identically equal Δt sph. The number of substeps is m ≈ int(Δt ext, mint sph, min) + 1, where the minimum is taken over all particles.

2.3.4. Individual particle timesteps

For simulations of stiff problems with a large range in timestep over the domain, it is more efficient to allow each particle to evolve on its own timestep independently (Bate Reference Bate1995; Springel Reference Springel2005; Saitoh & Makino Reference Saitoh and Makino2010). This violates all of the conservation properties of the Leapfrog integrator [see Makino et al. (Reference Makino, Hut, Kaplan and Saygın2006) for an attempt to solve this], but can speed up the calculation by an order of magnitude or more. We implement this in the usual blockstepped manner by assigning timesteps in factor-of-two decrements from some maximum timestep Δt max, which for convenience is set equal to the time between output files.

We implement a timestep limiter where the timestep for an active particle is constrained to be within a factor of 2 of its neighbours, similar to condition employed by Saitoh & Makino (Reference Saitoh and Makino2009). Additionally, inactive particles will be woken up as required to ensure that their timestep is within a factor of 2 of its neighbours.

The practical side of individual timestepping is described in Appendix A.6.

2.4. External forces

2.4.1. Point-mass potential

The simplest external force describes a point mass, M, at the origin, which yields gravitational potential and acceleration:

(86) $$\begin{equation} \Phi _{a} = -\frac{GM}{r_{a}}, \hspace{28.45274pt} {\bm a}_{{\rm ext},a} = -\nabla \Phi _{a} = -\frac{GM}{\vert {\bm r}_{a} \vert ^{3}} {\bm r}_{a}, \end{equation}$$

where $r_{a} \equiv \vert {\bm r}_{a} \vert \equiv \sqrt{{\bm r}_{a}\cdot {\bm r}_{a}}$. When this potential is used, we allow for particles within a certain radius, R acc, from the origin to be accreted. This allows for a simple treatment of accretion discs where the mass of the disc is assumed to be negligible compared to the mass of the central object. The accreted mass in this case is recorded but not added to the central mass. For more massive discs, or when the accreted mass is significant with respect to the central mass, it is better to model the central star using a sink particle (Section 2.8) where there are mutual gravitational forces between the star and the disc, and any mass accreted is added to the point mass (Section 2.8.2).

2.4.2. Binary potential

We provide the option to model motion in binary systems where the mass of the disc is negligible. In this case, the binary motion is prescribed using

(87) $$\begin{eqnarray} {\bm r}_{1} = [(1 - M) \cos (t),(1- M) \sin (t), 0], \end{eqnarray}$$
(88) $$\begin{eqnarray} {\bm r}_{2} = [- M \cos (t), - M \sin (t), 0], \end{eqnarray}$$

where M is the mass ratio in units of the total mass (which is therefore unity). For this potential, G and Ω are set to unity in computational units, where Ω is the angular velocity of the binary. Thus, only M needs to be specified to fix both m 1 and m 2. Hence, the binary remains fixed on a circular orbit at r = 1. The binary potential is therefore

(89) $$\begin{equation} \Phi _{a} = -\frac{M}{\vert {\bm r}_{a} - {\bm r}_{1} \vert } - \frac{(1- M) }{\vert {\bm r}_{a} - {\bm r}_{2} \vert }, \end{equation}$$

such that the external acceleration is given by

(90) $$\begin{equation} {\bm a}_{{\rm ext},a} = -\nabla \Phi _{a} = - M \frac{({\bm r}_{a} - {\bm r}_{1})}{\vert {\bm r}_{a} - {\bm r}_{1} \vert ^3} - (1- M) \frac{({\bm r}_{a} - {\bm r}_{2})}{\vert {\bm r}_{a} - {\bm r}_{2} \vert ^3}. \end{equation}$$

Again, there is an option to accrete particles that fall within a certain radius from either star (R acc, 1 or R acc, 2, respectively). For most binary accretion disc simulations (e.g. planet migration), it is better to use ‘live’ sink particles to represent the binary so that there is feedback between the binary and the disc (we have used a live binary in all of our simulations to date, e.g. Nixon, King, & Price Reference Nixon, King and Price2013; Facchini, Lodato, & Price Reference Facchini, Lodato and Price2013; Martin et al. Reference Martin, Nixon, Armitage, Lubow and Price2014a, Reference Martin, Nixon, Lubow, Armitage, Price, Doğan and King2014b; Doğan et al. Reference Doğan, Nixon, King and Price2015; Ragusa, Lodato, & Price Reference Ragusa, Lodato and Price2016; Ragusa et al. Reference Ragusa, Dipierro, Lodato, Laibe and Price2017), but the binary potential remains useful under limited circumstances—in particular, when one wishes to turn off the feedback between the disc and the binary.

Given that the binary potential is time-dependent, for efficiency, we compute the position of the binary only once at the start of each timestep, and use these stored positions to compute the accelerations of the SPH particles via (90).

2.4.3. Binary potential with gravitational wave decay

An alternative binary potential including the effects of gravitational wave decay was used by Cerioli, Lodato, & Price (Reference Cerioli, Lodato and Price2016) to study the squeezing of discs during the merger of supermassive black holes. Here the motion of the binary is prescribed according to

(91) $$\begin{eqnarray} {\bm r}_{1} & =& \left[-\frac{m_{2}}{m_{1} + m_{2}} a \cos (\theta ), -\frac{m_{2}}{m_{1} + m_{2}} a \sin (\theta ), 0 \right], \nonumber \\ {\bm r}_{2} & =& \left[\frac{m_{1}}{m_{1} + m_{2}} a \cos (\theta ), \frac{m_{1}}{m_{1} + m_{2}} a \sin (\theta ), 0 \right], \end{eqnarray}$$

where the semi-major axis, a, decays according to

(92) $$\begin{equation} a(t) = a_{0} \left(1 - \frac{t}{\tau } \right)^{\frac{1}{4}} . \end{equation}$$

The initial separation is a 0, with τ defined as the time to merger, given by the usual expression (e.g. Lodato et al. Reference Lodato, Nayakshin, King and Pringle2009)

(93) $$\begin{equation} \tau \equiv \frac{5}{256} \frac{a_{0}^{4}}{\mu _{12} (m_{1} + m_{2})^{2}}, \end{equation}$$

where

(94) $$\begin{equation} \mu _{12} \equiv \frac{m_1 m_2}{m_1 + m_2}. \end{equation}$$

The angle θ is defined using

(95) $$\begin{equation} \Omega \equiv \frac{{\rm d}\theta }{{\rm d}t} = \sqrt{\frac{G(m_{1} + m_{2})}{a^{3}}}. \end{equation}$$

Inserting the expression for a and integrating gives (Cerioli et al. Reference Cerioli, Lodato and Price2016)

(96) $$\begin{equation} \theta (t) = -\frac{8\tau }{5} \sqrt{\frac{G(m_{1} + m_{2})}{a_{0}^{3}}} \left( 1 - \frac{t}{\tau } \right). \end{equation}$$

The positions of the binary, ${\bm r}_1$ and ${\bm r}_2$, can be inserted into (89) to obtain the binary potential, with the acceleration as given in (90). The above can be used as a simple example of a time-dependent external potential.

2.4.4. Galactic potentials

We implement a range of external forces representing various galactic potentials, as used in Pettitt et al. (Reference Pettitt, Dobbs, Acreman and Price2014). These include arm, bar, halo, disc, and spheroidal components. We refer the reader to the paper above for the actual forms of the potentials.

For the non-axisymmetric potentials, a few important parameters that determine the morphology can be changed at run time rather than compile time. These include the pattern speed, arm number, arm pitch angle, and bar axis lengths (where applicable). In the case of non-axisymmetric components, the user should be aware that some will add mass to the system, whereas others simply perturb the galactic disc. These potentials can be used for any galactic system, but the various default scale lengths and masses are chosen to match the Milky Way’s rotation curve (Sofue Reference Sofue2012).

The most basic potential in Phantom is a simple logarithmic potential from Binney & Tremaine (Reference Binney and Tremaine1987), which allows for the reproduction of a purely flat rotation curve with steep decrease at the galactic centre, and approximates the halo, bulge, and disc contributions. Also included is the standard flattened disc potential of Miyamoto–Nagai (Miyamoto & Nagai Reference Miyamoto and Nagai1975) and an exponential profile disc, specifically the form from Khoperskov et al. (Reference Khoperskov, Vasiliev, Sobolev and Khoperskov2013). Several spheroidal components are available, including the potentials of Plummer (Reference Plummer1911), Hernquist (Reference Hernquist1990), and Hubble (Reference Hubble1930). These can be used generally for bulges and halos if given suitable mass and scale lengths. We also include a few halo-specific profiles; the NFW (Navarro, Frenk, & White Reference Navarro, Frenk and White1996), Begeman, Broeils, & Sanders (Reference Begeman, Broeils and Sanders1991), Caldwell & Ostriker (Reference Caldwell and Ostriker1981), and the Allen & Santillan (Reference Allen and Santillan1991) potentials.

The arm potentials include some of the more complicated profiles. The first is the potential of Cox & Gómez (Reference Cox and Gómez2002), which is a relatively straightforward superposition of three sinusoidal-based spiral components to damp the potential ‘troughs’ in the inter-arm minima. The other spiral potential is from Pichardo et al. (Reference Pichardo, Martos, Moreno and Espresate2003), and is more complicated. Here, the arms are constructed from a superposition of oblate spheroids whose loci are placed along a standard logarithmic spiral. As the force from this potential is computationally expensive it is prudent to pre-compute a grid of potential/force and read it at run time. The python code to generate the appropriate grid files is distributed with the code.

Finally, the bar components: We include the bar potentials of Dehnen (Reference Dehnen2000a), Wada & Koda (Reference Wada and Koda2001), the ‘S’ shaped bar of Vogt & Letelier (Reference Vogt and Letelier2011), both biaxial and triaxial versions provided in Long & Murali (Reference Long and Murali1992), and the boxy-bulge bar of Wang et al. (Reference Wang, Zhao, Mao and Rich2012). This final bar contains both a small inner non-axisymmetric bulge and longer bar component, with the forces calculated by use of Hernquist–Ostriker expansion coefficients of the bar density field. Phantom contains the coefficients for several different forms of this bar potential.

2.4.5. Lense–Thirring precession

Lense–Thirring precession (Lense & Thirring Reference Lense and Thirring1918) from a spinning black hole is implemented in a post-Newtonian approximation following Nelson & Papaloizou (Reference Nelson and Papaloizou2000), which has been used in Nixon et al. (Reference Nixon, King, Price and Frank2012b), Nealon, Price, & Nixon (Reference Nealon, Price and Nixon2015), and Nealon et al. (Reference Nealon, Nixon, Price and King2016). In this case, the external acceleration consists of a point-mass potential (Section 2.4.1) and the Lense–Thirring term:

(97) $$\begin{equation} {\bm a}_{{\rm ext},a} = -\nabla \Phi _a + {\bm v}_a \times {\bm \Omega }_{p,a}, \end{equation}$$

where Φa is given by (86) and ${\bm v}_a \times {\bm \Omega }_{p,a}$ is the gravitomagnetic acceleration. A dipole approximation is used, yielding

(98) $$\begin{equation} {\bm \Omega }_{p,a} \equiv \frac{2{\bm S}}{\vert {\bm r}_a \vert ^{3}} - \frac{6({\bm S}\cdot {\bm r}_a) {\bm r}_a}{\vert {\bm r}_a \vert ^{5}}, \end{equation}$$

with ${\bm S} = a_{\rm spin} (GM)^{2} {\bm k}/c^{3}$, where ${\bm k}$ is a unit vector in the direction of the black hole spin. When using the Lense–Thirring force, geometric units are assumed such that G = M = c = 1, as described in Section 2.2.3, but with the additional constraints on the unit system from M and c.

Since in this case the external force depends on velocity, it cannot be implemented directly into Leapfrog. The method we employ to achieve this is simpler than those proposed elsewhere [c.f. attempts by Quinn et al. (Reference Quinn, Perrine, Richardson and Barnes2010) and Rein & Tremaine (Reference Rein and Tremaine2011) to adapt the Leapfrog integrator to Hill’s equations]. Our approach is to split the acceleration into position and velocity-dependent parts, i.e.

(99) $$\begin{equation} {\bm a}_{\rm ext} = {\bm a}_{\rm ext,x}({\bm r}) + {\bm a}_{\rm ext,v}({\bm r}, {\bm v}). \end{equation}$$

The position-dependent part [i.e. $-\nabla \Phi ({\bm r})$] is integrated as normal. The velocity-dependent Lense–Thirring term is added to the predictor step, (66)–(67), as usual, but the corrector step, (69), is written in the form

(100) $$\begin{equation} {\bm v}^{n + 1} = {\bm v}^{n + \frac{1}{2}} + \frac{1}{2} \Delta t \left[ {\bm a}_{\rm sph}^{n+1} + {\bm a}_{\rm ext,x}^{n+1} + {\bm a}_{\rm ext,v}({\bm r}^{n+1}, {\bm v}^{n+1})\right], \end{equation}$$

where ${\bm v}^{n + \frac{1}{2}} \equiv {\bm v}^{n} + \frac{1}{2} \Delta t {\bm a}^{n}$ as in (65). This equation is implicit but the trick is to notice that it can be solved analytically for simple forcesFootnote 4. In the case of Lense–Thirring precession, we have

(101) $$\begin{equation} {\bm v}^{n + 1} = \tilde{\bm v} + \frac{1}{2} \Delta t \left[{\bm v}^{n+1} \times {\bm \Omega }_{p}({\bm r}^{n+1}) \right], \end{equation}$$

where $\tilde{\bm v} \equiv {\bm v}^{n + \frac{1}{2}} + \frac{1}{2} \Delta t ({\bm a}_{\rm sph}^{n+1} + {\bm a}_{\rm ext,x}^{n+1})$. We therefore have a matrix equation in the form

(102) $$\begin{equation} {\bm R} {\bm v}^{n+1} = \tilde{\bm v}, \end{equation}$$

where ${\bm R}$ is the 3 × 3 matrix given by

(103) $$\begin{equation} {\bm R} \equiv \left[ \arraycolsep5pt\begin{array}{ccc}1 & - \displaystyle\frac{\Delta t}{2} \Omega ^{z}_{p} & \displaystyle\frac{\Delta t}{2} \Omega ^{y}_{p} \\[3pt] \displaystyle\frac{\Delta t}{2} \Omega ^{z}_{p} & 1 & -\displaystyle\frac{\Delta t}{2} \Omega ^{x}_{p} \\[3pt] -\displaystyle\frac{\Delta t}{2} \Omega ^{y}_{p} & \displaystyle\frac{\Delta t}{2} \Omega ^{x}_{p} & 1 \\ \end{array} \right]. \end{equation}$$

Rearranging (102), ${\bm v}^{n+1}$ is obtained by using

(104) $$\begin{equation} {\bm v}^{n+1} = {\bm R}^{-1}\tilde{\bm v}, \end{equation}$$

where ${\bm R}^{-1}$ is the inverse of ${\bm R}$, which we invert using the analytic solution.

2.4.6. Generalised Newtonian potential

The generalised Newtonian potential described by Tejeda & Rosswog (Reference Tejeda and Rosswog2013) is implemented, where the acceleration terms are given by

(105) $$\begin{equation} {\bm a}_{{\rm ext},a} = -\frac{GM{\bm r}_a}{\vert {\bm r}_a \vert ^{3}}f^{2} + \frac{2R_{\rm g} {\bm v}_a ({\bm v}_a \cdot {\bm r}_a)}{\vert {\bm r}_a \vert ^{3} f} - \frac{3 R_{\rm g} {\bm r}_a ({\bm v}_a \times {\bm r}_a)^{2}}{\vert {\bm r}_a \vert ^{5}}, \end{equation}$$

with R gGM/c 3 and $f \equiv \left(1 - 2R_{\rm g}/\vert {\bm r}_a \vert \right)$. See Bonnerot et al. (Reference Bonnerot, Rossi, Lodato and Price2016) for a recent application. This potential reproduces several features of the Schwarzschild (Reference Schwarzschild1916) spacetime, in particular, reproducing the orbital and epicyclic frequencies to better than 7% (Tejeda & Rosswog Reference Tejeda and Rosswog2013). As the acceleration involves velocity-dependent terms, it requires a semi-implicit solution like Lense–Thirring precession. Since the matrix equation is rather involved for this case, the corrector step is iterated using fixed point iterations until the velocity of each particle is converged to a tolerance of 1%.

2.4.7. Poynting–Robertson drag

The radiation drag from a central point-like, gravitating, radiating, and non-rotating object may be applied as an external force. The implementation is intended to be as general as possible. The acceleration of a particle subject to these external forces is

(106) $$\begin{eqnarray} {\bm a}_{{\rm ext},a} &=& \frac{(k_0\beta _{\rm PR}-1)GM}{\vert {\bm r}_a \vert ^3} \bm r_a \nonumber \\ && -\beta _{\rm PR} \left(k_1\frac{GM}{\vert {\bm r}_a \vert ^3}\frac{v_r}{c}\bm r_a -k_2\frac{GM}{\vert {\bm r}_a \vert ^2}\frac{\bm v_a}{c}\right), \end{eqnarray}$$

where v r is the component of the velocity in the radial direction. The parameter βPR is the ratio of radiation to gravitational forces, supplied by a separate user-written module. Relativistic effects are neglected because these are thought to be less important than radiation forces for low (βPR < 0.01) luminosities, even in accreting neutron star systems where a strong gravitational field is present (e.g., Miller & Lamb Reference Miller and Lamb1993).

The three terms on the right side of (106) correspond, respectively, to gravity (reduced by outward radiation pressure), redshift-related modification to radiation pressure caused by radial motion, and Poynting–Robertson drag against the direction of motion. These three terms can be scaled independently by changing the three parameters k 0, k 1, and k 2, whose default values are unity. Rotation of the central object can be crudely emulated by changing k 2.

As for Lense–Thirring precession, the ${\bm a}^{n+1}$ term of the Leapfrog integration scheme can be expanded into velocity-dependent and non-velocity-dependent component. We obtain, after some algebra,

(107) $$\begin{equation} \bm v^{n+1} =-\frac{\bm T-Qk_1(\bm v^{n+1}\cdot \hat{\bm r})\hat{\bm r}}{1+Qk_2}, \end{equation}$$

where

(108) $$\begin{equation} \bm T=\bm v^n+\frac{1}{2}\Delta t\bm a^n-\frac{(1-k_0\beta _{\rm PR})GM\Delta t}{2r^3}\bm r \end{equation}$$

and

(109) $$\begin{equation} Q=\frac{GM\beta _{\rm PR}\Delta t}{2cr^2}. \end{equation}$$

Equation (107) yields a set of simultaneous equations for the three vector components that can be solved analytically. A detailed derivation is given in Worpel (Reference Worpel2015).

2.4.8. Coriolis and centrifugal forces

Under certain circumstances, it is useful to perform calculations in a co-rotating reference frame (e.g. for damping binary stars into equilibrium with each other). The resulting acceleration terms are given by

(110) $$\begin{equation} {\bm a}_{{\rm ext},a} = -{\bm \Omega } \times ({\bm \Omega } \times {\bm r}_a) - 2 ({\bm \Omega } \times {\bm v}_a), \end{equation}$$

which are the centrifugal and Coriolis terms, respectively, with ${\bm \Omega }$ the angular rotation vector. The timestepping algorithm is as described above for Lense–Thirring precession, with the velocity-dependent term handled by solving the 3 × 3 matrix in the Leapfrog corrector step.

2.5. Driven turbulence

Phantom implements turbulence driving in periodic domains via an Ornstein–Uhlenbeck stochastic driving of the acceleration field, as first suggested by Eswaran & Pope (Reference Eswaran and Pope1988). This is an SPH adaptation of the module used in the grid-based simulations by Schmidt, Hillebrandt, & Niemeyer (Reference Schmidt, Hillebrandt and Niemeyer2006) and Federrath, Klessen, & Schmidt (Reference Federrath, Klessen and Schmidt2008) and many subsequent works. This module was first used in Phantom by Price & Federrath (Reference Price and Federrath2010) to compare the statistics of isothermal, supersonic turbulence between SPH, and grid methods. Subsequent applications have been to the density variance–Mach number relation (Price et al. Reference Price, Federrath and Brunt2011), subsonic turbulence (Price Reference Price2012b), supersonic MHD turbulence (Tricco et al. Reference Tricco, Price and Federrath2016b), and supersonic turbulence in a dust–gas mixture (Tricco et al. Reference Tricco, Price and Laibe2017). Adaptations of this module have also been incorporated into other SPH codes (Bauer & Springel Reference Bauer and Springel2012; Valdarnini Reference Valdarnini2016).

The amplitude and phase of each Fourier mode is initialised by creating a set of six random numbers, z n, drawn from a random Gaussian distribution with unit variance. These are generated by the usual Box–Muller transformation (e.g. Press et al. Reference Press, Teukolsky, Vetterling and Flannery1992) by selecting two uniform random deviates u 1, u 2 ∈ [0, 1] and constructing the amplitude according to

(111) $$\begin{equation} z = \sqrt{2 \log (1 /u_{1})} \cos (2 \pi u_{2}). \end{equation}$$

The six Gaussian random numbers are set up according to

(112) $$\begin{equation} {\bm x}_{n} = \sigma {\bm z}_{n}, \end{equation}$$

where the standard deviation, σ, is set to the square root of the energy per mode divided by the correlation time, $\sigma = \sqrt{E_{\rm m}/t_{\rm decay}}$, where both E m and t decay are user-specified parameters.

The ‘red noise’ sequence (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930) is generated for each mode at each timestep according to (Bartosch Reference Bartosch2001)

(113) $$\begin{equation} {\bm x}_{n+1}= f {\bm x}_{n} + \sigma \sqrt{(1 - f^{2})} {\bm z}_n, \end{equation}$$

where f = exp(− Δt/t decay) is the damping factor. The resulting sequence has zero mean with root-mean-square equal to the variance. The power spectrum in the time domain can vary from white noise [P(f) constant] to ‘brown noise’ [P(f) = 1/f 2].

The amplitudes and phases of each mode are constructed by splitting ${\bm x}_{n}$ into two vectors, ${\bm \Phi }_{a}$ and ${\bm \Phi }_{b}$ of length 3, employed to construct divergence- and curl-free fields according to

(114) $$\begin{eqnarray} {\bm A}_{m} = w [{\bm \Phi }_{a} - ( {\bm \Phi }_{a} \cdot \hat{\bm k} )\hat{\bm k}] + (1 - w)[({\bm \Phi }_{b} \cdot \hat{\bm k} )\hat{\bm k}] , \end{eqnarray}$$\vspace*{-15pt}
(115) $$\begin{eqnarray} {\bm B}_{m} = w [{\bm \Phi }_{b} - ( {\bm \Phi }_{b} \cdot \hat{\bm k} )\hat{\bm k}] + (1 - w)[({\bm \Phi }_{a} \cdot \hat{\bm k} )\hat{\bm k}] , \end{eqnarray}$$

where ${\bm k} = [k_{x}, k_{y}, k_{z}]$ is the mode frequency. The parameter w ∈ [0, 1] is the ‘solenoidal weight’, specifying whether the driving should be purely solenoidal (w = 1) or purely compressive (w = 0) (see Federrath et al. Reference Federrath, Klessen and Schmidt2008, Reference Federrath, Roman-Duval, Klessen, Schmidt and Mac Low2010a).

The spectral form of the driving is defined in Fourier space, with options for three possible spectral forms

(116) $$\begin{equation} C_{m} = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}1 & \textrm {uniform} \\[6pt] 4 (a_{\rm min} - 1) \displaystyle\frac{(k - k_{c})^{2}}{(k_{\max } - k_{\min })^{2}} + 1 & \textrm {parabolic} \\[3pt] k/k_{\min }^{-5/3} & \textrm {Kolmogorov,} \end{array}\right. \end{equation}$$

where $k = \sqrt{k_{x}^{2} + k_{y}^{2} + k_{z}^{2}}$ is the wavenumber, with non-zero amplitudes defined only for wavenumbers where k minkk max, and a min is the amplitude of the modes at k min and k max in the parabolic case (we use a min = 0 in the code). The frequency grid is defined using frequencies from k x = n x2π/L x in the x direction, where n x ∈ [0, 20] is an integer and L x is the box size in the x-direction, while k y = n y2π/L y and k z = n z2π/L z with n y ∈ [0, 8] and n z ∈ [0, 8]. We then set up four modes for each combination of n x, n y, and n z, corresponding to [k x, k y, k z], [k x, −k y, k z], [k x, k y, −k z], and [k x, −k y, −k z]. That is, we effectively sum from [− (N − 1)/2, (N − 1)/2] in the k y and k z directions in the discrete Fourier transform, where N = max(n x) is the number of frequencies. The default values for k min and k max are 2π and 6π, respectively, corresponding to large-scale driving of only the first few Fourier modes, so with default settings there are 112 non-zero Fourier modes. The maximum number of modes, defining the array sizes needed to store the stirring information, is currently set to 1,000.

We apply the forcing to the particles by computing the discrete Fourier transform over the stirring modes directly, i.e.

(117) $$\begin{equation} {\bm a}_{{\rm forcing}, a} = f_{\rm sol} \sum _{m=1}^{n_{\rm modes}} C_{m} \left[ {\bm A}_{m} \cos ({\bm k}\cdot {\bm r}_{a}) - {\bm B}_{m}\sin ({\bm k}\cdot {\bm r}_{a}) \right], \end{equation}$$

where the factor f sol is defined from the solenoidal weight, w, according to

(118) $$\begin{equation} f_{\rm sol} = \sqrt{\frac{3}{n_{\rm dim}}} \sqrt{ \frac{3}{1 - 2w + n_{\rm dim} w^{2}}}, \end{equation}$$

such that the rms acceleration is the same irrespective of the value of w. We default to purely solenoidal forcing (w = 1), with the factor f sol thus equal to $\sqrt{3/2}$ by default. For individual timesteps, we update the driving force only when a particle is active.

To aid reproducibility, it is often desirable to pre-generate the entire driving sequence prior to performing a calculation, which can then be written to a file and read back at runtime. This was the procedure used in Price & Federrath (Reference Price and Federrath2010), Tricco et al. (Reference Tricco, Price and Federrath2016b), and Tricco et al. (Reference Tricco, Price and Laibe2017).

2.6. Accretion disc viscosity

Accretion disc viscosity is implemented in Phantom via two different approaches, as described by Lodato & Price (Reference Lodato and Price2010).

2.6.1. Disc viscosity using the shock viscosity term

The default approach is to adapt the shock viscosity term to represent a Shakura & Sunyaev (Reference Shakura and Sunyaev1973) α-viscosity, as originally proposed by Artymowicz & Lubow (Reference Artymowicz and Lubow1994) and Murray (Reference Murray1996). The key is to note that (39) and (40) represent a Navier–Stokes viscosity term with a fixed ratio between the bulk and shear viscosity terms (e.g. Murray Reference Murray1996; Jubelgas, Springel, & Dolag Reference Jubelgas, Springel and Dolag2004; Lodato & Price Reference Lodato and Price2010; Price Reference Price2012b; Meru & Bate Reference Meru and Bate2012). In particular, it can be shown (e.g. Español & Revenga Reference Español and Revenga2003) that

(119) $$\begin{equation} \sum _{b} \frac{m_{b}}{\overline{\rho }_{ab}} ({\bm v}_{ab}\cdot \hat{\bm r}_{ab} )\frac{\nabla _{a} W_{ab}}{\vert r_{ab}\vert } \approx \frac{1}{5} \nabla \left( \nabla \cdot {\bm v}\right) + \frac{1}{10}\nabla ^{2} {\bm v}, \end{equation}$$

where $\overline{\rho }_{ab}$ is some appropriate average of the density. This enables the artificial viscosity term, (39), to be translated into the equivalent Navier–Stokes terms. In order for the artificial viscosity to represent a disc viscosity, we make the following modifications (Lodato & Price Reference Lodato and Price2010):

  1. 1. The viscosity term is applied for both approaching and receding particles.

  2. 2. The speed v sig is set equal to c s.

  3. 3. A constant αAV is adopted, turning off shock detection switches (Section 2.2.9).

  4. 4. The viscosity term is multiplied by a factor h/|r ab|.

The net result is that (40) becomes

(120) $$\begin{equation} q^{a}_{ab} = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}- \displaystyle\frac{\rho _{a} h_{a}}{2 \vert r_{ab}\vert } \left(\alpha ^{\rm AV} c_{{\rm s}, a} + \beta _{\rm AV} \vert {\bm v}_{ab} \cdot \hat{\bm r}_{ab}\vert \right) {\bm v}_{ab} \cdot \hat{\bm r}_{ab}, & {\bm v}_{ab} \cdot \hat{\bm r}_{ab} < 0 \\[3pt] - \frac{\rho _{a} h_{a}}{2 \vert r_{ab}\vert } \alpha ^{\rm AV} c_{{\rm s}, a} {\bm v}_{ab} \cdot \hat{\bm r}_{ab}. & \text{otherwise.} \end{array}\right. \end{equation}$$

With the above modifications, the shear and bulk coefficients can be translated using (119) to give (e.g. Monaghan Reference Monaghan2005; Lodato & Price Reference Lodato and Price2010; Meru & Bate Reference Meru and Bate2012)

(121) $$\begin{eqnarray} \nu _{\rm AV} \approx \frac{1}{10} \alpha ^{\rm AV} c_{\rm s} h, \end{eqnarray}$$
(122) $$\begin{eqnarray} \zeta _{\rm AV} = \frac{5}{3} \nu _{\rm AV} \approx \frac{1}{6} \alpha ^{\rm AV} c_{\rm s} h. \end{eqnarray}$$

The Shakura–Sunyaev prescription is

(123) $$\begin{equation} \nu = \alpha _{\rm SS} c_{\rm s} H, \end{equation}$$

where H is the scale height. This implies that αSS may be determined from αAV using

(124) $$\begin{equation} \alpha _{\rm SS} \approx \frac{\alpha ^{\rm AV}}{10} \frac{\langle h \rangle }{H}, \end{equation}$$

where 〈h〉 is the mean smoothing length on particles in a cylindrical ring at a given radius.

In practice, this means that one must uniformly resolve the scale height in order to obtain a constant αSS in the disc. We have achieved this in simulations to date by choosing the initial surface density profile and the power-law index of the temperature profile (when using a locally isothermal EOS) to ensure that this is the case (Lodato & Pringle Reference Lodato and Pringle2007). Confirmation that the scaling provided by (124) is correct is shown in Figure 4 of Lodato & Price (Reference Lodato and Price2010) and is checked automatically in the Phantom test suite.

In the original implementation (Lodato & Price Reference Lodato and Price2010), we also set the βAV to zero, but this is dangerous if the disc dynamics are complex as there is nothing to prevent particle penetration (see Section 2.2.7). Hence, in the current version of the code, βAV = 2 by default even if disc viscosity is set, but is only applied to approaching particles (c.f. 120). Applying any component of viscosity to only approaching particles can affect the sign of the precession induced in warped discs (Lodato & Pringle Reference Lodato and Pringle2007), but in general retaining the βAV term is safer with no noticeable effect on the overall dissipation due to the second-order dependence of this term on resolution.

Using αAV to set the disc viscosity has two useful consequences. First, it forces one to consider whether or not the scale height, H, is resolved. Second, knowing the value of αAV is helpful, as αAV ≈ 0.1 represents the lower bound below which a physical viscosity is not resolved in SPH (that is, viscosities smaller than this produce disc spreading independent of the value of αAV, see Bate Reference Bate1995; Meru & Bate Reference Meru and Bate2012), while αAV > 1 constrains the timestep (Section 2.3.2).

2.6.2. Disc viscosity using the Navier–Stokes viscosity terms

An alternative approach is to compute viscous terms directly from the Navier–Stokes equation. Details of how the Navier–Stokes terms are represented are given below (Section 2.7), but for disc viscosity a method for determining the kinematic viscosity is needed, which in turn requires specifying the scale height as a function of radius. We use

(125) $$\begin{equation} H_{a} \equiv \frac{c^{a}_{\rm s}}{\Omega (R_{a})}, \end{equation}$$

where we assume Keplerian rotation $\Omega = \sqrt{GM/R^{3}}$ and c s is obtained for a given particle from the EOS (which for consistency must be either isothermal or locally isothermal). It is important to note that this restricts the application of this approach only to discs where R can be meaningfully defined, excluding, for example, discs around binary stars.

The shear viscosity is then set using

(126) $$\begin{equation} \nu _{a} = \alpha _{\rm SS} c_{\rm s}^{a} H_{a}, \end{equation}$$

where αSS is a pre-defined/input parameter. The timestep is constrained using C visch 2/ν as described in Section 2.7. The advantage to this approach is that the shear viscosity is set directly and does not depend on the smoothing length. However, as found by Lodato & Price (Reference Lodato and Price2010), it remains necessary to apply some bulk viscosity to capture shocks and prevent particle penetration of the disc midplane, so one should apply the shock viscosity as usual. Using a shock-detection switch (Section 2.2.9) means that this is usually not problematic. This formulation of viscosity was used in Facchini et al. (Reference Facchini, Lodato and Price2013).

2.7. Navier–Stokes viscosity

Physical viscosity is implemented as described in Lodato & Price (Reference Lodato and Price2010). Here, (23) and (24) are replaced by the compressible Navier–Stokes equations, i.e.

(127) $$\begin{eqnarray} \frac{{\rm d}v^{i}}{{\rm d}t} &=& -\frac{1}{\rho }\frac{\partial S_{\rm NS}^{ij}}{\partial x^{j}} + \Pi _{\rm shock} + {\bm a}_{\rm ext}({\bm r}, t) \nonumber \\ && +\, {\bm a}_{\rm sink-gas} + {\bm a}_{\rm selfgrav}, \end{eqnarray}$$
(128) $$\begin{eqnarray} \frac{{\rm d}u}{{\rm d}t} = -\frac{P}{\rho } \left(\nabla \cdot {\bm v}\right) + \Lambda _{\rm visc} + \Lambda _{\rm shock} - \frac{\Lambda _{\rm cool}}{\rho }, \end{eqnarray}$$

with the stress tensor given by

(129) $$\begin{equation} S_{\rm NS}^{ij} = \left[ P - \left(\zeta - \frac{2}{3}\eta \right) \frac{\partial v^{k}}{\partial x^{k}} \right] \delta ^{ij} - \eta \left( \frac{\partial v^{i}}{\partial x^{j}} + \frac{\partial v^{j}}{\partial x^{i}}\right), \end{equation}$$

where δij is the Kronecker delta, and ζ and η are the bulk and shear viscosity coefficients, related to the volume and kinematic shear viscosity coefficients by ζv ≡ ζ/ρ and ν ≡ η/ρ.

2.7.1. Physical viscosity using two first derivatives

As there is no clear consensus on the best way to implement physical viscosity in SPH, Phantom currently contains two implementations. The simplest is to use two first derivatives, which is essentially that proposed by Flebbe et al. (Reference Flebbe, Muenzel, Herold, Riffert and Ruder1994), Watkins et al. (Reference Watkins, Bhattal, Francis, Turner and Whitworth1996), and Sijacki & Springel (Reference Sijacki and Springel2006). In this case, (127) is discretised in the standard manner using

(130) $$\begin{eqnarray} \frac{{\rm d}v^{i}_{a}}{{\rm d}t} &=& -\sum _{b} m_{b}\left[ \frac{S^{ij}_{{\rm NS}, a}}{\Omega _{a}\rho _{a}^{2}} \frac{\partial W_{ab}(h_{a})}{\partial x^{j}_{a}} + \frac{S^{ij}_{{\rm NS}, b}}{\Omega _{b}\rho _{b}^{2}} \frac{\partial W_{ab}(h_{b})}{\partial x^{j}_{a}}\right] \nonumber \\ && +\, \Pi ^{i}_{\rm shock} + a^{i}_{\rm ext}({\bm r}, t) + a^{i}_{\rm sink-gas} + a^{i}_{\rm selfgrav}, \end{eqnarray}$$

where the velocity gradients are computed during the density loop using

(131) $$\begin{equation} \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}} = -\frac{1}{\Omega _{a}\rho _{a}} \sum _{b} m_{b} v_{ab}^{i} \nabla _a^{j} W_{ab} (h_{a}). \end{equation}$$

Importantly, the differenced SPH operator is used in (131), whereas (130) uses the symmetric gradient operator. The use of conjugate operatorsFootnote 5 is a common requirement in SPH in order to guarantee energy conservation and a positive definite entropy increase from dissipative terms (e.g. Price Reference Price2010; Tricco & Price Reference Tricco and Price2012). Total energy conservation means that

(132) $$\begin{equation} \frac{{\rm d}E}{{\rm d}t} = \sum _{a} m_{a} \left( \frac{{\rm d}u_{a}}{{\rm d}t} + v^{i}_{a} \frac{{\rm d}v^{i}_{a}}{{\rm d}t} \right) = 0. \end{equation}$$

This implies a contribution to the thermal energy equation given by

(133) $$\begin{equation} \frac{{\rm d}u_{a}}{{\rm d}t} = \frac{S_{{\rm NS}, a}^{ij}}{\Omega _{a} \rho _{a}^{2}} \sum _{b} m_{b} v_{ab}^{i} \nabla _a^{j} W_{ab} (h_{a}), \end{equation}$$

which can be seen to reduce to (24) in the inviscid case (Sij NS = Pδij), but in general is an SPH expression for

(134) $$\begin{equation} \frac{{\rm d}u_{a}}{{\rm d}t} = -\frac{S_{{\rm NS}, a}^{ij}}{\rho _{a}} \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}}. \end{equation}$$

Using Sij NS = S NSji, we have

(135) $$\begin{equation} \frac{{\rm d}u_{a}}{{\rm d}t} = -\frac{1}{2} \frac{S_{{\rm NS}, a}^{ij}}{\rho _{a}} \left( \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}} + \frac{\partial v_{a}^{j}}{\partial x_{a}^{i}} \right), \end{equation}$$

which, using (129), gives

(136) $$\begin{eqnarray} \Lambda _{\rm visc} = \left(\zeta _{v, a} - \frac{2}{3} \nu _{a} \right) (\nabla \cdot {\bm v})_{a}^{2} + \frac{\nu _{a}}{2} \left( \frac{\partial v_{a}^{i}}{\partial x_{a}^{j}} + \frac{\partial v_{a}^{j}}{\partial x_{a}^{i}} \right)^{2}. \end{eqnarray}$$

By the square in the last term, we mean the tensor summation ∑jiT ijT ij, where Tij ≡ ∂via/∂xja + ∂va j/∂xia. The heating term is therefore positive definite provided that the velocity gradients and divergence are computed using the difference operator (131), both in (136) and when computing the stress tensor (129).

The main disadvantage of the first derivatives approach is that it requires storing the strain tensor for each particle, i.e. six additional quantities when taking account of symmetries.

2.7.2. Physical viscosity with direct second derivatives

The second approach is to use SPH second derivative operators directly. Here we use modified versions of the identities given by Español & Revenga (Reference Español and Revenga2003) (see also Monaghan Reference Monaghan2005; Price Reference Price2012a), namely

(137) $$\begin{eqnarray} \hspace*{-38pt}&&\nabla \left[A (\nabla \cdot {\bm v})\right] \approx \nonumber \\ \hspace*{-38pt}&&\quad-\,\sum _{b} m_{b} \left[\frac{A_{a}}{\rho _{a}} G_{ab}(h_{a})+ \frac{A_{b}}{\rho _{b}} G_{ab}(h_{b})\right] ({\bm v}_{ab}\cdot \hat{\bm r}_{ab}) \hat{\bm r}_{ab} , \end{eqnarray}$$
(138) $$\begin{eqnarray} \hspace*{-20pt}\nabla \cdot (C \nabla {\bm v}) \approx -\,\sum _{b} m_{b} \left[\frac{C_{a}}{\rho _{a}} G_{ab}(h_{a}) + \frac{C_{b}}{\rho _{b}} G_{ab}(h_{b})\right]{\bm v}_{ab}, \end{eqnarray}$$

where G ab ≡ −2F ab/|r ab|, i.e. the scalar part of the kernel gradient divided by the particle separation, which can be thought of as equivalent to defining a new ‘second derivative kernel’ (Brookshaw Reference Brookshaw1985, Reference Brookshaw1994; Price Reference Price2012a; Price & Laibe Reference Price and Laibe2015a).

From the compressible Navier–Stokes equations, (127) with (129), the coefficients in these two terms are

(139) $$\begin{eqnarray} A \equiv \frac{1}{2} \left(\zeta + \frac{\eta }{3}\right), \end{eqnarray}$$
(140) $$\begin{eqnarray} C \equiv \frac{1}{2} \eta , \end{eqnarray}$$

so that we can simply use

(141) $$\begin{eqnarray} \left(\frac{{\rm d}v^{i}_{a}}{{\rm d}t}\right)_{\rm visc} &=& \sum_b m_b \left[ \frac{\tau_a}{\rho_a} G_{ab}(h_a) + \frac{\tau_b}{\rho_b} G_{ab}(h_b) \right] ({\bm v}_{ab}\cdot \hat{\bm r}_{ab})\hat{r}^i_{ab} \nonumber\\ &&+\, \sum_b m_b \left[ \frac{\kappa_a}{\rho_a} G_{ab}(h_a) + \frac{\kappa_b}{\rho_b} G_{ab}(h_b)\right] v_{ab}^i \end{eqnarray}$$

where

(142) $$\begin{eqnarray} \tau = \frac{5}{2} A, \end{eqnarray}$$
(143) $$\begin{eqnarray} \kappa = \left(C - \frac{A}{2}\right). \end{eqnarray}$$

The corresponding heating terms in the thermal energy equation are given by

(144) $$\begin{eqnarray} \Lambda _{\rm visc} &=& \frac{\tau _{a}}{\rho _{a}} \sum _{b} m_{b} ({\bm v}_{ab}\cdot \hat{\bm r}_{ab})^{2} G_{ab}(h_{a}) \nonumber \\ && +\, \frac{\kappa _{a}}{\rho _{a}} \sum _{b} m_{b} ({\bm v}_{ab})^{2} G_{ab}(h_{a}). \end{eqnarray}$$

This is the default formulation of Navier–Stokes viscosity in the code since it does not require additional storage. In practice, we have found little difference between the two formulations of physical viscosity, but this would benefit from a detailed study. In general one might expect the two first derivatives formulation to offer a less noisy estimate at the cost of additional storage. However, direct second derivatives are the method used in ‘Smoothed Dissipative Particle Dynamics’ (Español & Revenga Reference Español and Revenga2003).

2.7.3. Timestep constraint

Both approaches to physical viscosity use explicit timestepping, and therefore imply a constraint on the timestep given by

(145) $$\begin{equation} \Delta t^{a}_{\rm visc} \equiv C_{\rm visc} \frac{h_{a}^{2}}{\nu _{a}}, \end{equation}$$

where C visc = 0.25 by default (Brookshaw Reference Brookshaw1994). When physical viscosity is turned on, this constraint is included with other timestep constraints according to (77).

2.7.4. Physical viscosity and the tensile instability

Caution is required in the use of physical viscosity at high Mach number, since negative stress can lead to the tensile instability (Morris Reference Morris1996b; Monaghan Reference Monaghan2000; Gray, Monaghan, & Swift Reference Gray, Monaghan and Swift2001). For subsonic applications, this is usually not a problem since the strain tensor and velocity divergence are small compared to the pressure. In the current code, we simply emit a warning if physical viscosity leads to negative stresses during the calculation, but this would benefit from a detailed study.

2.7.5. Physical viscosity and angular momentum conservation

Neither method for physical viscosity exactly conserves angular momentum because the force is no longer directed along the line of sight joining the particles. However, the error is usually small (see discussion in Bonet & Lok Reference Bonet and Lok1999, Section 5 of Price & Monaghan Reference Price and Monaghan2004b or Hu & Adams Reference Hu and Adams2006). Recently, Müller, Fedosov, & Gompper (Reference Müller, Fedosov and Gompper2015) have proposed an algorithm for physical viscosity in SPH that explicitly conserves angular momentum by tracking particle spin, which may be worth investigating.

2.8. Sink particles

Sink particles were introduced into SPH by Bate, Bonnell, & Price (Reference Bate, Bonnell and Price1995) in order to follow star formation simulations beyond the point of fragmentation. In Phantom, these are treated separately to the SPH particles, and interact with other particles, including other sink particles, only via gravity. The differences with other point-mass particles implemented in the code (e.g. dust, stars, and dark matter) are that (i) the gravitational interaction is computed using a direct N 2 summation which is not softened by default (i.e., the N-body algorithm is collisional); (ii) they are allowed to accrete gas; and (iii) they store the accreted angular momentum and other extended properties, such as the accreted mass. Sink particles are evolved in time using the RESPA algorithm (Section 2.3.3), which is second-order accurate, symplectic, and allows sink particles to evolve on shorter timesteps compared to SPH particles.

2.8.1. Sink particle accelerations

The equations of motion for a given sink particle, i, are

(146) $$\begin{eqnarray} \frac{{\rm d}{\bm v}_{i}}{{\rm d}t} &= & -\sum ^{N_{\rm sink}}_{j=1} G M_{j} \phi ^{\prime }_{ij} (\epsilon ) \hat{\bm r}_{ij} \nonumber \\ && -\, \sum ^{N_{\rm part}}_{b=1} G m_{b} \phi ^{\prime }_{ib} (\epsilon _{ib}) \hat{\bm r}_{ib}, \end{eqnarray}$$

where ϕ′ab is the softening kernel (Section 2.12.2), N part is the total number of gas particles, and N sink is the total number of sink particles. The sink–gas softening length, εib, is defined as the maximum of the (fixed) softening length defined for the sink particles, ε, and the softening length of the gas particle, εb. That is, εib ≡ max(ε, εb). SPH particles receive a corresponding acceleration

(147) $$\begin{equation} {\bm a}^{a}_{\rm sink-gas} = -\sum ^{N_{\rm sink}}_{j=1} G M_{j} \phi ^{\prime }_{aj} (\epsilon _{aj}) \hat{\bm r}_{aj}. \end{equation}$$

Softening of sink–gas interactions is not applied if the softening length for sink particles is set to zero, in which case the sink–gas accelerations reduce simply to

(148) $$\begin{equation} {\bm a}^{a}_{\rm sink-gas} = -\sum ^{N_{\rm sink}}_{j=1} \frac{G M_{j}}{\vert {\bm r}_{a} - {\bm r}_{j}\vert ^{3}} {\bm r}_{aj}. \end{equation}$$

This is the default behaviour when sink particles are used in the code. Softening of sink–gas interactions is useful to model a point-mass particle that does not accrete gas (e.g. by setting the accretion radius to zero). For example, we used a softened sink particle to simulate the core of the red giant in Iaconi et al. (Reference Iaconi, Reichardt, Staff, De Marco, Passy, Price, Wurster and Herwig2017). The sink–sink interaction is unsoftened by default (ε = 0), giving the usual

(149) $$\begin{equation} {\bm a}^{i}_{\rm sink-sink} = -\sum ^{N_{\rm sink}}_{j=1} \frac{G M_{j}}{ \vert {\bm r}_{i} - {\bm r}_{j}\vert ^{3}} {\bm r}_{ij}. \end{equation}$$

Caution is required when integrating tight binary or multiple systems when ε = 0 to ensure that the timestep conditions (Section 2.3.2) are strict enough.

2.8.2. Accretion onto sink particles

Accretion of gas particles onto a sink particle occurs when a gas particle passes a series of accretion checks within the accretion radius r acc of a sink particle (set in the initial conditions or when the sink particle is created, see Section 2.8.4). First, a gas particle is indiscriminately accreted without undergoing any additional checks if it falls within f accr acc, where 0 ⩽ f acc ⩽ 1 (default f acc = 0.8). In the boundary region f accr acc < r < r acc, a gas particle undergoes accretion if

  1. 1. $\vert {\bm L}_{ai} \vert < \vert {\bm L}_{\rm acc} \vert$, that is, its specific angular momentum is less than that of a Keplerian orbit at r acc,

  2. 2. $e = \frac{v_{ai}^2}{2} - \frac{GM_i}{r_{ai}} < 0$, i.e., it is gravitationally bound to the sink particle, and

  3. 3. e for this gas–sink pair is smaller than e with any other sink particle, that is, out of all sink particles, the gas particle is most bound to this one.

In the above conditions, ${\bm L}_{ai}$ is the relative specific angular momentum of the gas–sink pair, ai, defined by

(150) $$\begin{eqnarray} \vert {\bm L}_{ai}^2\vert & \equiv& \vert {\bm r}_{ai} \times {\bm v}_{ai}\vert ^2 \nonumber \\ & = &r_{ai}^2v_{ai}^2 - \left({\bm r}_{ai}\cdot {\bm v}_{ai}\right)^2, \end{eqnarray}$$

while $\vert {\bm L}_{\rm acc} \vert = r_{\rm acc}^2 \Omega _{\rm acc}$ is the angular momentum at r acc, where $\Omega _{\rm acc} = \sqrt{GM_i / r_{ai}^{3}}$ is the Keplerian angular speed at r acc, v ai, and r ai are the relative velocity and position, respectively, and M i is the mass of the sink particle.

When a particle, a, passes the accretion checks, then the mass, position, velocity, acceleration, and spin angular momentum of the sink particle are updated according to

(151) $$\begin{eqnarray} {\bm r}_{i} = \frac{({\bm r}_{a} m_{a} + {\bm r}_{i} M_{i})}{M_{i} + m_{a}}, \end{eqnarray}$$
(152) $$\begin{eqnarray} {\bm v}_{i} = \frac{({\bm v}_{a} m_{a} + {\bm v}_{i} M_{i})}{M_{i} + m_{a}}, \end{eqnarray}$$
(153) $$\begin{eqnarray} {\bm a}_{i} = \frac{({\bm a}_{a} m_{a} + {\bm a}_{i} M_{i})}{M_{i} + m_{a}}, \end{eqnarray}$$
(154) $$\begin{eqnarray} {\bm S}_{i} = {\bm S}_{i} +\frac{m_{a} M_{i}}{M_{i} + m_{a}} \left[ \left({\bm r}_{a} - {\bm r}_{i}\right) \times \left({\bm v}_{a} - {\bm v}_{i}\right)\right], \end{eqnarray}$$
(155) $$\begin{eqnarray} M_{i} = M_{i} + m_{a}. \end{eqnarray}$$

This ensures that mass, linear momentum, and angular momentum (but not energy) are conserved by the accretion process. The accreted mass as well as the total mass for each sink particle is stored to avoid problems with round-off error in the case where the particle masses are much smaller than the sink mass. Accreted particles are tagged by setting their smoothing lengths negative. Those particles with h ⩽ 0 are subsequently excluded when the kd-tree is built.

2.8.3. Sink particle boundary conditions

No special sink particle boundary conditions are implemented in Phantom at present. More advanced boundary conditions to regularise the density, pressure, and angular momentum near a sink have been proposed by Bate et al. (Reference Bate, Bonnell and Price1995) and used in Bate & Bonnell (Reference Bate and Bonnell1997), and proposed again more recently by Hubber, Walch, & Whitworth (Reference Hubber, Walch and Whitworth2013b). While these conditions help to regularise the flow near the sink particle, they can also cause problems—particularly the angular momentum boundary condition if the disc near the sink particle has complicated structure such as spiral density waves (Bate 2014, private communication). Often it is more cost effective to simply reduce the accretion radius of the sink. This may change in future code versions.

2.8.4. Dynamic sink particle creation

As described in Bate et al. (Reference Bate, Bonnell and Price1995), it is also possible to create sink particles on-the-fly provided certain physical conditions are met and self-gravity is turned on (Section 2.12). The primary conditions required for sink particle formation are that the density of a given particle exceeds some threshold physical density somewhere in the domain, and that this density peak occurs more than a critical distance r crit from an existing sink. Once these conditions are met on a particular particle, a, the creation of a new sink particle occurs by passing the following conditions (Bate et al. Reference Bate, Bonnell and Price1995):

  1. 1. The particle is a gas particle.

  2. 2. $\nabla \cdot {\bm v}_{a} \le 0$, that is, gas surrounding the particle is at rest or collapsing.

  3. 3. h a < r acc/2, i.e., the smoothing length of the particle is less than half of the accretion radius.

  4. 4. All neighbours within r acc are currently active.

  5. 5. The ratio of thermal to gravitational energy of particles within r acc, αJ, satisfies αJ ⩽ 1/2.

  6. 6. αJ + βrot ⩽ 1, where βrot = |e rot|/|e grav| is the ratio of rotational energy to the magnitude of the gravitational energy for particles within r acc.

  7. 7. e tot < 0, that is, the total energy of particles within r acc is negative (i.e. the clump is gravitationally bound).

  8. 8. The particle is at a local potential minimum, i.e. Φ is less than Φ computed on all other particles within r acc (Federrath et al. Reference Federrath, Banerjee, Clark and Klessen2010b).

A new sink particle is created at the position of particle a if these checks are passed, and immediately the particles within r acc are accreted by calling the routine described in Section 2.8.2. The checks above are the same as those in Bate et al. (Reference Bate, Bonnell and Price1995), with the additional check from Federrath et al. (Reference Federrath, Banerjee, Clark and Klessen2010b) to ensure that sink particles are only created in a local minimum of the gravitational potential.

The various energies used to evaluate the criteria above are computed according to

(156) $$\begin{eqnarray} \hspace*{-26.2pt}e_{\rm kin} = \frac{1}{2} \sum ^{N<r_{\rm acc}}_{b=1} m_{b} ({\bm v}_{b} - {\bm v}_{a})^{2},\hspace*{26.2pt} \end{eqnarray}$$
(157) $$\begin{eqnarray} \hspace*{-49.14pt}e_{\rm therm} = \sum ^{N<r_{\rm acc}}_{b=1} m_{b} u_{b},\hspace*{49.14pt} \end{eqnarray}$$
(158) $$\begin{eqnarray} e_{\rm grav} = -\frac{1}{2}\sum ^{N<r_{\rm acc}}_{b=1}\sum ^{N<r_{\rm acc}}_{c=b} G m_{b} m_{c} , \nonumber \\ && \times\, \left[\phi (\vert {\bm r}_{b} - {\bm r}_{c}\vert , h_b) + \phi (\vert {\bm r}_{b} - {\bm r}_{c}\vert , h_c)\right], \end{eqnarray}$$
(159) $$\begin{eqnarray} e_{\rm tot} = e_{\rm kin} + e_{\rm therm} + e_{\rm grav}, \end{eqnarray}$$
(160) $$\begin{eqnarray} e_{\rm rot} \equiv \sqrt{e^{2}_{{\rm rot}, x} + e^{2}_{{\rm rot}, y} + e^{2}_{{\rm rot}, z}}, \end{eqnarray}$$
(161) $$\begin{eqnarray} e_{{\rm rot}, x} \equiv \frac{1}{2} \sum ^{N<r_{\rm acc}}_{b=1} m_{b} \frac{L^{2}_{{ab},x}}{\sqrt{(y_{a} - y_{b})^{2} + (z_{a} - z_{b})^{2}}}, \end{eqnarray}$$
(162) $$\begin{eqnarray} e_{{\rm rot}, y} \equiv \frac{1}{2} \sum ^{N<r_{\rm acc}}_{b=1} m_{b} \frac{L^{2}_{{ab},y}}{\sqrt{(x_{a} - x_{b})^{2} + (z_{a} - z_{b})^{2}}}, \end{eqnarray}$$
(163) $$\begin{eqnarray} e_{{\rm rot}, z} \equiv \frac{1}{2} \sum ^{N<r_{\rm acc}}_{b=1} m_{b} \frac{L^{2}_{{ab},z}}{\sqrt{(x_{a} - x_{b})^{2} + (y_{a} - y_{b})^{2}}}, \end{eqnarray}$$

where ${\bm L}_{ab} \equiv ({\bm r}_{a} - {\bm r}_{b}) \times ({\bm v}_{a} - {\bm v}_{b})$ is the specific angular momentum between a pair of particles, and ϕ is the gravitational softening kernel (defined in Section 2.12), which has units of inverse length. Adding the contribution from all pairs, bc, within the clump is required to obtain the total potential of the clump.

2.8.5. Sink particle timesteps

Sink particles are integrated together with a global, but adaptive, timestep, following the inner loop of the RESPA algorithm given in (80)–(83) corresponding to a second-order Leapfrog integration. The timestep is controlled by the minimum of the sink–gas timestep, (75), and a sink–sink timestep (Dehnen & Read Reference Dehnen and Read2011)

(164) $$\begin{equation} \Delta t_{\rm sink-sink} \equiv \min _{i} \left( C_{\rm force} \eta _\Phi \sqrt{\frac{\vert \Phi ^{\rm sink-sink}_{i}\vert }{\vert \nabla \Phi ^{\rm sink-sink}_i \vert ^{2}}} \right), \end{equation}$$

where the potential and gradient include other sink particles, plus any external potentials applied to sink particles except the sink–gas potential. We set $\eta _\Phi = 0.05$ by default, resulting in ~300–500 steps per orbit for a binary orbit with the default C force = 0.25 (see Section 5.5.1).

More accurate integrators such as the fourth-order Hermite scheme (Makino & Aarseth Reference Makino and Aarseth1992) or the fourth-order symplectic schemes proposed by Omelyan, Mryglod, & Folk (Reference Omelyan, Mryglod and Folk2002) or Chin & Chen (Reference Chin and Chen2005) are not yet implemented in Phantom, but it would be a worthwhile effort to incorporate one of these in a future code version. See Hubber et al. (Reference Hubber, Allison, Smith and Goodwin2013a) for a recent implementation of a fourth-order Hermite scheme for sink particles in SPH.

2.9. Stellar physics

A tabulated EOS can be used to take account of the departure from an ideal gas, for example, due to changes in ionisation or molecular dissociation and recombination. This tabulated EOS in Phantom is adapted from the logP gasT EOS tables provided with the open source package Modules for Experiments in Stellar Astrophysics mesa (Paxton et al. Reference Paxton, Bildsten, Dotter, Herwig, Lesaffre and Timmes2011). Details of the data, originally compiled from blends of equations of state from Saumon, Chabrier, & van Horn (Reference Saumon, Chabrier and van Horn1995) (SCVH), Timmes & Swesty (Reference Timmes and Swesty2000), Rogers & Nayfonov (Reference Rogers and Nayfonov2002, also the 2005 update), Potekhin & Chabrier (Reference Potekhin and Chabrier2010) and for an ideal gas, are outlined by Paxton et al. (Reference Paxton, Bildsten, Dotter, Herwig, Lesaffre and Timmes2011).

In our implementation (adapted from original routines for the Music code; Goffrey et al. Reference Goffrey2017), we compute the pressure and other required EOS variables for a particular composition by interpolation between sets of tables for different hydrogen abundance X = 0.0, 0.2, 0.4, 0.6, 0.8 and metallicity Z = 0.0, 0.02, 0.04. Pressure is calculated with bicubic interpolation, and Γ1 ≡ ∂lnP/∂lnρ|s with bilinear interpolation, in logu and logV ≡ logρ − 0.7logu + 20. The tables are currently valid in the ranges 10.5 ⩽ logu ⩽ 17.5 and 0.0 ⩽ logV ⩽ 14.0. Values requested outside the tables are currently computed by linear extrapolation. This triggers a warning to the user.

We have not tested the thermodynamic consistency of our interpolation scheme from the tables (Timmes & Arnett Reference Timmes and Arnett1999).

2.10. Magnetohydrodynamics

Phantom implements the smoothed particle magnetohydrodynamics (SPMHD) algorithm described in Price (Reference Price2012a) and Tricco & Price (Reference Tricco and Price2012, Reference Tricco and Price2013), based on the original work by Phillips & Monaghan (Reference Phillips and Monaghan1985) and Price & Monaghan (Reference Price and Monaghan2004a, Reference Price and Monaghan2004b, Reference Price and Monaghan2005). Phantom was previously used to test a vector potential formulation (Price Reference Price2010), but this algorithm has been subsequently removed from the code due to its poor performance (see Price Reference Price2010).

The important difference between Phantom and the gadget implementation of SPMHD (Dolag & Stasyszyn Reference Dolag and Stasyszyn2009; Bürzle et al. Reference Bürzle, Clark, Stasyszyn, Greif, Dolag, Klessen and Nielaba2011a, Reference Bürzle, Clark, Stasyszyn, Dolag and Klessen2011b), which also implements the core algorithms from Price & Monaghan (Reference Price and Monaghan2004a, Reference Price and Monaghan2004b, Reference Price and Monaghan2005), is our use of the divergence-cleaning algorithm from Tricco & Price (Reference Tricco and Price2012, Reference Tricco and Price2013) and Tricco, Price, & Bate (Reference Tricco, Price and Bate2016a). This is vital for preserving the divergence-free (no monopoles) condition on the magnetic field.

For recent applications of Phantom to MHD problems, see e.g. Tricco et al. (Reference Tricco, Price and Federrath2016b), Dobbs et al. (Reference Dobbs, Price, Pettitt, Bate and Tricco2016), Bonnerot et al. (Reference Bonnerot, Price, Lodato and Rossi2017), Forgan, Price, & Bonnell (Reference Forgan, Price and Bonnell2017), and Wurster et al. (Reference Wurster, Price and Bate2016, Reference Wurster, Price and Bate2017).

2.10.1. Equations of magnetohydrodynamics

Phantom solves the equations of MHD in the form

(165) $$\begin{eqnarray} \frac{{\rm d}v^{i}}{{\rm d}t} &=& -\frac{1}{\rho }\frac{\partial M^{ij}}{\partial x^{j}} + \Pi _{\rm shock} + f^{i}_{\rm divB} + a^{i}_{\rm ext} \nonumber \\ && +\, a^{i}_{\rm sink-gas} + a^{i}_{\rm selfgrav}, \end{eqnarray}$$
(166) $$\begin{eqnarray} \frac{{\rm d}u}{{\rm d}t} = -\frac{P}{\rho } \left(\nabla \cdot {\bm v}\right) + \Lambda _{\rm shock} - \frac{\Lambda _{\rm cool}}{\rho }, \end{eqnarray}$$
(167) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left(\frac{\bm B}{\rho } \right) = \frac{1}{\rho } \left[ \left({\bm B}\cdot \nabla \right){\bm v} - \nabla \psi + \mathcal {\bm D}_{\rm diss} \right], \end{eqnarray}$$
(168) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left( \frac{\psi }{c_{\rm h}} \right) = -c_{\rm h} \left(\nabla \cdot {\bm B}\right) - \frac{1}{2} \frac{\psi }{c_{\rm h}} \left(\nabla \cdot {\bm v}\right) - \frac{\psi }{c_{\rm h}\tau _{\rm c}}, \end{eqnarray}$$

where ${\bm B}$ is the magnetic field, ψ is a scalar used to control the divergence error in the magnetic field (see Section 2.10.8), and $\mathcal {D}_{\rm diss}$ represents magnetic dissipation terms (Sections 2.10.5 and 2.11). The Maxwell stress tensor, M ij, is given by

(169) $$\begin{equation} M^{ij} = \left( P + \frac{1}{2} \frac{B^{2}}{\mu _{0}} \right) \delta ^{ij} - \frac{B^{i} B^{j}}{\mu _{0}}, \end{equation}$$

where δij is the Kronecker delta and μ0 is the permeability of free space. A source term related to the numerically induced divergence of the magnetic field, given by

(170) $$\begin{equation} f^{i}_{\rm divB} \equiv - \frac{B^{i}}{\rho }\left(\nabla \cdot {\bm B}\right) \end{equation}$$

is necessary to prevent the tensile instability in SPMHD (Phillips & Monaghan Reference Phillips and Monaghan1985; Monaghan Reference Monaghan2000; Børve, Omang, & Trulsen Reference Børve, Omang and Trulsen2001; Price Reference Price2012a). With this source term, the equation set for ideal MHD in the absence of the divergence cleaning field, ψ, is formally the same as in the Powell et al. (Reference Powell, Roe, Linde, Gombosi and de Zeeuw1999) eight-wave scheme (Price Reference Price2012a), meaning that the divergence errors in the magnetic field are advected by the flow, but not dissipated, unless cleaning is used.

2.10.2. Discrete equations

The discrete version of (165) follows the same procedure as for physical viscosity (Section 2.7), i.e.

(171) $$\begin{eqnarray} \frac{{\rm d}v^{i}_{a}}{{\rm d}t} &=& -\sum _{b} m_{b}\left[ \frac{M^{ij}_{a}}{\Omega _{a}\rho _{a}^{2}} \frac{\partial W_{ab}(h_{a})}{\partial x^{j}_{a}} + \frac{M^{ij}_{b}}{\Omega _{b}\rho _{b}^{2}} \frac{\partial W_{ab}(h_{b})}{\partial x^{j}_{a}}\right] \nonumber \\ && + \,\Pi _{\rm shock}^{a} + f^{i}_{{\rm divB}, a} + a^{i}_{{\rm ext}, a} + a^i_{\rm sink-gas} \nonumber \\ && +\, a^i_{\rm selfgrav}, \end{eqnarray}$$

where Mija is defined according to (169), f divB is a correction term for stability (discussed below), and accelerations due to external forces are as described in Section 2.4.

Equations (167) and (168) are discretised according to (Price & Monaghan Reference Price and Monaghan2005; Tricco & Price Reference Tricco and Price2012; Tricco et al. Reference Tricco, Price and Bate2016a)

(172) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left( \frac{{\bm B}}{\rho } \right)_a &=& -\frac{1}{\Omega _{a} \rho _{a}^2} \sum _{b} m_{b} {\bm v}_{ab} \left[ {\bm B}_{a} \cdot \nabla _a W_{ab}(h_{a})\right] \nonumber \\ && -\, \sum _{b} m_{b} \left[ \frac{\psi _{a}}{\Omega _{a}\rho _{a}^{2}} \nabla _a W_{ab}(h_{a})+ \frac{\psi _{b}}{\Omega _{b}\rho _{b}^{2}} \nabla _a W_{ab}(h_{b})\right] \nonumber \\ && +\, \frac{1}{\rho _a} \mathcal {D}^{a}_{\rm diss}, \end{eqnarray}$$
(173) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left(\frac{\psi }{c_{\rm h}} \right)_a &=& \frac{c_{\rm h}^a}{\Omega _{a}\rho _{a}} \sum _{b} m_{b} {\bm B}_{ab} \cdot \nabla _a W_{ab}(h_{a}) \nonumber \\ && +\, \frac{\psi _{a}}{2c_{\rm h}^a\Omega _{a}\rho _{a}} \sum _{b} m_{b} {\bm v}_{ab} \cdot \nabla _a W_{ab}(h_{a}) - \frac{\psi _{a}}{c_{\rm h}^a\tau ^{a}_{\rm c}}.\qquad \end{eqnarray}$$

The first term in (173) uses the divergence of the magnetic field discretised according to

(174) $$\begin{equation} (\nabla \cdot {\bm B})_{a} = - \frac{1}{\Omega _{a}\rho _{a}} \sum _{b} m_{b} \left({\bm B}_{a} - {\bm B}_{b} \right) \cdot \nabla _a W_{ab}(h_{a}), \end{equation}$$

which is therefore the operator we use when measuring the divergence error (c.f. Tricco & Price Reference Tricco and Price2012).

2.10.3. Code units

An additional unit is required when magnetic fields are included to describe the unit of magnetic field. We adopt code units such that μ0 = 1, as is common practice. The unit scalings for the magnetic field can be determined from the definition of the Alfvén speed,

(175) $$\begin{equation} v_{\rm A} \equiv \sqrt{\frac{B^{2}}{\mu _{0} \rho }}. \end{equation}$$

Since the Alfvén speed has dimensions of length per unit time, this implies a unit for the magnetic field, u mag, given by

(176) $$\begin{equation} u_{\rm mag} = \left(\frac{\mu _{0} u_{\rm mass}}{u_{\rm dist} u_{\rm time}} \right)^{\frac{1}{2}}. \end{equation}$$

Converting the magnetic field in the code to physical units therefore only requires specifying μ0 in the relevant unit system. In particular, it avoids the differences between SI and cgs units in how the charge unit is defined, since μ0 is dimensionless and equal to 4π in cgs units but has dimensions that involve charge in SI units.

2.10.4. Tensile instability correction

The correction term f divB is necessary to avoid the tensile instability—a numerical instability where particles attract each other along field lines—in the regime where the magnetic pressure exceeds the gas pressure, that is, when plasma $\beta \equiv P / \frac{1}{2} B^2 < 1$ (Phillips & Monaghan Reference Phillips and Monaghan1985). The correction term is computed using the symmetric divergence operator (Børve et al. Reference Børve, Omang and Trulsen2001; Price Reference Price2012a; Tricco & Price Reference Tricco and Price2012)

(177) $$\begin{eqnarray} f^{i}_{{\rm divB}, a} &=& -\hat{B}_{a}^{i} \sum _{b} m_{b} \left[ \frac{{\bm B}_{a} \cdot \nabla _{a} W_{ab} (h_{a})}{\Omega _{a} \rho _{a}^{2}} \right. \nonumber \\ && \left. + \,\frac{{\bm B}_{b} \cdot \nabla _{a} W_{ab} (h_{b})}{\Omega _{b} \rho _{b}^{2}} \right]. \end{eqnarray}$$

Since this term violates momentum conservation to the extent that the $\nabla \cdot {\bm B}$ term is non-zero, several authors have proposed ways to minimise its effect. Børve et al. (Reference Børve, Omang and Trulsen2004) showed that stability could be achieved with $\hat{B}^i = \frac{1}{2} B^i$ and also proposed a scheme for scaling this term to zero for β > 1. Barnes, Kawata, & Wu (Reference Barnes, Kawata and Wu2012) similarly advocated using a factor of $\frac{1}{2}$ in this term. However, Tricco & Price (Reference Tricco and Price2012) showed that this could lead to problematic effects (their Figure 12). In Phantom, we use

(178) $$\begin{equation} \hat{B}^i = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}B^i & \beta < 2, \\ {[}(10 - \beta ) B^i]/8 & 2 < \beta < 10, \\ 0 & \textrm {otherwise} \end{array}\right. \end{equation}$$

to provide numerical stability in the strong field regime while maintaining conservation of momentum when β > 10. This also helps to reduce errors in the MHD wave speed caused by the correction term (Iwasaki Reference Iwasaki2015).

2.10.5. Shock capturing

The shock capturing term in the momentum equation for MHD is identical to (39) and (40) except that the signal speed becomes (Price & Monaghan Reference Price and Monaghan2004a, Reference Price and Monaghan2005; Price Reference Price2012a)

(179) $$\begin{equation} v_{{\rm sig}, a} = \alpha ^{\rm AV}_{a} v_{a} + \beta ^{\rm AV} \vert {\bm v}_{ab}\cdot \hat{\bm r}_{ab} \vert , \end{equation}$$

where

(180) $$\begin{equation} v_{a} = \sqrt{c^{2}_{{\rm s}, a} + v_{{\rm A}, a}^{2}} \end{equation}$$

is the fast magnetosonic speed. Apart from this, the major difference to the hydrodynamic case is the addition of an artificial resistivity term to capture shocks and discontinuities in the magnetic field (i.e. current sheets). This is given by

(181) $$\begin{equation} \mathcal {D}^{a}_{{\rm diss}} = \frac{\rho _{a}}{2} \sum _{b} m_{b} \left[ \frac{v^{B}_{{\rm sig}, a}}{\rho _{a}^{2}} \frac{F_{ab}(h_{a})}{\Omega _{a}} + \frac{v^{B}_{{\rm sig}, b}}{\rho _{b}^{2}} \frac{F_{ab}(h_{b})}{\Omega _{b}} \right]{\bm B}_{ab}, \end{equation}$$

where vB sig, a is an appropriate signal speed (see below) multiplied by a dimensionless coefficient, αB. The corresponding contribution to the thermal energy from the resistivity term in (42) is given by

(182) $$\begin{eqnarray} \Lambda _{\rm artres} &= & - \frac{1}{4} \sum _{b} m_{b} \left[ \frac{v^{B}_{{\rm sig}, a}}{\rho _{a}^{2}} \frac{F_{ab}(h_{a})}{\Omega _{a}} \right. \nonumber \\ && \left. \phantom{-\frac{1}{2}\sum _{b} m_{b}}+ \frac{v^{B}_{{\rm sig}, b}}{\rho _{b}^{2}} \frac{F_{ab}(h_{b})}{\Omega _{b}} \right] ({\bm B}_{ab})^{2}. \end{eqnarray}$$

As with the artificial viscosity, (181) and (182) are the SPH representation of a physical resistivity term, $\eta \nabla ^2 {\bm B}$, but with a coefficient that is proportional to resolution (Price & Monaghan Reference Price and Monaghan2004a). The resistive dissipation rate from the shock capturing term is given by

(183) $$\begin{equation} \eta \approx \frac{1}{2} \alpha ^{B} v_{\rm sig}^{B} \vert r_{ab} \vert , \end{equation}$$

where |r ab|∝h.

2.10.6. Switch to reduce resistivity

Phantom previously used the method proposed by Tricco & Price (Reference Tricco and Price2013) to reduce the dissipation in the magnetic field away from shocks and discontinuities. The signal velocity, vB sig, was set equal to the magnetosonic speed [Equation (180)] multiplied by the dimensionless coefficient αB, which was set according to

(184) $$\begin{equation} \alpha ^{B}_{a} = \min \left( h_{a} \frac{\vert \nabla {\bm B}_{a} \vert }{\vert {\bm B}_{a}\vert }, \alpha ^{B}_{\max }\right), \end{equation}$$

where αBmax = 1.0 by default and $\vert \nabla {\bm B}_{a} \vert$ is the two-norm of the gradient tensor, i.e. the root mean square of all nine components of this tensor. Unlike the viscous dissipation, this is set based on the instantaneous values of h and ${\bm B}$ and there is no source/decay equation involved, as Tricco & Price (Reference Tricco and Price2013) found it to be unnecessary. Since αB is proportional to resolution, from (183), we see that this results in dissipation that is effectively second order (∝h 2). When individual particle timesteps were used, inactive particles retained their value of αB from the last timestep they were active.

More recently, we have found that a better approach, similar to that used for artificial conductivity, is to simply set αB = 1 for all particles and set the signal speed for artificial resistivity according to

(185) $$\begin{equation} v_{\rm sig}^{B} = \vert {\bm v}_{ab} \times \hat{\bm r}_{ab} \vert . \end{equation}$$

We find that this produces better results on all of our tests (Section 5.6), in particular, producing zero numerical dissipation on the current loop advection test (Section 5.6.5). As with the Tricco & Price (Reference Tricco and Price2013) switch, it gives second-order dissipation in the magnetic field (demonstrated in Section 5.6.1, Figure 26). This is now the default treatment for artificial resistivity in Phantom.

2.10.7. Conservation properties

The total energy when using MHD is given by

(186) $$\begin{equation} E = \sum _{a} m_{a} \left( \frac{1}{2} v_{a}^{2} + u_{a} + \Phi _{a} + \frac{1}{2} \frac{B^{2}_{a}}{\mu _{0}\rho _{a}} \right). \end{equation}$$

Hence, total energy conservation, in the absence of divergence cleaning, corresponds to

(187) $$\begin{eqnarray} \frac{{\rm d}E}{{\rm d}t} &=& \sum _{a} m_{a} \left[ {\bm v}_{a} \cdot \frac{{\rm d}{\bm v}_{a}}{{\rm d}t} + \frac{{\rm d}u_{a}}{{\rm d}t} + \frac{{\rm d} \Phi _{a}}{{\rm d}t} \right. \nonumber \\ && +\, \left. \frac{{B}^{2}_{a}}{2 \mu _{0}\rho _{a}^{2}} \frac{{\rm d}\rho _{a}}{{\rm d}t} + \frac{{\bm B}_{a}}{ \mu _{0}} \cdot \frac{{\rm d}}{{\rm d}t} \left( \frac{{\bm B}}{\rho } \right)_a \right] = 0. \end{eqnarray}$$

Neglecting the f divB correction term for the moment, substituting (171), (35), and (172) into (187) with the ideal MHD and shock capturing terms included demonstrates that the total energy is exactly conserved, using the same argument as the one given in Section 2.2.6 (detailed workings can be found in Price & Monaghan Reference Price and Monaghan2004b). The total linear momentum is also exactly conserved following a similar argument as in Section 2.2.11. However, the presence of the f divB correction term, though necessary for numerical stability, violates the conservation of both momentum and energy in the strong field regime (in the weak field regime, it is switched off and conservation is preserved). The severity of this non-conservation is related to the degree in which divergence errors are present in the magnetic field, hence inadequate divergence control (see below) usually manifests as a loss of momentum conservation in the code (see Tricco & Price Reference Tricco and Price2012, for details).

2.10.8. Divergence cleaning

We adopt the ‘constrained’ hyperbolic/parabolic divergence cleaning algorithm described by Tricco & Price (Reference Tricco and Price2012) and Tricco et al. (Reference Tricco, Price and Bate2016a) to preserve the divergence-free condition on the magnetic field. This formulation addresses a number of issues with earlier formulations by Dedner et al. (Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002) and Price & Monaghan (Reference Price and Monaghan2005).

The main idea of the scheme is to propagate divergence errors according to a damped wave equation (Dedner et al. Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002; Price & Monaghan Reference Price and Monaghan2005). This is facilitated by introducing a new scalar field, ψ, which is coupled to the magnetic field in (167) and evolved according to (168).

Tricco et al. (Reference Tricco, Price and Bate2016a) generalised the method of Dedner et al. (Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002) to include the case where the hyperbolic wave speed, c h, varies in time and space. This is the approach we use in Phantom. The resulting ‘generalised wave equation’ may be derived by combining the relevant term in (167) with (168) to give (Tricco et al. Reference Tricco, Price and Bate2016a)

(188) $$\begin{eqnarray} &&\frac{{\rm d}}{{\rm d}t} \left[ \frac{1}{\sqrt{\rho } c_{\rm h}}\frac{{\rm d}}{{\rm d}t} \left( \frac{\psi }{\sqrt{\rho } c_{\rm h}} \right) \right] - \frac{\nabla ^{2} \psi }{\rho } \nonumber \\ &&\quad +\, \frac{{\rm d}}{{\rm d}t} \left[ \frac{1}{\sqrt{\rho } c_{\rm h}} \left(\frac{\psi }{\sqrt{\rho }c_{\rm h}\tau _{\rm c}}\right) \right] = 0. \end{eqnarray}$$

When c h, ρ, τc, and the fluid velocity are constant, this reduces to the usual damped wave equation in the form

(189) $$\begin{equation} \frac{\partial ^{2} \psi }{\partial t^{2}} - c_{\rm h}^{2} \nabla ^{2} \psi + \frac{1}{\tau _{\rm c}} \frac{\partial \psi }{\partial t} = 0. \end{equation}$$

The same equation holds for the evolution of ${\nabla \cdot {\bm B}}$ itself, i.e.

(190) $$\begin{equation} \frac{\partial ^{2} (\nabla \cdot {\bm B})}{\partial t^{2}} - c_{\rm h}^{2} \nabla ^{2} (\nabla \cdot {\bm B}) + \frac{1}{\tau _{\rm c}} \frac{\partial (\nabla \cdot {\bm B})}{\partial t} = 0, \end{equation}$$

from which it is clear that c h represents the speed at which divergence errors are propagated and τc is the decay timescale over which divergence errors are removed.

Tricco & Price (Reference Tricco and Price2012) formulated a ‘constrained’ SPMHD implementation of hyperbolic/parabolic cleaning which guarantees numerical stability of the cleaning. The constraint imposed by Tricco & Price (Reference Tricco and Price2012) is that, in the absence of damping, any energy removed from the magnetic field during cleaning must be conserved by the scalar field, ψ. This enforces particular choices of numerical operators for $\nabla \cdot {\bm B}$ and ∇ψ in (172) and (173), respectively, in particular that they form a conjugate pair of difference and symmetric derivative operators. This guarantees that the change of magnetic energy is negative definite in the presence of the parabolic term (see below).

In Phantom, we set the cleaning speed, c h, equal to the fast magnetosonic speed [Equation (180)] so that its timestep constraint is satisfied already by (72), as recommended by Tricco & Price (Reference Tricco and Price2013). The decay timescale is set according to

(191) $$\begin{equation} \tau ^a_{\rm c} \equiv \frac{h_a}{\sigma _{\rm c} c_{{\rm h},a}}, \end{equation}$$

where the dimensionless factor σc sets the ratio of parabolic to hyperbolic cleaning. This is set to σc = 1.0 by default, which was empirically determined by Tricco & Price (Reference Tricco and Price2012) to provide optimal reduction of divergence errors in three dimensions.

The divergence cleaning dissipates energy from the magnetic field at a rate given by (Tricco & Price Reference Tricco and Price2012)

(192) $$\begin{equation} \left(\frac{{\rm d} E}{{\rm d}t}\right)_{\rm cleaning} = -\sum _{a} m_{a} \frac{\psi ^{2}_{a}}{\mu _{0} \rho _{a} c_{{\rm h},a}^{2} \tau ^{a}_{\rm c}}. \end{equation}$$

In general, this is so small compared to other dissipation terms (e.g. resistivity for shock capturing) that it is not worth monitoring (Tricco et al. Reference Tricco, Price and Bate2016a). This energy is not added as heat, but simply removed from the calculation.

2.10.9. Monitoring of divergence errors and over-cleaning

The divergence cleaning algorithm is guaranteed to either conserve or dissipate magnetic energy, and cleans divergence errors to a sufficient degree for most applications. However, the onus is on the user to ensure that divergence errors are not affecting simulation results. This may be monitored by the dimensionless quantity

(193) $$\begin{equation} \epsilon _{\rm divB} \equiv \frac{h\vert \nabla \cdot {\bm B}\vert }{\vert {\bm B} \vert }. \end{equation}$$

The maximum and mean values of this quantity should be used to check the fidelity of simulations that include magnetic fields. A good rule-of-thumb is that the mean should remain ≲ 10−2 for the simulation to remain qualitatively unaffected by divergence errors.

The cleaning wave speed can be arbitrarily increased to improve the effectiveness of the divergence cleaning according to

(194) $$\begin{equation} c_{{\rm h}, a} = f_{\rm clean} v_a , \end{equation}$$

where f clean is an ‘over-cleaning’ factor (by default, f clean = 1, i.e. no ‘over-cleaning’). Tricco et al. (Reference Tricco, Price and Bate2016a) showed that increasing f clean leads to further reduction in divergence errors, without affecting the quality of obtained results, but with an accompanying computational expense associated with a reduction in the timestep size.

2.11. Non-ideal magnetohydrodynamics

Phantom implements non-ideal MHD including terms for Ohmic resistivity, ambipolar (ion-neutral) diffusion and the Hall effect. Algorithms and tests are taken from Wurster, Price, & Ayliffe (Reference Wurster, Price and Ayliffe2014) and Wurster, Price, & Bate (Reference Wurster, Price and Bate2016). See Wurster et al. (Reference Wurster, Price and Bate2016, Reference Wurster, Price and Bate2017) and Wurster, Bate, & Price (Reference Wurster, Bate and Price2018) for recent applications. Our formulation of non-ideal SPMHD in Phantom is simpler than the earlier formulation proposed by Hosking & Whitworth (Reference Hosking and Whitworth2004) because we consider only one set of particles, representing a mixture of charged and uncharged species. Ours is similar to the implementation by Tsukamoto, Iwasaki, & Inutsuka (Reference Tsukamoto, Iwasaki and Inutsuka2013) and Tsukamoto et al. (Reference Tsukamoto, Iwasaki, Okuzumi, Machida and Inutsuka2015).

2.11.1. Equations of non-ideal MHD

We assume the strong coupling or ‘heavy ion’ approximation (see e.g. Wardle & Ng Reference Wardle and Ng1999; Shu et al. Reference Shu, Galli, Lizano and Cai2006; Pandey & Wardle Reference Pandey and Wardle2008), which neglects ion pressure and assumes ρi ≪ ρn, where the subscripts i and n refer to the ionised and neutral fluids, respectively. In this case, (167) contains three additional terms in the form

(195) $$\begin{equation} \frac{{\rm d}}{{\rm d}t} \left( \frac{{\bm B}}{\rho } \right)_{\rm nimhd} = -\frac{1}{\rho } \nabla \times \left[ \frac{\bm J}{\sigma _e} + \frac{{\bm J} \times {\bm B}}{e n_e} - \frac{({\bm J} \times {\bm B}) \times {\bm B}}{\gamma _{\rm AD} \rho _i}\right], \end{equation}$$

where σe is the electrical conductivity, n e is the number density of electrons, e is the charge on an electron, and γAD is the collisional coupling constant between ions and neutrals (Pandey & Wardle Reference Pandey and Wardle2008). We write this in the form

(196) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left( \frac{{\bm B}}{\rho } \right)_{\rm nimhd} &=& -\frac{1}{\rho } \nabla \times \bigg [ \eta _{\rm O} {\bm J} + \eta _{\rm H} {\bm J} \times \hat{\bm B} \nonumber \\ && -\, \left. \eta _{\rm AD} ({\bm J} \times \hat{\bm B}) \times \hat{\bm B} \right], \end{eqnarray}$$

where $\hat{\bm B}$ is the unit vector in the direction of ${\bm B}$ such that ηO, ηAD, and ηHall are the coefficients for resistive and ambipolar diffusion and the Hall effect, respectively, each with dimensions of area per unit time.

To conserve energy, we require the corresponding resistive and ambipolar heating terms in the thermal energy equation in the form

(197) $$\begin{equation} \left( \frac{{\rm d}u}{{\rm d}t}\right)_{\rm nimhd} = \frac{ \eta _{\rm O}}{\rho } J^2 + \frac{\eta _{\rm AD}}{\rho } \left[ J^2 - ({\bm J}\cdot {\hat{\bm B}})^2\right]. \end{equation}$$

The Hall term is non-dissipative, being dispersive rather than diffusive, so does not enter the energy equation.

We currently neglect the ‘Biermann battery’ term (Biermann Reference Biermann1950) proportional to ∇P e/(en e) in our non-ideal MHD implementation, both because it is thought to be negligible in the interstellar medium (Pandey & Wardle Reference Pandey and Wardle2008) and because numerical implementations can produce incorrect results (Graziani et al. Reference Graziani, Tzeferacos, Lee, Lamb, Weide, Fatenejad and Miller2015). This term is mainly important in generating seed magnetic fields for subsequent dynamo processes (e.g. Khomenko et al. Reference Khomenko, Vitas, Collados and de Vicente2017).

2.11.2. Discrete equations

Our main constraint is that the numerical implementation of the non-ideal MHD terms should exactly conserve energy, which is achieved by discretising in the form (Wurster et al. Reference Wurster, Price and Ayliffe2014)

(198) $$\begin{eqnarray} \frac{{\rm d}}{{\rm d}t} \left( \frac{{\bm B}}{\rho } \right)_{{\rm nimhd},a} &=& - \sum _b \Bigg [ \frac{{\bm D}_a}{\Omega _a \rho _a^2} \times \nabla _a W_{ab} (h_a) \nonumber \\ &&+\, \frac{{\bm D}_b}{\Omega _b \rho _b^2} \times \nabla _a W_{ab} (h_b) \Bigg ], \end{eqnarray}$$

where

(199) $$\begin{equation} {\bm D} = -\eta _{\rm O} {\bm J} - \eta _{\rm H} ({\bm J} \times \hat{\bm B}) + \eta _{\rm AD} [({\bm J} \times \hat{\bm B}) \times \hat{\bm B}]. \end{equation}$$

The corresponding term in the energy equation is given by

(200) $$\begin{equation} \left( \frac{{\rm d}u_a}{{\rm d}t}\right)_{\rm nimhd} = - \frac{{\bm D}_a \cdot {\bm J}_a}{\rho _a}, \end{equation}$$

where the magnetic current density is computed alongside the density evaluation according to

(201) $$\begin{equation} {\bm J} = \frac{1}{\Omega _a\rho _a} \sum _b m_b ({\bm B}_a - {\bm B}_b) \times \nabla _a W_{ab} (h_a). \end{equation}$$

Non-ideal MHD therefore utilises a ‘two first derivatives’ approach, similar to the formulation of physical viscosity described in Section 2.7.1. This differs from the ‘direct second derivatives’ approach used for our artificial resistivity term, and in previous formulations of physical resistivity (Bonafede et al. Reference Bonafede, Dolag, Stasyszyn, Murante and Borgani2011). In practice, the differences are usually minor. Our main reason for using two first derivatives for non-ideal MHD is that it is easier to incorporate the Hall effect and ambipolar diffusion terms.

2.11.3. Computing the non-ideal MHD coefficients

To self-consistently compute the coefficients ηO, ηH, and ηAD from the local gas properties, we use the nicil library (Wurster Reference Wurster2016) for cosmic ray ionisation chemistry and thermal ionisation. We refer the reader to Wurster (Reference Wurster2016) and Wurster et al. (Reference Wurster, Price and Bate2016) for details, since this is maintained and versioned as a separate package.

2.11.4. Timestepping

With explicit timesteps, the timestep is constrained in a similar manner to other terms, using

(202) $$\begin{equation} \Delta t = \frac{C_{\rm nimhd} h^2}{\max (\eta _{\rm O},\eta _{\rm AD}, \vert \eta _{\rm H}\vert )}, \end{equation}$$

where C nimhd = 1/(2π) by default. This can prove prohibitive, so we employ the so-called ‘super-timestepping’ algorithm from Alexiades, Amiez, & Gremaud (Reference Alexiades, Amiez and Gremaud1996) to relax the stability criterion for the Ohmic and ambipolar terms (only). The implementation is described in detail in Wurster et al. (Reference Wurster, Price and Bate2016). Currently, the Hall effect is always timestepped explicitly in the code.

2.12. Self-gravity

Phantom includes self-gravity between particles. By self-gravity, we mean a solution to Poisson’s equation

(203) $$\begin{equation} \nabla ^{2} \Phi = 4\pi G \rho ({\bm r}), \end{equation}$$

where Φ is the gravitational potential and ρ represents a continuous fluid density. The corresponding acceleration term in the equation of motion is

(204) $$\begin{equation} {\bm a}_{\rm selfgrav} = -\nabla \Phi . \end{equation}$$

Since (203) is an elliptic equation, implying instant action, it requires a global solution. This solution is obtained in Phantom by splitting the acceleration into ‘short-range’ and ‘long-range’ contributions.

(205) $$\begin{equation} {\bm a}^{a}_{\rm selfgrav} = {\bm a}^{a}_{\rm short} + {\bm a}^{a}_{\rm long} , \end{equation}$$

where the ‘short-range’ interaction is computed by direct summation over nearby particles, and the ‘long-range’ interaction is computed by hierarchical grouping of particles using the kd-tree.

The distance at which the gravitational acceleration is treated as ‘short-’ or ‘long-range’ is determined for each node–node pair, nm, either by the tree opening criterion,

(206) $$\begin{equation} \theta ^{2} < \left(\frac{s_{m}}{r_{nm}}\right)^{2}, \end{equation}$$

where 0 ⩽ θ ⩽ 1 is the tree opening parameter, or by nodes whose smoothing spheres intersect,

(207) $$\begin{equation} r_{nm}^{2} < \left[s_{n} + s_{m} + \max (R_{\rm kern} h_{\max }^{n}, R_{\rm kern} h_{\max }^{m})\right]^{2}. \end{equation}$$

Here, s is the node size, which is the minimum radius about the centre of mass that contains all the particles in the node, and h max is the maximum smoothing length of the particles within the node. Node pairs satisfying either of these criteria have the particles contained within them added to a trial neighbour list, used for computing the short-range gravitational acceleration. Setting θ = 0, therefore, leads to the gravitational acceleration computed entirely as ‘short-range’, that is, only via direction summation, while θ = 1 gives the fastest possible, but least accurate, gravitational force evaluation. The code default is θ = 0.5.

2.12.1. Short-range interaction

How to solve (203) on particles is one of the most widely misunderstood problems in astrophysics. In SPH or collisionless N-body simulations (i.e. stars and dark matter), the particles do not represent physical point-mass particles, but rather interpolation points in a density field that is assumed to be continuous. Hence, one needs to first reconstruct the density field on the particles, then solve (203) in a manner which is consistent with this (e.g. Hernquist & Barnes Reference Hernquist and Barnes1990).

How to do this consistently using a spatially adaptive softening length was demonstrated by Price & Monaghan (Reference Price and Monaghan2007), since an obvious choice is to use the iterative kernel summation in (3) to both reconstruct the density field and set the softening length, i.e.Footnote 6

(208) $$\begin{equation} \rho _{a} = \sum _{b} m_{b} W_{ab} (\epsilon _{a}); \hspace{14.22636pt} \epsilon _{a} = h_{\rm fac} (m_{a}/\rho _{a})^{1/3}. \end{equation}$$

It can then be shown that the gravitational potential consistent with this choice is given by

(209) $$\begin{equation} \Phi _{a} = -G \sum _{b} m_{b} \phi (\vert {\bm r}_{a} - {\bm r}_{b} \vert , \epsilon _{a}), \end{equation}$$

where ϕ is the softening kernel derived from the density kernel via Poisson’s equation (Section 2.12.2). For a variable smoothing length, energy conservation necessitates taking an average of the kernels, i.e.

(210) $$\begin{equation} \Phi _{a} = -G \sum _{b} m_{b} \left[ \frac{\phi _{ab}(\epsilon _{a}) + \phi _{ab}(\epsilon _{b})}{2} \right]. \end{equation}$$

Price & Monaghan (Reference Price and Monaghan2007) showed how the equations of motion could then be derived from a Lagrangian in order to take account of the softening length gradient terms, giving equations of motion in the form

(211) $$\begin{eqnarray} \hspace*{-15pt}{\bm a}^{a}_{\rm selfgrav} &= & -\nabla \Phi _{a}, \nonumber \\ \hspace*{-15pt}&=& -G \sum _{b} m_{b} \left[ \frac{\phi ^{\prime }_{ab}(\epsilon _{a}) + \phi ^{\prime }_{ab}(\epsilon _{b})}{2}\right] \hat{\bm r}_{ab} \nonumber \\ \hspace*{-15pt}&& -\,\frac{G}{2} \sum _{b} m_{b} \left[ \frac{\zeta _{a}}{\Omega ^{\epsilon }_{a}} \nabla _{a} W_{ab} (\epsilon _{a}) + \frac{\zeta _{b}}{\Omega ^{\epsilon }_{b}} \nabla _{a} W_{ab} (\epsilon _{b})\right], \end{eqnarray}$$

where Ωε and ζ are correction terms necessary for energy conservation, with Ω as in (5) but with h replaced by ε, and ζ defined as

(212) $$\begin{equation} \zeta _{a} \equiv \frac{\partial \epsilon _{a}}{\partial \rho _{a}} \sum _{b} m_{b} \frac{\partial \phi _{ab}(\epsilon _{a})}{\partial \epsilon _{a}}. \end{equation}$$

The above formulation satisfies all of the conservation properties, namely conservation of linear momentum, angular momentum, and energy.

The short-range acceleration is evaluated for each particle in the leaf node n by summing (211) over the trial neighbour list obtained by node–node pairs that satisfy either of the criteria in (206) or (207). For particle pairs separated outside the softening radius of either particle, the short range interaction reduces to

(213) $$\begin{equation} {\bm a}^{a}_{\rm short, r > R_{\rm kern}\epsilon } = -G \sum _{b} m_{b} \frac{ {\bm r}_{a} - {\bm r}_{b}}{\vert {\bm r}_{a} - {\bm r}_{b}\vert ^{3}}. \end{equation}$$

We use this directly for such pairs to avoid unnecessary evaluations of the softening kernel.

It is natural in SPH to set the gravitational softening length equal to the smoothing length ε = h, since both derive from the same density estimate. Indeed, Bate & Burkert (Reference Bate and Burkert1997) showed that this is a crucial requirement to resolve gas fragmentation correctly. For pure N-body calculations, Price & Monaghan (Reference Price and Monaghan2007) also showed that setting the (variable) softening length in the same manner as the SPH smoothing length (Sections 2.1.32.1.4) results in a softening that is always close to the ‘optimal’ fixed softening length (Merritt Reference Merritt1996; Athanassoula et al. Reference Athanassoula, Fady, Lambert and Bosma2000; Dehnen Reference Dehnen2001). In collisionless N-body simulations, this has been found to increase the resolving power, giving results at lower resolution comparable to those obtained at higher resolution with a fixed softening length (Bagla & Khandai Reference Bagla and Khandai2009; Iannuzzi & Dolag Reference Iannuzzi and Dolag2011). It also avoids the problem of how to ‘choose’ the softening length, negating the need for ‘rules of thumb’ such as the one given by Springel (Reference Springel2005) where the softening length is chosen to be 1/40 of the average particle spacing in the initial conditions.

2.12.2. Functional form of the softening kernel

The density kernel and softening potential kernel are related using Poisson’s equation (203), i.e.

(214) $$\begin{equation} W(r, \epsilon ) = \frac{1}{4\pi r^{2}} \frac{\partial }{\partial r} \left( r^{2} \frac{\partial \phi }{\partial r} \right), \end{equation}$$

where $r \equiv \vert {\bm r} - {\bm r}^{\prime }\vert$. Integrating this equation gives the softening kernels used for the force and gravitational potential. As with the standard kernel functions (Section 2.1.5), we define the potential and force kernels in terms of dimensionless functions of the scaled interparticle separation, qr/h, according to

(215) $$\begin{eqnarray} \phi (r, \epsilon ) \equiv \frac{1}{\epsilon } \tilde{\phi } (q), \end{eqnarray}$$
(216) $$\begin{eqnarray} \phi ^{\prime }(r, \epsilon ) \equiv \frac{1}{\epsilon ^{2}} \tilde{\phi }^{\prime } (q), \end{eqnarray}$$

where the dimensionless force kernel is obtained from the density kernel f(q) (Section 2.1.52.1.6) using

(217) $$\begin{equation} \tilde{\phi }^{\prime } (q) = \frac{4\pi }{q^{2}} C_{\rm norm} \int f(q^{\prime }) q^{\prime 2} {\rm d}q^{\prime }, \end{equation}$$

with the integration constant set such that $\tilde{\phi }^{\prime }(q) = 1/q^{2}$ at q = R kern. The potential function is

(218) $$\begin{equation} \tilde{\phi }(q) = \int \tilde{\phi }^{\prime }(q^{\prime }) dq^{\prime }. \end{equation}$$

The derivative of ϕ with respect to ε required in (212) is also written in terms of a dimensionless function, i.e.

(219) $$\begin{equation} \frac{\partial \phi (r,\epsilon )}{\partial \epsilon } \equiv \frac{1}{\epsilon ^{2}} h(q), \end{equation}$$

where from differentiating (215), we have

(220) $$\begin{equation} h(q) = -\tilde{\phi }(q) - q\tilde{\phi }^{\prime }(q). \end{equation}$$

Included in the Phantom source code is a Python script using the sympy library to solve (217), (218), and (220) using symbolic integration to obtain the softening kernel from the density kernel. This makes it straightforward to obtain the otherwise laborious expressions needed for each kernel [expressions for the cubic spline kernel are in Appendix A of Price & Monaghan (Reference Price and Monaghan2007) and for the quintic spline in Appendix A of Hubber et al. (Reference Hubber, Batty, McLeod and Whitworth2011)]. Figure 3 shows the functional form of the softening kernels for each of the kernels available in Phantom. The kernel function f(q) is shown for comparison (top row in each case).

Figure 3. Functional form of the softening kernel functions −ϕ(r, h) and ϕ′(r, h) used to compute the gravitational force in Phantom, shown for each of the available kernel functions w(r, h) (see Figure 1). Dotted lines show the functional form of the unsoftened potential (−1/r) and force (1/r 2) for comparison.

2.12.3. Softening of gravitational potential due to stars, dark matter, and dust

In the presence of collisionless components (e.g. stars, dark matter, and dust), we require estimates of the density in order to set the softening lengths for each component. We follow the generalisation of the SPH density estimate to multi-fluid systems described by Laibe & Price (Reference Laibe and Price2012a) where the density (and hence the softening length and Ωε) for each component is computed in the usual iterative manner (Section 2.1.4), but using only neighbours of the same type (c.f. Section 2.13.3). That is, the softening length for stars is determined based on the local density of stars, the softening length for dark matter is based on the local density of dark matter particles, and so on. The gravitational interaction both within each type and between different types is then computed using (211). This naturally symmetrises the softening between different types, ensuring that momentum, angular momentum, and energy are conserved.

2.12.4. Long-range interaction

At long range, that is r > R kernεa and r > R kernεb, the second term in (211) tends to zero since ζ = 0 for qR kern, while the first term simply becomes 1/r 2. Computing this via direct summation would have an associated $\mathcal {O}(N^{2})$ computational cost, thus we adopt the usual practice of using the kd-tree to reduce this cost to $\mathcal {O}(N\log N)$.

The main optimisation in Phantom compared to a standard tree code (e.g. Hernquist & Katz Reference Hernquist and Katz1989; Benz et al. Reference Benz, Cameron, Press and Bowers1990) is that we compute the long-range gravitational interaction once per leaf-node rather than once per-particle and that we use Cartesian rather than spherical multipole expansions to enable this (Dehnen Reference Dehnen2000b; Gafton & Rosswog Reference Gafton and Rosswog2011).

The long-range acceleration on a given leaf node n consists of a sum over distant nodes m that satisfy neither (206) nor (207).

(221) $$\begin{equation} {\bm a}_{n} = \sum _{m} {\bm a}_{nm}. \end{equation}$$

The acceleration along the node centres, between a given pair n and m, is given (using index notation) by

(222) $$\begin{equation} a^{i}_{nm} = -\frac{GM_{m}}{r^{3}} r^{i} + \frac{1}{r^{4}} \left(\hat{r}^{k} Q^{m}_{ik} - \frac{5}{2} \hat{r}^{i} \mathcal {Q}^{m} \right), \end{equation}$$

where rixinxim is the relative position vector, $\hat{r}$ is the corresponding unit vector, M m is the total mass in node m, Qmij is the quadrupole moment of node m, and repeated indices imply summation. We define the following scalar and vector quantities for convenience:

(223) $$\begin{eqnarray} \mathcal {Q} \equiv \hat{r}^{i} \hat{r}^{j} Q_{ij},\vspace*{-15pt} \end{eqnarray}$$
(224) $$\begin{eqnarray} \mathcal {Q}_{i} \equiv \hat{r}^{j} Q_{ij}. \end{eqnarray}$$

Alongside the acceleration, we compute six independent components of the first derivative matrix:

(225) $$\begin{equation} \frac{\partial a^{i}_{n}}{\partial r^{j}} = \sum _{m} \frac{\partial a^{i}_{nm}}{\partial r^{j}}, \end{equation}$$

where

(226) $$\begin{eqnarray} \frac{\partial a_{nm}^{i}}{\partial r^{j}} &= & \frac{GM_{m}}{r^{3}} \bigg [ 3\hat{r}^{i} \hat{r}^{j} - \delta ^{ij} \bigg ] \nonumber \\ & &+ \,\frac{1}{r^{5}} \left[ Q_{ij}^{m} + \left(\frac{35}{2} \hat{r}^{i} \hat{r}^{j} - \frac{5}{2} \delta _{ij}\right) \mathcal {Q}^{m} \right. \nonumber \\ && -\, 5 \hat{r}^{i} \mathcal {Q}_{j}^{m} - 5 \hat{r}^{j} \mathcal {Q}_{i}^{m} \bigg ], \end{eqnarray}$$

and ten independent components of the second derivatives, given by

(227) $$\begin{eqnarray} \frac{\partial ^{2} a^{i}_{nm}}{\partial r^{j} \partial r^{k}} &= & -\frac{3GM_{m}}{r^{4}} \bigg [ 5 \hat{r}^{i} \hat{r}^{j} \hat{r}^{k} - \delta _{jk} \hat{r}^{i} - \delta _{ik} \hat{r}^{j} - \delta _{ij} \hat{r}^{k} \bigg ] \nonumber \\ && +\, \frac{1}{r^{6}} \bigg [ -5(\hat{r}^{k} Q^{m}_{ij} + \hat{r}^{i} Q^{m}_{jk} + \hat{r}^{j} Q^{m}_{ik}) \nonumber \\ && -\, \frac{315}{2}\hat{r}^{i}\hat{r}^{j}\hat{r}^{k} \mathcal {Q}^{m} \nonumber \\ && + \,\frac{35}{2} \left(\delta _{ij} \hat{r}^{k} + \delta _{ik} \hat{r}^{j} + \delta _{jk} \hat{r}^{i} \right) \mathcal {Q}^{m} \nonumber \\ && +\, 35 \left( \hat{r}^{j} \hat{r}^{k} \mathcal {Q}^{m}_{i} + \hat{r}^{i} \hat{r}^{k} \mathcal {Q}^{m}_{j} + \hat{r}^{i} \hat{r}^{j} \mathcal {Q}^{m}_{k} \right) \nonumber \\ && -\, 5 (\delta _{ij} \mathcal {Q}^{m}_{k} + \delta _{ik} \mathcal {Q}^{m}_{j} + \delta _{jk} \mathcal {Q}^{m}_{i}) \bigg ]. \end{eqnarray}$$

The acceleration on each individual particle inside the leaf node n is then computed using a second-order Taylor series expansion of ${\bm a}_{\rm node}^{n}$ about the node centre, i.e.

(228) $$\begin{equation} a^{i}_{{\rm long}, a} = a^{i}_{n} + \Delta x^{j} \frac{\partial a_{n}^{i}}{\partial r^{j}} + \frac{1}{2} \Delta x^{j} \Delta x^{k} \frac{\partial ^{2} a_{n}^{i}}{\partial r^{j}\partial r^{k}}, \end{equation}$$

where Δxiaxiaxi 0 is the relative distance of each particle from the node centre of mass. Pseudo-code for the resulting force evaluation is shown in Figure A5.

The quadrupole moments are computed during the tree build using

(229) $$\begin{equation} Q_{ij} = \sum _{a} m_{a} \left[ 3 \Delta x^{i} \Delta x^{j} - (\Delta x)^{2} \delta ^{ij}\right], \end{equation}$$

where the sum is over all particles in the node. Since Q is a symmetric tensor, only six independent quantities need to be stored (Q xx, Q xy, Q xz, Q yy, Q yz and Q zz).

The current implementation in Phantom is $\mathcal {O}(N_{\rm {leafnodes}} \log N_{\rm part})$ rather than the $\mathcal {O}(N)$ treecode implementation proposed by Dehnen (Reference Dehnen2000b), since we do not currently implement the symmetric node–node interactions required for $\mathcal {O}(N)$ scaling. Neither does our treecode conserve linear momentum to machine precision, except when θ = 0. Implementing these additional features would be desirable.

2.13. Dust–gas mixtures

Modelling dust–gas mixtures is the first step in the ‘grand challenge’ of protoplanetary disc modelling (Haworth et al. Reference Haworth2016). The public version of Phantom implements dust–gas mixtures using two approaches. One models the dust and gas as two separate types of particles (two-fluid), as presented in Laibe & Price (Reference Laibe and Price2012a, Reference Laibe and Price2012b), and the other, for small grains, using a single type of particle that represents the combination of dust and gas together (one-fluid), as described in Price & Laibe (Reference Price and Laibe2015a). Various combinations of these algorithms have been used in our recent papers using Phantom, including Dipierro et al. (Reference Dipierro, Price, Laibe, Hirsh and Cerioli2015, Reference Dipierro, Laibe, Price and Lodato2016), Ragusa et al. (Reference Ragusa, Dipierro, Lodato, Laibe and Price2017), and Tricco et al. (Reference Tricco, Price and Laibe2017) (see also Hutchison et al. Reference Hutchison, Price, Laibe and Maddison2016).

In the two-fluid implementation, the dust and gas are treated as two separate fluids coupled by a drag term with only explicit timestepping. In the one-fluid implementation, the dust is treated as part of the mixture, with an evolution equation for the dust fraction.

2.13.1. Continuum equations

The two-fluid formulation is based on the continuum equations in the form:

$$\begin{eqnarray*} \hspace*{30pt}\frac{\partial \rho _{\rm g}}{\partial t} + ({\bm v}_{\rm g}\cdot \nabla ) \rho _{\rm g} &=& -\rho _{\rm g} (\nabla \cdot {\bm v}_{\rm g}),\hspace*{64pt}(230)\\[10pt] \hspace*{30pt}\frac{\partial \rho _{\rm d}}{\partial t} + ({\bm v}_{\rm d}\cdot \nabla ) \rho _{\rm d} &=& -\rho _{\rm d} (\nabla \cdot {\bm v}_{\rm d}),\hspace*{64pt}(231)\\[10pt] \hspace*{20pt}\frac{\partial {\bm v}_{\rm g}}{\partial t} + ({\bm v}_{\rm g}\cdot \nabla ) {\bm v}_{\rm g} &=& -\frac{\nabla P}{\rho _{\rm g}} + \frac{K}{\rho _{\rm g}} ({\bm v}_{\rm d}- {\bm v}_{\rm g}),\hspace*{31pt}(232)\\[10pt] \hspace*{30pt}\frac{\partial {\bm v}_{\rm d}}{\partial t} + ({\bm v}_{\rm d}\cdot \nabla ) {\bm v}_{\rm d} &=& - \frac{K}{\rho _{\rm d}} ({\bm v}_{\rm d}- {\bm v}_{\rm g}),\hspace*{55pt}(233)\\[10pt] \hspace*{20pt}\frac{\partial u}{\partial t} + ({\bm v}_{\rm g}\cdot \nabla ) u &=& -\frac{P}{\rho _{\rm g}} (\nabla \cdot {\bm v}_{\rm g}) + \frac{\Lambda_{\rm drag}}{\rho_{\rm g}} ,\hspace*{32pt}(234) \end{eqnarray*}$$

where the subscripts g and d refer to gas and dust properties, K is a drag coefficient and the drag heating is

(235) $$\begin{equation} \Lambda _{\rm drag} \equiv K ({\bm v}_{\rm d}- {\bm v}_{\rm g})^{2}. \end{equation}$$

The implementation in Phantom currently neglects any thermal coupling between the dust and the gas (see Laibe & Price Reference Laibe and Price2012a), aside from the drag heating. Thermal effects are important for smaller grains since they control the heating and cooling of the gas (e.g. Dopcke et al. Reference Dopcke, Glover, Clark and Klessen2011). Currently, the internal energy (u) of dust particles is simply set to zero.

2.13.2. Stopping time

The stopping time

(236) $$\begin{equation} t_{\rm s} = \frac{\rho _{\rm g} \rho _{\rm d}}{K(\rho _{\rm g} + \rho _{\rm d})} \end{equation}$$

is the characteristic timescale for the exponential decay of the differential velocity between the two phases caused by the drag. In the code, t s is by default specified in physical units, which means that code units need to be defined appropriately when simulating dust–gas mixtures.

2.13.3. Two-fluid dust–gas in SPH

In the two-fluid method, the mixture is discretised into two distinct sets of ‘dust’ and ‘gas’ particles. In the following, we adopt the convention from Monaghan & Kocharyan (Reference Monaghan and Kocharyan1995) that subscripts a, b, and c refer to gas particles while i, j, and k refer to dust particles. Hence, (230)–(231) are discretised with a density summation over neighbours of the same type (c.f. Section 2.12.3), giving

(237) $$\begin{equation} \rho _{a} = \sum _{b} m_{b} W_{ab} (h_{a}); \hspace{14.22636pt} h_{a} = h_{\rm fact} \left( \frac{m_{a}}{\rho _{a}} \right) ^{1/3}, \end{equation}$$

for a gas particle, and

(238) $$\begin{equation} \rho _{i} = \sum _{j} m_{j} W_{ij} (h_{i}); \hspace{14.22636pt} h_{i} = h_{\rm fact} \left( \frac{m_{i}}{\rho _{i}} \right) ^{1/3}, \end{equation}$$

for a dust particle. The kernel used for density is the same as usual (Section 2.1.6). We discretise the equations of motion for the gas particles, (232), using

(239) $$\begin{equation} \left(\frac{{\rm d}{\bm v}_{a}}{{\rm d}t}\right)_{\rm drag} = -3 \sum _{j} m_{j} \frac{{\bm v}_{aj} \cdot \hat{\bm r}_{aj}}{(\rho _{a} + \rho _{j}) t^{\rm s}_{aj}} \hat{\bm r}_{aj} D_{aj} (h_{a}), \end{equation}$$

and for dust, (233), using

(240) $$\begin{equation} \left(\frac{{\rm d}{\bm v}_{i}}{{\rm d}t}\right)_{\rm drag} = -3 \sum _{b} m_{b} \frac{{\bm v}_{ib} \cdot \hat{\bm r}_{ib}}{(\rho _{i} + \rho _{b}) t^{\rm s}_{ib}} \hat{\bm r}_{ib} D_{ib} (h_{b}), \end{equation}$$

where D is a ‘double hump’ kernel, defined in Section 2.13.4. The drag heating term in the energy equation, (234), is discretised using

(241) $$\begin{equation} \left( \frac{{\rm d}u}{{\rm d}t} \right)_{\rm drag} = \frac{\Lambda_{\rm drag}}{\rho_{\rm g}} = 3 \sum _{j} m_{j} \frac{({\bm v}_{aj} \cdot \hat{\bm r}_{aj})^{2}}{(\rho _{a} + \rho _{j}) t^{\rm s}_{aj}} D_{aj} (h_{a}). \end{equation}$$

Notice that gas properties are only defined on gas particles and dust properties are defined only on dust particles, greatly simplifying the algorithm. Buoyancy terms caused by dust particles occupying a finite volume (Monaghan & Kocharyan Reference Monaghan and Kocharyan1995; Laibe & Price Reference Laibe and Price2012a) are negligible in astrophysics because typical grain sizes (μm) are negligible compared to simulation scales of ~au or larger.

2.13.4. Drag kernels

Importantly, we use a ‘double-hump’ shaped kernel function D (Fulk & Quinn Reference Fulk and Quinn1996) instead of the usual bell-shaped kernel W when computing the drag terms. Defining D in terms of a dimensionless kernel function as previously (c.f. Section 2.1.5).

(242) $$\begin{equation} D(r,h) = \frac{\sigma }{h^{3}} g(q), \end{equation}$$

then the double hump kernels are defined from the usual kernels according to

(243) $$\begin{equation} g(q) = q^{2} f(q), \end{equation}$$

where the normalisation constant σ is computed by enforcing the usual normalisation condition

(244) $$\begin{equation} \int D (r, h) {\rm d}V = 1. \end{equation}$$

Figure 4 shows the functional forms of the double hump kernels used in Phantom. Using double hump kernels for the drag terms was found by Laibe & Price (Reference Laibe and Price2012a) to give a factor of 10 better accuracy at no additional cost. The key feature is that these kernels are zero at the origin putting more weight in the outer parts of the kernel where the velocity difference is large. This also solves the problem of how to define the unit vector in the drag terms (239)–(241)—it does not matter since D is also zero at the origin.

Figure 4. Double hump smoothing kernels D(r, h) available in Phantom, used in the computation of the dust–gas drag force.

2.13.5. Stopping time in SPH

The stopping time is defined between a pair of particles, using the properties of gas and dust defined on the particle of the respective type, i.e.

(245) $$\begin{equation} t^{\rm s}_{aj} = \frac{\rho _{a} \rho _{j}}{K_{aj} (\rho _{a} + \rho _{j})}. \end{equation}$$

The default prescription for the stopping time in Phantom automatically selects a physical drag regime appropriate to the grain size, as described below and in Laibe & Price (Reference Laibe and Price2012b). Options to use either a constant K or a constant t s between pairs are also implemented, useful for testing and debugging (c.f. Section 5.9).

2.13.6. Epstein drag

To determine the appropriate physical drag regime, we use the procedure suggested by Stepinski & Valageas (Reference Stepinski and Valageas1996) where we first evaluate the Knudsen number

(246) $$\begin{equation} K_{\rm n} = \frac{9\lambda _{\rm g}}{4 s_{\rm grain}}, \end{equation}$$

where s grain is the grain size and λg is the gas mean free path (see Section 2.13.8, below, for how this is evaluated). For K n ⩾ 1, the drag between a particle pair is computed using the generalised formula for Epstein drag from Kwok (Reference Kwok1975), as described in Paardekooper & Mellema (Reference Paardekooper and Mellema2006) and Laibe & Price (Reference Laibe and Price2012b), giving

(247) $$\begin{equation} K_{aj} = \rho ^{a}_{\rm g} \rho ^{j}_{\rm d} \frac{4}{3} \sqrt{\frac{8\pi }{\gamma }} \frac{s_{\rm grain}^{2}}{m_{\rm grain}} c^{a}_{\rm s} f, \end{equation}$$

where

(248) $$\begin{equation} m_{\rm grain} \equiv \frac{4}{3} \pi \rho _{\rm grain} s_{\rm grain}^{3}, \end{equation}$$

and ρgrain is the intrinsic grain density, which is 3 g/cm3 by default. The variable f is a correction for supersonic drift velocities given by (Kwok Reference Kwok1975)

(249) $$\begin{equation} f \equiv \sqrt{1 + \frac{9\pi }{128} \frac{\Delta v^{2}}{c_{\rm s}^{2}}}, \end{equation}$$

where $\Delta v \equiv \vert {\bm v}_{\rm d}- {\bm v}_{\rm g}\vert = v_{\rm d}^{j} - v_{\rm g}^{a}$. The stopping time is therefore

(250) $$\begin{equation} t_{\rm s} = \frac{ \rho _{\rm grain} s_{\rm grain}}{\rho c_{\rm s}f} \sqrt{\frac{\pi \gamma }{8}}, \end{equation}$$

where ρ ≡ ρd + ρg. This formula, (250), reduces to the standard expression for the linear Epstein regime in the limit where the drift velocity is small compared to the sound speed (i.e. f → 1). Figure 5 shows the difference between the above simplified prescription and the exact expression for Epstein drag (Epstein Reference Epstein1924; c.f. equations (11) and (38) in Laibe & Price Reference Laibe and Price2012b) as a function of Δv/c s, which is less than 1% everywhere.

Figure 5. Dependence of the drag stopping time t s on differential Mach number, showing the increased drag (decrease in stopping time) as the velocity difference between dust and gas increases. The black line shows the analytic approximation we employ [Equation (250)] which may be compared to the red line showing the exact expression from Epstein (Reference Epstein1924). The difference is less than 1% everywhere.

2.13.7. Stokes drag

For K n < 1, we adopt a Stokes drag prescription, describing the drag on a sphere with size larger than the mean free path of the gas (Fan & Zhu Reference Fan and Zhu1998). Here, we use (Laibe & Price Reference Laibe and Price2012b)

(251) $$\begin{equation} K_{aj} = \rho ^{a}_{\rm g} \rho ^{j}_{\rm d} \frac{1}{2} C_{\rm D} \frac{\pi s_{\rm grain}^{2}}{m_{\rm grain}} \vert \Delta v\vert , \end{equation}$$

where the coefficient C D is given by (Fassio & Probstein Reference Fassio and Probstein1970) (see Whipple Reference Whipple and Elvius1972; Weidenschilling Reference Weidenschilling1977)

(252) $$\begin{equation} C_{\rm D} = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}24 R_{e}^{-1}, & R_{\rm e} < 1, \\ 24 R_{e}^{-0.6}, & 1 < R_{\rm e} < 800, \\ 0.44, & R_{\rm e} > 800, \end{array}\right. \end{equation}$$

where R e is the local Reynolds number around the grain

(253) $$\begin{equation} R_{\rm e} \equiv \frac{2 s_{\rm grain} \vert \Delta v \vert }{\nu }, \end{equation}$$

and ν is the microscopic viscosity of the gas (see below, not to be confused with the disc viscosity). Similar formulations of Stokes drag can be found elsewhere (see e.g. discussion in Woitke & Helling Reference Woitke and Helling2003 and references therein). The stopping time in the Stokes regime is therefore given by

(254) $$\begin{equation} t_{\rm s} = \frac{8 \rho _{\rm grain} s_{\rm grain}}{3\rho \vert \Delta v\vert C_{\rm D}}, \end{equation}$$

where it remains to evaluate ν and λg.

2.13.8. Kinematic viscosity and mean free path

We evaluate the microscopic kinematic viscosity, ν, assuming gas molecules interact as hard spheres, following Chapman & Cowling (Reference Chapman and Cowling1970). The viscosity is computed from the mean free path and sound speed according to

(255) $$\begin{equation} \nu = \sqrt{\frac{2}{\pi \gamma }} c_{\rm s} \lambda _{\rm g}, \end{equation}$$

with the mean free path defined by relating this expression to the expression for the dynamic viscosity of the gas (Chapman & Cowling Reference Chapman and Cowling1970) given by

(256) $$\begin{equation} \mu _{\nu } = \frac{5m}{64\sigma _{\rm s}} \sqrt{\frac{\pi }{\gamma }} c_{\rm s}, \end{equation}$$

with μν = ρgν, giving

(257) $$\begin{equation} \lambda _{\rm g} = \frac{5\pi }{64 \sqrt{2}} \frac{1}{n_{\rm g} \sigma _{s}}, \end{equation}$$

where n g = ρg/m is the number density of molecules and σs is the collisional cross section. To compute this, Phantom currently assumes the gas is molecular Hydrogen, such that the mass of each molecule and the collisional cross section are given by

(258) $$\begin{equation} m = 2 m_{\rm H}, \\ \end{equation}$$
(259) $$\begin{equation} \sigma _{\rm s} = 2.367 \times 10^{-15} {\rm cm}^{2}. \end{equation}$$

2.13.9. Stokes/Epstein transition

At the transition between the two regimes, assuming R e < 1, (254) reduces to

(260) $$\begin{equation} t_{\rm s} = \frac{2 \rho _{\rm grain} s_{\rm grain}^{2}}{9 \rho c_{\rm s} \lambda _{\rm g}} \sqrt{\frac{\pi \gamma }{2}}, \end{equation}$$

which is the same as the Epstein drag in the subsonic regime when λg = 4s grain/9, i.e. K n = 1. That this transition is indeed continuous in the code is demonstrated in Figure 6, which shows the transition from Epstein to Stokes drag and also through each of the Stokes regimes in (252) by varying the grain size while keeping the other parameters fixed. For simplicity, we assumed a fixed Δv = 0.01c s in this plot, even though in general one would expect Δv to increase with stopping time [see Equation (270)].

Figure 6. Drag stopping time t s (in years) as a function of grain size, showing the continuous transition between the Epstein and Stokes drag regimes. The example shown assumes fixed density ρ = 10−13g/cm3 and sound speed c s = 6 × 104 cm/s with subsonic drag Δv = 0.01c s and material density ρgrain = 3g/cm3.

2.13.10. Self-gravity of dust

With self-gravity turned on, dust particles interact in the same way as stars or dark matter (Section 2.12.3), with a softening length equal to the smoothing length determined from the density of neighbouring dust particles. Dust particles can be accreted by sink particles (Section 2.8.2), but a sink cannot currently be created from dust particles (Section 2.8.4). There is currently no mechanism in the code to handle the collapse of dust to form a self-gravitating object independent of the gas.

2.13.11. Timestep constraint

For the two-fluid method, the timestep is constrained by the stopping time according to

(261) $$\begin{equation} \Delta t^{a}_{\rm drag} = \min _{j} ( t_{\rm s}^{aj}). \end{equation}$$

This requirement, alongside the spatial resolution requirement hc st s (Laibe & Price Reference Laibe and Price2012a), means the two-fluid method becomes both prohibitively slow and increasingly inaccurate for small grains. In this regime, one should switch to the one-fluid method, as described below.

2.13.12. Continuum equations: One-fluid

In Laibe & Price (Reference Laibe and Price2014a), we showed that the two-fluid equations, (230)–(234), can be rewritten as a single fluid mixture using a change of variables given by

(262) $$\begin{eqnarray} \rho \equiv \rho _{\rm g} + \rho _{\rm d}, \end{eqnarray}$$
(263) $$\begin{eqnarray} \epsilon \equiv \rho _{\rm d}/\rho , \end{eqnarray}$$
(264) $$\begin{eqnarray} {\bm v} \equiv \frac{\rho _{\rm g} {\bm v}_{\rm g}+ \rho _{\rm d} {\bm v}_{\rm d}}{ \rho _{\rm g} + \rho _{\rm d} }, \end{eqnarray}$$
(265) $$\begin{eqnarray} \Delta {\bm v} \equiv {\bm v}_{\rm d}- {\bm v}_{\rm g}, \end{eqnarray}$$

where ρ is the combined density, ε is the dust fraction, ${\bm v}$ is the barycentric velocity, and $\Delta {\bm v}$ is the differential velocity between the dust and gas.

In Laibe & Price (Reference Laibe and Price2014a), we derived the full set of evolution equations in these variables, and in Laibe & Price (Reference Laibe and Price2014c), implemented and tested an algorithm to solve these equations in SPH. However, using a fluid approximation cannot properly capture the velocity dispersion of large grains, as occurs for example when large planetesimals stream simultaneously in both directions through the midplane of a protoplanetary disc. For this reason, the one-fluid equations are better suited to treating small grains, where the stopping time is shorter than the computational timestep. In this limit, we employ the ‘terminal velocity approximation’ (e.g. Youdin & Goodman Reference Youdin and Goodman2005) and the evolution equations reduce to (Laibe & Price Reference Laibe and Price2014a; Price & Laibe Reference Price and Laibe2015a)

(266) $$\begin{eqnarray} \frac{{\rm d}\rho }{{\rm d}t} = -\rho (\nabla \cdot {\bm v}), \end{eqnarray}$$
(267) $$\begin{eqnarray} \frac{{\rm d}{\bm v}}{{\rm d}t} = -\frac{\nabla P}{\rho } + {\bm a}_{\rm ext}, \end{eqnarray}$$
(268) $$\begin{eqnarray} \frac{{\rm d}\epsilon }{{\rm d}t} = -\frac{1}{\rho } \nabla \cdot \left[{\epsilon (1 - \epsilon ) \rho \Delta {\bm v}}\right], \end{eqnarray}$$
(269) $$\begin{eqnarray} \frac{{\rm d}u}{{\rm d}t} = -\frac{P}{\rho } (\nabla \cdot {\bm v}) + \epsilon (\Delta {\bm v} \cdot \nabla ) u, \end{eqnarray}$$

where

(270) $$\begin{equation} \Delta {\bm v} \equiv t_{\rm s} \left({\bm a}_{\rm d} - {\bm a}_{\rm g}\right), \end{equation}$$

where ${\bm a}_{\rm d}$ and ${\bm a}_{\rm g}$ refers to any acceleration acting only on the dust or gas phase, respectively. For the simple case of pure hydrodynamics, the only difference is the pressure gradient, giving

(271) $$\begin{equation} \Delta {\bm v} \equiv t_{\rm s} \frac{\nabla P}{\rho _{\rm g}} = \frac{t_{\rm s}}{(1 - \epsilon )} \frac{\nabla P}{\rho }, \end{equation}$$

such that (268) becomes

(272) $$\begin{equation} \frac{{\rm d}\epsilon }{{\rm d}t} = -\frac{1}{\rho } \nabla \cdot \left({\epsilon t_{\rm s}\nabla P}\right). \end{equation}$$

Importantly, the one-fluid dust algorithm does not result in any heating term in du/dt due to drag, because this term is $\mathcal {O}(\Delta {\bm v}^{2})$ and thus negligible (Laibe & Price Reference Laibe and Price2014a).

2.13.13. Visualisation of one-fluid results

Finally, when visualising results of one-fluid calculations, one must reconstruct the properties of the dust and gas in post-processing. We use

(273) $$\begin{eqnarray} \rho _{\rm g} = (1 - \epsilon ) \rho , \end{eqnarray}$$
(274) $$\begin{eqnarray} \rho _{\rm d} = \epsilon \rho , \end{eqnarray}$$
(275) $$\begin{eqnarray} {\bm v}_{\rm g} = {\bm v} - \epsilon \Delta {\bm v}, \end{eqnarray}$$
(276) $$\begin{eqnarray} {\bm v}_{\rm d} = {\bm v} + (1 - \epsilon ) \Delta {\bm v}. \end{eqnarray}$$

To visualise the one-fluid results in a similar manner to those from the two-fluid method, we reconstruct a set of ‘dust’ and ‘gas’ particles with the same positions but with the respective properties of each type of fluid particle copied onto them. We compute $\Delta {\bm v}$ from (271) using the pressure gradient computed using (34), multiplied by the stopping time and the gas fraction. See Price & Laibe (Reference Price and Laibe2015a) for more details.

2.13.14. One-fluid dust–gas implementation

Our implementation of the one-fluid method in Phantom follows Price & Laibe (Reference Price and Laibe2015a) with a few minor changes and correctionsFootnote 7. In particular, we use the variable $s = \sqrt{\epsilon \rho }$ described in Appendix B of Price & Laibe (Reference Price and Laibe2015a) to avoid problems with negative dust fractions. The evolution equation (268) expressed in terms of s is given by

(277) $$\begin{equation} \frac{{\rm d}s}{{\rm d}t} = -\frac{1}{2s} \nabla \cdot \left( \frac{\rho _{\rm g} \rho _{\rm d}}{\rho } \Delta {\bm v} \right) - \frac{s}{2} (\nabla \cdot {\bm v}), \end{equation}$$

which for the case of hydrodynamics becomes

(278) $$\begin{eqnarray} \frac{{\rm d}s}{{\rm d}t} & =& -\frac{1}{2s} \nabla \cdot \left( \frac{s^{2}}{\rho } t_{\rm s} \nabla P \right) - \frac{s}{2} (\nabla \cdot {\bm v}), \nonumber \\ & =& -\frac{1}{2} \nabla \cdot \left( \frac{s}{\rho } t_{\rm s} \nabla P \right) - \frac{t_{\rm s}}{2\rho } \nabla P \cdot \nabla s - \frac{s}{2} (\nabla \cdot {\bm v}). \end{eqnarray}$$

The SPH discretisation of this equation is implemented in the form

(279) $$\begin{eqnarray} \frac{{\rm d}s_{a}}{{\rm d}t} &= & - \frac{1}{2} \sum _{b} \frac{m_{b} s_{b}}{\rho _{b}} \left(\frac{t_{{\rm s}, a}}{\rho _{a}} + \frac{t_{{\rm s}, b}}{\rho _b} \right) (P_{a} - P_{b}) \frac{\overline{F}_{ab}}{\vert r_{ab} \vert } \nonumber \\ && +\, \frac{s_{a}}{2\rho _{a} \Omega _{a}} \sum _{b} m_{b} {\bm v}_{ab} \cdot \nabla _{a} W_{ab} (h_{a}), \end{eqnarray}$$

where $\overline{F}_{ab} \equiv \frac{1}{2} [ F_{ab} (h_{a}) + F_{ab}(h_{b}) ]$. The thermal energy equation, (269), takes the form

(280) $$\begin{equation} \frac{{\rm d}u}{{\rm d}t} = - \frac{P}{\rho } (\nabla \cdot {\bm v}) + \frac{s^{2} t_{\rm s}}{\rho \rho _{\rm g}}\nabla P \cdot \nabla u, \end{equation}$$

the first term of which is identical to (35) and the second term of which is discretised in Phantom according to

(281) $$\begin{equation} -\frac{\rho _{a}}{2\rho ^{\rm g}_{a}} \sum _{b} m_{b}\frac{s_{a} s_{b}}{\rho _{a}\rho _{b}} \left(\frac{t_{{\rm s}, a}}{\rho _{a}} + \frac{t_{{\rm s}, b}}{\rho _b} \right) (P_{a} - P_{b}) (u_{a} - u_{b})\frac{\overline{F}_{ab}}{\vert r_{ab} \vert } . \end{equation}$$

2.13.15. Conservation of dust mass

Conservation of dust mass with the one-fluid scheme is in principle exact because (Price & Laibe Reference Price and Laibe2015a)

(282) $$\begin{equation} \frac{{\rm d}}{{\rm d}t} \left( \sum _{a} \frac{m_{a}}{\rho _{a}} s_{a}^{2}\right) = \sum _{a} m_{a} \left( \frac{2s_{a}}{\rho _{a}} \frac{{\rm d}s_{a}}{{\rm d}t} - \frac{s_{a}^{2}}{\rho _{a}^{2}} \frac{{\rm d}\rho _{a}}{{\rm d}t}\right) = 0. \end{equation}$$

In practice, some violation of this can occur because although the above algorithm guarantees positivity of the dust fraction, it does not guarantee that ε remains less than unity. Under this circumstance, which occurs only rarely, we set ε = max(ε, 1) in the code. However, this violates the conservation of dust mass. This specific issue has been recently addressed in detail in the study by Ballabio et al. (Reference Ballabio, Dipierro, Veronesi, Lodato, Hutchison, Laibe and Price2018). Therefore, in the latest code, there are two main changes as follows:

  • Rather than evolve $s = \sqrt{\epsilon \rho }$, in the most recent code, we instead evolve a new variable $s^{\prime } = \sqrt{\rho _{\rm d}/\rho _{\rm g}}$. This prevents the possibility of ε > 1.

  • We limit the stopping time such that the timestep from the one fluid algorithm does not severely restrict the computational performance.

For details of these changes, we refer the reader to Ballabio et al. (Reference Ballabio, Dipierro, Veronesi, Lodato, Hutchison, Laibe and Price2018).

2.13.16. Conservation of energy and momentum

Total energy with the one-fluid scheme can be expressed via

(283) $$\begin{equation} E =\sum _{a} m_{a} \left[ \frac{1}{2} v^{2}_{a} + (1 - \epsilon _{a}) u_{a} \right], \end{equation}$$

which is conserved exactly by the discretisation since

(284) $$\begin{eqnarray} &&\sum _{a} m_{a} \left[ {\bm v}_{a}\cdot \frac{{\rm d} {\bm v}_{a}}{{\rm d}t} + (1 - \epsilon _{a}) \frac{{\rm d}u_{a}}{{\rm d}t} \right. \nonumber \\ &&\left. -\, u_{a} \left(\frac{2s_{a}}{\rho _{a}}\frac{{\rm d}s_{a}}{{\rm d}t} - \frac{s^{2}_{a}}{\rho _{a}^{2}} \frac{{\rm d}\rho _{a}}{{\rm d}t}\right)\right] = 0. \end{eqnarray}$$

Conservation of linear and angular momentum also hold since the discretisation of the momentum equation is identical to the hydrodynamics algorithm.

2.13.17. Timestep constraint

For the one-fluid method, the timestep is constrained by the inverse of the stopping time according to

(285) $$\begin{equation} \Delta t^{a}_{\rm drag} = C_{\rm force} \frac{h^{2}}{ \epsilon t_{\rm s} c_{\rm s}^{2}}. \end{equation}$$

This becomes prohibitive precisely in the regime where the one-fluid method is no longer applicable (t s > t Cour; see Laibe & Price Reference Laibe and Price2014a), in which case one should switch to the two-fluid method instead. There is currently no mechanism to automatically switch regimes, though this is possible in principle and may be implemented in a future code version.

2.14. Interstellar medium (ISM) physics

2.14.1. Cooling function

The cooling function in Phantom is based on a set of Fortran modules written by Glover & Mac Low (Reference Glover and Mac Low2007), updated further by Glover et al. (Reference Glover, Federrath, Mac Low and Klessen2010). It includes cooling from atomic lines (H I), molecular excitation (H2), fine structure metal lines (Si I, Si II, O I, C I, and C II), gas–dust collisions, and polycyclic aromatic hydrocarbon (PAH) recombination (see Glover & Mac Low Reference Glover and Mac Low2007 for references on each process). Heating is provided by cosmic rays and the photoelectric effect. The strength of the cooling is a function of the temperature, density, and various chemical abundances. Table 2 summarises these various heating and cooling processes. Figure 7 shows the resulting emissivity ΛE(T) for temperatures between 104 and 108 K. The cooling rate per unit volume is related to the emissivity according to

(286) $$\begin{equation} \Lambda _{\rm cool} \equiv n^2 \Lambda _{\rm E} (T) \, {\rm erg} \phantom{s}{\rm s}^{-1} {\rm cm}^{-3}, \end{equation}$$

where n is the number density. The cooling in the energy equation corresponds to

(287) $$\begin{equation} \left(\frac{{\rm d}u}{{\rm d}t}\right)_{\rm cool} = -\frac{\Lambda _{\rm cool}}{\rho }. \end{equation}$$

These routines were originally adapted for use in sphng (Dobbs et al. Reference Dobbs, Glover, Clark and Klessen2008) and result in an approximate two-phase ISM with temperatures of 100 K and 10 000 K. Note that the cooling depends on a range of parameters (dust-to-gas ratio, cosmic ray ionisation rate, etc.), many of which can be specified at runtime. Table 3 lists the default abundances, which are set to values appropriate for the Warm Neutral Medium taken from Sembach et al. (Reference Sembach, Howk, Ryans and Keenan2000). The abundances listed in Table 3, along with the dust-to-gas ratio, are the only metallicity-dependent parameters that can be changed at runtime. An option for cooling appropriate to the zero metallicity early universe is also available.

Figure 7. Emissivity, ΛE (erg s−1 cm3) as a function of temperature for the ISM cooling assuming default abundances appropriate for the warm neutral medium (WNM). Note that as we treat cooling from atomic hydrogen using a full non-equilibrium treatment, the behaviour of ΛE close to 104 K is highly sensitive to the electron fraction, which in the case shown here is much smaller than it would be in collisional ionisation equilibrium. Values of ΛE below 104 K depend strongly on the current chemical state of the gas and are not shown in this plot.

Table 2. Heating and cooling processes in the Phantom cooling module.

Table 3. Default fractional abundances for C, O, Si, and e in the ISM cooling and chemistry modules. Abundances are taken from Sembach et al. (Reference Sembach, Howk, Ryans and Keenan2000) appropriate for the warm neutral medium (WNM). These are lower than solar because it is assumed some fraction of the metals are locked up in dust rather than being available in the gas phase.

2.14.2. Timestep constraint from cooling

When cooling is used, we apply an additional timestep constraint in the form

(288) $$\begin{equation} \Delta t _{\rm cool}^a = C_{\rm cool} \left| \frac{u}{({\rm d}u/{\rm d}t)_{\rm cool}} \right|, \end{equation}$$

where C cool = 0.3 following Glover & Mac Low (Reference Glover and Mac Low2007). The motivation for this additional constraint is to not allow the cooling to completely decouple from the other equations in the code (Suttner et al. Reference Suttner, Smith, Yorke and Zinnecker1997), and to avoid cooling instabilities that may be generated by numerical errors (Ziegler, Yorke, & Kaisig Reference Ziegler, Yorke and Kaisig1996).

Cooling is currently implemented only with explicit timestepping of (287), where u is evolved alongside velocity in our leapfrog timestepping scheme. However, the substepping scheme described below for the chemical network (Section 2.14.4) is also used to update the cooling on the chemical timestep, meaning that the cooling can evolve on a much shorter timestep than the hydrodynamics when it is used in combination with chemistry, which it is by default. Implementation of an implicit cooling method, namely the one proposed by Townsend (Reference Townsend2009), is under development.

2.14.3. Chemical network

A basic chemical network is included for ISM gas that evolves the abundances of H, H+, e, and the molecules H2 and CO. The number density of each species, n X, is evolved using simple rate equations of the form

(289) $$\begin{equation} \frac{dn_X}{dt} = C_X - D_X n_X, \end{equation}$$

where C X and D X are creation and destruction coefficients for each species. In general, C X and D X are functions of density, temperature, and abundances of other species. The number density of each species, X, is time integrated according to

(290) $$\begin{equation} n_X(t+\Delta t) = n_X(t) + \frac{dn_X}{dt} \Delta t. \end{equation}$$

There are in effect only three species to evolve (H, H2, and CO), as the H+ and e abundances are directly related to the H abundance.

The chemistry of atomic hydrogen is effectively the same as in Glover & Mac Low (Reference Glover and Mac Low2007). H is created by recombination in the gas phase and on grain surfaces, and destroyed by cosmic ray ionisation and free electron collisional ionisation. H2 chemistry is based on the model of Bergin et al. (Reference Bergin, Hartmann, Raymond and Ballesteros-Paredes2004), with H2 created on grain surfaces and destroyed by photo-dissociation and cosmic rays (see Dobbs et al. Reference Dobbs, Glover, Clark and Klessen2008 for computational details).

The underlying processes behind CO chemistry are more complicated, and involve many intermediate species in creating CO from C and O by interactions with H species. Instead of following every intermediate step, we use the model of Nelson & Langer (Reference Nelson and Langer1997) (see Pettitt et al. Reference Pettitt, Dobbs, Acreman and Price2014 for computational details). CO is formed by a gas phase reaction from an intermediate CHZ step after an initial reaction of C+ and H2 (where Z encompasses many similar type species). CO and CHZ are subject to independent photo-destruction, which far outweighs the destruction by cosmic rays. Abundances of C+ and O are used in the CO chemistry, and their abundance is simply the initial value excluding what has been used up in CO formation. Glover & Clark (Reference Glover and Clark2012) test this and a range of simpler and more complicated models, and show that the model adopted here is sufficient for large-scale simulations of CO formation, although it tends to over-produce CO compared to more sophisticated models.

The details for each reaction in the H, H2, and CO chemistry are given in Table 4, with relevant references for each process.

Table 4. Processes and references for the Phantom ISM chemistry module tracing the evolution of H, H2, and CO.

* H + 2 ions produced by cosmic ray ionisation of H2 are assumed to dissociatively recombine to H + H, so that the effective reaction in the code is actually H2 + c.r. → H + H.

† Process is intermediate and is assumed rather than fully represented.

†† C is not present in our chemistry, but is assumed to rapidly photoionise to C+.

2.14.4. Timestep constraint from chemistry

H chemistry is evolved on the cooling timestep, since the timescale on which the H+ abundance changes significantly is generally comparable to or longer than the cooling time. This is not true for H2 and CO. Instead, these species are evolved using a chemical timestepping criterion, where (290) is subcycled during the main time step at the interval Δt chem. If the abundance is decreasing, then the chemical timestep is

(291) $$\begin{equation} \Delta t_{\rm chem} = -\frac{1}{10} \frac{n_X}{(C_X - D_X n_X)}, \end{equation}$$

i.e., 10% of the time needed to completely deplete the species. If the abundance is increasing,

(292) $$\begin{equation} \Delta t_{\rm chem} = \frac{\Delta t_{\rm hydro}}{200} , \end{equation}$$

where Δt hydro is the timestep size for the hydrodynamics, and was found to be an appropriate value by test simulations. These updated abundances feed directly into the relevant cooling functions. Although the cooling function includes Si I and C I, the abundances of these elements are set to zero in the current chemical model.

2.15. Particle injection

We implement several algorithms for modelling inflow boundary conditions (see Toupin et al. Reference Toupin, Braun, Siess, Jorissen, Gail and Price2015a, Reference Toupin, Braun, Siess, Jorissen, Price, Kerschbaum, Wing and Hron2015b for recent applications) This includes injecting SPH particles in spherical winds from sink particles (both steady and time dependent), in a steady Cartesian flow and for injection at the L 1 point between a pair of sink particles to simulate the formation of accretion discs by Roche Lobe overflow in binary systems.

3 INITIAL CONDITIONS

3.1. Uniform distributions

The simplest method for initialising the particles is to set them on a uniform Cartesian distribution. The lattice arrangement can be cubic (equal particle spacing in each direction, Δx = Δy = Δz), close-packed ($\Delta y = \sqrt{3/4}\Delta x$, $\Delta z = \sqrt{6}/3 \Delta x$, repeated every three layers in z), hexagonal close-packed (as for close-packed but repeated every two layers in z), or uniform random. The close-packed arrangements are the closest to a ‘relaxed’ particle distribution, but require care with periodic boundary conditions due to the aspect ratio of the lattice. The cubic lattice is not a relaxed arrangement for the particles, but is convenient and sufficient for problems where the initial conditions are quickly erased (e.g. driven turbulence). For problems where initial conditions matter, it is usually best to relax the particle distribution by evolving the simulation for a period of time with a damping term (Section 3.6). This is the approach used, for example, in setting up stars in hydrostatic equilibrium (Section 3.4).

3.2. Stretch mapping

General non-uniform density profiles may be set up using ‘stretch mapping’ (Herant Reference Herant1994). The procedure for spherical distributions is the most commonly used (e.g. Fryer, Hungerford, & Rockefeller Reference Fryer, Hungerford and Rockefeller2007; Rosswog & Price Reference Rosswog and Price2007; Rosswog, Ramirez-Ruiz, & Hix Reference Rosswog, Ramirez-Ruiz and Hix2009), but we have generalised the method for any density profile that is along one coordinate direction (e.g. Price & Monaghan Reference Price and Monaghan2004b). Starting with particles placed in a uniform distribution, the key is that a particle should keep the same relative position in the mass distribution. For each particle with initial coordinate x 0 in the relevant coordinate system, we solve the equation

(293) $$\begin{equation} f(x) = \frac{M(x)}{M(x_{\rm max})} - \frac{x_{0} - x_{\rm min}}{(x_{\rm max} - x_{\rm min})} = 0, \end{equation}$$

where M(x) is the desired density profile integrated along the relevant coordinate direction, i.e.

(294) $$\begin{equation} M(x) \equiv \int _{x_{\rm min}}^{x} \rho (x^{\prime }) \text{d}S(x^{\prime }) \text{d}x^{\prime }, \end{equation}$$

where the area element dS(x′) depends on the geometry and the chosen direction, given by

(295) $$\begin{equation} dS(x) = \left\lbrace \arraycolsep5pt\begin{array}{ll}1 & \textrm {Cartesian or cyl./sph. along}\ \phi , \theta \ \text{or}\ z, \\ 2\pi x & \textrm {cylindrical along r}, \\ 4\pi x^{2} & \textrm {spherical along r}. \end{array}\right. \end{equation}$$

We solve (293) for each particle using Newton–Raphson iterations

(296) $$\begin{equation} x = x - \frac{f(x)}{f^{\prime }(x)}, \end{equation}$$

where

(297) $$\begin{equation} f^{\prime }(x) = \frac{\rho (x)\text{d}S(x)}{M(x_{\rm max})}, \end{equation}$$

iterating until |f(x)| < 10−9. The Newton–Raphson iterations have second-order convergence, but may fail in extreme cases. Therefore, if the root finding has failed to converge after 30 iterations, we revert to a bisection method, which is only first-order but guaranteed to converge.

Stretch mapping is implemented in such a way that only the desired density function need be specified, either as through an analytic expression (implemented as a function call) or as a tabulated dataset. Since the mass integral in (294) may not be known analytically, we compute this numerically by integrating the density function using the trapezoidal rule.

The disadvantage of stretch mapping is that in spherical or cylindrical geometry it produces defects in the particle distribution arising from the initial Cartesian distribution of the particles. In this case, the particle distribution should be relaxed into a more isotropic equilibrium state before starting the simulation. For stars, for example, this may be done by simulating the star in isolation with artificial damping added (Section 3.6). Alternative approaches are to relax the simulation using an external potential chosen to produce the desired density profile in equilibrium (e.g. Zurek & Benz Reference Zurek and Benz1986; Nagasawa, Nakamura, & Miyama Reference Nagasawa, Nakamura and Miyama1988) or to iteratively ‘cool’ the particle distribution to generate ‘optimal’ initial conditions (Diehl et al. Reference Diehl, Rockefeller, Fryer, Riethmiller and Statler2015).

3.3. Accretion discs

3.3.1. Density field

The accretion disc setup module uses a Monte Carlo particle placement (details in Section A.7) in cylindrical geometry to construct density profiles of the form

(298) $$\begin{equation} \rho (x,y,z) = \Sigma _{0} f_{\rm s} \left(\frac{R}{R_{\rm in}}\right)^{-p} \exp {\left(\frac{-z^{2}}{2H^{2}}\right)} , \end{equation}$$

where Σ0 is the surface density at R = R in (if f s = 1), Hc s/Ω is the scale height (with $\Omega \equiv \sqrt{GM/R^{3}}$), p is the power-law index (where p = 3/2 by default following Lodato & Pringle Reference Lodato and Pringle2007), and $f_{\rm s} \equiv (1 - \sqrt{R_{\rm in}/R} )$ is an optional factor to smooth the surface density near the inner disc edge.

Several authors have argued that a more uniform particle placement is preferable for setting up discs in SPH (Cartwright, Stamatellos, & Whitworth Reference Cartwright, Stamatellos and Whitworth2009; Vanaverbeke et al. Reference Vanaverbeke, Keppens, Poedts and Boffin2009). This may be important if one is interested in transient phenomena at the start of a simulation, but otherwise the particle distribution settles to a smooth density distribution within a few orbits (c.f. Lodato & Price Reference Lodato and Price2010).

3.3.2. Velocity field

The orbital velocities of particles in the disc are set by solving the equation for centrifugal equilibrium, i.e.

(299) $$\begin{equation} v_{\phi }^{2} = \frac{GM}{R} - f_{p} - 2v_{\phi } f_{\rm BH}, \end{equation}$$

where the correction from radial pressure gradients is given by

(300) $$\begin{equation} f_{p} = -c_{\rm s}^{2}(R) \left(\frac{3}{2} + p + q + \frac{1}{2f_{\rm s}}\right), \end{equation}$$

where q is the index of the sound speed profile such that c s(R) = c s, in(R/R in)q and f BH is a correction used for discs around a spinning black hole (Nealon et al. Reference Nealon, Price and Nixon2015)

(301) $$\begin{equation} f_{\rm BH} = - \frac{2a}{c^{3}} \left(\frac{GM}{R}\right)^{2}, \end{equation}$$

where a is the black hole spin parameter. The latter assumes Lense–Thirring precession is implemented as in Section 2.4.5. Where self-gravity is used, M is the enclosed mass at a given radius M( < R), otherwise it is simply the mass of the central object. Using the enclosed mass for the self-gravitating case is an approximation since the disc is not spherically symmetric, but the difference is small and the disc relaxes quickly into the true equilibrium. Equation (299) is a quadratic for v ϕ which we solve analytically.

3.3.3. Warps

Warps are applied to the disc (e.g. Lodato & Price Reference Lodato and Price2010; Nealon et al. Reference Nealon, Price and Nixon2015) by rotating the particles about the y-axis by the inclination angle i [in general, a function of radius ii(R)], according to

(302) $$\begin{equation} x^{\prime } = x\cos (i) + z\sin (i), \\ \end{equation}$$
(303) $$\begin{equation} y^{\prime } = y, \\ \end{equation}$$
(304) $$\begin{equation} z^{\prime } = -x\sin (i) + z\cos (i), \end{equation}$$
with the velocities similarly adjusted using
(305) $$\begin{equation} v^{\prime }_{x} = v_{x}\cos (i) + v_{z}\sin (i), \\ \end{equation}$$
(306) $$\begin{equation} v^{\prime }_{y} = v_{y}, \\ \end{equation}$$
(307) $$\begin{equation} v^{\prime }_{z} = -v_{x}\sin (i) + v_{z} \cos (i). \end{equation}$$

3.3.4. Setting an α-disc viscosity

The simplest approach to mimicking an α-disc viscosity in SPH is to employ a modified shock viscosity term, setting the desired αSS according to (124) as described in more detail in Section 2.6.1. Since the factor 〈h〉/H is dependent both on resolution and temperature profile (i.e. the q-index), it is computed in the initial setup by taking the desired αSS as input in order to give the required αAV. Although this does not guarantee that αSS is constant with radius and time (this is only true with specific choices of p and q and if the disc is approximately steady), it provides a simple way to prescribe the disc viscosity.

3.4. Stars and binary stars

We implement a general method for setting up ‘realistic’ stellar density profiles, based on either analytic functions (e.g. polytropes) or tabulated data files output from stellar evolution codes [see Iaconi et al. (Reference Iaconi, Reichardt, Staff, De Marco, Passy, Price, Wurster and Herwig2017) for a recent application of this to common envelope evolution].

The basic idea is to set up a uniform density sphere of particles and set the density profile by stretch mapping (see above). The thermal energy of the particles is set so that the pressure gradient is in hydrostatic equilibrium with the self-gravity of the star for the chosen EOS. We then relax the star into equilibrium for several dynamical times using a damping parameter (Section 3.6), before re-launching the simulation with the star on an appropriate orbit.

For simulating red giants, it is preferable to replace the core of the star by a sink particle (see Passy et al. Reference Passy2012; Iaconi et al. Reference Iaconi, Reichardt, Staff, De Marco, Passy, Price, Wurster and Herwig2017). When doing so, one should set the accretion radius of the sink to zero and set a softening length for the sink particle consistent with the original core radius (see Section 2.8.1).

3.5. Galactic initial conditions

In addition to simulating ISM gas in galactic discs with analytic stellar potentials, one may represent bulge-halo-disc components by collisionless N-body particles (see Section 2.12.3). To replace a potential with a resolved system requires care with the initial conditions (i.e. position, velocity, and mass). If setup incorrectly the system will experience radial oscillations and undulations in the rotation curve, which will have strong adverse effects on the gas embedded within. We include algorithms for initialising the static-halo models of Pettitt et al. (Reference Pettitt, Dobbs, Acreman and Bate2015) (which used the sphng SPH code). These initial conditions require the NFW profile to be active and care must be taken to ensure the mass and scale lengths correspond to the rotation curve used to generate the initial conditions. Other codes may alternatively be used to seed multi-component N-body disc galaxies (e.g. Kuijken & Dubinski Reference Kuijken and Dubinski1995a; Boily, Kroupa, & Peñarrubia-Garrido Reference Boily, Kroupa and Peñarrubia-Garrido2001; McMillan & Dehnen Reference McMillan and Dehnen2007; Yurin & Springel Reference Yurin and Springel2014), including magalie (Boily et al. Reference Boily, Kroupa and Peñarrubia-Garrido2001) and galic (Yurin & Springel Reference Yurin and Springel2014) for which we have implemented format readers.

The gas in galactic scale simulations can be setup either in a uniform surface density disc, or according to the Milky Way’s specific surface density. The latter is based on the radial functions given in Wolfire et al. (Reference Wolfire, Hollenbach, McKee, Tielens and Bakes1995). As of yet, we have not implemented a routine for enforcing hydrostatic equilibrium (Springel, Di Matteo, & Hernquist Reference Springel, Di Matteo and Hernquist2005; Wang et al. Reference Wang, Klessen, Dullemond, van den Bosch and Fuchs2010), this may be included in a future update.

3.6. Damping

To relax a particle distribution into equilibrium, we adopt the standard approach (e.g. Gingold & Monaghan Reference Gingold and Monaghan1977) of adding an external acceleration in the form

(308) $$\begin{equation} {\bm a}_{\rm ext, damp}^{a} = - f_{\rm d} {\bm v}, \end{equation}$$

such that a percentage of the kinetic energy is removed each timestep. The damping parameter, f d, is specified by the user. A sensible value for f d is of order a few percent (e.g. f d = 0.03) such that a small fraction of the kinetic energy is removed over a Courant timescale.

4 SOFTWARE ENGINEERING

No code is completely bug free (experience is the name everyone gives to their mistakes; Wilde Reference Wilde1892). However, we have endeavoured to apply the principles of good software engineering to Phantom. These include

  1. 1. a modular structure,

  2. 2. unit tests of important modules,

  3. 3. nightly builds,

  4. 4. automated nightly tests,

  5. 5. automated code maintenance scripts,

  6. 6. version control with git,

  7. 7. wiki documentation, and

  8. 8. a bug database and issue tracker.

Together these simplify the maintenance, stability, and usability of the code, meaning that Phantom can be used direct from the development repository without fear of regression, build failures, or major bugs.

Specific details of how the algorithms described in Section 2 are implemented are given in the Appendix. Details of the test suite are given in Appendix A.9.

5 NUMERICAL TESTS

Unless otherwise stated, we use the M 6 quintic spline kernel with h fac = 1.0 by default, giving a mean neighbour number of 113 in 3D. Almost all of the test results are similar when adopting the cubic spline kernel with h fac = 1.2 (requiring ≈58 neighbours), apart from the tests with the one-fluid dust method where the quintic is required. Since most of the algorithms used in Phantom have been extensively tested elsewhere, our aim is merely to demonstrate that the implementation in the code is correct, and to illustrate the typical results that should be achieved on these tests when run by the user. The input files used to run the entire test suite shown in the paper are available on the website, so it should be straightforward for even a novice user to reproduce our results.

Unless otherwise indicated, we refer to dimensionless ${\cal L}_{1}$ and ${\cal L}_{2}$ norms when referencing errors, computed according to

(309) $$\begin{eqnarray} {\cal L}_{1} \equiv \frac{1}{N C_0}\sum _{i=1}^{N} \vert y_{i} - y_{\rm exact} \vert , \end{eqnarray}$$
(310) $$\begin{eqnarray} {\cal L}_{2} \equiv \sqrt{\frac{1}{N C_0}\sum _{i=1}^{N} \vert y_{i} - y_{\rm exact} \vert ^{2}}, \end{eqnarray}$$

where y exact is the exact or comparison solution interpolated or computed at the location of each particle i and N is the number of points. The norms are the standard error norms divided by a constant, which we set to the maximum value of the exact solution within the domain, C 0 = max(y exact), in order to give a dimensionless quantity. Otherwise quoting the absolute value of ${\cal L}_{1}$ or ${\cal L}_{2}$ is meaningless. Dividing by a constant has no effect when considering convergence properties. These are the default error norms computed by splash (Price Reference Price2007).

Achieving formal convergence in SPH is more complicated than in mesh-based codes where linear consistency is guaranteed (see Price Reference Price2012a). The best that can be achieved with regular (positive, symmetric) kernels is second-order accuracy away from shocks provided the particles remain well ordered (Monaghan Reference Monaghan1992). The degree to which this remains true depends on the smoothing kernel and the number of neighbours. Our tests demonstrate that formal second-order convergence can be achieved with Phantom on certain problems (e.g. Section 5.6.1). More typically, one obtains something between first- and second-order convergence in smooth flow depending on the degree of spatial disordering of the particles. The other important difference compared to mesh-based codes is that there is no intrinsic numerical dissipation in SPH due to its Hamiltonian nature—numerical dissipation terms must be explicitly added. We perform all tests with these terms included.

We use timestep factors of C Cour = 0.3 and C force = 0.25 by default for all tests (Section 2.3.2).

5.1. Hydrodynamics

5.1.1. Sod shock tube

Figure 8 shows the results of the standard Sod (Reference Sod1978) shock tube test, performed in 3D using [ρ, P] = [1, 1] in the ‘left state’ (x ⩽ 0) and [ρ, P] = [0.125, 0.1] for the ‘right state’ (x > 0) with a discontinuity initially at x = 0 and zero initial velocity and magnetic field. We perform the test using an adiabatic EOS with γ = 5/3 and periodic boundaries in y and z. While many 1D solutions appear in the literature, only a handful of results on this test have been published for SPH in 3D (e.g. Dolag et al. Reference Dolag, Vazza, Brunetti and Tormen2005; Hubber et al. Reference Hubber, Batty, McLeod and Whitworth2011; Beck et al. Reference Beck2016; a 2D version is shown in Price Reference Price2012a). The tricky part in a 3D SPH code is how to set up the density contrast. Setting particles on a cubic lattice is a poor choice of initial condition since this is not a stable arrangement for the particles (Morris Reference Morris1996b, Reference Morris1996a; Lombardi et al. Reference Lombardi, Sills, Rasio and Shapiro1999; Børve et al. Reference Børve, Omang and Trulsen2004). The approach taken in Springel (Reference Springel2005) (where only the density was shown, being the easiest to get right) was to relax the two halves of the box into a stable arrangement using the gravitational force run with a minus sign, but this is time consuming.

Figure 8. Results of the Sod shock tube test in 3D, showing projection of all particles (black dots) compared to the analytic solution (red line). The problem is set up with [ρ, P] = [1, 1] for x ⩽ 0 and [ρ, P] = [0.125, 0.1] for x > 0 with γ = 5/3, with zero velocities and no magnetic field. The density contrast is initialised using equal mass particles placed on a close packed lattice with 256 × 24 × 24 particles initially in x ∈ [− 0.5, 0] and 128 × 12 × 12 particles initially in x ∈ [0, 0.5]. The results are shown with constant αAV = 1.

Here we take a simpler approach which is to set the particles initially on a close-packed lattice (Section 3.1), since this is close to the relaxed arrangement (e.g. Lombardi et al. Reference Lombardi, Sills, Rasio and Shapiro1999). To ensure continuity of the lattice across periodic boundaries, we fix the number of particles in the y (z) direction to the nearest multiple of 2 (3) and adjust the spacing in the x-direction accordingly to give the correct density in each subdomain. We implement the boundary condition in the x-direction by tagging the first and last few rows of particles in the x direction as boundary particles, meaning that their particle properties are held constant. The results shown in Figure 8 use 256 × 24 × 24 particles initially in x ∈ [− 0.5, 0] and 128 × 12 × 12 particles in x ∈ [0, 0.5] with code defaults for the artificial conductivity [αu = 1 with vu sig given by equation (43)] and artificial viscosity (αAV = 1, βAV = 2). The results are identical whether global or individual particle timesteps (Section 2.3.4) are used. Figure 9 shows the results when code defaults for viscosity are also employed, resulting in a time-dependent αAV (see lower right panel). There is more noise in the velocity field around the contact discontinuity in this case, but the results are otherwise identical.

Figure 9. As in Figure 8, but with code defaults for all dissipation terms (αAV ∈ [0, 1]; αu = 1). These leave more noise in the velocity field behind the shock but provide second-order convergence in smooth flow. The lower right panel in this case shows the resultant values for the viscosity parameter αAV.

The ${\cal L}_2$ errors for the solutions shown in Figure 8 are 0.0090, 0.0022, 0.0018, and 0.0045 for the density, velocity, thermal energy, and pressure, respectively. With dissipation switches turned on (Figure 9), the corresponding errors are 0.009, 0.0021, 0.0019, and 0.0044, respectively. That is, our solutions are within 1% of the analytic solution for all four quantities in both cases, with the density profile showing the strongest difference (mainly due to smoothing of the contact discontinuity).

Puri & Ramachandran (Reference Puri and Ramachandran2014) compared our shock capturing algorithm with other SPH variants, including Godunov SPH. Our scheme was found to be the most robust of those tested.

5.1.2. Blast wave

As a more extreme version of the shock test, Figures 10 and 11 show the results of the blast wave problem from Monaghan (Reference Monaghan1997), set up initially with [ρ, P] = [1, 1000] for x ⩽ 0 and [ρ, P] = [1.0, 0.1] for x > 0 and with γ = 1.4 (appropriate to air). As previously, we set the particles on a close-packed lattice with a discontinuous initial pressure profile. We employ 800 × 12 × 12 particles in the domain x ∈ [− 0.5, 0.5]. Results are shown at t = 0.01. This is a more challenging problem than the Sod test due to the higher Mach number. As previously, we show results with both αAV = 1 (Figure 10) and with the viscosity switch α ∈ [0, 1]. Both calculations use αu = 1 with (43) for the signal speed in the artificial conductivity. For the solution shown in Figure 10, we find normalised ${\cal L}_2$ errors of 0.057, 0.063, 0.051, and 0.018 in the density, velocity, thermal energy, and pressure, respectively, compared to the analytic solution. Employing switches (Figure 11), we find corresponding ${\cal L}_2$ errors of 0.056, 0.059, 0.052, and 0.017. That is, our solutions are within 6% of the analytic solution at this resolution.

Figure 10. Results of the 3D blast wave test, showing projection of all particles (black dots) compared to the analytic solution (red line). The problem is set up with [ρ, P] = [1, 1000] for x ⩽ 0 and [ρ, P] = [1, 0.1] for x > 0 with γ = 7/5, with zero velocities and no magnetic field. We use equal mass particles placed on a close packed lattice with 800 × 12 × 12 particles initially in x ∈ [− 0.5, 0.5]. Results are shown with constant αAV = 1.

Figure 11. As in Figure 10, but with code defaults for dissipation switches (αAV ∈ [0, 1]; αu = 1). As in Figure 9, the velocity field behind the shock is more noisy with switches applied, but the switches reduce the numerical dissipation away from shocks.

The main source of error is that the contact discontinuity is over-smoothed due to the artificial conductivity, while the velocity profile shows a spurious jump at the location of the contact discontinuity. This glitch is a startup error caused by our use of purely discontinuous initial conditions—it could be removed by adopting smoothed initial conditions but we prefer to perform the more difficult version of this test. There is also noise in the post-shock velocity profile because the viscosity switch does not apply enough dissipation here. As in the previous test, this noise can be reduced by increasing the numerical viscosity, e.g. by adopting a constant αAV (compare Figures 10 and 11).

Figure 12 quantifies the error in energy conservation, showing the error in the total energy as a function of time, i.e. |EE 0|/|E 0|. Energy is conserved to a relative accuracy of better than 2 × 10−6.

Figure 12. Relative error in energy conservation for the 3D blast wave test. Energy is conserved to an error of 10−6 in this test.

5.1.3. Sedov blast wave

The Sedov–Taylor blast wave (Taylor Reference Taylor1950a, Reference Taylor1950b; Sedov Reference Sedov1959) is a similar test to the previous but with a spherical geometry. This test is an excellent test for the individual timestepping algorithm, since it involves propagating a blast wave into an ambient medium of ‘inactive’ or ‘asleep’ particles, which can cause severe loss of energy conservation if they are not carefully awoken (Saitoh & Makino Reference Saitoh and Makino2009). For previous uses of this test with SPH, see e.g. Springel & Hernquist (Reference Springel and Hernquist2002) and Rosswog & Price (Reference Rosswog and Price2007), and for a comparison between SPH and mesh codes on this test, see Tasker et al. (Reference Tasker, Brunino, Mitchell, Michielsen, Hopton, Pearce, Bryan and Theuns2008).

We set up the problem in a uniform periodic box x, y, z ∈ [− 0.5, 0.5], setting the thermal energy on the particles to be non-zero in a barely resolved sphere around the origin. We assume an adiabatic EOS with γ = 5/3. The total energy is normalised such that the total thermal energy in the blast is E 0 = ∑am au a = 1, distributed on the particles within r < R kernh 0 using the smoothing kernel, i.e.

(311) $$\begin{equation} u_{a} = \left\lbrace \arraycolsep5pt\begin{array}{@{}ll@{}}E_{0}W(r, h_{0}), & r/h_{0} \le R_{\rm kern}\\ 0 & r/h_{0} > R_{\rm kern}, \end{array}\right. \end{equation}$$

where $r = \sqrt{x^{2} + y^{2} + z^{2}}$ is the radius of the particle and we set h 0 to be twice the particle smoothing length.

We simulate the Sedov blast wave using both global and individual timesteps at several different resolutions. Figure 13 shows the evolution of the relative error in total energy for our suite, while Figure 14 shows the density at t = 0.1 compared to the analytical solution given by the solid line. Energy is conserved to better than 1% in all cases. Using higher spatial resolution results in a better match of the post-shock density with the analytic solution. The scatter at the leading edge of the shock is a result of the default artificial conductivity algorithm. Given the initial strong gradient in thermal energy, using more artificial conductivity than our switch permits would reduce the noise on this problem (Rosswog & Price Reference Rosswog and Price2007).

Figure 13. Evolution of the relative error in total energy for the Sedov blast wave problem. Resolutions are given in the legend; solid lines use individual timestepping while dotted lines show global timestepping.

Figure 14. Density as a function of radius in the Sedov blast wave problem at three resolutions. All particles are placed initially on a closepacked lattice, are evolved using individual timesteps and we use the quintic kernel. The analytic solution is given by the solid line, and the bottom panels show the residuals compared to the analytical solution.

5.1.4. Kelvin–Helmholtz instability

Much has been written about Kelvin–Helmholtz instabilities with SPH (e.g. Agertz et al. Reference Agertz2007; Price Reference Price2008; Abel Reference Abel2011; Valdarnini Reference Valdarnini2012; Read & Hayfield Reference Read and Hayfield2012; Hubber, Falle, & Goodwin Reference Hubber, Falle and Goodwin2013c). For the present purpose, it suffices to say that the test problems considered by Agertz et al. (Reference Agertz2007) and Price (Reference Price2008) are not well posed. That is, the number of small-scale secondary instabilities will always increase with numerical resolution because high wavenumber modes grow fastest in the absence of physical dissipation or other regularising forces such as magnetic fields or surface tension. The ill-posed nature of the test problem has been pointed out by several authors (Robertson et al. Reference Robertson, Kravtsov, Gnedin, Abel and Rudd2010; McNally, Lyra, & Passy Reference McNally, Lyra and Passy2012; Lecoanet et al. Reference Lecoanet2016), who have each proposed well-posed alternatives.

We adopt the setup from Robertson et al. (Reference Robertson, Kravtsov, Gnedin, Abel and Rudd2010), similar to the approach by McNally et al. (Reference McNally, Lyra and Passy2012), where the initial density contrast is smoothed using a ramp function. This should suppress the formation of secondary instabilities long enough to allow a single large-scale mode to grow. The density and shear velocity in the y direction are given by

(312) $$\begin{equation} \rho (y) = \rho _{1} + R(y)[\rho _{2} - \rho _{1}], \end{equation}$$

and

(313) $$\begin{equation} v_{x}(y) = v_{1} + R(y)[v_{2} - v_{1}], \end{equation}$$

where ρ1 = 1, ρ2 = 2, v 1 = −0.5, and v 2 = 0.5 with constant pressure P = 2.5, γ = 5/3. The ramp function is given by

(314) $$\begin{equation} R(y) \equiv \left[1 - f(y)\right]\left[1 -g(y)\right], \end{equation}$$

where

(315) $$\begin{eqnarray} f & \equiv& \frac{1}{1 + \exp \left[2(y-0.25)