© ESO, 2012

1. Introduction

Debris discs have been detected surrounding a significant fraction (between  ~5 and  ~30%) of main-sequence stars of all spectral types (e.g. Su et al. 2006; Trilling et al. 2008). They consist of collisionally interacting solid bodies spanning a wide size range, from  ~1–100 kilometre-sized parent bodies, probably leftovers from the planet formation process, down to micron-sized debris. As they constitute most of the geometrical cross-section, only the smallest,  ≤1 cm, particles of this collisional cascade are detectable, in thermal emission or scattered light (see review by Wyatt 2008).

For most discs that have been resolved, pronounced structures have been observed, such as two-sided asymmetries, spirals, warps, clumps or rings (e.g Kalas et al. 2005; Golimovski et al. 2006; Schneider et al. 2009). These features are indicative of complex dynamical environments, and numerous studies have have attempted to explain their origin (e.g. Krivov 2010). Most of these studies are based on N-body numerical codes that follow the dynamical evolution of the system, exploring, depending on each system’s characteristics, different scenarios: the perturbing influence of (often undetected) planets or stellar companions, violent transient events, or interaction with gas remnants. In such a collisionless N-body approach, the disc is sampled by a population of N_num test particles, whose trajectories can be precisely tracked. The forces that control their evolution are: the central star’s gravity, the gravitational pull of planets and other stars, radiation pressure, Poynting-Robertson (PR) drag, and stellar wind drag. The precision of the N-body method allows one to study fine effects such as orbital resonances or transient spatial inhomogeneities. However, this comes at the price that collisions are neglected, leading to several shortcomings, which can be more or less worrying depending on the issues that are being addressed. These possible shortcomings are:

• 1) No size-distribution evolution. Since collisions control the particle size distribution and its evolution (by merger,erosion, or shattering), collisionless N-body codes cannot predict how the distribution varies with time and space. In particular they cannot assess the extent to which the specific dynamics that they predict for a given system can affect the evolution of its population. This is a serious problem since the local dynamics of a given region of the disc, in particular the local value of encounter velocities, should directly control both the collision rates and their outcomes.

• 2) No feedback of collisions on the dynamics. Impacts, be they accreting, bouncing, or fragmenting ones, necessarily alter the orbits of colliding objects, because angular momentum is redistributed and kinetic energy is lost to heat. These effects can significantly change the results of pure N-body simulations by damping the high eccentricities induced by a planetary or stellar perturbations, or by ejecting collisional fragments far from their progenitors.

• 3) Particle lifetimes / spatial structures. Dynamical processes have characteristic timescales, which can be relatively long, as for instance in the case of either mean motion resonances or secular effects. Depending on the disc’s surface density, collision timescales could be shorter than the dynamical ones, thus interfering with, or even preventing the development of the corresponding spatial structures. This problem gets even more complex because collisional timescales are in most cases size-dependent. For the smallest grains, which are of specific interest because they usually dominate the disc luminosity (Thébault & Augereau 2007), the dynamical timescales also become size-dependent owing to the effect of radiation pressure. These competing timescale effects, in addition to affecting the development of spatial structures, can have direct consequences on the size distribution, which can, at some locations in the disc, strongly depart from any classical equilibrium power-law.

Incorporating collisions into an N-body scheme is extremely difficult in the context of debris discs, where typical impact velocities are large and lead to destructive collisions producing a vast amount of small fragments. Monitoring the evolution of all these fragments would lead to an exponential increase in the number of particles considered and a simulation that would be unmanageable.

The collisional evolution of debris discs is usually investigated separately with a different type of model, based on statistical particle-in-a-box codes, where the solid body population is distributed into logarithmically-spaced size bins, whose number density is followed using detailed prescriptions for collision outcomes (fragmentation, cratering, accretion). The price to pay for this approach is its low spatial resolution and very limited modelling of the dynamics: collision rates are estimated using spatially averaged estimates of orbital parameters that are fixed or follow very simplified evolution laws (e.g. Krivov et al. 2006; Thébault & Augereau 2007). A purely statistical approach thus cannot provide any precise information about the creation and evolution of spatial structures. For this, an N-body scheme is unavoidable.

1.1. N-body and collisions: previous studies

