Collisions and drag in debris discs with eccentric parent belts

TL;DR

This study investigates how collisions and drag forces influence asymmetries in debris discs with eccentric belts, revealing complex interactions that affect observable features and spectral energy distributions over long timescales.

Contribution

It extends a numerical model to include azimuthal variations, analyzing the combined effects of eccentricity, collisions, and drag on debris disc morphology and observables.

Findings

Radiation pressure causes size-dependent asymmetries in the halo.

Collisional equilibrium favors smaller grains on the apastron side.

Poynting-Robertson drag creates a more tenuous interior region on the apastron side.

Abstract

Context: High-resolution images of circumstellar debris discs reveal off-centred rings that indicate past or ongoing perturbation, possibly caused by secular gravitational interaction with unseen stellar or substellar companions. The purely dynamical aspects of this departure from radial symmetry are well understood. However, the observed dust is subject to additional forces and effects, most notably collisions and drag. Aims: To complement the studies of dynamics, we therefore aim to understand how new asymmetries are created by the addition of collisional evolution and drag forces, and existing ones strengthened or overridden. Methods: We augmented our existing numerical code "Analysis of Collisional Evolution" (ACE) by an azimuthal dimension, the longitude of periapse. A set of fiducial discs with global eccentricities ranging from 0 to 0.4 is evolved over giga-year timescales. Size distribution and spatial variation of dust are analysed and interpreted. The basic impact of belt eccentricity on spectral energy distributions (SEDs) and images is discussed. Results: We find features imposed on characteristic timescales. First, radiation pressure defines size cutoffs that differ between periapse and apoapse, resulting in an asymmetric halo. The differences in size distribution make the observable asymmetry of the halo depend on wavelength. Second, collisional equilibrium prefers smaller grains on the apastron side of the parent belt, reducing the effect of pericentre glow and the overall asymmetry. Third, Poynting-Robertson drag fills the region interior to an eccentric belt such that the apastron side is more tenuous. Interpretation and prediction of the appearance in scattered light is problematic when spatial and size distribution are coupled.

Equations

$$e_{\text{f}}=e_{\text{planet}}\frac{b_{3/2}^{(2)}(\alpha )}{b_{3/2}^{(1)}(\alpha )}=\left[\frac{5}{4}\alpha +o(\alpha )\right]e_{\text{planet}},$$

$$A=\sqrt{\frac{G{M}_{\ast }(1-\beta )}{a_{\text{b}}^3}}\frac{M_{\text{planet}}}{4M_{\ast }(1-\beta )}\alpha b_{3/2}^{(1)}(\alpha ),$$

$$T_{\text{prec}}^{\text{full}}\equiv \frac{2\pi }{A}=\frac{M_{\ast }(1-\beta )}{M_{\text{planet}}}\frac{4P_{\text{b}}}{\alpha b_{3/2}^{(1)}(\alpha )},$$

$$T_{\text{prec}}\approx T_{\text{prec}}^{\text{full}}\frac{e_{\text{b}}}{2\pi e_{\text{f}}}$$

$$\Delta T_{\text{prec}}\equiv \frac{\Delta e_{\text{b}}}{\Delta A},$$

$$\Delta A=\frac{\mathrm{d}(A{e}_{\text{f}})}{\mathrm{d}\alpha }\frac{\mathrm{d}\alpha }{\mathrm{d}{a}_{\text{b}}}\Delta a_{\text{b}}=A{e}_{\text{f}}\frac{\Delta a_{\text{b}}}{a_{\text{b}}}\left[\frac{5}{2}+\frac{\alpha }{b_{3/2}^{(2)}}\frac{\mathrm{d}{b}_{3/2}^{(2)}}{\mathrm{d}\alpha }\right].$$

$$\Delta T_{\text{prec}}=\frac{\Delta e_{\text{b}}}{A{e}_{\text{f}}}\frac{a_{\text{b}}}{\Delta a_{\text{b}}}{\left[\frac{9}{2}+\frac{7}{2}{\alpha }^2+\frac{119}{32}{\alpha }^4+o({\alpha }^4)\right]}^{-1},$$

