A&A 455, 509-519 (2006)
DOI: 10.1051/0004-6361:20064907

Dust distributions in debris disks: effects of gravity, radiation pressure and collisions[*]

A. V. Krivov1 - T. Löhne1 - M. Sremcevic2


1 - Astrophysikalisches Institute, Friedrich-Schiller-Universität Jena, Schillergäßchen  2-3, 07745 Jena, Germany
2 - LASP, University of Colorado, 1234 Innovation Drive, Boulder, CO 80303, USA

Received 24 January 2006 / Accepted 22 March 2006

Abstract
We model a typical debris disk, treated as an idealized ensemble of dust particles, exposed to stellar gravity and direct radiation pressure and experiencing fragmenting collisions. Applying the kinetic method of statistical physics, written in orbital elements, we calculate size and spatial distibutions expected in a steady-state disk, investigate timescales needed to reach the steady state, and calculate mass loss rates. Particular numerical examples are given for the debris disk around Vega. The disk should comprise a population of larger grains in bound orbits and a population of smaller particles in hyperbolic orbits. The cross section area is dominated by the smallest grains that still can stay in bound orbits, for Vega about 10  ${\rm\mu m}$ in radius. The size distribution is wavy, implying secondary peaks in the size distribution at larger sizes. The radial profile of the pole-on surface density or the optical depth in the steady-state disk has a power-law index between about -1 and -2. It cannot be much steeper even if dust production is confined to a narrow planetesimal belt, because collisional grinding produces smaller and smaller grains, and radiation pressure pumps up their orbital eccentricities and spreads them outward, which flattens the radial profile. The timescales to reach a steady state depend on grain sizes and distance from the star. For Vega, they are about 1 Myr for grains up to some hundred ${\rm\mu m}$ at 100 AU. The total mass of the Vega disk needed to produce the observed amount of micron and submillimeter-sized dust does not exceed several earth masses for an upper size limit of parent bodies of about 1 km. The collisional depletion of the disk occurs on Gyr timescales.

Key words: planetary systems: formation - circumstellar matter - meteors, meteoroids - celestial mechanics - stars: individual: Vega

1 Introduction

Debris disks are solar system-sized, gas-poor dust disks found around hundreds of main-sequence stars through infrared excesses and, in a dozen cases, directly resolved at different wavelengths, from the visual to submillimeter (see, e.g., Greaves 2005, for recent review). Simple estimates show that the lifetime of dust grains in these systems cannot exceed $\sim $1 Myr, implying that the dust material in debris disks is not primordial and must be steadily replenished by internal sources. It is widely believed therefore that debris disks are created and maintained by mutual collisions (Weissman 1984) and, possibly, comet-type activity of small bodies (Beust et al. 1989), similar to asteroids, comets, and Kuiper belt objects in the solar system. The dust material evolves then through a collisional cascade under the action of stellar gravity and radiation pressure forces.

Debris disks must be clearly distinguished from protoplanetary disks - denser disks with a high gas content around young T Tauri and Herbig Ae/Be stars. These disks are thought to be made of primordial material from which the central star has formed and planets should form, but the planet formation process is still ongoing. The physics of the dust component in protoplanetary disks are completely different (gas-driven dust dynamics, grain growth, etc.). These disks are not the subject of this paper.

Stimulated by a growing bulk of observational data on debris disks, substantial effort has been invested into their theoretical studies. Many authors focused on modeling of inner gaps and substructure found in most of the debris disks and attributed to the gravity of presumed planets (Moro-Martín & Malhotra 2002; Quillen & Thorndike 2002,Kuchner & Holman 2003; Wilner et al. 2002; Ozernoy et al. 2000; Wyatt et al. 1999; Wyatt & Dent 2002; Deller & Maddison 2005; Moran et al. 2004, among others). A number of studies dealt with dynamics of disks at transitional phases between protoplanetary and debris disks, exemplified by $\beta $ Pic, AU Mic, HR4796A and similar objects, which may have retained gas in moderate amounts, still sufficient to influence the dust dynamics (e.g., Takeuchi & Artymowicz 2001; Thébault & Augereau 2005). Kenyon & Bromley (2005) studied non-stationary processes in debris disks - consequences of major planetesimal collisions for the dust distributions.

Ironically, these - in fact more complicated - cases receive more attention than the underlying "regular'' debris disks. In this paper, we treat a debris disk as an ensemble of dust particles, exposed to stellar gravity and radiation and experiencing fragmenting collisions and try to identify essential, generic properties of such an idealized debris disk. The goal is to find out what kind of size/mass and spatial distributions of the disk material can be expected from theory, and how these distributions might depend on distributions of the parent bodies, parameters of the central star, and grain properties. We therefore do not try to construct realistic, but very sophisticated models of specific objects, and seek more general models. Studies of this kind were done, e.g., by Krivov et al. (2000), Dominik & Decin (2003), Thébault et al. (2003), Wyatt (2005).

Specific application is made to the disk around Vega. This disk is observed pole-on and does not show any clear substructure in the infrared images (Su et al. 2005), which makes it an ideal application of the model.

Section 2 describes the basic physical processes acting on the dust grains: stellar gravity, direct radiation pressure, and destructive intraparticle collisions. In Sect. 3, we present the kinetic model used to calculate the resulting dust distributions. Numerical results for the Vega disk are presented and discussed in Sect. 4. Section 5 lists our conclusions.

2 Physical processes and grain dynamics

2.1 Stellar gravity and radiation pressure

The main force acting upon macroscopic material and keeping it on closed orbits is the central star's gravity:

$\displaystyle \vec{F}\ensuremath{_{\rm grav}} = - \frac{GMm}{r^3} \vec{r},$     (1)

where G is the gravitational constant, M the stellar mass, and m the grain mass. Another force considered here is radiation pressure caused by the central star. At its simplest, it shows the same first-order r-2-behavior, but is pointed outward along $\vec{r}$. Because of this proportionality between the two forces the total photogravitational force can be written as
$\displaystyle \vec{F}\ensuremath{_{\rm pg}} = -\frac{G M (1-\beta) m }{r^3}\vec{r}.$     (2)

Here, $\beta $ is the ratio between radiation pressure and gravitational pull and is a constant depending on the grain size and optical properties.

In old systems with the optical depth $\tau$ less than roughly 10-5, exemplified by the presumed solar system's debris disk in the Kuiper belt region, the Poynting-Robertson (P-R) effect causes migration of smaller grains toward the primary star where they evaporate, while larger grains are typically lost to mutual collisions (Grün et al. 1985). If $\tau > 10^{-5}$ (all observable extrasolar disks), P-R drag is inefficient, as the collisional lifetimes are much shorter than the P-R times (Wyatt 2005; Krivov et al. 2000).

To evaluate the direct radiation pressure, one has to specify the properties of the grains to be treated. Here we assume compact grains of spherical shape that are characterized by their size or radius s, mass density $\rho$, and pressure efficiency $Q\ensuremath{_{\rm pr}} $. The last coefficient controls the fraction of the momentum transferred from the infalling radiation to the grain. The theoretical values range from 0 for perfect transmitters to 2 for perfect backscatterers. An ideal absorber's radiation pressure efficiency equals unity, the value adopted here. Then, the $\beta $-ratio is given by (Burns et al. 1979)

$\displaystyle \beta
=
\frac{3LQ\ensuremath{_{\rm pr}} }{16\pi GMcs\rho}
=
0.574...
...\cdot \frac{1 \: {\rm g} / {\rm cm}^3}{\rho}\cdot\frac{ 1\:\mu{\rm m}}{\rm s} ,$     (3)

where $L/L_\odot$ and $M/M_\odot$ are luminosity and mass of the star in solar units.

The smaller or the less dense the grains are, the more the radiation pressure they experience compensates the central star's gravity. Below the critical size where $\beta = 1$ the effective force is repelling, and the grains can no longer be held on bound orbits. Larger grains with $\beta
< 1$ orbit the star on Keplerian trajectories at velocities reduced by a factor of $\sqrt{1 - \beta}$ compared to macroscopic bodies ( $\beta \ll 1$), which are purely under the influence of gravitation. If, by fragmentation or any other erosive process, smaller particles are released from larger ones, in the first instance, they move along with their source due to a relative velocity close to zero. But the actual elements of their new orbits depend on their response to radiation pressure, which leads to separation of the fragments. A parent body on a circular orbit, for example, releases fragments into bound elliptic orbits with larger semimajor axes and eccentricities up to $\beta = 0.5$. Beyond this limit the grains are unbound and leave the system on normal ($\beta
< 1$) or "anomalous'' (i.e., open outward, $\beta > 1$) hyperbolic orbits (Fig. 1). For parent bodies in elliptic orbits, the boundaries between different types of orbits are somewhat smeared. Figure 2 depicts critical grain radii that separate particles in bound and hyperbolic orbits, as a function of the star's luminosity, for parent bodies with different orbital eccentricities. Borrowing the terminology from solar system studies (Zook & Berg 1975), we will call dust grains in bound and unbound orbits $\alpha$-meteoroids and $\beta $-meteoroids, respectively.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{4907f1.eps}
\end{figure} Figure 1: Three possible types of orbits of dust particles under the combined action of stellar gravity and direct radiation pressure. For illustrative purposes, grains are assumed to be released from a circular orbit.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8cm,clip]{4907f2.eps}
\end{figure} Figure 2: Grain radius that separates particles in bound and hyperbolic orbits, as a function of the star's luminosity (assuming dust bulk density of $2~\hbox{g}~\hbox{cm}^{-3}$, a unit radiation pressure efficiency, and a standard mass-luminosity $L\propto M^{3.8}$ relation for MS stars). Different linestyles are for different typical eccentricities of parent bodies. Grains between two lines of the same style may be in both types of orbits, depending on where between the periastron and apastron they are ejected. The upper lines separate bound and unbound orbits for ejection at the periastron and the lower ones correspond to ejection at the apastron.
Open with DEXTER

