© ESO, 2017

1. Introduction

Observed dust in debris discs is produced in collisions amongst orbiting planetesimals. Resolved images at submillimeter wavelengths, which trace large grains, suggest that the dust parent bodies are arranged in narrow belts, similar to the classical Kuiper belt in the solar system. These narrow planetesimal belts are also evident in scattered light images in the optical and near-infrared. The small dust grains visible at these wavelengths are most abundant at the same locations where their parent bodies reside.

In some of the debris discs, these narrow planetesimal belts appear eccentric and show a global offset between the belt centre and the star. The disc in the Fomalhaut A system (Kalas et al. 2005, 2013) is perhaps the most prominent example. Another example is the HD 202628 disc (Krist et al. 2012). Many other discs, such as HD 32297 (Kalas 2005), HD 61005 (also known as the Moth, Hines et al. 2007), and HD 15115 (the Blue Needle, Kalas et al. 2007), exhibit global asymmetries between the two wings. It is possible that these asymmetries also derive from the offsets in the underlying belts, which may not be seen because of the edge-on orientation of the discs and/or insufficient spatial resolution of the submillimeter facilities. Belt offsets can naturally be explained by as yet undiscovered planets in eccentric orbits interior to (Lee & Chiang 2016; Esposito et al. 2016) or substellar companions exterior to the belts (Thébault et al. 2010; Thébault 2012; Nesvold et al. 2016). Alternatively, the wing asymmetries may also be caused by recent giant collisions (e.g. Kral et al. 2015; Olofsson et al. 2016) or displacement of the dust by the surrounding interstellar gas (Debes et al. 2009) or dust (Artymowicz & Clampin 1997).

Interpretation of asymmetries in discs in terms of potential perturbers requires models to predict how exactly such planets would shape the distribution of the disc material. Whereas the influence of planets on the parent belts is easily understood with the Laplace–Lagrange secular perturbation theory, the task gets more complicated for small dust grains. These dominate the cross section and thus also the observable appearance of extrasolar debris discs. However, these are not direct tracers of the underlying distribution of parent bodies from which they are produced because they are subject to an additional array of forces and effects, including collisional production and removal, radiation pressure, and drag forces (e.g. Wyatt et al. 1999).

Previous work (Stark & Kuchner 2008, 2009; Kuchner & Stark 2010; Thébault et al. 2012, 2014; Vitense et al. 2012; Kral et al. 2013; Nesvold et al. 2013; Vitense et al. 2014; Kral et al. 2015; Nesvold & Kuchner 2015b,a; Lee & Chiang 2016; Esposito et al. 2016, among others) extended purely gravitational models of planet-disc interactions by including these effects and forces acting on dust grains. In this paper, we tackle the problem with a novel approach that is based on modelling of the evolution of the phase-space distributions of the material, rather than N-body integrations as was done previously. We show that a combination of grain-grain collisions and Poynting–Robertson (PR) drag with the gravitational perturbations by massive bodies in the system creates inseparable size-spatial distributions of solids. These are used to explore the distributions of dust grain sizes in different disc locations and expected observable signatures in discs.

We start with a discussion of secular perturbations by planets as a potential cause for the narrow eccentric belts in Sect. 2. In Sect. 3 we introduce a new version of our collisional code ACE which can now treat azimuthal asymmetry. Sections 4 and 5 present the outcomes of simulations of fiducial discs in terms of dust distributions and observables, respectively. We summarize our findings in Sect. 6.

2. On possible origins of the offsets

The presence of eccentric belts is evident from observed images, but the origin of these offsets is yet unclear. Mechanisms other than a planet in eccentric orbit in the cavity of the disc are conceivable. For instance, Shannon et al. (2014) have proposed that the eccentricity of the belt around Fomalhaut A was set by dynamical interactions with the other two companions of this triple system, Fomalhaut B and C. Yet the planetary scenario is considered the most generic, and we now address it in more detail.

2.1. Secular perturbations

Where both short-period perturbations and resonances are unimportant, secular Laplace–Lagrange theory provides the means to compute the perturbing influence of planets or substellar companions (see e.g. Murray & Dermott 2000). Differences between mean anomalies of perturber and perturbee are assumed to be random in that approximation. While no energy is exchanged, orbital eccentricities and orientations of the orbits change. In a space spanned radially by eccentricity e and azimuthally by longitude of periapse ϖ (Fig. 1), secular perturbation makes the eccentricity vector (ecosϖ,esinϖ) precess uniformly along a circle with a radius called proper eccentricity, e_p, centred around a forced eccentricity e_f (Hirayama 1918). The closer the combination of e_p and e_f gets to unity, the higher are the deviations from perfect circles (e.g. Beust et al. 2014). The forced eccentricity vector is aligned parallel to the eccentricity vector of the planet, where its absolute value is (Murray & Dermott 2000): (1)where e_planet is the absolute eccentricity of the orbit of the planet. The bs are Laplace coefficients, which only depend on the ratio of semi-major axes of the (interior) planet and a perturbed belt object, α ≡ a_planet/a_b< 1. A given forced eccentricity can be caused by a nearby planet of the same eccentricity or a closer-in planet of higher eccentricity. Figure 2 depicts the corresponding planetary orbital eccentricity as a function of its relative distance to the perturbed belt. Although better approximations exist for perturber eccentricities e_planet> 0.2 (see e.g. Mustill & Wyatt 2009, and references therein), we use Laplace–Lagrange theory for the broad analysis in this section.

While the mass of the perturber does not influence e_f (as long as M_planet ≪ M_∗), this mass determines the timescale on which this precession occurs. Based on the (approximate) angular precession frequency for a single perturber (Murray & Dermott 2000), (2)the time for a full precession cycle can be estimated from (3)where is the orbital period in the belt. The additional factor 1 − β with β ≡ F_pr/F_grav, the ratio of the forces due to radiation pressure and gravitational pull, accounts for the reduction of effective stellar mass because of radiation pressure.

