© ESO 2020

1 Introduction

Planets grow in protoplanetary discs by the gradual accumulation of material, starting with micrometre-sized dust grains. Sticking collisions lead to the formation of millimetre-sized pebbles; growth barriers such as radial drift, bouncing, and fragmentation prevent the bodies from growing to larger sizes (Weidenschilling 1977a; Blum & Wurm 2008; Brauer et al. 2008; Güttler et al. 2010; Zsom et al. 2010; Birnstiel et al. 2012). Therefore, other mechanisms are necessary to form planetesimals. Streaming instability is currently thought to be the main channel that leads to the formation of planetesimals, which typically have sizes from a hundred kilometres up to about Ceres size (Youdin & Goodman 2005; Johansen et al. 2007, 2014; Carrera et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Yang et al. 2017). These planetesimals grow to protoplanets and the cores of giant planets through mutual collisions and the accretion of the remnant millimetre-sized pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Ormel & Kobayashi 2012; Ormel 2017; Lambrechts et al. 2019; Liu et al. 2019). Lastly, collisions between protoplanets and the accretion of gas from the protoplanetary disc by protoplanets more massive than ~10 M_⊕ lead to the formation of Earth-like planets and gas giants (Pollack et al. 1996; Raymond et al. 2009; Bitsch et al. 2019; Schulik et al. 2019).

In the intermediate phase where planetesimals grow by collisions and pebble accretion, their dynamics plays a crucial role. The velocity distribution of planetesimals directly affects the accretion rate. However, the total number of planetesimals in protoplanetary discs is high and, in general, the gravitational interaction between all of them needs to be taken into account. Recent progress in the development and optimisation of N-body integrators, which use an efficient splitting of the N-body Hamiltonian (Wisdom & Holman 1991; Duncan et al. 1998; Chambers 1999; Rein & Tamayo 2015), have made it possible to perform dynamical studies of Solar System formation and evolution. However, when included, numerous planetesimal-sized bodies are typically treated either as non-interacting test particles that are only scattered and accreted by planets, or as super particles that represent a subset of the planetesimals.

Because of the low density of planetesimal discs, stirring of eccentricity and inclination comes mainly from the cumulative effect of the scattering of two planetesimals that undergo close encounters in the gravitational field of the Sun. Therefore, a different approach to study planetesimal dynamics is to use statistical methods. Based on three-body calculations for the outcome of planetesimal scattering and averaging over an assumed velocity distribution, semi-analytic models (Stewart & Kaula 1980; Stewart & Wetherill 1988; Ida 1990; Ohtsuki 1999; Stewart & Ida 2000; Ohtsuki et al. 2002) successfully reproduce N-body simulations (Ida & Makino 1992a,b, 1993; Palmer et al. 1993; Aarseth et al. 1993) and provide a way to account for planetesimal dynamics in toy models of planet formation (e.g. Chambers 2006a,b; Ormel & Kobayashi 2012). A third approach to modelling planetesimal dynamics is based on the fact that the dominant contribution to orbital changes arises from close encounters between two planetesimals (Henon & Petit 1986; Petit & Henon 1986; Weidenschilling 1989; Hasegawa & Nakazawa 1990; Ida 1990). Therefore, as in the statistical approach, planetesimal dynamics is simulated as the cumulative effect of close encounters. These method types were already used by Cox & Lewis (1980), Lecar & Aarseth (1986), and Beauge & Aarseth (1990). Cox & Lewis (1980) evolved planetesimals along Keplerian orbits as long as their distance was larger than the Tisserand sphere of influence. Inside this sphere the close encounter was resolved by assuming hyperbolic trajectories for both planetesimals. A different approach was chosen by Lecar & Aarseth (1986). Instead of resolving the close encounter only when inside the sphere of influence, they divided the orbital plane into radial and azimuthal zones to include the gravitational perturbations of the nearest neighbours, given in terms of the grid, when integrating the equations of motion of the planetesimals. Lecar & Aarseth (1986) pointed out that the method of Cox & Lewis (1980) misses out the close encounters that dominate the stirring because the Tisserand sphere of influence is smaller than the relevant encounter distance they derived, which is comparable to the size of the Hill radius.

Here we use a hybrid approach combining the ideas of symplectic N-body methods (e.g. Chambers 1999) and close encounters (Cox & Lewis 1980; Lecar & Aarseth 1986). Instead of a statistical treatment of planetesimal dynamics, we develop a method that efficiently evolves a large number of planetesimals on their orbits taking mutual gravitational interaction into account. This method has the advantage that we can simulate many bodies in reasonable computational time, in contrast to conventional N-body methods, and does not rely on averaging over distribution functions, as done in the statistical methods. We achieve this by separating the dynamics of the planetesimals into unperturbed Keplerian motion and two-body scattering in the solar gravitational field and by using an efficient method for the detection of close encounters. To make the method even more versatile and applicable to studies of planet formation, we treat planets in the conventional N-body sense.

We outline the method in Sect. 2; we start by introducing the concept of symplectic N-body integrators and continue with a technical description of the different components of our method. In Sect. 3, we present the results of different test cases, performance, and scaling tests. Limitations are discussed in Sect. 4 before we summarise and conclude in Sect. 5. Algorithms, our prescription for pebble accretion, gas drag, and migration can be found in Appendices A–C.

2 Method

We construct our close-encounter method based on the theory of symplectic integrators for the gravitational N-body problem. Symplectic integrators have become popular for this because they use the Hamiltonian to construct efficient and accurate numerical schemes. Symplectic integrators preserve the Hamiltonian structure by conserving the phase space volume and by solving a Hamiltonian that is close to the true Hamiltonian of the system, which means energy is conserved (Yoshida 1993; Hernandez & Dehnen 2017).

2.1 Theory of symplectic integrators

The Hamiltonian of the gravitational N-body problem is the sum of kinetic and potential energy, $$\mathcal{H}={\sum }_{i=1}^N\frac{p_i^2}{2m_i}-{\sum }_{i=1}^N{\sum }_{j=i+1}^N\frac{G{m}_i{m}_j}{|x_i-x_j|},$$(1)

where p_i = m_iv_i and x_i are the momentum and the position of body i, and the equations of motion are given by Hamilton’s equations: $$\stackrel{.}{x_i}=\frac{\partial \mathcal{H}}{\partial p_i},\stackrel{.}{p_i}=-\frac{\partial \mathcal{H}}{\partial x_i},i\in (1,\dots ,N).$$(2)

The time derivative of any quantity f that is a function of the x_i and p_i can then be written as $$\dot{f}={\sum }_{i=1}^N\left(\frac{\partial f}{\partial x_i}\frac{\partial \mathcal{H}}{\partial p_i}-\frac{\partial f}{\partial p_i}\frac{\partial \mathcal{H}}{\partial x_i}\right)=\widehat{D}f,$$(3)

where $\widehat{D}$ is a differential operator (Hanslmeier & Dvorak 1984; Chambers 1999). We can now formally integrate Eq. (3) to obtain the general solution of f at time t + Δt: $$f(t+\Delta t)={\text{e}}^{\widehat{D}\Delta t}f(t)$$(4)

(e.g. Hanslmeier & Dvorak 1984; Yoshida 1993; Hernandez & Dehnen 2017; Wisdom 2018; Tamayo et al. 2019). Therefore, the solution of Eq. (3) is formally given by an ideal operator $\mathcal{D}={\text{e}}^{\widehat{D}\Delta t}$, which depends on the Hamiltonian $\mathcal{H}$ and advances f by one time step Δt. The exponential of an operator is defined as $${\text{e}}^{\widehat{D}\Delta t}={\sum }_{k=0}^{\infty }\frac{1}{k!}{\left(\widehat{D}\Delta t\right)}^k=1+\Delta t\widehat{D}+\frac{1}{2}\Delta t^2{\widehat{D}}^2+\dots$$(5)

(Hanslmeier & Dvorak 1984). However, $\mathcal{D}$ has no closed-form solution in general. Therefore, we make use of splitting the Hamiltonian into parts, each of which can be solved easily on its own and applied to approximate the true solution of the problem (e.g. Wisdom & Holman 1991; Duncan et al. 1998; Chambers 1999).

For example, if we split the Hamiltonian into two parts, such that $\mathcal{H}={\mathcal{H}}_{\text{A}}+{\mathcal{H}}_{\text{B}}$, we would get two operators  and $\widehat{B}$, $$\widehat{A}={\sum }_{i=1}^N\left(\frac{\partial }{\partial x_i}\frac{\partial {\mathcal{H}}_{\text{A}}}{\partial p_i}-\frac{\partial }{\partial p_i}\frac{\partial {\mathcal{H}}_{\text{A}}}{\partial x_i}\right),$$(6)