2.2 Binary collisions

In contrast to protoplanetary disks, collisions in optically thin, gas-poor debris disks occur at high relative velocities ($\ga$1 km s-1) and are, therefore, destructive and create smaller fragments. Removal of fine debris by stellar radiation pressure is the main loss "channel'' of material in such systems. So the questions to be addressed are: when or how often do grains collide? and: what happens if they collide?

The first question deals with the probability that two (spherical) objects on Keplerian orbits come close enough to touch. This probability can be expressed as a product of a purely geometric factor $\Delta $, describing the overlap of two orbits, the density of grains at the desired location in the space of orbital elements, the relative or impact velocity \ensuremath {\overline {v_{\rm imp}}}, and the cross section $\sigma$ for the collision, which is basically $\pi(s_{\rm p}+ s_{\rm t})^2$, where $s_{\rm p}$ and $s_{\rm t}$ are the radii of target and projectile, respectively (Krivov et al. 2005).

The answer to the second question depends on the physics of the impact and fragmentation process. As said above, the impact velocities easily exceed 1 km s-1, involving too much kinetic energy to allow the grains to stick together. Instead, the particles tend to destroy each other, producing smaller fragments, whose size distribution can be approximated by a classical power law $\ensuremath{{\rm d}} n/\ensuremath{{\rm d}} s \propto s^{-3.5}$ (Fujiwara 1986).

3 Disk model

3.1 Kinetic approach

Knowledge of single-particle dynamics and effects of individual binary collisions has now to be translated into the properties of the whole disk, considered as an ensemble of grains.

A straightforward, N-body approach - to follow dynamics of many individual objects and to perform true collision simulations - remains important for studying "difficult'' cases where many other methods fail, such as the final stages of planet formation (e.g., Charnoz et al. 2001; Ida & Makino 1993; Charnoz & Brahic 2001; Kokubo & Ida 1998). It can also be useful when the dynamics are complex, whereas any collisional event is assumed to simply eliminate both colliders (Lecavelier des Etangs et al. 1996). However, this method cannot treat a large number of objects sufficient to cover a broad range of particle masses. Further, it has an intrinsic problem in detecting collisions during the integration, which restricts its applicability. Hydrodynamics or SPH are not suitable for collision-dominated systems either.

The kinetic approach of statistical physics is more suitable. It is based on the continuity-like equation for the distribution of dust in an appropriate phase space. The idea is to introduce a multi-dimensional "phase-space distribution'' of dust, for instance a distribution of grain sizes, coordinates, and velocites, and to write down terms that describe the supply, loss and transport of dust grains due to different mechanisms. The resulting equation can be solved for the phase space distribution as a function of time, from which usual size and spatial distributions can easily be calculated. The Kinetic approach was used to model debris disks, including the interplanetary dust cloud and presumed Kuiper belt dust disk, by many authors (e.g., Gor'kavyi et al. 1997a; Rhee 1976; Ishimoto 1999; Dohnanyi 1978; Dohnanyi 1973; Ishimoto 2000; Gor'kavyi et al. 1997b; Southworth & Sekanina 1973; Leinert et al. 1983). However, in most of these papers the goal was to solve for either size or spatial distribution - not both. Also no realistic treatment of collisions was involved. An exception is our earlier paper (Krivov et al. 2000) in which, however, another major assumption has been made: small orbital eccentricities of dust grains. The velocity evolution of the dust material was not included either.

Our most recent model (Krivov et al. 2005) is free of the shortcomings listed above. However, it was developed for macroscopic objects rather than dust. As a specific application we addressed the collisional evolution of the classical population in the Edgeworth-Kuiper belt. In this paper, the model is generalized to include radiation pressure, which makes it applicable to dust disks.

3.2 Kinetic theory in orbital elements

Following Krivov et al. (2005), we choose the grain mass m and a vector of its orbital elements ${\vec k}$ as phase space variables. This is a distinct feature of our approach: usually, coordinates and velocities in a corotating "box'' are used instead (e.g. Greenberg et al. 1978; Weidenschilling et al. 1997; Spaute et al. 1991). To accommodate systems where large eccentricites and/or inclination cause large radial and/or vertical excursions of particles, standard codes employ multiannulus schemes, whereas our approach treats such cases in a straightforward way. Roughly speaking, we translate from the single-particle dynamics to an ensemble of particles using the idea known for centuries in classical celestial mechanics: to use slowly changing osculating elements instead of rapidly changing coordinates and velocities.

The disk is fully characterized by the distribution function $n(m,{\vec k})$, defined in such a way that $n(m,{\vec k}) \:\:\: \ensuremath{{\rm d}} m \; \ensuremath{{\rm d}} {\vec k}$is the number of particles with $[m, m+\ensuremath{{\rm d}} m]$, $[{\vec k}, {\vec k}+\ensuremath{{\rm d}} {\vec k}]$present in the disk at a certain instant of time t. The equation for the distribution function n has the form (Krivov et al. 2005)

 
$\displaystyle {\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t}
(m, {\vec k...
...}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm loss}} (m, {\vec k})$     (4)

