Online material

Appendix A: Numerical techniques

Our model is computationally very demanding, primarily because of the large number of calculations needed for collisions, which scales with the number of phase-space bins cubed. To ensure that the model can produce results in a reasonable amount of time, an effort was made to reduce the amount of calculations needed per run. To achieve this, we used an integration technique that allows large time-step sizes, which reduces the number of steps needed per run (Appendix A.1). Furthermore, we calculated as many time-independent numerical factors as possible before the actual integration, so these factors do not have to be determined at every time step (Appendix A.2).

Appendix A.1: Integration method

After discretization, the distribution function n(m,k,t) is replaced by the vector of population levels n, and Eq. (33) can be treated as the system of linear equations (A.1)where A is a matrix of coefficients, and b is a constant vector containing the artificial source terms that replenish dust in the parent belt¹⁶. This system of equations suffers from stiffness: the population levels of some bins change very rapidly compared to others mainly due to large differences in collisional timescales, and the time-step size is therefore determined by the stability of the integration method rather than by its accuracy. When using standard explicit integrators, this leads to an impractically small step size that prevents long integrations.

Following Löhne (2008), we resolve the stiffness by writing the differentials as (A.2)where A^′ is a diagonal matrix that contains only the diagonal elements of A, while the terms marked “const” contain the off-diagonal parts of A, as well as the constant terms b. Our integration scheme for the jth component of n, for timestep m + 1, reads as(A.3)where Δt denotes the time-step size¹⁷. This scheme is known as the exponential Euler method (see Hochbruck & Ostermann 2010 for an introductory review). It is suitable for time integration of semi-linear problems, which consist of a stiff linear part and a non-stiff non-linear part. In short, the strategy is to solve the linear part exactly and to approximate the non-linear part using an explicit integration scheme.

The time-step size Δt is determined dynamically from the condition that population levels can never become negative. We use a scheme similar to that of Krivov et al. (2005), but adapted for the exponential Euler method.

Appendix A.2: Precalculation

Only considering collisions (i.e., ignoring diffusion terms), the evolution of the jth component of n can be written as (Löhne 2008) (A.4)Here, t and p are the target and projectile bin indices, respectively, and B_jtp is an entry in the time-independent tensor of collisional coefficients B. Specifically, coefficient B_jtp is the rate at which the population level of bin j changes, per particle in bin t, per particle in bin p, combining all considered relative orbit orientations. If j = t or j = p, this is the loss rate of the target or projectile bin, respectively, caused by collisional destruction (assuming that the mass grid resolution is high enough that fragments do not end up in their parent bin). Otherwise, it is a rate at which fragments are created.

For the phase-space grids we use, the entire tensor B is too large to be stored in memory. However, it is very sparse, because (1) not all orbits that are part of the phase-space grid have mutual overlap; (2) overlap may occur for a limited range in relative orbit orientation; (3) impact velocities are not always high enough to cause catastrophic collisions; and (4) only a fraction of all possible masses are created as fragments.

By only storing the non-zero entries of B, it becomes possible to keep it in memory. We store the source and sink terms separately. For the source terms, which are by far the most memory intensive, each non-zero entry requires (1) the index of the target bin; (2) the index of the projectile bin; (3) the index of the bin fragment bin; and (4) the rate at which particles are created in the fragment bin, per particle in the target bin, per particle in the projectile bin. By further manual compression, only (3) and (4) need to be stored for each entry separately. The target bin index (1) only needs to be stored once for each possible target bin, along with the number of possible projectile bins for that target bin. Then, the projectile bin index (2) only needs to be stored once for each possible collisional pair of bins, along with the number of possible fragment bins for that pair. A similar scheme is used for the sink terms.

The compressed version of B is still very large (several gigabytes for the runs presented in this work). A disadvantage of the precalculation technique is therefore that the size of the grid that can be used is restricted by the available amount of memory, whereas a code that recalculates (parts of) B at each time step is only limited by CPU power. A cubic dependence of the non-zero entries of B on the number of mass bins (as opposed to quadratic dependencies on the orbital element variables) motivates the choice of a relatively small mass grid, representing only a part of the collisional cascade.

Appendix B: Model verification

We performed several tests to verify our numerical debris disk model. The P–R drag and sublimation modules of the code were tested by comparing the behavior of the model with independent numerical solutions of the equations of motion and sublimation for individual particles. For this purpose, the collisional module of the code was switched off, and a single bin was filled as the initial setup of the model, corresponding to the orbital elements of the particle. For all these tests, the resulting evolution of the dust distribution (not shown here) matches that of the independent solution. The accuracy of the predictions is limited by numerical diffusion and becomes better with higher resolution (i.e., larger phase-space grids).

