© ESO, 2013

1. Introduction

Debris discs are circumstellar optically thin dusty discs around main-sequence stars with little or no gas¹. Hundreds of such discs have been detected through the IR-excess associated with the warm and/or cold circumstellar dust in the micron- to centimetre-size range (see reviews by Wyatt 2008 and Krivov 2010).

The observed dust cannot be primordial, because its removal timescale due to Poynting-Robertson drag and destructive collisions is much shorter than the system’s age. It must thus be steadily replenished from a mass reservoir contained in larger, unseen bodies, most probably through a collisional cascade starting at bodies exceeding the kilometre-size range. Studying the collisional evolution of such systems is thus of fundamental importance, as it allows one to address crucial questions such as the size-distribution of particles in the disc, the link between the observed dust and its hidden larger progenitors, the total mass in the system, and its mass loss with time. All these questions are usually numerically investigated with particle-in-a-box codes, with no or poor spatial resolution, dividing the disc into size bins whose mutual collisional interactions are treated in a statistical way (Kenyon & Bromley 2002; Thebault et al. 2003; Krivov et al. 2006; Thebault & Augereau 2007; Löhne et al. 2008; Gáspár et al. 2012). Such models have greatly improved our understanding of debris discs in the past decade and are essential for understanding the formation and evolution of planetary systems.

For a minority of debris discs (40 to date²), resolved images have also been obtained. These images are usually obtained at visible and near-IR wavelengths, but new observational facilities such as the Herschel Space Observatory are now also providing increasingly more images in the mid- to far-IR (Löhne et al. 2012; Wyatt et al. 2012; Lestrade et al. 2012; Donaldson et al. 2012; Ertel et al. 2012). Almost all of these resolved images show pronounced structures, such as bright clumps, two-side asymmetries, spiral arms, or warps (e.g. Kalas et al. 2005; Golimovski et al. 2006; Acke et al. 2012; Booth et al. 2013). Understanding the origin of these structures has been a major objective of debris disc studies, which have explored several possible scenarios mostly through numerical modelling. In many cases, these scenarios involve a perturbing planet, for instance for the warp of the β Pictoris disc (Mouillet et al. 1997; Augereau et al. 2001), the brightness asymmetries in the Epsilon Eridani system (Kuchner & Holman 2003), or the confined rings around HR4796A and Fomalhaut (Wyatt et al. 1999; Chiang et al. 2009). However, alternative scenarios, such as stellar perturbations (e.g., Augereau & Papaloizou 2004; Thebault et al. 2010), interaction with gas (Takeuchi & Artymowicz 2001; Besla & Wu 2007), interaction with the interstellar medium (Artymowicz & Clampin 1997; Debes et al. 2009; Marzari & Thebault 2011), or transient violent events (Kenyon & Bromley 2005; Grigorieva et al. 2007) have also been investigated.

The numerical exploration of these different scenarios has mostly been made using N-body codes, which are designed to accurately follow the development of dynamical structures such as resonances and migrations (e.g., Reche et al. 2009; Kuchner & Holman 2003). The price to pay for the spatial precision allowed by the N-body approach is that collisions are usually neglected in these models, which follow the evolution of collisionless test particles.

As can be seen, the collisional and dynamical evolution of debris discs are usually studied separately, the focus of collisional models being mostly global characteristics such as grain-size distributions or mass loss, while dynamical models focus on local or transient phenomenae that can leave signatures in resolved images. While such separate studies can (and have) produce(d) important results, they suffer from unavoidable limitations.

For the purely statistical models these limitations are obviously the absence of, or only poor, spatial resolution, but also the likelihood that dynamical processes might strongly affect impact rates and velocities and hence the collisional evolution. This might lead to strong discrepancies in the way some parts of the disc collisionally evolve with respect to others.

The N-body collisionless models also suffer from severe handicaps. The absence of collisions can indeed strongly bias or even invalidate results obtained in such codes, and this for several reasons. As an example, if collisional timescales are shorter than dynamical ones, collisions can hinder or even prevent the build-up of dynamical structures. Likewise, the identification of dynamically stable and unstable regions in a perturbed debris disc (be it by a planet or a stellar companion) can also be strongly affected by collisional activity, as a collision cascade will steadily produce small grains that can be launched by radiation pressure on highly eccentric or unbound orbits that can populate regions that are in principle dynamically forbidden.

Unfortunately, the complex interplay between the system’s dynamical and collisional evolutions is very difficult to handle numerically, mainly because the particle-in-a-box and the N-body approaches are radically different in their principles and structures. However, the first attempts to partially couple dynamics and collisions have been published recently, most of them taking the N-body approach as a basis into which some collision-imposed properties are injected.

The most reliable way of including destructive collisions into an N-body scheme would in principle be by brute-force methods, where bodies are effectively broken into fragments whose evolution is then dynamically followed (Beauge & Aarseth 1990). However, such codes can only follow a very limited number of collisional fragments and lead in any case to an exponential increase of the number of particles that very quickly becomes unmanageable.

An alternative option is to run collisionless N-body runs, and to post-process them assuming that each test particle stands for a dust-producing collisional cascade of solids (Booth et al. 2009). This approach allows one to explore very long timescales, up to several Gyr. It has been used to study the signature, in terms of spectral energy distribution (SED), of the late stages of terrestrial planet formation (Raymond et al. 2011). However, it has some important limitations, because it implicitly supposes that collisions and dynamics are fully decoupled, that is, it assumes that collision processes are unaffected by the dynamics, but also that collisions have no influence on the formation and fate of spatial structures.

An important recent improvement is the LIPAD code developed by Levison et al. (2012). LIPAD follows the dynamical evolution of tracers, which represent populations of single-sized planetesimals, and tracks down their mutual collisions. Collisions change tracer velocities as well as increase or decrease the planetesimal sizes they represent, depending on the collision outcome (accretion or fragmentation). However, this code is designed to study the earlier stage of kilometre-sized planetesimal accretion and is not adapted, in its present form, to debris disc studies. Indeed it mainly focuses on the fate of large growing bodies and has a very simplified prescription for the dust: all material below the smallest planetesimal size is converted into dust tracers that have a single physical size and no longer collisionally interact with each other. Furthermore, it neglects the effect of radiation pressure, an effect that is crucial for debris discs because it creates a strongly size-dependent behaviour of dust particles, which makes the numerical treatment of the dust population much more complex.

Amongst the most sophisticated debris disc models that have been developed to date are the CGA (Stark & Kuchner 2009) and DyCoSS (Thebault 2012) algorithms. These codes are different in their principles, but both implement similar levels of collisional processes into their N-body main structure. They are specifically designed to estimate how collisional lifetimes of grains are affected by dynamical perturbations, and how this variety of collisional lifetimes in turn affects the development of dynamical structures (resonances, PR-drag migrations, etc.). Both codes allow studying very fine spatial structures and have been used for some specific cases such as the Neptune-Kuiper-Belt system (Kuchner & Stark 2010), debris discs in binaries (Thebault 2012) or discs with an embedded or exterior planet (Thebault et al. 2012). However, they are both restricted to specific set-ups, that is, systems at collisional and dynamical steady-state under the influence of one perturbing body only. Moreover, the coupling between dynamics and collisions is only partial, and both codes assume that collisions are fully destructive, that is, particles that have exceeded their collisional lifetimes are simply removed, and the fate of small collisional fragments is not followed.

The first attempt at truly coupling the collisional and dynamical modelling of a debris disc was performed with the hybrid code of Grigorieva et al. (2007). The principle of this approach is that each particle of an N-body simulation is a tracer (or a super-particle) that represents a vast population of particles of a given size, and that the code tracks mutual collisions between these tracers. These collisions are then treated with a classical particle-in-a-box approach, estimating the amount of mass lost by the impacting tracers and the corresponding mass that is injected into new tracers that carry the collisional fragments created. This code was successfully used to identify and quantify the mechanism of collisional avalanches produced by the shattering of a large planetesimal deep inside a dense dust belt. It had one major limitation, however, which is that the number of tracers is constantly increasing and rapidly becomes too much for the code to handle. As such it was restricted to short timescales, typically a few orbital periods.

This major limitation has recently been overcome by an algorithm aimed at studying the very different astrophysical case of young, massive, and gaseous protoplanetary discs (Charnoz & Taillifet 2012). This code, named LIDT3D, is based on a similar tracer approach but integrates a routine to identify superfluous and redundant tracers and to reassign these tracers to regions where they are needed more. This procedure prevents the number of tracers from increasing and allows studies spanning much longer times. However, LIDT3D cannot be used to study debris discs because of the very different physical processes at play in a protoplanetary disc, in particular, the strong homogenizing effect of gas drag that greatly simplifies the dynamics, as well as the very different collisional regime that prevails, that is, mostly low-velocity impacts instead of high-Δv shattering-and-fragment-producing ones in debris discs.