and a similar expression for $\widehat{B}$, such that $\mathcal{D}={\text{e}}^{(\widehat{A}+\widehat{B})\Delta t}$. However, Â and $\widehat{B}$ do not necessarily commute, which means that the relation $[\widehat{A},\widehat{B}]=\widehat{A}\widehat{B}-\widehat{B}\widehat{A}\ne 0$ in general. With this in mind, we can easily show that by expanding the exponential expression ${\text{e}}^{(\widehat{A}+\widehat{B})\Delta t}={\text{e}}^{\widehat{A}\Delta t}{\text{e}}^{\widehat{B}\Delta t}$ is only true to first order in Δt. Higher-order solutions can be constructed in the following way $${\text{e}}^{(\widehat{A}+\widehat{B})\Delta t}={\displaystyle \prod _{i=1}^k{\text{e}}^{c_i\widehat{A}\Delta t}}{\text{e}}^{d_i\widehat{B}\Delta t}+\mathcal{O}(\Delta t^{n+1})$$(7)

for a suitable set of coefficients (c_i, d_i) with i ∈ (1, …, k) (Yoshida 1993). For example, a second-order accurate scheme is obtained by setting c₁ = 1∕2 = c₂, d₁ = 1, and d₂ = 0, which results in $${\text{e}}^{(\widehat{A}+\widehat{B})\Delta t}={\text{e}}^{\frac{1}{2}\widehat{A}\Delta t}{\text{e}}^{\widehat{B}\Delta t}{\text{e}}^{\frac{1}{2}\widehat{A}\Delta t}+\mathcal{O}(\Delta t^3)$$(8)

(Yoshida 1993).

Splitting the Hamiltonian in kinetic and potential energy is often used for the general N-body problem, which allows us to construct the second-order accurate leap-frog algorithm. For the Solar System, however, where the Sun is the dominating body forcing the other bodies on Keplerian orbits, a different splitting is more beneficial. Using so-called democratic–heliocentric coordinates, which means heliocentric position and barycentric momentum, the Hamiltonian can be transformedto $$\mathcal{H}={\mathcal{H}}_{\text{A}}+{\mathcal{H}}_{\text{B}}+{\mathcal{H}}_{\text{C}},$$(9)

with $$\begin{array}{ll}{\mathcal{H}}_{\text{A}}\hfill & ={\sum }_{i=1}^N\left(\frac{p_i^2}{2m_i}-\frac{G{m}_i{m}_{\odot }}{\left|r_i\right|}\right),\hfill \\ {\mathcal{H}}_{\text{B}}\hfill & =-{\sum }_{i=1}^N{\sum }_{j=i+1}^N\frac{G{m}_i{m}_j}{|r_i-r_j|},\hfill \\ {\mathcal{H}}_{\text{C}}\hfill & =\frac{1}{2m_{\odot }}{\left({\sum }_{i=1}^N{p}_i\right)}^2.\hfill \end{array}$$

Here, r_i = x_i −x_⊙ is the heliocentric position of body i and p_i is its barycentric momentum; ${\mathcal{H}}_{\text{A}}$ is the Hamiltonian of the Kepler problem; ${\mathcal{H}}_{\text{B}}$ takes the interactions between the bodies, excluding the Sun, into account; and ${\mathcal{H}}_{\text{C}}$ is the kinetic energy of the Sun. Each of the Hamiltonians can be easily solved in the absence of the others. Efficient methods exist for solving the Kepler problem (Danby 1992; Rein & Tamayo 2015; Wisdom & Hernandez 2015), the interaction part reduces to calculating forces between the bodies, and ${\mathcal{H}}_{\text{C}}$ does not depend on the position, which makes it easy to solve as well. Therefore, we can use this splitting to construct a second-order accurate scheme in the form of Eq. (8) as $${\text{e}}^{(\widehat{A}+\widehat{B}+\widehat{C})\Delta t}={\text{e}}^{\frac{1}{2}\widehat{B}\Delta t}{\text{e}}^{\frac{1}{2}\widehat{C}\Delta t}{\text{e}}^{\widehat{A}\Delta t}{\text{e}}^{\frac{1}{2}\widehat{C}\Delta t}{\text{e}}^{\frac{1}{2}\widehat{B}\Delta t}+\text{O}(\Delta t^3)$$(11)

to integrate the motion of a body i (Chambers 1999), that is for f = (r_i, p_i). The first step is a kick for half of a time step owing to the interaction with the other bodies. Then, the body drifts from the change in kinetic energy of the Sun. Next, the body moves on a Keplerian orbit for Δt, followed by another drift and a kick.

2.2 Keplerian motion

Without mutual gravitational perturbation, bodies move on Keplerian orbits. The equivalent one-body Hamiltonian of Keplerian motion for body i is $${\mathcal{H}}_{\text{A}}^i=\frac{p_i}{2m_i}-\frac{G{m}_i{m}_{\odot }}{\left|r_i\right|},$$(12)

We use the method described in Mikkola & Innanen (1999) and Rein & Tamayo (2015), which uses the Gauss f- and g-functions to solve for it. The Gauss f- and g-function formalism is advantageous for efficiently solving the Kepler problem because it avoids computationally expensive transformations to orbital elements and problems arising with circular orbits (Rein & Tamayo 2015).

The key idea is that the position r_i and velocity v_i = p_i∕m_i of a body at time t = t₀ + Δt can be expressed in terms of the initial values at time t₀, r_i,0, and v_i,0 as $$\begin{array}{ll}{r}_i\hfill & =f{r}_{i,0}+g{v}_{i,0},\hfill \\ v_i\hfill & =\dot{f}{r}_{i,0}+\dot{g}{v}_{i,0},\hfill \end{array}$$

where f and g are scalar functions that only depend on time and ḟ and ġ are the time derivatives of these functions. A detailed description of how to obtain f and g, their time derivatives, and how to implement and use them to calculate the orbit can be found, for instance, in Danby (1992), Mikkola & Innanen (1999), Rein & Tamayo (2015), Wisdom & Hernandez (2015), and references therein.

2.3 Interactions between the bodies

The integration scheme of Eq. (11) constructed from the Hamiltonian splitting of Eq. (9) can be used as long as bodies are well separated. Because the Keplerian part dominates in that case, the error for integrating one time step is $~\mathcal{O}(\varepsilon \Delta t^3)$ with ɛ being the mass ratio of the bodies with respect to the Sun (Chambers 1999). However, when two bodies have a close encounter, the interaction term becomes comparable to ${\mathcal{H}}_{\text{A}}$ and the error per time step increases, which prevents us from obtaining accurate results. To avoid this problem, hybrid methods were developed to make sure that ${\mathcal{H}}_{\text{B}}$ remains small compared to the Keplerian part all the time. The key idea is to move the interacting bodies from ${\mathcal{H}}_{\text{B}}$ to ${\mathcal{H}}_{\text{A}}$ and to integrate those bodies numerically with a higher-order method instead of using a Kepler solver. We can move bodies between ${\mathcal{H}}_{\text{A}}$ and ${\mathcal{H}}_{\text{B}}$ by replacing Eqs. (10a) and (10b) with $$\begin{array}{ll}{\mathcal{H}}_{\text{A}}=\hfill & {\sum }_{i=1}^N\left(\frac{p_i^2}{2m_i}-\frac{G{m}_i{m}_{\odot }}{\left|r_i\right|}\right)\hfill \\ \hfill & +{\sum }_{i=1}^N{\sum }_{j=i+1}^N\frac{G{m}_i{m}_j}{|r_i-r_j|}\left[1-K(r_{ij})\right],\hfill \\ {\mathcal{H}}_{\text{B}}=\hfill & -{\sum }_{i=1}^N{\sum }_{j=i+1}^N\frac{G{m}_i{m}_j}{|r_i-r_j|}K(r_{ij}),\hfill \end{array}$$

where K is a function of the mutual distance of two planetesimals, r_ij = |r_i −r_j|, which acts as a switch to move terms between the two Hamiltonians (e.g. Duncan et al. 1998; Chambers 1999; Rein et al. 2019). The corresponding operators then become $$\begin{array}{ll}\widehat{A}=\hfill & {\sum }_{i=1}^N\{{\sum }_{j\ne i}^N\frac{\partial }{\partial r_i}\frac{p_i}{m_i}-\frac{\partial }{\partial p_i}G{m}_i{m}_{\odot }\frac{r_i}{|r_i{|}^3}\hfill \\ \hfill & -\frac{\partial }{\partial p_i}G{m}_i{\sum }_{j\ne i}^N{m}_j\frac{r_i-r_j}{|r_i-r_j{|}^3}\left[1-L(r_{ij})\right]\},\hfill \\ \widehat{B}=\hfill & {\sum }_{i=1}^N\left\{\frac{\partial }{\partial p_i}G{m}_i{\sum }_{j\ne i}^N{m}_j\frac{r_i-r_j}{|r_i-r_j{|}^3}L(r_{ij})\right\},\hfill \end{array}$$