To test the collisional module of our code, we benchmarked the model against the code of Krivov et al. (2006), by replicating one of their model runs for the debris disk of Vega as accurately as possible. Since the code of Krivov et al. (2006) uses semi-major axes as the distance dimension of the phase-space grid, the benchmark runs were performed with a version of our code that also uses the semi-major axis (as opposed to periastron distance, used in the rest of this paper). These runs do not include the effects of P–R drag or sublimation, so they can be used to separately test the collisional module of our code by switching off P–R drag and sublimation. Of the various runs presented by Krivov et al. (2006), the specific one that was reproduced is characterized by an initial optical depth profile in the outer disk (beyond 120 AU) that scales as τ_geo(r) ∝ r^-4, an initial eccentricity distribution between 0 and 0.375, and material properties for “rocky” grains. We refer the reader to Krivov et al. (2006) for the specific values used for parameters describing the phase-space grid, the initial setup of the disk, material properties, etc.

	Fig. B.1 Evolution of the radial geometrical optical depth profile of the benchmark run, to be compared with Fig. 10 of Krivov et al. (2006).
Open with DEXTER

	Fig. B.2 Evolution of the size distribution at r = 100 AU of the benchmark run, to be compared with Fig. 6 of Krivov et al. (2006).
Open with DEXTER

The evolution of the radial and size distribution predicted by the benchmark run are shown in Figs. B.1 and B.2, respectively. The corresponding results of Krivov et al. (2006) are their Figs. 10 and 6, respectively. Comparison of the results reveals good agreement between the two codes, and the remaining discrepancies can be accounted for. Relative to the benchmark, our model predicts (1) lower unbound particle populations and (2) shorter evolutionary timescales. Point (1) is to be expected since the unbound particle densities of Krivov et al. (2006) are computed using the product of their production rate and their disk crossing time, rather than from Eqs. (C.6) (Löhne, priv. comm.). We attribute point (2) to a mislabeling of the time steps in Figs. 6 and 10 of Krivov et al. (2006). Further discrepancies can be due to small differences in the input parameters and numerical techniques used.

Appendix C: Post-processing of model output

The output of a model run are the population levels n of all bins in the three-dimensional phase space of particle mass and orbital elements, as well as their change rates ṅ (which should equal zero once steady state is reached, except for bins corresponding to unbound orbits). While this data is useful for analyzing the orbital characteristics of different classes of particles in the debris disk, it is often more convenient to know about the state of the disk as a function of distance from the star. This is essential, for example, if the results of the model are to be compared with observations. Here, we describe the processing steps that are applied to the model output to derive the quantities used in Sect. 3.3.

Appendix C.1: Conversion from orbital elements to radial distance

To find the radial distribution of matter in the debris disk, the orbital element phase-space distribution function n(m,q,e) (dimension: [g^-1 cm^-1]) needs to be converted to the configuration space distribution function , which denotes the vertically averaged number density of particles with masses [m,m + dm] at distance r (dimension: [g^-1 cm^-3]). This problem was first solved by Haug (1958) for a rotationally symmetric ensemble of particles on Keplerian orbits. Here, we give a brief derivation under the additional assumption that the distribution of inclinations is independent of the distribution of the other orbital elements.

Consider an individual particle on a bound Keplerian orbit that spends dt time to cross a radial annulus with width dr at distance r from the star. The contribution of this particle to is (C.1)where P is the particle’s orbital period, and h = 2r sin ε is the disk height, where ε is the semi-opening angle of the disk. The explicit factor 2 in the numerator accounts for the fact that the particle passes through this radial annulus twice during each orbit. In terms of orbital elements q and e, the orbital period P, accounting for direct radiation pressure, is given by (C.2)The radial velocity ṙ of the particle is given by (C.3)Combining Eqs. (C.1)–(C.3), and considering all particles on bound orbits, gives (cf. Krivov et al. 2005) (C.4)with the integration domain (C.5)The contribution of unbound particles to is calculated using their production rate ṅ(m,q,e) and the radial velocity with which they leave the system (Löhne, priv. comm.). Assuming all unbound particles are created at the periastron of their orbits, their vertically-averaged particle number density is given by with the integration domain (C.8)Negative eccentricities correspond to “anomalous” hyperbolic orbits (see Krivov et al. 2006).

In applying this theory to the raw data, we replace the integrals in Eqs. (C.4) and (C.7) with sums; replace n(m,q,e) and ṅ(m,q,e) with n and ṅ, respectively; and evaluate the resulting equations at discrete points of r. For each bin, we sample the phase space it represents using a Monte Carlo method.

Appendix C.2: Derived quantities

For radial dust distribution profiles, we use the vertical geometrical optical depth τ_geo, defined as the surface density of collective (i.e., combining particles of all sizes) cross-section. It is computed from as To characterize size distributions, we use the quantity A(s,r), defined as the cross-section density per base-10 logarithmic unit of size. It is given by This quantity allows for an easy comparison between the relative contributions of particles with different sizes to the total cross-section of the disk. In the size distributions plots (Figs. 5 and 6), a horizontal line means particles of all sizes contribute equally to the cross-section, and equal areas under the curve correspond to equal contributions to the total cross-section.

© ESO, 2014