Our aim is here to create a new version of the LIDT code, called LIDT-DD, which keeps the basic tracer architecture of LIDT3D, but is able to handle the very constraining demands of modelling debris discs. We present in Sect. 2 a brief summary of the main LIDT architecture that will be kept in the new algorithm. In Sect. 3, we first describe the challenging aspects of debris disc physics and their implications in terms of numerical handling. We then present the LIDT-DD code itself as well as an illustrative pedagogical run. Sect. 4 is devoted to different tests that have been performed to determine the code’s reliability. In Sect. 5, we illustrate the potential of this new code by showing some preliminary results for an example set-up of a massive catastrophic collision within a debris disc.

2. LIDT basic architecture

The basic architecture of the LIDT-DD debris disc model takes its roots in the LIDT3D code (OPEN-MP) developed by Charnoz & Taillifet (2012) for the study of primordial protoplanetary discs. Because of the radical differences between these two astrophysical cases in terms of dominant physical mechanisms but also regarding the strong constraints imposed by high-velocity destructive collisions in debris discs (see Sect. 3.1), LIDT-DD strongly departs from its predecessor, incorporating several new key procedures, and can be considered as an independent, stand-alone code. However, before presenting it in detail in the next section, we first briefly recall here the basic features of the LIDT code that remains as the backbone of this new model (for a more detailed description, see Charnoz & Taillifet 2012).

2.1. Tracers

The building blocks of the LIDT scheme are tracers, each representing a whole population of particles of a given size at a given location in the system. At each time step, these tracers are dynamically evolved with a deterministic N-body integrator and do collisionally interact following a statistical particle-in-a-box scheme.

Tracers are defined by their position R_t, velocity V_t, the physical size s_t of the particles they stand for, and the individual mass ${\mathit{m}}_{\mathrm{t}}\mathrm{=}\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi \rho }{\mathit{s}}_{\mathrm{t}}^{\mathrm{3}}$ of these particles. A tracer represents a large number N_t of such individual particles, and thus represents a total mass M_t such as M_t = N_t   m_t. The whole system is then composed of N such tracers, a number that can evolve with time depending on the number of tracers that are needed in every region of the disc. For our LIDT-DD debris disc simulations, N is of the order of ~10⁵, which is the typical maximum total number of tracers that are needed in our simulations to handle two tracers per dynamical category per size bin per spatial cell per time step (see Sect. 3.7). N will be worked out automatically by the code at every time step because the number of tracers directly depends on the dynamical and collisional activity within the disc (see Sect. 3.8 for more details).

2.2. Dynamical evolution

The dynamical evolution of a tracer representing a cloud of particles of size s is followed using a Lagrangian N-body approach, integrating at each time step the equation of motion corresponding to one particle of size s at position R_t with velocity V_t. Forces that are taken into account are the central star’s gravity, stellar radiation pressure, and friction generated by turbulent gas. The integrator is a Burlisch-Stoer with a semi-implicit solver (Bader & Deuflhard 1983).

2.3. Collision treatment

Once all tracers have been dynamically evolved for a time step, their collisional evolution is computed in the following way: the system is divided into a 2D grid, in (r,z) or (r,θ), of spatial cells, into which the collisional evolution will be followed. Within each cell, mutual collisions between tracers of each size group are considered. The first step is to compute average impact velocities ΔV_i,j between all pairs of sizes s_i and s_j. Then, from the values of ΔV_i,j and the particles number N_i and N_j, the number of collisions Nc_i,j between all particles contained in tracer i and those contained in tracer j within a time step Δt is derived, following: $\mathit{N}{\mathrm{c}}_{\mathit{i,j}}\mathrm{~}\frac{\mathrm{\Delta }{\mathit{V}}_{\mathit{i,j}}\mathit{.}{\mathit{N}}_{\mathit{i}}\mathit{.}{\mathit{N}}_{\mathit{j}}\mathit{.}{\mathit{\sigma }}_{\mathit{i,j}}}{{\mathrm{Vol}}_{\mathit{i,j}}}\mathrm{\Delta }\mathit{t,}$(1)where Vol_i,j is the cell’s volume and σ_i,j = π(s_i + s_j)² is the total cross-section during the impact.

Then, for each impacting pair (i,j), a standard collisional outcome prescription is plugged in to estimate wether the collision results in accretion, rebound, or erosion. In the latter case, the number of collisional fragments produced for each size is estimated using simple collisional laws (see Charnoz & Taillifet 2012). For each fragment size, the number of produced fragments is then added to the tracers present in the cell corresponding to this size. If no corresponding tracers are present, a new one is created.

A key feature of the LIDT scheme is a tracer reassignment procedure that avoids an unmanageable increase of the number of tracers. Its basic principle is that each redundant tracer, for instance when there are too many tracers of a given size in a given cell than would be required to produce a statistically significant behaviour, is taken away (its mass being given to other neighbouring tracers of the same size) and stored in a buffer of freed tracers that can be used whenever new tracers are created by collisions in other cells.

To avoid too fast movements of the tracers in and out of the cells, which would put an artificial constraint on the time step, the (r,θ) grid is rotating at the Keplerian velocity calculated at the centre of each cell.

Note that the collisional spatial grid necessarily has a finite spatial extension and cannot extend over the whole space over which the dynamical evolution of tracers is followed. There is thus a region beyond the outer limit of the grid where collisions are implicitly not taken into account, but this is an acceptable simplification if the outer limit of the grid is located in regions that are sparsely populated and not very collisionally active.

In summary, the basic principle of the LIDT scheme can be schematically presented as follows:

• Step 1: Evolve dynamics.

• Step 2: Create a 2D grid superimposed onto the tracers to divide the system into different spatial cells.

• Step 3: Compute relative velocities between tracers in each cell.

• Step 4: Compute collisional outcomes and produce collisional fragments in each cell.

• Step 5: Create tracers and reorganize them within the system.

3. LIDT-DD

3.1. Specificities of debris discs physics

The main changes that have been implemented in LIDT-DD with respect to the initial protoplanetary-disc version were motivated by the specificities of debris discs physics, which impose very strong constraints on the way collisions and dynamics are to be treated. We present here a brief description of what these specificities are and how they affect the numerical treatment of the system within the LIDT frame.

One crucial difference between protoplanetary and debris discs is the absence (or strong depletion) of gas in the latter. While this absence is a simplification in terms of the physical processes at play, in particular allowing one to dispense with the problematic handling of turbulence and gas-disc parameterizing (see Charnoz & Taillifet 2012), it does in fact add great complexity to the system’s dynamical and collisional evolution. Indeed, in a protoplanetary disc, small grains are strongly coupled to the gas, and gas drag very quickly smoothes out any disparities between orbits of similarly sized particles. As a consequence, the way these particles are produced and the orbits into which they are released is of little importance because any initial conditions are very quickly relaxed. The welcome consequence is that at any given location in the disc, there is only one possible dynamics for a given particle size. In numerical terms, this means that in a given cell of the LIDT code, only one category of tracer is needed to represent one size bin. This is unfortunately no longer the case for debris discs, where the absence of gas drag means that initial conditions are never relaxed. So that, depending on where and how they are produced, particles of the same size can have very different dynamics even though they are in the same region of the system. This puts very strong constraints on the numerical treatment because several categories of tracers might be needed for each size bin in each cell. A collateral problem is that this number of categories cannot be known in advance. Another consequence is that the way new particles (or their tracers) are produced suddenly becomes a problem because this will control their fate and the fate of the new particles they will in turn spawn because of later collisions, whereas in LIDT3D tracers are simply automatically given the only possible gas-imposed dynamics in the region they are produced.

Another major difference to the protoplanetary-disc case is that impact velocities are much higher, partly because the damping effect of gas drag is no longer present, but also because debris discs are expected to be dynamically stirred by massive (and often unseen) bodies (see, e.g., Thebault 2009, for a discussion on this matter). These high velocities have dramatic consequences on the treatment of collisions, because they will lead to very destructive impacts, producing numerous small fragments of all sizes. In the protoplanetary-disc case, no such high-dv impacts are expected, so that the treatment of collisions was handled in a very simplified way, that is, by setting a threshold velocity of 1 m/s, regardless of the impacting bodies sizes and compositions, beyond which all impacts were considered as erosive, and all erosive impacts resulted in similar outcomes in terms of mass loss and fragment distribution. Such simplified laws cannot be used for debris discs, because the outcome of fragmenting impacts strongly depends on parameters such as colliding bodies’ mass ratio and velocities (e.g., Benz & Asphaug 1999). The consequences of, for instance, a 100 m s^-1 impact and a 1 km s^-1 one are radically different in terms of mass loss, size of the largest remaining fragment, or size distribution of the produced debris.

