© ESO 2019

1. Introduction

The study of protoplanetary disks in nonisolated system has become a relevant topic in numerical astrophysics because several planetary systems and disks around stars inside open clusters have recently been observed. The recent discovery of a Neptune-sized planet that is hosted by a binary star in the nearby Hyades cluster (Ciardi et al. 2018) has opened new perspectives for the study of the evolution of primordial disks that interact with binary stars in a nonisolated environment. In the context of the study of isolated binary systems, we note that most of the classical models have studied low-mass disks, with the result that even considering their self-gravity, no appreciable change is observed in the time evolution of the stellar orbital parameters. Additionally, low-mass disks give a poor feedback on the hosting stars, which means that the timescale for the variation in their orbital parameters is larger than the other dynamical timescales that are involved.

To investigate such systems, we built our own smoothed particle hydrodynamics (SPH) code to integrate the evolution of the composite star+gas system. Following the scheme described in Capuzzo-Dolcetta & Miocchi (1998), our code treats gravity by means of a classical tree-based scheme (see also Barnes 1986; Barnes & Hut 1986; Miocchi & Capuzzo-Dolcetta 2002) and adds a proper treatment of the close gravitational interactions of the gas particles. The evolution of a limited number of point-mass objects (which may represent either stars or planets) is treated with a high-order explicit method.

This preliminary work provides an instrument that is not only suited to modeling heavy protoplanetary disks that interact with single and binary stars. It can also be used to study these systems not in isolation, but in stellar systems such as open star clusters.

In the next section we describe our code after preliminarily introducing the numerical framework, and in Sect. 3 we present and discuss some physical and performance tests. In Sect. 4 we describe the disk model we adopted and the application of our code to the study of heavy Keplerian disks. Section 5 is dedicated to the conclusions.

2. Numerical algorithm

In order to study a selg-gravitating gas, we developed an SPH code, coupled with a tree-based scheme for the Newtonian force integration. In this section, after a breaf recall of some basic theory of gravitational tree codes and of the SPH formalism (Sect. 2.1), we describe the main features of our numerical algorithm (Sect. 2.2). For further technical details, we refer to Appendices A.1 and A.2.

2.1. Basic theory

Introduced for the first time by Lucy (1977) and Gingold & Monaghan (1977), the SPH scheme has been widely adopted to investigate a huge set of astronomical problems that involve fluid systems. An SPH scheme allows integrating the fluid dynamical equations in a Lagrangian approach by representing the system through a set of points, or so-called pseudo-particles. For each particle, a set of fundamental quantities (such as density ρ, pressure P, internal energy u, and velocity v) are calculated by means of an interpolation with a proper kernel function over a suitable neighbor. For an exhaustive explanation of the method, we refer to various papers in the literature, such as those by Monaghan & Lattanzio (1985), Monaghan (1988, 2005). Here we recall some basic aspects.

Interpolations are performed with a continuous kernel function W(r, h), whose spread scale is determined by a characteristic length h, called the smoothing length. It can easily be shown (see, e.g., Hernquist & Katz 1989) that under some additional constraints, interpolation errors are limited to the order O(h²). We used as kernel function the cubic spline that was for the first time adopted by Monaghan & Lattanzio (1985), who developed a formalism that had been introduced by Hockney & Eastwood (1981). This kernel function has the following form:

$$\begin{array}{c}\hfill W(r,h)=\frac{1}{\pi h^3}·\{\begin{array}{cc}{\displaystyle 1-{\displaystyle \frac{3}{2}}{\left({\displaystyle {\displaystyle \frac{r}{h}}}\right)}^2+{\displaystyle \frac{3}{4}}{\left({\displaystyle \frac{r}{h}}\right)}^3,}\hfill & 0\le r<h\hfill \\ {\displaystyle {\displaystyle \frac{1}{4}}{\left({\displaystyle 2-{\displaystyle \frac{r}{h}}}\right)}^3,}\hfill & h\le r<2h\hfill \\ 0,\hfill & r\ge 2h\hfill \end{array}.\end{array}$$(1)

The SPH interpolation involves only a limited set of N′ neighboring particles that are enclosed within the range 2h, therefore the computational effort is expected to scale linearly with the total particle number N. On the other hand, when long-range interactions such as gravity are considered, the computational effort increases because each particle interacts with the whole system. A classical direct N-body code would therefore require a computational weight scaling as N². However, a suitable gravitational tree-based scheme allows us to efficiently evaluate the Newtonian force by approximating the potential with a harmonic expansion (see Barnes & Hut 1986; Barnes 1986, for a full explanation). For each particle, only the contribution from a local neighborhood is calculated through a direct particle–particle coupling, while the contribution from farther particles is suitably approximated. The following expressions (Eqs. (2) and (3)) represent the approximated potential Φ(r) and the force (per unit mass) a(r)= − ∇Φ(r) given by a far cluster of particles,

$$\begin{array}{c}\hfill \mathrm{\Phi }(\mathit{r})=-{\displaystyle \frac{\mathit{GM}}{r}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{G\mathit{r}\overline{\overline{\mathit{Q}}}\mathit{r}}{r^5}}\end{array}$$(2)

$$\begin{array}{c}\hfill \mathit{a}(\mathit{r})=-{\displaystyle \frac{\mathit{GM}}{r^3}}\mathit{r}+{\displaystyle \frac{G\overline{\overline{\mathit{Q}}}·\mathit{r}}{r^5}}\mathit{r}-{\displaystyle \frac{5}{2}}{\displaystyle \frac{G\mathit{r}\overline{\overline{\mathit{Q}}}\mathit{r}}{r^7}}\mathit{r}.\end{array}$$(3)

M is the total mass of this ensemble, r = |r| is the distance of the particle under study to the center of mass of the cluster. The symbol $\overline{\overline{\mathit{Q}}}$ represents the so-called quadrupole tensor, which is associated with the specific cluster. In indexed form, it is given by

$$\begin{array}{c}\hfill Q_{\mathit{ij}}={\displaystyle \sum _{k=1}^{N_C}(3x_i^{(k)}{x}_j^{(k)}-r_k^2{\delta }_{\mathit{ij}})m_k,}\end{array}$$(4)

where $x_i^{(k)}$ and $x_i^{(k)}$ (i, j = 1, 2, 3) refer to the Cartesian coordinates of the kth particle of mass m_k. The summation is performed over all the N_C particles that are included in the cluster.

In the following section we describe the main structure and formalism used by GaSPH. Further computational details related to the implementation of the algorithm can be found in Appendices A and B.

2.2. Main structure of the algorithm

A single step to compute the acceleration contains two preliminary phases. One cycle is dedicated to map the particles into an octal grid domain. A further cycle, linear in N, is needed to evaluate some key parameters such as the density, ρ, the smoothing length, h, and the pressure, P. Then a third set of operations, the most complicated operation, is that of evaluating the gravitational and hydrodynamical forces in addition to the internal energy rate $\dot{u}$ of the gas.

GaSPH can also easily treat a system that consists of a set of point masses. To do this, the part of the SPH computations is turned off and the tree scheme alone is used for gravity interactions. On the other hand, a gas can be treated with pure hydrodynamics by turning off the gravitational field and using the SPH formalism alone.

After the main computations of the acceleration a and the energy rate $\dot{u}$, the algorithm updates in time the velocity, position, and internal energy of the gas with a second-order Verlet method. Because of the structure of the second-order technique, the three main computational cycles should be performed twice into a single time iteration to obtain two estimates of a and $\dot{u}$.

In addition, the smooth particles may interact with a small number N_ob of additional objects, an ensemble of point masses that mimic stars and/or planets. Differently from the other particles, the motion of these few objects is integrated with a 14th-order Runge–Kutta method by direct particle-to-particle N-body interactions without any approximation for the gravitational field. When N_ob is sufficiently small, these operations request a little additional computational effort that scales roughly linearly with respect to the total number of points (including both the SPH particles number N and the objects number N_ob). For the specific purpose of our investigation, where we have N_ob ≤ 2, the general efficiency of the code is not affected.

2.2.1. Particle mapping and density computation

For a set of N equal-mass points, we preliminarily need to subdivide the system into a hierarchical series of subgroups of points in order to apply the multipole approximation for the Newtonian field contribution given by a “cluster”. To do this, we use a classical Barnes-Hut tree-code to map the particles into an octal grid space, according to their positions. We follow in particular the technique adopted by Miocchi & Capuzzo-Dolcetta (2002) by mapping the points through a 3-bit-based codification (see Sect. 2.2.6 for further details).

Before the accelerations are computed, SPH particles need a preliminary stage in which densities and smoothing lengths are computed. To perform a good interpolation, we need to keep a fixed number of neighbors for each point. For inhomogeneous fluids, we must therefore use a smoothing length h ≡ h(r, t) that varies in space and in time.

Individual smoothing lengths should be chosen in such a way that the higher the local number density n = ρ/m, the smaller the interpolation kernel radius: h ∝ n^−1/3, in order to have a roughly constant number of neighbors of the given particle. For this purpose, we adopt a commonly used prescription (Hernquist & Katz 1989; Monaghan 2005). For each particle, we start from an initial guess for h, then we vary it until the number of particles that lie within the kernel dominion reaches a fixed value N₀. We iterate a process in which each time the number of neighbor points, N′, is counted using a certain smoothing length h_prev, then we update this length to a new value h_new according to the following formula:

$$\begin{array}{c}\hfill h_{\mathrm{new}}=h_{\mathrm{prev}}{\displaystyle \frac{1}{2}}[1+{\left({\displaystyle \frac{N_0}{N\prime }}\right)}^{1/3}].\end{array}$$(5)

If the fluid were homogeneous, h_prev(N₀/N′)^1/3 would immediately provide the correct value of the smoothing length, without any further iteration. The addend 1 lets the program perform an average of the old smoothing length, for which any excessive oscillation error due to non-homogeneities in the spatial distribution of particles is damped. The iteration is stopped when convergence is reached according to the criterion |N′−N₀|≤ΔN, where ΔN is a tolerance number. In this regard, Attwood et al. (2007) investigated the acoustic oscillations of some models of polytrope around the equilibrium by imposing a constant neighbor number N′ and letting ΔN vary. They found that the fluctuation of N′±ΔN introduced an additional numerical noise that was able to break the stability of this system, giving rise to errors. To prevent errors, the authors found that ΔN should be set to zero, which is the choice we adopt in this paper. Moreover, they showed that the calculation of h according to the iterative process illustrated above and with ΔN = 0 is equivalent to solving for all the particles the 2N-equations system described by the following two equations:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\rho }_i=\sum _{j=1}^{N\prime }{m}_{j}W(r_{\mathit{ij}},h_i)\hfill \\ h_i=\delta {\left({\displaystyle \frac{m_i}{{\rho }_i}}\right)}^{1/3}\hfill \end{array},\end{array}$$(6)

and to finding the exact solutions of density and smoothing length {ρ_i, h_i | i ∈ [1, N]}, with $\delta \approx 0.31 N_0^{1/3}$ a suitable constant, and r_ij the mutual distance between the ith and the jth particles. We typically use a number of neighbors N′=60, such as δ ≈ 1.2.

When the density ρ_i is evaluated, the corresponding pressure P_i can be computed by means of a suitable equation of state. Appendix B.1 illustrates further technical details about the neighbor-search procedure.

2.2.2. Force calculation and softened interactions

For a generic ith particle, the acceleration a_i is computed by adding both the SPH terms and the Newtonian terms in the same iteration. Together with the acceleration, the internal energy rate of the particle is also computed.

To treat a self-gravitating gas with an SPH scheme, a proper treatment of the gravitational potential is necessary to avoid overestimating the gravity field. Particles can be considered as point sources of the Newtonian field when their mutual distance is larger than 2h. Otherwise, their Newtonian interaction is, in consistency with the assumed kernel function (Gingold & Monaghan 1977), such that it vanishes at an interparticle distance approaching zero.