$$\Delta T_{\text{prec}}\approx T_{\text{prec}}^{\text{full}}\frac{\Delta e_{\text{b}}}{2\pi e_{\text{f}}}\frac{a_{\text{b}}}{\Delta a_{\text{b}}}\frac{1-{\alpha }^2}{4}$$

$$\frac{\Delta e_{\text{b}}=e_{\text{b}}}{2\pi e_{\text{f}}}\frac{a_{\text{b}}}{\Delta a_{\text{b}}}\frac{1-{\alpha }^2}{4}>\frac{e_{\text{b}}}{2\pi e_{\text{f}}}+N$$

$$\frac{\Delta a_{\text{b}}}{a_{\text{b}}}<\frac{1-{\alpha }^2}{4}{\left(1+\frac{2\pi N{e}_{\text{f}}}{e_{\text{b}}}\right)}^{-1}.$$

$${\dot{n}}_i=\sum _{jk}{G}_{ijk}{n}_j{n}_k-\sum _j{L}_{ij}{n}_i{n}_j+\sum _j{T}_{ij}{n}_j,$$

$$Q_{ij}^{\ast }$$

$$m_{\text{frag}}=\frac{Q}{Q_{\ast }}\left(\frac{m_i}{2}+m_j\right).$$

$$m_{\text{lf}}\equiv \frac{2.0+\eta }{-\eta -1.0}{m}_{\text{frag}}=\frac{1}{5}{m}_{\text{frag}}.$$

$$\frac{r}{a}$$

$$p$$

$$e$$

$$a=a_{\text{t}}\frac{1-\beta }{1-2\beta a_{\text{t}}/r}.$$

$${\beta }_{\text{lim}}=\frac{r}{2a_{\text{t}}}.$$

$$e\sin \theta =\frac{m_{\text{t}}^{\prime }}{m_{\text{sum}}^{\prime }}\sqrt{\frac{p}{p_{\text{t}}}}{e}_{\text{t}}\sin {\theta }_{\text{t}}+\frac{m_{\text{p}}^{\prime }}{m_{\text{sum}}^{\prime }}\sqrt{\frac{p}{p_{\text{p}}}}{e}_{\text{p}}\sin {\theta }_{\text{p}}.$$

$$e\cos \theta =\frac{p}{r}-1.$$

$${\theta }^{\prime }-\theta =\varpi -{\varpi }^{\prime }.$$

$${\theta }_{\text{t}}=\arcsin \frac{A}{\sqrt{C}}\pm \arccos \frac{B}{\sqrt{C}}$$

$$\sin {\theta }_{\text{t}}=\frac{AB\pm D\sqrt{C-B^2}}{C},\cos {\theta }_{\text{t}}=\frac{BD\mp A\sqrt{C-B^2}}{C},$$

$$A$$

$$B$$

$$C$$

$$D$$

$$\frac{\Delta Q_{\text{b}}}{Q_{\text{b}}}=\sqrt{{\left(\frac{\Delta a_{\text{b}}}{a_{\text{b}}}\right)}^2+{\left(\frac{\Delta e_{\text{b}}}{1+e_{\text{b}}}\right)}^2}$$

$$\frac{\mathrm{d}i}{\mathrm{d}\ln e}={\left\{{\left(\ln 1.25\right)}^{-3}+{\left[{\left(0.1/e\right)}^3+{\left(0.1/1.5\right)}^3\right]}^{-1}\right\}}^{1/3}.$$

$$\dot{a}$$

$$\dot{e}$$

$${\tau }_q=\frac{{\dot{\sigma }}_q(r)}{|{\dot{r}}_q|P_q{v}_q},$$

Full text

*Aims. *To complement the studies of dynamics, we therefore aim to understand how new asymmetries are created by the addition of collisional evolution and drag forces, and existing ones strengthened or overridden.