An important criterion that the dynamical evolution must fulfil is to allow for narrow belts. We consider two principle possibilities (Fig. 1), which is a simplification of the scheme of four classes discussed by Thilliez & Maddison (2015). One is that planetesimals in the belt have low proper eccentricities, e_p. In equilibrium, over timescales longer than the precession period, the width of the belt is then set by these proper eccentricities. If they are low, the belt remains narrow at all times (Fig. 1 left). However, this raises the question of how the orbits of the parent bodies came close to their forced value in the e-ϖ plane. It would be more natural to assume a second possibility, in which planetesimals were born in nearly circular orbits, before the planetary perturber emerged in the protoplanetary disc. If so, we must require the precession circle to cross the centre (i.e. the point e = 0) and the currently observed average belt eccentricity e = e_b. The minimum time for the belt eccentricity to precess to e_b is a fraction, (4)of a full cycle (see Fig. 1 right), potentially increased by if N precession cycles have already passed. However, precession also smears out the orbits around the forced eccentricity in a wide circle if the belt had no offset prior to the perturbation. The resulting spread in e–ϖ would be twice as large as e_f itself, leading to wide instead of narrow belts (Beust et al. 2014). Thus, for this scenario to produce a narrow belt, the differential precession timescales must be longer than the time since perturbation started, to prevent this smearing out. We can define such a differential timescale as (5)where Δe_b is the spread in eccentricities and ΔA denotes the spread in linear change rates in eccentricity space as (6)Inserting (6) into (5) and expanding the term in brackets in a series in α, we obtain (7)or with Eq. (3), (8)because the coefficients in the series approach a value of 4. The spread in eccentricities (Δe_b) remains lower than the average eccentricity (e_b) and the belt remains narrow as long as ΔT_prec>T_prec. We find (9)and (10)This constraint on the belt width is mild for N = 0 and low values of α, but is strong for α close to 1 and N> 0. That is, it is more likely that an observed narrow belt is still in its first precession cycle. Assume, for example, that the belt is distant from the planet, α ≪ 1, and that the observed belt eccentricity is approximately the same as the forced eccentricity. We find Δa_b/a_b< 1 / 4 for N = 0 and Δa_b/a_b ≲ 1 / 29 for N = 1. While the first belt can be broad, the second belt needs to be narrow.

2.2. Constraints on perturbing planets

Not excluding either of the two possibilities described in Sect. 2.1, we now briefly discuss what both would mean for the unseen perturbing planet. In the modelling described in the rest of the paper, we use the mean eccentricity of the parent belt, e_b, as a key paremeter. However, its interpretation in these two cases is different. In the low-e_p scenario, e_b is equal to the forced eccentricity e_f (Fig. 1 left). In that case, the apsidal line of the orbit of the planet is aligned with the major axis of the belt. In the slow-precession scenario, e_b is not equal to e_f. Instead, it represents the instantaneous value of the complex eccentricity e (Fig. 1 right). In that case, the planetary orbit is misaligned with the major axis of the belt (Beust et al. 2014).

3. Collisional model

The number of particles in debris discs is orders of magnitude beyond the scope of pure N-body simulations. Hence, statistical representations for particle distributions and/or collisions are used. Collision rates and outcomes are calculated for whole groups of similar particles, called super-particles, bins, tracers, or streamlines. A major difficulty common to all approaches is the sampling; the number of groups needs to be high enough to properly represent the modelled distribution and low enough to be computationally tractable. There are two main approaches to this grouping: (A) time-resolved and (B) orbit-averaged.

In approach A, particles in close spatial proximity and with similar velocity vectors are grouped into so-called super-particles (Grigorieva et al. 2007), which can be viewed as more or less coherent clouds of particles that move in parallel. When two clouds collide, collision rates among individual particles are calculated based on the local particle-in-a-box principle. The collision cross section – or the volume of interaction – of the super-particles can either be defined as a co-moving sphere (e.g. Grigorieva et al. 2007) or a (revolving) grid element in polar coordinates (e.g. Kral et al. 2013). Smaller interaction volumes reduce the rate of collisions per super-particle, but increase the rate of individual collisions per super-collision. Smaller super-particles allow for higher spatial (and temporal) resolution, but require more super-particles to reduce noise artefacts. The biggest advantage of this approach is the ability to model short-term effects, such as collisional avalanches (Grigorieva et al. 2007), major break-ups of planetesimals (Kral et al. 2015), or close stellar flybys (Nesvold et al. 2017). Owing to the underlying N-body integration, additional forces are easily implemented in these codes, including radiation pressure, drag-induced spiralling, resonant capture, and scattering. Collisional grooming, the algorithm presented by Stark & Kuchner (2009), can be considered a variant in which collisions do not create new dust. Once released at a (constant) production rate, grains stream along their trajectories, while their number densities are gradually reduced in collisions as their trajectories cross others (or themselves). The model settles towards an equilibrium. Levison et al. (2012) have introduced a further example for approach A.

In approach B, which our code ACE follows, particles are grouped according to their orbits. Instead of local clouds, each group populates a given ellipse (or hyperbola) with particle density uniform accross mean anomalies, i.e. uniform in time. Discrete orbits are fixed throughout the simulations, parameterized either by orbital elements (ACE), or again, by location and velocity vector (Thébault et al. 2003). This orbit-averaging makes the models meaningful only on timescales longer than the orbital period. It excludes application to very dense discs with short collision timescales and to the short-term perturbations that can be analyzed with approach A. On the other hand, the phase space of orbital elements eases long-term computations with time steps of millions of years, i.e. orders of magnitude longer than orbital periods.