With the cubic spline kernel, we can obtain a different form of the Newtonian interaction between two particles, such that the classical term is softened if the particles approach within a distance of the order of a softening length ϵ = 2h. See the appendix in Hernquist & Katz (1989) for more details, and the appendix in Price & Monaghan (2007) for an explicit expression of the force and the potential. When SPH interaction is turned off, a constant value ϵ is generally used in place of 2h for the softening length. In this case, the total energy is conserved within the numerical error. On the other hand, the Hamiltonian becomes time dependent with SPH systems because the softening length varies in time, and so the energy is no longer conserved. To solve this problem, the equations of motion must be rewritten in a conservative form that takes the variation in h into account. We followed the Hamiltonian formalism as adopted for the first time by Springel & Hernquist (2002) for the hydrodynamical interactions, which was further developed by Price & Monaghan (2007) for the gravitational field. The SPH equation assume the form

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{\mathit{v}}_i}{\mathrm{d}t}}& =-\sum _j\frac{1}{2}(g_{\mathrm{soft}}(r_{\mathit{ij}},h_i)+g_{\mathrm{soft}}(r_{\mathit{ij}},h_j))\frac{{\mathit{r}}_{\mathit{ij}}}{r_{\mathit{ij}}}\hfill \\ \hfill & -\sum _j{m}_j\frac{G}{2}(\frac{{\zeta }_i}{{\mathrm{\Omega }}_i}{\mathbf{\nabla }}_{i}W(r_{\mathit{ij}},h_i)+\frac{{\zeta }_j}{{\mathrm{\Omega }}_j}{\mathbf{\nabla }}_{i}W(r_{\mathit{ij}},h_j))\hfill \\ \hfill & -\sum _j{m}_j({\displaystyle \frac{P_i}{{\rho }_i^2{\mathrm{\Omega }}_i}}{\mathbf{\nabla }}_{i}W(r_{\mathit{ij}},h_i)+{\displaystyle \frac{P_j}{{\rho }_j^2{\mathrm{\Omega }}_j}}{\mathbf{\nabla }}_{i}W(r_{\mathit{ij}},h_j))\hfill \\ \hfill & -\sum _j{m}_j{\mathrm{\Pi }}_{\mathit{ij}}\left[\frac{{\mathbf{\nabla }}_{\mathit{i}}W(r_{\mathit{ij}},h_i)+{\mathbf{\nabla }}_{i}W(r_{\mathit{ij}},h_j)}{2}\right]+{\displaystyle \frac{\mathrm{d}{\mathit{v}}_i^{[\mathrm{stars}]}}{\mathrm{d}t}}\hfill \end{array}$$(7)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{u}_i}{\mathrm{d}t}}=\sum _j{m}_j({\displaystyle \frac{P_i}{{\rho }_i^2{\mathrm{\Omega }}_i}}+{\displaystyle \frac{1}{2}}{\mathrm{\Pi }}_{\mathit{ij}}){\mathit{v}}_{\mathit{ij}}·{\mathbf{\nabla }}_{i}W(h_i)\end{array},\end{array}$$(8)

where the index i refers to a generic ith particle and the index j in the sums refers to the jth particle that is enclosed within the range 2h_M = 2 ⋅ max(h_i, h_j). The term g_soft represents the softened gravitational force per unit mass mentioned above: it is a function only of the mutual particle distance r_ij and of the smoothing length h. It tends to zero as r_ij → 0 and assumes the classical Newtonian form $m_{j}G{r}_{\mathit{ij}}^{-2}$ for r_ij ≥ 2h. The operator ∇_i represents the gradient with respect to the coordinates of the ith particle. The gradient is performed over two different expressions of the kernel W, with two different lengths h_i and h_j. The terms ζ_i and Ω_i are suitable functions that account for the variation in the smoothed Newtonian potential with respect to the softening length and for the non-uniformity of the softening length itself, respectively. They assume this form for a generic particle of index i:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_i=1+{\displaystyle \frac{h_i}{3{\rho }_i}}\sum _j{m}_j{\mathbf{\nabla }}_{\mathit{i}}W(r_{\mathit{ij}},h_i)\end{array}$$(9)

$$\begin{array}{c}\hfill {\zeta }_i=-{\displaystyle \frac{h_i}{3{\rho }_i}}\sum _j{m}_j{\displaystyle \frac{\partial {\varphi }_{\mathrm{soft}}(r_{\mathit{ij}},h_i)}{\partial h_i}},\end{array}$$(10)

where in the same way for the system in Eq. (6), the sum extends over the particles that are enclosed within the range 2h_i. In Eq. (10), the function ϕ_soft represents the softenend gravitational potential, such that ∇ϕ_soft = −g_softr_ij/r_ij. The potential reaches a constant value as r_ij → 0 and becomes equal to the Newtonian potential for r_ij ≥ 2h (for an explicit expression, see, e.g., Price & Monaghan 2007). Terms Ω and ζ are computed in the same neighbor-search iterative loop where ρ and h are worked out.

Only if the gas interacts with stars does the last term in Eq. (7), $\mathrm{d}{\mathit{v}}_i^{[\mathrm{stars}]}/\mathrm{d}t$ (discussed in Sect. 2.2.4), represent a non-null acceleration that accounts for the Newtonian interaction between ith particle and the point masses. The function Π_ij, which we discuss in the following section, characterizes the well-known artificial viscosity. The expression of the equation of motion, Eq. (7), guarantees a symmetric exchange of linear momentum between the particles.

2.2.3. Artificial viscosity

In high-compression regions such as shock wavefronts, the velocity gradient may be so strong that two layers interpenetrate and the hydrodynamical equations may not be integrated correctly and generate unphysical effects. Additional artificial pressure terms are a possible solution for this problem. In our code, we added an artificial term by adopting the same classical schematization as Monaghan (1989), which corresponds to introducing a suitable artificial viscosity that dampens the velocity gradient when two particles approach. Practically, a viscous-pressure term Π_ij is included in Eqs. (7) and (8). It assumes the expression

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\mathit{ij}}=\{\begin{array}{cc}{\displaystyle \frac{-\alpha \overline{c}{\mu }_{\mathit{ij}}+\beta {\mu }_{\mathit{ij}}^2}{{\overline{\rho }}_{\mathit{ij}}},}\hfill & \text{if}{\mathit{v}}_{\mathit{ij}}·{\mathit{r}}_{\mathit{ij}}<0,\hfill \\ 0,\hfill & \text{if}{\mathit{v}}_{\mathit{ij}}·{\mathit{r}}_{\mathit{ij}}>0,\hfill \end{array}\end{array}$$(11)

where ${\mu }_{\mathit{ij}}={\displaystyle \frac{h{\mathit{v}}_{\mathit{ij}}·{\mathit{r}}_{\mathit{ij}}}{r_{\mathit{ij}}^2+{\eta }^2{\overline{h}}^2}}$. The dot product v_ij ⋅ r_ij involves the relative velocity and the distance of a pair of particles i − j. Only the particles that move in, for which v_ij ⋅ r_ij <  0, contribute to the artificial viscosity. The parameter η is a suitable term to prevent singularities when two particles approach very closely (we use the typical value of η = 0.1). The terms $\overline{h},\overline{\rho ,}$ and $\overline{c}$ represent the average values of the smoothing length $\frac{1}{2}(h_i+h_j)$, the density $\frac{1}{2}({\rho }_i+{\rho }_j),$ and the speed of sound $\frac{1}{2}(c_{\mathrm{s}i}+c_{\mathrm{s}j})$, respectively. We set β = 2α. In this simple formulation, the artificial viscosity is activated throughout the fluid; but in two circumstances it should be damped to prevent unphysical effects. Artificial viscosity must be damped in regions where shear dominates, and where the velocity gradient is low.

For two shearing layers of fluid, the relative velocity between the particles leads to an approach that is interpreted by the artificial viscosity (11) as a compression. This incorrect interpretation causes the code to overestimate the strength of the viscous interaction. To prevent false compressions, Balsara (1995) multiplied the term μ_ij by a proper switching coefficient:

$$\begin{array}{c}\hfill f=\frac{|\mathbf{\nabla }·\mathit{v}|}{|\mathbf{\nabla }·\mathit{v}|+|\mathbf{\nabla }\times \mathit{v}|+{10}^{-4}{c}_{\mathrm{s}}{h}^{-1}},\end{array}$$(12)

in which the divergence of velocity and the velocity curl are evaluated for a particle of index i as

$$\begin{array}{c}\hfill \{\begin{array}{c}{(\mathbf{\nabla }·\mathit{v})}_i={\rho }_i^{-1}\sum _j{m}_j{\mathit{v}}_{\mathit{ij}}·{\mathbf{\nabla }}_{i}W(r,h_i)\hfill \\ {(\mathbf{\nabla }\times \mathit{v})}_i={\rho }_i^{-1}\sum _j{m}_j{\mathit{v}}_{\mathit{ij}}\times {\mathbf{\nabla }}_{i}W(r,h_i)\hfill \end{array}.\end{array}$$(13)

We implemented the term f by multiplying μ_ij for an average value $\overline{f}=\frac{1}{2}(f_i+f_j)$. Further problems may arise far away from high-compression regions. In the classical formulation of Π_ij, α = 1 = cost. (e.g., in Monaghan 1992). In this scheme, the viscosity acts with the same effectiveness in every region, while we would expect the artificial term to be efficient only where it is needed, that is, close to the shock fronts. To solve this problem, we used the same formalism as was introduced by Morris & Monaghan (1997) and further developed by Rosswog et al. (2000) by considering an individual α_i for each particle that follows the time-variation equation,

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{\alpha }_i}{\mathrm{d}t}}=-{\displaystyle \frac{({\alpha }_i-{\alpha }_{\min })}{{\tau }_{\alpha }}}+S_i,\end{array}$$(14)

where S_i = max(−(∇ ⋅ v)_i, 0) (α_max − α_min) represents a “source” term that increases in the proximity of the shock front; α_min represents a minimum threshold value for α, and α_max represents its maximum. The (increasing) rate of the viscosity coefficient is driven by a characteristic timescale τ_α = h_i/bc_s that depends on how the fluid allows the perturbations to propagate through the resolution length. The individual viscosity coefficients α_i and α_j, when referred to a generic i − j particle pairing, are averaged in the same way as was done with the other quantities.

For a gas with γ = 5/3, a good value for the b coefficient can be set such that 5 ≤ b⁻¹ ≤ 10 (Morris & Monaghan 1997). For our tests, we set α_max = 2, α_min = 0.1 and b⁻¹ = 5. These are the most commonly adopted values in literature for a wide class of problems involving collapse, merging stars or protoplanetary disks (see, e.g., Rosswog & Price 2007; Stamatellos et al. 2011; Hosono et al. 2016). The implementation of the artificial viscosity term (Eqs. (12)–(14)), together with its form implemented in Eqs. (7) and (8), may affect the accuracy of the code in preserving the total angular momentum. In Sect. 3.1.5 we discuss how this form of viscosity, with different choices of the coefficients α_min and b, guarantees the conservation of the angular momentum.

2.2.4. Additional stellar objects

We calculated a direct point-to-point interaction for the mutual interaction between stars and to couple stars with SPH particles. The equation of motion of a generic p-th star takes the following form:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{\mathit{v}}_{\mathrm{p}}}{\mathrm{d}t}}=-\sum _j\frac{1}{2}(g_{\mathrm{soft}}(r_{\mathrm{p}j},{ϵ}_{\mathrm{p}})+g_{\mathrm{soft}}(r_{\mathrm{p}j},h_j))\frac{{\mathit{r}}_{\mathrm{p}j}}{r_{\mathrm{p}j}}\\ \hfill -\sum _{\mathrm{s}}\frac{1}{2}(g_{\mathrm{soft}}(r_{\mathrm{ps}},{ϵ}_{\mathrm{p}})+g_{\mathrm{soft}}(r_{\mathrm{ps}},{ϵ}_{\mathrm{s}}))\frac{{\mathit{r}}_{\mathrm{ps}}}{r_{\mathrm{ps}}},\end{array}\end{array}$$(15)

where g_soft(r, ϵ) represents the Newtonian acceleration, which takes the form discussed above in Sect. 2.2.2. The force softening is accounted for the stars as well according to a constant softening length ϵ_s = cost. The gravity is thus softened when the mutual distance approaches ϵ_s. The first summation is extended over all SPH particles, while the index s in the second sum refers to the generic stars.

Similarly, the equation of motion, Eq. (7), when referred to a gas ith particle contains the following sum:

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{\mathit{v}}_i^{[\mathrm{stars}]}}{\mathrm{d}t}}=-{\displaystyle \sum _s\frac{1}{2}(g_{\mathrm{soft}}(r_{i\mathrm{s}},{ϵ}_{\mathrm{s}})+g_{\mathrm{soft}}(r_{i\mathrm{s}},h_i))\frac{{\mathit{r}}_{i\mathrm{s}}}{r_{i\mathrm{s}}},}\end{array}$$(16)

where the index s again refers to the stars, and r_is is the distance vector between a gas particle and a star.

2.2.5. Time integration and time-stepping

