A&A 461, 537-549 (2007)
DOI: 10.1051/0004-6361:20065210

Collisional dust avalanches in debris discs

A. Grigorieva¹ - P. Artymowicz^1,2 - Ph. Thébault^1,3

1 - Stockholm Observatory, SCFAB, 10691 Stockholm, Sweden
2 - University of Toronto at Scarborough, 1265 Military Trail, Toronto, Ontario, M1C 1A4, Canada
3 - Observatoire de Paris, Section de Meudon, 92195 Meudon Principal Cedex, France

Received 15 March 2006 / Accepted 29 August 2006

Abstract
We quantitatively investigate how collisional avalanches may develop in debris discs as the result of the initial breakup of a planetesimal or comet-like object, triggering a collisional chain reaction due to outward escaping small dust grains. We use a specifically developed numerical code that follows both the spatial distribution of the dust grains and the evolution of their size-frequency distribution due to collisions. We investigate how strongly avalanche propagation depends on different parameters (e.g., amount of dust released in the initial breakup, collisional properties of dust grains, and their distribution in the disc). Our simulations show that avalanches evolve on timescales of $\sim$1000 years, propagating outwards following a spiral-like pattern, and that their amplitude exponentially depends on the number density of dust grains in the system. We estimate the probability of witnessing an avalanche event as a function of disc densities, for a gas-free case around an A-type star, and find that features created by avalanche propagation can lead to observable asymmetries for dusty systems with a $\beta$ Pictoris-like dust content or higher. Characteristic observable features include: (i) a brightness asymmetry of the two sides for a disc viewed edge-on, and (ii) a one-armed open spiral or a lumpy structure in the case of face-on orientation. A possible system in which avalanche-induced structures might have been observed is the edge-on seen debris disc around HD 32297, which displays a strong luminosity difference between its two sides.

Key words: stars: circumstellar matter - planetary systems: formation - planetary systems: protoplanetary disks - stars: individual: $\beta$ Pictoris - stars: individual: HD 32297

1 Introduction

Direct imaging of circumstellar discs (e.g., Heap et al. 2000; Schneider et al. 2005; Liu 2004; Clampin et al. 2003) have provided resolved disc morphologies for several systems (e.g., $\beta$ Pic, HD 141569A, HD 100546, HD 32297) and have shown that dust distribution is not always smooth and axisymmetric. Warps, spirals, and other types of asymmetries are commonly observed (e.g., Kalas & Jewitt 1995, for the $\beta$ Pic system). These morphological features can provide hints on important ongoing processes in the discs and improve our understanding of the evolution of circumstellar discs and of planetary formation.

The usual explanation proposed for most of these asymmetries is the perturbing influence of an embedded planet. As an example, the warp in the $\beta$ Pic disc has been interpreted as induced by a jovian planet on an inclined orbit (Mouillet et al. 1997; Augereau et al. 2001). Likewise, for annulus-like discs with sharp inner or outer edges, the most commonly proposed explanation is truncation or gap opening due to planets or bound stellar companions (e.g., Augereau & Papaloizou 2004), although alternative mechanisms such as gas drag on dust grains within a gas disc of limited extent have also been proposed (Takeuchi & Artymowicz 2001). For spiral structures, authors have also been speculating on gravitational instabilities (Fukagawa et al. 2004), as well as on a bound stellar companion (Augereau & Papaloizou 2004).

The catastrophic breakup of one single large object releasing a substantial amount of dust fragments could be an alternative explanation for some of the observed asymmetries. Wyatt & Dent (2002) have examined how such collisionally produced bright dust clumps could be observed in Fomalhaut's debris disc. Likewise, such clumps have been proposed by Telesco et al. (2005) as a possible explanation for mid-infrared brightness asymmetries in the central $\beta$ Pictoris disc, but only based on preliminary order of magnitude estimates. More recently, the detailed study of Kenyon & Bromley (2005) investigated the possibility of detecting catastrophic two-body collisions in debris discs and found that such a detection would require the breakup of 100-1000 km objects. The common point between these different studies is that they focus on global luminosity changes due to the debris cloud directly produced by the shattering events themselves.

In the present paper, we re-examine the consequences of isolated shattering impacts from a different perspective, i.e., by considering the collisional evolution of the produced dust cloud after its release by the shattering event. The main goal here is to study one possibly very efficient process, first proposed by Artymowicz (1997), but never quantitatively studied so far, i.e., the so-called collisional avalanche mechanism. The basic principle of this process is simple. After a localized disruptive event, such as the collisional breakup of a large cometary or planetesimal-like object, a fraction of the dust then produced is driven out by radiation pressure on highly eccentric or even unbound orbits. These grains moving away from the star with significant radial velocities can breakup or microcrater other particles farther out in the disc, creating in turn even more small particles propagating outwards and colliding with other grains. Should this collisional chain reaction be efficient enough, then a significant increase in the number of dust grains could be achieved. In this case, the consequences of a single shattering event, in terms of induced dust production, could strongly exceed that of the sole initially released dust population. The outward propagation of the dusty grains could then induce observable asymmetric features in the disc, even if the initially released dust cloud is undetectable.

The goal of this work is to perform the first quantitative study of the avalanche process and investigate the morphology of avalanches in debris discs, under the assumption that dust dynamics is not controlled by gas (Lagrange et al. 2000). For this purpose we have created a numerical code, described in Sect. 3, that enables us to simulate the coupled evolution of dynamics and size-frequency distribution of dusty grains. The results of our simulations, which explore the effect of several parameters (total mass and radial distribution of dust in the disc, mass and size distribution of the planetesimal debris, physical properties of the grains and the prescription for collisional outcome for grain-grain collisions) are presented in Sect. 4. In Sect. 5 we examine under which conditions avalanche-induced features might become observable. We end with a discussion of the probability of witnessing an avalanche (Sect. 6) and finally a summary (Sect. 7).

   
2 Simplified theory of dust avalanches

A dust avalanche is a chain reaction of outflowing debris impacting disc particles and creating even more debris accelerated outwards by the star's radiation pressure. The basic principle of this mechanism can be illustrated by a set of analytical equations. We present here a simplified theory of avalanches based on the order-of-magnitude approach of Artymowicz (1997), firstly for its pedagogical virtues, but also because it can serve as a reference that facilitates the understanding of the main results derived from our extensive numerical exploration.

Let us assume that N particles of size $s_{\rm gr}$ (radius) move through a cloud of dust grains of size s at a relatively high velocity. Let us further assume that each collision produces a constant number  $N_{\rm\beta}$ of such debris, which are quickly accelerated to velocities leading to further destructive collisions. To derive the total number of debris produced by the avalanche, we define the optical depth as

$$% {\rm d}\tau= n(s) \sigma(s) {\rm d}l,$$ (1)

where n(s) is the number density of dust particles of size s in the system, $\sigma(s) \approx \pi (s+s_{\rm gr})^2$ is the cross-section for collisional interaction between grains, and ${\rm d}l$ is the length measured along the grain path. The number of debris produced in the interval  ${\rm d}\tau$ is then

$$% {\rm d}N= N N_{\rm\beta}~ {\rm d}\tau.$$ (2)

Integration over the whole path of the grains gives the total number of debris produced by the avalanche

$$% N_{\rm tot}= N_0{} \exp (N_{\rm\beta} \tau_{\rm }),$$ (3)

where N₀ is the number of outflowing grains initially released.

In a disc, $\tau$ can be approximated by the optical thickness in the disc midplane,

$$% \tau_{\rm\Vert}=\int\!\!\int \pi s^2 {\rm d}n(s) {\rm d}R,$$ (4)

where R is the radial cylindrical coordinate. We replace $N_{\rm\beta}$ by its average value $\langle N_{\rm\beta} \rangle$ to emphasize the fact that in reality  $N_{\rm\beta}$ depends on the details of each collision. Equation (3) then takes the form

$$% N_{\rm tot} \sim N_0{} \exp( \langle N_{\rm\beta} \rangle \tau_{\rm\Vert}).$$ (5)

This equation gives an estimate of an avalanche efficiency in a disc through the total number of grains  $N_{\rm tot}$it produces. However, one should keep in mind that the relevance of this set of equations is limited to global, order of magnitude estimates. Furthermore, these equations do not give any insight into the temporal development and spatial structure of a given avalanche. For these crucial issues, numerical modeling is clearly required.

   
3 The model

The number of dust grains in a circumstellar disc is far too large to follow every grain individually during the calculation; some kind of statistical approach must therefore be used. Models of dust disc evolution developed to date fall into two main categories. On the one hand, "particle in a box'' models divide the dust grains into statistical bins according to their size and enable us to compute the evolution of the size distribution within a given spatially homogeneous region (e.g., Thébault et al. 2003). While it is possible to mimic a spatially inhomogeneous system by integrating a set of coupled particle-in-a-box models, this can become unwieldy in the absence of strong simplifying symmetries. Kenyon & Bromley (2004) use a multiannulus code for example, but their model is one-dimensional in space. On the other hand, direct N-body simulations (treating the dust as test particles in the potential of a 2 or multi-body system) are used to accurately follow the spatial evolution of dynamical structures such as planet induced gaps or resonances (e.g., Augereau & Papaloizou 2004; Wyatt 2003). In this case, however, the sizes of the dust grains are either not taken into account or assumed to be equal.

For the present problem, however, we need to follow both the spatial distribution of the grains and their size distribution with reasonable accuracy. To do this, we developed a new code in which all grains with similar parameters (size, chemical composition, spatial coordinates, and velocity) are represented by a single superparticle (hereafter SP). We follow the dynamical evolution of these SPs and compute the collisional destruction and production of grains as SPs pass through each other. We represent newly created grains as new SPs. The maximum number of SPs our code can handle is about one million.

3.1 Superparticles

A detailed description of our SP modeling is given in the appendix. Here we briefly outline its main characteristics: a SP is described by the position and velocity of its center of mass (which coincides with its geometrical center), by its size, shape, and internal density profile, and by the number of dust grains it contains. For the present work all SPs are treated as cylinders and their geometrical centers are constrained to lie in the midplane of the disc. The cylinders have constant radii  $r_{\rm sp}$ and variable heights  $h_{\rm sp}(R)$, where R is the distance from the star (see Appendix A). All grains inside a given SP are assumed to have the same physical properties. We assume that all grains in our simulation are spheres with identical densities, chemical compositions, and porosities. The grains (and thus the SPs) are distributed into mass bins separated by a factor 2 logarithmic mass increment (i.e., a factor of 1.26 in size).