Given the sheer difficulty of incorporating fragmenting collisions into N-body codes, there have only been a few attempts at this task. The conceptually most natural, but in practice most difficult way to do this is by a “brute force” method that follows the fragments produced after each impact. The pioneering study of that method is that of Beauge & Aarseth (1990), whose code followed 4 fragments after each shattering collision, but had to be restricted to only 200 initial bodies, and was limited to relatively short timescales before reaching a critical number of particles. More recently, Leinhardt & Richardson (2005) considered the breakup of rubble piles of gravitationally bound hard spheres, but did not study sizes smaller than the initially defined “hard” sphere units.

Another solution is to use a mix of the N-body and the statistical approaches. To our knowledge, the only published attempt in this direction is that of Grigorieva et al. (2007), using a population of “super particles” (SPs), whose orbits are deterministically followed, each representing a cloud of monosize grains. When SP paths cross, collisions are considered between their respective grain populations in a statistical way, and new SPs are created to account for the fragments. However, this model was restricted to very short timescales of a few orbital periods. This combined N-body/statistical approach is thus far from being able to address long-term phenomenon, but it is probably the most promising one in terms of its potentialities in solving all three main problems listed in the previous section.

As of today, the most sophisticated available model is probably that of Stark & Kuchner (2009), whose main ambition is to address shortcoming number 3¹ regarding the competing effects of dynamical and collisional lifetimes on the development of spatial structures. Its principle is to release a population of collisionless test particles and record, at regularly spaced time intervals, their position to construct a stream of positions and velocities for each particle. All streams are then recombined to create a surface density map on which the same particles are released once again and removed following collision probabilities derived from the density map. A new map of streams is then created and the process is iterated until it converges. This pioneering model has provided impressive results for the specific case of the Kuiper Belt (Kuchner & Stark 2010). In its current version, this code is designed for the specific case of one perturber on a circular orbit. Although the case of multi-perturbers and of an eccentric perturber could in principle be implemented, these issues have not yet been addressed and the code might not be able to handle rapidly evolving spatial structures (Stark, priv. comm.). In the case of an eccentric perturber, this code requires that the time between two consecutive records must equal a multiple of the orbital timescale, which must be less than the collision timescale. As a result, this code may be unable to model highly collisional discs with distant perturbers on eccentric orbits, as in the case of a circumprimary debris disc in a binary.

We present here our new code, which is also designed to study how collisional lifetimes affect the development of perturbed spatial structures (shortcoming number 3 in the list given in Sect. 1), which can be used in the generic case of a collisional disc perturbed by one massive body, planet, or stellar companion, with any possible orbit, circular or eccentric, internal to, embedded in, or external to the disc. This code can handle size-dependent effects, in particular for small grains placed on high-eccentricity orbits by radiation pressure. The only required assumption is that a dynamical and collisional steady state has been reached in the system. We illustrate this new model with the specific case of a circumprimary debris disc perturbed by an external stellar companion.

2. Model

Fig. 1 Schematic presentation of the numerical model. nsav = 10 “collisional” runs are performed, each taking as a starting point one of the nsav parent body (PB) disc configurations, corresponding to nsav different positions of the companion star on its orbit separated by a constant time interval Δtsav = torb/nsav, obtained at the end (i.e., steady state) of the parent body runs. Nnum = 2 × 105 particles are released from the PB’s positions, following a size distribution dN ∝ sqds down to the radiation-pressure blow-out size sRP. Particles positions are recorded at time intervals separated by Δtsav. Each simulation is run until all particles have been removed, either by dynamical encounters with the companion star or by collisions. The collected data is then used to reconstruct the dust distribution at dynamical and collisional steady state by means of the following procedure (indicated by the red arrows on the graph): the dust distribution at a given time t0, corresponding to a given position iφ of the companion star, is the combination of the dust released at t0 for the iφ run, plus the dust particles released at t0 − Δtsav for the iφ − 1 run that have not been removed at t0, plus the dust released at t0 − 2Δtsav for the iφ − 2 run that has not been removed at t0, etc. The procedure is iterated until we reach the jfinal record (particles released at t0 − jfinalΔtsav), for which no particle has survived until t0. Given that the binary orbit is divided into nsav = 10 positions, all iφ − 10 × j runs correspond to the iφ one. (The different disc profiles displayed in the figure are not simulation results but only illustrative sketches.)