To evolve the gas system in time, we adopted a second-order integration method that is similar to a classical second-order Runge–Kutta scheme but is at the same time very similar to a leap-frog integrator: the well-known velocity-Verlet method (see Andersen 1983; Allen & Tildesley 1989, for detailed references). The Verlet method is based on a trapezoidal scheme coupled with a predictor-corrector technique for the estimation of v and u. The structure of this scheme is very similar to that of classical symplectic leap-frog algorithms, although it requires two computations of the force at every time iteration (see Appendix A.1). Nevertheless, the general velocity-Verlet method applied to gas evolution shows some advantages compared to the symplectic algorithm of the same order. Like a standard Runge–Kutta method, velocity and positions are updated in synchronized steps, without the Δt/2 shift. This feature provides a good flexibility in problems that are approached with nonuniform time-steps that involve the interaction of the gas component with other components that are integrated with different methods, as in our case. Various applications of velocity-Verlet methods in SPH schemes are found in the literature, for example, in Hubber et al. (2013) or in Hosono et al. (2016).

The additional point masses are ballistic elements, whose equation of motion needs to be integrated with very high precision to avoid secular trends that are typical of few-body gravitational problems. Although the SPH precision is at only second order, we decided to integrate the Newtonian motion of the (few) stars and planets in the system with a 14th-order Runge–Kutta method that was recently developed by Feagin (2012) through the so-called m-symmetry formalism. The method consists of 35 force computations per time-step, and in analogy with the well-known second- and fourth-order RK methods, it updates the velocities and the positions by suitable linear combinations of 35 different Kr and Kv coefficients (see Appendix A.2 for further details).

For the gas, we chose the time-step Δt following a criterion similar to the standard Courant–Friedrichs–Lewy (CFL) criterion that is commonly adopted for SPH systems (see, e.g., Monaghan 1992), together with some additional criteria. A global time-step Δt_min can be determined by taking the minimum between the following two quantities:

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{term}}=\underset{i}{min}({\displaystyle \frac{Ch}{c_{\mathrm{s}i}+h_i{|\mathbf{\nabla }·\mathit{v}|}_i+\phi {\alpha }_i[c_{\mathrm{s}i}+2\underset{j}{max}({\mu }_{\mathit{ij}})]}},C_u{\displaystyle \frac{u_i}{{\dot{u}}_i}})\end{array}$$(17)

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{dyn}}=\underset{i}{min}(C_a\sqrt{{\displaystyle \frac{h_i}{a_i}}},C_d{\displaystyle \frac{v_i}{a_i}}),\end{array}$$(18)

where c_s is the sound speed, C is a coefficient whose typical value lies between 0.1 and 0.4; we usually choose 0.15. Moreover, C_u, C_a, and C_d are coefficients to be set < 1. We chose C_u = 0.04, C_a = 0.15, and C_d = 0.02. Finally, the coefficient φ typically ranges from 0.6 to 1.2 (we adopt φ = 1.2 throughout). Similarly to the control of the variation in kinetic energy, we control the time variation of the thermal energy $u/\dot{u}$ in a single time-step by limiting it to a certain fraction C_u = 0.04. The index i refers to an individual time-step Δt_i that is related to a specific particle.

For the point-particle phase in the system (i.e., stars or planets), we chose a characteristic time-step, Δt, that is defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{stars}}=\underset{\mathrm{s}}{min}(C_{\mathrm{ob}}\sqrt{{\displaystyle \frac{{ϵ}_{\mathrm{s}}}{a}}},C_d{\displaystyle \frac{v_{\mathrm{s}}}{a_{\mathrm{s}}}}),\end{array}$$(19)

where we use C_ob = 0.15. The various quantities with the index s of course characterize a specific star particle.

For a homogeneous medium, the integration can be performed with a global time-step, that is, the lowest value of gas and stars. Generally, the particles have different resolutions h_i and different accelerations, which leads to a wide class of typical evolution timescales. For some particles, the integration might therefore be made with different Δt_i, which avoids the explicit force calculation at every time iteration and saves some computing time. We adopted a technique that was implemented in several N-body algorithms, such as the classical TREESPH (Hernquist & Katz 1989) or in the multi-GPU-parallelized N-body code HiGPUs (Capuzzo-Dolcetta et al. 2013). We assigned to each point a time-step as a negative two-power fraction of a reference time $\mathrm{\Delta }{t}_{\max }=\underset{i}{max}(\mathrm{\Delta }{t}_i)$ (it can be a fixed quantity, or it may change periodically during the simulation). The particle motion is updated periodically according to their Δt_i in such a way that after an integration time Δt_max, all of them are synchronized (further details are explained in Appendix A.1). Particle mapping and sorting are performed every time for every particle as well, independently of their individual time-step. Thus, the configuration of the tree grid, together with total mass and quadrupole momentum of the boxes, are computed every single step Δt_min. Similarly, at each minimum time-step iteration, the gravitational interactions between gas and star are computed even during “inactive” stages of the gas. Thus, the acceleration of the inactive SPH particles is split into two terms: one is given by a fixed non-updated hydrodynamical and self-gravitaty term, and the other is given by a constantly updated gas-star gravitational force.

In our scheme, the stars and planets do not follow an individual time-step scheme, and their mutual interactions are computed for every single step Δt_min, even when Δt_stars ≠ Δt_min. Furthermore, we force the particles that lie close to the stars within a tolerance distance to be integrated at every time iteration. Practically, we compute the distance for a generic ith particle from the stars, and we furthermore predict this distance at the following time iteration. When these values are lower than a tolerance of κϵ_ob (with a constant κ ≥ 2), the particle time-step drops to Δt_min. For our practical purposes, we used a small number of objects (in the current investigation, N_ob ≤ 2), thus, the 35-stage RK scheme requires a relatively short CPU-time (less than 2% of the total).

In gas problems that involve strong shocks, the use of individual time-steps may lead to strong errors. Even though CFL conditions are satisfied, the strong velocity gradients may determine a great discrepancy in time-steps between close particles. Consequently, close particles may evolve with very strongly different timescales. This may create too many asymmetries in the mutual hydrodynamical interactions, which causes unphysical discontinuities in velocity and pressure. Following the idea of Saitoh & Makino (2009), we limit for each pair of neighboring ith and jth particles the ratio of time-steps ${\displaystyle \frac{\mathrm{\Delta }{t}_i}{\mathrm{\Delta }{t}_j}}\le A$. These investigators have shown that a good compromise is reached by the choice of A = 4, which gives good results without abruptly affecting the efficiency of the code.

2.2.6. Approximation of the gravitational field: opening criterion

The decomposition of the system into a series of clusters is performed by the tree algorithm through a recursive octal cube eight sub-boxes, each one subdivided into a further eight cubes of order L = 2, and so on. A tree structure is thus constituted, made of several nested boxes, each of which contains a group of particles. To calculate the acceleration of an ith particle, the algorithm walks along the tree, starting from the low-order cubes toward the highest order cubes (that contain just one particle), and evaluates the distance between the particle and the center of mass of the boxes. Each time a box is probed, the code decides to open it and probes its internal cubes only if the well-known opening criterion is satisfied,

$$\begin{array}{c}\hfill {\displaystyle \frac{D_{\mathrm{L}}}{r}}>\theta ,\end{array}$$(20)

where θ is the so-called opening angle parameter for which reasonable values range among 0.3 and 1 (see Sect. 3.2, dedicated to performance tests), and D_L = D₀ × 2^−L is the side length of the box. In the opposite case, the algorithm decides to approximate the gravitational field by adding for the acceleration of the ith particle just the contribution of the box (given by Eq. (3)). With this scheme, the net amount of computation scales down to N log N, which is far shorter than N² for large N.

With r_a the geometrical center of a certain cube and r_CM its center of mass position, we may find a very large offset Δ_CM ≡ |r_a − r_CM| under specific circumstances. With a center of mass far away from the box center, some errors may arise in the force approximation because a cube can be considered “far enough” from a particle according to the opening criterion even though some of the points enclosed in the box may still be very close to the particle (see Fig. 1). These close particles are therefore ignored, and the whole box gives the multi-polar approximated contribution to the particle acceleration. The acceleration is therefore calculated with less accuracy than might be expected. A first key to avoid these errors should be adopted by checking whether the particle lies very close to a box (as was done, e.g., by Springel 2005). If the test particle is inside a cube or close to its borders according to a certain tolerance, the box is always opened, independently of the truthfulness of Eq. (20).

Fig. 1. Schematic 2D example of the lack of accuracy in the field computation caused by a large offset ΔCM. The center of gravity is far away from the ith particle, and the cube is not opened. Nevertheless, some particles that lie near the edge of the box, such as the jth and the kth points, are very close to the ith point, but their direct contribution is missed, which results in a loss of accuracy.

We can furthermore optionally modify the opening criterion in our code by taking the offset term into account. We may use the following rule to open a box:

$$\begin{array}{c}\hfill r<{\displaystyle \frac{D_{\mathrm{L}}}{\theta }}+{\mathrm{\Delta }}_{\mathrm{CM}}.\end{array}$$(21)

This prescription is equivalent to the classic opening criterion, but with an effective opening angle θ′< θ, to guarantee that every close box is opened. In some peculiar cases in which Δ_CM is large (i.e., comparable with the length of the semidiagonal of the box), as in the example of Fig. 1, the effective opening angle is considerably smaller than θ. In Sect. 3.2 we show that for a typical value of θ = 0.6, the adoption of the new criterion does not require too much additional computational effort, especially when many particles are involved.

3. Code testing

We illustrate here some basic physics tests (Sect. 3.1) and a series of performance tests (Sect. 3.2). In Sect. 3.1 we apply GaSPH to two basic problems: (i) a non-hydrodynamical system, characterized by a cluster of point-mass particles distributed according to a Plummer profile, and (ii) a classical shock-wave problem. These quality tests are followed by some applications to hydrodynamical systems at equilibrium. First, we treat some polytropes with finite radius. Then, we compare our algorithm with a well-known hydrodynamical tree-based code (Gadget-2) in the case of a gaseous Plummer sphere. In Sect. 3.2 we analyze the computational efficiency and the accuracy of our code in different contexts.

3.1. Tests with gas and pressureless systems

3.1.1. Turning off the SPH: the evolution of a pressureless system

In order to test the stability of our numerical method, we performed a series of simulations in which we placed a set of points according to the standard Plummer configuration (Plummer 1911), which is often adopted to study the distribution of stars in globular clusters. The Plummer sphere is pressureless, so that the particles interact only though gravity, and there is no SPH interaction. As units of measurement, we chose the total mass M and the gravitational constant G, and we placed an ensemble of N = 10⁵ particles in a Plummer distribution with core radius R = 1 and cutoff radius R_out = 10R. The particles had equal masses m = N⁻¹ and equal softening length ϵ, chosen as a fraction of the central mean interparticle distance: $ϵ={\alpha }_{\mathrm{s}}{\left(\frac{m}{{\rho }_0}\right)}^{1/3}={\alpha }_{\mathrm{s}}{\left(\frac{4\pi }{3N}\right)}^{1/3}$ (with α_s ∈ [0.2, 1.0]). Starting the Plummer distribution at the virial equilibrium, we integrated its time evolution for 50 mean crossing-times τ_c. This parameter is defined as the initial ratio between the half-mass radius and the mean dispersion velocity $\frac{R_{1/2}}{\sqrt{<v^2>}}$. Figure 2 shows the virial ratio $\frac{2T}{|\mathrm{\Omega }|}$ as a function of time and compares four runs made with different combinations of the opening angle θ (0.6; 1.0) and ϵ (0.2; 0.5). The four results illustrated in Fig. 2 do not show any relevant difference: the virial ratios oscillate within a small fraction < 0.5%, especially for the configuration with θ = 1 and α_s = 0.5, which was expected to be the worst case.

Fig. 2. Oscillation of the virial ratio as a function of time for different choices of code configuration parameters for the simulation of a Plummer distribution of 105 equal-mass particles. The continuous line shows θ = 0.6, ϵ = 0.2; the dashed line represents θ = 1.0, ϵ = 0.2; the dotted line shows θ = 0.6, ϵ = 0.5; and the dash-dotted line denotes θ = 1.0, ϵ = 0.5.

We note that we do not work with a classical high-precision N-body code such as Nbody-6 (for example see Aarseth 1999). The Newtonian force is approximated by means of both the multipolar expansion, which occurs when particles are sufficiently far, and the softening length damping, which occurs when the particles approach within a distance of about ϵ. Despite these approximations, acceptable results can be obtained in a noncollisional system like our Plummer distribution. The results lie within reasonable errors.