3.1. Phase space and master equation

While previous versions of ACE could only treat axisymmetric debris discs (e.g. Krivov et al. 2006, 2013; Reidemeister et al. 2011; Löhne et al. 2012), the azimuthal distribution is now allowed to be non-uniform. The phase space spans an additional dimension that covers the orientations of object orbits, parameterized by longitude of periapsis ϖ ≡ Ω + ω, i.e. the sum of the longitude of the ascending node (Ω) and the argument of periapsis (ω). The vertical dimension is still averaged over (cf. Krivov et al. 2006).

In total, the discretized distribution of material is represented by bins with four dimensions: (1) object masses m; (2) orbital pericentres q; (3) eccentricities e; and (4) longitudes of periapse ϖ. Collisions among pair-wise combinations of bins are possible at up to two distinct points defined by the qs, es, and ϖs of the colliders. Collision velocities, rates, and outcomes follow directly (Krivov et al. 2006).

The discretized master equation to be integrated over time reads (11)where n_i is the number (or mass) of objects in the bin specified by multi-index i ≡ (i_m,i_q,i_e,i_ϖ). Coefficients G_ijk denote the gain, that is the specific rates at which objects of type i are formed in collisions among objects of types j and k; L_ij denote the loss, that is the specific rate at which objects of type i are removed in collisions with j. The ACE code models drag forces by advection from grid cell to grid cell (since Reidemeister et al. 2011). The coefficients T_ij hence denote the transport to cell i from (neighbouring) cell j. For example, PR drag reduces q in bin j at a rate . Given a bin width Δq_j, the contents of bin j are moved towards i at a rate .

With collision physics depending on the orientations of the orbits only through the difference ϖ_i − ϖ_j, the dependence of G_ijk and L_ij on i_ϖ, j_ϖ, and k_ϖ is only through pair-wise differences i_ϖ − j_ϖ, etc. The relation G_ijk = G_i′j′k′ for i′ = (i_m,i_q,i_e,i_ϖ + o), j′ = (j_m,j_q,j_e,j_ϖ + o), and k′ = (k_m,k_q,k_e,k_ϖ + o) with o ∈ N can be employed to speed up the calculations; instead of looping colliders over j_ϖ and k_ϖ, k_ϖ − j_ϖ and j_ϖ are used, making the innermost loop over j_ϖ trivial.

3.2. Collision outcomes

Depending on the masses of the colliders and the impact energy, we consider three main outcomes. First, disruption and dispersal occurs if the energy suffices to overcome both the material strength and the combined gravitational potential of the colliders. For that threshold specific energy we follow Löhne et al. (2012) and assume the size dependence described by Benz & Asphaug (1999) together with the modification from Stewart & Leinhardt (2009), i.e. (12)where the two terms in brackets on the right-hand-side represent (1) shock disruption in strength regime and (2) in gravity regime, scaled by impact velocity v_imp. By we denote the equivalent radius of a sphere with the combined volume of the colliders. The last line in Eq. (12) approximates the specific energy required to overcome self-gravity. It is important only for radii s ≳ 30 km and has not been taken into account in previous ACE versions.

Second, below that threshold, we call collisions cratering if the target retains at least half of its original mass, but half the impact energy is enough to disrupt the projectile. Gravitational accretion also falls into this category. Third, if both colliders stay intact, they are assumed to separate again, unless impact velocities are below 10 m/s. In all three cases, a cloud of smaller fragments is produced in addition to the remnants of the colliders. In the model, the total mass in escaping fragments is proportional to impact energy, (13)The fragment mass distribution is assumed to follow a power law with exponent η = − 11 / 6 up to a limiting mass, above which the power-law distribution would accumulate to exactly one further particle. The mass of the largest fragment is thus given by (14)where m_frag is the total mass in escaping fragments. The same prescription was used in previous ACE versions.

We treated the remnants from erosive collisions somewhat differently from new fragments. If their new combination of mass and orbit still has them in the same bin, the total mass in that bin is reduced by how much is transformed into fragments. The mass of the remnants is only added to the loss of one bin and the gain of another if the remnants move to a different bin. However, in what follows we simply refer to both fragments and remnants as fragments.

3.3. Fragment orbits

Momentum conservation requires that the cloud of fragments produced in a collision has the same centre of mass as the original colliders. Neglecting relative velocities of the fragments, i.e. assuming full energy dissipation, all share a common initial velocity. As soon as the cloud becomes optically thin, radiation (and wind) pressure segregate the fragments according to their size-dependent β ratios. In the spirit of Krivov et al. (2006) the following orbital semi-major axis a, semilatus rectum p, and eccentricity e can be derived for a fragment produced in a collision between a target (subscript t) and a projectile (subscript p) at a distance r, i.e. where and denote the effective masses of target and projectile, and the sum of the collider masses, weighed by the β ratio of the fragment. The signs of e and p are equal to that of 1 − β, corresponding to anomalous hyperbolae with e< 0 for β> 1.

For grains that are released from a macroscopic body with β_t = 0 on impact of a projectile with m_p ≪ m_t, Eq. (15) is reduced to (17)The blowout limit is reached where a → ∞, corresponding to (Kresák 1976) (18)That limit varies between (1 + e) / 2 and (1 − e) / 2 for a grain release at apastron (r = Q = a(1 + e)) and periastron (r = q = a(1 − e)), respectively. When launched from circular orbits, grains become unbound for β ≥ 1 / 2.

If the material distribution is non-axisymmetric, the relative orientations of fragment orbits with respect to those of the initial colliders need to be considered in the modelling. The true anomaly θ of a freshly released fragment can be determined from(19)and (20)At the mutual crossing points of two orbits, the difference of true anomalies equals the (negative) difference between their longitudes of periapse (21)For collisions between projectiles and targets, the true anomalies of the latter at these points are given by (cf. Krivov et al. 2006, Eq. (8)) (22)or (23)where