We assume that the studied system has reached a steady state, both dynamically and collisionally. We consider the forces of gravity and radiation pressure from the central star, and the gravitational pull of a perturbing planet and/or stellar companion. Collisions are assumed to produce bodies following a size distribution of the form dN ∝ s^qds. We divide particles into two categories: parent bodies (PB), i.e., bodies that are large enough not to be affected by radiation pressure, and small fragments steadily produced from these parent bodies by collisions and placed on eccentric or unbound orbits by radiation pressure. The magnitude of the radiation pressure force is parameterized, as usual, by the parameter β, corresponding to the size-dependent ratio of this force to that of star’s gravity. The dynamical evolution of the parent bodies and the smaller grains is computed using an adaptive-step fourth order Runge-Kutta integrator.

For the sake of simplicity, we consider a reference case of particles orbiting a central star and one external perturber on an eccentric orbit. Our numerical procedure is divided into three steps that can be schematically described as follows.

2.1. Parent body runs

A first run is performed with only the parent bodies and thus only the central and perturbing gravitational forces. This run is stopped once a dynamical steady state has been reached, i.e., when no changes are observed between two consecutive passages of the perturber at the same orbital phase. In this steady state, the position of the parent body disc is then recorded for a series of n_sav different positions of the perturber on its orbit, separated by a fixed time interval Δt_sav = t_Porb/n_sav (where t_Porb is the orbital period of the perturber)².

2.2. Collisional runs

A series of “collisional” runs are then performed, each taking as an initial condition the positions of the parent bodies for one of the n_sav phases of the perturber, as denoted by the number i_φ, where 1 ≤ i_φ ≤ n_sav, from the sample of recorded PB configurations. For each collisional run, at t = 0, N_num particles are released following a size distribution given by dN ∝ s^qds, where q = −3.5 (Dohnanyi 1969). The positions of these particles are then recorded, at regularly space time intervals Δt_sav corresponding to the same sample of orbital phases of the perturber as in the PB run.

At each integration timestep, each of the released particles is assigned a collisional destruction probability f_Dcoll depending on its size s, velocity, and the local geometrical optical depth τ_r,θ of the system: (1)where dt is the time increment of the code, t_orb the orbital period of a β = 0 body at this location of the disc, s₀ is a reference size, and δV and δV₀ are the departures from the local Keplerian velocity V_Kep of the considered particle and for a parent body with β = 0, respectively (see Thébault & Augereau 2005, for more details of the procedure). The optical depth values τ_(r,θ) are derived from the n_sav parent body runs, under the assumption that collisions occur mainly in the parent body region. This assumption agrees with previous debris disc modelling results (e.g. Krivov et al. 2006; Thébault & Augereau 2007) and has been adopted in all studies investigating the outer regions of debris discs (e.g. Strubbe & Chiang 2006; Thébault & Wu 2008). In practice, we record the (steady-state) final positions of all parent body runs and transform them into (r,θ) maps of the vertical optical depth. These maps only give relative values and have to be normalized by ⟨τ⟩, the average normal optical depth assumed for the PB disc. The collisional activity can thus be tuned by setting the initial value of ⟨τ⟩. We note that this procedure requires that dt is shorter than the collision time for the grain to properly resolve the collisional destruction probability, a criteria that is always met for the study of perturbed debris discs, where dt is a fraction of the orbital period and t_coll is always  ≥  t_orb.

The collisional runs were stopped when all particles had been eliminated, either by collisions or by dynamical ejection after an encounter with the perturber. At the end of this step, to each initial phase i_φ of the perturber we assigned a series of surface density maps σ′(i_φ,0),σ′(i_φ,1),...,σ′(i_φ,j), recorded at the same regularly spaced time intervals Δt_sav, following the fate of all particles released when the perturber was at the i_φ position. We note that all snapshots σ′(i_φ,j), σ′(i_φ,j + n_sav), σ′(i_φ,j + 2n_sav), etc., correspond to successive passages of the perturber at the same orbital phase (each perturber orbit being divided into n_sav positions).

2.3. Recombining