3.1.2. Sedov–Taylor blast wave

To test the code with strong shock waves, we simulated the effects of a point explosion on a homogeneous infinite hydrodynamical medium with constant density ρ₀ and null pressure. If an amount of energy E₀ is injected at a certain point r₀, an explosion occurs and then a radial symmetric shock wave propagates outward. Sedov (1959) investigated this problem and found a simple analytical law for the time evolution of the shock front:

$$\begin{array}{c}\hfill r_{\mathrm{s}}(t)={\left({\displaystyle \frac{E_0}{{\rho }_{0}a}}\right)}^{1/5}{t}^{2/5},\end{array}$$(22)

where r_s is the radial position of the front relative to the point of the explosion r₀, while a is a function of the adiabatic constant γ (it is close to 0.5 for γ = 5/3, and it approaches 1 for γ = 7/5). Furthermore, the fluid density immediately behind the shock front (r ≤ r_s) has the following radial profile:

$$\begin{array}{c}\hfill \rho (r,t)={\displaystyle \frac{\gamma +1}{\gamma -1}}{\rho }_0{G}_{\gamma }\left({\displaystyle \frac{r}{r_{\mathrm{s}}}}\right),\end{array}$$(23)

where G_γ is an analytical function of the relative radial coordinate r/r_s. Similarly as in many previous works (see, e.g., Rosswog & Price 2007; Tasker et al. 2008), we set the initial conditions for a homogeneous and static medium (ρ₀ = 1, v = 0) by placing 10⁶ equal-mass particles in a cubic lattice structure, confined in a box with x, y, and z coordinates each ranging from −1 to 1. γ was set to 5/3, and the explosion was simulated by giving an amount of energy E₀ = 1 to the origin of the system. We were unable to reproduce a point explosion with an SPH system because its spatial resolution is determined by the kernel support. We therefore needed to inject the energy in a small region with the same scale as 2h. We thus gave at a time t_* the energy E₀ to the particles that were enclosed in a sphere with radius R = 2h. Figure 3 shows three different radial average density profiles ρ(r, t′) that correspond to the times t = 0.05, t = 0.1, and t = 0.2. The results are compared with the analytical solution. Although the position of the front follows the expected law of Eq. (22), the peak does not reach the expected value $\frac{\gamma +1}{\gamma -1}{\rho }_0=4{\rho }_0$.

Fig. 3. Radial density profiles of the Sedov–Taylor blast wave at several times (increasing rightward) t = 0.05 (circles), t = 0.1 (triangles), and t = 0.2 (squares) in a simulation with N = 106. The results obtained using a higher resolution (N = 3 375 000) are plotted with dotted lines. The full lines represent the classical Sedov–Taylor auto-similar solutions.

Intrinsic errors in approximating the physical quantities, given by the smoothing kernel, allow the density to spread out and follow a wider distribution than the true profile. This corresponds to a smoothing of the vertical discontinuity and so to a lower peak of the density. The same figure shows a comparison with results obtained from a further test that was made with the same system, but using a better resolution (N = 3 375 000). The peak of the curve clearly reaches a higher value.

3.1.3. Polytropes at equilibrium

We tested our code in the case of hydrodynamic self-gravitational systems by building static polytropes with different indexes (n = 1, n = 3/2, and n = 2). A generic polytrope of index n constitutes a radially symmetric system whose equation of state follows the expression

$$\begin{array}{c}\hfill P(r)=K_n\rho {(r)}^{1+\frac{1}{n}},\end{array}$$(24)

where the density is parametrized as ρ(r)/ρ₀ = θⁿ(r), and ρ₀ is the central value. The static radial solution θ(r) can be found by writing an equilibrium condition between the hydrostatic pressure gradient and the gravitational forces, from which the well-known Lane–Emden equation can be obtained (an exhaustive treatment can be found, e.g., in Chandrasekhar 1958):