with
  
                     $\displaystyle \left(
\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm gain}} (m, {\vec k})$ = $\displaystyle \int\limits_{m_{\rm p}}
\int\limits_{m_{\rm t}}
\int\limits_{{\ve...
...m_{\rm p},{\vec k}_{\rm p}; \;
m_{\rm t},{\vec k}_{\rm t}; \;
m, {\vec k}) \;\;$  
    $\displaystyle \times~
n(m_{\rm p},{\vec k}_{\rm p}) \;
n(m_{\rm t},{\vec k}_{\r...
...th{_{\rm imp}} (m_{\rm p},{\vec k}_{\rm p}; \; m_{\rm t},{\vec k}_{\rm t}) \;\;$  
    $\displaystyle \times~
\delta [{\vec r}({\vec k}_{\rm p}) - {\vec r}({\vec k}_{\rm t})] \;
\sigma (m_{\rm p},m_{\rm t}) \;$  
    $\displaystyle \times~
\ensuremath{{\rm d}} m_{\rm p}\;
\ensuremath{{\rm d}} m_{...
...\ensuremath{{\rm d}} {\vec k}_{\rm p}\;
\ensuremath{{\rm d}} {\vec k}_{\rm t}\;$ (5)
$\displaystyle \left(
\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm loss}} (m, {\vec k})$ = $\displaystyle n(m, {\vec k}) \;
\int\limits_{m_{\rm p}}
\int\limits_{{\vec k}_{...
... p}) \;
v\ensuremath{_{\rm imp}} (m_{\rm p},{\vec k}_{\rm p}, \; m,{\vec k}) \;$  
    $\displaystyle \times~
\delta [{\vec r}({\vec k}_{\rm p}) - {\vec r}({\vec k})] ...
...},m) \;
\ensuremath{{\rm d}} m_{\rm p}\;
\ensuremath{{\rm d}} {\vec k}_{\rm p}.$ (6)

Subscripts p and t mark two colliding particles, a projectile (the smaller of the two masses) and a target. Other quantities in Eqs. (4)-(6) are: $\sigma$ the collisional cross section, $\sigma(m_{\rm p}, m_{\rm t}) = \pi (s_{\rm p}+ s_{\rm t})^2$, where s is the radius of a particle with mass m; $v\ensuremath{_{\rm imp}} $ is the relative velocity of two particles $(m_{\rm p},{\vec k}_{\rm p})$ and $(m_{\rm t},{\vec k}_{\rm t})$ at collision; f is the fragment-generating function: ${f}(m_{\rm p},{\vec k}_{\rm p};\;
m_{\rm t},{\vec k}_{\rm t};\;
m, {\vec k})
\ensuremath{{\rm d}} m \ensuremath{{\rm d}} {\vec k}
$is the number of fragments with $[m, m+\ensuremath{{\rm d}} m]$, $[{\vec k}, {\vec k}+\ensuremath{{\rm d}} {\vec k}]$, produced by a collision of particles with $(m_{\rm p},{\vec k}_{\rm p})$ and $(m_{\rm t},{\vec k}_{\rm t})$. Finally, Dirac's $\delta [{\vec r}({\vec k}_{\rm p}) - {\vec r}({\vec k}_{\rm t})]$ensures that integrands are evaluated at collision.

3.3 Averaging

As discussed in Krivov et al. (2005), ${\vec k}$ must not necessarily consist of all six Keplerian elements. We can split all six into a vector of "important'' elements ${\vec p}$ and a vector of "unimportant'' elements ${\vec q}$. The letter ${\vec k}$ in Eq. (4)-(6) can simply be replaced with the letter ${\vec p}$. The equations will still be valid, only the interpretation of quantities slightly changes. The distribution n should now be understood as a distribution averaged over the remaining elements ${\vec q}$. Other quantities in integrands of (5) and (6), too, are averages over qs. In this case, we use $\overline{f}$ instead of f, \ensuremath {\overline {v_{\rm imp}}} instead of $v_{\rm imp}$, and $\Delta $ instead of $\delta$. They are formally defined as integrals over qs. For instance, Dirac's $\delta$ is replaced by

 
$\displaystyle \Delta ({\vec p}_{\rm p}, {\vec p}_{\rm t}) \equiv
\int\int_{{\ve...
...\; \varphi({\vec q}_{\rm t}) \;
{\rm d}{\vec q}_{\rm p}{\rm d}{\vec q}_{\rm t},$     (7)

where ${\vec q}_{\rm p}$ and ${\vec q}_{\rm t}$ are the variables to average over and $\varphi({\vec q}_{\rm p})$ and $\varphi({\vec q}_{\rm t})$ are their distributions normalized to unity.

An important question is which Keplerian elements of the six in total should be included in ${\vec p}$. The more elements that are included in ${\vec p}$, the more accurate the treatment, but the required computational resources grow rapidly with the increasing dimension of the phase space. In this paper, we choose an idealized debris disk that is rotationally-symmetric (uniformly distributed angular elements) and geometrically thin (small orbital inclinations). As a reasonable compromise, we therefore use ${\vec p}$ with two elements: semimajor axis a and eccentricity e. Altogether, we have a three-dimensional (m, a, e) phase space.

In the original model (Krivov et al. 2005), averaging over ${\vec q}$s (inclination i, longitude of ascending node $\Omega$, longitude of pericenter $\varpi$ and an anomaly) was "pre-done'' analytically before we started actual calculations with the model. For instance, we found an approximate expression for (7).

Here, we use a slightly different, and more accurate, method of averaging over $\varpi$s and anomalies. Consider two intersecting elliptic orbits (Fig. A.2). In the 2D approximation such two ellipses always cross twice or graze once. The relative orientation of the two orbits is given by the difference of their longitudes of pericenters: $\Lambda = \varpi _{\rm p}- \varpi _{\rm t}$. For this reason, we simply add one more argument, $\Lambda$, to the functions $\overline{f}$, $\Delta $, and $\ensuremath{\overline{v_{\rm imp}}} $, and perform one additional (innermost) integration over $\Lambda$, in Eqs. (5) and (6). Having $\Lambda$ as an argument of functions is convenient because $\Lambda$ immediately determines the true anomaly

 
$\displaystyle \vartheta_{\rm t}
=
\arcsin\frac{p_{\rm t}e_{\rm p}\sin\Lambda}{\...
... p_{\rm t}^2 e_{\rm p}^2 - 2p_{\rm p}p_{\rm t}e_{\rm p}e_{\rm t}
\cos\Lambda}},$     (8)

which, in turn, determines the distance r, through equation of conic section (14). Then, \ensuremath {\overline {v_{\rm imp}}} has a simple interpretation: it is now the average relative velocity of two particles colliding at distance r to the star. The Interpretation of functions $\overline{f}$ and $\Delta $ changes in a similar way.

The final equations have the form

 
$\displaystyle {\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t}
(m, {\vec p...
...}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm loss}} (m, {\vec p})$     (9)

with
  
                   $\displaystyle \left(
\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm gain}} (m, {\vec p})$ = $\displaystyle \int\limits_{m_{\rm p}}
\int\limits_{m_{\rm t}}
\int\limits_{{\ve...
...{\vec p}_{\rm p}; \;
m_{\rm t},{\vec p}_{\rm t}; \;
m, {\vec p};\;\Lambda) \;\;$  
    $\displaystyle \times~
n(m_{\rm p},{\vec p}_{\rm p}) \;
n(m_{\rm t},{\vec p}_{\rm t}) \;\;$  
    $\displaystyle \times~
\ensuremath{\overline{v_{\rm imp}}} (m_{\rm p}, {\vec p}_{\rm p}; \; m_{\rm t}, {\vec p}_{\rm t};\;\Lambda) \;\;$  
    $\displaystyle \times~
\Delta ({\vec p}_{\rm p}, {\vec p}_{\rm t};\;\Lambda) \;
\sigma (m_{\rm p},m_{\rm t})$  
    $\displaystyle \times~
\varphi (\Lambda) \;
\ensuremath{{\rm d}} m_{\rm p}\;
\en...
...{\rm p}\;
\ensuremath{{\rm d}} {\vec p}_{\rm t}\;
\ensuremath{{\rm d}}\Lambda\;$ (10)
$\displaystyle \left(
\ensuremath{{\rm d}} n \over \ensuremath{{\rm d}} t
\right)\ensuremath{_{\rm loss}} (m, {\vec p})$ = $\displaystyle n(m, {\vec p}) \;
\int\limits_{m_{\rm p}}
\int\limits_{{\vec p}_{\rm p}}
\int\limits_{\Lambda}
n(m_{\rm p},{\vec p}_{\rm p}) \;$  
    $\displaystyle \times~
\ensuremath{\overline{v_{\rm imp}}} (m_{\rm p}, {\vec p}_{\rm p}; \; m, {\vec p} ;\;\Lambda) \;\;$  
    $\displaystyle \times~
\Delta ({\vec p}_{\rm p}, {\vec p}; \;\Lambda) \;
\sigma (m_{\rm p},m) \;$  
    $\displaystyle \times~
\varphi (\Lambda) \;
\ensuremath{{\rm d}} m_{\rm p}\;
\ensuremath{{\rm d}} {\vec p}_{\rm p}\;
\ensuremath{{\rm d}}\Lambda .$ (11)

Due to the assumed axisymmetry, the distribution over $\Lambda$ is uniform, so that $\varphi (\Lambda) = 1/(2\pi)$. In the subsequent subsections, we discuss in turn the quantities that appear in the integrands in Eqs. (9)-(11): orbital elements and mass of fragments generated in a binary collision, needed to specify the function $\overline{f}$; collisional probabilities that determine $\Delta $; and impact velocities $\ensuremath{\overline{v_{\rm imp}}} $.


  \begin{figure}
\par\includegraphics[width=14cm]{fi3.eps} \end{figure} Figure 3: Fragments' orbits produced by collision of two particles ( $m_{\rm t}= 2m_{\rm p}= 2$, $a_{\rm t}=1$, $e_{\rm t}=0.9$, $\beta _{\rm t}=0$, $a_{\rm p}=1.5$, $e_{\rm p}=0.26$, $\beta _{\rm p}=0$) whose locations are shown by the two boxes. Lines of equal $\beta $-ratio (at steps of 0.1 from 0.0 to 2.0, $\beta = 1$ is excluded) and a running ejection point are shown, as are lines of equal ejection point (at true anomaly steps of 45$^\circ $ starting at the pericenter) and running $\beta $. The upper-left part contains anomalous hyperbolas, the ellipses are located in the middle, and the normal hyperbolas are found at the lower-right. The constant GM is set to unity.
Open with DEXTER

3.4 Orbital elements of collisional fragments

Initially, a cloud of fragments, produced by a disruptive collision, carries the sum of the momenta of both colliders in one direction. The fragments have zero relative velocities, as we assume maximal collisional damping. The resulting momentum is that of the center of mass and, in the 2D case, is described by two conservation laws, one for the radial component and one for the angular:

$\displaystyle m\ensuremath{\ensuremath{_{\rm sum}}}\dot r = m_{\rm p}\dot r_{\rm p}+ m_{\rm t}\dot r_{\rm t},$     (12)
$\displaystyle m\ensuremath{\ensuremath{_{\rm sum}}} r\dot\vartheta= m_{\rm p}r \dot\vartheta_{\rm p}+ m_{\rm t}r\dot\vartheta_{\rm t},$     (13)

where $m\ensuremath{\ensuremath{_{\rm sum}}} = m_{\rm p}+ m_{\rm t}$. Krivov et al. (2005) considered the angular momentum only. The fragments were all placed on the orbit of the center of mass and remained there due to their being too large to show significant response to radiation pressure. Now that we concentrate on micron-sized dust, we extend the fragment-generating function to consider the distance $r=r_{\rm p}=r_{\rm t}$, at which the parent particles p and t collide, and to include radiation pressure described by the $\beta $-ratios. Using the equation of conic section
 
$\displaystyle r = {p \over {1 + e\cos\vartheta}}$     (14)

with
p = a(1-e2),     (15)

and the derivatives
$\displaystyle r\dot\vartheta= {L \over mr},\qquad
\dot r = {L \over mr}\left({1 \over r}{\partial
r \over \partial\vartheta}\right)$     (16)

with
 
$\displaystyle L = m\sqrt{GM(1-\beta)p},$     (17)
$\displaystyle {1 \over r}{\partial r \over \partial\vartheta} = \pm \sqrt{{r \over
p}\left(2 - {r\over a} - {p\over r}\right)}$     (18)

for the colliders and the resulting fragments, a lengthy, but straightforward algebra results in
  
                          $\displaystyle \frac{r}{a}$ = $\displaystyle 2 - \frac{m_{\rm p}^2}{m\ensuremath{\ensuremath{_{\rm sum}}} ^2}\...
...} ^2}\cdot \frac{1-\beta_{\rm t}}{1-\beta}
\left[2 - \frac{r}{a_{\rm t}}\right]$  
    $\displaystyle - 2\frac{m_{\rm p}m_{\rm t}}{m\ensuremath{\ensuremath{_{\rm sum}}...
...{\rm p})(1-\beta_{\rm t})}}{1-\beta}
\Bigg[\frac{1}{r}\sqrt{p_{\rm p}p_{\rm t}}$  
    $\displaystyle \pm \sqrt{\left(2 - \frac{r}{a_{\rm p}} - \frac{p_{\rm p}}{r}\right)
\left(2 - \frac{r}{a_{\rm t}} - \frac{p_{\rm t}}{r}\right)}\Bigg],$ (19)
e = $\displaystyle \pm \left[1 - \frac{1}{a}\left(\frac{m_{\rm p}^2}{m\ensuremath{\e...
...suremath{_{\rm sum}}} ^2}\frac{1-\beta_{\rm t}}{1-\beta}p_{\rm t}\right.\right.$  
    $\displaystyle \left.\left. + 2\frac{m_{\rm p}m_{\rm t}}{m\ensuremath{\ensuremat...
...\beta_{\rm t})}}
{1-\beta} \sqrt{p_{\rm p}p_{\rm t}}\right)\right]^\frac{1}{2}.$ (20)

Here, the sign of e is equal to that of $(1-\beta)$, yielding anomalous hyperbolas with e < 0 for $\beta > 1$. Given $a\ensuremath{_{\rm p,t}} $, $e\ensuremath{_{\rm p,t}} $, $\beta\ensuremath{_{\rm p,t}} $, and $m\ensuremath{_{\rm p,t}} $ for the colliders, we are now able to calculate a and e for all fragments (all $\beta $) and positions r. The influence of $\vartheta$ (or r) and $\beta $ on the fragments' orbital elements is demonstrated in Fig. 3.

3.5 Masses of collisional fragments