In the final stage, we use all data collected in step 2 to reconstruct the dust distribution, at steady state, for each possible orbital phase of the perturber. The principle is that, at a given time t_iφ corresponding to the perturber’s position i_φ, the total dust population regroups grains that have just been produced (for the present position i_φ of the perturber), as well as survivor grains that have been produced earlier (when the perturber was at a different orbital position), and that have not yet been collisionally destroyed or dynamically ejected. The procedure is the following: we start with the most recent dust particles, produced at t_iφ, whose spatial distribution is given by the saved record σ′(i_φ,0). We then add the previous generation, produced at t_iφ − Δt_sav when the perturber was at an angular position index i_φ − 1, whose spatial distribution at time t_iφ (i.e., at time Δt_sav after their release) is given by the saved record σ′(i_φ − 1,1). This procedure is then iterated, working our way back in time and piling up all the surviving grains from the successive records σ′(i_φ − j,j), to produce the total geometrical optical depth at time t_iφ(2)The index j_max corresponds to the most ancient record, for which all initially released particles have been eliminated after the time j_max × Δt_sav differentiating their release from the present time. We note that all the snapshots σ′(i_φ − j,j), σ′(i_φ − j − n_sav,j + n_sav), σ′(i_φ − j − 2n_sav,j + 2n_sav), etc., are taken from the same source collisional run, with the perturber position at release being (i_φ − j), but separated by an integer number of orbital periods. In other words, the σ′(i_φ − j − kn_sav,j + kn_sav) record corresponds in practice to the σ′(i_φ − j,j + kn_sav) one.

An illustrated summary of the main steps of our numerical procedure can be found in Fig. 1.

2.4. Tests

Fig. 2 Test runs: left panel: azimutally averaged radial profile of the normalized surface density for a fiducial case with no companion. The expected theoretical asymptotic behaviour in r-1.5 is plotted as a reference. The initial extension of the parent body ring is set to r = 1. Right panel: azimutally averaged radial profile of the normalized surface density for a companion of eccentricity eb = 0.5 as compared to the results of Thébault et al. (2010) with a different code. For this with-companion run, distances have been normalized to the theoretical limit for orbital stability rcrit (see Sect. 3.1).

As there is no analytical solution to the full problem of a collisionally active disc perturbed by a massive body, we follow an approach similar to that of Stark & Kuchner (2009) and first check the correctness of our algorithm’s solution for a simplified case with no companion. In this case, Strubbe & Chiang (2006) and Thébault & Wu (2008) demonstrated that, for a collisional debris ring whose parent body distribution has a sharp outer edge r_out, the radial profile of the surface density beyond r_out should tend towards an aymptotic slope in r^-1.5 (owing to the small grains produced within the ring that radiation pressure has placed on eccentric orbits entering these outer regions). We note that we do not try to reproduce the numerical simulations of Strubbe & Chiang (2006) or Thébault & Wu (2008), which addressed slightly different issues, but the robust conclusions of their analytical derivations that a collisional ring always tends towards a profile with a slope in −1.5 outside the collisional active region. We thus run a test simulation with no companion and a parent body ring having a sharp outer edge at r_out = 1 (arbitrary unit). As can be clearly seen in Fig. 2a, we obtain a profile that tends towards the standard slope in −1.5 beyond r ~ 1.5. We note that this test is not a validation of the whole collisional procedure, as the −1.5 slope is a result that depends only weakly on the collisional rate within the PB ring. However, it validates the recombination procedure, which is rather complex, piling up maps after maps coming from different runs and at different times, none of these individual maps having a radial profile tending towards −1.5.

In the case of a massive external companion, we perform an additional test by comparing our results to those obtained by Thébault et al. (2010) who, using a completely different algorithm with only one dimensional (1D) resolution, estimated the azimuthally averaged radial profile of the surface density. We take as a reference the case of a M₂ = 0.25M₁ companion on a e_B = 0.5 eccentric orbit and a parent body ring of optical depth τ = 2 × 10^-3 extending out to the orbital stability limit r_crit (see Sect. 3.1). We then azimuthally average the spatial distribution and display it in Fig. 2. As can be seen in Fig. 2b, the obtained averaged radial profile is a very good match to the one obtained by Thébault et al. (2010) with their different method (see Fig. 2 of that paper), i.e. a steeper profile than in the previous no-perturber case owing to the steady removal of high-β particles by the companion. This second test is a robust validation of the collisional and dynamical prescriptions, as the shape of the radial profile depends directly on the balance between the collisional activity in the PB disc and the dynamical ejection of particles perturbed by the companion.

3. A case study: debris disc in a binary