3.4. Model and grid parameters

The aim of our study is to identify the characteristic influence of collisions and drag on perturbed discs. In this first paper, we refrain from covering the complexity of the wide parameter space spanned by observed discs. Instead we present and discuss results for a small set of fiducial, typical debris systems. For the central star, we choose an A3 V main-sequence star, with mass, luminosity, and effective temperature adopted from the values reported for Fomalhaut, i.e. M = 1.92 M_Sun, L = 16.6 L_Sun, T = 8590 K (Mamajek 2012). The stellar photospheric emission is modelled with the nearest point in the PHOENIX/NextGen grid of models (Hauschildt et al. 1999). We model the material with a homogeneous mix (Bruggeman 1935) of astrosilicate (Draine 2003) and water ice (Li & Greenberg 1998) in equal volume fractions with a bulk density of 2.35 g cm^-3. The combination of this material with the high luminosity and radiation pressure of the early-type star sets a clear blowout limit at grain sizes of a few microns.

Asymmetries are most easily identified where narrow belts are resolved, which correspond to narrow volumes in the space of orbital elements. The relative radial width ΔQ_b/Q_b at the apastron Q_b = a_b(1 + e_b) of a belt can be estimated from (28)if a and e vary independently. A given radial HWHM of, for instance, δr/r = 10 % does not only limit Δa/a ≤ 10 %, but also Δe_b = e_p ≤ 0.1. For the parent bodies in our main simulation runs, we therefore assumed initial distributions that are confined to circular regions in the (ecosϖ,esinϖ) plane, centred on (e_b,0), with radii e_p = 0.1. The values used for e_b in the individual runs range from 0.0 to 0.4.

The eccentricity grid that we use spans [ 0.015,1.5 ] and is logarithmically spaced for e ≲ 0.4 and for e> 1, with factors Δe/e ≈ 0.25. In order to preserve accuracy for barely bound grains, the step size is additionally limited to Δe ≤ 0.1 for e< 1, resulting in nearly linear steps in the range 0.4 ≲ e< 1. The transitions between linear and logarithmic regimes are smooth. In total, we model the number of bins per unit logarithm of eccentricity with (29)The resulting number of eccentricity bins is 26. The grid of orbit orientations has 32 linear steps in ϖ, covering [ 0,2π). The resolution elements Δe × eΔϖ thus measure 0.25e × 0.2e for e ≲ 0.4 and 0.1 × 0.2e for 0.4 ≲ e< 1. Figure 3 illustrates the sub-grid of e and ϖ in a polar plot.

In the radial direction, a grid of logarithmically spaced periastron distances q ≡ a(1 − e) is used. A number of 85 bins spans a distance range from 30 au to 600 au, corresponding to relative bin widths of 0.036. In combination with an eccentricity spread e_p = 0.1, this width ensures that orbits from neighbouring q annuli cross. In all runs, the q-and-e sub-grid is filled initially such that a range of semi-major axes from 95 to 110 au is covered.

Finally, masses are binned logarithmically from a minimum grain radius of 0.26 ?m to a maximum of 49 km with factors of 12 in mass (or 2.3 in radius) between adjacent bins. For grains with radii s ≲ 30 ?m (corresponding to β ≳ 0.1), where radiation pressure is important, the spacing is refined to factors of 12^{1 / 4} in mass (or 1.23 in radius) with a smooth transition between these regimes. A number of 48 mass bins result. The total number of bins in the grid that are actually filled increases over time and then saturates at about 10⁶ in the runs presented here.

The initial size distribution is assumed to follow a power law n(s) ∝ s^νT_orbit/T₀ with ν = − 3.66, normalized to a total mass of two Earth masses. We follow Strubbe & Chiang (2006) and Lee & Chiang (2016) in scaling the initial abundances of grains with their orbital timescales T_orbit with respect to that of large grains (β = 0), T₀. This scaling is meant to account for the increased lifetimes of smaller grains as a result of their spending most of the time close to their apocentres, i.e. far from the star. We are free to choose this initial set-up to ease comparison with other work. That choice does not influence the dust distribution towards which the subsequent collisional evolution quickly converges.

4. Resulting size and radial distributions

In this section, we identify the impact that radiation pressure (Sect. 4.1), collisions (Sect. 4.2), and drag forces (Sect. 4.3) have on the spatial and size distribution of dust. Figure 4 shows the different timescales for these effects. Initially, grains populate the elliptic orbits induced by radiation pressure on short, orbital timescales. At an intermediate stage (t = 10⁷ yr), enough time has passed to bring grains with sizes s ≲ 10 cm to collisional equilibrium. Later, at an evolved stage (t = 2 × 10⁸ yr), grains with sizes s ≲ 100 ?m are in PR drag equilibrium, which means that the system is older than the time these grains need to spiral from the belt to the star. In what follows, we use these three stages to illustrate the different effects. A more detailed comparison of timescales can be found in Sect. 4.4.

4.1. Dynamical consequences

For the large parent bodies that are unaffected by radiation pressure, a non-vanishing mean eccentricity vector (ecosϖ,esinϖ) corresponds to a global offset. The typical rates and velocities at which collisions occur are dictated by the spread in orbital elements. Discs with different mean eccentricities still have similar erosion rates as long as this spread is comparable. The parallel evolution of total mass and dust mass in our simulation runs for different e_b, shown in Fig. 5, confirms this expectation. When looking at parent bodies and larger grains, a global eccentricity results in literally just an offset, where the relative widths of periastron and apastron side are equal.