The masses of fragments produced by the collision are supposed to range from arbitrarily small vapor particles to a certain upper limit $m\ensuremath{_{\rm x}} $, which is determined by the properties of the parent bodies and their impact velocity (e.g., Paolicchi et al. 1996):

 
$\displaystyle m\ensuremath{_{\rm x}} =
\frac{m_{\rm t}}{2}\left[2\frac{m_{\rm t...
...math{_{\rm D}} ^*(m_{\rm t})}{\ensuremath{\overline{v_{\rm imp}}} ^2}\right]^c,$     (21)

with c close to unity. The higher the velocity, the smaller the fragments. Following Krivov et al. (2005) the critical specific energy $Q\ensuremath{_{\rm D}} ^*$ is composed of two power laws, one for the strength regime and one for the gravitational regime:
 
$\displaystyle Q\ensuremath{_{\rm D}} ^* = A\ensuremath{_{\rm s}} s_{\rm t}^{b\e...
...emath{_{\rm s}} } + A\ensuremath{_{\rm g}} s_{\rm t}^{b\ensuremath{_{\rm g}} }.$     (22)

Its minimum and its transition zone are reported to lie at around 100 m, and thus, the typical dust grains and parent bodies would be sufficiently described by the strength regime alone. The energy $m_{\rm t}Q\ensuremath{_{\rm D}} ^*(m_{\rm t}) + m_{\rm p}Q\ensuremath{_{\rm D}} ^*(m_{\rm p})$ is needed to break up both target and projectile. If the kinetic energy of the inelastic impact is higher than this critical energy, the collision is treated as fragmenting. Thus, at a given impact velocity, only projectiles of a mass $m_{\rm p}\ga m\ensuremath{_{\rm cr}} $ are considered, and as usually $m\ensuremath{_{\rm cr}}\ll m_{\rm t}$, this reduces to
$\displaystyle m\ensuremath{_{\rm cr}}\approx {2 m_{\rm t}\over \ensuremath{\overline{v_{\rm imp}}} ^2}Q\ensuremath{_{\rm D}} ^*(m_{\rm t}).$     (23)

In agreement with previous works, the mass distribution function in this range is a power law $m^{-\eta}$ or $s^{-\nu}$. While experimental results give a range $1.5 < \eta < 2.0$ (Fujiwara 1986, and references therein), the "classical'' value is 11/6, corresponding to $\nu = 3.5$.

3.6 The fragment-generating function

The fragment-generating function $\overline{f}$ that enters Eq. (10) describes the number of fragments with $(m, {\vec p})$ that are produced by a destructive collision of parent bodies with $(m_{\rm p}, {\vec p}_{\rm p})$and $(m_{\rm t}, {\vec p}_{\rm t})$. It is a product of two distributions:

$\displaystyle \overline{f}(m_{\rm p},{\vec p}_{\rm p};\;
m_{\rm t},{\vec p}_{\r...
...,{\vec p};\;m_{\rm p},{\vec p}_{\rm p};\;m_{\rm t},{\vec p}_{\rm t};\;\Lambda).$     (24)

Here, $\overline g$ represents the mass distribution mentioned above,
$\displaystyle \overline g = (2 - \eta)\cdot(m_{\rm p}+ m_{\rm t})\cdot\left({m\...
...math{_{\rm x}}\over m}\right)^\eta\cdot{1 \over m\ensuremath{_{\rm x}} ^2}\cdot$     (25)

The second function, $\overline h$, is the orbital element distribution. This function evaluates whether or not a fragment with the mass m, produced in a collision between $(m_{\rm p}, {\vec p}_{\rm p})$ and $(m_{\rm t}, {\vec p}_{\rm t})$at the apsidal angle $\Lambda$, can have orbital elements ${\vec p}=(a,e)$. Therefore, $\overline{h} (m,a,e;\;m_{\rm p},{\vec p}_{\rm p};\;m_{\rm t},{\vec p}_{\rm t};\;\Lambda)$is proportional to

\begin{displaymath}\delta[a - a(m_{\rm p},{\vec p}_{\rm p};\:m_{\rm t},{\vec p}_{_{\rm t}};\:m;\:\Lambda)]
\end{displaymath} (26)

and

\begin{displaymath}\delta[e - e(m_{\rm p},{\vec p}_{\rm p};\:m_{\rm t},{\vec p}_{_{\rm t}};\:m;\:\Lambda)],
\end{displaymath} (27)

where $a(m_{\rm p},{\vec p}_{\rm p};\;
m_{\rm t},{\vec p}_{\rm t};\;
m;\;\Lambda)$and $e(m_{\rm p},{\vec p}_{\rm p};\;
m_{\rm t},{\vec p}_{\rm t};\;
m;\;\Lambda)$are given by Eqs. (19) and (20). The function $\overline h$ is normalized to unity:

\begin{displaymath}\int \overline{h} (m,a,e;\;m_{\rm p},{\vec p}_{\rm p};\;m_{\r...
...ambda)\: \ensuremath{{\rm d}} a \: \ensuremath{{\rm d}} e = 1.
\end{displaymath} (28)

3.7 Collisional probabilities and impact velocities

The "geometric probability'' of collision $\Delta ({\vec p}_{\rm p}, {\vec p}_{\rm t};~\Lambda)$(dimension: reciprocal of volume) and the impact velocity $\ensuremath{\overline{v_{\rm imp}}} (m_{\rm p}, {\vec p}_{\rm p}; \; m_{\rm t}, {\vec p}_{\rm t};\;\Lambda)$for two given elliptic orbits with a given angle $\Lambda$ between their apsidal lines are

  
$\displaystyle \Delta = \frac{1}{v_{\rm p}v_{\rm t}P_{\rm p}P_{\rm t}2r \sin\varepsilon
\vert\sin\gamma\vert},$     (29)
$\displaystyle \ensuremath{\overline{v_{\rm imp}}} = \sqrt{v_{\rm p}^2 + v_{\rm t}^2 - 2v_{\rm p}v_{\rm t}
\cos\overline\gamma}.$     (30)

The variables are: the orbital velocities $v_{\rm p}$ and $v_{\rm t}$, the orbital periods $P_{\rm p}$ and $P_{\rm t}$, the angle $\gamma$ at which the orbits cross, and the disk's semi-opening angle $\varepsilon$. The distance of collision r is calculated with Eqs. (8) and (14).

To account for the third dimension, the impact velocity is calculated using an average crossing angle:

$\displaystyle \cos\overline{\gamma} = \frac{\cos\gamma}{\sqrt{1+(B\sin\varepsilon)^2}},$     (31)

where B is an empirical weighting constant set to 2/3.

For a derivation of these equations see the Appendix.

3.8 Numerical implementation

To solve Eqs. (9)-(11) numerically, we discretize them by introducing a 3D-grid of m, a, and e and by replacing integrations with summations, and evolve the system in time with a first-order Euler routine with an adaptive stepsize control, based on analysis of increments of the function n.

The collisional probabilities $\Delta $ and velocities \ensuremath {\overline {v_{\rm imp}}} are, as far as possible, computed only once, and then used to calculate the gain and loss terms for each time step, describing the destruction of parent bodies and the creation of fragments. Explicitly, for every bound orbit ( $a_{\rm t},e_{\rm t}$) the distance r at which a possible collision with another orbit ( $a_{\rm p},e_{\rm p}$) would take place is calculated at equidistant steps of the true anomaly $\vartheta_{\rm t}$. If the projectile p reaches this distance during its revolution around the star, the following quantities are calculated: the true anomalies $\vartheta_{\rm p}{\ensuremath{^+} }$ and $\vartheta_{\rm p}{\ensuremath{^-} }$ (Fig. A.2), together with $\Delta\ensuremath{^{\rm\pm}} $, the instantaneous orbital velocities $v_{\rm t}$ and $v_{\rm p}$ for the case where $\beta_{\rm t}= \beta_{\rm p}= 0$, and the effective crossing angles $\overline\gamma\ensuremath{^{\rm\pm}} $. These values, which are described in the Appendix, are stored. Then at each time step for each crossing combination of orbits and each combination of colliders' masses the discrete range of anomalies $\vartheta_{\rm t}$ with stored values is processed. The orbital velocities are properly scaled by $(1 - \beta_{\rm t})^{1/2}$ and $(1 - \beta_{\rm p})^{1/2}$ to determine the relative impact velocity, and momentum conservation is applied to obtain the orbits of all possible fragments up to the limiting mass $m\ensuremath{_{\rm x}} $ which is calculated using \ensuremath {\overline {v_{\rm imp}}}. A linear interpolation between the $\vartheta_{\rm t}$-steps is used to distribute the fragments over the (m,a,e)-grid. This selection of fragment orbits and the distribution of mass together give the implementation of the fragment-generating function $\overline{f}$.

The abundance of grains on unbound orbits, $\beta $-meteoroids, is calculated as the product of their production rate and the average time it takes for them to leave the disk, starting from the pericenter of their orbits.

4 Application to the Vega disk

4.1 Observations

The so-called Vega phenomenon for main-sequence stars refers to a mid-infrared excess in the observed spectrum over the purely stellar emission. It was discovered in 1983 by the IRAS mission, at first for Vega itself (Aumann et al. 1984), and was soon attributed to a population of cold dust surrounding the star. From the analysis of the spectral energy distribution in the spectral range from 12  ${\rm\mu m}$ (Aumann et al. 1984; Beichman et al. 1988) to 1.3 mm (Chini et al. 1990), conclusions about the properties of the assumed disk were made, for instance the existence of an inner gap was inferred from the absence of hot grains. The first spatially resolved submillimeter images of the Vega disk were obtained with the SCUBA camera on the JCMT on Hawaii (Holland et al. 1998), followed by other groups (Koerner et al. 2001; Wilner et al. 2002), who concentrated on finding asymmetries. Recently, Su et al. (2005) observed the Vega system at 24, 70, and 160  ${\rm\mu m}$ with the Spitzer space telescope and found radial profiles of the Vega disk to be nearly rotationally symmetric and featureless. The absence of asymmetries and substructure in the infrared makes the Vega disk an ideal application of our approach.

4.2 Input data

We set the initial size distribution to the well-known relation $n(s) \propto s^{-3.5}$, starting at a minimum radius of ${\approx} 0.1$  ${\rm\mu m}$. The total mass strongly depends on the upper mass cut-off (proportional to $\sqrt{m\ensuremath{_{\rm max}} }$). To conform with previous mass estimates (Holland et al. 1998; Su et al. 2005), the total observable mass was set to be made up of grains up to a limiting radius of 1.5 mm. Grains larger than this - the population of parent bodies - follow the same power law, but deliver additional mass. The mass range used was 10-14 g to $3.56\times 10^{14}$ g. Thus, setting an observable mass of $\approx $0.5  $M\ensuremath{_{\rm lunar}} $ (lunar masses) leads to a total initial mass of close to 300  $M\ensuremath{_{\rm lunar}} $ or 3.5 $M_\oplus$ (earth masses). This range was covered by a grid of 70 bins with logarithmic steps with a step factor of 2.6. The equivalent size grid has steps with factors of 1.37.

To model the distribution of orbital elements, from which then the spatial distribution of dust material was calculated, we used a mesh of the semimajor axes from -400 to 400 AU with steps of 20 AU, and eccentricity bins from -3 to 3 with width of 0.125, as depicted in Fig. 8. According to Su et al. (2005) and Dent et al. (2000), for the Vega disk we adopted an initial distribution of semimajor axes between 80 AU to 120 AU, chosen in such a way that for zero eccentricities the normal optical depth $\tau$ in this distance range would be constant. The distribution of semimajor axes outside a = 120 AU was taken to be a power law, corresponding to the optical depth $\tau(a) \propto a^{-\alpha}$. Different values of the index $\alpha$ were tested, see below. No material was placed initially inside a = 80 AU. The initial distribution of eccentricities both in the "ring'' between 80 AU and 120 AU and in the outer disk outside 120 AU was set uniform from zero to an upper limit of $e\ensuremath{_{\rm max}} = 1/8$ or 3/8, depending on the run. The vertical extension of the disk is defined by a full-opening angle, which we arbitrarily set to $2\varepsilon = 17^\circ $.

The optical and mechanical properties of dust in the disk are wildly unknown, and we use here two types of grains: "rocky'' and "icy'' ones. Assuming a mixture of 70% astronomical silicate of 3.5 g/cm3(Laor & Draine 1993), 30% amorphous carbon of 1.85 g/cm3(Zubko et al. 1996) and porosity close to zero, the effective grain bulk density used here for "rocky'' grains is 3.0 g/cm3. The mean radiation pressure efficiency $\langle Q\ensuremath{_{\rm pr}}\rangle$, describing the momentum transfer from radiation to grain, is set to unity. The dependence of the critical energy for disruption on the grain size is given by Eq. (22). The value of $Q\ensuremath{_{\rm D}} ^*$ at a size of 1 m was set to 106 erg/g, and the slope $b\ensuremath{_{\rm s}} $ for this strength regime was set to -0.24, although these values are not well defined by theoretical and laboratory work. To check how strongly the results depend on the adopted material, one run was done for "icy'' grains of low bulk density (1 g/cm3) and mechanical strength ( $Q\ensuremath{_{\rm D}} ^* = 2\times 10^5$ erg/g at a size of 1 m). In Eq. (21), we used c=1.24 for rock (Paolicchi et al. 1996) and c=0.91 for ice (Arakawa 1999).

4.3 Results

The resulting dust distributions show strong fingerprints of the radiation pressure. The most noticeable is the peak at the particle size where $\beta = 1$, so that the radiation pressure equals gravity (cf. Fig. 2). Below this size bound orbits are impossible (Fig. 4), and the grains are blown away on a timescale of the order of 102 to 103 years. Due to this lack of possible impactors, grains slightly above this size limit are overabundant, thereby reducing the number of grains of the next larger population, and so on. This dependence produces a well-known wavy pattern in the size or mass distribution (e.g., Campo Bagatin et al. 1994; Durda & Dermott 1997; Thébault et al. 2003), whose "wavelength'' depends on the ratio of the average impact energy available and the impact energy $Q\ensuremath{_{\rm D}} ^*$ needed to disrupt a given target. The first is controled by impact velocities, which depend on the disk's layout. The range of inclinations and especially the range of eccentricities are important, because they determine the orbits' crossing angles. The critical energy $Q\ensuremath{_{\rm D}} ^*$, on the other hand, is determined by material properties.


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{$\frac{\mbox{cross section density}...
...ain radius [$\rm\mu m$ ]}
\includegraphics[width=8.8cm]{4907f4.eps} \end{figure} Figure 4: Grain size distribution at 100 AU after 10 Myr for a rocky disk with an initial outer profile ${\propto } r^{-4}$ and eccentricities from 0 to 0.375 (solid line). The contributions from the three different types of orbits, shifted down by one order of magnitude for better visibility, are shown by long-dashed (ellipses), short-dashed (normal hyperbolas), and dotted (anomalous hyperbolas) lines.
Open with DEXTER


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{$\frac{\mbox{cross section density}...
...ain radius [$\rm\mu m$ ]}
\includegraphics[width=8.8cm]{4907f5.eps} \end{figure} Figure 5: Grain size distribution at 100 AU and 10 Myr (4.5 Myr for "ice'') for different eccentricity ranges: 0.0 to 0.375 for "rock'' and "ice'', and only 0.0 to 0.125 for "circular'' orbits. The dashed lines are rescaled to coincide with the solid at large radii. The horizontal shift of the maximum for icy grains is largely due to a different bulk density.
Open with DEXTER


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{$\frac{\mbox{cross section
density}...
...ain radius [$\rm\mu m$ ]}
\includegraphics[width=8.8cm]{4907f6.eps} \end{figure} Figure 6: Evolution of the size distribution of the disk of Fig. 4 at a distance of 100 AU.
Open with DEXTER


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{$\frac{\mbox{cross section
density}...
...ain radius [$\rm\mu m$ ]}
\includegraphics[width=8.8cm]{4907f7.eps} \end{figure} Figure 7: Grain size distribution at distances from 50 AU to 350 AU in steps of 50 AU after 30 Myr. The initial disk is the same as in Fig. 4. The inner boundary of semimajor axes was placed at 80 AU. Therefore, only a small fraction of grains on sufficiently eccentric orbits contribute to the density at 50 AU.
Open with DEXTER

As is shown in Fig. 5 the disk's and the dust's set-up can have a noticeable influence on the size distribution. Low maximum eccentricities shorten the pattern's wavelength due to lower impact velocities leading to the critical impactor's size coming closer to the given target's. In contrast, the fluffy constitution of icy grains enlarges the gap between the target mass and the mass needed to disrupt it at a given impact velocity. Furthermore, it shifts the lower cut-off of the particle radius to higher values, according to the ratio of bulk densities, which is here 3:1. We note that a realistic disk cannot be built up of a perfectly homogeneous material, which implies a dispersion of densities, fragmentation energies etc. This could weaken or smear the waviness of the size distribution.

The time evolution of a particle size distribution of a modeled Vega disk is shown in Fig. 6. The larger the particles, the longer the collisional timescales, and the more time it takes for the size distribution to deviate significantly from the initial power law and to converge towards a quasi-steady state distribution whose absolute values are subject to change due to the amount of material delivered by collisions of larger bodies. Similarly, the comparison of different radial distances at the same time (Fig. 7) shows slower evolution at regions that are farther out with lower number densities and decreased impact velocities, where the latter also result in a reduced wavelength of the pattern.

Because we use the semimajor axis a and the eccentricity e as variables, all quantities involving the radial distance r are computed by integrating the distributions (cf. Krivov et al. 2005, Eq. (2.18)) over all the orbits crossing or grazing this distance. The actual distribution of the cross sectional area over the phase space of a and e is plotted in Fig. 8 for rocky material with an initial maximum eccentricity of 0.375 and a radial slope $\alpha = 2$. Near the inner edge of the disk this surface area is dominated by eccentricities in the range of the parent population's, whereas in the outer region the maximum contribution comes from small grains on more eccentric orbits, because a higher eccentricity at the same semimajor axis implies an origin closer to the star, where the density of parent bodies is higher. Thus, the fragments created near the disk's inner edge produce the main portion of optical depth out to a considerable distance, setting an upper limit to the radial slope $\alpha$. Even a ring-like or toroidal parent population with a sharp outer cut-off would lead to a slope not steeper than this.


  \begin{figure}
\par\psfrag{eccentricity}[bc][bc]{Eccentricity $e$ }
\psfrag{sem...
...ajor axis $a$\space [AU]}
\includegraphics[width=8.8cm]{4907f8.eps} \end{figure} Figure 8: The phase space distribution for the Vega disk of Fig. 4 at a quasi-steady state after 10 Myr. Each bin in the e-a-grid represents the surface area (in cm2 per phase space bin-volume) covered by the grains of all masses belonging to it.
Open with DEXTER


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{normal optical depth}
\psfrag{XAxi...
...c]{distance to Vega [AU]}
\includegraphics[width=8.8cm]{4907f9.eps} \end{figure} Figure 9: Radial profiles of the normal optical depth after 10 Myr: one for an initially ring-like disk only (solid), one for the ring followed by the semimajor axis distribution ${\propto } a^{-4}$ (dashed), and two for the ring followed by the ${\propto } a^{-2}$ profile (dotted for a maximum eccentricity $e\ensuremath {_{\rm max}} = 0.375$ and dash-dotted for $e\ensuremath {_{\rm max}} = 0.125$). The ring boundaries at 80-120 AU are shown by vertical lines. Note that these values are the minimum and maximum semimajor axes, so that the ring in space extends somewhat inside r=80 AU and outside r=120 AU.
Open with DEXTER

The radial distribution of the optical depth at the end of the integration is shown in Fig. 9 for four different runs, corresponding to the source ring of constant surface density with and without adjacent outer part with different slopes and assuming different eccentricity distributions. As can be seen, the average slopes of all four profiles are comparable inside and outside the ring, which is located at 80-120 AU. The average outer slope is 1.5, but the ring-only configuration drops more sharply at the edge of the ring (1.7) and flattens more at larger distances (1.2), thereby deviating more strongly from a single power law. This flattening is caused partly by the fact that the outermost semimajor axis bin is overpopulated, because it represents orbits from 380 AU to infinity. The two power-law profiles make a smoother turn-over.

The convergence of one such radial profile towards a steady state is illustrated in Fig. 10. This plot demonstrates the transition from a steep profile $r^{-\alpha}$ with $\alpha \ga 3$(corresponds to $\tau(a) \propto a^{-4}$) to a flattened one with $\alpha \approx 1.5$. At the beginning of the simulation, the optical depth grows, reflecting overproduction of the smallest $\alpha$-meteoroids, i.e. development of the primary maximum in the wavy size distibution (Fig. 6). Then the optical depth starts to decrease, due to ongoing collisional depletion of large bodies and therefore decreased production of collisional fragments at dust sizes. The figure also shows how the radial location of the peak optical depth is shifted outward due to higher rates of mass loss in the inner regions.


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{normal optical depth}
\psfrag{XAxi...
...]{distance to Vega [AU]}
\includegraphics[width=8.8cm]{4907f10.eps} \end{figure} Figure 10: Convergence of the radial profile of the normal optical depth for a rocky disk with an outer slope of -4.
Open with DEXTER

The resulting steady-state slopes of surface density that we find, between 1 and 2, have to be compared to observations. For large grains in thermal equilibrium the temperature is proportional to r-1/2. For wavelengths much larger than the peak flux wavelength for grains of 50-100 K (60-30  ${\rm\mu m}$), the emitted black-body energy is proportional to the temperature itself. In this case, the total surface brightness I of the thermal emission, which is a product of the surface density or the optical depth and the intensity in the observed spectral range, varies as $I \propto r^{-\alpha - 1/2}$. At wavelengths close to the maximum of the Planck function, the maximum intensity goes as the fifth power of the temperature, giving $I \propto r^{-\alpha - 5/2}$. Therefore, a radial profile of surface brightness with indexes in a broad range between about -2 and -4 may be consistent with our model, depending on the temperature distribution and the wavelength of the observations. Su et al. (2005) report that the observed surface brightness follows a power law with index $\approx $4 for wavelengths of 24 and 70  ${\rm\mu m}$, and $\approx $3.5 for 160  ${\rm\mu m}$. They suggest that small particles, $\approx $ ${\rm\mu m}$, which sweep out throught the disk with an ideal $\alpha = 1$ and, according to their calculations, with $I \propto r^{-3}$, are responsible for the observed flux. However, the latter contradicts our model: the contribution of grains of 2  ${\rm\mu m}$ is by two orders of magnitude less than that of grains of 10  ${\rm\mu m}$ (Fig. 6). Collisions of $\beta $-meteoroids with $\alpha$-meteoroids, which are neglected here, could slow the escape or produce more small fragments, thereby prolonging the effective residence time, but probably not efficiently enough.

So, we expect the main contribution to come from $\alpha$-meteoroids, i.e. 10  ${\rm\mu m}$ and larger grains. As mentioned above, this may be compatible with the observations. For a more conclusive comparison to the observational data, realistic calculations of grain temperatures and thermal fluxes need to be applied to our results, which is beyond the scope of this paper. There are many other physical mechanisms that could affect the radial profiles. For instance, the disk may not be in a steady state. A single catastrophic event in the past may have occurred, whose remnants produce the main fraction of the currently visible fragmented material. The primordial population of fragments from this event is expected to spread soon over a circumstellar ring, which resembles our initial situation. If the time elapsed after the event were too short to reach the steady state, a steeper profile would be expected (Fig. 10). The upper limit (from Fig. 10) for this time is of the order of 106 to 107 yr.

The radial profile may be steepened by a mechanism that continuously removes grains in orbits with higher eccentricities, for instance a planet interior or exterior to the belt of parent bodies. As the radial distribution for a dynamically cold disk with reduced eccentricities shows (Fig. 9), the absence of eccentric parent bodies, which could be explained by an inner planet, does not influence the final result. A simple cut-off criterion forbidding planet crossing orbits only for parent bodies would have no impact on the results, either, but a continuous removal of small fragments on orbits of high eccentricity and short pericentric distances, which exist (Fig. 7, 50 AU; Fig. 8), could steepen the profile. The influence of an outer planet could be similar.

Table 1: Masses of the modeled disks producing $0.5 M\ensuremath {_{\rm lunar}} $ of dust.

The total masses of the model disks are listed in Table 1 for different runs. The "rocky'' disks are made up of grains from 0.1  ${\rm\mu m}$ to 300 m, the "icy'' particles span a size range from 0.14  ${\rm\mu m}$ to 430 m, so that the same mass range is covered. The given masses and loss rates for ice are probably less realistic, as the real dust is assumed to be a mixture of silicates and carbonaceous material. As the mass-per-size-decade distribution (Fig. 11) shows, the derived total mass, which determines the absolute values of ordinates, is very sensitive to where the upper cut-off is set, especially when considering the wave pattern. The relative mass loss, $\dot M / M$, shows that all disks provide enough material to sustain the current debris production over at least 80 Myr up to 2.8 Gyr. A higher mass-loss rate of the icy disk stems from the fact that ice is easier to destroy than rock (smaller $Q\ensuremath{_{\rm D}} ^*$, smaller mcr, and therefore shorter collisional lifetimes). The existence of even larger parent bodies would further prolong the disk's half-life and increase the total mass, according to $M \propto
m\ensuremath{_{\rm max}} ^{0.5}$, if the mass distribution's power law still holds for that regime. To analyse the evolution of total mass with time, we chose a rocky disk with a ring-only distribution of parent bodies with a maximum radius of 5 km. Figure 12 depicts the resulting moderate mass loss over a period of 2 Gyr. From the mass loss rate of a system in a steady state, $\dot M = -c_M M^2$, where cM is a constant, we can deduce the dependence of mass on time:

 
$\displaystyle M(t) = M_0 / [1 + c_M M_0 \cdot (t - t_0)],$     (32)

where M0 = M(t0) (cf. Dominik & Decin 2003).

While Su et al. (2005) estimate the current mass loss due to blow-out of small grains as $6\times 10^{14}$ g/s or 3200 $M_\oplus$/Gyr, we conclude that mass-loss rates of 1 to 20 $M_\oplus$/Gyr (Table 1) are consistent with an observed amount of 0.5  $M\ensuremath{_{\rm lunar}} $ of submillimeter-sized dust. A linear extrapolation by 350 Myr back in time would lead to initital disk masses several orders of magnitude below the 330-3300 Jupiter masses given by Su et al. (2005). Their estimation is based on the short lifetime of unbound $\beta $-meteoroids, which they assume to be responsible for the observable flux. The lifetime of small $\alpha$-meteoroids on bound orbits is determined by the rate of destructive collisions and is much longer for optically thin disks. Given the observed amount of material and the longer lifetime, the necessary production rate and the total mass decrease to more plausible values.

  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{$\frac{\mbox{mass
density}}{\mbox{u...
...in radius [$\rm\mu m$ ]}
\includegraphics[width=8.8cm]{4907f11.eps} \end{figure} Figure 11: Mass distribution at 100 AU after 10 Myr for the same disk as in Fig. 4.
Open with DEXTER


  \begin{figure}
\par\psfrag{YAxisLabel}[c][c]{total mass [$M_\oplus$ ]}
\psfrag{...
...Label}[c][c]{time [Myr]}
\includegraphics[width=8.8cm]{4907f12.eps} \end{figure} Figure 12: Decay of the total mass of a rocky, initially ring-only disk (thick line) and its fits through Eq. (32) with t0 = 200 (dotted), 600 (thin solid), and 2000 Myr (dashed). The largest bodies in this simulation measured 5 km in radius.
Open with DEXTER

5 Discussion and conclusions

The subject of this paper is a "typical circumstellar debris disk''. Instead of including a large array of forces and effects in a sophisticated numerical code to get as realistic model as possible, we have chosen here to concentrate on a few fundamental physical effects and to develop a kind of a reference model. Our main interest was to identify several essential, generic properties of such a debris disk, driven by these effects.

We thus treat the disk as an idealized ensemble of dust particles, exposed to stellar gravity and direct radiation pressure and experiencing fragmenting collisions. Two other effects, often thought to be of primary importance - gas drag and Poynting-Robertson drag - intentionally, are not included, for the following reasons. In debris disks - as opposed to protoplanetary disks and possibly transitional disks - gas drag does not play a significant role, simply because there is not enough gas (Thébault & Augereau 2005). The Poynting-Roberston effect is usually not important either, since collisional timescales are typically much shorter than the Poynting-Robertson decay time (Wyatt 2005; Krivov et al. 2000).

We use a method that can be referred to as a "kinetic theory in orbital elements''. Solving a balance equation for material gain and loss, we obtain a three-dimensional distribution of mass, semimajor axes, and eccentricities, averaged over the other elements. This distribution can easily be transformed to "usual'' distributions, such as size distribution or radial profiles.

We have investigated what kind of size/mass and spatial distributions of the disk material can be expected from theory, and how these distributions might depend on the distributions of the parent bodies, parameters of the central star, and grain properties. With a number of tests, we have checked that the properties we found are generic. The model designed here may serve as a skeleton for the development of more complete models with a number of additional effects included: physical sources, drag mechanisms, more accurate impact physics with cratering collisions and restitution, etc.

An application is made to the disk around Vega. This disk is seen edge-on and appears nearly uniform, representing an ideal application of our model.

Our main findings are as follows:

1.
Size distribution. The radiation pressure strongly influences the grain size distribution in a debris disk. The cross section area of the grains, and therefore the disk brightness, are dominated by particles somewhat above the blowout limit, i.e. by smallest grains that can barely stay in bound orbits against the radiation pressure. For the Vega disk, this corresponds to a grain radius of roughly 10  ${\rm\mu m}$. The cross sectional density of smaller particles, sweeping through the disk in hyperbolic orbits, is by about two orders of magnitude lower.

Radiation pressure affects the distribution towards larger sizes as well. Overabundance or underabundance of smaller projectiles are repeatedly reflected at larger sizes, and thus, a kind of resonance comes into play, generating a wavy pattern in the size distribution. This might influence modeling of disk's spectral properties and disk mass estimates. The pattern is well understood (Campo Bagatin et al. 1994) as the result of the lower size cut-off produced by the blowout due to radiation pressure and was obtained by other models (Thébault et al. 2003), too. The reported peak at a grain size twice as large as the smallest bound grains is found as well as the lack of grains 20 to 100 times as large. However, the quantitative results differ, because the pattern's wavelength and amplitude strongly depend on the grains' mechanical properties and average impact velocities.

2.
Spatial distribution. The outer slope of the radial profile of the normal optical depth, ${\propto} r^{-\alpha}$, converges towards an upper limit of $\alpha \approx 1.5$. This would correspond to a power law index of the surface brightness profile between 2 and 4 if the thermal emission comes mainly from $\sim $10  ${\rm\mu m}$-sized grains, as suggested by the size distribution derived here. This may be consistent with observations of the Vega disk that reveal profiles with a surface brightness index of 3 to 4 (Su et al. 2005). Even if it is not, there may be a variety of other physical mechanisms not included in our model that affect the spatial distribution. The bright component of the disk may not have reached a steady state yet. For instance, this can happen in a disk where "large'' collisions between planetesimals have recently produced copious amounts of material that has not had enough time to be fully reprocessed by collisional grinding, so that the populations of solids at dust sizes are not "heated'' by radiation pressure yet. Alternatively, a planet - either interior or exterior to the location of parent bodies - could systematically remove grains in eccentric orbits, increasing the slope.

3.
Timescales and steady state. The size distribution converges towards a quasi-steady state on timescales depending on the collisional rates. The latter, in turn, depend on grain sizes and the distance from the star. For grains up to some hundred ${\rm\mu m}$ at 100 AU, about 1 Myr is needed to reach a steady state.

4.
The disk mass and mass loss. The total disk masses needed for the model disks to produce the observed amount of micron and submillimeter-sized dust do not exceed several Earth masses for an upper size limit of 300 m for rocky bodies. Even if the size distribution holds for the 2 orders of magnitude larger bodies (30 km), the total mass of 10 to 30 $M_\oplus$ will still be reasonable, although higher than the expected mass of 1 to 2 $M_\oplus$ of the Solar system's Edgeworth-Kuiper belt. Because the derived lifetimes of the rocky disks are of the order of 1 Gyr, and the smallest particles that are still in bound orbits are assumed to dominate the optical depth, an extrapolation 350 Myr back to the early times of Vega's disk would not lead to a disk as massive as the 330 to 3300 Jupiter masses concluded by Su et al. (2005). Thus, the results given here do not support their idea of a recent large collisonal event as the dust source, but do not rule it out either.

Acknowledgements
We wish to thank the anonymous referee for a thorough and helpful review. Miodrag Sremcevic is supported by the Cassini UVIS project.

References

 

  
Online Material

Appendix A: Calculation of the collisional probabilities and impact velocities

Consider first a single particle flying a certain distance through an ensemble of targets distributed at a given number density $\rho_{\rm t}$ with an interaction cross-section $\sigma$. The differential hit probability with respect to x, the distance flown, is given by $\sigma\rho_{\rm t}\:\ensuremath{{\rm d}} x$. If, instead of a single projectile, we consider a second ensemble of a given density $\rho_{\rm p}$, interacting with the target ensemble, we get the symmetric expression

 
$\displaystyle \ensuremath{{\rm d}} W = \sigma \rho_{\rm p}\rho_{\rm t}V \ensuremath{{\rm d}} x$     (A.1)

for the number of collisions, with V being the effective volume of interaction. The actual rate of collisions involving an impact velocity $v\ensuremath{_{\rm imp}} = \ensuremath{{\rm d}} x / \ensuremath{{\rm d}} t$ is then
 
$\displaystyle \frac{\ensuremath{{\rm d}} W}{\ensuremath{{\rm d}} t} = \sigma \rho_{\rm p}\rho_{\rm t}v\ensuremath{_{\rm imp}} V.$     (A.2)

In 2D, a single bound orbit is fully represented by a triplet of elements: the semimajor axis a, the eccentricity e , and the longitude of the pericenter (or: angular position of the apsidal line) $\varpi = \Omega + \omega$. Under the axisymmetric condition assumed here, the distribution over $\varpi$is uniform from 0 to $2\pi$. Then, collisions between two orbits $(a_{\rm t}, e_{\rm t}, \varpi_{\rm t})$ and $(a_{\rm p},
e_{\rm p}, \varpi_{\rm p})$ occur in one or two volumes of (infinitely) small extent (or: points) if the orbits are grazing or intersecting. The actual size of the volume is given by the two "thicknesses'' $d_{\rm t}$ and $d_{\rm p}$ of the ellipses, the crossing angle $\gamma$, and the vertical thickness or height h of the disk at distance r (Fig. A.1). This height depends on the disk's geometry and is given by $h = 2r\sin\varepsilon$ for a standard case with semi-opening angle $\varepsilon$.


  \begin{figure}
\includegraphics[width=8.8cm]{fiA1.eps}\end{figure} Figure A.1: Collision of two groups of particles in a small volume of interaction $V = d_{\rm t}d_{\rm p}h / \vert\sin\gamma\vert$, where h is the height perpendicular to the plane of projection, i.e. the disk's thickness at that point.

The linear number density $\hat\rho$ of projectiles and targets along their orbits is determined by their phase space density n together with the elementary phase space volume $\ensuremath{{\rm d}} {\vec p} = \ensuremath{{\rm d}} a\: \ensuremath{{\rm d}} e\: \ensuremath{{\rm d}} m$ and the Kepler equation, which gives the individual orbital velocity v and the orbital period P:

$\displaystyle \ensuremath{{\rm d}}\hat\rho = \frac{n}{vP}\:\ensuremath{{\rm d}} {\vec p}.$     (A.3)

Thus, the usual number densities at the point of collision and the interaction volume can be written as
$\displaystyle \ensuremath{{\rm d}}\rho_{\rm t}= \frac{n_{\rm t}\;\ensuremath{{\...
...\rm p}d_{\rm p}h},
\qquad
V = \frac{d_{\rm t}d_{\rm p}h}{\vert\sin\gamma\vert},$     (A.4)

containing phase space densities $n(a,e,\varpi) = n(a,e)/(2\pi)$. The involved orbital velocities, the orbital periods, and the crossing angle $\gamma$ can be expressed in terms of orbital elements:
  
                      $\displaystyle v\ensuremath{_{\rm p,t}} $ = $\displaystyle \sqrt{GM(1-\beta\ensuremath{_{\rm p,t}} )\left[\frac{2}{r} -
\frac{1}{a\ensuremath{_{\rm p,t}} }\right]},$ (A.5)
$\displaystyle P\ensuremath{_{\rm p,t}} $ = $\displaystyle 2\pi\sqrt{\frac{a\ensuremath{_{\rm p,t}} ^3}{GM(1-\beta\ensuremath{_{\rm p,t}} )}},$ (A.6)
$\displaystyle \cos\gamma$ = $\displaystyle \frac{1 + \left(\frac{1}{r}\frac{\partial r}{\partial\vartheta_{\...
...left(\frac{1}{r}\frac{\partial r}{\partial\vartheta_{\rm t}}\right)^2\right]}},$ (A.7)
$\displaystyle \vert\! \sin\gamma\vert$ = $\displaystyle \frac{\left\vert\frac{1}{r}\frac{\partial r}{\partial\vartheta_{\...
...(\frac{1}{r}\frac{\partial r}{\partial\vartheta_{\rm t}}\right)^2\right]}}\cdot$ (A.8)

The resulting variant of Eq. (A.2) contains the phase space densities n instead of number densities $\rho$:
 
                        $\displaystyle \ensuremath{{\rm d}} {\ensuremath{{\rm d}} W \over \ensuremath{{\rm d}} t}$ = $\displaystyle \frac{\sigma n(a_{\rm p},e_{\rm p},\varpi_{\rm p}) \: n(a_{\rm t}...
...p}\; \ensuremath{{\rm d}} {\vec p}_{\rm t}\; \ensuremath{{\rm d}}\varpi_{\rm t}$  
  = $\displaystyle \frac{\sigma n(a_{\rm p},e_{\rm p}) \: n(a_{\rm t},e_{\rm t}) \: ...
..._{\rm p}\; \ensuremath{{\rm d}} {\vec p}_{\rm t}\; \ensuremath{{\rm d}}\Lambda,$ (A.9)

where $\varphi(\Lambda)$ represents the distribution of pairs of colliders over $\Lambda$, which is uniform, as the system is axisymmetric. As all variables in Eqs. (A.5)-(A.8) depend on $\varpi_{\rm p}$ and $\varpi_{\rm t}$ through $\Lambda = \varpi _{\rm p}- \varpi _{\rm t}$ only, $\ensuremath{{\rm d}}\varpi_{\rm p}\;\ensuremath{{\rm d}}\varpi_{\rm t}/(4\pi^2)$ is transformed to $\ensuremath{{\rm d}} (\varpi_{\rm p}-\varpi_{\rm t})\;\ensuremath{{\rm d}} (\varpi_{\rm p}+\varpi_{\rm t})/(8\pi^2)$, and integrated once over $\varpi_{\rm p}+ \varpi_{\rm t}$ from 0 to $4\pi$, is giving $\ensuremath{{\rm d}}\Lambda/(2\pi)$ or $\varphi(\Lambda)\ensuremath{{\rm d}}\Lambda$.


  \begin{figure}
\includegraphics[width=8.8cm]{fiA2.eps}\end{figure} Figure A.2: Two elliptic orbits crossing in 2D: (a) two cases for fixed difference $\Lambda = \varpi _{\rm p}- \varpi _{\rm t}$, (b) two cases for fixed $r = r_{\rm t}= r_{\rm p}$.

After separating the number densities, the impact velocity, the cross section $\sigma$, and the differentials, the remaining term describes the geometric probability of a collision of two such orbits, but averaged over all variables not involved in the above calculations: the inclinations and the absolute orientation of (one of) the orbits. This probability is of unit one per volume and is very similar to the $\Delta $-integral in Krivov et al. (2005), but without averaging over the relative positions of both orbits - there is still a dependence on the distance to the star r at which the collision occurs:

 
$\displaystyle \Delta({\vec p}_{\rm p},{\vec p}_{\rm t}, \Lambda)
=
\frac{1}{v_{...
...rm p}P_{\rm t}2 r(\Lambda) \sin\varepsilon
\vert\sin\gamma (\Lambda)\vert}\cdot$     (A.10)

Except for the special case of grazing orbits with $\gamma = 0$, two orbits in the same plane always cross at two points (Fig. A.2), or given the target's true anomaly, there are two ways in which the projectile can cross the target's orbit: moving inward or moving outward. Accordingly, both signs have to be treated for the derivative $\partial r /
\partial\vartheta_{\rm p}$, giving $\gamma{\ensuremath{^+} }$ and $\gamma{\ensuremath{^-} }$. From the absolute values of $v_{\rm p}$, $v_{\rm t}$ and $\gamma$ we derive the relative velocity
 
$\displaystyle v\ensuremath{_{\rm imp}} ({\vec p}_{\rm p},{\vec p}_{\rm t}, \Lam...
...
=
\sqrt{ v_{\rm p}^2 + v_{\rm t}^2 - 2v_{\rm p}v_{\rm t}\cos\gamma(\Lambda) }.$     (A.11)

We now dicsuss 3D corrections necessary to account for (small) orbital inclinations. As the gain and loss terms (10) and (11) contain integrals of $\Delta({\vec p}_{\rm p},{\vec p}_{\rm t}, \Lambda)$ and $v\ensuremath{_{\rm imp}} ({\vec p}_{\rm p},{\vec p}_{\rm t}, \Lambda)$over $\Lambda$, we calculated averages
  
                   $\displaystyle \Delta ({\vec p}_{\rm p},{\vec p}_{\rm t})$ = $\displaystyle {1 \over 2\pi}
\int_0^{2\pi}\limits
\Delta ({\vec p_{\rm p}},{\vec p_{\rm t}};\;\Lambda) \:\ensuremath{{\rm d}}\Lambda,$ (A.12)
$\displaystyle \ensuremath{\overline{v_{\rm imp}}} ({\vec p_{\rm p}},{\vec p_{\rm t}})$ = $\displaystyle {
\int_0^{2\pi}\limits
\ensuremath{\overline{v_{\rm imp}}} ({\vec...
...a ({\vec p_{\rm p}},{\vec p_{\rm t}};\;\Lambda) \:\ensuremath{{\rm d}}\Lambda
}$ (A.13)

and compared the results with those yielded by strict Monte-Carlo integrations in 3D (cf. Krivov et al. 2005). We found that Eq. (A.10) does not show any noticeable dependence on the vertical extension of the disk beyond the term $\sin\varepsilon$ in the denominator. The impact velocity turned out to be more sensitive to the possible range of inclinations, and Eq. (A.11) is corrected to
 
$\displaystyle v\ensuremath{_{\rm imp}} ({\vec p}_{\rm p},{\vec p}_{\rm t}, \Lam...
...
= \sqrt{v_{\rm p}^2 + v_{\rm t}^2 - 2v_{\rm p}v_{\rm t}
\cos\overline{\gamma}}$     (A.14)

with
$\displaystyle \cos\overline{\gamma} = \frac{\cos\gamma}{\sqrt{1+(B\sin\varepsilon)^2}}\cdot$     (A.15)

The correction constant B was found empirically from the comparison with the Monte-Carlo integration: B = 2/3. Figure A.3 shows examples of both quantities, $\Delta ({\vec p}_{\rm p},{\vec p}_{\rm t})$ (Eq. (A.12)) and $\ensuremath{\overline{v_{\rm imp}}} ({\vec p}_{\rm p},{\vec p}_{\rm t})$ (Eq. (A.13)), with integrands given by Eqs. (A.10) and (A.14), as well as their comparison with the Monte-Carlo results.


  \begin{figure}
\psfrag{YAxisLabelDel}[bc][bc]{$\Delta(a_{\rm p},e_{\rm p},a_{\rm...
...th=8.8cm]{4907fA3a.eps} \includegraphics[width=8.8cm]{4907fA3b.eps} \end{figure} Figure A.3: The $\Delta $ integral ( top) and the impact velocities \ensuremath {\overline {v_{\rm imp}}} ( bottom) for different combinations of semimajor axes and eccentricities of two colliding particles with $(a_{\rm t},e_{\rm t})$ and $(a_{\rm p},e_{\rm p})$. The disk's full-opening angle was set to $2\varepsilon = 17^\circ $. For each set of orbits the "exact'' values (symbols), obtained with a time-consuming Monte Carlo evaluation, are compared to the 2D-based approximations, Eqs. (A.10) and (A.14) as solid lines. The semimajor axes are: $a_{\rm t}= 1.6$, $a_{\rm p}= 1.0$. The projectile's eccentricities are $e_{\rm p}= 0.2$and $e_{\rm p}= 0.5$.

At this point we are able to evaluate the gain and loss terms (10) and (11). For example, the gain term can be written as

$\displaystyle \left(
\frac{\ensuremath{{\rm d}} n}{\ensuremath{{\rm d}} t}
\rig...
...-} }v\ensuremath{_{\rm imp}} {\ensuremath{^-} }\Delta{\ensuremath{^-} }
\right]$      
$\displaystyle \times~\: n_{\rm p}n_{\rm t}\: \varphi(\Lambda) \: \ensuremath{{\...
...m d}} e\ensuremath{_{\rm p,t}}\;\ensuremath{{\rm d}} m\ensuremath{_{\rm p,t}} .$     (A.16)

To avoid some ambiguities and singularities during the calculation, it is convenient to switch from integration over $\Lambda$ to integration over $\vartheta_{\rm t}$. The Jacobian gives
 
$\displaystyle \frac{\ensuremath{{\rm d}}\Lambda}{\ensuremath{{\rm d}}\vartheta_...
...l\vartheta_{\rm t}}}{\frac{1}{r}\frac{\partial r}
{\partial\vartheta_{\rm p}}},$     (A.17)

fitting quite well to $\sin\gamma$ (Eq. (A.8)), because their numerators cancel out. The resulting probability can be written as
 
$\displaystyle {\Delta ({\vec p}_{\rm p},{\vec p}_{\rm t};\; \vartheta_{\rm t})}...
...) \over a_{\rm p}} - {p_{\rm p}\over
r(\vartheta_{\rm t})}\right)}
\right]^{-1}$     (A.18)

and is no longer symmetric with respect to the target and the projectile. In the code, the Jacobian was not applied in its differential form, but used together with discretized steps of $\vartheta_{\rm t}$ as a quotient of differences. Near singularities, that still appear, the program slightly shifts the eccentricity of both colliders to mimic the averaging over a whole bin.

One gets the same result as in Eq. (A.18) by starting from Eq. (7), reformulating the Dirac functions to depend on one collider's longitude of the pericenter, $\varpi_{\rm t}$, as well as $\vartheta_{\rm t}$ and $\vartheta_{\rm p}$. The Jacobian, transforming in 2D from the Cartesian case ${\vec
r}=(x,y)$ to $r=r(\vartheta_{\rm t})=r(\vartheta_{\rm p})$ and $\Lambda = \vartheta_{\rm t}^* - \vartheta_{\rm p}^*$, is given by

                        J-1 = $\displaystyle \left\vert\frac{\partial(x_{\rm t}- x_{\rm p}, y_{\rm t}- y_{\rm ...
...rm t}- \vartheta_{\rm t}^*,\vartheta_{\rm p}-
\vartheta_{\rm p}^*)}
\right\vert$  
  = $\displaystyle \frac{p_{\rm t}p_{\rm p}\:\vert e_{\rm t}e_{\rm p}\sin\Lambda + e...
...{(1+e_{\rm t}\cos\vartheta_{\rm t})^2(1+e_{\rm p}\cos\vartheta_{\rm p})^2}\cdot$ (A.19)

The distributions in Eq. (7) are
           $\displaystyle \varphi(\vartheta_{\rm t})$ = $\displaystyle \frac{(1-e_{\rm t}^2)^{3/2}}{(1+e_{\rm t}\cos\vartheta_{\rm t})^2},$  
$\displaystyle \varphi(\vartheta_{\rm p})$ = $\displaystyle \frac{(1-e_{\rm p}^2)^{3/2}}{(1+e_{\rm p}\cos\vartheta_{\rm p})^2},$  
$\displaystyle \varphi(\varpi_{\rm t})$ = $\displaystyle {1 \over 2\pi}\cdot$ (A.20)

Finally, the 2D probability
$\displaystyle \Delta\ensuremath{_{\rm 2D}} =
\int_{\varpi_{\rm t},\vartheta_{\r...
...; \ensuremath{{\rm d}}\vartheta_{\rm t}\; \ensuremath{{\rm d}}\vartheta_{\rm p}$     (A.21)

leads again to Eq. (A.18), after division by $2r\sin\varepsilon$ to account for 3D.



Copyright ESO 2006