To illustrate the potential of our code, we consider a case study of a circumprimary debris disc in a binary system, of particular interest in light of the increasing amount of research regarding planet formation in binaries (for a review, see Thébault 2011). This case has been investigated in several recent studies, both observational (Trilling et al. 2007; Plavchan et al. 2009; Duchêne 2010) and theoretical (Thébault et al. 2010). The main issue investigated has been the extent of circumprimary discs, in particular whether the companion star can induce a truncation that can be detectable when looking at either the infra-red excesses or the radial profile of the resolved disc. The numerical study of Thébault et al. (2010) showed that the coupling of collisional activity and radiation pressure plays a crucial role, steadily placing small dust grains in regions that are in principle dynamically unstable. Since these grains dominate the flux at all wavelengths up to mid-IR, debris discs can thus appear to extend far beyond the theoretical radial distance r_crit for orbital stability around the primary. However, Thébault et al. (2010) used a collisional code that had only 1-D spatial resolution (for which all azimuthal information had therefore been lost in phase averaging), and were thus unable to study how binarity affects the shape of circumstellar discs. Such a study would be invaluable, as the presence of a companion star has been invoked as a possible explanation of the characteristics of several systems, such as HR4796 (Augereau et al. 1999; Schneider et al. 2009) or HD 141569 (Augereau & Papaloizou 2004; Quillen et al. 2005; Ardila et al. 2005).

3.1. Setup

We consider a binary of separation a_B, eccentricity e_B and mass ratio μ. We make the not-so-unreasonable assumption that the disc of large parent bodies has been shaped and truncated by the companion star. To ensure that this is true, we allow the initial disc to extend slightly beyond the empirical (1-D) outer limit r_crit necessary for orbital stability derived by Holman & Wiegert (1999), and leave the code naturally remove all unstable particles.

To allow the straightforward comparison of different cases, all distances are normalized so that r_crit = 1. In order to avoid the huge computational cost of calculating orbits very close to the primary, and since these regions are in any case the least affected by binarity, we consider a ring-like configuration for the initial disc: 0.6 ≤ r ≤ 1.1 (in units of r_crit).

For the binary, we consider a fixed mass ratio μ = 0.25, corresponding to estimates of the mean mass ratio in binaries derived from extensive surveys (Kroupa et al. 1990; Duquennoy & Mayor 1991), and explore four different orbital configurations: e_B = 0, e_B = 0.2, e_B = 0.5, and e_B = 0.75 (the values of a_B are then automatically obtained from the requirement that r_crit = 1).

We consider an average vertical optical depth of ⟨τ⟩ = 2 × 10^-3, which is typical of dense debris discs such as β Pictoris.

For the parent body runs, all particles start on circular orbits (e = 0) in the binary’s midplane (i = 0). For the collisional runs, particles are released from the parent bodies seeds following a size distribution of the form dN(s) ∝ s^-3.5ds. Particle sizes are parameterized by the value of their β parameter. Sizes are distributed continuously between β_max = 0.5 and β_min = 0.012. We take the minimum particle size to be β_max = 0.5, corresponding to the smallest grains in bound orbits. The contribution of smaller grains, below the blow-out size limit, is negligible for systems with optical depths in the ⟨τ⟩ ~ 10^-3 range (see Fig. 3 of Thébault et al. 2010). The number of sampled positions within a binary orbit is n_sav = 10. Convergence test runs have shown that this is the optimal value: larger values do not lead to significant changes in the final recombined density maps, while reducing the value of n_sav leads to increasingly divergent results. We note that for the first binary orbit of the collisional run, the sample value is higher, n_sav = 100, to enable us to accurately track the rapid initial blow-out of high-β grains. All initial parameters are summarized in Table 1.

3.2. Results