The trajectory of a SP corresponds to the trajectory of a test particle (with dynamical properties identical to the SP's grains) located at the SP's center of mass (see Sect. 3.2). SPs can overlap and freely pass through each other. In this event, collisional interactions between their respective grain populations is considered. This process is treated as a passage of two clouds of grains through each other (see Appendix A.3). It results in the loss, by destructive collisions, of a fraction of the initial grain populations and the production of smaller collisional fragments. These newly produced debris are placed into newly created SPs in accordance with the grain sizes and velocities. In the current version of the code, the centers of all SPs move in the same plane and the dust distribution is symmetric with respect to this midplane. However, the SPs representation method could in principle be used to model systems with vertical deformations (e.g., warps).

The size of a SP is fairly large ( $r_{\rm sp}=5$ AU). This puts unavoidable constraints on the spatial resolution of our simulations and prevents us from modeling processes occurring on scales smaller than the SP radius. It would, for example, be difficult to model fine resonant structures induced by disc-planet interaction. Moreover, the current version of the method with a constant value of the SP radius is not applicable to collisional evolution in the inner regions ($\la$20 AU) of debris discs. Although this limitation could be overcome by introducing a dependence of the size of a SP on the distance to the star (e.g., $r_{\rm sp} \propto R$), we have not implemented it in the current version of the code, since our main goal here is to model collisional avalanches that propagate outwards, inducing observationally significant features in the outer ($\ga$100 AU) regions of the disc.

The grains inside a SP do not have explicit vertical velocity components. To check the validity of this assumption, we have performed test runs, for which an artificial vertical velocity dispersion term was added to the planar velocity, which showed no significant departure from the in-plane velocities case. Note that a vertical velocity component is, however, indirectly taken into account by the fact that SP heights increase with distance from the star (see Appendix A.3), accounting for the geometrical dilution of grain spatial densities.

   
3.2 SP trajectories

As has been mentioned earlier, the trajectory of a SP is identical to the trajectory of a test particle (with mass, size, and chemical composition identical to those of the SP's grains) located at the SP's center of mass. Test particles move in the gravitational field of a star under the influence of the stellar radiation force. The equation of motion reads:

$$% m\frac{{\rm d}^2\vec{r}}{{\rm d}t^2}= - \frac{GMm}{r^3}\vec{r}+\vec{F}_{\rm rad}+\vec{F}_{\rm PR},$$ (6)

where m and $\vec{r}$ are the mass and position of the test particle, G is the gravitational constant, M is the star mass, and $\vec{F}_{\rm rad}$ and $\vec{F}_{\rm PR}$ are the radiation pressure and Poynting-Robertson drag, respectively. In our simulations we can neglect the Poynting-Robertson drag, since it acts on a timescale much longer than the time intervals considered here.

The radiation pressure force is expressed as a function of the gravitational force through the radiation pressure coefficient, $\beta$, as

$$% \vec{F}_{\rm rad} = - \beta \vec{F}_{\rm grav} = \beta \frac{GMm}{r^3}\vec{r}.$$ (7)

The parameter $\beta$ is a function of the stellar luminosity, grain size, and optical properties of the grain material (Burns et al. 1979). We use the Mie theory code developed by Artymowicz (1988) to calculate $\beta$ (Fig. 1).

A 7th-8th order Runge-Kutta method is used for integrating test particles trajectories. Although in the simulations presented in this paper the dynamics of the SPs is purely Keplerian, we have decided not to use analytical solutions since the Runge-Kutta integrator allows for an easy inclusion of any additional gravitational (due to planetary or stellar perturbers) or dissipative forces (such as PR and gas drag).

Figure 1: Ratio of the radiation pressure force to the gravitational force vs. grain size for different grain materials and porosities, P, calculated for an A5V ($\beta$ Pictoris-like) star. The thick lines represent silicate Mg $_{\rm0.95}$Fe $_{\rm0.05}$SiO$_{\rm 3}$. The thin lines are used for water ice. The solid lines are for solid grains, the dashed lines are for 50% porous grains, and the dash-dotted lines are for 80% porous grains.

Open with DEXTER

3.3 Initial dusty disc structure

The SP representation method is used to model the initial dusty disc structure. The total number of SPs for each size bin is chosen so that, at any given location in the disc, there are at least 2-5 overlapping SPs to account for different dynamical characteristics of grains of this size at this location. Each of these overlapping SPs thus differs from the others by its local velocity. To model the initial dust distribution in the disc we use $\simeq$5 $\times$ 10⁴ SPs (test runs with larger number of SPs do not lead to significant changes in the results). The number density of dust grains at a given location in the disc is calculated as the sum of the grain densities of the overlapping SPs.

The archetypical, and still by far the best known, debris disc of $\beta$ Pictoris is taken as a reference system for the initial dusty disc structure. In the present study we do not aim to model this particular system and just adopt its global properties for the dust distribution. Alternative dust distributions are also explored in Sect. 4.5. For the dust profile in $\beta$ Pictoris, we take the results of Augereau et al. (2001), who numerically derived the dust distribution giving the best fit to the resolved scattered light images as well as the long-wavelength photometric data, as a reference. We assume here that all grains are produced from parent bodies on circular orbits following the best-fit parent body distribution given in Augereau et al. (2001), where most of the bodies are located within an extended annulus between 80 and 120 AU, with a depletion in the inner <50 AU region and a sharp drop of the density distribution outside 120 AU (see for example Fig. 1 of Thébault & Augereau 2005). Grains with small $\beta$ have almost the same orbits as their parent bodies (the biggest grains), while smaller grains (i.e., with higher $\beta$) have more elliptic orbits depending on their $\beta$ value. The initial number of grains as a function of their size follows a classical single power law size-frequency distribution

$$% {\rm d}n=N_{\rm0s} (s/s_0)^{-p_{\rm0s}}{\rm d}s$$ (8)

where the power-law coefficient, $p_{\rm0s}=3.5$, corresponds to an idealized collisionally evolved system (Dohnanyi 1969). The minimum size  $s_{\rm min}$ for the disc grains is given by the radiation pressure blow-out cutoff and corresponds to $\simeq$$2~\mu$m for the compact grains considered in the nominal case. The maximum size  $s_{\rm max}$ for the disc grains is taken to be 1 cm. Runs with an order of magnitude higher maximum grain size give very similar results while being more computationally demanding. At the same time we cannot lower  $s_{\rm max}$ since millimeter particles make up a few percent of the disc's optical thickness and their contribution starts to be significant for the avalanche development.

The vertical structure of the disc is expressed in terms of the vertical geometrical optical thickness, $\tau _{\rm\perp }$, per unit length, z, as

$$% \frac{{\rm d}\tau_{\rm\bot}}{{\rm d}z} =C_{\rm\tau} \frac{\... ...eft[ - \left(\frac{\vert z\vert}{w}\right)^{p_{\rm z}}\right],$$ (9)

where $C_{\rm\tau}$ is a normalizing constant, $\tau_{\rm\bot}(R)=\int\!\!\int \pi s^2 {\rm d}n(R,s) {\rm d}z$ is the vertical optical thickness of the disc at distance Rfrom the star, n(R,s) is the number density of dust grains of size s, w(R) is the disc width, and $p_{\rm z}=0.7$ is a parameter, that determines the shape of the vertical profile. The disc width changes with radius as

$$% w_{\rm }(R)=0.055 R_{\rm m} \left(\frac{R}{R_{\rm m}}\right)^{p_w},$$ (10)

where $R_{\rm m}=117$ AU and p_w=0.75 for most of the runs. (Alternative values of p_w have been explored in test runs, which have shown that results only weakly depend on it.)

   
3.4 Collisional outcomes

Collisions are the crucial mechanism for the development of the avalanche phenomenon. The result of a collision, in terms of the size-frequency distribution of the debris, depends on several parameters: projectile and target materials and structures, sizes, impact velocities, and angle of incidence. Since it is not possible to model every collision in such detail, we have to adopt a simplified algorithm. We assume that the impact energy of colliding bodies, $E_{\rm col}$, is equally shared between them. Laboratory experiments show that this is the case when both bodies are made from identical material regardless of their sizes (Ryan et al. 1991). $E_{\rm col}$ is:

$$\begin{eqnarray*}E_{\rm col}=\frac{M_1 M_2 v_{\rm rel}^2 }{ 2(M_1+M_2)}, \end{eqnarray*}$$

where $v_{\rm rel}$ and $M_{\rm 1,2}$ are the relative velocity and the masses of the colliding bodies. The relative velocity between grains is derived from relative velocities between SPs after compensating for the artificial Keplerian shear induced by the SP finite radius.

Collision outcomes are traditionally divided into two classes: (1) catastrophic fragmentation, when the largest remaining fragment, $M_{\rm lf}$, is less than half of the parent body mass, M, and (2) cratering, when $M_{\rm lf} > 0.5M$. The energy per unit mass that is needed to get $M_{\rm lf} = 0.5M$ is called the threshold specific energy Q^*. If the specific energy $Q=0.5E_{\rm col}/M$ received by a body is more than Q^*, then the collision leads to catastrophic breakup, whereas cratering occurs if Q < Q^* (Fujiwara et al. 1977; Benz & Asphaug 1999; Petit & Farinella 1993). Q^* is a function of size for which we adopt a classical power law dependence (e.g., Ryan & Melosh 1998; Housen & Holsapple 1999). The collisional response of the small objects considered in the present work falls into the so-called strength regime, where the target's internal strength is the dominant factor, for which

 

Q^*=Q^*₀(s/s₀)^-p_Q, (11)

where Q^*₀ corresponds to the value of the threshold energy for size s₀. In this regime, Q^* decreases with size. Housen & Holsapple (1990) and Ryan & Melosh (1998) present a wide range of values for p_Q. We cannot directly apply their result because we are dealing with much smaller sizes. As pointed out in Tielens et al. (1994) "for submicron-sized bodies, cracks play little role, and the strength of a material approaches the ultimate yield strength of the material'', which corresponds to Q^*=2 $\times$ 10⁸ ergs/g (Tielens et al. 1994, and references therein). The threshold energy size-dependence most probably has a knee in the micron-submillimeter size range, but since we do not have any information about the slope change we decided not to introduce two additional unknown parameters into the code, but simply to adopt a slightly shallower slope, p_Q=0.2, for our calculations.

To account for the effect of different incidence angles, we correct the value of Q^* by a correction factor  $x_{\rm cr}$corresponding to an average over all incidence angles

$$\begin{eqnarray*}Q^*=Q^*_{\rm head~on}/x_{\rm cr}, \end{eqnarray*}$$

where $x_{\rm cr}=0.327$ (Petit & Farinella 1993). For both catastrophic shattering and cratering prescription, we use the approach and the algorithm presented in Petit & Farinella (1993). However, we refine this model by assuming that the fragment mass distribution produced follows a broken power-law instead of a single-index one:

 

$$% {\rm d}n=N~ {\rm d}m\left\{ \begin{array}{ll} \left( {m}/{m_{\r... ...}/{m_{\rm s}}\right)^{-q_2} &\mbox{if $m \geq m_{\rm s}$ .} \end{array} \right.$$ (12)

Such a change of slope between the small and large fragments, always corresponding to a flattening of the slope in the small particle range, is indeed supported by experimental results (e.g., Davis & Ryan 1990). For the choice of values for the power-law indexes (q₁, q₂) and the transition mass $m_{\rm s}$ for the slope change, we take the experimental studies of Davis & Ryan (1990) as a reference and explore values within the range of possible values obtained by these authors. The minimum fragment size is assumed to be $0.1~\mu$m, unless otherwise explicitly specified.

In our calculations we do not consider changes in the orbital parameters of the colliding bodies (i.e., SP), since this effect is not important for the present study. There are 2 reasons for this: (i) the lifetime of an avalanche (typically $\sim$10³ years) is very short from the point of view of the global disc evolution, thus we can neglect any changes in the disc dynamics caused by mutual collisions between the disc particles (i.e., "field SPs'' in our simulations); (ii) the dynamics for the majority of the avalanche SPs are controlled by the radiation pressure. Their orbital parameters are thus determined mostly by their $\beta$ values and only weakly depend on the velocities at which these SPs are born (as is verified in Sect. 4.2.1 for the first generation of avalanche grains). In this respect, taking the velocity of the center of mass of the colliding grains as the initial velocity for the produced debris is a good approximation within the frame of our simulations.

   
3.5 Initial planetesimal breakup

As previously mentioned, we assume that the initial source of the collisional avalanche is the breakup of a large, at least kilometer-sized object. We do not perform a simulation of the initial shattering event itself, but implement a simple parametric prescription for the dust released in the breakup. In most runs, we consider a "nominal'' case, in which M₀=10²⁰ g of dust is released in the $0.1~\mu$m to 1 cm range at $R_{\rm0}=20$ AU from the star, unless otherwise explicitly specified. It should be noted that the released dust mass M₀ is the only relevant parameter for our simulations. In this respect, the exact process leading to the initial release is not crucial. However, when it comes to estimating the probability for such a dust-release event to occur (as will be done in Sect. 6.1), one has to consider the mass  $M_{\rm PB}$ of the parent body whose shattering produces a mass M₀ of dust. The ratio  $M_0/M_{\rm PB}$ is obviously <1, but strongly depends on several poorly constrained parameters, mainly related to the physics of the shattering event. For an idealized case when the largest fragment produced has mass $M_{\rm lf}=0.5M_{\rm PB}$ and smaller fragments follow the Dohnanyi "equilibrium'' size distribution ( ${\rm d}n \propto s^{-3.5}{\rm d}s$), one gets $M_{\rm PB}\simeq 10^{24}$ g. However, laboratory and numerical studies as well as observations of asteroid families (e.g., Tanga et al. 1999; Davis & Ryan 1990) all point towards smaller  $M_{\rm lf}$ and steeper size distributions for highly disruptive impacts of large objects, with indexes typically in the -3.7 to -4 range for the largest $\geq$ $0.01M_{\rm lf}$ fragments and closer to -3.5for the smallest ones. Using for example the fragmentation prescription for large objects of Thébault et al. (2003), we determine that for a typical shattering at 1 km s^-1, $M_{\rm PB}\simeq 10^{21}$ g = 10M₀, which corresponds to an object of radius $\sim$40 km. We shall thus assume a nominal $M_0/M_{\rm PB}$ ratio of 0.1 for the discussion in Sect. 6.1. For the size spectrum of the dust particles released in the $0.1~\mu$m to 1 cm range, we assume a single power law (Eq. (8) with p_s0=3.5) for our nominal case. The dependence of the results on M₀, p_s0, and other parameters related to the planetesimal debris is explored in Sect. 4.2.

   
4 Results

For the sake of the readability of the results, it is convenient to divide the system into two populations: 1) the avalanche particles, representing all bodies initially released by the planetesimal breakup plus all grains later created by collisions between the avalanche particles and the disc material; and 2) the field particles, i.e., the population of grains in the disc unaffected by the avalanche mechanism. To quantify the magnitude of an avalanche we introduce the area amplification factor, $F_{\rm }$, which is the ratio of the total cross-sectional area of the avalanche grains, within 500 AU from the star, to the initial cross-sectional area of planetesimal debris released. The maximum value  $F_{\rm max}$ reached by the amplification factor while the avalanche is propagating is used to measure the amplitude of a given avalanche and to compare avalanches obtained for different initial conditions. Time is expressed in orbital periods at 20 AU ($\sim$70 yr), unless otherwise explicitly specified. Table 1 summarizes the set of initial parameters chosen for our "nominal'' case. All free parameters of the simulations are then explored in separate runs.

4.1 Nominal case (NC)

Table 1: Main model parameters for the nominal case.

Figure 2: Nominal case. Color-coded maps (log-scale) of the vertical optical thickness of avalanche grains, $\tau _{\rm\perp ,{\rm av}}$, at different stages of the avalanche evolution (t= 0.6, 5, 10, 40 orbital periods at 20 AU). The planetesimal debris are released at t=0 at 20 AU from the star. Field particles are not included in the plots. The position of the star is marked by the white cross.

Open with DEXTER

Figure 2 shows the temporal evolution of an avalanche, for the nominal case, in terms of the vertical optical thickness, $\tau _{\rm\perp ,{\rm av}}$, of the avalanche particles ( $\tau_{\rm\perp, av}= \int\!\!\int \pi s^2 {\rm d}n_{\rm av}(s) {\rm d}z$). As expected, the first stages correspond to a fast development and multiplication of the avalanche grains. In this early expansion phase the surface density is dominated by the smallest high-$\beta$ ($\ga$0.5) particles, which contribute to $\simeq$$85\%$ of  $\tau _{\rm\perp ,{\rm av}}$. The maximum value of the amplification factor is $F_{\rm max}=210$ and is reached at $t \simeq 5$ ($\simeq$350 yr, see Fig. 3). After that, the loss of small grains on unbound orbits dominates over the collisional production of new dust particles, and the avalanche begins to fade. In these later stages, the total cross sectional area of the avalanche grains (within 500 AU) is increasingly dominated by the bigger grains on bound orbits. It is important to point out that the timescale for the avalanche propagation is short in comparison with orbital periods in the outer part of the disc (e.g., only $\sim$1/5 of the orbital period at 200 AU).

The amplification achieved by the avalanche mechanism is impressive, i.e., an increase in grain cross-sectional surface density by two orders of magnitude compared to the particles initially released by the planetesimal breakup (Fig. 3). However, absolute values of  $\tau _{\rm\perp ,{\rm av}}$ are still very small compared to those of the field particles, with $\tau _{\rm\perp ,{\rm av}}$/ $\tau_{\rm\perp,{\rm field}}$ never exceeding 10^-2(see Sect. 5 for a more detailed discussion of this crucial parameter).

To compare the results of our simulation with the simplified theory of Sect. 2, we plot the ratio of the total number of grains  $N_{\rm tot}$ produced by the avalanche until time t to the initial number $N_{\rm0}$ of released planetesimal debris (Fig. 4). As can be clearly seen, $N_{\rm tot}/N_{\rm0}$ quickly reaches a plateau, and we take $N_{\rm tot}/N_{\rm0} \simeq 200$ at t=15 as a reference value. Plugging values for the average number of particles produced by each grain-grain collision, $\langle N_{\rm\beta} \rangle \approx 150$, and $\tau_{\rm\Vert}=0.022$ into Eq. (5), we get $N_{\rm tot, theory}/N_{\rm0} \approx 30$, i.e., a factor of $\sim$7 difference with the result of our simulation. This is mainly due to the fact that  $\tau_{\rm\Vert}$ underestimates the real value of $\tau$, firstly because the real path of a grain is curved rather than parallel to a disc radius, and secondly because in $\tau_{\rm\Vert}$ the size of the avalanche grains $s_{\rm gr}$ is neglected. From our simulations, we were able to estimate the discrepancy between $\tau$ and $\tau_{\rm\Vert}$ to be roughly of a factor of 1.6. We thus get

$$% \frac{N_{\rm tot}}{N_{\rm0}}\approx \mbox{e}^{1.6\tau_{\rm\Vert}\langle N_{\rm\beta}\rangle} \simeq200,$$ (13)

a value close to the numerical results.

Figure 3: Time evolution of the cross-sectional area amplification factor (the ratio of the total cross-sectional area of the avalanche grains within 500 AU to its initial value at t=0). Initial increase is due to dust production by outflowing planetesimal debris colliding with the disc material. When the grain removal (due to star radiation pressure) rate exceeds the grain production, the value of F drops (see text for more details).

Open with DEXTER

	Figure 4: Ratio of the total number of grains produced by the avalanche by the time t, $N_{\rm tot}$, to the initial number of released planetesimal debris, N0.
Open with DEXTER

It is also interesting to link $N_{\rm tot}/N_0$ to the amplification factor parameter  $F_{\rm max}$. By definition,

$$% F(t)=\frac{\int s^2 {\rm d}N_{\rm in}(s,t)}{\int s^2 {\rm d... ..._{\rm in}(t) \langle s(t)^2\rangle}{N_0 \langle s_0^2\rangle},$$ (14)

where $N_{\rm in}(s,t)$ is the total number of avalanche grains inside 500 AU at time t, and $\langle s(t)^2\rangle$ and  $\langle s_0^2\rangle$ are the averaged cross-sectional areas of the avalanche grains at time t and of the initially released planetesimal debris, respectively. Since the avalanche's dust production comes mainly from a rather short peak of activity (see Figs. 3 and 4), $N_{\rm tot}$ is close to  $N_{\rm in}(t_*)$, where t_* is the time at which  $N_{\rm in}$ reaches its maximum. A numerical check showed that indeed $N_{\rm in}(t_*) \approx 2/3 N_{\rm tot}$, so that

$$% F_{\rm max} \approx \frac{2}{3}~\frac{N_{\rm tot}}{N_0} \fr... ... } ~ \mbox{e}^{1.6 \langle N_{\rm\beta} \rangle \tau_{\Vert}}.$$ (15)

The validity of this relation is easily numerically verified, with $2/3 N_{\rm tot}/N_0 \times \langle s(t_*)^2 \rangle/ \langle s_0^2 \rangle \simeq 200$, a value that is indeed relatively close to the $F_{\rm max}$ value obtained in the simulation.

   
4.2 Dependence of ${F_{\max}}$ on the initial planetesimal debris parameters

   
4.2.1 Initial mass and velocity of the planetesimal dust cloud

The initial mass M₀ of dust released has been explored as a free parameter. The simulations show that the maximum amplification factor, $F_{\rm max}$, does not vary with M₀, at least in the 10¹² g-10²¹ g range, a result that is in good agreement with Eqs. (15) and (3). Likewise, $F_{\rm max}$ does not change much when varying the initial speed $v_{\rm0}$ of the center of mass of the planetesimal dust cloud. There is only a 20% increase of  $F_{\rm max}$ when $v_{\rm0}$ is increased from  $v_{\rm kep}$ to  $1.41v_{\rm kep}$. This weak dependence on the initial velocity of the debris confirms the fact that avalanches are driven mostly by the smallest particles, which are accelerated to high speeds weakly correlated to the initial release velocity.

   
4.2.2 Size distribution of planetesimal debris

For the nominal case we choose $s_{\rm min,pl}=0.1~\mu$m. This value is compatible with the lower limit for the size of the interplanetary dust particles (Fraundorf et al. 1982). It is also in good agreement with the size distribution deduced from studies of cometary comas that show that the smallest particles are about $0.08{-}0.28~\mu$m in diameter (Kolokolova et al. 2001). McDonnell et al. (1991) observed smaller grains in comas, but their contribution to the total dust population remained marginal. Even if grains smaller than $0.1~\mu$m are produced abundantly in the planetesimal breakup, they are not expected to contribute significantly to the avalanche process since they are in the size range where $\beta$ decreases for smaller grains (see Fig. 1). As a consequence, they have lower outgoing velocities which, together with their smaller masses, lead to a marginal contribution in terms of impacting kinetic energy. We thus believe $0.1~\mu$m to be a reliable minimum value for  $s_{\rm min,pl}$and explored $s_{\rm min,pl}$ values in the $0.1~\mu$m to $1~\mu$m range, the latter value being the one considered by Kenyon & Bromley (2005). Although the impacting kinetic energy per grain is increasing in the $0.1~\mu$m to $1~\mu$m size range (leading to an increase of the $N_{\rm tot}/N_0$ ratio), $F_{\rm max} \propto N_{\rm tot}/N_0 \times \langle s(t_*)^2 \rangle / \langle s_0^2 \rangle$ decreases with increasing $s_{\rm min,pl}$ (Fig. 5) because of the decreasing value of the $\langle s(t_*)^2 \rangle / \langle s_0^2 \rangle$ factor.

	Figure 5: Maximum amplification factor as a function of the minimum size assumed for the initial planetesimal debris. The power-law index for the size distribution is equal to its nominal value $p_{\rm0,pl}=3.5$. For smaller $p_{\rm0,pl}$ the variation with  $s_{\rm min,c}$ is less pronounced.
Open with DEXTER

	Figure 6: Dependence of the maximum amplification factor on the power-law index of the initial size-frequency distribution of the planetesimal debris (Eq. (8)), $s_{\rm min,pl}=0.1~\mu$m.
Open with DEXTER

Test runs have also been performed to check the $s_{\rm max,pl}$ dependence. This exploration has shown that results do not depend on this parameter for values higher than 1 cm. This is due to the fact that grains bigger than this size have very small $\beta$ values, as well as a low total cross-sectional area, which do not allow them to significantly contribute to the avalanche propagation.

The dependence of $F_{\rm max}$ on the power-law index for the initial planetesimal debris size distribution, $p_{\rm0,pl}$, is shown in Fig. 6. It can be noted that increasing $p_{\rm0,pl}$ above 3.5 (value for the nominal case) does not lead to a significant increase of  $F_{\rm max}$, whereas less steep power laws lead to a significant decrease in  $F_{\rm max}$.

   
4.2.3 Position of the planetesimal breakup

We perform a set of runs in which the position of the planetesimal breakup, R₀, is varied between 20 and 100 AU (Fig. 7), but all the other parameters remain identical to the nominal case. The maximum amplification factor decreases with increasing R₀ for two reasons: (i) the total amount of disc material through which the outflowing grains propagate is higher when the grains are released close to the star; (ii) the unbound grains ( $\beta > 0.5$) have time to reach higher radial velocities if they are released closer to the star (see Fig. 8), which leads to more violent collisions and hence higher dust production per collision.

	Figure 7: Dependence of the maximum amplification factor on the location of the primary planetesimal breakup.
Open with DEXTER

	Figure 8: Radial velocities vs. distance from the star for grains with different $\beta$ values released by parent bodies on circular orbits for a $\beta$ Pic-like star. The release distances are 20 AU and 70 AU.
Open with DEXTER

4.3 Dependence of ${F_{\max}}$ on the prescription for collisional outcome

  
4.3.1 Minimum size of avalanche produced debris

In the nominal case we assume that the minimum size, $s_{\rm min,col}$, for the debris produced by collisions and the minimum size of the initial planetesimal debris, $s_{\rm min,pl}$, are both equal to $0.1~\mu$m. We have seen in Sect. 4.2.2 that we do not expect significant changes when $s_{\rm min,pl}<0.1~\mu$m. However, the situation is slightly different for avalanche grains, since these grains are continuously produced through collisions. We investigate this parameter's effect in test runs exploring different values for $s_{\rm min,col}$ (Fig. 9). As can easily be seen, $F_{\rm max}$ does not strongly vary with $s_{\rm min,col}$ for $s_{\rm min,col}<0.1~\mu$m. There are two reasons for that: i) $\beta$ decreases with decreasing sizes for grains smaller then $\sim$$0.1~\mu$m (Fig. 1), thus preventing them from significantly contributing to the avalanche propagation; ii) the broken power law for the debris size distribution, and especially the flatter index q₁=1.5 for the smallest grains, prevents them from taking up most of the cross-sectional area of the avalanche grains.

	Figure 9: Dependence of the maximum amplification factor on the minimum size for debris produced in avalanche collisions (Eq. (12) with $q_1=1.5, q_2=1.83, m_s=M_{\rm lf}/3$).
Open with DEXTER

4.3.2 Size distribution of the debris grains

Numerical exploration of the position of the slope change m_s (with $s_{\rm min}=0.1~\mu$m) and of the m>m_s power-law index q₂ shows that the resulting amplification factor only weakly depends on these parameters. Changing $m_s/M_{\rm lf}$ from 1 to 10 leads to changes in $F_{\rm max}$ by only a factor $\sim$2. On the contrary, variation in the q₁ index can significantly affect  $F_{\rm max}$, especially for q₁>5/3, when most of the produced cross-sectional area resides in the smaller grains (Fig. 10).

	Figure 10: Dependence of the maximum amplification factor on the value q1 in Eq. (12) for collisionally produced grains (q2=1.83).
Open with DEXTER

   
4.4 Grain chemical composition and impact strength

The exact chemical composition of circumstellar disc material is not well constrained and might in any case vary from one system to the other. There is observational evidence for silicates, ices, and metals (e.g., Pantin et al. 1997; Bouwman et al. 2003), but their exact proportions in individual grains are difficult to estimate. Several detailed studies have addressed this issue for the specific $\beta$ Pic case (e.g., Pantin et al. 1997; Li & Greenberg 1998), but the estimates remain model dependent. Changes in grain compositions might affect the results in two ways: (i) compositional changes can lead to different values of $\beta$, so that the grains experience different radiation pressure and as a consequence reach different outgoing (and impacting) velocities; and (ii) their collisional response properties can be significantly different.

We first explore the role of grain porosities by varying this parameter between P=0 (compact grains, nominal case) and P=0.8 (highly porous grains), with Q^*₀ remaining constant. This constant Q^*₀ prescription might seem counter-intuitive at first, since more porous grains should be expected to be more fragile, but it is, in fact, supported by numerical experiments showing that porous targets often prove more resistant than non-porous ones (e.g., Ryan et al. 1991; Flynn & Durda 2004; Love et al. 1993), the reason being that impact shock waves are effectively dissipated by the pores. Figure 11 shows that avalanche strength is maximum for the nominal case of compact grains ( $F_{{\rm max}(P=0)}\simeq210$) and decreases for porous grains ( $F_{{\rm max}(P=0.8)}\simeq70$).

We numerically explore the importance of chemical composition by performing runs for the 2 extreme cases of pure (compact) silicates (Mg $_{\rm0.95}$Fe $_{\rm0.05}$SiO₃) and pure (compact) water ices. Here again, we take the possibly counter-intuitive constant Q^*₀ assumption, which is here again supported by experimental results showing that for target-projectile pairs of the same material, ices can be as strong as silicates (e.g., Ryan et al. 1999). Furthermore, compact ices and silicates of equivalent sizes have similar $\beta$ values in the $s>0.1~\mu$m range (see Fig. 1). It is thus not surprising that our results show no significant difference between the pure-ice ( $F_{\rm max}=300$) and pure-silicate runs ( $F_{\rm max}=210$).

In a third set of runs we separately explore the Q^* parameter, whose values for given grain compositions and dynamical conditions are still not well constrained by experiments or theoretical studies. The threshold energy estimates for silicate-silicate and ice-ice collisions might vary between $\simeq$10⁶ and a few $\times 10^7$ erg/g (e.g., Ryan et al. 1999; Holsapple et al. 2002; Benz & Asphaug 1999). We explore Q₀^* values between 10⁶ and 10⁸ erg/g and obtain strong variations in  $F_{\rm max}$ (Fig. 12). For the lowest explored Q₀^* value of 10⁶ erg/g, we get $F_{\rm max}=10^4$, which is about 50 times higher than in the nominal case ( Q₀^* = 10⁷ erg/g).

	Figure 11: Maximum amplification factor as a function of porosity for pure silicate grains. The value of the threshold energy, Q*, is assumed to be the same as in the nominal case.
Open with DEXTER

	Figure 12: Maximum amplification factor vs. value of the threshold energy, Q*. Values for s0=1 cm grains are denoted on the axis. For the other sizes the threshold energy is given by Eq. (11).
Open with DEXTER

   
4.5 Field particle population

	Figure 13: Different test radial distributions for $\tau _{\rm\perp }$ ($\tau_{\Vert}$ being constant). The thick solid line is the distribution for the nominal case, taken from Augereau et al. (2001).
Open with DEXTER

As mentioned earlier, our reference field particle disc was assumed to be similar to the $\beta$ Pictoris system, for which the dust profile derived by Augereau et al. (2001) has been taken. Here we explore alternative profiles (Fig. 13). Results show that  $F_{\rm max}$ does not strongly depend on the shape of the density distribution profile as long as the total radial optical depth of the system (within 500 AU) $\tau_{\rm\Vert}$remains the same. This result is in agreement with the simplified theory presented in Sect. 2.

On the other hand, we get drastic $F_{\rm max}$ variations when changing the value of $\tau_{\rm\Vert}$ (regardless of the radial profile). Figure 14 shows for example that increasing the number density by a factor of 5 leads to a value of  $F_{\rm max}$, which is a factor of $\sim$1000 higher. This strong increase in  $F_{\rm max}$ is in agreement with Eq. (15), which predicts a strictly exponential growth with  $\tau_{\rm\Vert}$, if $\langle N_{\rm\beta} \rangle$ is constant. However, in the simulations we find that $\langle N_{\rm\beta} \rangle$ weakly varies with  $\tau_{\rm\Vert}$ through the empirical relation $\langle N_{\rm\beta}\rangle \approx 150~ (\tau_{\rm\Vert}/\tau_{\rm\Vert,nom})^{-0.45}$. Plugging this expression for $\langle N_{\rm\beta} \rangle$ and ${\langle s(t_*)^2\rangle }/{\langle s_0^2\rangle }\simeq 1.5$ into Eq. (15) we get:

$$% \begin{array}{ll} F_{\rm max} & \simeq \exp\left[240~ \tau... ...}}{\tau_{\rm\Vert,nom}} \right)^{0.55}\right]\cdot \end{array}$$ (16)

Thus $\tau_{\rm\Vert}$ proves to be the most efficient parameter for increasing $F_{\rm max}$ by several orders of magnitude. This is a point of crucial importance when considering the absolute strength of an avalanche and its possible observability.

	Figure 14: Maximum amplification factor as a function of $\tau_{\Vert}$. The open circles are the results of our simulations. The dashed line is the theoretical prediction (Eq. (16)).
Open with DEXTER

   
5 Avalanche observability

So far we have been concerned with the way an avalanche develops in a disc, in particular with how much dust can be created compared to the initial amount of released grains. This was quantified by the amplification factor  $F_{\rm max}$. However, the crucial issue is under which conditions such an avalanche might become observable. In this respect, looking at  $F_{\rm max}$ is not enough. What matters here is the ratio between the luminosity of the avalanche particles and that of the "field'' particles, $L_{\rm av}/L_{\rm d}$.

The value of $L_{\rm av}/L_{\rm d}$ corresponding to the observability limit depends on several factors, such as the physical parameters of the system, observational conditions, and the observing devices' characteristics. Since our current study is not dedicated to a specific system, it is impossible to give a precise criterion for avalanche observability. We shall thus adopt a simple and probably conservative criteria in which an avalanche is deemed observable when $L_{\rm av}/L_{\rm d} \ga 1$ is reached at a given location in the disc.

Most of the resolved debris disc images have been obtained in the visual or near-infrared (NIR) domains, dominated by scattered starlight. We shall thus focus here more specifically on scattered light luminosity. The amount of light scattered towards an observer coming from a given region of the disc is proportional to

$$% L{\rm }(\lambda_0) \propto \int_{V}\!\!\int_{s} F_{*}(\lam... ...\rm sca}(s,\lambda_0) f(\theta) {\rm d}n(s,\vec{r}) {\rm d}V,$$ (17)

where the integration is done over the spatial volume of the considered region and over the whole grain size range, s, with the grain number density n, the scattering coefficient  $Q_{\rm sca}$, and the scattering function $f(\theta)$. $F_{*}(\lambda_0,\vec{r})$ is the monochromatic star flux at the distance $\vec{r}$ from the star.

Figure 15: Top panel: face-on case: color-coded map of the ratio between the geometric surface densities of the avalanche grains and that of the "field'' population for the nominal case. Bottom panels: edge-on case: midplane fluxes (arbitrary units) for the nominal case at avalanche maximum, edge-on orientation. The 2 solid lines indicates the total midplane fluxes ("field''+avalanche) for each side of the disc (differences between the 2 sides are so small here that the 2 lines are almost indistinguishable). The dashed lines show the midplane fluxes for just the avalanche particles. Plot  a) corresponds to the forward scattering function and  b) to the isotropic case.