The picture is different for smaller grains. The additional action of radiation pressure induces a lower blowout limit that depends on the birth location, which changes the size distribution in the parent belt. Blowout occurs for lower values of β for grains that are launched near periastron of a parent orbit. The excess velocity there helps these grains overcome the gravitational bond, increasing the maximum size of unbound grains. Vice versa, when released near the belt apastron, smaller grains can stay bound. Figure 6 illustrates this shift in blowout size with grain size distributions near the periastra and apastra of parent belts with different e_b. According to Eq. (18), the ratio between the blowout sizes at periastron and apastron is given by (1 + e_b) / (1 − e_b). This ratio reaches a value of 7:3 for e_b = 0.4, which is consistent with the ACE results.

In the left panel of Fig. 7, radial profiles of the normal optical depth, τ, are plotted for t ≈ 0 yr. At large radii, these profiles drop almost as τ ∝ r^{− 3 / 2}, the behaviour expected for discs in equilibrium (Krivov et al. 2006). This is by design because we employed the initial set-up of Strubbe & Chiang (2006), who derived that slope analytically. Wiggles in the radial profiles are artefacts of the narrow size distribution in the halo reaching the resolution limit of our mass grid.

The blowout limit then defines the typical sizes of the barely bound grains that form the halo. The different blowout limits induce different grain sizes, depending on the side of the halo. Barely bound grains are typically produced near their periastra, but spend most of their time near the apastra of their orbits, i.e. on the opposite side. Therefore the part of the halo that extends beyond the apastron side of the belt is produced at its periastron, and vice versa. The two upper rows in the left column of Fig. 3 show how larger grains have orbits that are still aligned azimuthally with the parent belt. Their periastra (as described by e and ϖ in that figure) remain close to the periastron of the belt. The small grains shown in the bottom row of Fig. 3 can only stay bound when they are created with their periastra at the apastron side of the belt. In consequence, the halo on the apastron side is formed by grains larger than those on the periastron side. Figure 8 shows this effect for a disc with e_b = 0.4. There, the β values of the bound halo grains on the periastron side approach 0.7 as distance increases. The halo on the apastron side is populated by larger grains, initially limited by radiation pressure to β< 0.3.

The widths of the size distributions on the two sides differ as illustrated in Fig. 9, where radial and azimuthal distributions for different grain sizes are compared. For all three grain sizes shown, the range of distances covered beyond the apastron side of the belt is wider because they are bound less in that direction. The regions populated by grains of different sizes overlap more strongly on the apastron side, and hence, grains of a wider range of sizes populate a given distance.

The main reason for the asymmetries seen in the halos in our simulations are the different conditions under which grains are launched from different sites along the parent belt. Initial orbital velocities differ along the belt, as do the resulting fragment orbits. Most of the effects described in this section can also be found in Lee & Chiang (2016). A prominent example is the bow or wing asymmetry seen in Fig. 9e as a dark arc on the right-hand side of that panel. Lee & Chiang have an equivalent arc in their Fig. 1 (and Fig. 7). It traces the outer boundary of the region populated by barely bound grains on the periastron side.

A similarly asymmetric halo is described by Kral et al. (2015) for dust produced in a giant breakup. However, the asymmetry there is mainly caused by the asymmetric distribution of launch sites, where the initial round of fragments are all produced from a single parent.

4.2. Effects of collisions

In collisional equilibrium a size distribution of infinite extent can be described well by a power law n(s) ∝ s^α, with α ~ − 3.5 (Dohnanyi 1969; Durda & Dermott 1997; O’Brien & Greenberg 2003; Wyatt et al. 2011; Pan & Schlichting 2012). However, notable ripples appear near physical breaks or cut-offs. Waves in the size distribution are induced for asteroids by the transition from strength to self-gravity (Durda & Dermott 1997) and for grains above the blowout limit by the radiation pressure cut-off (Campo Bagatin et al. 1994). At this lower size end, barely bound grains become overabundant because of a lack of smaller projectiles. In turn, this overabundance leads to a depletion of somewhat larger grains. Wavelengths and amplitudes of these waves are determined by impact energies relative to disruption thresholds (Krivov et al. 2006). More realistic impact physics, such as cratering collisions, quickly damp the waves towards larger grain sizes (Thébault et al. 2003; Thébault & Augereau 2007; Müller et al. 2010), when compared to simulations where only disruptive collisions are considered (Löhne et al. 2008).

This wave near the blowout limit is overlaid by an effect first described in Thébault & Wu (2008). For larger grains, the typical orbital eccentricities, and hence the typical relative velocities and collision timescales, are determined by that of their parent bodies. For smaller grains, radiation pressure is more important. An additional break in the size distribution occurs where the two effects are equal, i.e. where e from Eq. (16) equals e_b. For grains produced in a parent belt with average proper eccentricities e_p = 0.05, this break is expected near β = e_p/ (1 + e_p) ≈ 0.05, corresponding to grain sizes s ≈ 45 ?m in our set-up. The right panel of Fig. 6 shows both the break around this size and the depletion and blowout-induced waviness below.

The closer grain sizes get to the blowout limit, the further size distributions near the parent-belt apastra and periastra deviate from one another. The peak just above the blowout limit is higher at the parent apastra than at the parent periastra. This difference increases with increasing e_b, reaching an order of magnitude for e_b = 0.4 (Fig. 6). This can be understood because the two sides are coupled; the larger grains coming from the periastra of the belts suffer from collisions with the smaller grains coming from the apastra as their orbits cross.

