Phantom: A Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code for Astrophysics

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 (h^d _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, new} − h _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⁻⁴ 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 ³ 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 ₄ 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 ²) 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 ³ 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 ₅ 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 ₆ 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 ²:

(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 ⁴ 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 ⁶ 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 ⁶ 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 ⁴ 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 ₆ or C ⁶ 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 x–y 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.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 q^a_ab and q^b_ab 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 α^AV_a ∈ [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 ²). 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 ² + y ² + z ²)^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⁻¹⁶ 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⁻³ and linear and angular momentum conservation to ~10⁻⁶ with default code settings. The code execution is aborted if conservation errors exceed 10%.

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 ⁿ 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 Δt ≡ t ^{n + 1} − t ⁿ. 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⁻². 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 ^dt_sig 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 − 1}_jΔ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, min}/Δt _{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.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 ₁ and m ₂. 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 ₀, 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 ⁴. 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.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 ₀, k ₁, and k ₂, whose default values are unity. Rotation of the central object can be crudely emulated by changing k ₂.

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.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 ²/ν 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.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 ⁵ 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 (S^ij _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 S^ij _NS = S _NS^ji, 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 ∑_j∑_iT _ijT _ij, where T_ij ≡ ∂vⁱ_a/∂x^j_a + ∂v_a ^j/∂xⁱ_a. 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.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, a − i, 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, b − c, 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.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 μ₀ 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 M^ij_a 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 μ₀ = 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 μ₀ in the relevant unit system. In particular, it avoids the differences between SI and cgs units in how the charge unit is defined, since μ₀ 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 v^B _{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, v^B _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 α^B_max = 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 ²). 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⁻² 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.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.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 ⁶

(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.3–2.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, q ≡ r/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.5–2.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 ²) 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 q ⩾ R _kern, while the first term simply becomes 1/r ². 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 rⁱ ≡ xⁱ_n − xⁱ_m is the relative position vector, $\hat{r}$ is the corresponding unit vector, M _m is the total mass in node m, Q^m_ij 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 Δxⁱ_a ≡ xⁱ_a − xⁱ ₀ 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.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/cm³ 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⁻¹³g/cm³ and sound speed c _s = 6 × 10⁴ cm/s with subsonic drag Δv = 0.01c _s and material density ρ_grain = 3g/cm³.

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 h ≲ c _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 ⁷. 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.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 (H₂), 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 10⁴ and 10⁸ 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⁻¹ cm³) 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 10⁴ 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 10⁴ 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 H₂ 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, H₂, 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. H₂ chemistry is based on the model of Bergin et al. (Reference Bergin, Hartmann, Raymond and Ballesteros-Paredes2004), with H₂ 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 CH_Z step after an initial reaction of C⁺ and H₂ (where Z encompasses many similar type species). CO and CH_Z 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, H₂, 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, H₂, and CO.

* H ⁺ ₂ ions produced by cosmic ray ionisation of H₂ are assumed to dissociatively recombine to H + H, so that the effective reaction in the code is actually H₂ + 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 H₂ 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.