*Methods. *We augmented our existing numerical code “Analysis of Collisional Evolution” (ACE) by an azimuthal dimension, the longitude of periapse. A set of fiducial discs with global eccentricities ranging from 0 to 0.4 is evolved over giga-year timescales. Size distribution and spatial variation of dust are analysed and interpreted. The basic impact of belt eccentricity on spectral energy distributions (SEDs) and images is discussed.

*Results. *We find features imposed on characteristic timescales. First, radiation pressure defines size cutoffs that differ between periapse and apoapse, resulting in an asymmetric halo. The differences in size distribution make the observable asymmetry of the halo depend on wavelength. Second, collisional equilibrium prefers smaller grains on the apastron side of the parent belt, reducing the effect of pericentre glow and the overall asymmetry. Third, Poynting–Robertson drag fills the region interior to an eccentric belt such that the apastron side is more tenuous. Interpretation and prediction of the appearance in scattered light is problematic when spatial and size distribution are coupled.

1 Introduction

Observed dust in debris discs is produced in collisions amongst orbiting planetesimals. Resolved images at submillimeter wavelengths, which are tracing large grains, suggest the dust parent bodies to be 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 a global asymmetry between the two wings. It is possible that these asymmetry also derives 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 the belts (Lee & Chiang, 2016; Esposito et al., 2016) or substellar companions exterior to them (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 no direct tracers of the underlying distribution of parent bodies from that 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; Vitense et al., 2012; Kral et al., 2013; Nesvold et al., 2013; Thébault et al., 2014; 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 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 Section 2. In Section 3 we introduce a new version of our collisional code ACE that 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 summarise our findings in Section 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 disc’s cavity are conceivable. For instance, Shannon et al. (2014) propose 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 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 $\varpi$ (Fig. 1), secular perturbation makes the eccentricity vector $(e\cos\varpi,e\sin\varpi)$ precess uniformly along a circle with a radius called proper eccentricity, $e_{\text{p}}$, centred around a forced eccentricity $e_{\text{f}}$ (Hirayama, 1918). The closer the combination of $e_{\text{p}}$ and $e_{\text{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 planet’s eccentricity vector, its absolute value being (Murray & Dermott, 2000):

[TABLE]

where $e_{\text{planet}}$ is the absolute eccentricity of the planet’s orbit. The $b$’s are Laplace coefficients, which only depend on the ratio of semi-major axes of the (interior) planet and a perturbed belt object, $\alpha\equiv{a_{\text{planet}}}/a_{\text{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_{\text{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 perturber’s mass does not influence $e_{\text{f}}$ (as long as $M_{\text{planet}}\ll M_{*}$), it determines the timescale on which this precession occurs. Based on the (approximate) angular precession frequency for a single perturber (Murray & Dermott, 2000),

[TABLE]

the time for a full precession cycle can be estimated from

[TABLE]

where $P_{\text{b}}=2\pi\sqrt{a_{\text{b}}^{3}/[GM_{*}(1-\beta)]}$ is the orbital period in the belt. The additional factor $1-\beta$ with $\beta\equiv F_{\text{pr}}/F_{\text{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), 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_{\text{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 will remain 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$-$\varpi$ 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_{\text{b}}$. The minimum time for the belt eccentricity to precess to $e_{\text{b}}$ is a fraction

[TABLE]

of a full cycle (see Fig. 1 right), potentially increased by $NT_{\text{prec}}^{\text{full}}$ if $N$ precession cycles have already passed. However, precession will also smear 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$–$\varpi$ would be twice as large as $e_{\text{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, in order to prevent this smearing out. We can define such a differential timescale as

[TABLE]

where $\Delta e_{\text{b}}$ is the spread in eccentricities and $\Delta A$ denotes the spread in linear change rates in eccentricity space:

[TABLE]

Inserting (6) into (5) and expanding the term in brackets in a series in $\alpha$, we obtain

[TABLE]

or with equation (3),

[TABLE]

because the coefficients in the series approach a value of $4$. The spread in eccentricities ($\Delta e_{\text{b}}$) remains lower than the average eccentricity ($e_{\text{b}}$) and the belt remains narrow as long as $\Delta T_{\text{prec}}>T_{\text{prec}}$. We find

[TABLE]

and

[TABLE]

This constraint on the belt width is mild for $N=0$ and small values of $\alpha$, but strong for $\alpha$ 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, $\alpha\ll 1$, and that the observed belt eccentricity is approximately the same as the forced eccentricity. We find $\Delta a_{\text{b}}/a_{\text{b}}<1/4$ for $N=0$ and $\Delta a_{\text{b}}/a_{\text{b}}\lesssim 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 Section 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_{\text{b}}$, as a key paremeter. However, its interpretation in these two cases is different. In the low-$e_{\text{p}}$ scenario, $e_{\text{b}}$ is equal to the forced eccentricity $e_{\text{f}}$ (Fig. 1 left). In that case, the planet orbit’s apsidal line is aligned with the major axis of the belt. In the slow-precession scenario, $e_{\text{b}}$ is not equal to $e_{\text{f}}$. Instead, it represents the instantaneous value of the complex eccentricity $e$ (Fig. 1 right). In that case, the planetary orbit will be 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 or bins or 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, with LIDT-DD) or close stellar flybys (Nesvold et al., 2017, with SMACK). Due 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. LIPAD (Levison et al., 2012) is 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, parameterised 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; Reidemeister et al., 2011; Löhne et al., 2012; Krivov et al., 2013), 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 $\varpi\equiv\Omega+\omega$, the sum of the longitude of the ascending node ($\Omega$) and the argument of periapsis ($\omega$). 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 $\varpi$. Collisions among pair-wise combinations of bins are possible at up to two distinct points defined by the colliders’ $q$s, $e$s, and $\varpi$s. Collision velocities, rates, and outcomes follow directly (Krivov et al., 2006).

The discretized master equation to be integrated over time reads

[TABLE]

where $n_{i}$ is the number (or mass) of objects in the bin specified by multi-index $i\equiv(i_{m},i_{q},i_{e},i_{\varpi})$. Coefficients $G_{ijk}$ denote the gain, 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, the specific rate at which objects of type $i$ are removed in collisions with $j$. ACE 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, Poynting–Robertson drag reduces $q$ in bin $j$ at a rate $\dot{q}_{j}$. Given a bin width $\Delta q_{j}$, the contents of bin $j$ are moved towards $i$ at a rate $T_{ij}n_{j}=n_{j}\dot{q}_{j}/\Delta q_{j}$.

With collision physics depending on the orbits’ orientations only through the difference $\varpi_{i}-\varpi_{j}$, the dependence of $G_{ijk}$ and $L_{ij}$ on $i_{\varpi}$, $j_{\varpi}$, and $k_{\varpi}$ is only through pair-wise differences $i_{\varpi}-j_{\varpi}$, etc. The relation $G_{ijk}=G_{i^{\prime}j^{\prime}k^{\prime}}$ for $i^{\prime}=(i_{m},i_{q},i_{e},i_{\varpi}+o)$, $j^{\prime}=(j_{m},j_{q},j_{e},j_{\varpi}+o)$, and $k^{\prime}=(k_{m},k_{q},k_{e},k_{\varpi}+o)$ with $o\in N$ can be employed to speed up the calculations; instead of looping colliders over $j_{\varpi}$ and $k_{\varpi}$, $k_{\varpi}-j_{\varpi}$ and $j_{\varpi}$ are used, making the innermost loop over $j_{\varpi}$ trivial.

3.2 Collision outcomes

Depending on the colliders’ masses 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):

[TABLE]

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_{\text{imp}}$. By $s_{ij}\equiv(s_{i}^{3}+s_{j}^{3})^{1/3}$ we denote the equivalent radius of a sphere with the combined volume of the colliders. The last line in equation (12) approximates the specific energy required to overcome self-gravity. It is important only for radii $s\gtrsim 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 colliders’ remnants. In the model, the total mass in escaping fragments is proportional to impact energy:

[TABLE]

The fragment mass distribution is assumed to follow a power law with exponent $\eta=-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

[TABLE]

where $m_{\text{frag}}$ is the total mass in escaping fragments. The same prescription was used in previous ACE versions.

The remnants from erosive collisions are treated somewhat different 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. Only if the remnants move to a different bin, is their mass added to the loss of one bin and to the gain of another. However, in what follows we will 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 $\beta$ 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$:

[TABLE]

where $m_{\text{t}}^{\prime}\equiv m_{\text{t}}\sqrt{1-\beta_{\text{t}}}$ and $m_{\text{p}}^{\prime}\equiv m_{\text{p}}\sqrt{1-\beta_{\text{p}}}$ denote the effective masses of target and projectile, and $m_{\text{sum}}^{\prime}\equiv(m_{\text{t}}+m_{\text{p}})\sqrt{1-\beta}$ the sum of the collider masses, weighed by the $\beta$ ratio of the fragment. The signs of $e$ and $p$ are equal to that of $1-\beta$, corresponding to anomalous hyperbolae with $e<0$ for $\beta>1$.

For grains that are released from a macroscopic body with $\beta_{t}=0$ on impact of a projectile with $m_{\text{p}}\ll m_{\text{t}}$, equation (15) is reduced to

[TABLE]

The blowout limit is reached where $a\rightarrow\infty$, corresponding to (Kresák, 1976)

[TABLE]

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 $\beta\geq 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 $\theta$ of a freshly released fragment can be determined from

[TABLE]

and

[TABLE]

At the mutual crossing points of two orbits, the difference of true anomalies equals the (negative) difference between their longitudes of periapse:

[TABLE]

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))

[TABLE]

or

[TABLE]

where

[TABLE]

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 will 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: $M=1.92\leavevmode\nobreak\ M_{\text{Sun}}$, $L=16.6\leavevmode\nobreak\ L_{\text{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 $\Delta Q_{\text{b}}/Q_{\text{b}}$ at the apastron $Q_{\text{b}}=a_{\text{b}}(1+e_{\text{b}})$ of a belt can be estimated from

[TABLE]

if $a$ and $e$ vary independently. A given radial HWHM of, for instance, $\delta r/r=10$ % does not only limit $\Delta a/a\leq 10$ %, but also $\Delta e_{\text{b}}=e_{\text{p}}\leq 0.1$. For the parent bodies in our main simulation runs, we therefore assume initial distributions that are confined to circular regions in the $(e\cos\varpi,e\sin\varpi)$ plane, centred on $(e_{\text{b}},0)$, with radii $e_{\text{p}}=0.1$. The values used for $e_{\text{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\lesssim 0.4$ and for $e>1$, with factors $\Delta e/e\approx 0.25$. In order to preserve accuracy for barely bound grains, the step size is additionally limited to $\Delta e\leq 0.1$ for $e<1$, resulting in nearly linear steps in the range $0.4\lesssim 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

[TABLE]

The resulting number of eccentricity bins is 26. The grid of orbit orientations has 32 linear steps in $\varpi$, covering $[0,2\pi)$. The resolution elements $\Delta e\times e\Delta\varpi$ thus measure $0.25e\times 0.2e$ for $e\lesssim 0.4$ and $0.1\times 0.2e$ for $0.4\lesssim e<1$. Figure 3 illustrates the sub-grid of $e$ and $\varpi$ in a polar plot.

In the radial direction, a grid of logarithmically spaced periastron distances $q\equiv 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_{\text{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\lesssim 30$ µm (corresponding to $\beta\gtrsim 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^{6}$ in the runs presented here.

The initial size distribution is assumed to follow a power law $n(s)\propto s^{\nu}T_{\text{orbit}}/T_{0}$ with $\nu=-3.66$, normalized to a total mass of 2 earth masses. We follow Strubbe & Chiang (2006) and Lee & Chiang (2016) in scaling the initial abundances of grains with their orbital timescales $T_{\text{orbit}}$ with respect to that of large grains ($\beta=0$), $T_{0}$. 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 setup to ease comparison with other work. That choice does not influence the dust distribution towards which the subsequent collisional evolution will quickly converge.

4 Resulting size and radial distributions

In this section, we will 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^{7}$ yr), enough time has passed to bring grains with sizes $s\lesssim 10$ cm to collisional equilibrium. Later, at an “evolved” stage ($t=2\times 10^{8}$ yr), grains with sizes $s\lesssim 100$ µm are in P–R 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 time scales can be found in Section 4.4.

4.1 Dynamical consequences

For the large parent bodies that are unaffected by radiation pressure, a non-vanishing mean eccentricity vector $(e\cos\varpi,e\sin\varpi)$ 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 will 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_{\text{b}}$, shown in Figure 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, changing the size distribution in the parent belt. Blowout occurs for smaller values of $\beta$ for grains that are launched near periastron of a parent orbit. The excess velocity there helps them 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_{\text{b}}$. According to equation (18), the ratio between the blowout sizes at periastron and apastron is given by $(1+e_{\text{b}})/(1-e_{\text{b}})$. This ratio reaches a value of $7/3$ for $e_{\text{b}}=0.4$, which is consistent with the ACE results.

In the left panel of Figure 7, radial profiles of the normal optical depth, $\tau$, are plotted for $t\approx 0$ yr. At large radii, these profiles drop almost as $\tau\propto r^{-3/2}$, the behaviour expected for discs in equilibrium (Krivov et al., 2006). This is by design because we employed the initial setup 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 their orbits’ apastra, 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 Figure 3 show how larger grains have orbits that are still aligned azimuthally with the parent belt. Their periastra (as described by $e$ and $\varpi$ in that figure) remain close to the belt’s periastron. The small grains shown in the bottom row of Figure 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 will be formed by grains larger than those on the periastron side. Figure 8 shows this effect for a disc with $e_{\text{b}}=0.4$. There, the $\beta$ 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 $\beta<0.3$.

The widths of the size distributions on the two sides differ as illustrated in Figure 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, and so 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 Figure 9e as a dark arc on the right-hand side of that panel. Lee & Chiang have an equivalent arc in their Figure 1 (and Figure 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, with the initial round of fragments all being produced from a single parent.

4.2 Effects of collisions

In collisional equilibrium a size distribution of infinite extent can be well-described by a power law $n(s)\propto s^{\alpha}$, with $\alpha\sim-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 cutoffs. 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 cutoff (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 equation (16) equals $e_{\text{b}}$. For grains produced in a parent belt with average proper eccentricities $e_{\text{p}}=0.05$, this break is expected near $\beta=e_{\text{p}}/(1+e_{\text{p}})\approx 0.05$, corresponding to grain sizes $s\approx 45$ µm in our setup. The right panel of Figure 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 do 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_{\text{b}}$, reaching an order of magnitude for $e_{\text{b}}=0.4$ (Fig. 6). This can be understood because the two sides are coupled; the larger grains coming from the belts’ periastra 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. In 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_{\text{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 wavelength. The detailed effects on the observable SEDs and images will be discussed in Section 5.

In Figure 3 the most notable difference between the initial stage and the intermediate stage is the filling of regions that could not be reached initially. In particular, grains with $s=4.3$ µm ($\beta=0.5$) appear on bound orbits ($e<1$) with periastra aligned with the periastron of the belt. These grains cannot be just produced in the parent belt. P–R drag cannot be responsible either because it cannot alter $\varpi$ and it cannot re-bind unbound grains in significant numbers (which is why this process is not modelled in ACE). Thus, they must stem from medium sized grains. An illustration of this process is given in Figure 10, where the orbits of three particle types actually represented in the ACE grid are shown. The target ($s=12$ µm, $\beta=0.17$, $q=98$ au, $e=0.22$, $\varpi=23^{\circ}$) and the projectile ($s=5.4$ µm, $\beta=0.4$, $q=85$ au, $e=0.75$, $\varpi=-56^{\circ}$) 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, $\beta=0.5$, $q=101$ au, $e=0.95$, $\varpi=0^{\circ}$). This symmetrization of the phase space of grains with $0.3<\beta<0.7$ leads to a symmetrization of their spatial distribution. The arc seen on the periastron side in Figure 9e is gone in Figure 9f. Figure 8b shows how these grains then contribute to and strengthen the halo on the apastron side.

Comparison of the panels in Figure 7 suggests that the halo on the periastron side is strengthened even more. While the circular belt produces a halo with a classical $\tau\propto r^{-3/2}$ behaviour (Strubbe & Chiang, 2006; Krivov et al., 2006), the halos on the periastron sides of eccentric belts are closer to $\tau\propto r^{-1\ldots-3/4}$. The apastron sides see a steeper falloff. The slopes of all curves converge beyond $\sim 400$ au.

4.3 Inward transport through drag forces

Poynting–Robertson 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 grains’ orbital semi-major axes and eccentricities (see, e. g, Robertson, 1937; Wyatt & Whipple, 1950; Burns et al., 1979). The classical results for the orbit-averaged reduction rates of $a$s and $e$s are

[TABLE]

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 $\tau(r)$ towards the star. Collisional sinks (Wyatt et al., 1999; Wyatt, 2005) make $\tau$ decrease further in. Additional sources, such as active comets (Leinert et al., 1983), increase $\tau$.

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 $\dot{\sigma}$ 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 $\dot{r}$ times orbital period $P$. All these quantities differ between periastron and apastron side. On the periastron side, we have

[TABLE]

and for $r=q=a_{q}(1-e_{q})$, the product of orbital period and orbital speed is

[TABLE]

The radial component of the orbital velocity vanishes and $\dot{r}$ is given by the P–R induced reduction of periastron distance:

[TABLE]

Inserting equations (30) and (31) into (34) results in

[TABLE]

and

[TABLE]

At apastron, where $r=Q=a_{Q}(1+e_{Q})$, we find

[TABLE]

For the flux of cross section, we assume

[TABLE]

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

[TABLE]

For the radial dependence of $\dot{\sigma}_{a}(r)$ under the action of drag and collisions, we adopt the analytical result of Wyatt et al. (1999):

[TABLE]

where $r_{0}=a_{\text{b}}=100$ au and $\eta\approx 100$ is the ratio of drag and collision time scales in the parent belt.

The combination of equations (32), (36), (37), (39), and (40) leads to

[TABLE]

The first term on the right-hand side accounts for collisional loss. It dominates close to the parent belt and for high $\eta$. The remainder describes pure drag, dominating close to the star and for low $\eta$.

P–R drag reduces orbital eccentricities as grains spiral in. At a given distance $r$ from the star, grains on the apastron side will 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

[TABLE]

Wyatt & Whipple (1950) obtain

[TABLE]

for P–R drag between States 1 and 2. For $q=q_{1}$, $Q=Q_{2}$, and $r=q=Q$ this leads to

[TABLE]

which can be solved numerically to find $e_{Q}(=e_{2})$ as a function of $e_{q}(=e_{1})$.

In Figure 11, we show the profiles of optical depth expected from this analytic approach. Comparison with Figure 7 shows that these profiles can well reproduce the results of our ACE runs outside of the parent region, i. e. for $r<q_{\text{min}}=(a_{\text{b}}-\Delta a_{\text{b}})(1-e_{\text{b}}-\Delta e_{\text{b}})$ on the periastron side and $r<Q_{\text{min}}=(a_{\text{b}}-\Delta a_{\text{b}})(1+e_{\text{b}}-\Delta e_{\text{b}})$ on the apastron side. The lower panel of Figure 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_{\text{b}}$, we find $\tau_{\text{Q}}/\tau_{\text{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_{\text{b}}\approx 0.1$ observed for the outer Fomalhaut disc, we predict a flux deficit of 20 % at $r=0.5a_{\text{b}}\approx 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 HST (Kalas et al., 2005) and ALMA (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.

Note that the above analysis assumes that grains in the P–R 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, Figure 8 shows that smaller grains with potentially higher initial eccentricities contribute as well. A higher typical eccentricity of P–R 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, both present from the beginning and static throughout the simulation. At the same time we show in Figure 4 that a modest Neptune at 30 au can already induce precession periods as short as $\sim 10^{8}$ 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 $\beta$ ratios will twist the disc. Dragged-in dust will precess faster, dust in the halo slower. In regions where collision timescales are longer, precession will smear the distribution of complex eccentricities, increasing the spread in $\Delta 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_{\text{b}}=0.2$ and $s_{\text{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 $\lambda\gtrsim 100$ µm, the SEDs reflect the differences in size distributions described in Section 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_{\text{b}}=0.4$ corresponds to a difference in power-law slopes of $\log(1.1)/\log(10)\approx 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 Figure 13 show the fiducial discs with belt eccentricities $e_{\text{b}}=0.0$, $0.2$, and $0.4$ in thermal emission at 24 µm, 160 µm, and 1.2 mm. Characteristic grain sizes $s_{\text{c}}$ at different wavelengths $\lambda$ can be estimated from $s_{\text{c}}\approx\lambda/2\pi$ (Backman & Paresce, 1993). At $\lambda=1.2$ mm, the halo is invisible and the discs appear as narrow belts because the dominant grains have $s_{\text{c}}\approx 200$ µm, and with $\beta\approx 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^{8}$ yr (Fig. 4). At 160 µm, correponding to $s_{\text{c}}\approx 25$ µm and $\beta\approx 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 $\beta$ of the grains observed makes the chosen initial distribution, which corresponds to the setup by Lee & Chiang (2016), a good proxy to the collisional steady state at these longer wavelengths.

At $\lambda=24$ µm, we find $s_{\text{c}}\approx 4$ µm ($\beta\approx 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, surface brightness drops by about 1.5–2 orders of magnitude, in agreement with the drops in optical depth and at 160 µm. However, Figure 14 shows that brightness follows $r^{-2.5\ldots-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, the behaviour can produce a significant total excess – if the dust is not intercepted by inner planets (e. g., Liou & Zook, 1999; Reidemeister et al., 2011). The wiggles seen at distances $r\gtrsim 200$ au in Figure 14 correspond to the wiggles in the radial profiles of optical depth in Figure 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 Figure 14 show ratios between the peaks for $e_{\text{b}}=0.4$ that amount to factors of 1.6 at 1.2 mm, 2.6 at 160 µm, and 3.1 at 24 µm. Despite the Wien side of the SED being very sensitive to temperature, the asymmetry is only slightly more pronounced at 24 µm than at 160 µm for two reason: (1) because typical grains near belt apastron are smaller than those near periastron, reducing the temperature difference between both sides; and (2) because 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 will strongly depend 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) analyse this effect and show 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 Figure 2, we smoothened the images with Gaussian kernels and plotted azimuthal profiles for $e_{\text{b}}=0.4$ in our Figure 15 for two stages. Our initial setup is represented by $t=0$ yr. After $10^{7}$ yr, collisional equilibrium is reached but P–R 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, roughly consistent with $1+e_{\text{b}}$ derived for small $e_{\text{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 Figure 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\ll\lambda$, it fails to reproduce the scattering phase function for irregularly shaped grains with $s\gtrsim\lambda$, where the morphology of surface and sub-surface layers becomes important. Compared to a polished sphere, a rougher surface can increase back-scattering 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 Saturn’s G and D68 rings, respectively.

In Figure 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 Saturn’s G ring. 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; Debes et al., 2008; Thalmann et al., 2011; Schneider et al., 2014).

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 have 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:

On dynamical timescales, the different radiation pressure blowout limits on opposing disc sides create an asymmetric halo. At shorter wavelengths, where small grains dominate, the halo appears more extended beyond the periastron side of the parent belt. For the larger grains seen at longer wavelengths, the apastron side is more extended. 2. 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. 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_{\text{b}}\approx 0.1$. Towards the star, the asymmetries reduce along with the eccentricities. 4. 4.

Belt eccentricity affects spectral energy distributions 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 ($\gtrsim 0.4$) eccentricities. 5. 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.

Acknowledgements.

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.