The last fundamental specificity of debris discs is the crucial role played by grains very close to the blowout size s_cutoff imposed by radiation pressure. Using the parameterization with the ratio β = F_PR/F_grav, which is ∝1/s in a wide size range, these grains correspond to values very close to β = 0.5 ³. As has been shown by most collisional evolution models, the total geometrical cross-section in debris discs, and hence their luminosity at all wavelengths short of the mid-to-far IR, is dominated by solids in the s_cutoff to ~2−3   s_cutoff range (e.g. Thebault & Augereau 2007). Unfortunately, in contrast to protoplanetary discs that are made to behave very similarly to grains of other sizes by gas coupling, in gas-poor discs these small grains have a very complex dynamical and collisional evolution. They are indeed placed by radiation pressure on high-eccentricity orbits, making them populate vast regions that extend very far from their production location. Conversely, small grains present at any given location in the disc can potentially originate from far away regions much closer to star. Furthermore, small grains have strongly varying orbital (and thus impact) velocities along their elongated orbits. Last but not least, close to the s ~ s_cutoff limit, this complex dynamical behaviour becomes extremely sensitive to very small size differences. Indeed, a grain with a high β produced from a circular orbit of semi-major axis a has an apastron $\mathit{Q}\mathrm{=}\mathit{a}\hspace{0.17em}\hspace{0.17em}\frac{\mathrm{1}\mathrm{-}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\beta }}\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\right)\mathrm{=}\frac{\mathit{a}}{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\beta }}\mathrm{·}$(2)This means that grains with β = 0.48 have an apastron of 25a, while β = 0.45 grains, which are just 7% bigger, reach only 10a. These extreme characteristics of the particles that happen to be the ones dominating disc luminosities put very strong constraints on the way the code has to handle the critical ~[s_cutoff,2s_cutoff] size range.

3.2. Dynamical evolution

The dynamical computation part of the code is very modular and could be adapted to the present case. One important update made to fit debris disc physics was to add Poynting-Robertson drag, following the equation (Robertson 1937) $\frac{{\mathrm{d}}^{\mathrm{2}}{r}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}^{\mathrm{2}}}\hspace{0.17em}{{e}}_{{r}}\mathrm{+}\frac{\mathit{\beta G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}^{\mathrm{2}}}\left[\left(\mathrm{1}\mathrm{-}\frac{{v}\mathrm{·}{{e}}_{{r}}}{\mathit{c}}\right){{e}}_{{r}}\mathrm{-}\frac{{v}}{\mathit{c}}\right]\mathit{,}$(3)where r is the position vector for a particle in the frame of its central star (of mass M_∗), r is its norm, v is the velocity vector, c is the speed of light, and e_r is the radial normalized vector used in spherical coordinates. Note, however, that PR-drag is never dominant in the highly collisional regime that is considered throughout this paper. Another important change was to take into account the possibility of handling dynamically perturbing bodies such as planets or stellar companions.

As in LIDT3D, the dynamical evolution of the tracers is followed beyond the outer limit of the collisional grid (see Sect. 2.3). However, to avoid computing the evolution of irrelevant tracers too far removed from the region of interest, an outer limit r_outDyn for following their orbits was also imposed, beyond which tracers were considered as lost and were removed from the system. In contrast to the protoplanetary-disc case, this external border has to be located relatively far out to follow the orbits of the population of small, high-β particles that can have their apastron very far away from the central star. For the rest of the paper the dynamical outer border was fixed to r = 300 AU so that the smallest bound particles taken into account in the simulations can return and collide into their birth ring.

3.3. Particle-size sampling

As has been seen in Sect. 3.1, small high-β grains are of extreme importance because they dominate the disc luminosity for most observations. Furthermore, we have seen that these grains’ behaviour is very sensitive to small size differences. As a consequence, the resolution between adjacent size bins must be small enough to have a correct sample of grain sizes in this high-β domain.

Another major constraint on the size distribution is that it should extend to bodies that are large enough to sustain collisional cascades on long timescales, that is, bodies of size s_max should have a collisional lifetime exceeding the simulation timescale. This constraint will vary depending on the disc’s stirring, particle composition, and of course the considered timescales. To be on the safe side, we assumed s_max = 50 km, which corresponds to the smallest primordial bodies found by Löhne et al. (2008) for their highly stirred simulation after 5 × 10⁹ years.

With a classical size sampling where size bins are separated by a constant increment in logarithmic scale, these two constraints put together would lead to a number N_b of size bins that is much too large to be handled in the simulations, for which this number can typically not exceed ~60. However, these constraints can be relaxed, because a very fine size-sampling is only needed in the small-size range, but not for larger bodies, for which dynamical behaviours no longer have extreme variations with size. For these reasons we considered two different scales, a fine sampling for high-β grains, and a coarser one for all other sizes. For the fine-sampling domain, the logarithmic increment of size ϵ_f was taken as a free parameter in the [1.05,1.25] range. This choice of ϵ_f, combined with the value for N_b, automatically constrained both the limiting size s_lim between the fine and coarse domains, as well as the size increment ϵ_c in the latter, which is always within the [1.8,2.1] range. Note that the coarse size increment does not reach too high values, so that the code does not loose too much precision in the size-dependent treatment of collisions. As a matter of fact, it is similar to the size increment considered for the ACE code (Krivov et al. 2006). For the rest of the paper we considered a nominal case with ϵ_f = 1.15 for the smallest size bins and ϵ_c = 2 for the coarse domain.

3.4. Impact velocity estimates and collision search

For protoplanetary discs, impact velocities Δv_i,j in a given cell were simply derived by comparing the time-and-space-averaged local excitation velocities (i.e., departure from the local circular Keplerian velocity) of all tracers of sizes s_i and s_j. This fast and linear procedure cannot work in the much more complex context of debris discs because of possible short timescale variations and the potential existence of different dynamical categories within the same size range. As a consequence, we performed a real-time estimate of all mutual Δv_i,j = ⁽(Vim_x − Vjm_x)² + (Vim_y − Vjm_y)² + (Vim_z − Vjm_z)^2)0.5 values, where Vim and Vjm are Vi and Vj velocity vectors corrected for biases due to finite cell sizes. These biases are due to the Keplerian shear between tracers at different radial distances r and the orientation of velocity vectors because of different azimuthal positions θ within the cell, which has to be corrected for virtually rotating the tracers to the same longitude at the centre of the cell. In practice, we switched to spherical coordinates to facilitate the correction process. We first estimated the departures δV_i and δV_j from the local Keplerian velocity for each tracer. We then corrected the magnitude of these two vectors by a factor accounting for the Keplerian shear to obtain the debiased vectors Vim and Vjm. This process corrects both for the difference in radial position within the cell and the azimuthal one.

3.5. Collision outcome prescription