Fig. 3 Normalized vertical optical depth, at collisional and dynamical steady state, for four binary configurations; from top to bottom: eB = 0, eB = 0.2, eB = 0.5, and eB = 0.75. For each case, the disc profile is shown at periastron (left) and apoastron (right) passages of the companion star. The dotted circle is the theoretical outer radius for orbital stability rcrit. The rectangle in the upper right corner shows the orbit of the binary relative to the rcrit radius (dotted circle) and the current position of the companion on its orbit. (An animated version of these results, displaying disc profiles for ten different orbital positions of the companion, can be found at http://lesia.obspm.fr/perso/philippe-thebault/bindeb.html)

Depending on the considered configuration, the parent body simulations are run for between 200 (e_B = 0 case) and 500 (e_B = 0.75) orbital periods of the binary in order to reach a steady state. This steady state is reached once all particles on unstable orbits have been removed and once all transient dynamical features have disappeared. As shown by Augereau & Papaloizou (2004), the most notable of these transient features are spiral arms, which develop owing to the sudden introduction of a perturbing body, and disappear on the scale of a few hundred orbital periods. We then perform the collisional runs following the procedure presented in Sect. 2.2. The final disc profiles, after recombination of all collisional runs following the method described in Sect. 2.3, are displayed in Fig. 3.

A first important point is that these results confirm the main conclusion of Thébault et al. (2010) that circumprimary debris discs extend far beyond the outer limit for orbital stability. For all four binary configurations the regions beyond r_crit = 1 are populated by small grains steadily produced by collisions in the parent body ring and placed on eccentric orbits by radiation pressure. The level of dustiness in these “forbidden” regions increases with e_B because, for highly eccentric binaries, the perturbing companion is far from the main ring between two passages at periastron, so that small grains placed by radiation pressure in the r ≥ r_crit region have a longer dynamical survival timescale before being removed.

Fig. 4 Steady-state surface density profiles for the same (from top to bottom) eb = 0, eb = 0.2, eb = 0.5 and eb = 0.75 cases as in Fig. 3. Left panels: radial profiles, along the binary semi-major axis direction, at periastron (full line) and apoastron (dotted line) passages of the companion star. Right panels: azimutal profiles, for a Δr = 0.1 wide annulus centred around r = 1 (in rcrit units), at periastron (full line) and apoastron (dotted line) passages of the companion star.

Regarding the spatial structure of the disc, we find that, for low e_B values, it consists of precessing spiral features directed towards the companion’s position that develop from the main ring³. These precessing spirals are the most pronounced for the e_B = 0 case, being clearly visible for the whole binary orbit with no significant shape evolution, as logically expected for this circular orbit case. As e_B increases, these spiral structures become fainter. For the e_B = 0.2 case, they are only visible during half of the binary orbit (from periastron to apoastron), while they only briefly appear at periastron passages in the e_B = 0.5 run. For the e_B = 0.75 case, spirals are no longer visible. This disappearance of the precessing spirals for high e_B reflects the more general trend of the disc towards having a time-invariant structure. In particular, no significant shape change is observed at different orbital phases of the binary. The disc assumes a fixed and strongly asymmetric shape, extending much further out in the apoastron direction than in the periastron one (see radial profiles in Fig. 4e and g).

Another important result is the presence of strong inhomogeneities in the main parent body ring itself. For low e_B cases, they consist of a pronounced dip in the opposite direction of the companion’s position (Fig. 4b), which precesses at the binary’s angular velocity. In the moderately eccentric case (e_B = 0.2), this dip is surrounded by two overdensitites, precessing at the same rate and symmetrically positioned at  ±70−90 degrees with respect to the companion’s antipolar position (Fig. 4d). The density contrast around the density dip can reach up to  ~60% (Fig. 4b). For the e_B = 0.2 case, in addition to the precessing features, a time-invariant structure appears, i.e., a dip in the binary orbit’s periastron direction (Fig. 4d). As e_B increases even more (e_B ≥ 0.5), this time-invariant asymmetry becomes more pronounced, while the precessing structures progressively fade away (Fig. 4f and h).

3.3. Discussion

The results of the previous section show that a companion star can significantly affect the shape of a circumprimary disc, inducing pronounced radial and azimutal asymmetries. Basically, two different regimes can be identified depending on the binary’s eccentricity.

For low e_B binaries, the disc’s shape varies with time and precesses with the binary’s orbit. Strong asymmetries are observed, whose position is always related to that of the companion on its orbit. Strong asymmetries are also observed for highly eccentric binaries, but they are in this case almost time-invariant and the disc has roughly the same shape regardless of the location of the companion on its orbit. In these high e_B cases, the main asymmetries are oriented with respect to the binary orbit’s geometry, not with the actual position of the companion on the orbit. Between these two regimes, there is an intermediate state (see the e_B = 0.5 runs) where the disc has a globally invariant asymmetric structure with transient features close to the companion’s periastron passages.

These fixed, precessing or transient features are caused by the combination of three distinct processes: 1) collisional activity within the PB ring, steadily producing vast quantities of small grains carrying a large fraction of the geometrical cross section; 2) radiation pressure on the these small fragments, which places them on high-eccentricity orbits; and 3) gravitational perturbations by the companion star that eventually, but not immediately, remove all grains beyond r_crit. Equilibrium between these three effects results in a steady-state where a large fraction of the small, high-β fragments produced in the PB ring remain in the dynamically forbidden regions long enough to cause the appearance of pronounced structures, whose shape and evolution strongly depends on e_B.