Open with DEXTER

For avalanche detection, the visible domain ($\sim$$0.5~\mu$m) is probably more favorable than the NIR ($1{-}2~\mu$m). This is because avalanches consist mostly of submicron grains, which scatter very inefficiently at 1-$2~\mu$m compared to bigger grains a few microns in size (which is the average size for the "field'' population). At the same time, in the visual domain  $Q_{\rm sca}$ is nearly the same (within a factor of 2, depending on the exact chemical composition) for submicron and micron grains. Thus the ratio $L_{\rm av}/L_{\rm d}$ is expected to be higher in the visual than in the NIR. For simplicity we assume that $f(\theta)$ is only a function of the scattering angle $\theta$ and that  $Q_{\rm sca}$ is independent of the grain size. Although it is not exactly the case, this simplification can be considered as a reasonable starting point.

Here we consider two extreme cases of disc orientation, i.e., discs seen exactly edge-on and exactly pole-on. For the pole-on case, we determine that $L_{\rm av}/L_{\rm d} \simeq$ $\tau _{\rm\perp ,{\rm av}}$/ $\tau_{\rm\perp,{\rm field}}$ and thus use maps of the ratio between these vertical optical depths. For the edge-on case, we consider the synthetic midplane flux in scattered light, computed for different scattering functions. Since we do not know the exact optical properties of the circumstellar grains, two bracket cases have been considered for our calculations: isotropic and forward scattering. For the forward scattering function we use an analytical approximation of the empirical  $f_{\rm scat}$ for cometary dust, obtained from measurements of solar system comets (Artymowicz 1997, and references therein):

$$% f_{\rm scat}(\theta) =f_0 \left[ \frac{0.3}{(0.2+\theta/2)^3} +1.4 \left(\frac{\theta}{3.3}\right)^4 +0.2\right].$$ (18)

In the following subsections, we investigate under which conditions avalanche-induced asymmetries might become observable for these two disc viewing angles. As appears from Figs. 15-18, these asymmetries consist of partial spiral or lumpy patterns in the face-on case, and of two-sided asymmetries, for which one side of the system becomes brighter than the other, in the edge-on case. We would like to point out that at other wavelengths, e.g., infrared, the observability criteria ( $L_{\rm av}/L_{\rm d} \ga 1$) might be reached for lower $\tau _{\rm\perp ,{\rm av}}$/ $\tau_{\rm\perp,{\rm field}}$ ratios due to the fact that avalanche grains are expected to be hotter than field particles, since their average size is about 10 times smaller then the average size of the initial disc population.

   
5.1 Nominal case

As can be clearly seen in Fig. 15, in the nominal case  $L_{\rm av}/L_{\rm d}$ never exceeds 10^-2, neither in the edge-on nor in the head-on configuration. This value is far below our observability criterion and the asymmetries induced by the corresponding avalanche would thus probably be undetectable by scattered light observations.

   
5.2 Larger amount of released dust M₀

The most straightforward way of getting a more prominent avalanche is to increase the initially released amount of dust. As shown in Sect. 4.2.1, $F_{\rm max}$ remains constant with varying M₀, so that the ratio  $L_{\rm av}/L_{\rm d}$ increases linearly with M₀. As a consequence, the release of $\approx$10²² g of dust would be required for the avalanche-induced luminosities to become comparable to that of the rest of the disc. One might wonder however if a planetesimal shattering releasing this large amount of dust is a common event (see discussion in Sect. 6.1).

5.3 Collisionally weaker grains

	Figure 16: Same as Fig. 15, but for collisionally weaker grains with $Q^*_0(s_0=1~{\rm cm})=10^6$ erg/g (see Sect. 4.4).
Open with DEXTER

As has been seen in Sect. 4.4, $F_{\rm max}$ increases strongly for grains with lower specific energy values Q_*. The lowest Q_* value explored in Fig. 11, $Q^*_{\rm0(s_0=1~{\rm cm})}=10^6$ erg/g, leads to $\tau _{\rm\perp ,{\rm av}}$/ $\tau_{\rm\perp,{\rm field}}$ $\simeq$ 0.4-0.5. Thus, observability might be marginally reached when assuming the minimum shattering resistance for dust grains.

   
5.4 Dustier systems

	Figure 17: Same as Fig. 15, but for the case of a disc 5 times more massive than in the nominal case.
Open with DEXTER

The parameter exploration of Sect. 4 has clearly shown that the most efficient parameter for reaching high $F_{\rm max}$ values is the field particles number density $\tau_{\Vert}$ (Eq. (16)). An obvious way of increasing $\tau_{\Vert}$ is to assume a more massive disc, as has been done in Sect. 4.5. In terms of avalanche observability, we find that the observability criteria, $L_{\rm av}/L_{\rm d} \simeq 1$, is reached for a disc that is 4-5 times more dusty than in the nominal case. In this case, azimuthal asymmetries become clearly visible in the face-on configuration and two-sided brightness asymmetries for the edge-on case (Fig. 17). A massive dusty disc thus looks very promising from the point of view of avalanche observation.

	Figure 18: Same as Fig. 15, but for a vertically thinner disc, described by Eq. (19). The total mass of the disc is the same as in the nominal case, but the asymmetries become prominent.
Open with DEXTER

Another way to reach higher values of $\tau_{\Vert}$ is to keep the same total amount of dust, but distributed in a vertically thinner disc. In the nominal case the dusty disc has the characteristic width $w_{\rm }$ (Eq. (10)), a superexponential vertical profile (Eq. (9)), and a corresponding $F_{\rm max}=210$. Assuming now a disc of thickness $w^*_{\rm }=0.25w_{\rm nom}$, with a constant vertical profile

 

$$% \frac{{\rm d}\tau}{{\rm d}z} = \left\{ \begin{array}{ll} \tau(R... ...$ } \\ 0, &\mbox{if $\vert z\vert > w^*_{\rm disc}$ ,} \\ \end{array} \right.$$ (19)

we get $F_{\rm max}=2 \times 10^4$ for the same total amount of dust as in the nominal case. In this case, the number density of avalanche grains can even largely exceed that of the field particles (see Fig. 18). This density enhancement is azimuthally asymmetric due to the spiral structure of an avalanche. If the disc is orientated face-on, then the azimuthal asymmetry persists for about 800 years. If the system is viewed edge-on then a two-sided asymmetry can be observed. Figure 18 displays midplane fluxes ("field''+avalanche) for different azimuthal angles and scattering functions, clearly showing that, for favorable system orientations, one side gets significantly brighter because of the avalanche.

   
6 Discussions

   
6.1 Probability of witnessing an avalanche event

The numerical investigation of the previous sections has shown that collisional avalanches are a powerful and efficient mechanism that naturally develops in debris discs after the breakup of a large planetesimal. However, in our nominal case of a $\beta$ Pic-like system and M₀=10²⁰ g of dust initially released, the asymmetric features produced by the avalanche probably remain too weak to be observable in scattered light (Sect. 5.1). This result should, however, be taken with great care since our parameter exploration has shown that avalanche strength strongly depends on several critical and often poorly constrained parameters. The first set of parameters is those linked to the initial breakup event. Here we obtain the intuitive result that higher amounts of initially released dust leads to more powerful avalanches (see Sect. 5.2), with the avalanche strength scaling linearly with M₀. This is not unique to the avalanche mechanism: Kenyon & Bromley (2005) find a similar dependence when only considering the signature of the cloud of primary debris produced immediately after the planetesimal breakup. What distinguishes our results from studies in which only dust released at impact is considered is that avalanches strongly depend on the number density of dust in the disc. Section 4.5 has indeed shown that the global optical depth of the dust disc $\tau_{\Vert}$ is the parameter avalanche development depends most on, the dependence being close to an exponential. We have seen that other parameters, mostly related to the way the physical response of grains to collisions is modeled, might also lead to observable events when stretched to the extreme values that were numerically explored here. This is in particular the case for Q_*, for which very low $\simeq$10⁶ erg/g values might lead to powerful avalanches.

We shall however leave these "technical'' parameters aside to focus on the 2 parameters directly related to the system's properties themselves, i.e., the optical depth, both  $\tau_{\bot}$ and $\tau_{\Vert}$, and the initial amount of dust released M₀, and derive an order-of-magnitude estimate for the probability of witnessing avalanche events as a function of these parameters. From the results of Sect. 4, the $L_{\rm av}/L_{\rm d}$ criterion for observability might be written

$$% \left( \frac{\tau_{\bot}}{\tau_{\rm\bot,nom}} \right)^{-1} ... ... max}}{F_{\rm max(nom)}}~ \frac{M_0}{10^{20}~{\rm g}} \ga 100,$$ (20)

which is equivalent to saying that the luminosity ratio between avalanche and field grains should be at least 100 times higher than in the nominal case (for which $L_{\rm av}/L_{\rm d} \sim 10^{-2}$). Section 4.2.1 has shown that  $F_{\rm max}$ is independent of M₀, so that in our approximation  $F_{\rm max}$ is only a function of $\tau_{\Vert}$, and this $\tau_{\Vert}$ dependence is given by Eq. (16). Thus, Eq. (20) reduces to

$$% \exp \left[ 5.3 \left( \frac{\tau_{\Vert}}{\tau_{\rm\Vert,n... ...rac{\tau_{\bot,{\rm nom}}}{\tau_{\rm\bot}}~ \ga 2\times10^{4},$$ (21)

which gives a direct link between a given disc density ( $\tau_{\bot}$ and $\tau_{\Vert}$) and the minimum mass of released dust able to produce a visible avalanche in such a disc (the denser the disc, the smaller the corresponding M₀ value). The other important issue affecting witnessing probabilities is the duration of an avalanche. Our simulations show that the typical lifetime of an avalanche-induced pattern is $t_{\rm av} \sim 10^3$ yr. With this value and Eq. (21), one can estimate the probability  $P_{\rm obs}$of witnessing an observable avalanche event in a given disc:

$$% P_{\rm obs} = \frac{t_{\rm av}}{t_{\rm shatt(M_0,\tau_{\bot})}},$$ (22)

where $t_{\rm shatt(M_0,\tau)}$ is the average time between 2 shatterings producing M₀ of dust in a disc of average optical depth $\tau_{\bot}$, with M₀, $\tau_{\Vert}$ and $\tau_{\bot}$ satisfying Eq. (21). As suggested in Sect. 3.5, we consider that the object releasing M₀ of dust has a mass $M_{\rm PB}\simeq 10M_0$. From unpublished results of the Thébault et al. (2003) simulations of collisional rates and outcomes in the inner $\beta$ Pic disc, we determine that the approximate timescale for the shattering of a  $M_{\rm PB}=10M_0$ object to occur in the innermost <50 AU (the typical location for the initial shattering events considered in our simulations) of a $\beta$ Pic like system is $t_{\rm shatt}\simeq 150[(10M_0)/10^{21}$ g]^1.25 yr. Since, for systems with similar spatial distributions, the frequency of collisional events is proportional to the square of a system's total mass, we get the empirical relation:

$$% t_{\rm shatt(M_0,\tau_\bot)} \simeq ~150~ \left( \frac{\tau... ...eft( \frac{M_0}{10^{20}~ \mbox{g}} \right)^{1.25} ~\mbox{yrs,}$$ (23)

where we implicitly assume that the system's spatial distribution is the same as in the nominal case, so that the ratio between two systems' total masses is equal to the ratio $\tau_\bot / \tau_{\rm\bot,nom}$ anywhere in the disc. This equation should of course be regarded as giving a very rough estimate, since $t_{\rm shatt(M_0,\tau_\bot)}$ depends on many poorly constrained parameters, like the number density of planetesimals and their average kinetic energy at impact. Equation (23), however, gives the global trend of the way $t_{\rm shatt(M_0,\tau_\bot)}$ increases with M₀. Taking the lowest M₀ value satisfying Eq. (21) and plugging it into Eq. (23), we get, from Eq. (22):

$$% P_{\rm obs} \approx 3 \times 10^{-5} \left( \frac{\tau_\bot... ...ac{\tau_\Vert}{\tau_{\rm\Vert,nom}} \right)^{0.55}\right)\cdot$$ (24)

Equation (24) indicates that $P_{\rm obs}\simeq 0.03$ for the nominal case field particle disc, which means that we have about a $3\%$ chance of witnessing the avalanche caused by the breakup of a  $M_{\rm PB}=10M_0 \sim 10^{23}$ g object ( M₀=10²² g being the smallest released dust mass able to trigger a visible avalanche for such a disc, as given by Eq. (21)). This makes it a rather unlikely event, although it cannot be completely ruled out. Nevertheless, slightly denser discs (i.e., higher $\tau_{\bot}$ and $\tau_{\Vert}$) can easily raise $P_{\rm obs}$ up to 1. As a matter of fact, the dependence on $\tau_{\Vert}$is so sharp that $P_{\rm obs}=1$ is obtained for $\tau_{\Vert}\simeq 2.1\tau_{\Vert,{\rm nom}} \simeq 0.046$. We thus see that a $\beta$ Pic-like system is below, but not too far from the limit for which chances of witnessing an avalanche are high, especially when considering the uncertainties regarding avalanche strength due to its dependence on several poorly constrained parameters related to the collision-outcome prescription (also keeping in mind that higher $\tau_{\Vert}$ values could alternatively be achieved for a thinner disc of the same dust mass (see Sect. 5.4)). Moreover, our $L_{\rm av}/L_{\rm d}>1$ criterion for observability is probably too conservative, and avalanche-induced patterns might be detectable for lower luminosity excess values. Taking, for example, $L_{\rm av}/L_{\rm d}>0.1$ would raise the detection probability to $\simeq$$45\%$ for a $\beta$ Pic-like system.

6.2 Avalanches in observed systems, perspectives

We defer a detailed application of our model to specific circumstellar discs to a future study. However, the present results can already give a good idea of the typical profile for a "good'' avalanche-system candidate.

Our numerical exploration has shown that structures that are the most likely to be associated with avalanche-events have two-sided asymmetry for discs viewed edge-on and open spiral patterns for discs viewed pole-on or at intermediate inclinations. An additional requirement is that these discs should be dust-rich systems, with a dustiness at least equal to, and preferably higher than that of $\beta$-Pic. Note also that our model makes an additional prediction, i.e., that avalanche affected regions should consist of grains significantly smaller than the "field'' particles in the rest of the disc. If the blow-out radius of grains is of the order of the wavelength of the observed light, then this should translate into color differences between avalanche (bluer) and non-avalanche (redder) regions.

In this respect, one good edge-on candidate might be the recently discovered HD 32297 system, which exhibits a strong two-sided asymmetry. As reported by Schneider et al. (2005) and Kalas (2005), this system is a $\beta$ Pic analog with its SW side significantly brighter than the NE one within $\simeq$100 AU (Schneider et al. 2005) and possibly outside 500 AU (Kalas 2005). Such a two-sided asymmetry would be compatible with the ones obtained in our simulations (as shown for example in the bottom panel of Fig. 18). Furthermore, Kalas (2005) also reported a color asymmetry between the two sides, with the brighter one (SW) being significantly bluer. This seems to indicate that this side is made of smaller, possibly submicron grains (Kalas 2005). As previously discussed, this is what should be expected for an avalanche-affected region. However, an alternative scenario, like the collision with a clump of interstellar medium proposed by (Kalas 2005), might also explain the HD 32287 disc structure. Future imaging and spectroscopic observations are probably needed before reaching any definitive conclusions.

Among all head-on observed systems, the one displaying the most avalanche-like structure is without doubt HD 141569 (Clampin et al. 2003), with its pronounced spiral pattern. Furthermore, the disc's mass, significantly higher than $\beta$-Pictoris, makes it a perfect candidate in terms of witnessing probabilities. Of course, avalanche is not the only possible scenario here, and several alternative explanations, like an eccentric bound planet or stellar companion, or a stellar flyby have already been proposed (e.g., Augereau & Papaloizou 2004; Ardila et al. 2005; Wyatt 2005). One should, however, be aware that this system is strictly speaking not a "standard'' debris disc as defined by Lagrange et al. (2000) and as considered in the present study. Indeed, several studies seem to suggest the presence of large amounts of primordial gas (Zuckerman et al. 1995; Ardila et al. 2005).

Gas drag effects have been left out of the present study on purpose, mainly because, in the strict sense of the term, debris discs are systems where dust dynamics is not dominated by gas friction (Lagrange et al. 2000). Moreover, the correct description of dust-gas coupling adds several additional free parameters (gas density and temperature distributions, etc.) and requires a full 2D or 3D treatment of gas by far exceeding the scope of the present paper. However, the issue of avalanches in a gaseous medium might be a crucial one for those systems that are most favorable for avalanches, i.e., discs more dusty than $\beta$-Pic, a system which is already at the upper end of debris-discs in terms of dustiness (e.g., Spangler et al. 2001). Such more massive systems should fall into a loosely defined category of "transition'' discs between T-Tauri or Herbig Ae protoplanetary systems and "proper'' debris discs (see for example Sect. 4 of Dutrey et al. 2004). For such systems (of which HD 141569 is a typical example), which are younger than the more evolved debris discs, risks (or chances) of encountering large amounts of remaining gas are high. This crucial issue will be the subject of a forthcoming paper.

   
7 Summary

This paper presents the first quantitative study of the collisional avalanche process in debris discs, i.e., the chain reaction of dust grain collisions triggered by the initial breakup of a planetesimal-like body. We have developed a code that allows us to simultaneously follow both spatial and size distributions of the dust grain population, which is collisionally evolving because of impacts caused by small particles blown out of the system by the star's radiation pressure. Our results can be summarized as follows:

1.

Collisional avalanches propagate outwards leaving a characteristic spiral-shaped pattern in the system. Depending on the system's orientation, these patterns might appear as open spirals or lumpy structures (face-on geometry) or a two-sided brightness asymmetry (edge-on case). In a $\beta$ Pic-like disc, an avalanche lasts for about 10³ years.

2.

The strength of an avalanche depends linearly on the mass of the initially shattered object, but nearly exponentially on the optical depth of the dust disc in which it propagates. The disc's dustiness is by far the most crucial parameter here, making dusty discs much more favorable cases for avalanche propagation than tenuous ones.

3.

We define a conservative criterion for avalanche observability, from which we infer a relation between a given disc density and the minimum mass of the object that has to be shattered to reach observability. When coupling this relation to estimates for catastrophic disruption probabilities among planetesimal-objects in debris discs, we are able to derive a first-order estimate for the probability of witnessing an observable avalanche event in a given debris disc. For our reference $\beta$ Pic-like system it is of a few percents, but probabilities rapidly increase for slightly denser systems.

4.

Modeling of dustier young transitional discs may require the inclusion of gas drag, which may change both the morphology and the strength of the avalanche.

References

• Ardila, D. R., Lubow, S. H., Golimowski, D. A., et al. 2005, ApJ, 627, 986 [NASA ADS] [CrossRef]
• Artymowicz, P. 1988, ApJ, 335, L79 [NASA ADS] [CrossRef] (In the text)
• Artymowicz, P. 1997, Annual Review of Earth and Planetary Sciences, 25, 175 [NASA ADS] [CrossRef] (In the text)
• Augereau, J. C., & Papaloizou, J. C. B. 2004, A&A, 414, 1153 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Augereau, J. C., Nelson, R. P., Lagrange, A. M., Papaloizou, J. C. B., & Mouillet, D. 2001, A&A, 370, 447 [NASA ADS] [CrossRef] [EDP Sciences]
• Benz, W., & Asphaug, E. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef]
• Bouwman, J., de Koter, A., Dominik, C., & Waters, L. B. F. M. 2003, A&A, 401, 577 [NASA ADS] [CrossRef] [EDP Sciences]
• Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1 [NASA ADS] [CrossRef] (In the text)
• Clampin, M., Krist, J. E., Ardila, D. R., et al. 2003, AJ, 126, 385 [NASA ADS] [CrossRef]
• Davis, D. R., & Ryan, E. V. 1990, Icarus, 83, 156 [NASA ADS] [CrossRef] (In the text)
• Dohnanyi, J. W. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] (In the text)
• Dutrey, A., Lecavelier Des Etangs, A., & Augereau, J.-C. 2004, Comets II, 81 (In the text)
• Flynn, G. J., & Durda, D. D. 2004, Planet. Space Sci., 52, 1129 [NASA ADS] [CrossRef]
• Fraundorf, P., Walker, R. M., & Brownlee, D. E. 1982, in Comet Discoveries, Statistics, and Observational Selection, ed. L. L. Wilkening, IAU Colloq., 61, 383 (In the text)
• Fujiwara, A., Kamimoto, G., & Tsukamoto, A. 1977, Icarus, 31, 277 [NASA ADS] [CrossRef]
• Fukagawa, M., Hayashi, M., Tamura, M., et al. 2004, ApJ, 605, L53 [NASA ADS] [CrossRef] (In the text)
• Heap, S. R., Lindler, D. J., Lanz, T. M., et al. 2000, ApJ, 539, 435 [NASA ADS] [CrossRef]
• Holsapple, K., Giblin, I., Housen, K., Nakamura, A., & Ryan, E. 2002, Asteroids III, 443
• Housen, K. R., & Holsapple, K. A. 1990, Icarus, 84, 226 [NASA ADS] [CrossRef] (In the text)
• Housen, K. R., & Holsapple, K. A. 1999, Icarus, 142, 21 [NASA ADS] [CrossRef]
• Kalas, P. 2005, ApJ, 635, L169 [NASA ADS] [CrossRef] (In the text)
• Kalas, P., & Jewitt, D. 1995, AJ, 110, 794 [NASA ADS] [CrossRef] (In the text)
• Kenyon, S. J., & Bromley, B. C. 2004, AJ, 127, 513 [NASA ADS] [CrossRef] (In the text)
• Kenyon, S. J., & Bromley, B. C. 2005, AJ, 130, 269 [NASA ADS] [CrossRef] (In the text)
• Kolokolova, L., Lara, L. M., Schulz, R., Stüwe, J. A., & Tozzi, G. P. 2001, Icarus, 153, 197 [NASA ADS] [CrossRef] (In the text)
• Lagrange, A.-M., Backman, D. E., & Artymowicz, P. 2000, Protostars and Planets IV, 639 (In the text)
• Li, A., & Greenberg, J. M. 1998, A&A, 331, 291 [NASA ADS]
• Liu, M. C. 2004, Science, 305, 1442 [NASA ADS] [CrossRef]
• Love, S. G., Hörz, F., & Brownlee, D. E. 1993, Icarus, 105, 216 [NASA ADS] [CrossRef]
• McDonnell, J. A. M., Lamy, P. L., & Pankiewicz, G. S. 1991, in Comets in the post-Halley era, ed. R. L. Newburn, Jr., M. Neugebauer, & J. Rahe, ASSL Vol. 167, IAU Colloq., 116, 1043 (In the text)
• Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896 [NASA ADS]
• Pantin, E., Lagage, P. O., & Artymowicz, P. 1997, A&A, 327, 1123 [NASA ADS]
• Petit, J.-M., & Farinella, P. 1993, Celestial Mechanics and Dynamical Astronomy, 57, 1 [NASA ADS] [CrossRef]
• Ryan, E. V., Davis, D. R., & Giblin, I. 1999, Icarus, 142, 56 [NASA ADS] [CrossRef] (In the text)
• Ryan, E. V., Hartmann, W. K., & Davis, D. R. 1991, Icarus, 94, 283 [NASA ADS] [CrossRef] (In the text)
• Ryan, E. V., & Melosh, H. J. 1998, Icarus, 133, 1 [NASA ADS] [CrossRef]
• Schneider, G., Silverstone, M. D., & Hines, D. C. 2005, ApJ, 629, L117 [NASA ADS] [CrossRef]
• Spangler, C., Sargent, A. I., Silverstone, M. D., Becklin, E. E., & Zuckerman, B. 2001, ApJ, 555, 932 [NASA ADS] [CrossRef] (In the text)
• Takeuchi, T., & Artymowicz, P. 2001, ApJ, 557, 990 [NASA ADS] [CrossRef] (In the text)
• Tanga, P., Cellino, A., Michel, P., et al. 1999, Icarus, 141, 65 [NASA ADS] [CrossRef]
• Telesco, C. M., Fisher, R. S., Wyatt, M. C., et al. 2005, Nature, 433, 133 [NASA ADS] [CrossRef] (In the text)
• Thébault, P., & Augereau, J.-C. 2005, A&A, 437, 141 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Thébault, P., Augereau, J. C., & Beust, H. 2003, A&A, 408, 775 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Tielens, A. G. G. M., McKee, C. F., Seab, C. G., & Hollenbach, D. J. 1994, ApJ, 431, 321 [NASA ADS] [CrossRef] (In the text)
• Wyatt, M. C. 2003, ApJ, 598, 1321 [NASA ADS] [CrossRef]
• Wyatt, M. C. 2005, A&A, 440, 937 [NASA ADS] [CrossRef] [EDP Sciences]
• Wyatt, M. C., & Dent, W. R. F. 2002, MNRAS, 334, 589 [NASA ADS] [CrossRef] (In the text)
• Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Nature, 373, 494 [NASA ADS] [CrossRef]

 

  
Online Material

   
Appendix A: Superparticles' structure

All SPs are cylinders with variable height $h_{\rm sp}$ and constant radius $r_{\rm sp}$. Their geometrical centers always stay in the disc midplane. For most of the runs we take $r_{\rm sp}=5$ AU. Test runs showed that an avalanche development weakly depends on the value for  $r_{\rm sp}$. The maximum amplification factor for runs with $r_{\rm sp}$ in the 4-8 AU range differs by less then $\sim$$10\%$. SPs are modeled in different ways depending on the physical origin of the grains they represent. We distinguish between 3 types of SP.

A.1 Field SP

All SPs that represent the initial structure of the dusty disc (i.e., non-avalanche SPs) have a superexponential vertical density profile

$$% n_{{{\rm gr}},i}= \frac{N_{{{\rm gr}},i}} {\pi r_{{{\rm sp... ...-0.5h_{{\rm sp},i}}^{0.5h_{{\rm sp},i}} f(z) {\rm d}z} f(z),$$ (A.1)

where $N_{{{\rm gr}},i}$ is the total number of grains in a given SP, $r_{\rm sp}~{\rm and}~h_{\rm sp}$ are the SP's radius and height, and f(z) is the SP profile function, which reads

$$% f(z)=\exp \left( - \left(\frac{\vert z\vert}{w_{{{\rm sp}},i}(R)}\right)^{p_{\rm z}}\right),$$ (A.2)

where $p_{\rm z}=0.7$ and $w_{{{\rm sp}},i}$ is the SP's width, which depends on the distance to the star, R, and is equal to the disc width (Eq. (10)) for field SPs. The height of a field SP is

$$% h_{{{\rm sp}},i}(R)=20w_{{{\rm sp}},i}(R).$$ (A.3)

All field SPs standing for the biggest grains (1 cm) have circular orbits. The SPs radial distribution and the number of grains in each SP is chosen in accordance with the best-fit parent body distribution from Augereau et al. (2001). All other grains are assumed to be produced from these biggest grains following the power law size-frequency distribution of Eq. (8). The number of SPs in each size bin is taken such that in the steady-state configuration there are 2-5 overlapping SPs of the same grain size at any given location in the system.

A.2 Initially released planetesimal debris

When a planetesimal is shattered the fragments are created with an initial spread in velocities. Ryan & Melosh (1998) show that the fragment velocities depend on grain sizes and that for small grains they are about a few percent of the impact velocity, with $v_{\rm fr} \approx 0.01{-}0.1 v_{\rm impact}$. Fragments produced from a body on a Keplerian orbit would spread in the vertical direction, and the height of the layer can be estimated as

$$% h_{\rm sp}\approx 2 \frac{v_{\rm fr}}{v_{\rm impact}} R.$$ (A.4)

For most of the runs we take $v_{\rm fr} /v_{\rm impact}=0.1$.

If the debris have velocities with isotropic distributions around the center of mass then the vertical distribution of grains of the same $\beta$ can be approximated as a constant. Thus these SPs are assumed to have no internal density structure and constant grain density.

$$% n_{{{\rm gr}},i}=\frac{N_{{{\rm gr}},i}}{\pi r_{{{\rm sp}},i}^2 h_{{{\rm sp}},i}},$$ (A.5)

which corresponds to f(z)= 1 (i.e., $w_{{{\rm sp}},i}= \infty$ in Eq. (A.1)).

   
A.3 SP created due to collisions

As mentioned in Sect. 3, when two SPs are passing through each other, a fraction of their grains can be destroyed producing new and smaller grains, which are then combined into new SPs. The number of grains that are destroyed in each SP, if the relative velocity is high enough to reach catastrophic fragmentation (see Sect. 3.4), is calculated as

$$% N_{\rm gr, -}= \int_0^{t_{\rm col}}\!\!\int_{\rm -z_{\rm0}}... ...i} n_{{{\rm gr}},j}v_{\rm rel} A_{\rm over}{\rm d}z {\rm d}t,$$ (A.6)

where $z_{\rm0}=0.5~{\rm min}(h_{{{\rm sp}},i},h_{{{\rm sp}},j})$, $\sigma=\pi (s_{{{\rm gr}},i} +s_{{{\rm gr}},j})^2/4$ is the collisional cross-section for grains with physical sizes $s_{{{\rm gr}}(i,j)}$, $v_{\rm rel}$ is the relative velocity of the grains, $A_{\rm over}$is the overlapping area of the two SPs viewed top-on, and  $t_{\rm col}$is the time while the SPs are passing through each other. The time dependences in Eq. (A.6) might be neglected if a reasonably small "collisional'' time step is taken for the calculations. We adopt $\Delta t_{\rm col}=0.02$ (in units of the orbital period at 20 AU) for most of the runs. Simulations with much smaller collisional time steps (e.g., $\Delta t_{\rm col}=0.005$) do not lead to a significant improvement of the results. The debris velocities and the size-frequency distribution for the newly created grains are identical to the values that one would get, considering a collision between $N_{\rm gr, -}$ pairs of grains with sizes $s_{\rm gr,(i,j)}$ and relative velocity  $v_{\rm rel}$(Sect. 3.4).

The structure of the new SPs (one SP for each size bin), which are produced after a collision, is obtained through the following equations. The initial height of the SP is equal to

$$% h_{{{\rm sp}},k}= min(h_{{{\rm sp}},i},h_{{{\rm sp}},j}).$$ (A.7)

The SP vertical profile f_k(z) is calculated as

f_k(z)=f_i(z) f_j(z), (A.8)

which corresponds to a SP of width (Eq. (A.2))

$$% w_{{{\rm sp}},k}=\left[ \left(\frac{1}{w_{{{\rm sp}},i}}\ri... ...{\rm sp}},j}}\right)^{-p_{\rm z}} \right]^{-1/p_{\rm z}}\cdot$$ (A.9)

In the case of dust production by SPs with superexponential profiles (Eq. (A.2)), the new SPs will have a smaller width than the two "colliding'' SPs. Keeping the same width would allow us to overlook the fact that "real'' grains in the SPs have vertical velocities and so will the produced debris. This component can be neglected if we only consider collisional outcomes, but it should be taken into account for calculations of the density structures. In our calculations this effect is taken into account by increasing the newly-created-SPs' width and height with time. The increase rate is taken to be

$$% \dot{w}_{\rm sp}=\frac{w_{\rm disc}-w_{\rm sp}}{\Delta t_w},$$ (A.10)

where $\Delta t_w$ is equal to  $\frac{1}{4}$ of the orbital period at the location where the SP is created. The same kind of time dependence is applied to the $h_{\rm sp}$ growth, namely

$$% \dot{h}_{\rm sp}=\left\{ \begin{array}{ll} \frac{h_{\rm max}-h_... ...m max}$ } \\ 0 &\mbox{if $ h_{\rm sp} \geq h_{\rm max},$ } \end{array} \right.$$ (A.11)

where $h_{\rm max}=h_0 w_{\rm disc}$, t_h=h₀ t_w, and h₀ is a constant parameter. The results for h₀=0.3 and h₀=3 differ by a factor of 2, which is acceptable for our order-of-magnitude calculations.

Appebdix B: Recombining SPs

To keep the total number of SPs manageable (we can trace the evolution of about one million of them) we have developed an algorithm that allows us to recombine SPs with similar parameters (grain size, velocities, positions in the disc). This allows us to avoid a too fast increase of the total number of SPs while keeping it to a value large enough from a statistical point of view.

Figure B.1: The total number of SPs as a function of time for a typical nominal case run.

The merging procedure is applied only between SPs standing for grains of the same size. The proximity condition is obtained by dividing the disc plane into a 2D space grid with the cell size equal to the SP's diameter. For each cell we list all SPs whose centers are located within the cell. Thus a given list contains SPs from the same spatial volume and with the same grain size. For each SP from a given list we calculate the "normalized'' velocity $v^*=v/v_{\rm kep}(R)$, where v is the SP velocity and R is the distance from the star to the SP center. If several SPs fall into one velocity bin (of width ${\rm d}v^*=10^{-2}$), they are combined into one SPs with the number of grains being equal to the sum of  $N_{\rm gr}$ of the combined SPs. The vertical structure of the new SP is obtained through the averaged values $w_{{\rm sp}}=\frac{\sum w_{{\rm sp},i} N_{{\rm gr},i}} {\sum N_{{\rm gr},i}}$ and $h_{{\rm sp}}=\frac{\sum h_{{\rm sp},i} N_{{\rm gr},i}} {\sum N_{{\rm gr},i}}$. The velocity and the position of the new SP are chosen to be equal to those of one SP, randomly chosenfrom the list of the recombined SPs. This SP is chosen randomly with the probability proportional to $\frac{ N_{{\rm gr},i}} {\sum N_{{\rm gr},i}}$. Choosing positions and velocities this way has the advantage of not introducing new trajectories for SPs, and this makes us more certain about the general shape of an avalanche. We have also tested another procedure for SPs recombination in which all SPs from the same velocity bin are recombined into one SP with position and velocity such that the total angular momentum and kinetic energy are preserved. In both cases we obtain similar results.

The recombination described above makes a significant reduction of the total number of avalanche SP spossible. However, an additional optimization procedure has been implemented to speed up the calculations. The idea is the following: the dust production rate is approximately proportional to the number density of avalanche grains (multiplied by the number density of the field population). As a consequence, regions with the lowest density of avalanche grains (of a given size) cannot make a significant contribution to avalanche dust production. Thus these regions are not very interesting for our simulations and there is no need to do very accurate representations of the avalanche grains' population there. A very rough representation is enough here. The question is to determine which number density values should be considered as "unproductive''. Since midplane densities for the field population vary by a factor $\sim$100 throughout the disc, regions where the density of avalanche grains (of a given size bin) are less than 1/1000that of the most dense regions will be considered as "not important''. In these regions all SPs from the same space cell (the proximity condition as defined above) representing grains of the same size bin are combined into one SP regardless of their velocities. Test runs have proved the efficiency of this optimization procedure: i) the general shape of an avalanche and the amplification factor evolution are not modified by it, since the "important'' regions are not affected; ii) the total number of SPs is significantly reduced, significantly speeding up the calculations. However, this procedure is limited to cases for which there is a significant density gradient. At later stages of avalanche propagation ($\ga$12 orbital periods) it is inefficient. Figure B.1 shows the evolution of the total number of SPs with time.