where we defined a new function L = K − r_ij∂K∕∂r (Rein et al. 2019). We chose L to be a Heaviside step function, which is defined such that $$L=\{\begin{array}{ll}0\hfill & r_{ij}<d_{\text{ce}},\hfill \\ 1\hfill & r_{ij}>d_{\text{ce}},\hfill \end{array}$$(19)

where d_ce is the critical distance between bodies for classifying it as a close encounter (see Sect. 2.3.1). It is now easy to see that without a close encounter, Â reduces to the operator for the Kepler problem, which can be easily solved using a Kepler solver, and $\widehat{B}$ adds small perturbations. For a close encounter, on the other hand, Â contains an additional term and we resort to numerically integrating a few-body problem using a Bulirsch–Stoer algorithm (Press et al. 1992).

So far, we have only summarised the basic ideas of hybrid–symplectic integrators for the planetary N-body problem. However, a bottleneck of simulating a large number of bodies is their mutual interaction (i.e. ${\mathcal{H}}_{\text{B}}$). Because forces between everybody pair must becalculated, which is an N² calculation, simulating many bodies would result in unfeasibly long computation times. We have approached this problem by dividing the bodies into two classes, namely planets and planetesimals, and by using an efficient close-encounter detection method.

Planets

We cannot ignore the long-range forces for planets, and therefore treat them in full N-body fashion. We solve the full Hamiltonian, which means that we not only resolve close encounters, but also take the distant perturbations into account. Therefore, our method is equivalent to a hybrid–symplectic integrator, such as mercury (Chambers 1999).

Planetesimals

We takeplanetesimals as low-mass bodies for which we chose to ignore the long-range interaction and only take close encounters into account. This approach is justified because of the low total mass of the planetesimals compared to the Sun and the fact that significant orbital changes only occur if the approach is very close, typically within a few times the Hill radius (Petit & Henon 1986; Greenberg et al. 1991). The strategy is thus to evolve planetesimals on their Keplerian orbit if there is no close encounter, as described in the previous section, and to evolve the few-body problem of the Sun plus nearby planetesimals otherwise. Because planetesimals grow by mutual collisions and pebble accretion, they naturally evolve into planets. To account for this, we promote planetesimals to planets, when they reach a minimum mass, which we take as 10⁻² M_⊕.

2.3.1 Distance for close encounter

If the distance r_ij between two bodies i and j is less than a critical distance d_ce, they have a close encounter. This critical distance is typically a few times the Hill radius: $$R_{\text{h}}=r{\left(\frac{m}{3m_{\odot }}\right)}^{1/3}.$$(20)

Here r is the radial heliocentric distance and m is the mass of the body. Within a sphere of radius R_h, the gravity of the body dominates over the solar gravity for the motion of a test particle in the co-rotating frame. As a consequence, the mutual gravitational interaction between the two bodies results in strong scattering and significant changes of their respective orbits (Henon & Petit 1986; Petit & Henon 1986; Ida 1990). Distant encounters, where r_ij ≫ R_h, only lead to weak scattering with small orbital changes (Weidenschilling 1989; Hasegawa & Nakazawa 1990). The same applies for encounters with r_ij ≪ R_h, which result in horseshoe orbits. In these types of encounters, eccentricity and inclination hardly change (Henon & Petit 1986; Hasegawa & Nakazawa 1990). Therefore, the encountering body can enter the Hill sphere of the other object and will be strongly scattered only if r_ij ~ R_h.

In our close-encounter method, d_ce in units of the Hill radius, is a parameter that is set for each body individually. Numerical simulations of planetesimal encounters showed that orbital changes are strong if 0.8 R_h≲r_ij≲2.5 R_h (Petit & Henon 1986; Greenberg et al. 1991). Therefore, a typical value to use in our method would be d_ce ~ 3 R_h. We do not use a lower value, which would result in missing close encounters. On the other hand, a higher value catches more distant encounters, and thus we typically choose a value of d_ce ~ 3…5 R_h, which reproduces the dynamics of the planetesimals to a high degree, as we show in Sect. 3.

2.3.2 Cell-list approach

A computationally efficient way to detect close encounters is with a cell list, an approach that is frequently used in molecular dynamics simulations (Frenkel & Smit 2002) and dust particle collisions (Johansen et al. 2012). A grid divides the simulation domain into n_cell individual cells. The grid size is chosen in a way that the cell dimensions are larger than the typical distance associated with a close encounter, that is ~d_ce. For each cell, we assign a body as head of the list and link it to the next one in this cell, which creates a linked list of nearby bodies (see Algorithm 1 for the algorithm in pseudo-code).

To identify a close encounter, we loop through the cell list. The criterion for a close encounter is that the distance between two bodies is less than d_ce. Bodies residing close to the cell boundaries potentially undergo a close encounter with an object from a neighbouring cell. Therefore, looping over the neighbour cells is necessary to account for these cases as well. There are at most 27 cells to loop over. To check for close encounters the orbits of the bodies are approximated as straight lines connecting their initial positions and the new positions they would have if advanced along their unperturbed Keplerian orbit. If the minimum separation between these lines is less than d_ce, the encounter is classified as a close encounter and the body pair is added to a list (see Algorithm 2 for the algorithm in pseudo-code) that is used afterwards to integrate the encounter (see Sect. 2.3.4). Because creating the cell list scales linearly with the number of bodies and because the average number of bodies per grid cell is lower than the total number in the simulation, the close-encounter detection is faster than a simple loop over each body pair, which scales as N², especially for a large number of bodies.

2.3.3 Grouping into close-encounter groups

If more than two bodies have a close encounter, which may occur for high surface density of planetesimals or for large d_ce, the bodies need to be grouped together. To do this we follow the grouping algorithm of Grimm & Stadel (2014) (see Algorithm 3 for the algorithm in pseudo-code).

2.3.4 Close-encounter integration

Having a list of close encounters, we numerically integrate the motion of these bodies. We use the Bulirsch–Stoer (BS) method (Press et al. 1992), which is also implemented, for example, in mercury (Chambers 1999). The basic idea of the BS method is to integrate the motion of the body for one time step Δt by sub-dividing Δt into smaller bits, and integrating the smaller bits with a modified midpoint rule until the results change by less than a given tolerance parameter.

The gravitational acceleration acting on body i according to Eq. (17) is $$a_i=-G{m}_{\odot }\frac{r_i}{|r_i{|}^3}-{\sum }_{j\ne i}^{N_{\text{ce}}}G{m}_j\frac{r_j-r_i}{{(r_{ij}^2+b^2)}^{3/2}}.$$(21)

The first term is the gravity of the Sun and the second term is the acceleration from the other bodies. The summation is over all bodies N_ce that are involved in the close encounter. We use a gravitational softening parameter b to avoid high acceleration in very close encounters. The softening parameter is set to be the sum of the physical radii of planetesimals i and j. This choice for the softening is reasonable because the two planetesimals would undergo a physical collision if r_ij ≤b. If additional forces are included, they are simply added to a_i. We provide two different ways to integrate the close encounters: a sequential approach, which uses the list of close-encounter pairs, and a groupapproach, which uses close-encounter groups.

Sequential approach

The sequential approach is appropriate if planetesimals have at most one encounter per time step. We therefore have N_ce = 1, which means that we solve the three-body problem of the Sun and two other bodies. To avoid correlations in rare cases when a body might encounter more than one other body, we randomise the close-encounter list.

Group approach

The group approach is useful if the bodies encounter more than one other body per time step. In this case we have to treat the encounter group as a small N-body system with N = N_ce > 1 bodies. For large groups with N_ce ≫ 1, the group approach becomes comparable to a standard hybrid–symplectic scheme. Therefore, large groups should be avoided to obtain the full potential of the close-encounter method.

2.4 Grid size and time step

Using a cell list requires a grid. We construct the grid in cylindrical coordinates (r, ϕ, z), which is especially useful for simulating rings of planetesimals. The grid cells need to be large enough to include the close-encounter region of a body. This sets the minimum size for the grid cells to be Δr > d_ce, Δϕ > d_ce∕r, and Δz > d_ce. If the grid cell were smaller, bodies would potentially undergo close encounters with bodies outside the neighbouring cells. On the other hand, the grid cells need to be small enough, such that the average number of planetesimals per grid cell is ≪ N. Otherwise, the cell list approach loses its scaling advantage over a conventional N² -pair search.

Using a cell list also imposes a constraint on the maximum time step that we can use in the simulation. As already pointed out, only objects in the current and neighbouring cells are checked for close encounters. Therefore, bodies should not move over more than one grid cell per time step to avoid missing close encounters with bodies outside the neighbouring cells. To adjust the time step, we calculate the crossing times in each direction for all bodies. The current allowed time step is then the minimum of these crossing times $$\Delta t={\min }_i\left(\frac{\Delta r}{|v_{r,i}},\frac{\Delta \varphi }{\left|v_{\varphi ,i}\right|},\frac{\Delta z}{\left|v_{z,i}\right|}\right),$$(22)

where i runs over all bodies and the radial, azimuthal, and vertical velocities of the bodies are v_r,i, v_ϕ,i, and v_z,i, respectively.Because the time step is not constant throughout the simulation, our method is not symplectic.

2.5 Additional forces

To be able to include gas drag or migration, we allow for additional forces in our method. We do so by adding the acceleration in the interaction part that comes from ${\mathcal{H}}_{\text{B}}$ (Chambers 1999; Tamayo et al. 2019).

3 Results

We test the close-encounter method by simulating the viscous stirring of a ring of planetesimals and compare the results to N-body simulations done with mercury (Chambers 1999) and to the semi-analytical toy model of Ohtsuki et al. (2002). We use the same setup as in Ohtsuki et al. (2002) to test the viscous stirring of planetesimals. A ring of 1000 planetesimals with surface density of 10 g cm⁻² is located at 1 au. For a single planetesimal mass, we use a mass of m = 10²⁴ g. For a bimodal mass distribution, we split the 1000 planetesimals into 800 with mass 10²⁴ g (component 1) and 200 with mass 4 × 10²⁴ g (component 2). The initial orbits have eccentricities and inclinations following Rayleigh-distributions with e_rms = 10⁻⁴ and i_rms = 5 × 10⁻⁵. The phase angles, argument of perihelion ω, longitude of ascending node Ω, and the mean anomaly M are drawn from a uniform distribution in (0, 2π). We then integrate the planetesimals for a total time of 10⁴ yr (Ohtsuki et al. 2002, only integrated for 10³ yr). In mercury, we use the hybrid-symplectic method, a time step of 8 days, and a Bulirsch-Stoer tolerance of 10⁻¹². For our model and for the comparison runs with mercury, we show averages over five runs with varying random seed for setting up the initial conditions.

3.1 Equal-mass planetesimals

Figure 1 shows the viscous stirring of the equal-mass planetesimals. For the close-encounter distance, we used 10 Hill radii in our model. Every encounter with a minimum distance less than that is integrated with the Bulirsch–Stoer part of our code. In this test case we use the group approach. Our method reproduces the results of the full N-body simulation and the toy model sufficiently well, showing that the close-encounter approach works well for this particular case. We have also tested the sequential approach for this specific test case using 5 Hill radii for the close-encounter distance and found no noticeable difference other than a speed-up of ~ 1.5.

3.2 Unequal-mass planetesimals

Figure 2 shows the viscous stirring and dynamical friction of the two planetesimal populations. The results of our model, the close-encounter method, match the mercury and toy-model results. Because the gravitational kicks on the lower mass bodies (component 1) from the more massive planetesimals (component 2) are stronger than the kicks on higher mass bodies from the less massive planetesimals, on average e_rms and i_rms are higher for component 1. This effect is called dynamical friction which drives the system towards energy equipartition and a mass-dependent velocity dispersion. However, equipartition is not exactly reached because viscous stirring also increases e_rms and i_rms of both components.

3.3 Embedded planet

In this test case shown in Fig. 3, we simulate the viscous stirring of a ring of planetesimals with an embedded planet. To see the effects of viscous stirring and dynamical friction, we initialise the planet with e = 10⁻² and i = 5 × 10⁻³, which are 1000 times larger than the rms values of the planetesimals, which we initialised here with e_rms = 10⁻⁵ and i_rms = 5 × 10⁻⁶. We see that theorbit of the planet is rapidly damped towards very low eccentricity and inclination. This is the effect ofdynamical friction through the swarm of planetesimals. On the other hand, the planetesimals are excited into highly eccentric and inclined orbits. The results of our method match the results of N-body simulations and the toy model in this case as well.

Fig. 1 Viscous stirring for equal-mass planetesimals. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of equal-mass planetesimals with mass 1 × 1024 g. The ring is centred at 1 au and the surface density is 10 g cm−2. Shown are the results of our model: the close-encounter approach (orange), a full N-body simulation with mercury (grey), and the semi-analytic toy model of Ohtsuki et al. (2002; black dashed).

3.4 Dependence on the close-encounter distance

The critical distance for close encounters d_ce controls the contribution of distant encounters to the stirring of the planetesimals. The stirring increases for larger d_ce because more distant encounters are included. We investigate how the choice of d_ce affects our method. Therefore, we run the equal-mass planetesimal setup for different values of d_ce ranging from 1 R_h to 10 R_h.

Figure 4 shows the results when using the sequential approach in our method. For d_ce = 1 R_h, we miss many close encounters and therefore obtain little stirring of the planetesimals. Increasing it to d_ce ~ 5−7 R_h produces a close match with the toy model and mercury. Increasing d_ce further results in increased stirring at times t≳3000 yr. The reason is that the sequential approach breaks down. Because d_ce is large, the close-encounter regions of planetesimals overlap and they encounter more than one other body per time step. In this type of situation, integrating all close encounters sequentially gives the wrong results.

Figure 5 shows the same experiment of increasing the close-encounter distance, but this time we used the group approach of our method. The results are similar for d_ce = 1 R_h, for which we miss too many encounters, but the solution converges towards the toy model and mercury for larger values. Therefore, the sequential approach is safe to use as long as there is only one encounter per planetesimal per time step. We discuss this limitation more in Sect. 4.

Fig. 2 Viscous stirring for bimodal mass distribution. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of planetesimals with bimodal mass distribution. There are 800 planetesimals with mass 1024 g (component 1, m1) and 200 with mass 4 × 1024 g (component 2, m2). The ring is centred at 1 au and the surface density is 10 g cm−2. Shown are the results of our model (orange), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black dashed). Component 1 and 2 are indicated with solid and dotted lines, respectively.