The simplified 1 m s^-1 velocity barrier between erosion and accretion used for the protoplanetary-disc case is clearly not adapted to the high-velocity destructive collisions regime of debris discs. We therefore implemented a collision outcome prescription as sophisticated as those commonly used in collisional evolution particle-in-a-box codes. We chose the energy scaling prescription described in Thebault & Augereau (2007), where erosive impacts are separated into two regimes, catastrophic or cratering, depending on the value of the specific impact energy per target mass unit Q_kin = E_col/M_t, where E_col = 1/2M_pM_tΔv²/(M_p + M_t) (kinetic energy of the impact in the barycentric frame) as compared to the critical specific energy Q^∗: ${\mathit{Q}}^{\mathrm{\ast }}\mathrm{=}{\mathit{\alpha }}_{\mathrm{1}}\hspace{0.17em}{\left(\frac{{\mathit{R}}_{\mathrm{tar}}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{\mathit{a}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{2}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\left(\frac{{\mathit{R}}_{\mathrm{tar}}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{\mathit{b}}\mathit{,}$(4)where R₀ = 1  m, M_p and M_t are the respective mass of the projectile and target, and Δv is the relative velocity between the two colliding bodies. The values for a,b,α₁ and α₂ can be found in the literature and depend on the physical composition of grains (e.g. Housen & Holsapple 1990; Benz & Asphaug 1999). The values used for this study were taken from Benz & Asphaug (1999): a = −0.38, b = 1.36, α₁ = 3.5 × 10³ J/kg and α₂ = 3 × 10^-8 (SI).

The first term in Eq. (4) represents the strength regime for which grains become harder to break as they become smaller (a < 0). The second term, dominant for larger bodies, is the gravitational regime, for which bodies become harder to break as they become bigger (b > 0). The code is able to cope with different chemical compositions of grains, but for the sake of clarity we assumed pure silicates throughout the rest of this paper.

The fragmenting (Q > Q^∗) and cratering (Q ≤ Q^∗) regimes were handled as described in the appendix of Thebault & Augereau (2007). For cratering impacts we in particular distinguished between large-scale craters and grain-hitting-a-wall small impacts, connecting these two regimes by a smooth polynomial transition. The only simplification that we implemented with respect to Thebault & Augereau (2007) is that we only considered one power law (instead of two) for the size distribution of fragmented or cratered fragments. The index p of the power law is constrained by the mass of the largest remaining fragment M_LF coupled to the mass-conservation condition and the constraint that there is only one body of mass larger than M_LF⁴. This gives $\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{-}\mathrm{1}\mathrm{-}\frac{{\mathit{N}}_{\max }\hspace{0.17em}{\mathit{M}}_{\mathrm{LF}}}{{\mathit{M}}_{\mathrm{tot}}}}\mathit{,}$(5)where N_max is the number of particles in the size bin of the largest fragment and M_tot is the total mass of fragments. M_LF is given by the prescription of Fujiwara et al. (1977) for fragmenting impacts and by that of Wetherill & Stewart (1993) for cratering ones.

3.6. Creating collision-fragment virtual tracers

As already mentioned, we tracked all mutual impacts between tracers within each cell. For each tracer-tracer collision, we followed the procedure described in the previous section and estimated the total collisional mass lost by the impacting tracers. We then distributed this produced mass of collisional fragments into every size bin spanned by the m < M_LF range, following the size distribution of index p. For every size bins that have received collisional matter, a new virtual tracer was created, carrying the amount of mass received. This tracer was produced at the position of the largest (that is, the one representing particles with the largest size) of the two impacting tracers and was given its velocity vector.

3.7. Dynamical category-sorting

At the end of the collision treatment routine, each spatial cell was thus populated with a vast number of virtual tracers, carrying the mass of the fragments created by all tracer-tracer impacts in addition to the remaining initial tracers that had lost some of their mass. The next and fundamental step of our procedure was to reduce this total number of tracers and regroup them into similar families for which only a few representative tracers were kept.

Identifying these families is a challenge in the context of the complex dynamics of debris discs, because LIDT-DD must be able to sort and regroup tracers not only as a function of physical size and spatial location (as in LIDT3D), but also according to their orbits, which can be very different for a given grain size at a given location (see Sect. 3.1). Note that for small, high-β grains this variety of dynamical behaviours for a given grain size is a problem even in non-perturbed quiet discs. Indeed, at any given radial distance r, such grains can either have been produced locally, and hence have their periastron q close to r, or much closer to the star, with q ≪ r and hence a very different orbit, and reach r because of the highly eccentric orbit imposed by radiation-pressure.

The sorting method has to be generic enough so that no special treatment will be required when applied to different cases (unperturbed discs, planet/disc interactions, transient massive break-ups, etc.). We chose the hierarchical cluster analysis of Ward (1963), which satisfies this requirement and has the advantage of not requiring the user to define in advance the number of different dynamical populations. This procedure searches for distinct groups, as defined by the distribution of mutual distances between tracers in a multi-parameter phase space. The two parameters that we considered for the cluster analysis are the tracers’ periastron q and apastron Q, which are enough to constrain a particle orbit present in a (r,θ) spatial cell. In practice, the dynamical family identification is performed in a q + a vs. Q − a plane, instead of simply q vs. Q one, to give the same weight to q and Q when searching for distances in the parameter phase space. Note that the longitude of the periastron ω is not needed for the sorting, because the limited azimuthal extension of the spatial cell implicitly confines the values for ω within a narrow range once the q and Q sorting has been concluded. Likewise, the particles’ angular position on their orbit is almost fully implicitly constrained by the narrow radial extension of the cell, the only two possible distinct configurations being particles on their way out (i.e., moving towards apastron) or in (moving towards periastron). To distinguish between these two categories, we simply sorted tracers into two subcategories depending on the sign of their radial velocity v_r. The only free parameter for the sorting was the precision criterion, which determines how different two populations need to be to be considered as dynamically distinct. To adjust this criterion we chose the trial-and-error approach by performing test runs with an unperturbed disc and chose the largest possible value that could reliably sort high-β particles as a function of their different production annuli. Note that this procedure bears some similarities with the closest-streamline-search of Stark & Kuchner (2009), except for the crucial difference that our procedure is able to create new, not pre-existing dynamical families as the system evolves.

The procedure is summed up in Fig. 1. This plot shows for a given spatial cell and a given size bin (corresponding to small grains with β = 0.44), the tracers’ orbital distribution in a q + a vs. Q − a plane, before (init tracers) and after (final tracers) a collisional time step. We also plot the virtual tracers that are temporarily created as a result of collisions amongst tracers corresponding to larger size bins. As can be seen, some of these virtual tracers will be kept as final tracers by the sorting procedure, while some will be discarded and their mass transferred to the selected final tracers from their dynamical family (see Sect. 3.8). To avoid loosing any dynamical information, tracers that were present at the beginning of the time step (init tracers) have priority over virtual tracers to be selected as final tracers. The set-up we considered corresponds to the pedagogical run presented in Sect. 3.10, with an unperturbed disc initially confined to a narrow ring around 12 UA. The spatial cell that is considered in Fig. 1 is located ~8  AU beyond the ring and is thus mostly populated with small grains whose orbits stretch to these regions because of their high-e imposed by radiation pressure. These grains, which originate from the inner ring, constitute one of the dynamical families (delimited by the thick black lines) identified by the sorting procedure, that is, the one in the lower-left region of the plot (smaller q and Q). The second dynamical family, with larger periastron, corresponds to grains that have been produced locally as fragments from collisions involving slightly larger particles.

3.8. Tracer reassignment

As soon as all dynamical categories were identified amongst all initial and virtual tracers present in a cell, we eliminated all redundant tracers within these families, that is, we kept only a small number of representative tracers that carried all the mass of their category.

Given the potentially large number of dynamical categories for each particle size at each given location, and to avoid an unmanageable number of tracers, we only keep a maximum number of two tracers per dynamical family per size bin per spatial cell. To select these two representative tracers, we gave priority to the initial tracers, that is, those that were already present at the beginning of the time step before the collisional stage, and either selected two tracers at random amongst them or, if there was only one such tracer, selected this tracer and chose the other at random amongst the virtual tracer pool. When there was no initial tracer at all within a given family, the two representative tracers were directly chosen at random amongst the virtual population.

After selecting the two representative tracers, all other tracers within the same family were discarded and their mass was evenly added to each of the two final tracers. All these discarded tracers were turned off and stored to be used at another time step in other cells when they were needed. As illustrated by Fig. 2, this sorting and selection procedure effectively limits the increase of the number of active tracers after an unavoidable initial adjustment period.

3.9. Collisional energy dissipation

High-velocity collisions necessarily alter the orbits of impacting bodies and their resulting fragments by redistributing and dissipating kinetic energy. In simplified collision models, where bodies are treated as inelastically bouncing hard spheres, the kinetic energy dissipation is usually modelled by a normal and a tangential restitution coefficient, ϵ_N and ϵ_T, for the relative velocity, whose standard values are ϵ_N = −0.3 and ϵ_T = 1 (e.g., Thebault & Brahic 1998; Marzari & Scholl 2000; Charnoz et al. 2001; Lithwick & Chiang 2007). This simple prescription is impossible to directly implement here given the complexity of the tracer creation, sorting and reattribution procedures. However, it is possible to post-process tracer orbits to account for the average energy dissipation induced by collisions.

In practice, at the end of our collision procedure, after selecting all final tracers within a spatial cell, we ran an additional procedure that treated these tracers as potentially bouncing hard spheres. For each pair of tracers i and j, we first computed mutual collision probabilities in the same way as in Sect. 2.3. Using these probabilities, we then randomly selected tracer pairs that effectively had their velocity vectors modified by inelastic collisions during time step dt. For these chosen pairs of tracers we then assumed that an inelastic collision occured between two particles with the physical sizes that these tracers represent. This collision was treated in the centre-of-mass frame and each body’s velocity was modified following the ϵ_N and ϵ_T prescription in the same way as in Thebault & Brahic (1998) and Charnoz et al. (2001). This procedure is of course not an exact estimation of the post-collision evolution of each individual tracer present in a cell, since a given tracer will have its velocity vector modified by this procedure typically every ~t_coll/dt time steps. However, on timescales exceeding a few t_coll it does accurately model the average energy dissipation induced by collisional activity and reproduces all the expected behaviours predicted for simple cases (see Sect. 4.3)

3.10. Illustrative test run

To better illustrate the way the code functions, we present a simple pedagogical run, considering the collisional evolution of an unperturbed debris disc. The system’s evolution is presented in Fig. 3, showing snapshots of the tracer’s distribution at different times.

We started from a narrow ring, between a_min = 11 AU and a_max = 12 AU, where all the matter was initially confined (Fig. 3a). The (r,θ) collisional grid where the collisional evolution of solids is treated (see Sect. 3.5) extended out to 150 AU. This grid consisted of 10 radial and 20 azimuthal cells from 10 to 150 AU, and a log-scale was used for the radial spacing between the cells. The grid was rotated differentially so that each ring (at same r) rotated at its local Keplerian velocity. The outer edge beyond which the dynamical evolution is not followed was set at r_outDyn = 300 AU. We considered a β Pic-like A5V star and pure compact silicates for all solids in the disc. The size distribution of this population of solids was followed from 50 km to the radiation-pressure cut-off size at 2   μm. This size range was divided into size bins for which the logarithmic spacing is 1.15 for the critical domain of small grains in the 2   μm to 0.1mm range, and 2 for larger solids (see Sect. 3.3). The total initial mass of the disc was 7 × 10²⁴ kg, corresponding to an average optical depth ⟨ τ ⟩  ~ 10^-3 in the birth ring, and was distributed following a differential power law in s^-3.5 from s_max to s_min. All the main parameters for the considered setup are summarized in Table 1.

In the earliest stage of the disc evolution, we see the tracers corresponding to the smallest grains moving out from the birth ring because of radiation pressure that places them on high-e orbits (Fig. 3b). This outward movement of the initial tracers that leave the ring is compensated for the creation of new tracers in the ring due to impacts amongst larger particles. As a consequence, the total number of tracers initially increases (see Fig. 2). The increase stops when the small-grain tracers produced in the ring have had enough time to return to the ring after travelling through their elongated orbits. The number of tracers thus reaches a plateau after a few dynamical timescales of the smallest grains, that is, a few 10³ years. After this point, no significant change is visible in the spatial distribution of the tracers, since each collisional cell is populated with approximately two tracers per dynamical category per size bin (see Sect. 3.7). The importance of these dynamical categories and of the dynamical class-sorting procedure is clearly illustrated in Fig. 4, which shows the spatial distribution of tracers corresponding to a very small size bin (β = 0.4) as well as to the thirtieth size bin (β = 5 × 10^-3), the latter corresponding to larger particles that are not affected by radiation pressure. As can easily be seen, the big tracers are logically confined to the innermost cells of the grid, while the small tracers populate the whole system. In addition, while the number of big tracers is on average close to four per collisional cell (two for those moving towards apastron and two for those on their way back), the number of small-grain tracers is much larger. This is because, for these small grains, there are, for the same given spatial location, several dynamical behaviours that are possible depending on where the grains have been produced (see Sect. 3.1).

It is important to stress that the system has not yet reached a steady state by the time the number of tracers has reached its plateau (a few 10³  years). It will continue to evolve because of collisions that progressively change and redistribute the amount of mass carried by each tracer. Therefore even if the global distribution of tracers no longer changes after ~2 × 10³ years (Fig. 3c and d), the corresponding 2D maps of the disc’s optical depth is still evolving long past this point (Fig. 5a and b). The steady state, for which the optical depth maps no longer evolves, is reached only after ~10⁵ years.

The computation time to reach 10⁵  years evolution using the OPEN-MP version of the code on eight CPUs is ~65 h. The collisional procedure is by far the most costly in terms of CPU time. More precisely, for one typical time step the dynamical evolution procedure (to work out new orbits) takes up only ~1% of the total computation time, whereas the collisional procedure uses the remaining 99%. Within this collisional procedure, ~1.5% of the time is used to compute all tracer-tracer encounters and impact velocities, ~73% to compute the outcome (fragments) of these collisions, ~24% to sort out the dynamical families for all tracers and ~1.5% to select all the final tracers that will be kept at the end of the time step and to reassign the mass to these final tracers.

Note that we chose this simple unperturbed example because it has the pedagogical virtue of allowing one to easily distinguish the dynamical and collisional evolution of the tracers, that is, the spatial distribution of tracers stabilizes long before their collisional evolution does. In more complex cases (the ones that are really interesting to investigate with this code) such an easy distinction cannot be made, because the spatial distribution and the collisional evolution of tracers always change simultaneously (see for example the case study considered in Sect. 5).

4. Tests

Testing LIDT-DD was a challenge, as this code is the first of its kind, and there are no reliable results regarding the coupled dynamical and collisional evolution of debris discs that can be used as references. There are simplified cases, however, for which robust results have been obtained in past studies, which can be used as a benchmark to test the different aspects of our code.

4.1. Conservation of angular momentum

In the absence of external perturbers, the disc’s angular momentum needs to be conserved in a problem where all the forces are central (we do not consider PR-drag in this section). Figure 6 presents the evolution of the angular momentum derivative (dLog(L)/dt) for a few million years of our illustrating test run of an unpertubed ring (see Sect. 3.10). As can be seen, there are unavoidable small stochastic variations on short timescales due to the tracer selection and reattribution procedure. These variations remain very limited, however, less than 2 × 10^-14  yr^-1 in relative amplitude. More importantly, these stochastic variations do not increase in amplitude and induce no drift of the angular momentum over long timescales.

4.2. Mass loss

The collisional grinding of a debris disc naturally removes mass from it because of the blow-out of the smallest grains by radiation pressure. For an unperturbed system left to itself, the expected temporal evolution of both the system’s total and dust masses has been investigated in numerous studies. For an idealized system where the collisional cascade has had enough time (t > t_max) to reach the largest bodies in the size distribution, the expected behaviour is a decrease of M_tot and M_dust that is ∝t^-1 (Dominik & Decin 2003; Wyatt et al. 2007). However, this asymptotic behaviour is only expected to be reached at very late times, which can be longer than a system’s age (e.g., Löhne et al. 2008). For a more realistic case where t < t_max, detailed numerical and analytical investigations have shown that the evolution of M_tot and M_dust is much more complex, in particular because of the size dependency of the Q^∗ parameter in both the strength and gravity regimes (Löhne et al. 2008; Wyatt et al. 2011; Gáspár et al. 2013). We computed both M_tot and M_dust for our test case of an unperturbed disc and compared our results with that of these earlier studies.

The normalized total mass evolution displayed in Fig. 7 does closely match the behaviour that was obtained and thoroughly analysed by Löhne et al. (2008): during an initial stage, for which the collisional cascade has reached a steady-state only for bodies in the strength regime, M_tot remains almost constant. This phase ends when t ≥ t_b, where t_b is the time at which the collisional cascade has reached objects that are large enough for Q^∗ to be in the gravity regime. After that, M_tot decreases faster, at a rate that closely matches that predicted by Eq. (39) of Löhne et al. (2008). If we indeed substitute our own parameters for the initial size distribution and the Q^∗ dependency in this equation, we obtain M_tot ~ M₀(1−0.0030   t^0.3). As can be seen, this behaviour fits our test run at later times very well. Note that our best fit is obtained for 1.05M₀ instead of M₀, but this small difference is expected, as Löhne’s formula neglects all mass evolution during the t < t_b period.

As for the dust mass (Fig. 8), it first increases during the initial t ≤ t_b phase. This is an expected result, because the equilibrium size distribution of particles in the strength regime is steeper than the initial PSD (Gáspár et al. 2013). Beyond t_b, M_dust falls off approximately as ∝t^-0.41, which is again fully compatible with the rate predicted by (Löhne et al. 2008), that is, between t^-0.3 and t^-0.5 (see Fig. 10 of that paper).

4.3. Collisional energy dissipation

There is to our knowledge no reference analytic expression for how the kinetic energy should dissipate for complex systems where collisions and dynamics are coupled. Even for an unperturbed disc, there is no available law for high-velocity fragment-producing collisions affecting an extended size distribution, especially when taking into account radiation pressure. But we can test our code against available results obtained by Lithwick & Chiang (2007) for the simplified case of a ring made of monosize particles large enough not to be affected by radiation pressure. To mimic this monosize-particle case with LIDT-DD while retaining our collision-outcome prescription and tracer creation routine, we considered a system with ten size bins, but confined within a relatively narrow spread in sizes (all material that goes below the smallest bin is labelled as “lost dust”).

We take the same set-up as the nominal case considered by Lithwick & Chiang (2007), that is, a narrow ring at 1 AU from the central star, a mean tracer eccentricity of 0.01, and an optical depth τ₀ ~ 10^-2. We also chose the standard values ϵ_N = −0.3 and ϵ_T = 1 for the collision restitution factors. As an indicator of how the energy evolves in the system we chose to follow the evolution of the average eccentricity of tracers in the disc.

The ⟨ e(t) ⟩ evolution we obtain (Fig. 9) is very close to the one obtained by Lithwick & Chiang (2007) in their Fig. 1. for their hot case, with for example an initial decay time of ~2.5 × 10³ years to reach e₀/2.

Another important result from this graph is that the tracers’ dynamical evolution does not suffer from an artificial heating due to the finite cell sizes. Indeed, the values of ⟨ e ⟩ that are reached are far below the ratio Δ   r_cell/r ~ 0.25 between the width of the spatial cells and their radial distance. This proves the validity of our Keplerian-shear and azimuthal correction procedures when estimating relative velocities and collision outcomes between tracers of the same cell.

4.4. Viscous spread

Another consequence of the kinetic energy redistribution after collisions is the viscous spread of the disc, dissipative collisions playing the role of viscosity here (e.g., Lynden-Bell & Pringle 1974). For an unperturbed narrow ring, analytical derivations predict that the radial width of the ring should expand as ∝t^1/3 (Lithwick & Chiang 2007). We checked this behaviour for the no-radiation-pressure and narrow-size-distribution case considered in the previous subsection. As the typical timescale for the ring’s expansion scales as Δ³ (where Δ is the ring width) and ${\mathit{\tau }}_{\mathrm{0}}^{-1}$, and can thus be very long for wide and/or tenuous discs, we followed Lithwick & Chiang (2007) and considered a very narrow and dense initial ring, centred at 10 AU, of initial width Δ₀ = 0.05 AU. We considered the same optimum configuration as in Lithwick & Chiang (2007), that is, an initial eccentricity distribution such that the mean radial excursion of the tracers is equal to Δ₀.

As can be seen in Fig. 10, for a fixed initial set-up, the evolution of the ring’s width does closely match the expected behaviour in t^1/3. This is exactly the behaviour expected for a narrow ring evolving under dissipative collisions. Moreover, we also verified that the magnitude of the spreading is fully compatible with that viscosity being physical (i.e., collisional) rather than numerical. To do so, we followed Lithwick & Chiang (2007) and considered the simplifying assumption that the typical time for a disc to diffuse the width Δ is ${\mathit{t}}_{\mathrm{spread}}\mathrm{~}{\mathit{t}}_{\mathrm{col}}\hspace{0.17em}\hspace{0.17em}{\left(\frac{\mathrm{\Delta }}{\mathrm{⟨}\mathit{e}\mathrm{⟩}\mathit{.}\hspace{0.17em}\mathit{r}}\right)}^{\mathrm{2}}\mathit{,}$(6)where ⟨ e ⟩ is the average eccentricity of particles in the ring, and r the ring’s radial distance to the star. This equation is valid as soon as Δ is greater than ⟨ e ⟩ .   r. Replacing t_col by its approximate value t_orb/(2πτ), where τ is the optical depth in the ring when it has reached width Δ, we obtain ${\mathit{t}}_{\mathrm{spread}}\mathrm{~}\frac{{\mathit{t}}_{\mathrm{orb}}{\mathrm{\Delta }}^{\mathrm{3}}}{\mathrm{2}\mathit{\pi }{\mathit{\tau }}_{\mathrm{0}}{\mathrm{\Delta }}_{\mathrm{0}}\mathrm{(}\mathrm{⟨}\mathit{e}\mathrm{⟩}\mathit{.}\hspace{0.17em}\mathit{r}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(7)where we used τ ~ τ₀(Δ₀/Δ). Plugging in the parameters considered for the present set-up (r = 10  UA, Δ₀ = 0.05  UA, τ₀ = 0.01), we compared the estimates given by Eq. (7) with our numerical results displayed in Fig. 10. We found, for example, that the typical time to reach Δ = 0.25  UA should be ~12   000  years, which is relatively close to the ~10   000   retrieved from Fig. 10. Likewise, the time to reach Δ = 0.37   UA should be ~34   000  years, once again relatively close to the ~40   000   obtained in our run.

Note that for this ultra-narrow/very dense ring case, the spreading effect is much stronger than it is for the nominal fiducial ring considered elsewhere in this study (see Sect. 3.10). This is fully logical and is due to the fact that our nominal ring is both much wider, ~1  AU instead of 0.05 UA, and more tenuous, with τ₀ = 10^-3 instead of 0.01⁵.

4.5. Outward and inward flow of mass

One of the most innovative aspects of the code is the identification and sorting of dynamical categories within tracer populations and the subsequent tracer-reassignment procedure. As described in Sects. 3.7 and 3.8, these two procedures are relatively complex and together constitute one of the most critical parts of LIDT-DD.

A test for the reliability of these procedures can be made with the simplified set-up considered in Sect. 3.10 of an unperturbed inner ring of parent bodies. In this case, in regions outside the main ring there should be a two-way flow of small, radiation-pressure-affected grains produced in the annulus and placed on high-e orbits by radiation pressure: one flow of particles moving outward towards their apastron and another flow of particles moving inwards on their way back to their periaston in the ring. For a given particle size s₀, there should be an almost perfect balance between these two flows at steady-state, with only a small excess of outwardbound grains due to the small amount of s₀ grains that have been collisionally destroyed during their stay in the outer regions on their elongated orbits. In a fully deterministic (and fully unachievable) code where all numerical particles correspond to real physical bodies, this behaviour should be obtained automatically. Here, however, this behaviour is far from being straightforward, because particles are regrouped into tracers that can in principle be reshuffled, created, destroyed and/or merged at every collisional time step, depending on the number of dynamical categories that are identified and the number of redundant tracers within each of them.

We tested wether this result is obtained by displaying in Fig. 11 the mass contained in all tracers corresponding to the β = 0.4 size bin, which are present at a given time in all the azimuthal cells located in the fifth radial ring between 29 and 39 AU. We plot these tracers as a function of the periastron of their orbit, which gives the region of origin of the grains they represent. For tracers that have their periastron in one given radial region, we see that although there are significant variations in the mass carried by individual tracers, the total masses carried by outwardbound and inwardbound tracers (green and red lines) are remarkably close, with only the expected small excess in favour of outbound grains. The destruction rate, numerically estimated to be of the order of 10¹⁰ kg per time step, is so low (in outer parts of the disc τ < 10^-4) in this simulation that the red and green lines should be, and are, almost superposed.

4.6. Shape of the size distribution

Over the past decade, statistical particle-in-a-box models have identified several reliable results regarding the shape and profile of the particle size distribution (PSD) of collisional debris discs, in particular, how it might depart from the standard equilibrium distribution in s^-3.5dr (Dohnanyi 1969). We verified that all these robust features were obtained when our code is applied to unperturbed discs similar to those considered in those past investigations.

One well-established feature is the waviness of the PSD in the small-grain domain, which is caused by the depletion of unbound grains with β > 0.5. This depletion triggers a chain reaction in the PSD: there is an overabundance of grains that should have been destroyed by the absent unbound grains, that is, those just below the β = 0.5 limit, which in turn causes a depletion of slightly bigger grains that are destroyed by them, in turn causing an overabundance farther up the PSD, etc. (Campo Bagatin et al. 1994; Thebault et al. 2003; Krivov et al. 2006). This waviness is clearly visible in Fig. 12, displaying the PSD in the birth ring of the simplified inner-parent-body ring case considered in Sect. 3.10. The first peak in the wavy PSD is located at ~1.5   s_cutoff and the first dip around 10   s_cutoff, which agrees relatively well with the results obtained by the statistical code of Thebault & Augereau (2007, see Fig. 18 of that paper) for a similar set-up. The waviness is then progressively damped at larger sizes and is absent in the s > 100   s_cutoff domain; this result also agrees well with previous particle-in-a-box studies.

Another important feature of the PSD is that its slope is steeper than −3.5 in the strength regime, that is, up to objects in the sub-km range. In the 10^-4  m to 10²  m range we found dm/ds ∝ s^-0.64, corresponding to dN ∝ s^-3.64   ds. This −3.64 slope agrees remarkably well with the −3.7 slope found by Thebault & Augereau (2007) and the standard slope of −3.65 recently found in the extensive study by Gáspár et al. (2012).

Beyond s ~ 100  m, we clearly see another generic feature of debris disc PSDs, namely the flattening of the PSD corresponding to the transition from the tensile-strength to the gravity-dominated regimes (e.g., Löhne et al. 2008).

4.7. Radial surface density profile

Another well-established result is that beyond a narrow collisionally active ring, the surface density naturally tends towards a radial profile in r^-1.5. This somehow counter-intuitive feature⁶ arises because small high-e grains produced in the ring spend most of their time in almost-empty collisionally inactive regions far outside the ring, where they can safely accumulate, while they can only be destroyed during the small fraction of their orbit spent in the production ring (Strubbe & Chiang 2006; Thebault & Wu 2008). The time t_ss to reach this radial profile can be relatively long, of the order of the collisional timescale divided by the fraction of time f the smallest bound particles spend in the birth ring, ${\mathit{t}}_{\mathrm{ss}}\mathrm{~}\frac{{\mathit{T}}_{\mathrm{orb}}}{\mathrm{2}\mathit{\pi \tau }}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{f}}\mathit{,}$(8)where t_orb is the typical orbital period at the birth ring’s location, τ is the optical depth within the birth ring, and f is calculated according to the formula of Strubbe & Chiang (2006).

Figure 13 shows the result found with LIDT-DD simulating the evolution of such an initially confined ring (the set-up considered in Sect. 3.10) over t = 5 × 10⁵  years. It can be seen that a steady-state is reached after ~10⁵  years. This value is very close to the theoretical value t_ss ~ 1.1 × 10⁵  years that can be derived from Eq. (8) when taking t_orb for the middle of the birth ring at 11.5 AU and calculating f for the β ~ 0.4 particles that dominate the density profile in the outer regions. Once steady-state is reached, the slope of the radial density profile in the outer regions is close to the theoretical −1.5 value, in accordance with the conclusions of Strubbe & Chiang (2006) and Thebault & Wu (2008). Note, however, that the slope is in fact ~−1.7, that is, slightly steeper than −1.5. This small difference is expected and arises because the −1.5 slope is valid for an idealized case where all collisions beyond the ring are neglected, whereas some collisional activity (handled by LIDT-DD) occurs in these outer regions, thus removing some small grains and steepening the profile. After the onset of the quasi steady-state, the shape of the density profile no longer evolves, the only change being the progressive decrease of the total optical depth, which is due to the steady erosion of parent bodies in the birth ring.

4.8. Dynamically cold system

Another important result of recent debris disc studies is the peculiar PSD of debris discs that are dynamically cold, that is, where the stirring of orbits does not exceed e ~ 0.01. For this specific category of discs, Thebault & Wu (2008) found that there is a strong depletion of small grains, which arises from an imbalance between the rate at which these grains are produced and the rate at which they are destroyed. Indeed, while they are mainly produced from collisions amongst larger parent bodies, which will happen at low velocities (and thus produce few small fragments) if ⟨ e ⟩ is low, they are destroyed by impacts involving themselves, which will always occur at high velocities because the smallest grains are always placed on high-e orbits by radiation pressure regardless of the low dynamical stirring of the system (see discussion in Thebault & Wu 2008).

In Fig. 14 we display the results obtained with LIDT-DD for a dynamically cold disc for which e = 0.005. As expected, the depletion of grains in the s < 0.1mm range is clearly visible, and the geometrical cross-section is dominated by large grains. Note that for excitation values even lower than e = 0.005, the disc might enter a very different dynamically very cold regime for which impacts are no longer fragmenting and mutual accretion is possible.

Another important characteristics of dynamically cold discs that directly follows from the dominance of large grains is that they are expected to display sharp outer edges. This is exactly what we obtain when displaying the system’s radial profile in optical depth (Fig. 15), which abruptly drops just beyond the production ring. The cut-off in the profile farther out (~100  AU) is due to the limited radial excursion of the smallest particles (β = 0.44) in the simulation, because they are here created from parent bodies on nearly-circular orbits, so that their maximum apastron will be ~a₀/(1−2   β) = 100 AU, where a₀ = 12 AU is the maximum semi-major axis for a parent body in the main ring.

5. Case study: break-up of a massive object inside a debris disc

The detailed investigation of concrete astrophysical cases exceeds the scope of the present numerical paper and will be the purpose of future studies. However, to illustrate the potential of this new code for investigating cases that were so far inaccessible to detailed numerical modelling, we present here preliminary results for the case study of the violent break-up of a large planetesimal in an extended debris disc. Such violent and transient events have often been invoked to explain the presence of pronounced clumps in some resolved debris discs (Wyatt & Dent 2002; Telesco et al. 2005) or the high fractional dust luminosities of some discs, which are either too high or have too rapid variations to be explained with steady collisional cascades (Wyatt et al. 2007; Lisse et al. 2009; Melis et al. 2012). Until now, such massive collisions could not be directly simulated. Instead, they were usually modelled with analytical and statistical order-of-magnitude estimates (e.g. Kenyon & Bromley 2005; Lisse et al. 2009). An alternative approach is that by Jackson & Wyatt (2012), using an N-body scheme, resolving the initial impact with test particles, each representing a collisional cascade whose mass is steadily decreased using analytical laws.

The most sophisticated attempt at modelling catastrophic collisional events is probably that by Grigorieva et al. (2007). However, this model was unable to handle the case of a planetesimal exploding inside a disc, and could not be used on long timescales (see Sect. 1). The investigated case was that of a massive explosion far outside (much closer to the central star) the disc, from which only the unbound fragments were followed as they entered an outer disc at high velocities. This set-up allowed the authors to use a bimodal dust population (unbound explosion fragments and larger bound field particles) and only required short timescales.

5.1. Set-up

Because the goal of this example run is not to investigate specific real systems, such as HD172555 (Lisse et al. 2009) or TYC 8241 2652 1 (Melis et al. 2012), whose detailed study is deferred to a forthcoming dedicated paper, we consider here only one fiducial case, with one given set of initial parameters, chosen for its simplicity and illustrating virtues⁷.

We first considered an unperturbed debris disc, similar to the one considered in the test run of Sect. 3.10, that is, produced from a parent-body ring confined between 11 and 12 AU. The only difference is a lower optical depth ⟨ τ ⟩  =    1 × 10^-4 to enhance the contrast between the background debris disc and the fragments created from the impact.

We then considered the catastrophic explosion of a large body releasing a mass M_Ftot = 2    ×    10²² kg into fragments of sizes 2   μm < s < 10  cm, following a size distribution in dN ∝ s^-3.5ds. This explosion was assumed to take place at 35  AU from the star, and we assumed that the exploding parent body was on a circular orbit. The mass M_Ftot corresponded to that of an object of size ~10³  km ⁸. We implicitly assumed that the collision was just energetic enough for all M_Ftot to escape the potential of the shattered body. To ensure this, we followed the procedure of Jackson & Wyatt (2012) and released all fragments from a small circle corresponding to 100 physical radii of a ~10³   km body with a velocity kick randomly distributed between the escape velocity from the shattered body at this distance, ~150 m s^-1, and 1/4 of the surface escape velocity of this body, ~375 m s^-1.

All relevant initial parameters can be found in Table 2.

5.2. Results

Figure 16 plots the spatial location of all tracers at four different epochs with a colour scale giving the size bin corresponding to each tracer, while Fig. 17 shows the resulting smoothed 2D map of the system’s geometrical optical depth.

As can be clearly seen, an outward-propagating spiral structure builds up. It corresponds mostly to the smallest explosion fragments, which are pushed on high-eccentricity or unbound orbits by radiation pressure. After ~300  years, the smallest grains in the spiral arm have reached the outer limit of the plotted region at 200  AU. In parallel, the differential precession of the arm progressively populates the whole r > 35  AU region with explosion fragments. Because the precession rate is faster closer to the star, this populating process proceeds inside out. Meanwhile, at the radial location of the impact around 35AU, the largest fragments progressively form a homogeneous ring because of Keplerian shear.

The development of an outward-propagating spiral after a massive, fragment-releasing impact is a well-known feature that has been witnessed in previous N-body studies (e.g., Jackson & Wyatt 2012). However, we are here for the first time able to follow the collisional evolution of the released fragments as they propagate through the disc, as well as the collisional fate of the fragments that stay in the radial region of the explosion. In particular, the code takes into account the production of new generations of collisional fragments from impacts involving the initially released primordial fragments. Such a second, third or nth generation debris can significantly contribute to the disc’s luminosity in some regions.

LIDT-DD’s handling of collisions gives us the possibility of estimating the timescale during which the signature of the explosion remains visible in the disc. We see that the system settles back to an unperturbed 2D profile in two steps. In the outer regions, the propagating spiral and its aftermath are no longer visible, meaning that their signature drops below the noise level in the 2D map after ~2000   years (Fig. 17g). The disappearence of the spiral pattern in ~2000   years is a purely geometrical effect, because this time is approximately the orbital period of the small, high-β grains that were released by the initial shattering (this geometrical effect has also been identified by Jackson & Wyatt 2012). But in contrast to the collisionless case of Jackson & Wyatt (2012), for which, even after the fading of the spiral, the presence of small released grains remained visible in the outer regions for a rather long time, we found here that collisional destruction of these particles prevents them from leaving an observable trace after these 2000 years have elapsed. The situation is very different at the radial location of the initial release, where the secondary ring made of the largest primordial explosion fragments is still clearly visible. This secondary ring at the explosion’s location is much more long-lived. It takes another few million years before it is resorbed in the background level (Fig. 17h). This slow fading of the secondary ring is due to the progressive grinding down of the largest explosion fragments⁹, settling into a collisional cascade that produces increasingly smaller debris, which eventually are blown out by radiation pressure.

Note that these two timescales significantly depend on our arbitrary choice of initial parameters, mainly the value for M_Ftot and the size of the largest explosion fragments, and we repeat that this fiducial set-up was not chosen so as to match a specific real astrophysical system. The important point here is that once a set-up has been chosen, depending on the case that is to be investigated, then LIDT-DD is able to provide a reliable estimate for the survival time of these transient signatures of massive breakups. It is also able to evaluate how the luminosity profiles evolve with time in the wake of the initial violent event.

5.3. Tracers vs. real disc

The release of the explosion fragments results in a significant increase of the number of tracers in the system (Fig. 16). Note that this increase does not reflect the amount of matter that is released in the disc, but the fact that the smallest explosion fragments, which rapidly populate the outer regions because of radiation pressure, have different dynamical characteristics than the matter at rest in the unperturbed debris disc. Since the code’s structure requires that each collisional cell must contain at least two tracers for each identified dynamical category (see Sect. 3.7), this means that additional tracers will populate the cells. These tracers can either be primordial ones released at the moment of the explosion, or can have been created later as the primordial fragments pass through the outer disc and spawn new subfragments that need to be represented by new tracers.

As such, the number of tracers that is needed in the system is uncorrelated to the amount of mass that is initially released. For example, the number of tracers and their locations would be exactly the same for a simulation with M_Ftot/10 instead of M_Ftot. The difference would be in the number of physical particles each tracer stands for. This is also why the number of tracers does not significantly decrease, even long after the effect of the explosion is no longer visible in the real outer disc (see Figs. 16c and d), simply because a dynamical category, once created in a cell, is maintained as long as even a little amount of matter corresponds to it. In the present case, because there is always some collisional production of new fragments from the excess matter remaining around 35 AU, there is always some matter that is produced with dynamical characteristics that roughly match those of the first initially released fragments. True, the number of these newly injected fragments will decrease with time, and they are no longer visible outside of 35 AU on the 2D maps after ~2000 years, but the code is still able to identify and handle them even when they do not show up in the global density maps, hence the still large number of tracers at this stage (see Fig. 16d).

6. Limitations

Despite its many new features and the improvement it brings to disc studies, LIDT-DD, as all numerical codes, has unavoidable limitations. The main ones are related to the complexity of the numerical treatment itself, and especially the tracer-tracer collision treatment procedure (and the tracer creation and reassignment that follows), which is very time-consuming. This limits the number of tracers that can be handled to typically a few 10⁵. While this number allows a relatively good spatial sampling, it is at least an order of magnitude lower than the number of test particles that can be typically handled in pure N-body codes, or more exactly, N-body codes where the restricted three-body case applies (e.g. Reche et al. 2009), so that the spatial resolution will always be lower than in these purely dynamical codes. This also limits the timescale that can be investigated with the present LIDT-DD version to a few 10⁶ years. While this is enough to study the evolution of young and bright discs such as β Pic or HR4796A, it does not allow one to follow the evolution of older systems over their entire lifetime. However, this limitation might not be too constraining, as many processes shaping the structure of even very old discs might develop and evolve on timescales shorter than a few million years.

The fact that only two tracers per dynamical family are kept and that virtual tracers are assigned the same velocity as the larger of the two initial colliding tracers necessarily introduces small errors in the velocity distribution with every time step. These small errors are clearly visible in the form of short-term variations of the system’s angular momentum (Fig. 6). However, they remain limited in their amplitude (lower than ~0.2%) and have no cumulative effect, since they do not lead to systematic variations of the angular momentum over long timescales.

The limited number of tracers also introduces a relatively high noise in the density maps that are obtained. This is because of the finite distances between neighbouring tracers, which can be relatively large in scarcely populated regions. This noise can be overcome by averaging the surface density over larger regions, but the price to pay in this case is a loss in spatial resolution. Note, however, that this problem mainly affects the regions that are the most dynamically quiet, whereas dynamically perturbed regions, which are usually of main interest in concrete studies, are much more densely populated in terms of tracers, thus allowing a much better spatial resolution (see Fig. 17).

Another limitation is that a 2D (r,θ) collisional grid is used to compute size distributions. This implicitly assumes a symmetry in the Z direction for collision rates, an approximation that might be faulty when considering non-planar cases such as an inclined planet. A 3D collisional grid extending in the Z direction as well would solve the problem, but would multiply all calculations by at least a factor N_z (N_z being the number of cells in the Z direction at a fixed r and θ). Note, however, that the dynamical evolution of the system is always 3D.

7. Summary and perspectives

LIDT-DD is a new code developed for studying debris discs, based on a global approach that couples an N-body handling of the dynamics to a particle-in-a-box treatment of collisions. It takes its basic architecture from the LIDT3D code developed by Charnoz & Taillifet (2012) for the very different context of protoplanetary discs, with tracers that each represent a whole population of particles of a given size at a given location in the system, whose collisional interactions are handled, with a statistical approach, within superimposed spatial cells. It can, however, be considered as an independent stand-alone code, because of the major modifications and upgrades introduced to account for the very constraining specificities of debris disc physics. The three major such constraints are that 1) impact velocities are much higher in debris discs, so that collisions are highly destructive and fragment-producing; 2) the dynamics is much more complex because, in contrast to protoplanetary discs, in the absence of the smoothing effecf of gas drag, initial conditions are not relaxed so that all grains retain the characteristics of where and how they have been produced. As such, there is a very wide range of possible dynamics for same-sized particles located in the same region of the disc; and 3) the collisional and dynamical evolution of small grains close to the blow-out limit, which usually dominate the disc’s luminosity, is extremely sensitive to small size variations.

The main specificities of the new code are the following:

• All mutual tracer-tracer impacts are handled individually, and the produced collisional fragments have orbits derived from those of their collisional progenitors.

• The collision outcome procedure is now able to handle high-velocity fragmenting impacts and is similar in sophistication to those implemented in classical particle-in-a-box codes (e.g., Krivov et al. 2006; Thebault & Augereau 2007).

• In addition to sorting tracers by size and spatial location, the code is able to sort and regroup them by dynamical families, using a specially developed hierarchical clustering procedure. This allows LIDT-DD to track down and keep a trace of all the possible dynamical origins for grains present in a given region of the disc.

• The procedure for tracer creation and reassignment is designed to preserve two tracers per dynamical family in each spatial cell and each size bin. This prevents both the accidental removal of important dynamical information and the undesired increase in the number of tracers used by the code.

In the absence of any reference result for the coupled dynamical and collisional evolution of discs that could be used as reliable comparison, LIDT-DD was tested on three simplified cases for which reliable results have been obtained in past studies:

• We were able to retrieve the main features for the particle size distribution of unperturbed collisional debris discs, such as its waviness in the small grain domain and its steeper-than s^-3.5 slope in the millimetre-to-subkilometre range. The value we found for this slope in our nominal test case was −3.64, very close to the standard value of −3.65 inferred by Gáspár et al. (2012).

• For an initially narrow ring-like disc, we reproduced the results of Thebault & Wu (2008), showing that the surface density profile beyond the initial ring converges towards a radial decrease in ~r^-1.5.

• For a dynamically cold disc with particles on very low orbital eccentricities, we confirmed the result obtained by Thebault & Wu (2008) and Löhne et al. (2012) that such a disc is depleted of small grains and dominated by particles bigger than ~10   s_cutoff.

In addition to these comparison tests, we also independently tested the reliability of our tracer creation and reassignment procedure by verifying that the mass fluxes of outgoing and ingoing particles in a collisional disc at steady-state do almost perfectly balance one-another. Moreover, for simplified unperturbed cases, we verified that the angular momentum of the disc is conserved and that collisional energy dissipation leads to the expected decrease of the rms eccentricity of tracers as well as to the expected viscous spread of the disc.

As an illustration of the code’s potential to tackle astrophysical cases that are as yet numerically unexplored, we considered the case study of the break-up of a massive planetesimal in an extended disc. LIDT-DD was able to follow the collisional fate of the smallest grains released by the break-up that are pushed by radiation pressure on high-e or unbound orbits as they propagate outward in the disc. We see the expected formation of outward propagating spiral arms, which are resorbed by collisions and dynamical dilution after a few dynamical periods. The code was also able to monitor the fate of the large amount of mass that remains at the location of the initial break-up. This matter rapidly formed a ring-like structure that was much longer-lived than the outbound spiral and was resorbed into the disc’s background after a few million years. This was the first time that a code was able to estimate the amplitude and survival time of the signatures such events can leave in a disc.

An exhaustive investigation of this crucial question of massive transient collisional events, exploring the wide range of possible set-ups (location of the break-up, released mass, disc density, etc.) exceeds the scope of the present code-introducing paper. This problem will be addressed in a forthcoming study, whose main objectives will be to explore individual systems with abnormal flux excesses, such as HD172555 (Lisse et al. 2009; Johnson et al. 2012) or TYC 8241 2652 1 (Melis et al. 2012), and also to investigate how generic the massive-break-up scenario can be for explaining all bright debris discs with luminosities that cannot be explained by classical collisional cascades.

Another area of interest for LIDT-DD are the planet-disc interactions and the extent to which planetary companions can sculpt debris discs. This question has recently been explored by both the CGA and DyCoSS algorithms (e.g. Kuchner & Stark 2010; Thebault et al. 2012). Although these studies have given promising results, these codes’ structures limit them to relatively restrictive set-ups, that is, systems at steady-state, with no more than one planet and with collisional prescriptions assuming that all impacts are fully destructive, thus neglecting all second-generation fragments. LIDT-DD is not bound by these limitations and will thus enable the exploration of transient events, multi-planet systems as well as the feedback of the planetary perturbations on the collisional evolution. We stress that the CGA and DyCoSS codes are by no means obsolete, as they enable a very high spatial resolution that cannot be attained with the present code for the time being (see Sect. 6).

Another potential application of the code is the puzzling case of bright exozodiacal discs that have been identified by interferometry very close to several main-sequence stars such as Vega or Fomalhaut (e.g. Absil et al. 2009; Defrère et al. 2011). The estimated dustiness of these hot exozodis is in most cases far too high to be explained by steady collisional cascade, and alternative scenarios have to be considered (Mennesson et al. 2013). One possible scenario could be the injection of material from a planetesimal belt farther out in the disc that is dynamically destabilized by a perturbing planet. This scenario has been proven to be dynamically viable, at least for young systems, by the N-body investigation of Bonsor et al. (2012). However, it remains to be seen how this inward-scattering of large planetesimals can translate into observable small dust. This crucial question directly depends on the collisional fate of these bodies and can be investigated by LIDT-DD. Alternative scenarios for explaining exozodis, such as transient massive impacts, falling evaporating bodies (Thebault & Beust 2001; Beust & Valiron 2007), or pile-up due to the complex interplay of PR drag, sublimation and collisions (Kobayashi et al. 2011), can in principle also be investigated with LIDT-DD. More generally, LIDT-DD enables the handling of all complex scenarios where both dynamics and collisions are expected to play an important role in a debris disc evolution.

---

Acknowledgments

We thank Jean-Charles Augereau and Amy Bonsor for fruitful discussions. We thank the anonymous referee for helping improving the quality of the paper. Q.K. acknowledges financial support from the French National Research Agency (ANR) through contract ANR-2010 BLAN-0505-01 (EXOZODI).

References

All Tables

All Figures