$$\begin{array}{c}\hfill {\displaystyle \frac{{\alpha }^2}{r^2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(r^2{\displaystyle \frac{\mathrm{d}\theta }{\mathrm{d}r}}\right)=-{\theta }^n,\end{array}$$(25)

with ${\alpha }^2=(n+1)K_n{\rho }_0^{\frac{1}{n}-1}/4\pi G$, and K_n a suitable normalization coefficient. For an index n ∈ (0, 5), the system has a finite radius and the coefficient K_n depends, through α², both on the radius R and on the total mass M. We set both of them to 1 in our tests, implying K₁ ≈ 0.637, K_3/2 ≈ 0.424 and K₂ ≈ 0.365.

We tested the ability of our code to let a system spontaneously relax in a polytrope configuration, following the prescription adopted in Price & Monaghan (2007). Starting from a homogeneous sphere of particles placed in a lattice structure, the system was let evolve by forcing the pressure to follow Eq. (24). We forced the SPH system to evolve by damping the velocities with an additional acceleration a_damp = −0.05 v, until the kinetic energy decreased to a small fraction (1%) of the total energy. A standard nonconstant α was chosen for the artificial SPH viscosity, with α₀ = 0.1, and a number of neighbors of 110 was set for the particles.

A correct treatment of self-gravity and hydrodynamic interactions among SPH particles, and the choice of the equation of state (24) allows the system to acquire the density profile ρ₀θⁿ, which is the solution of Eq. (25). Figure 4 shows the three radial density profiles we obtained for the different polytropic indexes. The resolution, related to the particles number, affects the accuracy of the code in correctly sampling the profile ρ(r), especially in the central denser regions. Mainly for the higher index n = 2, a higher particle number is needed to let the numerical density approach the theoretical expected value at a specific accuracy level. 10 000 particles were used for the models with n = 1 and n = 3/2, while the polytrope with index n = 2 was built with 20 000 particles.

Fig. 4. Equilibrium analytical solutions of the density profiles ρ(r) related to three different models of polytropes with indexes n = 1, n = 3/2, and n = 2, drawn as solid curves. Comparison with the computed profiles (dots).

3.1.4. Gaseous Plummer distribution

We tested the equilibrium of a static gas density distribution according to the Plummer function:

$$\begin{array}{c}\hfill \rho (r)={\rho }_0{[1+{\left({\displaystyle \frac{r}{a}}\right)}^2]}^{-5/2},\end{array}$$(26)

with the central density defined by ρ₀ = 3M/4πa³, where M and a are the total mass of the system and a characteristic length, respectively. We set M = 1 and a = 1 (G = 1 in internal code units) so that the half-mass radius of the system is r(50%) ≈ 1.3. In a static configuration with a null velocity field, the gas SPH particles compensate for the mutual self-gravity with a pressure gradient that results from a temperature distribution T(r)=κρ(r)^1/5. κ represents a constant that is calibrated by taking the equation state of a perfect gas P = (γ − 1)ρu into account and by imposing the virial equilibrium between gravitational energy W and total thermal energy $U=\sum _i{u}_i{m}_i$, that is, |W|=2U. We placed 50 000 particles according to a Monte Carlo sampling of the distribution in Eq. (26). A realistic distribution has an infinite radius, therefore we used a cutoff at a proper radial distance r ≈ 22, such that the distribution contained 99.8 % of the mass of a realistic infinitely extended Plummer sphere. Figure 5 shows the time variation of some Lagrangian radii that contain 5%, 15%, 30%, 50%, 75%, and 90% of the total system mass for an integration time of 90 central free-fall timescales ( τ₀ = (3π / 32Gρ₀)^1/2 ≈ 1 in our code units).

Fig. 5. Lagrangian radii as a function of time for a hydro-Plummer distribution at equilibrium (see Sect. 3.1.4). The radii are normalized to their respective initial values.

We compared the results with the well-known gravitational SPH code Gadget-2 (Springel 2005). Figures 6 and 7 show a comparison between the radial density profiles obtained by the two algorithms with the same choice of the main parameters. The α viscosity coefficient was set constant and equal to 1. The reported density was computed at t = 90 τ₀, even though the system reached an acceptable equilibrium state already within a few units of τ₀ after several slight oscillations.

Fig. 6. Radial density profile for a Plummer distribution with 50 000 SPH particles (solid line). The results obtained with Gadget-2 are also shown (dashed line). The analytical Plummer profile is plotted with a dotted line. The density is in units such that ρ0 = 3/(4π).

Fig. 7. Logarithmic ratio ρN/ρA of the numerical to analytical radial density profile for a Plummer distribution with 50 000 SPH particles (solid line). The results obtained with Gadget-2 are also shown (dashed line).

In Fig. 7 we can distinguish three main radial zones for r ≤ 0.3, 0.3 <  r <  2, and r ≥ 2. In the middle zone, the codes are in good agreement and provide a density profile with an accuracy lower than 2% with respect to the analytical model. For r <  0.3, Gadget-2 describes a density that deviates by up to the 6% from the expected value, while our program has a maximum deviation of 11%. For both models, these higher errors can be ascribed to the fact that the system contains only about 2% of the total mass (and thus 2% of the total particles) within the radial distance r = 0.3, which causes a poor sampling of the potential inside the sphere. Consequently, the system tends to shrink slightly. In the outward zone, the deviations can reach significantly higher values with both codes because the density values are far lower than in the central zone. We therefore conclude that in the context of a standard physical environment, the two codes show a satisfactory agreement overall.

3.1.5. Artificial viscosity and angular momentum conservation

A nonconstant artificial viscosity may lead to nonconservation of the angular momentum, L. The actual conservation of this quantity was tested with different settings of the artificial viscosity parameters in Eq. (14), integrating the time evolution of a system that was very similar to the one described in the previous section. We used the same Plummer distribution (see Eq. (26)), with M = 1 and a = 1, made of 50 000 SPH particles. The same thermal energy profile was adopted, but scaled down by a factor 1/2, so that $T(r)=\frac{\kappa }{2}\rho {(r)}^{1/5}$. We converted the (subtracted) thermal energy into kinetic energy by assigning to each ith particle a clockwise azimuthal velocity, with absolute value ${\mathit{v}}_i=\sqrt{u_i}$ (where u_i was the original specific thermal energy characterizing the Plummer system used in Sect. 3.1.4), and a direction parallel to the X, Y plane. The system thus acquired a nonzero vertical component of the angular momentum, $L_z=\sum _i{m}_i(x_i{\mathit{v}}_{\mathit{y}i}-{\mathit{y}}_i{\mathit{v}}_{\mathit{xi}})$. The virial equilibrium was still formally preserved because gravitational potential energy and thermokinetic energy were the same such that |W|=2(K + U), but the (new) angular rotation triggered changes in the density distribution. We integrated in time for about 100 initial central free-fall timescales τ₀ (which is on the same order as the azimuthal dynamical timescale, taken as the ratio r_c/v(r_c)≈1.4, where r_c is the initial radius at which the density drops by a factor 1/2). We performed three different simulations by varying in the α rate Eq. (14) the parameters α_min and b. We set α_min = 0.1,  b⁻¹  =  5, α_min  =  0.02,  b⁻¹  =  5, and α_min  =  0.1,  b⁻¹  =  7. The angular rotation changes the configuration of the system, which causes the initial Plummer density distribution to become flatter perpendicularly to the z-axis, while the whole system expands. During an initial phase of about 20τ₀, the distribution underwent some rapid variation followed by a slow secular evolution. Figure 8 shows the quantity (L_z−L_z0)/|L_z0|, which represents the variation as a function of time of the component L_z compared to its initial value L_z0. The three lines refer to the different choices of the parameters α_min and b. The curves show a conservation of the angular momentum within 10⁻³ up to 100 evolution timescales. In particular, the choice of α_min = 0.1,  b⁻¹ = 7, compared with the other configurations, gives a stronger variation of L_z during the first phases, while it shows a lower change rate during the secular evolution of the system. On the other hand, a low value α_min = 0.02 gives rise to a better conservation in the initial phases and a higher deviation during later stages. We observed in all three simulations that the two components L_x and L_y maintain values compared to L_z within a relative error of 10⁻³.

Fig. 8. Fractional variation in the z-component of the angular momentum for our simulated rotating Plummer model (see Sect. 3.1.5). We plot the results obtained by varying some configuration parameters: αmin = 0.1,  b−1 = 5 (full line), αmin = 0.02,  b−1 = 5 (dotted line), and αmin = 0.1,  b−1 = 7 (dash-dotted line).

3.2. Code performance

In order to analyze the computational efficiency of our algorithm as a function of the particle number, we performed several tests by measuring the average CPU time that was spent for a single run by the main routines. We therefore studied the performances of the GaSPH code in three different contexts:

• (1)

  A system with pure self-gravity and zero pressure, adding a comparison with the results of Gadget-2.

• (2)

  A system with self-gravity and SPH pressure.

• (3)

  A system similar to that of case 2, but with the addition of 20 point star-like external objects.

We performed the tests by placing a set of N particles with the same Plummer density profile distribution as adopted in Sect. 3.1.1. The program was tested on an Intel®Core^TM i7-4710HQ architecture with 6MB of cache memory and with 16GB of RAM memory DDR3L with a data-transferring speed of 1600 MHz.

For a standard tree code without SPH, the computational time per particle is expected to be linear in log N because the overall time scales as N log N. To increase the efficiency and save considerable memory resources, we can also use a simple formalism made by considering only the first “monopole” term −MGr⁻³r that appears in the right member of Eq. (3). This is a technique that was also adopted in Gadget-2 and simplifies the complexity of the algorithm by neglecting the efforts for the quadrupole tensor computation. The suppression of the quadrupole term decreases the computational time, with a minor cost in terms of accuracy. For a pure self-gravitating system, Fig. 9 shows the CPU time needed for a single particle force calculation as a function of the ten-based logarithm of the particle number N (ranging from 10⁴ to 5 × 10⁶). Choosing an opening angle θ = 0.6, we performed a series of force evaluation by considering a simple pressureless system, with particles interacting only with the Newtonian field. The computational times measured using our code (averaged over a reasonable number ≥30 of equal tests) are comparable with the average CPU times measured using Gadget-2. The figure also shows the results based on a second series of runs with GaSPH performed with the quadrupole term included in the gravitational field. Including this term, an additional CPU time of about 30% is requested. Figure 10 compares the previous CPU times with the times needed by GaSPH for a full self-gravitating SPH system, with the quadrupole term included in the computation and the same value of θ = 0.6. The times per particle for the density computation routine are also shown. In computing the acceleration, the additional time per particle is fairly independent of N, as can be seen in the figure, because the close SPH interactions are always made over a fixed number of neighboring points, which we set to 60 in this example. For the same reasons, the average time per particle needed to calculate the density is also expected to be constant, as Fig. 10 shows. The calculus of ρ and h requests an iterative process in which the routine for each particle is called several times. The CPU times illustrated in the figure are the average values per single iteration. Determining the optimum value of h the code typically requires no more than two iterations.

Fig. 9. CPU time per particle for the pure gravitational force calculation in monopole and quardupole approximations at different N. Gadget-2 results (empty circles) are compared to our code results (squares connected with dashed lines) in the same monopole approximation. The continuous line refers to the performance of GaSPH, and the quadrupole term is included in the field.

Fig. 10. GaSPH average CPU times per particle as a function of N. Comparison of different routines: SPH neighbor-search routine averaged for a single iteration (circles), pure gravity computation with monopole term (empty squares) and with quadrupole term (filled squares), self-gravity computation up to quadrupole and including the SPH terms (triangles). The units are the same as in Fig. 9.

The optional introduction of the offset Δ_CM term in the opening criterion (as discussed in Sect. 2.2.6) causes in some cases a considerable reduction of the effective angle θ. Consequently, the number of direct particle-to-particle interactions increases, which decreases the code performance. Figure 11 illustrates the code efficiency in terms of number of particles processed in a second. The results, related to the two acceleration routines (pure self-gravity and self-gravity with SPH) shown in the previous graph, are compared with other result obtained by including the offset term Δ_CM in the opening criterion of Eq. (21). A substantial but not drastic worsening in performance can be observed. For instance, using 5 × 10⁶ SPH particles and including Δ_CM, the code computes the accelerations at a rate of ≈24 000 particles per second (about 17% slower than the case without Δ_CM). Computations were made with θ = 0.6 and including the quadrupole terms. The performance of the same two force subroutines (pure gravity and gravity plus hydrodynamics) were also studied at different values of θ (CPU times per particle as a function of N are shown in Fig. 12).

Fig. 11. Number of processed particles per second. The continuous lines and the dashed lines refer to results with and without the correction term ΔCM, respectively. We show the simple gravity field calculation (empty circles) and the full self-gravity routine with hydrodynamics (squares). The quadrupole term is considered for the gravity field. θ = 0.6.

Fig. 12. Computational time per particle vs. log N for various values of the opening angle θ. The pure tree gravitational algorithm (dashed line) and full SPH+gravity algorithm are shown as dashed and solid lines, respectively. The quadrupole term is included in the force evaluation. The ΔrCM offset term is not considered by the opening criterion.

Smaller angles should provide higher precision at the cost of a longer computational time. On the other hand, with larger angles there are fewer direct point-to-point interactions and we gain in efficiency, but we expect a lower accuracy. We evaluated the accuracy of our tree code by measuring a “mean relative error” in computing the accelerations according to the prescription suggested by Hernquist (1987), with different conditions of particle number N and opening angle θ. The prescription consists of a comparison of the three components of the acceleration vector as computed by means of the tree scheme, $a_k^{\text{TREE}},k=1,2,3$, with the “exact” value $a_k^{\text{NBODY}}$ computed by direct summation. A mean error $⟨\delta a_k⟩=\frac{1}{N}\sum _i(a_k^{(i)\mathrm{TREE}}-a_k^{(i)\mathrm{NBODY}})$ is computed by averaging over all the N particles. Then, the relative error is computed as follows:

$$\begin{array}{c}\hfill \mathrm{Err}(a_k)={\displaystyle \frac{\sum _{i=1}^N|a_k^{(i)\mathrm{TREE}}-a_k^{(i)\mathrm{NBODY}}-⟨\delta a_k⟩|}{\sum _i|a_k^{(i)\mathrm{NBODY}}|}}·\end{array}$$(27)

Figure 13 shows these relative errors, obtained with GaSPH, as a function of the CPU time. The figure does not show the relative errors for each single component, but shows the mean values, computed by the simple average $\frac{1}{3}\sum _{k=1}^3\mathrm{Err}(a_k)$. Results for several setup configurations are illustrated. The figure shows the results in three different panels, according to the value of N (N  =  10⁴, N  =  10⁵, and N  =  10⁶). For each value of N, we used different combinations of the parameter θ (0.4, 0.6, 0.8) with different opening criteria (Eq. (20) or Eq. (21)) and different multipole approximations (only the monopole term, or also including the quadrupole term). As the data show, the approximation of the field with the quadrupole moment always represents the best choice in terms of performance because at the same error, it requests a shorter CPU time than the monopole approximation. On the other hand, the choice of the new opening criterion gives a smaller improvement of the error with respect to the benefits that are obtained by switching from monopole to quadrupole term.

Fig. 13. Tree code relative errors Err(ak) (averaged over all the three Cartesian coordinates) for the gravitational field computation as a function of the CPU time. The data are illustrated for different particle numbers (N = 104, N = 105, and N = 106) in panels a, b, and c. In each panel, the full line connects the points related to a computation of the gravitational field made with the quadrupole approximation, while the dashed line refers to computations made by using just the monopole term. The shape of the “void” markers distinguishes different choices of θ with the “standard” criterion (Eq. (20)): θ = 0.8 (void squares), θ = 0.6 (void triangles), and θ = 0.4 (void circles). Results for different opening angles with the “opening law” (Eq. (21)) criterion are also marked: θ = 0.8 (solid squares), θ = 0.6 (solid triangles), and θ = 0.4 (solid circles).

For lower particle numbers, a better computational performance without loss of accuracy can be obtained by including the quadrupole approximation and the (more expensive) opening criterion given in Eq. (21), together with a suitable change of the θ angle. We focus, for example, on the simulation setups characterized by θ = 0.6 with any possible opening criterion, monopole approximation, and N = 10⁴ or N = 10⁵. The change θ  =  0.6 → θ  =  0.8, together with the use of criterion (21), represents a good choice and provides better performing simulations without degrading the precision of the algorithm. We obtain the same advantage when we pass, similarly, from θ = 0.4 to θ = 0.6. On the other hand, for N  =  10⁶ (panel c in Fig. 13), the results related to the approximation with monopole and with quadrupole have smaller differences than the other cases with different N. The choice of a larger opening angle (passing from θ = 0.6 to θ = 0.8 or passing from θ = 0.4 to θ = 0.6) together with the use of the new opening criterion can give a better performance even though the accuracy slightly decreases.

If we wish to preserve the high efficiency of the tree code by keeping the CPU time to scale as N log N, an angle θ ≥ 0.3 must be chosen (Hernquist 1987). The choice of θ = 0.6 in quadrupole approximation or the choice of θ = 0.4 in monopole approximation together with criterion (21) represents a satisfying option because it provides relative errors of at most about 10⁻³.

We now added N_ob = 20 stars to the SPH distribution. As explained in the previous section, a 14th-order explicit method was applied to calculate the evolution of these objects, and giving their low number in comparison to N, very little additional CPU time is expected. Table 1 reports the percentage or workload related to the main relevant subroutines for the new gas+stars system: tree-building + particle-sorting routine, density computation routine, acceleration routine, and the star evolution routine. In addition, we report the rest of the time needed for the basic operations (such as v, u, and r updating, energy computation, and time-step computation). Different work balances are shown for several values of N. The percentage of workload related to the tree-building routine is stable to the order of 6% at different N. In contrast, the work needed by the density routine becomes less and less relevant as N increases, while the gravity+SPH computation becomes increasingly essential. The computational effort to treat the evolution of stars, together with their interaction with the gas, is due both to the pure N-body RK coupling, which is expected to scale as $N_{\text{ob}}^2$, and to the time for coupling each star with each SPH particle, which is expected to be linear in N. Nevertheless, Table 1 shows that 20 stars contribute very little to the total CPU time. For the specific purposes of our current work, we used fewer than two stars or two stars, and their contribution to the code effort is therefore far lower than the 2%÷3%.

4. Protoplanetary disks

This section is dedicated to the description of some tests performed on two main problems involving protoplanetary disks. The first test is to compare the numerical integration of a disk around one star with the analytical prediction (Sect. 4.1). The second test involves a self-gravitating disk that interacts with a binary star (Sect. 4.2).

4.1. Protoplanetary disks around one star

4.1.1. Disk model

Here we illustrate the general setup we used to model a protoplanetary disk in equilibrium around a star of mass M_s = 1 M_⊙. According to the classical flared-disk model (see, e.g., Garcia 2011; Armitage 2011), we let the disk revolve around the central object with a roughly Keplerian frequency ${\mathrm{\Omega }}_k\approx \sqrt{M_{\mathrm{s}}G/R^3}$ (where R is the cylindrical coordinate $R=\sqrt{x^2+{\mathit{y}}^2}$ in the reference frame centered on the central object). The disk evolution is essentially driven by secular viscous dissipation. According to the well-known α-disk model (Shakura & Sunyaev 1973), the turbulence in the internal disk is schematized by means of a pseudo-viscosity of the following form:

$$\begin{array}{c}\hfill \nu ={\alpha }_{\mathit{SS}}{c}_{\mathrm{s}}H.\end{array}$$(28)

This kinematic viscosity perturbs the fluid equations by leading to a net transport of matter inward and an outward flux of angular momentum. α_SS represents a characteristic efficiency coefficient for the momentum transport, while H = c_s/Ω_k represents a characteristic vertical pressure of the disk scale height. The viscous evolution is usually much slower than the dynamical evolution (the characteristic secular timescale is ∝r²/ν, which typically is 2 or 3 orders of magnitude larger than ${\mathrm{\Omega }}_k^{-1}$). This modeling of the turbulence is basically dimensional and is made by mainly taking dynamical turbulence processes into account. Thus, α_SS extends over a wide range of variability (typically 10⁻⁴ and 10⁻²). When the disk self-gravity is strong enough, another important effect arises that is due to the gravitational perturbations. Several works (see, e.g., Mayer et al. 2002; Boss 1998, 2003) have numerically estimated the gravitational timescales in a protoplanetary disk, which is on the same order as its dynamical time. They showed that under certain conditions, matter can undergo instabilities and eventually condense to form clumps in 10³ ÷ 10⁴ yr, which may give rise to gaseous planets. It can been shown that the disks maintain their equilibrium state against collapse according to the Toomre criterion,

$$\begin{array}{c}\hfill Q={\displaystyle \frac{c_{\mathrm{s}}{\mathrm{\Omega }}_{\mathrm{e}}}{\pi G\mathrm{\Sigma }}}>1.5,\end{array}$$(29)

where Ω_e represents the epicyclic frequency, which is approximatively equivalent to Ω_k for Keplerian disks (see Binney & Tremaine 1987; Toomre 1964, for a detailed study). The Toomre factor is a general coefficient that quantifies the predominance of the gravitational processes over the typical thermal and dynamical actions.

We initially let our disk revolve with an azimuthal velocity ${\mathit{v}}_k\equiv {\mathit{v}}_{\varphi }(R)=\sqrt{\frac{G(M_{\mathrm{s}}+M(R))}{R}}$ that depends both on the mass of the central star M_s and on the internal mass of the disk itself $M(R)=\underset{0}{\overset{R}{\int }}\mathrm{\Sigma }(R)2\pi R\mathrm{d}R$. The cumulative mass M(R) can be neglected only for low disk masses M_D ≪ M_s. The shape of the disk in the direction perpendicular to the revolving midplane depends on the vertical pressure scale height H, such that pressure and density scale with a Gaussian profile exp(−z²/2H²). Here a local vertically isothermal approximation was used because we assumed that any radiative input energy from the star is efficiently dissipated away: the cooling times are far shorter than the dynamical timescales. The disk is thus vertically isothermal, and the temperature depends only on the radial distance from the central star.

We set the thermal disk profile according to the well-known flared-disk model, for which the ratio H/R increases with R (see Garcia 2011; Dullemond et al. 2007, for a full clarification). The disk temperature thus follows the profile

$$\begin{array}{c}\hfill T=T_0{\left({\displaystyle \frac{R}{R_0}}\right)}^{-q},\end{array}$$(30)

which is commonly used by setting q = 1/2, while R₀ represents a scale length. We used a slightly different slope q = 3/7, adopted by D’Alessio et al. (1999) by assuming that the thermal processes in the inner layers of the disk do not affect its dynamical stability.

With this temperature profile (independent of t and z), the gas pressure follows a barotropic equation of state $P=c_{\text{s}}^2 \rho$. This choice represents a rough approximation of the cooling processes and allows us to model self-gravitating disks in equilibrium only when Q >  2, which excludes disk models with a state of marginal stability (Q ≈ 1). In a realistic model of a disk without the isothermal approximation, when the Toomre parameter approaches unity, the loss of thermal energy due to radiative cooling processes leads to matter aggregation, which in turn causes shock waves that heat the gas again. If the disk is capable of retaining a sufficient amount of the additional thermal energy that is generated, the collapse only causes some spiral instabilities that do not increase exponentially. The collapse process is thus arrested and the disk reaches a meta-stable state. Every time that a gravitational instability occurs, it is further dissipated by the heat back-production in this meta-state. For a good treatment, see for example Kratter & Lodato (2016). Conversely, the isothermal equation adopted by our model forces the system to cool down at an infinitely high efficiency rate, and to expel all the additional thermal energy that is generated by the compression of matter. Thus, in regions where 1 ≤ Q <  2, the density increases and the collapse process is not halted by production of heat. Our model of a disk in equilibrium is thus limited to masses M_D for which the self-gravity guarantees the condition Q ≥ 2.

We used $\overline{\mu }=2.33$ as the mean molecular weight for the gas. This parameter is commonly adopted to model protoplanetary disks (see, e.g., Kratter et al. 2008; Liu et al. 2017) because it is an average mean weight for a gas composed of H₂ and He, based on the observed cosmic abundance of the elements.

The effects of the Shakura-Sunyaev viscosity associated with Keplerian disks can be emulated by means of the SPH artificial viscosity. Meglicki et al. (1993) found that the SPH viscosity coefficient α provides a viscous acceleration with an effective kinematic viscosity that contains a similar form with a shear component plus a bulk viscosity. When a cubic spline function is used for the kernel function, the authors showed that ν assumes the following form:

$$\begin{array}{c}\hfill \nu \propto \alpha c_{\mathrm{s}}h.\end{array}$$(31)

In several works (e.g., Artymowicz & Lubow 1994; Nelson et al. 1998) the viscosity term in Eq. (11) is used under peculiar conditions: it acts not only for approaching particles, but also for points that move out (with r_ij ⋅ v_ij >  0). It turns out that ν = 0.1 αc_sh (for a comprehensive explanation, see Meru & Bate 2012). We thus have the following law that connects the Shakura–Sunyaev viscosity coefficient to the α parameter used in SPH:

$$\begin{array}{c}\hfill {\alpha }_{\mathit{SS}}={\displaystyle \frac{1}{10}}\alpha {\displaystyle \frac{h}{H}}·\end{array}$$(32)

This modification of the SPH formalism provides a more realistic prediction of the effect given by a kinematic viscosity because it acts under compression and under gas expansion. This prescription is reliable except for strong velocity gradients, that is, shock waves due to strong compressions. In this case, the classical Morris & Monaghan (1997) amplification law (Eq. (14) in this paper) would generate high dissipative forces even in expansion regions. Protoplanetary disks are usually modeled as quiet systems and are not expected to undergo such huge compressions to let strong shock waves arise.

Using Eq. (11) for the viscosity, and consequently, activating dissipation only for particles approaching each other, the law in Eq. (32) can be modified and improved by also considering the effects of the β coefficient on the kinematic viscosity (Meru & Bate 2012; Picogna & Marzari 2013). It follows that for Keplerian disks,

$$\begin{array}{c}\hfill {\alpha }_{\mathit{SS}}={\displaystyle \frac{31}{525}}\alpha {\displaystyle \frac{h}{H}}+{\displaystyle \frac{9}{70\pi }}\beta {\displaystyle \frac{h^2}{H^2}}·\end{array}$$(33)

Equations (32) and (33) formally do not contain any effect of the Balsara switch to compensate for the false sharing attenuation (Eq. (12)). As was done in Picogna & Marzari (2013), we used the artificial viscosity term of the disk in Eq. (11), and multiplied the factor μ_ij by the term of Balsara (1995) $\overline{f}$.

4.1.2. Evolution of the viscous disk

We tested the response of our model with respect to dissipative processes with long timescales as are characteristic of viscous turbulent disks. We started from the well-known disk model developed by Lynden-Bell & Pringle (1974) (see also Pringle 1981; Hartmann et al. 1998). It consists of a thin (H/R ≪ 1) non-self-gravitating disk that is subjected to a power-law dissipative turbulent viscosity. The surface density evolution is described by the following equation (see Pringle 1981; Hartmann et al. 1998):

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathrm{\Sigma }(R,t)}{\partial t}}={\displaystyle \frac{3}{R}}{\displaystyle \frac{\partial }{\partial R}}\left[R^{1/2}{\displaystyle \frac{\partial }{\partial R}}\left(R^{1/2}\nu \mathrm{\Sigma }(R,t)\right)\right].\end{array}$$(34)