3.5 Computational time and scaling

Our aim is to develop a computationally fast method. Therefore, we look at the computational time and the scaling of our close-encounter method with number of bodies N for the setup of a ring of equal-mass planetesimals. Figure 6 shows the result. We measured the computational time for the sequential approach (sequential; orange), the group approach (grouping; blue), and, for comparison, a full N-body run with 1000 bodies with mercury.

We can see the advantage of our method in Fig. 6. Simulating the dynamical evolution of 1000 equal-mass planetesimals for 10⁴ yr takes about 15 minutes with the close-encounter method. The same setup with full N-body with mercury requires 10 hr, and more for a larger number of bodies. However, in 10 hr, the close-encounter method (both approaches, sequential and group) can integrate the dynamical evolution of 10 000 bodies. The close-encounter method is hence faster without losing the ability to obtain correct results. Full N-body methods typically scale with the number of bodies as N². Our sequential approach scales linearly, which is expected because N² operations are avoided. Setting up the cell list is an operation of order N and the integration of close encounters reduces to a sequence of three-body integrations. The group approach does not scale linearly with N, but better than quadratic as ∝ N^1.8. The execution time is much faster than for conventional N-body methods in both cases. This allows simulating a larger number of bodies than with full N-body methods.

Fig. 3 Viscous stirring and dynamical friction for a planet embedded in a ring of planetesimals. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of planetesimals with an embedded planet. There are 1000 planetesimals with mass 1024 g and 1 planet with mass 1026 g. The ring is centred at 1 au and the surface density is 10 g cm−2. The planet is initially at 1 au. Shown are the results of our model (orange), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed). Planetesimals and the planet are indicated with solid and dotted lines, respectively.

3.6 Energy and angular-momentum conservation

We look at the energy and angular-momentum conservation of our close-encounter method by studying the relative change of energy and angular momentum with time and compare it to a full N-body simulation with mercury. We again use the setup for the equal-mass planetesimal ring. The relative change of energy E is |(E − E₀)∕E₀| (and likewise for the angular momentum L), where E₀ is the total energy at the beginning and E is the energy at time t.

Figure 7 shows the results. Energy and angular momentum show a steady increase over the simulation time of 10⁴ yr. We first look at the angular momentum. Regardless of the approach (sequential or group), the change in angular momentum is the same. Compared to mercury, we achieve the same accuracy. The situation looks slightly different for the energy. The hybrid-symplectic integrator mercury conserves energy because mercury is a symplectic method and solves a Hamiltonian that is close to the true N-body Hamiltonian of the system at any time (Chambers 1999). Our close-encounter method, however, has two significant differences. Firstly, we ignore the interaction terms for planetesimals unless they undergo a close encounter, and thus only for planets is our method comparable to mercury. Therefore, each time there is a close encounter between planetesimals, our method changes the Hamiltonian (i.e. the energy of the system) by adding a small perturbation owing to the interaction. Secondly, our method does not use a constant time step, which breaks the symplecticity. These errors accumulate over time, as seen in Fig. 7, which results in the steady increase in |(E − E₀)∕E₀|.

Fig. 4 Viscous stirring of equal-mass planetesimals for increasing close-encounter distance using the sequential approach. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) for the same setupas for the equal-mass planetesimal run. Shown are the results of our model for different values of dce (coloured lines), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed).

3.7 Collisions