For grains larger than around 20 ?m, the situation seems reversed. In collisional equilibrium, grains in this size range contribute more to the optical depth at the belt periastra. This effect is mainly caused by these grains being spread out more widely on the apastron side, leading to lower densities there. However, part of this asymmetry extends to grains sizes where radiation pressure and the resulting radial spread are neglible. As a consequence, the size distributions are shallower overall for grains between 20 ?m and 1 mm, which translates to shallower spectral energy distributions (SEDs) in the corresponding range of wavelengths (Draine 2006). When combined with the increased abundance of barely bound grains at the apastra of discs with higher e_b, the effective grain sizes shift towards smaller radii. This would imply higher temperatures near the apastra, reducing the brightness asymmetry due to pericentre glow at shorter wavelengths. The detailed effects on the observable SEDs and images are discussed in Sect. 5.

In Fig. 3 the most notable difference between the initial stage and intermediate stage is the filling of regions that could not be reached initially. In particular, grains with s = 4.3 ?m (β = 0.5) appear on bound orbits (e< 1) with periastra aligned with the periastron of the belt. These grains cannot just be produced in the parent belt. Drag cannot be responsible either because it cannot alter ϖ and it cannot re-bind unbound grains in significant numbers, which is why this process is not modelled in ACE. Thus, these grains must stem from medium-sized grains. An illustration of this process is given in Fig. 10, where the orbits of three particle types actually represented in the ACE grid are shown. The target (s = 12 ?m, β = 0.17, q = 98 au, e = 0.22, ϖ = 23°) and the projectile (s = 5.4 ?m, β = 0.4, q = 85 au, e = 0.75, ϖ = − 56°) are produced at different locations in the parent belt. At one of their mutual collision points, they produce a barely bound fragment that is aligned with the parent belt (s = 4.3 ?m, β = 0.5, q = 101 au, e = 0.95, ϖ = 0°). This symmetrization of the phase space of grains with 0.3 <β< 0.7 leads to a symmetrization of their spatial distribution. The arc seen on the periastron side in Fig. 9e is gone in Fig. 9f. Figure 8b shows how these grains then contribute to and strengthen the halo on the apastron side.

Comparison of the panels in Fig. 7 suggests that the halo on the periastron side is strengthened even more. While the circular belt produces a halo with a classical τ ∝ r^{− 3 / 2} behaviour (Strubbe & Chiang 2006; Krivov et al. 2006), the halos on the periastron sides of eccentric belts are closer to τ ∝ r^{− 1... − 3 / 4}. Optical depth falls off more steeply on the apastron sides. The slopes of all curves converge beyond ~400 au.

4.3. Inward transport through drag forces

The PR drag and stellar wind drag cause grains to spiral towards the star on timescales set by their size-dependent susceptibility to radiation and wind pressure, stellar luminosity and mass loss rate, and the orbital semi-major axes and eccentricities of the grains (see, e.g., Robertson 1937; Wyatt & Whipple 1950; Burns et al. 1979). The classical results for the orbit-averaged reduction rates of as and es are Other orbital elements are not affected secularly and non-relativistically. When grains are released from circular orbits in a source belt, drag alone produces constant optical depth τ(r) towards the star. Collisional sinks (Wyatt et al. 1999; Wyatt 2005) make τ decrease further in. Additional sources, such as active comets (Leinert et al. 1983), increase τ.

Figure 7 shows how drag and collisions shape the optical depth profiles in the inner regions of our model runs. From the peak in the parent belt to its inner edge, optical depth drops by about 2 orders of magnitude. Closer to the star, the slope flattens out as collisions become less important.

In the runs with eccentric parent belts, optical depths differ between apastron and periastron sides. Values on the apastron sides are systematically lower because drag rates are higher there. This can be explained in more detail with the following analytic model. At a given distance r from the star, the optical depth is determined by three factors: (1) the rate at which cross section gets dragged accross that distance; (2) the azimuthal spread of that cross section, caused by orbital speed v; and (3) the radial spread, caused by radial drift speed ṙ times orbital period P. All these quantities differ between periastron and apastron side. On the periastron side, we have (32)and for r = q = a_q(1 − e_q), the product of orbital period and orbital speed is (33)The radial component of the orbital velocity vanishes and ṙ is given by the PR induced reduction of periastron distance (34)Inserting Eqs. (30) and (31) into (34) results in (35)and (36)At apastron, where r = Q = a_Q(1 + e_Q), we find (37)For the flux of cross section, we assume (38)i.e. material flux does not change from periastron to apastron of a single orbit. The ratio of fluxes at the same distance on both sides is thus given by (39)For the radial dependence of under the action of drag and collisions, we adopt the analytical result of Wyatt et al. (1999), i.e. (40)where r₀ = a_b = 100 au and η ≈ 100 is the ratio of drag and collision timescales in the parent belt.

The combination of Eqs. (32), (36), (37), (39), and (40) leads to (41)The first term on the right-hand side accounts for collisional loss. It dominates close to the parent belt and for high η. The remainder describes pure drag, dominating close to the star and for low η.

The PR drag reduces orbital eccentricities as grains spiral in. At a given distance r from the star, grains on the apastron side therefore have an eccentricity e_Q(r) that is lower than the corresponding eccentricity e_q(r) of grains on the periastron side because r is reached later on the apastron side. The exact relation between the two eccentricities can be deduced from integrating the evolutions of a and e simultaneously. From (42)Wyatt & Whipple (1950) obtain (43)for PR drag between states 1 and 2. For q = q₁, Q = Q₂ and r = q = Q, we find (44)which can be solved numerically to find e_Q( = e₂) as a function of e_q( = e₁).

In Fig. 11, we show the profiles of optical depth expected from this analytic approach. The comparison with Fig. 7 shows that these profiles can well reproduce the results of our ACE runs outside of the parent region, i.e. for r<q_min = (a_b − Δa_b)(1 − e_b − Δe_b) on the periastron side and r<Q_min = (a_b − Δa_b)(1 + e_b − Δe_b) on the apastron side. The lower panel of Fig. 11 shows the resulting asymmetry ratios for a range of eccentricities e_b and distances r. For e_b = 0.2 and r = 0.5a_b, we find τ_Q/τ_q = 0.6. The asymmetry is vanishing slowly with decreasing distance from the star because typical eccentricities also decrease. This trend is seen both in the ACE output and the analytic model.