It has been shown that when the disk is perturbed by a radial power-scaling kinematic viscosity ν ∝ R^ϕ, the differential equation admits the following similarity solution (see Lynden-Bell & Pringle 1974; Hartmann et al. 1998):

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R,t)={\mathrm{\Sigma }}_0{\left({\displaystyle \frac{R}{R_1}}\right)}^{-\varphi }{({\displaystyle \frac{t}{{\tau }_{\nu }}}+1)}^{-\gamma }exp[-{({\displaystyle \frac{t}{{\tau }_{\nu }}}+1)}^{-1}{\left({\displaystyle \frac{R}{R_1}}\right)}^{2-\varphi }],\end{array}$$(35)

where $\gamma =\frac{5/2-\varphi }{2-\varphi }$. R₁ is a characteristic radial scale that contains about the 68% of the total disk mass, while Σ₀ is the normalization scale density. In Eq. (35), ${\tau }_{\nu }=\frac{R_1^2}{3{(2-\varphi )}^2{\nu }_1}$ represents the characteristic viscous timescale of the disk, and it is proportional to the inverse of the viscosity evaluated in correspondence to the scale radius (ν₁ = ν(R₁)).

We sampled a disk, made of 20 000 particles, which followed the thermal law of Eq. (30) described in the previous section, where we assumed R₀ = 10 AU for scaling and T₀ = 25 K. The SPH particle distribution was achieved by sampling the initial (t = 0) radial density profile as from Eq. (35),

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R,0)={\mathrm{\Sigma }}_0{\left({\displaystyle \frac{R}{R_1}}\right)}^{-\varphi }exp[-{\left({\displaystyle \frac{R}{R_1}}\right)}^{2-\varphi }],\end{array}$$(36)

where R₁ was set here to 50 AU. We used a constant value for the viscosity coefficient, α_SS  =  10⁻². The kinematic viscosity is proportional to the sound speed and the vertical scale height (Eq. (28)), and when we assume that the exponent q = 3/7 in the power-law Eq. (30), ν follows a radial power law with a positive exponent ϕ = 15/14 ≈ 1.07. In order to reproduce the effects of such a radial viscosity law, we imposed a specific prescription for the artificial viscosity in the code. Taking the law of Eq. (33) into account and considering that we used β = 2α, we obtain the following expression for the coefficient α:

$$\begin{array}{c}\hfill \alpha ={\alpha }_{\mathit{SS}}{\displaystyle \frac{H}{h}}{\displaystyle \frac{1}{31/525+(9/35\pi )(h/H)}},\end{array}$$(37)

with constant α_SS  =  10⁻². To apply this form, we inserted the expression above in the artificial viscosity term of Eq. (11), which we used in all our disk simulations, without considering the Morris & Monaghan (1997) variation law.

The (infinite) disk has been truncated at R_out  =  8R₁  =  400 AU. The distribution was also truncated at an inner cutoff R_in  =  R₁/5  =  10 AU. The gravitational softening radius of the central star together with its sink radius were set equal to R_in. The inner border condition we set is that the gas particles that cross the R_in radius are absorbed by the star and are therefore excluded from time integration. The mass of the sunk SPH particles was considered to increase the mass of the central star.

The ratio R_out/R_in  ≈  40 is high enough to match the infinite extension of the analytical density profile as accurately as possible. This reduces the external radial boundary discontinuity. Hartmann et al. (1998) indeed pointed out that the inward flux of matter (which constitutes the most important process in guiding the evolution of an α disk) considerably depends on the outward angular momentum transportation and thus on the disk expansion through the external shells.

The boundary conditions at R_in >  0 may affect the disk evolution throughout the whole spatial extension. As pointed out in Lynden-Bell & Pringle (1974) and remarked by Hartmann et al. (1998), a viscous disk may formally extend to R → 0, but at a critical radius that is on the same order as the radius of the star, both the torque and the viscosity ν go to zero. Computational efficiency purposes prevent us from modeling a disk by introducing a very small cutoff radius. We imposed zero viscosity ν by varying the coefficient α within R = 3R_in down to zero, assumed at R = R_in.

Parameter α in Eq. (37) varies with time because the disk evolution processes lead to a time variation of the ratio between h and the height scale H. For R = R₁/2 the disk initially has an average h/H ≈ 2, while this ratio reaches a minimum value of about 1.5, in correspondence to R = 2R₁, which provides a value of α ≈ 0.02 and α ≈ 0.04, respectively. As the gas is captured by the central star, we expect the density to decrease and the h-to-H ratio to increase. We show below that the surface density varies quickly during an initial phase and later evolves at a relatively lower rate; this causes h/H to follow the same cadence. At later stages (after about 1.6 Myr), the disk is indeed integrated with lower values of α that range from approximately 0.01 to 0.03.

Our disk is virtually non-selfgravitating because Q  ≫  2 throughout its surface (the minimum value is about 20), even though we formally took the disk mass into account in setting the azimuthal velocity v_ϕ(R). For R = R₁, it has a vertical aspect ratio H/R ≈ 0.05. With this setup, the disk should have a viscous evolution timescale τ_ν ≈ 840 000 yr according to the analytical model.