We apply the close-encounter method to study the collisional growth of planetesimals located in a ring at 1 au. The setup is similar to the one studied in Kokubo & Ida (1998) and Kokubo & Ida (2000). Four thousand planetesimals with mass 3 × 10²³ g are placed at 1 au in a ring that is 0.085 au in width. Initial eccentricities and inclinations are Rayleigh-distributed with e_rms = 2i_rms = 2 × 10⁻³. Gas drag is included as described in Appendix C.1. Bodies that reach a mass of 10⁻² M_⊕ are promoted to planets. We integrate for 0.5 Myr. Figure 8 shows snapshots of eccentricity, inclination, and mass as a function of semi-major axis at three different times.

Starting with equal masses, collisions result in the formation of runaway bodies within ~ 5 × 10⁴ yr. In the later evolution, these bodies grow the fastest and separate from the bulk of the planetesimals. This can be seen in the bottom panels of Fig. 8. Figure 9 shows the mean and maximum mass of the system. There the runaway growth can be seen more clearly. The mean mass increases by about a factor of two. The maximum mass, however, increases more rapidly withtime and separates from the mean mass, as is characteristic for runaway growth. By the end of the simulation at 0.5 Myr, the maximum mass increased by a factor of ~100, whereas the mean mass increased by a factor of ~2. This is consistent with Kokubo & Ida (1998) and Kokubo & Ida (2000). The few massive bodies also start to stir the surrounding planetesimals, which is visible in the a–e and a–i plots (Fig. 8). The lower mass planetesimals are excited to eccentricities and inclinations of ~ 0.02. The more massive bodies, on the other hand, retain low eccentricities and inclinations because of dynamical friction. At 0.5 Myr, a total number of 1484 planetesimals (m < 10⁻² M_⊕) and 2 protoplanets (m > 10⁻² M_⊕) remain. Simulating 0.5 Myr of evolution took about 2.3 days. In comparison, with mercury it took us about 3 weeks to simulate the same setup for 0.1 Myr.

Fig. 5 Viscous stirring of equal-mass planetesimals for increasing close-encounter distance using the group approach. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) for the same setupas for the equal-mass planetesimal run. Shown are the results of our model for different values of dce (coloured lines), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed).

Fig. 6 Scaling of computational time with number of bodies. The setup is the same as in Sect. 3.1. The surface density is kept constant at 10 g cm−2. The number of planetesimals increases from 1000 to 10 000. Shown are the total computational time to integrate 104 yr with thesequential approach (sequential; orange), the group approach (grouping; blue), and, for comparison, with mercury (grey diamond). Also shown is a power-law fitting to find the scaling with N for our model (blue and orange lines).

Fig. 7 Energy and angular momentum. The figure shows the time evolution of the relative change of energy (top panel) and angular momentum (bottom panel) for equal-mass planetesimals. Shown are the results obtained with the sequential approach (sequential; orange) and the group approach (grouping; blue) of our method and compared with mercury (grey).

3.8 Pebble accretion

As another test application, we apply the close-encounter method to planet formation with collisions and pebble accretion in the giant planet region. Therefore, we distribute six embryos of mass 10⁻² M_⊕ and 444 planetesimals of mass 10⁻⁴ M_⊕ in a ring between 5 au and 10 au. A total flux of pebbles of 100 M_⊕ Myr⁻¹ enters the system from the outer disc and is potentially accreted by the bodies. The pebbles have a constant Stokes number of 10⁻² and the ratio of the pebble to the gas scale-height is 0.1. We use a prescription for pebble accretion that is based on the work by Ormel & Klahr (2010), Lambrechts & Johansen (2012), Johansen & Lambrechts (2017), Ormel (2017), and Lambrechts et al. (2019). Details can be found in the Appendix B. As bodies grow in mass they start to migrate. We therefore implemented an additional force for type I migration following Cresswell & Nelson (2008) (see Appendix C.2). The gas disc is a minimum mass solar nebula (Weidenschilling 1977b; Hayashi 1981, see Appendix C.1). We note that this is only a demonstration of the capability of the close-encounter method without the goal of providing a quantitative study of pebble accretion.

Figure 10 shows snapshots of eccentricity, inclination, and mass at times 0 Myr, 1.5 Myr, and the final outcome after 3 Myr. The six embryos grow efficiently by pebble accretion to a few Earth masses within 3 Myr (see Fig. 11). Pebble accretion proceeds in a runaway fashion, as seen in Fig. 12. The embryos grow very quickly, but the majority of the planetesimals remain at low masses. The exponential growth comes from the fact that the embryos grow in the 3D regime (Lambrechts et al. 2019). Collisions with planetesimals do not play a significant role. As the embryos grow bigger, they start to migrate inward by type I migration. The most massive body reaches pebble isolation mass of ~ 10 M_⊕ and stops migrating at the inner edge of the disc, which in this case is located at 0.3 au. We do not form cold gas giants in our simulations. The assumed pebble flux of 100 M_⊕ Myr⁻¹ corresponds to a nominal gas accretion-rate of the disc of ~3 × 10⁻⁸ M_⊙ yr⁻¹ (Hartmann et al. 2016) for a dust-to-gas ratio of 0.01 and a nominal metallicity. Slightly higher accretion rates and metallicities are needed to form cold gas giants within the disc lifetime (Johansen et al. 2019; Bitsch & Battistini 2020).

Figure 10 also shows that the massive embryos efficiently stir the surrounding planetesimals, which renders collisional growth inefficient. The combined effect of eccentricity and inclination damping due to the gas and dynamical friction with the lower mass planetesimals keeps the embryos on orbits with low eccentricities and inclinations throughout the simulation.

One simulation took about 3 days to complete 3 Myr of evolution.This is considerably longer than the runs with only collisions. However, with pebble accretion more bodies potentially reach the threshold mass of 10⁻² M_⊕ and need to be considered planets in our close-encounter method, which is computationally more expensive. Additionally, type I migration reduces the maximum time step allowed in the simulation because of the grid which is necessary for creating the cell list, which requires more single integration steps.

Fig. 8 Collisional growth of planetesimals at 1 au. Shown are scatter plots for eccentricity (first row), inclination (second row), and mass (third row). 4000 planetesimals are located at 1 au in a ring of width 0.085 au. The three columns show snapshots at different times. Planetesimals that grow more massive than 10−3 M⊕ are shown as coloured filled circles. All other bodies are grey open circles. The colour and size of the symbols scale with the instantaneous mass.

Fig. 9 Mean and maximum mass of planetesimals. Shown are the mean (orange line) and maximum (blue line) mass of the planetesimals as a function of time. Both lines are averages of five runs with different random realisations of the initial conditions.

Fig. 10 Pebble accretion and collisions in the giant planet region. Shown are the scatter plots for eccentricity (first row), inclination (second row), and mass (third row). A total mass of 0.1 M⊕ is distributed amongst 450 planetesimals between 5 au and 10 au. The three columns show snapshots at 0 Myr, 1.5 Myr, and at the final time of 3 Myr. Bodies that are more massive than 10−2 M⊕ are shown as coloured filled circles, their growth tracks are shown by the coloured lines, and the centre of the circle is indicated by a black dot. All other bodies are grey open circles. The colour and size of the symbols scale with the instantaneous mass.

Fig. 11 Cumulative number of bodies for pebble accretion in the giant planet region. Shown are N(> m), the cumulative number of bodies with mass higher than m, for the initial distribution and for the distribution at times 0, Myr, 1 Myr, 2 Myr, and 3 Myr. Each line is the average of five runs with randomised initial configuration of the bodies.

Fig. 12 Mass evolution of planetesimals for collisions and pebble accretion in the giant planet region. Each line shows the growth history of a planetesimal or embryo. Bodies that grow more massive than 10−2 M⊕ are shown as coloured lines. The colour scales with mass and is the same as in Fig. 10. The grey lines show all other bodies. The endpoint of the simulation for each body is indicated with a circle. Because the bodies do not grow significantly before 104 yr, only the later times are shown.

4 Discussion

As shown in Fig. 6, the sequential close-encounter approach scales almost linearly with the number of bodies. The group approach scales as N^1.8, which is only slightly better than quadratically. However, for both methods the computational time is significantly lower than for mercury. Therefore, it is possible to simulate a large number of planetesimals in reasonable times. The speed-up occurs not only because we ignored the long-range interactions for planetesimals, and thus the N² process of calculating the forces between all planetesimals, but also because we chose a more efficient detection method for close encounters. The test cases, which were compared toN-body simulations with mercury and semi-analytic toy models, show that our approach sufficiently reproduces the dynamics of the planetesimals. When applied to more physical cases (i.e. including planets, collisions, and pebble accretion) our method becomes a faster alternative to conventionalN-body methods.

