A&A 461, 537-549 (2007)
A. Grigorieva1 - P. Artymowicz1,2 - Ph. Thébault1,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
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 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 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: Pictoris - stars: individual: HD 32297
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., 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 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 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 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).
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
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
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
In a disc,
can be approximated by the optical thickness in the disc midplane,
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.
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 and variable heights , 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 ( 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 (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., ), 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 (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.
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:
The radiation pressure force is expressed as a function of the
gravitational force through the radiation pressure coefficient,
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 ( Pictoris-like) star. The thick lines represent silicate Mg Fe SiO. 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.|
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 5 104 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 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 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
have almost the same orbits as their parent bodies (the biggest grains), while smaller grains (i.e., with higher )
have more elliptic orbits
depending on their value. The initial number of grains as a function of their
size follows a classical single power law size-frequency distribution
The vertical structure of the disc is expressed in terms of the vertical
geometrical optical thickness,
length, z, as
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,
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).
Collision outcomes are traditionally divided into two classes:
(1) catastrophic fragmentation, when the largest remaining fragment,
is less than half of the parent body mass, M,
and (2) cratering, when
The energy per unit mass that is
needed to get
is called the threshold specific energy Q*.
If the specific energy
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
To account for the effect of different incidence angles, we
correct the value of Q* by a correction factor
corresponding to an average over all incidence angles
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 103 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 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.
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 M0=1020 g of dust is released in the m to 1 cm range at AU from the star, unless otherwise explicitly specified. It should be noted that the released dust mass M0 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 of the parent body whose shattering produces a mass M0 of dust. The ratio 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 and smaller fragments follow the Dohnanyi "equilibrium'' size distribution ( ), one gets 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 and steeper size distributions for highly disruptive impacts of large objects, with indexes typically in the -3.7 to -4 range for the largest 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, g = 10M0, which corresponds to an object of radius 40 km. We shall thus assume a nominal ratio of 0.1 for the discussion in Sect. 6.1. For the size spectrum of the dust particles released in the m to 1 cm range, we assume a single power law (Eq. (8) with ps0=3.5) for our nominal case. The dependence of the results on M0, ps0, and other parameters related to the planetesimal debris is explored in Sect. 4.2.
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, , 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 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 (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.
|material||Mg Fe SiO|
|porosity||compact grains (P=0)|
|grain density||3.5 g cm-3|
|Disc ("field'' population)|
|radial distribution||Augereau et al. (2001)|
|optical thickness along radius|
|in the midplane|
|disc extension||[20, 500] AU|
|Initial planetesimal debris|
|size distribution (see Eq. (8))|
|initial mass of dust released||M0=1020 g|
|distance from the star for|
|the planetesimal breakup||R0=20 AU|
|initial velocity of|
|the center of the mass|
|threshold energy, s0=1 cm||Q*0=107 erg g-1|
|power-law index (Eq. (11))||pQ=-0.2|
|size distribution of debris (see Eq. (12))|
|power-law indexes for m<ms||q1=1.5|
|power-law indexes for||q2=1.83|
|position of the slope change|
|Figure 2: Nominal case. Color-coded maps (log-scale) of the vertical optical thickness of avalanche grains, , 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.|
Figure 2 shows the temporal evolution of an avalanche, for the nominal case, in terms of the vertical optical thickness, , of the avalanche particles ( ). 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- (0.5) particles, which contribute to of . The maximum value of the amplification factor is and is reached at (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 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 are still very small compared to those of the field particles, with / 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
produced by the avalanche until time t to the initial number
of released planetesimal debris
(Fig. 4). As can be clearly seen,
a plateau, and we take
at t=15 as a reference value. Plugging values for the average number of particles produced by each grain-grain collision,
into Eq. (5), we get
i.e., a factor of 7 difference with the result of our simulation. This is mainly due to the fact that
underestimates the real value of ,
because the real path of a grain is curved rather
than parallel to a disc radius,
and secondly because in
the size of the avalanche grains
From our simulations, we were able to estimate the discrepancy
to be roughly of
a factor of 1.6. We thus get
|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).|
|Figure 4: Ratio of the total number of grains produced by the avalanche by the time t, , to the initial number of released planetesimal debris, N0.|
It is also interesting to link
to the amplification factor parameter
The initial mass M0 of dust released has been explored as a free parameter. The simulations show that the maximum amplification factor, , does not vary with M0, at least in the 1012 g-1021 g range, a result that is in good agreement with Eqs. (15) and (3). Likewise, does not change much when varying the initial speed of the center of mass of the planetesimal dust cloud. There is only a 20% increase of when is increased from to . 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.
For the nominal case we choose 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 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 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 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 m to be a reliable minimum value for and explored values in the m to m range, the latter value being the one considered by Kenyon & Bromley (2005). Although the impacting kinetic energy per grain is increasing in the m to m size range (leading to an increase of the ratio), decreases with increasing (Fig. 5) because of the decreasing value of the 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 . For smaller the variation with is less pronounced.|
|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)), m.|
Test runs have also been performed to check the 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 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 on the power-law index for the initial planetesimal debris size distribution, , is shown in Fig. 6. It can be noted that increasing above 3.5 (value for the nominal case) does not lead to a significant increase of , whereas less steep power laws lead to a significant decrease in .
We perform a set of runs in which the position of the planetesimal breakup, R0, 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 R0 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 ( ) 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.|
|Figure 8: Radial velocities vs. distance from the star for grains with different values released by parent bodies on circular orbits for a Pic-like star. The release distances are 20 AU and 70 AU.|
In the nominal case we assume that the minimum size, , for the debris produced by collisions and the minimum size of the initial planetesimal debris, , are both equal to m. We have seen in Sect. 4.2.2 that we do not expect significant changes when 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 (Fig. 9). As can easily be seen, does not strongly vary with for m. There are two reasons for that: i) decreases with decreasing sizes for grains smaller then 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 q1=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 ).|
Numerical exploration of the position of the slope change ms (with m) and of the m>ms power-law index q2 shows that the resulting amplification factor only weakly depends on these parameters. Changing from 1 to 10 leads to changes in by only a factor 2. On the contrary, variation in the q1 index can significantly affect , especially for q1>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).|
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 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 , 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*0 remaining constant. This constant Q*0 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 ( ) and decreases for porous grains ( ).
We numerically explore the importance of chemical composition by performing runs for the 2 extreme cases of pure (compact) silicates (Mg Fe SiO3) and pure (compact) water ices. Here again, we take the possibly counter-intuitive constant Q*0 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 values in the m range (see Fig. 1). It is thus not surprising that our results show no significant difference between the pure-ice ( ) and pure-silicate runs ( ).
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 106 and a few erg/g (e.g., Ryan et al. 1999; Holsapple et al. 2002; Benz & Asphaug 1999). We explore Q0* values between 106 and 108 erg/g and obtain strong variations in (Fig. 12). For the lowest explored Q0* value of 106 erg/g, we get , which is about 50 times higher than in the nominal case ( Q0* = 107 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.|
|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).|
|Figure 13: Different test radial distributions for ( being constant). The thick solid line is the distribution for the nominal case, taken from Augereau et al. (2001).|
As mentioned earlier, our reference field particle disc was assumed to be similar to the 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 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) remains the same. This result is in agreement with the simplified theory presented in Sect. 2.
On the other hand, we get drastic
variations when changing the value of
(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
which is a factor of 1000 higher. This strong increase in
is in agreement with
Eq. (15), which predicts a strictly exponential growth with
is constant. However, in the simulations we find that
weakly varies with
through the empirical relation
Plugging this expression for
into Eq. (15) we get:
|Figure 14: Maximum amplification factor as a function of . The open circles are the results of our simulations. The dashed line is the theoretical prediction (Eq. (16)).|
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 . However, the crucial issue is under which conditions such an avalanche might become observable. In this respect, looking at is not enough. What matters here is the ratio between the luminosity of the avalanche particles and that of the "field'' particles, .
The value of 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 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
The amount of light scattered towards an observer
coming from a given region of the disc is proportional to
|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.|
For avalanche detection, the visible domain (m) is probably more favorable than the NIR (m). This is because avalanches consist mostly of submicron grains, which scatter very inefficiently at 1-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 is nearly the same (within a factor of 2, depending on the exact chemical composition) for submicron and micron grains. Thus the ratio is expected to be higher in the visual than in the NIR. For simplicity we assume that is only a function of the scattering angle and that 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
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
for cometary dust, obtained from measurements
of solar system comets (Artymowicz 1997, and references therein):
As can be clearly seen in Fig. 15, in the nominal case 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.
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, remains constant with varying M0, so that the ratio increases linearly with M0. As a consequence, the release of 1022 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).
As has been seen in Sect. 4.4, increases strongly for grains with lower specific energy values Q*. The lowest Q* value explored in Fig. 11, erg/g, leads to / 0.4-0.5. Thus, observability might be marginally reached when assuming the minimum shattering resistance for dust grains.
|Figure 17: Same as Fig. 15, but for the case of a disc 5 times more massive than in the nominal case.|
The parameter exploration of Sect. 4 has clearly shown that the most efficient parameter for reaching high values is the field particles number density (Eq. (16)). An obvious way of increasing 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, , 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.|
Another way to reach higher values of
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
(Eq. (10)), a superexponential vertical profile
(Eq. (9)), and a corresponding
Assuming now a disc of thickness
with a constant vertical profile
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 Pic-like system and M0=1020 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 M0. 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 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 106 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
and the initial amount
of dust released M0, 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
criterion for observability might be written
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 -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 Pic analog with its SW side significantly brighter than the NE one within 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 -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 -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.
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:
All SPs are cylinders with variable height and constant radius . Their geometrical centers always stay in the disc midplane. For most of the runs we take AU. Test runs showed that an avalanche development weakly depends on the value for . The maximum amplification factor for runs with in the 4-8 AU range differs by less then . SPs are modeled in different ways depending on the physical origin of the grains they represent. We distinguish between 3 types of SP.
All SPs that represent the initial structure of the dusty disc
(i.e., non-avalanche SPs)
have a superexponential vertical density profile
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
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
If the debris have velocities with isotropic distributions around the center of mass then
the vertical distribution of grains of the same
can be approximated
as a constant. Thus these SPs are assumed to have no internal density structure and constant
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
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
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.
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 , 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 ), they are combined into one SPs with the number of grains being equal to the sum of of the combined SPs. The vertical structure of the new SP is obtained through the averaged values and . 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 . 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 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 (12 orbital periods) it is inefficient. Figure B.1 shows the evolution of the total number of SPs with time.