For e = e_b ≈ 0.1 observed for the outer Fomalhaut disc, we predict a flux deficit of 20% at r = 0.5a_b ≈ 70 au on the apastron side compared to the periastron side. This asymmetry is not seen in currently available observational data because these are limited by either low resolution in the case of Herschel/PACS (Acke et al. 2012) or sensitivity in the cases of the Hubble Space Telescope (Kalas et al. 2005) and the Atacama Large Millimetre Array (Boley et al. 2012). If the dust in the region interior to the outer Fomalhaut belt does not exhibit such an asymmetry, this would speak against inward drag as the dominating mechanism for replenishment.

The above analysis assumes that grains in the PR region start from orbits that follow those of the parent bodies, i.e. they are large enough not to be affected strongly by radiation pressure. However, Fig. 8 shows that smaller grains with potentially higher initial eccentricities contribute as well. A higher typical eccentricity of PR grains would further strengthen the asymmetry. Despite uncertainties in e and the simplified collisional depletion, the effect is robust and provides a testable prediction for the asymmetry of the drag-filled inner regions of eccentric belts.

4.4. Comparison with precession timescales

In the presented ACE runs we assumed a fixed average belt eccentricity and orientation, which are both present from the beginning and static throughout the simulation. At the same time we show in Fig. 4 that a modest Neptune at 30 au can already induce precession periods as short as ~10⁸ yr for a distant belt at 100 au. In the regime of observable dust, precession can thus act on timescales longer than those for collisions, but shorter than those of drag. If precession periods and forced eccentricities were equal for all objects, the effects on dust distribution would be negligible. The disc would precess as a whole while collisions and drag take place. However, differential precession due to different semi-major axes and β ratios twists the disc. Dragged-in dust precesses faster and dust in the halo precesses slower. In regions where collision timescales are longer, precession smears the distribution of complex eccentricities, increasing the spread in Δe, collision velocities, and potentially, depletion rates. An updated model that accounts for this process is in preparation.

5. Spectral energy distributions and images

Figure 12 illustrates the mild influence that the combination of geometrical offset, collisions, and drag has on the overall SED. Even for the rather eccentric cases e_b = 0.2 and s_b = 0.4, the differences in the wavelength range above a few tens of microns do not exceed 10%. It is only at shorter wavelengths that fluxes from eccentric discs become significantly higher, which is mainly due to pericentre glow, i.e. grains near the belt periastron having higher temperatures.

At wavelengths λ ≳ 100 ?m, the SEDs reflect the differences in size distributions described in Sect. 4.2. Shallower size distributions in more eccentric belts result in shallower SEDs. The discrepancy of 10% over a factor of ten in wavelength for e_b = 0.4 corresponds to a difference in power-law slopes of log (1.1) / log (10) ≈ 0.04. There exist several effects that have stronger impacts on SED slopes, but an unknown eccentricity adds to the uncertainty in the derivation of these other parameters. SED slopes are commonly used to infer underlying grain size distributions, which in turn are related to, e.g., collisional physics. Assuming that a difference of 0.04 in the SED slopes translates to no more than 0.04 in the inferred slopes of the size distributions (cf. Draine 2006), we conclude that other observational and modelling uncertainties dominate.

The panels in Fig. 13 show the fiducial discs with belt eccentricities e_b = 0.0, 0.2, and 0.4 in thermal emission at 24 ?m, 160 ?m, and 1.2 mm. Characteristic grain sizes s_c at different wavelengths λ can be estimated from s_c ≈ λ/ 2π (Backman & Paresce 1993). At λ = 1.2 mm, the halo is invisible and the discs appear as narrow belts because the dominant grains have s_c ≈ 200 ?m, and with β ≈ 0.01, are only weakly affected by radiation pressure. The drag timescales are such that these large grains just start to fill the inner gap after a few times 10⁸ yr (Fig. 4). At 160 ?m, correponding to s_c ≈ 25 ?m and β ≈ 0.08, the belt is wider and the halo and the inner region start to become visible. As long as drag is not important, though, the low β of the grains observed makes the chosen initial distribution, which corresponds to the set-up by Lee & Chiang (2016), which is a good proxy to the collisional steady state at these longer wavelengths.

At λ = 24 ?m, we find s_c ≈ 4 ?m (β ≈ 0.5). Observations at this wavelength are thus sensitive to the distribution of barely bound grains around the A3 V star assumed in our simulations. Grains of this size only stay bound when released from the apastron side of the parent belt, strengthening the halo on the periastron side. As a result and in contrast to the longer wavelengths, emission on the periastron extends further away from the belt. Images and radial profiles become more symmetric with increasing distance.

At the inner edges of all belts, the surface brightness drops by about 1.5–2 orders of magnitude, in agreement with the drops in optical depth and at 160 ?m. However, Fig. 14 shows that brightness follows r^{− 2.5... − 3} at 24 ?min the drag-filled region, increasing strongly towards the star because of the increasing temperature. Although the profile flattens off further towards the star, where 24 ?m is no longer on the Wien side, if the dust is not intercepted by inner planets the behaviour can produce a significant total excess (e.g. Liou & Zook 1999; Reidemeister et al. 2011). The wiggles seen at distances r ≳ 200 au in Fig. 14 correspond to the wiggles in the radial profiles of optical depth in Fig. 7, which are artefacts of the mass binning.