During the earlier integration phases, the disk experiences a fast relaxation caused by the inner cutoff. In this stage the internal density discontinuity is smoothed out and fades out near the inner border. We therefore allowed the system to relax for a time of about 50 000 yr (about 900 Keplerian orbits at R = R₁), which is far shorter than τ_ν and sufficient to obtain a steady state. From this time, the only expected changes in the disk density profile are those that are due to the secular viscous evolution. Thus, after 50 000 yr, the disk is quickly arranged in a configuration that is slightly different from the initial one: it is characterized by a density distribution that follows the density law of Eq. (36), but now with a larger radius R₁ ≈ 60 AU. With this new radius, the disk has an evolution timescale τ_ν ≈ 10⁶ yr. The properties of self-similarity according to Eq. (35) indeed guarantee that the surface density profile maintains the same analytical form for every time. Thus, the surface density of a disk can at every time be considered an initial solution of a new disk that is described by Eq. (36), with a different parameter R₁ and thus a different viscosity timescale τ_ν. The initial profile and the density profile after 50 000 yr (which we set conventionally as the instant t = 0) are shown in the top panel of Fig. 15. The initial t = 0 state fits (full line) the disk profile of Eq. (36) with R₁ ≈ 60 AU. We considered the disk evolution to start from this configuration and plot the density profile for various times (t = τ_ν/3 ≃ 330 000 yr and t ≃ 1 000 000 yr). We also compare this with the analytical predictions obtained by Eq. (35) (dashed lines). Results and model do not match because the numerical disk appears to evolve faster than the disk predicted by the analytical theory. For t ≃ 330 000 yr the discrepancy between the analytical density was about 45% higher than the numerical result at R = R₁, while for t ≃ 1 000 000 yr the discrepancy increased up to 95%. The disk thus looses mass more quickly than would be expected based on analytical theory, and thus its density decreases too rapidly. This excessive loss can be ascribed to the inner border, which is absent in the analytical theory, where the disk matter flows onto the star in a single point. To check how incorrectly the SPH disk depleted its mass, we considered the radial-cumulative disk mass, predicted at a given time by the analytical theory, which is strictly connected to the viscosity timescale,

$$\begin{array}{c}\hfill M(\le R,t)=M_{\mathrm{D}}(0){({\displaystyle \frac{t}{{\tau }_{\nu }}}+1)}^{\frac{-1}{2(2-\varphi )}}[1-exp\left({\left({\displaystyle \frac{R}{R_1}}\right)}^{2-\varphi }{({\displaystyle \frac{t}{{\tau }_{\nu }}}+1)}^{-1}\right)],\end{array}$$(38)

where M_D(0) is the starting disk mass at t = 0. In Fig. 15 we illustrate the fractional disk mass-loss rate $-{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{M(t)}{M(0)}}\right)$ as a function of time. The result obtained with our model is compared with the theoretical rate, where the mass as a function of time was considered as the integral of the analytical surface density in Eq. (35), from R = 10 AU and R = ∞, that is, $M(t)=\underset{10\mathrm{AU}}{\overset{\infty }{\int }}\mathrm{\Sigma }(R,t)2\pi R\mathrm{d}R=M_{\mathrm{D}}(t)-M(\le 10\mathrm{AU})$, where M_D(t) is the total disk mass at a given time. The SPH-disk mass-loss rate tends to reach the analytical curve only for t → τ_ν), while for earlier stages, we have a huge mass-loss rate that decreases as time approaches τ_ν.

To make a more substantial comparison, we moved out from the early stages and focused on the later phases at t′=850 000, where the mass-loss rate illustrated in Fig. 15 approached the analytic value to within 10%. For this time, we considered our model density distribution as the starting state of a new disk and studied its subsequent evolution. We summarize this evolution in Fig. 16. The figure shows that at t′, the disk density surface corresponds to a Lynden-Bell solution of the same form as the starting one, but with a different characteristic radius $R_1^{\prime }\approx 153\text{AU}\simeq 2.5R_1$. This radius corresponds to a viscous timescale ${\tau }_{\nu }^{\prime }\approx 2.4\times {10}^6$ yr. In the same figure we show the evolution of the new surface density after a time of ${\tau }_{\nu }^{\prime }/3\approx \mathrm{790\; 000}$ yr, together with the theoretically expected value (dashed line). This was made by analogy with Fig. 14, where the result at t = τ_ν/3 ≈ 330 000 yr is shown (plotted in the middle). Comparing the two curves, we note that during the later stages the density of the SPH disk shows less deviation from the analytical prediction (within an error of 20% at $R=R_1^{\prime }$), although the discrepancies are not negligible. The discrepancies are also smaller in internal regions at $R=20\text{AU}\ll R_1^{\prime }$, close to the inner border. To relate the discrepancy between the numerical density and its analytical prediction, we further verified this by adopting a smaller inner border radius, that is, by reducing the star sink from 10 AU to 5 AU, and we indeed observed a substantial reduction of the inner mass-loss rate of the disk in the initial integration phases. We also note that the gap of the surface density within the star sink radius is reduced in amplitude when a smaller inner boundary is chosen.

Fig. 14. Surface density profiles Σ(R) of the disk at different times (t = 0 corresponds to the disk state after it has been relaxed up to a time of 50 000 yr). The dots represent the numerical results, while the lines refer to the analytical model (Eq. (35)). The actual initial disk setting (black triangles) is reported for the sake of comparison with the relaxed state. The abscissa is in AU and the density is expressed in 10−2 M⊙ AU−2. The evolved density profiles at t >  0 (t = 330 000 yr and t = 1 000 000 yr) are shifted down by −3 and −6 in logarithm, respectively, for clarity. The fitting curve of the density profile at t = 0 is also plotted (full lines).

Fig. 15. Relative (percentage) mass-loss rate of the disk $-{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{M(t)}{M_0}}\right)$, expressed in units of Myr−1, and the time is in Myr. The results of our numerical simulation (dots) are compared with the theoretical behavior (full line).

Fig. 16. Logarithm of the numerical radial profile of the disk density at t′=850 000 yr (dots). The density profile at the beginning of the simulation is also plotted (black triangles). The density at t = t′ matches the self-similar solution (Eq. (36)) with a parameter $R_1^{\prime }\approx 153$ AU and a viscosity timescale ${\tau }_{\nu }^{\prime }=2.4\times {10}^6$ yr. The evolution of this disk after a time ${\tau }_{\nu }^{\prime }/3$ is also plotted (bottom, artificially shifted by −3) together with the analytical prediction (dashed line). The fit density profile at t = t′ is plotted twice (in the top and bottom curves) as a full line for clarity.

4.2. Self-gravitating disk in a binary system

We tested our code on a more complex dynamical system by treating the evolution of a self-gravitating disk that interacts with a binary star. The system is characterized by a circumprimary disk around a 1 M_⊙ star, truncated by the gravitational field of a 0.4 M_⊙ external companion star. This topic has been treated in a relevant work by Marzari et al. (2009), who integrated the time evolution of this configuration and studied the effects of the stars on the parameters of the orbital disk: eccentricity and argument of periastron. The authors used the well-known Eulerian code FARGO, implemented with a full scheme for the self-gravity (see Masset 2000; Baruteau & Masset 2008), and performed a 2D simulation. We investigated the evolution of a system like this with GaSPH with a 3D model for the particular case of a binary with eccentricity 0.4.

4.2.1. System setup

The orbit of the binary system is characterized by an eccentricity e_b = 0.4 and a semimajor axis of 30 AU. As in Marzari et al. (2009), we kept the orbit of the two stars fixed during the integration because we focus on the gravitational effects of the two stars on the disk. This means that the dynamics of the binary star is not affected by the gas feedback. The initial configuration of the disk that was adopted in the original model is characterized by a radial surface density Σ ∝ R^−1/2 that extended from 0.5 AU to 11 AU from the central primary star. Beyond 11 AU, the density quickly fades out; the total disk mass is M_D = 0.04 M_⊙. Instead of using a typical flared disk, we adopted a flat disk by setting the linear vertical scale height to H = 0.05 R.

The choice adopted by Marzari et al. (2009) for the disk shape and the viscosity law leads to a slightly different schematization than in the model we used in Sect. 4.1.2, for which we needed to set the parameters inside our SPH 3D code. The authors used a constant kinematic viscosity ν, which leads α_SS to vary. Furthermore, the constant value of the aspect ratio H/R causes the speed of sound to scale as c_s ∝ R^−1/2. When ν = cost. and given the alpha-disk law of Eq. (28), we have that ${\alpha }_{\mathit{SS}}\propto H^{-1}{c}_{\mathrm{s}}^{-1}=H^{-2}{\mathrm{\Omega }}^{-1}\propto R^{-1/2}$. The coefficient α_SS was indeed set by the above authors by calibrating it to correspond to the value α_SS = 2.5 × 10⁻³ in the central regions about 5 AU within the disk. Thus, we can deduce it as follows:

$$\begin{array}{c}\hfill {\alpha }_{\mathit{SS}}=2.5\times {10}^{-3}{\left(\frac{R}{R_{\mathrm{ref}}}\right)}^{-1/2},\end{array}$$(39)

with R_ref = 5 AU. We applied the artificial viscosity term by setting the SPH α parameter according to the same expression of Eq. (37) as we used in the previous section, but now using the nonconstant coefficient α_SS with the radial profile of Eq. (39) illustrated above.

The disk is coplanar with the star orbit, and we built it by confining a set of SPH particles between R = 0.5 AU and R = 11 AU. We integrated the system for about 3000 yr. All the gas particles that flew across the inner border were excluded from the integration. Three runs for the same model were made with different particle numbers, N = 20 000, N = 50 000, and N = 100 000.

As in the case of protoplanetary disk around one star we discussed in the previous section, the ratio h/H, and consequently α, are not constant in time. For N = 20 000, the disk acquires a steady state configuration after an initial quick relaxation phase, where h/H has a rather slow evolution along the timescale of a binary rotation period. After about 500 yr, the disk acquires an average ratio h/H ≈ 6.5 (α ≈ 0.001) in the inner regions (R = 1 AU), which reaches a minimum of about 1 (α ≈ 0.011) in the middle regions (R = 5.5 AU). Similarly, the disk with N = 50 000 has an average ratio h/H ≈ 3.7 in the inner regions, with a minimum of h/H ≈ 0.7, corresponding to α ≈ 0.004 and α ≈ 0.02. The disk with N = 100 000 has α ≈ 0.008 in the internal regions and a maximum α ≈ 0.03 in the intermediate radial regions. Regardless of the resolution, after an integration time of 3000 yr, the h-to-H ratio slightly changes in the inner regions while its minimum substantially increases, giving rise to a maximum value of α of about 0.008, 0.015, and 0.02 for the disk with the lower, medium, and higher particle numbers, respectively.

4.2.2. Disk deformation and gravitational feedback to the stars

The gravitational field of the stars affects the disk configuration by altering its average orbital parameters such as the eccentricity and the argument of periastron. In order to describe the disk evolution, we calculated the mean eccentricity and the mean argument of periastron by averaging over the disk surface, with the same prescription as was adopted by Pierens & Nelson (2007),

$$\begin{array}{cc}& e_{\mathrm{disk}}={\displaystyle \frac{1}{M_{\mathrm{D}}}}\sum _i{e}_i(R,\varphi )m_i,R_A\le R\le R_B\hfill \\ \hfill & {\omega }_{\mathrm{disk}}={\displaystyle \frac{1}{M_{\mathrm{D}}}}\sum _i{\omega }_i(R,\varphi )m_{\mathrm{i}},R_A\le R\le R_B.\hfill \end{array}$$(40)

Here M_D is the disk mass included within R_A and R_B (the latter being a quantity on the order of the effective radius R_D). The disk radius is defined by the following expression:

$$\begin{array}{c}\hfill R_{\mathrm{D}}\equiv {⟨R_{\mathrm{D}}⟩}_{\mathrm{L}}\propto {\left({\displaystyle \frac{L}{M_{\mathrm{D}}}}\right)}^2\end{array}$$(41)

and is computed as the radial distance containing the total angular momentum L of the disk, with M_D its total mass. R_D is not very different from the half-mass radius.

The local orbital parameters are evaluated using the eccentricity vector,

$$\begin{array}{c}\hfill \mathit{e}={\displaystyle \frac{\mathit{r}x\mathit{l}}{G{M}_{\mathrm{c}}}}-\widehat{\mathit{r}},\end{array}$$(42)

where l = rxv represents the angular momentum per unit mass, and M_c is the mass of the central star. The vector e characterizes the orbit described by the position r of a point about the center of mass of the binary, assuming that it corresponds to an elliptical trajectory. The absolute value e = |e| corresponds to the orbital eccentricity. Moreover, e is always parallel to the semimajor axis, thus, its normalized components e_x/e = cos(ω) and e_y/e = sin(ω) provide the local argument of periastron ω, and finally, the local semimajor axis. In Eq. (40), only the particles on orbits that are tied to the central primary star are considered part of the disk and thus are included in the summation.

The choice of a 3D model introduces some new degrees of freedom with respect to a 2D scheme because the self-gravity of the disk and the gravity of the stars act even along the vertical direction. In particular, the angular momentum, whose flux across the disk plays a crucial role on the disk evolution itself, is also spread out in the vertical direction.