The sequential treatment of close encounters increases the speed of the method significantly; however, one has to carefully check the validity of this approximation. The sequential approach breaks down when planetesimals have more than one encounter per time step. This can occur either in a very dense system or when the close-encounter regions of many planetesimals overlap. Because the physical sizes of planetesimals are small, the latter situation is much more likely. The group approach and conventional N-body codes do not suffer from this limitation, but integrating the encounters sequentially will give incorrect results.

We now derive an approximate criterion for the validity of the sequential approach. Therefore, we calculate the ratio of the mean-free path $\lambda ={(n\sigma )}^{-1}$ of the planetesimals, where n is their number density and σ the relevant close-encounter cross section, to the close-encounter distance d_ce. Substituting the number density by Σ∕2mH, where 2H is the thickness of the planetesimal disc, the mean-free path is given by $$\lambda =\frac{2mH}{\sigma \Sigma },$$(23)

where m is the mass of the planetesimal, H is their scale height, and Σ is the surface density. The scale height can be expressed in terms of the root mean square inclination as H = r i_rms, the cross section is $\sigma =\pi d_{\text{ce}}^2$, and we obtain the ratio $$\begin{array}{ll}\frac{\lambda }{d_{\text{ce}}}\hfill & =\frac{6m_{\odot }{i}_{\text{rms}}}{\pi \Sigma r^2{(d_{\text{ce}}/R_{\text{h}})}^3},\hfill \\ \hfill & =3{\left(\frac{\Sigma }{10\text{g}{\text{cm}}^{-2}}\right)}^{-1}{\left(\frac{r}{1\text{au}}\right)}^{-2}\left(\frac{i_{\text{rms}}}{5\times {10}^{-5}}\right){\left(\frac{d_{\text{ce}}}{3R_{\text{h}}}\right)}^{-3}.\hfill \end{array}$$(24)

Figure 13 shows the ratio of mean-free path and close-encounter distance for different values of d_ce assuming that the surface density varies as $\Sigma =10\text{g}{\text{cm}}^{-2}{(r/\text{au})}^{-1}$. Here, we can see why Fig. 4 shows the increased stirring for large d_ce. The ratio λ∕d_ce is much lower than unity for d_ce≳5 R_h. This means that in our simulations most of the planetesimal close-encounter regions, which are defined as the spheres of radius d_ce around the bodies, start to overlap. In this case, however, integrating the close encounters sequentially results in stronger stirring because our method evolves the dynamics of the planetesimals as a sequence of kicks without taking the gravitational pull from the other bodies into account. This results in a random-walk-like increase in the eccentricities and inclinations. This can be seen in Fig. 4 where at times ≳3000 yr and for large d_ce, e_rms scales with the square root of time.

Fig. 13 Ratio of mean-free path to close-encounter distance as a function of heliocentric distance. The solid lines represent λ∕dce for different values of dce. For the dashed line, this ratio is unity, indicating where the sequential approach breaks down. A value of Σ0 = 10 g cm−2 at 1 au and irms = 5 × 10−5 was used. Surface density is assumed to be a power law $\Sigma ={\Sigma }_0{(r/\text{au})}^{-1}$.

5 Summary

The dynamics of planetesimals is important for planet formation because it affects how efficiently they grow to bigger bodies. However, the total number of planetesimals in protoplanetary discs is high, and to resolve the dynamics the gravity between all of the bodies would need to be taken into account, which makes it computationally expensive. Therefore, we developed a close-encounter method for simulating a large number of planetesimals. Compared to conventional N-body methods which are designed to simulate planetary dynamics with very high accuracy but low N, our method focuses on efficiently evolving many interacting small bodies under the assumption that distant encounters can be ignored. Forrings of planetesimals this is justified because significant changes in their orbits only occur during close encounters. With this assumption, the close-encounter method provides the following advantages compared to more detailed N-body methods:

• Our method is conceptually simple. We employ the same ideas as were developed for symplectic N-body integrators. Thanks to this, the numerical integration reduces to solving the Kepler problem and a sequence of few-body problems;

• The method can handle planets as well, which makes it applicable for planet formation simulations;

• Our method is fast. We achieve a significant speed-up through an optimised close-encounter detection using a cell list (see Sect. 2.3);

• It scales linearly when encounters can be integrated sequentially and better than N² when encounters are grouped together (see Sect. 3.5).

However, there are also limitations to this method:

• The Hamiltonian that is used here describes single-star systems. Therefore, our method is not suitable for simulating binary systems. We plan to generalise the method to handle binaries in future work;

• The sequential treatment of close encounters is only valid if the close-encounter regions are not overlapping (see Sect. 4);

• If planets are included, their number should remain low. Otherwise our method becomes less efficient, comparable to a conventional N-body method, because the code will spend most of the time in calculation interactions between planets (i.e. ${\mathcal{H}}_{\text{B}}$). For the ring of equal-mass planetesimals (see Sect. 3.1) the interactions between planets and all other bodies becomes the most expensive process if we turn more than 25 of the 1000 planetesimals into planets. For ≳100 planets, the interaction part accounts for more than 50% of the total computational time. We therefore recommend limiting the number of planets to ≲100;

• The cell list requires a spatial grid, which imposes a time step constraint for the code (see Sect. 2.4). The time step is not constant, but depends on the fastest particle in the smallest grid cell;

• The method is not symplectic because we do not use a constant time step and because we ignore interactions other than close encounters between planetesimals. Therefore, energy is not conserved which might be a problem for long-term integrations (see Sect. 3.6).

Acknowledgements

We thank the anonymous referee for a thorough review that helped to significantly improve the original manuscript. AJ thanks the Swedish Research Council (grant 2018-04867), the Knut and Alice Wallenberg Foundation (grant 2017.0287) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) for research support.

Appendix A Algorithms

Algorithm 1 describes the algorithm for setting up the list in pseudo-code. The array head of length n_cell stores (for each grid cell) the index of the head. The array list of length N and initialised with 0 contains (for each body) the index of the next body in that cell.

Algorithm 1 Construct cell list

Algorithm 2 shows how to loop through the cell list. The body pairs that match the criterion for a close encounter are added to a list.

Algorithm 2 Detect close encounter

Algorithm 3 Grouping of close encounters

Algorithm 3 shows how to group bodies in close-encounter groups. The array link is initialised as link[i] = i for i = 1, …, N. Having the list of close-encounter pairs, a new array of links of length N is created which contains (for each body) the lowest index of the interacting planetesimals. The indices of the encountering bodies are i and j. We set the corresponding value of link via $$\begin{array}{ll}link[i]\hfill & =\min (link[i],link[j]),\hfill \\ link[j]\hfill & =\min (link[j],link[i])\hfill \end{array}$$

(Grimm & Stadel 2014). Bodies that undergo a close encounter are now linked to the body with the lowest index in their group. Next we assign a consecutive group index to the close-encounter groups. We first initialise an array group such that group[i] = i for i = 1, …, N and a counter for the total number of groups n_groups. In the second step we loop through all bodies, increase the counter if i = link[i], and set the array entry to $$group[i]=\min (group[link[i]],n_{\text{groups}}).$$(A.3)

In this way we map each body to its group and have a consecutive group index. Finally, we construct a matrix cegroup that contains (for each group) the indices of the bodies in that group. Another array cenumb contains the total number of bodies within each group.

Appendix B Pebble accretion

In this section, we briefly describe our prescription for pebble accretion. Pebble accretion is the process by which a planetesimal grows through the accretion of pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Johansen & Lambrechts 2017; Ormel 2017). Pebbles are aerodynamically coupled to and carried along by the gas. When the gravitational pull by the planetesimal is strong enough, pebbles are accreted. This translates into a timescale criterion: $$t_{\text{def}}=t_{\text{fric}}$$(B.1)

(Johansen & Lambrechts 2017; Ormel 2017; Lambrechts et al. 2019). Here, t_def is the deflection timescale and t_fric is the friction timescale of the pebble.

A pebble approaching the planetesimal at a distance r_acc with velocity Δv is pulled towards the planetesimal on a timescale of $$t_{\text{def}}=\frac{\Delta v{r}_{\text{acc}}^2}{Gm},$$(B.2)

where m is the mass of the planetesimal. The friction timescale is the time it takes for the pebble to adjust to the gas velocity. The relative velocity between pebble and planetesimal, Δv, can be decomposed into a radial, azimuthal, and vertical component: $$\begin{array}{ll}\Delta v_r\hfill & =v_r+\frac{2\eta {\tau }_{\text{f}}{v}_{\text{K}}}{1+{\tau }_{\text{f}}^2},\hfill \\ \Delta v_{\varphi }\hfill & =v_{\varphi }-v_{\text{K}}+\frac{\eta v_{\text{K}}}{1+{\tau }_{\text{f}}^2}+\frac{3}{2}{\Omega }_{\text{K}}{r}_{\text{acc}},\hfill \\ \Delta v_z\hfill & =v_z.\hfill \end{array}$$