A careful examination of the variety of structures obtained for the different cases displayed in Fig. 3 shows that the action of a companion star could be a potential explanation of some features that have been routinely observed in debris discs: spiral arms, azimuthal inhomogeneties in ring-like structures, blobs, etc. We stress that all of these structures are obtained at steady-state, and thus do not appear during a brief specific period in the system’s history, as would be the case for a fly-by with a passing star (Ardila et al. 2005; Reche et al. 2009).

This explanation should of course be considered with some caution. One reason is that other mechanisms have already been found to be able to produce asymmetries and inhomogeneities in discs, the most commonly invoked being the gravitational pull of a planet (Krivov 2010). Another important point is that, in contrast to the planet-perturbation scenario for which an undetected planet can often be freely postulated given the observational constraints on the system, the presence or absence of a companion star is usually far more tightly constrained. Of all the 25 resolved debris discs⁴, only three are known to inhabit a confirmed binary system, namely HR4796A, HD 181296 and ζ Ret (HD 141569A is more a “transitional” disc than a bona fide debris disc according to the definition given by Lagrange et al. 2000). A detailed analysis of these systems exceeds the scope of the present numerical paper, as it would require exploring several additional free parameters, such as τ or the mass ration μ. This analysis will be undertaken in a forthcoming study.

4. Limitations and perspectives

We note that this code is only a first, albeit important, step towards a fully integrated dynamics + collisions model. We now list its main limitations and the possible improvements that can be expected:

• The most obvious limitation is that collisions are notmodelled in a self-consistent way, but through an empiricallydefined collision destruction probability assigned to eachparticle, and that no feedback of the collisions on thedynamics is included. These two problematic issues canprobably not be handled by an N-body approach, for which the present code is maybe the upper limit to what can be achieved, at least for the case of high-velocity fragment-producing collisions. In this case of fragmenting impacts, a fully integrated and self-consistent dynamics+collisions model is still beyond our reach today, and can probably only be achieved by the next generation of “hybrid” codes in the spirit of the pioneering model explored by Grigorieva et al. (2007) .

• Our procedure is suited to systems where a level of symmetry has been reached, where the location of the parent bodies is symmetric with respect to the perturber’s orbit, i.e., the global shape of the PB disc is identical between two passages of the perturber at the same orbital position. We note that this symmetry is not required for each individual PB’s orbit, which would be impossible because of the different periodicities between PBs and the perturber. What is required is a global symmetry for orbital streamlines, provided that each parent body’s orbit is populated by a large number of bodies distributed uniformly in mean anomaly. This assumption of a phasing between global PB locations and pertuber orbit might not be valid for all perturbed disc configurations, in particular for systems strongly shaped by mean motion resonances. This case, which usually occurs for discs interacting with embedded planets, will be explored in a forthcoming study. We note however that, for the present case of a circumprimary disc in a binary, this global symmetry is neither assumed nor postulated but is a verified behaviour of the PB disc, which always reaches this state after a transitory period of 100–300 binary periods (see Sect.  3.2).

• The equation governing the collisional probability does not resolve vertical structure in the disc since it is based on the vertically integrated optical depth. This restrics the current version of the model to two dimensional problems, but those have been by far the most widely investigated in all past studies. A three-dimensional version of the collisional probability prescription is however fully implementable, at the cost of an increase in the size of the saved density maps at the end of the parent body runs.

• Our model assumes that all collisions take place in the region of the parent body belt. This is not a major limitation for the present case of highly collisional massive debris discs, but the model does not perform well for fainter, drag-dominated discs where collisions can occur well within the belt. Our code is thus unable to include significant drag effects in its current version.

5. Conclusions

We have developed a new code to study the spatial structure of dynamically perturbed debris discs (be it by planets or companion stars) taking into account the effect of collisions. It is especially designed to quantify the competing effects of the steady collisional production of small grains, which have been placed by radiation pressure in dynamically unstable regions, and their removal by gravitational perturbations. The principle of our numerical method is the following: we assume that the system has reached a steady state and we perform a series of distinct N-body runs, each corresponding to a different initial position of the perturber in its orbit, for which each particle has a collisional destruction probability depending on its size and location. For each run, all particle positions are recorded at regularly spaced time intervals. This database of “collisional” runs is then used to reconstruct the steady state profile of the disc at different orbital phases of the perturber.