Figures 17 and 18 show the evolution of the disk eccentricity and of the angle ω_disk, respectively. After a few orbital periods, the secondary star initially truncates the disk and a chaotic phase arises during which the mean eccentricity increases abruptly. In a second phase, the eccentricity stabilizes around an average value that is similar to the value obtained by Marzari et al. (2009), which is e_disk ≈ 0.075. Moreover, in the same way as the other investigation, Fig. 17 shows that e_disk slightly oscillates. The oscillation is modulated with the binary period (P_bin ≈ 134 yr), and can be ascribed to the strong variation of the gravitational field of the companion star at periastron.

Fig. 17. Disk eccentricity, edisk, evolution under the perturbation of a binary system of e = 0.4. Values referring to different simulations are plotted: N = 20 000 (dotted line), N = 50 000 (dashed line), and N = 100 000 (full line).

Fig. 18. Time evolution of argument of periastron, ωdisk, in an eccentric (e = 0.4) binary system. Like in Fig. 17, results for different simulations are illustrated: N = 20 000 (dotted line), N = 50 000 (dashed line), and N = 100 000 (full line).

Similarly, Fig. 18 shows the mean inclination of the oscialltions in the disk semimajor axis around the initial value (which conventionally was taken as π). Even for the argument of periastron of the disk, a convergence can be observed near the simulations with higher resolution.

5. Summary and conclusions

The primary intent of this paper is the presentation, testing, and preliminary application of our new SPH code, GaSPH, which is designed to be a multi-purpose code that is applicable to a variety of astrophysical, multi-phase, and self-gravitating environments. We briefly summarize the main points below.

• We presented and discussed the characteristics of our code in some detail. At the moment, the code does not treat radiative transfer, but properly takes internal gas gravity and the gas-star mutual gravity into account.

• The code fully passes the classic tests in slowly varying situations, which assess its stability, and also in violently varying cases that reproduce the Sedov–Taylor blast wave well.

• The code performs well numerically (speed) and has a good stability and quality, as we discussed in Sect. 3.2.

• The capability of the code to treat the evolution of a protoplanetary disk when it interacts with a single star and with a binary star was tested.

In a near future, we aim to reach a much better resolution, which is achievable with an MPI parallel version of our code. This will allow us to study disks on the smaller, planetary, scale. Further scientific applications of our code will include the evolution of protoplanetary disks in a star cluster environment in order to study the star-to-disk feedback.

Acknowledgments

We express our gratitude to Fabrizio Capaccioni of INAF-IAPS-Istituto di Astrofisica e Planetologia Spaziali (Rome, Italy) for the valuable support and useful discussions throughout the preparation of this work. We also wish to thank the anonymous referee for the suggestions and comments that greatly helped to improve the paper.

References

Appendix A: Details of the numerical method

A.1. Velocity-Verlet method

At each time iteration, a and $\dot{u}$ are evaluated twice, in correspondence to the current step n and the next n + 1. Starting from a generic nth time iteration, we use a^[n] and ${\dot{u}}^{[n]}$ to predict the velocity and the energy:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{v}}^{\ast [n+1]}={\mathit{v}}^{[n]}+{\mathit{a}}^{[n]}\mathrm{\Delta }t\hfill \\ u^{\ast [n+1]}=u^{[n]}+{\dot{u}}^{[n]}\mathrm{\Delta }t\hfill \end{array},\end{array}$$(A.1)

while the position is directly updated to the next step n + 1, without any prediction:

$$\begin{array}{c}\hfill {\mathit{r}}^{[n+1]}={\mathit{r}}^{[n]}+{\mathit{v}}^{[n]}\mathrm{\Delta }t+{\displaystyle \frac{1}{2}}{\mathit{a}}^{[n]}\mathrm{\Delta }{t}^2.\end{array}$$(A.2)

With these new quantities, a new calculation is performed for a and $\dot{u}$; we thus have

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{a}}^{\ast [n+1]}=\mathit{a}({\mathit{r}}^{[n+1]},{\mathit{v}}^{\ast [n+1]},u^{\ast [n+1]})\hfill \\ {\dot{u}}^{\ast [n+1]}=\dot{u}({\mathit{r}}^{[n+1]},{\mathit{v}}^{\ast [n+1]},u^{\ast [n+1]})\hfill \end{array},\end{array}$$(A.3)

which we can use to correct velocity and energy,

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{v}}^{[n+1]}={\mathit{v}}^{[n]}+({\mathit{a}}^{[n]}+{\mathit{a}}^{\ast [n+1]}){\displaystyle \frac{\mathrm{\Delta }t}{2}}\hfill \\ u^{[n+1]}=u^{[n]}+({\dot{u}}^{[n]}+{\dot{u}}^{\ast [n+1]}){\displaystyle \frac{\mathrm{\Delta }t}{2}}\hfill \end{array}.\end{array}$$(A.4)

It can be straightforwardly shown that when acceleration only depends on the positions, that is, in a Newtonian problem without hydrodynamics, the Verlet method described above is equivalent to a standard second-order kick-drift-kick (KDK) leap-frog method. When the acceleration depends neither on the velocity field nor on the internal energy, the quantity a^*[n + 1] corresponds to the actual acceleration a^[n + 1] that is related to the next step. The numerical method can thus be rewritten in the following way:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{r}}^{[n+1]}={\mathit{r}}^{[n]}+({\mathit{v}}^{[n]}+{\displaystyle \frac{1}{2}}{\mathit{a}}^{[n]}\mathrm{\Delta }t)\mathrm{\Delta }t={\mathit{r}}^{[n]}+{\mathit{v}}^{[n+1/2]}\mathrm{\Delta }t\hfill \\ {\mathit{v}}^{[n+1]}={\mathit{v}}^{[n]}+{\displaystyle \frac{1}{2}}{\mathit{a}}^{[n]}\mathrm{\Delta }t+{\displaystyle \frac{1}{2}}{\mathit{a}}^{[n]}\mathrm{\Delta }t={\mathit{v}}^{[n+1/2]}+{\mathit{a}}^{[n]}{\displaystyle \frac{\mathrm{\Delta }t}{2}}\hfill \end{array},\end{array}$$(A.5)

which represents indeed the standard expression of a KDK leap-frog integrator that requires just one force calculation per time step (see Hockney & Eastwood 1988; Hut et al. 1995; Quinn et al. 1997, for leap-frog methods and further improvements).

Each particles has its own individual time steps, sorted by the code as submultiples of the maximum time step Δt_max. For a generic particle of index i, the time step is calculated according to the criteria expressed by Eqs. (17) and (18) and approximated to the nearest value Δt_i = 2^−P ⋅ Δt_max, where P is a positive integer number. The simple scheme in Fig. A.1 shows that a single time iteration between two consecutive steps n and n + 1 is performed between the time t^[n] and t^[n + 1] = t^[n] + Δt_min, while a generic particle of index i is updated between two characteristic times: $t_i^{[\text{PREV}]}$ and $t_i^{\text{NEXT}}=t_i^{[\text{PREV}]}+\mathrm{\Delta }{t}_i$. The routines dedicated to the neighbor search, to the hydrodynamic forces, and to the gas self-gravity are activated for the particle i only when its time step is synchronized, that is, if the conditions $t^{[n]}=t_i^{[\text{PREV}]}$ or $t^{[n+1]}=t_i^{[\text{NEXT}]}$ are satisfied. Figure A.2 shows a scheme of a single time iteration. First, the code sorts and maps all the particles by building the tree, calculating the quadrupole momentum and all the other key quantities related to the cubes. This operation is thus independent on the particle time steps. Then, the code runs the main cycles of neighbor search and acceleration computation. For each ith particle, it verifies whether the particle is synchronized at t = t^[n], that is, whether the condition $t^{[n]}=t_i^{[\text{PREV}]}$ is satisfied. In that case, the algorithm computes the density ${\rho }_i^{[\text{PREV}]}$, the pressure $P_i^{[\text{PREV}]}$, the velocity gradient $\nabla \cdot {\mathit{v}}_i^{[\text{PREV}]}$ and the velocity rotor $\nabla \times {\mathit{v}}_i^{[\text{PREV}]}$, the switching coefficient of Eq. (12), the quantities ${\omega }_i^{[\text{PREV}]},{\zeta }_i^{[\text{PREV}]}$, the hydrogravitational acceleration ${\mathit{a}}_i^{[\text{PREV}]}$, and the time variation of the internal energy ${\dot{u}}_i^{[\mathrm{PREV}]}$. These quantities are suitably stored in memory, to be used in further phases of the integration and in further iterations. For the remaining nonsynchronized particles, the algorithm indeed uses the quantities that were calculated in previous stages.

Fig. A.1. Sketch of the hierarchical time-step subdivision. A generic individual time step, Δti is a power of two multiples of the minimum time step Δtmin, and a power of two submultiple of the maximum one, Δtmax. The integration is performed by means of elementary iterations from t[n] to t[n + 1] = t[n] + Δtmin, while the particle is updated from t[PREV] to t[NEXT] = t[PREV] + Δti. The hydrodynamics and force routines are activated only for particles that are synchronized at t[n] = t[PREV] or at t[n + 1] = t[NEXT].

Fig. A.2. Flow chart showing the main scheme of a single time iteration.

If some stars are included in the simulation, $a_i^{[\text{PREV}]}$ is incremented by a contribution from the star-gas interaction for every ith particle.

After this, we update the time of v, u, and r. For particles synchronized at t = t^[n + 1], so that $t_i^{[\text{NEXT}]}=t^{[n+1]}$, the velocity and the energy are updated with the same predictor scheme as expressed by Eq. (A.1), thus, we have

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{v}}_i^{\ast [n+1]}={\mathit{v}}_i^{\ast [\mathrm{NEXT}]}={\mathit{v}}_i^{[\mathrm{PREV}]}+{\mathit{a}}_i^{[\mathrm{PREV}]}\mathrm{\Delta }{t}_i\hfill \\ u_i^{\ast [n+1]}=u_i^{\ast [\mathrm{NEXT}]}=u_i^{[\mathrm{PREV}]}+{\dot{u_i}}^{[\mathrm{PREV}]}\mathrm{\Delta }{t}_i\hfill \end{array}.\end{array}$$(A.6)

The position is updated as well, in the same manner as indicated by Eq. (A.2):

$$\begin{array}{c}\hfill {{\mathit{r}}_{\mathit{i}}}^{[n+1]}={\mathit{r}}_i^{[\mathrm{NEXT}]}={\mathit{r}}_i^{[\mathrm{PREV}]}+{\mathit{v}}_i^{[\mathrm{PREV}]}\mathrm{\Delta }{t}_i+{\displaystyle \frac{1}{2}}{\mathit{a}}_i^{[\mathrm{PREV}]}\mathrm{\Delta }{t}_i^2.\end{array}$$(A.7)

On the other hand, for nonsynchronized points, we do not update, but estimate the quantities v^*[n + 1], u^*[n + 1] and r^[n + 1] by means of the following predictor scheme:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathit{v}}_i^{\ast [n+1]}={\mathit{v}}_i^{[\mathrm{PREV}]}+{\mathit{a}}_i^{[\mathrm{PREDICT}]}\delta t_i\hfill \\ u_i^{\ast [n+1]}=u_i^{[\mathrm{PREV}]}+{\dot{u}}_i^{[\mathrm{PREV}]}\delta t_i\hfill \\ {\mathit{r}}_i^{[n+1]}={{\mathit{r}}_{\mathit{i}}}^{[\mathrm{PREV}]}+{\mathit{v}}_i^{[\mathrm{PREDICT}]}\delta t_i+{\displaystyle \frac{1}{2}}{\mathit{a}}_i^{[\mathrm{PREDICT}]}\delta t_i^2,\hfill \end{array}\end{array}$$(A.8)

which contains the following quantities:

$$\begin{array}{cc}& \delta t_i=t^{[n+1]}-t_i^{[\mathrm{PREV}]}\hfill \\ \hfill & {\mathit{a}}_i^{[\mathrm{PREDICT}]}=\frac{1}{2}({\mathit{a}}_i^{[\mathrm{PREV}]}+{\mathit{a}}_i^{[n]})\hfill \\ \hfill & {\mathit{v}}_i^{[\mathrm{PREDICT}]}=\frac{1}{2}({\mathit{v}}_i^{[\mathrm{PREV}]}+{\mathit{v}}_i^{\ast [n+1]})\hfill \end{array}$$

${\mathit{a}}_i^{[n]}$ is exactly equal to the old value ${\mathit{a}}_i^{[\text{PREV}]}$ in simulations with pure gas because the algorithm does not compute the acceleration at the current iteration. When one or more stars interact with the gas, however, the two accelerations are different because they contain different contributions from the gas-star interaction, which is independently calculated at every step in the particle synchronization.