Pericentre glow (Wyatt et al. 1999) and broadening of the belt towards its apastron induce an asymmetry between the peak brightnesses in these two loci. The radial cuts in Fig. 14 show ratios between the peaks for e_b = 0.4, which amount to factors of 1.6 at 1.2 mm, 2.6 at 160 ?m, and 3.1 at 24 ?m. Although the Wien side of the SED is very sensitive to temperature, the asymmetry is only slightly more pronounced at 24 ?m than at 160 ?m for two reasons. First, typical grains near belt apastron are smaller than those near periastron, reducing the temperature difference between both sides; and second, the periastron side of the belt is wider at 24 ?m. In the initial disc, where only the second reason applies, the brightness ratio is 4.4 at 24 ?m.

The actually observed contrast strongly depends on how well these peaks are resolved. If the narrow periastron side is PSF broadened to the width of the apastron side, the observable difference is drastically reduced. Pan et al. (2016) have analysed this effect and have shown that pericentre glow can turn into apocentre glow at long wavelengths, where thermal emission depends less on temperature and resolution is typically lower. Adopting the idea behind their Fig. 2, we smoothened the images with Gaussian kernels and plotted azimuthal profiles for e_b = 0.4 in our Fig. 15 for two stages. Our initial set-up is represented by t = 0 yr. After 10⁷ yr, collisional equilibrium is reached but PR drag has not brightened the innermost region yet. The PSF broadened discs clearly show reduced pericentre glow at 24 and 160 ?m. At 1.2 mm, the apocentre is brighter than the pericentre by a factor of 1.3, which is roughly consistent with 1 + e_b that has been derived for small e_b by Pan et al. (2016). The bump seen around 180° at 24 ?m on the right panel reflects the tighter radial confinement of the small grains near the belt periastron seen in Fig. 13n.

Images in scattered light are similar to 24 ?m emission because both trace the small, barely bound grains. The major difference lies in the radial slopes, which are steeper in thermal emission because there the exponential dependence on grain temperature factors in.

The scattering cross section and its angular dependence are very sensitive to model assumptions on grain morphology. While Mie theory for perfect spheres works well for s ≪ λ, it fails to reproduce the scattering phase function for irregularly shaped grains with s ≳ λ, where the morphology of surface and sub-surface layers becomes important. Compared to a polished sphere, a rougher surface can increase backscattering and reduce absorption (Pollack & Cuzzi 1980). At scattering angles far away from the strong forward diffraction peaks of large grains, the resulting phase functions are flatter. In models of debris discs, these more symmetric phase functions are often associated with the presence of smaller, spherical Mie grains. This degeneracy between small grains and grains with small structures is discussed by Min et al. (2010) and Hedman & Stark (2015) in their models for the Fomalhaut disc and the G and D68 rings of Saturn, respectively.

In Fig. 16, we illustrate this problem with a set of images for different scattering models. Results from Mie calculations are compared to empirical Henyey–Greenstein (HG) phase functions, ranging from mild forward-scattering with anisotropy parameter g = 0.3, to strong forward scattering with g = 0.94, and the three-component best-fit model that Hedman & Stark (2015) derive for the G ring of Saturn. In the latter, the strongest component has g = 0.995. As expected from the rather large grains in our collisional models, the Mie results are best matched by g close to unity. However, the model based on the Hedman & Stark fit shows the degree to which the total brightness of the disc in scattered light may be underestimated. Accordingly, non-Mie fits to observed discs find weaker anisotropy, with g< 0.5 (e.g. Kalas et al. 2005; Schneider et al. 2006, 2014; Debes et al. 2008; Thalmann et al. 2011).

6. Conclusions

With a new version of our collisional code ACE, we have studied debris discs that are sustained by eccentric belts of parent bodies. We identified a set of features and asymmetries that are caused by the combined effects of global eccentricity, radiation pressure, collisional evolution, and drag forces. The features most easily observed are as follows:

• 1.

  On dynamical timescales, the different radiation pressureblowout limits on opposing disc sides create an asymmetric halo.At shorter wavelengths, where small grains dominate, the haloappears more extended beyond the periastron side of the parentbelt. For the larger grains seen at longer wavelengths, the apastronside is more extended.

• 2.

  In collisional equilibrium, the abundance ratios between barely bound grains and grains that are around an order of magnitude larger are different for opposing sides of the belt. On the periastron side, average grains are larger, while on the apastron side, grains are smaller. This size difference reduces the temperature difference between the two sides and weakens the brightness asymmetry expected from pericentre glow.

• 3.

  Poynting–Robertson and stellar wind drag induce an additional asymmetry because they reduce apocentre distances at a higher rate than pericentre distances. Apocentre sides of the drag-filled regions inside of eccentric belts are therefore populated more tenuously. The relative difference between the two sides is comparable with the belt eccentricity, i.e. 10 % for e_b ≈ 0.1. Towards the star, the asymmetries reduce along with the eccentricities.

• 4.

  Belt eccentricity affects SEDs only weakly. The azimuthal variation in size distribution and the pericentre glow result in an SED that is broader overall. The effect is most notable in the mid-IR, but significant only for high (≳0.4) eccentricities.

• 5.

  Interpretation of near-infrared images crucially depends on the scattering model on which it is based. Empirical models valid for the larger grains in the parent belt will be inadequate for the smaller grains that form the outer halo. Mie theory has the benefit of retaining the dependence on grain size, but on the other hand, it does not approximate the scattering phase functions of larger grains well.

A more detailed analysis of the influence of dust optical properties and the disc viewing geometry on the observables is on its way. Another update to ACE will allow us treat the secular precession of orbits by a perturbing planet in parallel with the collisional evolution.

Acknowledgments

We thank an anonymous referee for constructive suggestions that helped improve the presentation. Part of this work was supported by the Deutsche Forschungsgemeinschaft through grants LO 1715/2-1, KR 2164/13-1, KR 2164/15-1, WO 857/13-1, and WO 857/15-1.

References

All Figures