To illustrate the potential of this numerical procedure, we have applied it to the case of a circumprimary debris disc in a binary. We have demonstrated that the complex interplay between collisional activity, radiation pressure, and gravitational perturbations can create pronounced structures inside and outside the dynamical stability region. Depending on the binary’s eccentricity, two regimes can be distinguished. For low e_B, the disc’s structure is time varying, with spiral arms forming in the dynamically “forbidden” region beyond r_crit and precessing at the binary’s angular velocity. For high e_B, spiral structures disappear, the disc adopts a time-invariant structure, extends far outside the stability region in the binary’s apoastron direction, and has a pronounced material depletion, inside r_crit, in the periastron direction.

Our model has the potential to be applied to many other configurations, in particular debris discs sculpted by planets, which will be the purpose of a future study.

Acknowledgments

We thank the reviewer, Christopher Stark, for very useful comments that helped improve the quality of the manuscript. We thank the French National Research Agency (ANR) for financial support through contract ANR-2010 BLAN-0505-01 (EXOZODI).

References

All Tables

All Figures

Fig. 1 Schematic presentation of the numerical model. nsav = 10 “collisional” runs are performed, each taking as a starting point one of the nsav parent body (PB) disc configurations, corresponding to nsav different positions of the companion star on its orbit separated by a constant time interval Δtsav = torb/nsav, obtained at the end (i.e., steady state) of the parent body runs. Nnum = 2 × 105 particles are released from the PB’s positions, following a size distribution dN ∝ sqds down to the radiation-pressure blow-out size sRP. Particles positions are recorded at time intervals separated by Δtsav. Each simulation is run until all particles have been removed, either by dynamical encounters with the companion star or by collisions. The collected data is then used to reconstruct the dust distribution at dynamical and collisional steady state by means of the following procedure (indicated by the red arrows on the graph): the dust distribution at a given time t0, corresponding to a given position iφ of the companion star, is the combination of the dust released at t0 for the iφ run, plus the dust particles released at t0 − Δtsav for the iφ − 1 run that have not been removed at t0, plus the dust released at t0 − 2Δtsav for the iφ − 2 run that has not been removed at t0, etc. The procedure is iterated until we reach the jfinal record (particles released at t0 − jfinalΔtsav), for which no particle has survived until t0. Given that the binary orbit is divided into nsav = 10 positions, all iφ − 10 × j runs correspond to the iφ one. (The different disc profiles displayed in the figure are not simulation results but only illustrative sketches.)

In the text

Fig. 2 Test runs: left panel: azimutally averaged radial profile of the normalized surface density for a fiducial case with no companion. The expected theoretical asymptotic behaviour in r-1.5 is plotted as a reference. The initial extension of the parent body ring is set to r = 1. Right panel: azimutally averaged radial profile of the normalized surface density for a companion of eccentricity eb = 0.5 as compared to the results of Thébault et al. (2010) with a different code. For this with-companion run, distances have been normalized to the theoretical limit for orbital stability rcrit (see Sect. 3.1).

In the text

Fig. 3 Normalized vertical optical depth, at collisional and dynamical steady state, for four binary configurations; from top to bottom: eB = 0, eB = 0.2, eB = 0.5, and eB = 0.75. For each case, the disc profile is shown at periastron (left) and apoastron (right) passages of the companion star. The dotted circle is the theoretical outer radius for orbital stability rcrit. The rectangle in the upper right corner shows the orbit of the binary relative to the rcrit radius (dotted circle) and the current position of the companion on its orbit. (An animated version of these results, displaying disc profiles for ten different orbital positions of the companion, can be found at http://lesia.obspm.fr/perso/philippe-thebault/bindeb.html)

In the text

Fig. 4 Steady-state surface density profiles for the same (from top to bottom) eb = 0, eb = 0.2, eb = 0.5 and eb = 0.75 cases as in Fig. 3. Left panels: radial profiles, along the binary semi-major axis direction, at periastron (full line) and apoastron (dotted line) passages of the companion star. Right panels: azimutal profiles, for a Δr = 0.1 wide annulus centred around r = 1 (in rcrit units), at periastron (full line) and apoastron (dotted line) passages of the companion star.

In the text