Here, τ_f = t_fricΩ_K is the Stokes number of the pebble and η is the pressure gradient of the gas. In the radial direction, the relative velocity is given by the difference between the radial velocity of the planetesimal and the radial drift velocity of the pebble. Ignoring settling, the relative velocity in vertical direction is simply the vertical velocity of the planetesimal. In the azimuthal direction, we have to take the difference between the azimuthal velocity of the planetesimal, the azimuthal drift of the pebble, and the Keplerian shear, which is expressed in the third term of Δv_phi. By using this expression for Δv, the eccentricity and inclination of the planetesimal are taken into account and we can write $$\Delta v=\sqrt{\Delta v_r^2+\Delta v_{\varphi }^2+\Delta v_z^2}$$(B.4)

(Lambrechts et al. 2019).

Because Δv depends on the accretion radius, we solve the criterion Eq. (B.1) iteratively to obtain the accretion radius (Lambrechts et al. 2019). The accretion radius according to Eq. (B.1) is obtained with the assumption that the encounter timescale is long enough for the pebble to settle towards the planetesimals due to the action of gas drag. Therefore, we have to check whether this is the case by comparing the encounter timescale to the friction timescale of the pebble (Ormel 2017). The encounter timescale, t_enc, is the time it takes for the pebble to pass the gravitational reach of the planetesimal and can be expressed as $$t_{\text{enc}}=\frac{2r_{\text{acc}}}{\Delta v}$$(B.5)

(Ormel 2017). We calculate t_enc and compare it to t_fric. If t_enc < t_fric, we set the accretion radius to $$r_{\text{acc}}=R\sqrt{1+{\left(\frac{v_{\text{esc}}}{\Delta v}\right)}^2},$$(B.6)

which is the accretion radius for gravitational focusing. The parameter $v_{\text{esc}}=\sqrt{2Gm/R}$ is the escape velocity for a body of radius R (Ormel & Klahr 2010; Ormel 2017).

With the accretion radius at hand, we can determine the mass accretion rate. If the accretion radius is larger than the scale height of the pebble layer, H_peb, accretion takes place in the 2D regime and the mass accretion rate becomes $${\dot{m}}_{2\text{D}}=2\pi r_{\text{acc}}\Delta v{\Sigma }_{\text{peb}}.$$(B.7)

If, on the other hand, r_acc < H_peb, accretion is in 3D and the mass accretion rate reads $${\dot{m}}_{\text{3D}}=\frac{{\Sigma }_{\text{peb}}}{\sqrt{2\pi }{H}_{\text{peb}}}\pi r_{\text{acc}}^2\Delta v$$(B.8)

(Lambrechts & Johansen 2012; Johansen & Lambrechts 2017; Ormel 2017; Lambrechts et al. 2019). Typically, planetesimals start in the 3D regime and make a transition into 2D when they become massive enough. Accretion of pebbles stops when the body reaches its pebble isolation mass. In that case the body opens a gap in the gas which creates a pressure bump and prevents pebbles from drifting in. We take the pebble isolation mass as $$m_{\text{iso}}\approx 20{\left(\frac{H_{\text{g}}/r}{0.05}\right)}^3{M}_{\oplus },$$(B.9)

where H∕r is the aspect ratio of the gas disc (Johansen & Lambrechts 2017). Accretion also stops when the inclined planetesimal reaches a height that is above the pebble scale-height. Lastly, we take mutual filtering into account. That means that the pebble flux decreases towards inner orbits because every planetesimal on an outer orbit takes out a certain portion of the pebble flux (Liu et al. 2019).

Fig. B.1 Accretion radius and mass accretion rate for pebble accretion at 5 au. Shown are racc (top panel) and ṁ (bottom panel) as a function of planetesimal mass. We compare our prescription (black solid line) to the models developed in Johansen & Lambrechts (2017; JL17, blue dashed line) and Lambrechts et al. (2019; L19, dashed orange line). The vertical dashed line indicates the initial mass of the planetesimals used in Sect. 3.8. The horizontal dashed line (top panel) is the pebble scale-height, which indicates the transition from 3D to 2D accretion.

Fig. B.2 Accretion radius and mass accretion rate for pebble accretion at 10 au. Shown are racc (top panel) and ṁ (bottom panel) as a function of planetesimal mass. We compare our prescription (black solid line) to the models developed in Johansen & Lambrechts (2017; JL17, blue dashed line) and Lambrechts et al. (2019; L19, dashed orange line). The vertical dashed line indicates the initial mass of the planetesimals used in Sect. 3.8. The horizontal dashed line (top panel) is the pebble scale-height, which indicates the transition from 3D to 2D accretion.

Figures B.1 and B.2 show the accretion radius and mass accretion rate for pebble accretion at 5 au and 10 au, respectively.They represent the inner and outer edge of the ring of planetesimals simulated in Sect. 3.8. We used the same parameters for the Stokes number, the pebble flux, and the ratio of pebble to gas scale-height as in our simulation to generate these plots. Our prescription is in general very similar to that of Johansen & Lambrechts (2017) and Lambrechts et al. (2019). Our prescription has some differences because we do not include the weak-coupling Bondi regime (Ormel & Klahr 2010; Lambrechts et al. 2019) and our mass accretion rates are slightly higher than in Lambrechts et al. (2019). The vertical dashed lines in Figs. B.1 and B.2 indicate the initial mass of the planetesimals (10⁻⁴ M_⊕) used in our simulation. The biggest differences compared to the other two models occur around that mass for heliocentric distances between 5 au and 10 au. We therefore overestimate the growth of the low-mass planetesimals, but capture the embryos fairly well. The transition from 3D to 2D accretion is indicated by the horizontal dashed line, which shows the pebble scale-height. Only for the most massive bodies, ≳5 M_⊕, is accretion in the 2D regime. The prescription used here is thus sufficient for testing purposes without producing entirely incorrect results.

Appendix C Gas drag and migration

C.1 Gas drag

We implement gas drag through the protoplanetary disc acting on the planetesimals. We use as an example the simple minimum mass solar nebula (MMSN) model (Weidenschilling 1977b; Hayashi 1981). The surface density and temperature of the gas decrease with increasing heliocentric distance as power laws with respective slopes α and β: $$\begin{array}{ll}\Sigma \hfill & =1700\text{g}{\text{cm}}^{-2}{\left(\frac{r}{\text{au}}\right)}^{-\alpha },\hfill \\ T\hfill & =280\text{K}{\left(\frac{r}{\text{au}}\right)}^{-\beta }.\hfill \end{array}$$

The drag force acting on the planetesimal is determined by a drag coefficient C_D, the gas density ρ, the aerodynamic cross-section of the body of radius R, and the relative velocity v_rel = v −v_g between the planetesimal and the gas: $$F_{\text{D}}=-\frac{1}{2}{C}_{\text{D}}\rho \pi R^2\left|v_{\text{rel}}\right|v_{\text{rel}}.$$(C.3)

We assume the disc to be in vertical hydrostatic equilibrium. Therefore, the density is $$\rho =\frac{\Sigma }{\sqrt{2\pi }{H}_{\text{g}}}\exp \left(-\frac{z^2}{2H_{\text{g}}^2}\right),$$(C.4)

where H_g = c_s∕Ω_K is the gas scale-height of the disc, which is the ratio of sound speed c_s and Keplerian angular frequency Ω_K. The gas velocity is sub-Keplerian and can be written as $$v_{\text{g}}=v_{\text{K}}(1-\eta ),$$(C.5)

where η is set by the pressure gradient $$\eta =-\frac{1}{2}{\left(\frac{c_{\text{s}}}{v_{\text{K}}}\right)}^2\frac{\partial \log P}{\partial \log r}.$$(C.6)

The drag coefficient C_D depends on the properties of the planetesimal and of the disc gas. If the size of the body is smaller than the mean-free path of the gas λ_g (i.e. R∕λ_g < 4∕9) the body is in the Epstein drag regime and C_d becomes (Epstein 1924) $$C_{\text{d}}=\frac{8}{3}\frac{v_{\text{thm}}}{\left|v_{\text{rel}}\right|},$$(C.7)

where $v_{\text{thm}}=\sqrt{8/\pi }{c}_{\text{s}}$ is the thermal velocity of the gas. On the other hand, if the body is larger, the drag coefficient is more complicated because the drag arises from the hydrodynamic flow across the object. Depending on the Reynolds number Re = 2R|v_rel|∕ν, where ν = (1∕3)λ_gv_thm is the molecular viscosity, the drag coefficient becomes $$C_{\text{d}}=\{\begin{array}{ll}24R{e}^{-1}\hfill & Re<1\hfill \\ 24R{e}^{-0.6}\hfill & 1\le Re<800\hfill \\ 0.44\hfill & Re\ge 800\hfill \end{array}$$(C.8)

(Whipple 1972; Weidenschilling 1977a). If the body moves supersonically through the gas, that is |v_rel | > c_s, the drag coefficient takes on a constant value of C_d = 2 (Brasser et al. 2007).

C.2 Type I migration

In addition to gas drag, we also implement type I migration and the damping of eccentricity and inclination by adding a force that mimics the results of hydrodynamic simulations (Cresswell & Nelson 2008). The fundamental damping timescale is $$t_{\text{wave}}=\frac{M_{\odot }}{M}\frac{M_{\odot }}{\Sigma r^2}{\left(\frac{H_{\text{g}}}{r}\right)}^4{\Omega }_{\text{K}}^{-1} ,$$(C.9)

which is the timescale for the damping of semi-major axis a, eccentricity e, and inclination i due to the excitation of density waves in the disc by the planetesimal (Tanaka & Ward 2004).

Cresswell & Nelson (2008) then provide timescales for the damping of a, e, and i that fit their hydrodynamical simulations most accurately. These are $$\begin{array}{ll}{t}_e\hfill & =\frac{t_{\text{wave}}}{0.780}\left[1-0.14{\tilde{e}}^2+0.06{\tilde{e}}^3+0.18\tilde{e}{\tilde{i}}^2\right],\hfill \\ t_i\hfill & =\frac{t_{\text{wave}}}{0.544}\left[1-0.30{\tilde{i}}^2+0.24{\tilde{i}}^3+0.14{\tilde{e}}^2\tilde{i}\right],\hfill \\ t_a\hfill & =\frac{2t_{\text{wave}}}{2.7+1.1\alpha }{\left(\frac{H_{\text{g}}}{r}\right)}^{-2}\{P(\tilde{e})+\frac{P(\tilde{e})}{|P(\tilde{e})|}\hfill \\ \hfill & \times \left[0.070\tilde{i}+0.085{\tilde{i}}^4-0.080\tilde{e}{\tilde{i}}^2\right]\},\hfill \end{array}$$

with $$P(\tilde{e})=\frac{1+{(\tilde{e}/2.25)}^{1.2}+{(\tilde{e}/2.84)}^6}{1-{(\tilde{e}/2.02)}^4}.$$(C.13)

Here α is the surface density slope of the gas, and eccentricity and inclination are divided by the aspect ratio of the disc as ẽ = e∕(H_g∕r) and ĩ = i∕(H_g∕r). The accelerations that act on the planetesimal are then $$\begin{array}{ll}{a}_a\hfill & =-\frac{v}{t_a},\hfill \\ a_e\hfill & =-2\frac{(v\cdot r)r}{r^2{t}_e},\hfill \\ a_i\hfill & =-\frac{v_z}{t_i}{e}_z,\hfill \end{array}$$

where e_z is the unit vector in z-direction.

References

All Figures

Fig. 1 Viscous stirring for equal-mass planetesimals. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of equal-mass planetesimals with mass 1 × 1024 g. The ring is centred at 1 au and the surface density is 10 g cm−2. Shown are the results of our model: the close-encounter approach (orange), a full N-body simulation with mercury (grey), and the semi-analytic toy model of Ohtsuki et al. (2002; black dashed).

In the text

Fig. 2 Viscous stirring for bimodal mass distribution. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of planetesimals with bimodal mass distribution. There are 800 planetesimals with mass 1024 g (component 1, m1) and 200 with mass 4 × 1024 g (component 2, m2). The ring is centred at 1 au and the surface density is 10 g cm−2. Shown are the results of our model (orange), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black dashed). Component 1 and 2 are indicated with solid and dotted lines, respectively.

In the text

Fig. 3 Viscous stirring and dynamical friction for a planet embedded in a ring of planetesimals. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) of a ring of planetesimals with an embedded planet. There are 1000 planetesimals with mass 1024 g and 1 planet with mass 1026 g. The ring is centred at 1 au and the surface density is 10 g cm−2. The planet is initially at 1 au. Shown are the results of our model (orange), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed). Planetesimals and the planet are indicated with solid and dotted lines, respectively.

In the text

Fig. 4 Viscous stirring of equal-mass planetesimals for increasing close-encounter distance using the sequential approach. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) for the same setupas for the equal-mass planetesimal run. Shown are the results of our model for different values of dce (coloured lines), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed).

In the text

Fig. 5 Viscous stirring of equal-mass planetesimals for increasing close-encounter distance using the group approach. The figure shows root mean square eccentricity (top panel) and inclination (bottom panel) for the same setupas for the equal-mass planetesimal run. Shown are the results of our model for different values of dce (coloured lines), a full N-body simulation with mercury (grey), and the semi-analytic model of Ohtsuki et al. (2002; black-dashed).

In the text

Fig. 6 Scaling of computational time with number of bodies. The setup is the same as in Sect. 3.1. The surface density is kept constant at 10 g cm−2. The number of planetesimals increases from 1000 to 10 000. Shown are the total computational time to integrate 104 yr with thesequential approach (sequential; orange), the group approach (grouping; blue), and, for comparison, with mercury (grey diamond). Also shown is a power-law fitting to find the scaling with N for our model (blue and orange lines).

In the text

	Fig. 7 Energy and angular momentum. The figure shows the time evolution of the relative change of energy (top panel) and angular momentum (bottom panel) for equal-mass planetesimals. Shown are the results obtained with the sequential approach (sequential; orange) and the group approach (grouping; blue) of our method and compared with mercury (grey).
In the text

Fig. 8 Collisional growth of planetesimals at 1 au. Shown are scatter plots for eccentricity (first row), inclination (second row), and mass (third row). 4000 planetesimals are located at 1 au in a ring of width 0.085 au. The three columns show snapshots at different times. Planetesimals that grow more massive than 10−3 M⊕ are shown as coloured filled circles. All other bodies are grey open circles. The colour and size of the symbols scale with the instantaneous mass.

In the text

	Fig. 9 Mean and maximum mass of planetesimals. Shown are the mean (orange line) and maximum (blue line) mass of the planetesimals as a function of time. Both lines are averages of five runs with different random realisations of the initial conditions.
In the text

Fig. 10 Pebble accretion and collisions in the giant planet region. Shown are the scatter plots for eccentricity (first row), inclination (second row), and mass (third row). A total mass of 0.1 M⊕ is distributed amongst 450 planetesimals between 5 au and 10 au. The three columns show snapshots at 0 Myr, 1.5 Myr, and at the final time of 3 Myr. Bodies that are more massive than 10−2 M⊕ are shown as coloured filled circles, their growth tracks are shown by the coloured lines, and the centre of the circle is indicated by a black dot. All other bodies are grey open circles. The colour and size of the symbols scale with the instantaneous mass.

In the text

	Fig. 11 Cumulative number of bodies for pebble accretion in the giant planet region. Shown are N(> m), the cumulative number of bodies with mass higher than m, for the initial distribution and for the distribution at times 0, Myr, 1 Myr, 2 Myr, and 3 Myr. Each line is the average of five runs with randomised initial configuration of the bodies.
In the text

Fig. 12 Mass evolution of planetesimals for collisions and pebble accretion in the giant planet region. Each line shows the growth history of a planetesimal or embryo. Bodies that grow more massive than 10−2 M⊕ are shown as coloured lines. The colour scales with mass and is the same as in Fig. 10. The grey lines show all other bodies. The endpoint of the simulation for each body is indicated with a circle. Because the bodies do not grow significantly before 104 yr, only the later times are shown.

In the text

Fig. 13 Ratio of mean-free path to close-encounter distance as a function of heliocentric distance. The solid lines represent λ∕dce for different values of dce. For the dashed line, this ratio is unity, indicating where the sequential approach breaks down. A value of Σ0 = 10 g cm−2 at 1 au and irms = 5 × 10−5 was used. Surface density is assumed to be a power law $\Sigma ={\Sigma }_0{(r/\text{au})}^{-1}$.

In the text

	Algorithm 1 Construct cell list
In the text

	Algorithm 2 Detect close encounter
In the text

	Algorithm 3 Grouping of close encounters
In the text

Fig. B.1 Accretion radius and mass accretion rate for pebble accretion at 5 au. Shown are racc (top panel) and ṁ (bottom panel) as a function of planetesimal mass. We compare our prescription (black solid line) to the models developed in Johansen & Lambrechts (2017; JL17, blue dashed line) and Lambrechts et al. (2019; L19, dashed orange line). The vertical dashed line indicates the initial mass of the planetesimals used in Sect. 3.8. The horizontal dashed line (top panel) is the pebble scale-height, which indicates the transition from 3D to 2D accretion.

In the text

Fig. B.2 Accretion radius and mass accretion rate for pebble accretion at 10 au. Shown are racc (top panel) and ṁ (bottom panel) as a function of planetesimal mass. We compare our prescription (black solid line) to the models developed in Johansen & Lambrechts (2017; JL17, blue dashed line) and Lambrechts et al. (2019; L19, dashed orange line). The vertical dashed line indicates the initial mass of the planetesimals used in Sect. 3.8. The horizontal dashed line (top panel) is the pebble scale-height, which indicates the transition from 3D to 2D accretion.

In the text