Dynamical evolution of an eccentric planet and a less massive debris disc
Abstract
We investigate the interaction between an eccentric planet and a less massive external debris disc. This scenario could occur after planetplanet scattering or merging events. We characterise the evolution over a wide range of initial conditions, using a suite of nbody integrations combined with theory. Planets near the disc midplane remove the inner debris region, and surviving particles form an eccentric disc apsidally aligned with the planet. The inner disc edge is elliptical and lies just beyond the planet’s orbit. Moderately inclined planets ( for ) may instead sculpt debris into a bellshaped structure enveloping the planet’s orbit. Finally some highly inclined planets () can maintain a disc orthogonal to the planet’s plane. In all cases disc particles undergo rapid evolution, whilst the overall structures evolve more slowly. The shapes of these structures and their density profiles are characterised. The width of the chaotic zone around the planet’s orbit is derived in the coplanar case using eccentric Hill radius arguments. This zone is cleared within approximately ten secular or diffusion times (whichever is longer), and debris assumes its final shape within a few secular times. We quantify the planet’s migration and show it will almost always be small in this mass regime. Our results may be used to characterise unseen eccentric planets using observed debris features.
keywords:
planetdisc interactions  planets and satellites: dynamical evolution and stability  circumstellar matter1 Introduction
Extrasolar planets exhibit a broad range of eccentricities, including values far higher than those typical of present day Solar System planets^{1}^{1}1See the exoplanet.eu database. However planets should form on circular orbits as eccentricity excitations are damped by the protoplanetary disc (Lissauer, 1993), so high eccentricities are generally attributed to postformation dynamical processes (Tremaine & Zakamska, 2004). Solar System dynamical theories suggest that the giant planets underwent significant migration and orbital evolution at early times (Tsiganis et al., 2005; Walsh et al., 2011), and such models could apply to other star systems too. During such violent periods eccentricity driving interactions, such as planetplanet scattering (Ford & Rasio, 2008) and mergers (Lin & Ida, 1997), could be commonplace. Stellar flybys (Malmberg & Davies, 2009), planetdisc interactions (Bitsch et al., 2013) and longterm interactions between planets (Wu & Lithwick, 2011) could also drive up planet eccentricities. So with large eccentricity changes possible, planets could come into contact with bodies previously too distant for interactions to occur.
In this paper, we investigate the scenario where such an eccentric planet encounters a debris disc. The latter describes a circumstellar disc of solids, containing bodies ranging in size from submicrometer dust grains all the way up to dwarf planets, and which is probably composed of material left over from the protoplanetary disc (see Wyatt 2008 for a comprehensive review). The Solar System’s asteroid and Kuiper belts fit this description. Debris material forms a collisional cascade, such that the collisional breakup of objects produces ever smaller bodies which are eventually blown out of the system by radiation pressure (Backman & Paresce, 1993; Wyatt & Dent, 2002). The dust reemits starlight at longer wavelengths, and of order a thousand extrasolar debris discs have been detected as mid to farinfrared excesses (e.g. Patel, Metchev & Heinze 2014). Observed discs typically contain Earth masses () in emitting grains (Wyatt, Dent & Greaves, 2003), and extrapolating these to include larger unseen bodies yields total masses several orders of magnitude greater (e.g. Wyatt & Dent 2002). The discs are observed at a few to hundreds of au from their host stars (Najita & Williams, 2005); for context, the total mass of the early Kuiper belt is thought to have been tens of spread between 30 and 50 au (Stern & Colwell, 1997).
Many debris discs have been resolved and show eccentricity and/or density asymmetries, which are often attributed to perturbations by unseen planets (e.g. Kalas, Graham & Clampin 2005; Krist et al. 2012; Schneider et al. 2009). The effect of a low eccentricity, low inclination planet on a cold debris disc is well studied, and invoking such planets successfully explains observed debris disc features in several systems such as offsets, warps and spirals (e.g. Wyatt et al. 1999; Augereau et al. 2001; Wyatt 2005). However discs with large asymmetries (e.g. Eiroa et al. 2010) or substantial misalignments (e.g. Kennedy et al. 2012) have also been observed which, combined with the discovery of Fomalhaut b (a highly eccentric and possibly disc crossing object: Kalas et al. 2013; Beust et al. 2014), motivate the study of more extreme interactions. Studies involving a highly eccentric planet have recently been conducted by Beust et al. (2014), Tamayo (2014) and Faramaz et al. (2014), specifically attempting to explain the discs of Fomalhaut and Reticuli. We wish to characterise this interaction generally, to produce quantitative predictions applicable over a broad region of parameter space. This includes planets with large initial inclinations and those with nonnegligible orbital evolution. In this paper we only examine planets which are more massive than the disc.
The layout of this paper is as follows. We describe our simulations in Section 2, and outline the three primary outcomes in Section 3. In Section 4 we provide an overview of dynamical theory, and apply this in Section 5 to make predictions about the postinteraction debris structure, planet evolution and important timescales. We discuss our results in Section 6 and conclude in Section 7. We use the subscript “plt” to denote parameters of the planet throughout the paper, and colour versions of all plots are available in the online version.
2 Simulations
We modelled the interaction using the Mercury 6.2 nbody integrator (Chambers, 1999), with a hybrid symplectic / BulirschStoer algorithm. This allowed the rapid computation of distant interactions whilst maintaining accuracy in close encounters. Our systems consisted of a massive particle representing the planet and to massive test particles to model the disc. The planet therefore perturbed the disc and viceversa however the disc had no selfgravity, which is valid as the more massive planet would dominate over debrisdebris interactions (e.g. Beust et al. 2014). Relatively small numbers of particles let us run many integrations, and identical systems simulated with and particles resulted in the same outcome in terms of disc and planet evolution. Furthermore the planets did not evolve stochastically, so we are confident that our particle numbers are sufficiently high. We used a wide range of disc and planet masses ( and respectively), with the former 0.1 to 100 and the latter 1 to 1000 . The planet was at least ten times as massive as the disc. The stellar masses () were equal to a solar mass () in almost all integrations. Simulations typically lasted for years, long enough for the major evolution to have occurred (see Section 5.4).
The discs were what we might consider early Kuiperbelt analogues, typically with inner and outer edges at and au respectively. Each disc particle was assigned an initial semimajor axis within this range. The particles themselves were of equal mass, and semimajor axes were populated in such a way that the resulting surface density profile (if each particle had negligible eccentricity and inclination , the latter relative to the initial disc midplane) could be written as , where is the radial distance from the star and is the surface density index. We used , that of the Minimum Mass Solar Nebula (Hayashi, 1981). The discs were relatively dynamically cool before the interaction, with inclinations and eccentricities typically drawn from the ranges and respectively. Finally, the initial argument of pericentre , longitude of ascending node and mean anomaly of each disc particle were uniformly drawn between and (note that due to the initial symmetry of the disc, the planet’s ascending node was defined to lie along the axis without loss of generality).
The eccentric planets typically had pericentres at 5 au, the distance of Jupiter in our solar system and the location at which massive planets could be expected to form based on snow line arguments (e.g. Ida & Lin 2004). This is therefore a reasonable estimate for the pericentre of a body made eccentric by planetplanet interactions. Their apocentre distances were typically between the disc’s initial inner and outer edge radii, with eccentricities from . Half of our simulations had planets with orbits initially located in the disc midplane, and the other half had planets inclined by up to . The latter had initial arguments of pericentre (the orientation of the planet’s pericentre relative to the ascending node, where the orbit crosses the disc midplane) equal to 0, 45 or . The planets were prograde with respect to the disc particles in all simulations (i.e. ).
3 Results
We observed three outcomes in our simulations. Firstly if the planet’s orbit lay roughly in the disc midplane, the inner disc particles were ejected and remaining debris formed an elliptical disc apsidally aligned with the planet (outcome Ia). Secondly if the planet’s orbit brought it far out of the midplane then surviving debris instead formed a hollow bellshaped structure, again aligned with the planet. This structure completely enveloped the planet’s orbit, apart from for two holes at its pericentre and apocentre ends (outcome Ib). An extreme case of this outcome occurred for very highly inclined planets, where the holes opened to the extent that debris formed a disc orthogonal to the planet’s orbit. Finally a planet could eject all nonresonant disc particles if it came sufficiently close to the disc outer edge (outcome Ic). Tables of individual simulation parameters and their final outcomes are given in Appendix A, and we now describe these results in more detail.
3.1 Outcome Ia  the coherent disc
In this outcome most of the inner disc particles were ejected, and remaining material formed an elliptical disc apsidally aligned with the planet. We name this structure the “coherent disc”. The particles themselves covered a wide range of eccentricities and orientations, and these values oscillated whilst the disc as a whole maintained its shape. This outcome was noted by Faramaz et al. (2014). The longitude of pericentre (defined as ) of both the planet and disc precessed at the same rate, and these objects remained aligned. We plot the final state of an outcome Ia system (simulation 21) with an initially coplanar planet on Figure 1.
We plot the final semimajor axes, eccentricities, inclinations and orientations of the simulation 21 particles on Figure 2, and see that coherent disc particles occupy a relatively well defined region in versus space. In this simulation their semimajor axes all lie between about 50 au and the initial outer disc edge at 60 au, and these did not change significantly throughout the interaction. Their eccentricities have clearly increased, and oscillate between roughly zero and some maximum ( on Figure 2) which declines with increasing semimajor axis. These imply a secular interaction, a long term angular momentum exchange leading the eccentricity oscillating whilst the semimajor axis remains constant (a detailed discription of secular theory is given in Section 4.1). Not all particles at smaller semimajor axes were ejected; some were longterm stable if in resonance with the planet from the start of the simulation. These particles form the additional ring structures superimposed on the coherent disc on Figure 1.
The coherent discs are also well defined on the , plane, implying coupling between particle eccentricity and orientation. Each moves in a circle centred on a point between the location of the planet’s orbit on this plane and the origin, with those at larger semimajor axes forming tighter circles (i.e. their maximum eccentricities are smaller, as apparent from the versus plot). These circles all pass close to the origin. Each particle moves at a uniform rate anticlockwise about its circle, and the planet and the centre of the circles also rotate in this direction (far more slowly than individual disc particles), with each circle remaining aligned with the planet. Again, this is the behaviour expected of secular particles in the low eccentricity regime (e.g Wyatt et al. 1999). However one might not necessarily expect this behaviour if the eccentricities are large (see Section 4.1).
Beust et al. (2014) find similar behaviour in their analysis; our particles are most akin to those below the red line on the first panel of their figure 7, which either librate or circulate in whilst oscillating in eccentricity. We also note that resonant particles move in a similar manner, but with much smaller maximum eccentricities (the green points on Figure 2). Coherent disc particles also move in closed loops on the , plane, only not in circles but ovals and in a clockwise direction. These ovals are centred on the planet and aligned with . As evolves with time, the tracks of the particles on this plane rotate. Resonant particles also behave in this manner.
The paths of secular particles on the , and , planes are qualitatively the same when the planet is moderately inclined to the disc midplane. For , this means planets with and for and respectively, and all inclinations for (note that means the planet’s pericentre and apocentre both lie in the initial disc plane, whilst means its apocentre occurs at the point on its orbit which is highest above this plane). Inclined planets cause the coherent disc to align with the planet’s orbital plane, and this disc is puffed up compared to the coplanar regime. The discs also show vertical structure, appearing shaped when viewed along the planet’s line of apsides. The second row of Figure 4 shows an outcome Ia interaction where the planet was initially inclined by with (simulation 54).
Note that for very high planet inclinations with , particles still behave in the above manner but their inclinations may reach very high values. In this case debris forms an ovular structure enveloping the planet’s orbit, but which does not have holes at the ends nor the same density pattern observed in Ib simulations. Whilst such a structure may appear similar to the Ib case, secular particle behaviour is completely akin to the Ia regime. Thus this is an extreme case of outcome Ia, with particle planes still oscillating about that of the planet, only with inclinations high enough that the “opening angle” of the disc (as seen fronton on the second row of Figure 4) may be very high.
The planet evolution is qualitatively the same for all outcome Ia interactions, and that from simulation 21 is shown on Figure 3. The semimajor axis and eccentricity typically decrease by 0.01–1 and 1–10 per cent respectively, with the latter occurring rapidly at the beginning of the interaction. The inclination change for coplanar planets is negligible. Once the planet’s semimajor axis and eccentricity stop evolving its pericentre precesses at a constant rate. Inclined outcome Ia planets evolve as in the coplanar case, only now their inclinations rapidly decrease by 120 per cent at the start of the simulation.
3.2 Outcome Ib  the bellshaped structure
If the planet’s inclination is above some critical value with respect to the disc midplane (which is a function its argument of pericentre), the debris behaviour differs significantly from the previous regime. These critical inclinations were around 20 and for and respectively in our simulations with ; outcome Ib never occurred for . In outcome Ib secular particles no longer move about ovals on the , plane, but rather a squashed “kidney” shape centred on a point orthogonal to the planet’s orbital plane. Again the orientation of the shape is related to ; we find it remains parallel to a line at angle to the axis which passes through the planet’s location, and which rotates as the planet precesses. Particle behaviour on the , is more complicated; a number still move in circles or ellipses aligned with the planet, however most now move in much more complex shapes. Figure 5 shows the particles from simulation 56 on both planes, at an instant (30 Myr) when structure is most clearly defined. Higher values result in tighter particle paths on the , plane.
The physical debris structure resulting from the superposition of these particles is no longer an eccentric disc but rather a hollow bellshaped bubble enveloping the planet’s orbit, as shown on the third row of Figure 4. The structure is broadest around planet apocentre and is roughly symmetric about its line of apsides, appearing ovular rather than shaped when viewed along this line. There are holes in both ends of the shape near planet pericentre and apocentre, and if is increased then these holes become larger and the bell structure is squashed when viewed topdown or sideon. Therefore extremely high inclinations () may result in a stable debris disc lying orthogonal to the planet’s orbital plane (bottom row of Figure 4). Particles in all outcome Ib simulations behave in the same manner on the , plane as in outcome Ia, however far fewer particles are trapped in resonance in the former outcome.
The planet’s eccentricity and semimajor axis undergo only minor evolution similar to that on Figure 3, with the eccentricity again decreasing rapidly at the start of the simulation. In most cases the planet’s inclination decreases by less than 20 per cent, although this occurs over the timescale of the semimajor axis change (rather than the eccentricity change as for case Ia). The exceptions are simulations 57 and 58, the very highly inclined planets with stable discs orthogonal to their orbital planes. The inclinations of these planets undergo sinusoidal oscillations with periods of Myr, over which their values change by about 30 per cent. Inclined planets also eject debris at a slower rate; 80 per cent of particles were removed after 50 Myr in simulation 21, whilst simulation 56 (identical parameters but with a planet inclined by ) ejected only 50 per cent of disc bodies in the same time.
3.3 Outcome Ic  ejection of the disc
Some simulations result in the ejection of almost all disc particles. This outcome is favoured if the planet comes close to the outer disc edge, and is more likely for massive planets. All stable particles in outcome Ic simulations with low inclination planets lie in resonances, and occupy similar trajectories on the , and , planes to the resonant particles in outcome Ia. The particles do not move into such configurations, but rather remain there from the beginning of the simulations. Again few particles are caught in resonances in the inclined case, and most of those still bound at the end of outcome Ic simulations with inclined planets are in the process of being removed from the system. The planets evolve in the same manner as on Figure 3.
4 Dynamical effects
There will be three processes causing the observed behaviour: the secular (longterm) effect of the planet on the disc, the scattering of debris by the planet, and resonances between the planet and disc particles. We provide an overview of these mechanisms here, and will use them to explain our results in Section 5.
4.1 Secular interactions
Secular effects cause bodies to exchange angular momentum, whilst their energies (and hence semimajor axes) remain unchanged. This mechanism may be modelled using two equivalent methods. The most common technique is to isolate secular terms in the disturbing function by ignoring those associated with the location of the bodies on their orbits. The function is then expanded in eccentricity and inclination (often to to second order, known as LaplaceLagrange theory: see Murray & Dermott 1999) or the ratio of semimajor axes. This approach is only really suited to hierarchical or low eccentricity systems where the expansion rapidly converges; for our scenario the eccentricities and semimajor axis ratios are of order unity, and hence ignoring higher order terms will be invalid. We may however use second order theory to predict behaviour to order of magnitude accuracy, and to identify how secular effects scale with different parameters.
The second approach, known as Gauss averaging, is perhaps a more intuitive description of secular perturbations. The effect of discarding terms in the disturbing function which depend on mean longitude is equivalent to extending the mass of each body around its orbit, with the density at each point weighted inversely by the body’s velocity there. The system may then be thought of as several rings, and the gravitational interaction of these rings causes them to change shape and orientation. This evolution may be calculated numerically using the method of Touma, Tremaine & Kazandjian (2009), whereby the force of all other rings on a single point is found analytically and then these points are summed numerically over the ring. We implemented this method, using softening lengths of 1 per cent of the body’s semimajor axis. This technique has the advantage that it may be applied to systems of arbitrary eccentricities and semimajor axes, however the numerical summation makes it difficult to predict outcomes without running the code. We therefore use LaplaceLagrange theory to isolate and predict secular effects where possible, but switch to numeric Gauss averaging in cases where higher order terms significantly affect the results.
Second order theory predicts how a massless particle undergoing secular perturbations from a single planet moves on the , plane. Its instantaneous location on this plane may be thought of as the sum of two vectors, of magnitudes and , at angles and respectively to the axis. Each particle moves at a uniform rate anticlockwise around a circle of radius , the centre of which is offset by in the direction . To second order in , the forced eccentricity decreases with semimajor axis such that
(1) 
(Murray & Dermott 1999; Mustill & Wyatt 2009). The value of is set by the initial conditions, so if a particle starts with then . Similar behaviour to the above is predicted for an object’s inclination and longitude of ascending node, such that the particle moves in circles on the , plane at the same rate as before, albeit in a clockwise direction.
We now compare these predictions to our simulations. The semimajor axes of coherent disc particles in outcome Ia interactions do not change significantly (typically by less than 15 per cent over the course of a simulation), and the same is true of the bell structure particles in outcome Ib. On the , plane, coherent disc particles move anticlockwise about circles so far as we can discern by eye. These circles pass close to the origin and their centres are offset by an amount scaling with . These results agree with the predictions of low order secular theory. However the behaviour of coherent disc particles on the , plane differs from these predictions; they move in ovals on this plane aligned with , and hence rotate as evolves (see Section 5.2 for a description of the evolution of particle orbital planes measured relative to the orbital plane of the planet). We confirm this behaviour using the Gauss averaging method, and hence conclude that the coherent disc is a secular phenomenon. Bellstructure particles in outcome Ib also move in closed loops on the , and , planes, in the manner predicted by the Gauss averaging method. Hence these structures are secular phenomena too. However the shapes of these loops are much more complicated than in outcome Ia, and differ significantly from those predicted by second order theory.
We see that for outcome Ia, most predictions of second order theory appear qualitatively correct even when all bodies are significantly eccentric. On Figure 6 we plot the forced eccentricity of the outermost coherent disc particles from our outcome Ia simulations, versus the value predicted by Equation 1 (the outermost particles are those least susceptible to nonsecular effects). We see that the second order theory prediction is reasonably accurate even up to high forcing eccentricities; these values are within 25 per cent of the simulation values in all but three of our simulations (note that the maximum eccentricity of a particle is if the disc’s initial eccentricity spread is small). Thus we conclude that in the high eccentricity regime, the quantitative predictions of low order secular theory are still roughly applicable to the , and evolution provided mutual inclinations are not too high. However its predictions for the and evolution quickly break down, so Gauss averaging must be used for these parameters. If the mutual inclination is high enough to enter the Ib regime, the quantitative predictions of second order theory fail with regard to the , , and evolution. This will be due to higher order mixed terms in the disturbing function becoming nonnegligible, and significantly altering the eccentricity evolution. Therefore Gauss averaging must always be used in the Ib regime.
4.2 Planetdebris scattering
Scattering exchanges both energy and angular momentum, and hence may eject bodies from the system (unlike secular interactions). Any change in a nonresonant body’s semimajor axis must be a result of this effect, and we use this to find the location of the innermost stable nonresonant particles in Section 5.1.2. We also use scattering arguments to constrain the planet’s migration in Section 5.3.
4.3 Resonances
A particle may be longterm stable if its orbital period is a simple fraction of the planet’s, hence some particles may survive strong scattering if located in resonance. The versus plot on Figure 2 shows a significant number of particles trapped in the strong 2:1 external resonance, and also some in the weaker 3:1 resonance. Their behaviour on the , plane is similar to that of the secular particles as seen on Figure 2; those in the 2:1 resonance (green) move in trajectories which resemble smaller versions of those of coherent disc particles. The 3:1 resonance particles (orange) are also visible on this plane, forming the shape aligned with the planet and located halfway between the 2:1 resonance particles and the coherent disc. Thus resonant particles form structures similar to the coherent disc in outcome Ia, which are superimposed on the latter. We observe this on Figure 1; particles in 2:1 resonance appear to form a low eccentricity ring passing over the planet’s apocentre.
However in generating this figure we assumed all particles had random mean anomalies, which is incorrect for resonant particles as their mean anomalies are not random but depend on the instantaneous position of the planet. On Figure 7 we replot Figure 1, showing the instantaneous positions of particles in 2:1 resonance to preserve resonant structure. These particles were selected if their semimajor axes remained within two per cent of the nominal 2:1 resonance location for the final 0.5 Myr of the simulation. They do not form a continual ring but rather two crescents apart from one another. Resonant particles initially located in the two gaps would pass close to the planet when the latter was at apocentre, and have hence been ejected from resonance. The other particles in the 2:1 resonance never come close to the planet so are longterm stable. In addition the resonant angles of these particles librate, so the overdensities associated with material on resonant orbits oscillate between both ends of the crescent. This motion is slowest at the two extrema, hence these regions are overdense relative to the rest of the crescent^{2}^{2}2Note that the resonant structure on Figure 7 shows some similarities with the clumpy structure observed in the Eridani debris disc (Greaves et al., 2005). However, more detailed characterisation of this structure is needed to determine whether there is any link. (e.g. Wyatt 2006). Thus these particles do not actually have random mean anomalies, however we are probably justified in plotting them as a continuous ring on Figure 1 as any gaps would be difficult to observe. We have not plotted 3:1 resonance particles separately as there are fewer of these bodies and they lie farther from the star, so they do not contribute significantly to the surface density of the disc.
Owing to the small number of particles in our simulations and the thin resonant widths, resolution of all resonant structures is difficult and probably incomplete. On Figure 2 seven per cent of the surviving disc particles inhabit the 2:1 resonance with a further two per cent in the 3:1, amounting to under two per cent of the initial disc mass in resonant particles. However resonant features were well resolved by Faramaz et al. (2014), who ran outcome Ia simulations but with larger numbers of particles, and we refer the reader to that paper for more detailed views of the resonant structures. We expect the coherent disc will be the dominant structure in outcome Ia systems, with less massive resonant features superimposed on top of it. Likewise particles which would otherwise be ejected may be trapped in resonance in the outcome Ib interactions, however our resolution is too low to pick out such structures on Figure 4. In Ic systems the only surviving material would be in resonance with the planet, and so a sufficiently massive initial disc could result in a potentially observable amount of material on such orbits.
5 Analysis and predictions
We now characterise our results in greater detail to produce quantitative predictions for this interaction. We first examine the debris behaviour when the planet is coplanar with the disc midplane, followed by the inclined case. We then show that the evolution of the planet will almost always be small if it is at least ten times as massive as the disc, and finally we list the important timescales involved in the interaction.
5.1 Coplanar planets
For coplanar planets, surviving particles are either caught in resonances or form an elliptical disc apsidally aligned with the planet. By considering the dynamical effects from Section 4 we predict the location of the disc edges, the width of the unstable zone around planet apocentre (and hence differentiate between coplanar interactions resulting in outcome Ia and Ic) and characterise the nonuniform surface density of this structure.
5.1.1 Disc edges
By considering the behaviour of debris particles, we may write down the locations of the coherent disc edges at late times. Assuming the initial eccentricity spread is small each particle’s orbit evolves between a near circle and an ellipse, the latter aligned with the planet’s orbit, and back again. Superimposing all intermediate orbits forms an eccentric disc, as shown on Figure 8. The edges of this disc are exactly elliptical, and we henceforth label the apocentre and pericentre of the inner and outer edges as , , and as shown on the figure. A disc formed from particles at a single semimajor axis will have edges
(2) 
(3) 
and
(4) 
Our discs are composed of particles at various semimajor axes. If the maximum eccentricity decreases sufficiently quickly as semimajor axis increases then the inner and outer disc edges will be determined by particles with the smallest and largest semimajor axes respectively (this is satisfied by ). Extending the above equations, a coplanar coherent disc formed of particles with semimajor axes between and will have
(5) 
(6) 
(7) 
and
(8) 
at late times. The outermost semimajor axis will be given by , and we derive in Section 5.1.2. If the initial eccentricity spread is nonzero then no longer equals ; the above equations are roughly modified such that becomes , and . Hence these are minor corrections for realistic . Note that the coherent disc and planet maintain their alignment, meaning the shape of the inner disc edge is constant. Were this not the case and the precession rates differed, the planet would eventually carve out a much larger region as it precessed.
5.1.2 Innermost stable particle
We now calculate , which determines the inner edge of the coherent disc, using scattering arguments. All nonresonant particles with initial semimajor axes smaller than are ejected, and each particle’s eccentricity is minimised when farthest from alignment with the planet (Figure 2). Therefore if the maximum eccentricity of the innermost particle is not much greater than that of the planet, these objects will be closest when the latter is at apocentre and the former is on a near circular orbit. Furthermore the planet moves slowest relative to disc particles when it is at apocentre, which increases scattering efficiency. Therefore the inner edge of the disc must be determined by scattering when the planet is at apocentre. We define a width such that (Figure 8), where is determined by the strength of scattering at apocentre, and we see from Equation 5 that the innermost stable semimajor axis also equals . Any nonresonant particles with smaller semimajor axes are unstable and will eventually be ejected from the system. Hence predicting would allow us to determine the innermost stable semimajor axis. We could also differentiate between outcomes Ia and Ic: if the innermost stable semimajor axis lies beyond the initial outer edge of the disc , all nonresonant particles are unstable and will be ejected.
Instability zones around planetary orbits is a wellstudied problem (e.g. Wisdom 1980; Mustill & Wyatt 2012), but most results are only valid for low planet eccentricities. We therefore establish our own criterion, based on Hill sphere arguments. We find the instantaneous Hill radius of an eccentric planet at apocentre is
(9) 
(see Appendix B for the derivation). The Hill radius at pericentre is obtained by substituting for in the above, and it is larger at apocentre than at pericentre for all eccentricities and mass ratios. We hypothesise that the width of the instability zone is , where accounts for initially stable particles which evolve until they come close enough to the planet to be ejected. Estimating requires to be determined from the simulations, which is difficult to measure by eye as apparent from Figure 1. Resonant particles and those being ejected are still present, and these must be ignored if we wish to measure the disc inner edge. There are two ways to do this: firstly we may simply state that we are not interested in particles which are either planet crossing, in strong resonance or have significant inclinations. Ignoring these leaves only the coherent disc, and may be ascertained from the separation of the disc inner edge and the planet’s apocentre. The second method is to examine the small bodies in the , plane; coherent disc particles inhabit a specific region of semimajor axis space, with a well defined maximum eccentricity and unchanged (Figure 2). Therefore omitting particles whose semimajor axes have changed by more than a few percent leaves only those in resonance and the coherent disc, hence we may read off the minimum coherent disc semimajor axis . We use both of these methods to calculate and the other edges, and both give the same answers.
We plot the values of versus for coplanar outcome Ia simulations on Figure 9, as well as lower limits from coplanar outcome Ic simulations. For the latter we take the minimum value of to be the distance between the planet’s apocentre and the initial outer edge of the disc. We find that provides a good fit to the data, hence
(10) 
This result is in line with several previous works on instability zones in the low eccentricity regime, which suggest that systems are stable for 10100 Myr (the timescale of our simulations) if they are separated by of order 10 Hill radii (e.g. Chambers, Wetherill & Boss 1996, Smith & Lissauer 2009). Note that not all particles in the unstable regions may have been ejected by the end of the simulations, so the points on Figure 9 are more likely to be underestimations of than overestimations.
The cleared region will be wider at pericentre than apocentre if . Were this not the case, the innermost stable semimajor axis could potentially be determined by scattering at planet pericentre rather than apocentre. However the planet’s Hill radius at pericentre is very small, so the disc would have to be very eccentric for this to occur. Therefore the inner edge of the disc will probably be determined by scattering at planet apocentre in all cases.
5.1.3 Disc overdensities
We now characterise the coherent disc overdensities. To circumvent our low simulation resolution, and increase the number of coherent discs examined, we use the following method to generate these structures without running nbody simulations. Firstly we create many particles with semimajor axes between and , following surface density profiles as before. We then generate the preinteraction disc by assigning each particle an eccentricity between and , and a randomised orientation. For a given forcing eccentricity, a particle’s value is the distance between its initial position on the , plane and that of the forced eccentricity vector. Hence we find the circles on this plane about which each particle moves. Randomising particle positions around these circles assigns each an instantaneous eccentricity and orientation, and we set the inclinations to zero. Finally we randomise the mean anomaly of each particle to find its instantaneous position, and superimpose these to generate a coherent disc. We plot an example disc on Figure 10, which contains particles with similar parameters to those in simulation 21 (Figures 1 and 2). This shows how the disc would look at high resolution, without resonant or unstable particles.
The apocentre side of the coherent disc is much denser than the pericentre side, because eccentric bodies spend more time near apocentre than pericentre. We characterise this overdensity by first considering discs composed of particles at a single semimajor axis . These discs are functions of , and , and we generated many discs using a broad range of these parameters. We then made two wedgeshaped cuts, radiating from the star and apart, through each disc. Finally we summed the particles in each wedge and divided the two, yielding the ratio of the total crosssectional area on opposing sides of the disc. The results are shown on Figure 11 and are independent of and , only depending on forcing eccentricity. Extending this to discs composed of particles at a range of semimajor axes (which have five parameters: inner and outermost semimajor axis, , , and such that from Equation 1) we find that the pericentreapocentre surface density ratio is still roughly equal to that on Figure 11, if the forced eccentricity is now taken to be that of a particle at the mean semimajor axis. Measuring this ratio (noting that observations must be corrected for the different temperatures on both sides of the disc) would hence allow the properties of an unseen perturber to be constrained.
The coherent discs also exhibit finer structure. Unlike the general apocentre overdensity, these structures may not always be observable depending on the instrument beam size. However we will briefly describe these features, as they could potentially be observed with a high resolution instrument such as ALMA (e.g. Moór et al. 2013; Dent et al. 2014) or the HST (e.g. Golimowski et al. 2011). We use the discs generated above to calculate surface density profiles over the pericentre and apocentre wedges, with an example shown on Figure 12. Again, we first consider a disc formed of particles at one semimajor axis. There are either one or two overdensities on both sides of the disc, and we label their radial locations as , , and as on Figure 10. These peaks arise from the overlap of many orbits at these locations, where a particle’s rate of eccentricity change goes from positive to negative or vice versa. If then these peaks are extremely sharp and are located at the disc edges, and if each peak is formed from many independent peaks at slightly different locations. Hence consideration of a peak’s inner and outer edges shows its width to be , and the same for all peaks. The two apocentre peaks will merge when the outer edge of the peak overlaps the inner edge of the peak, which occurs if . The same result also arises on the pericentre side, so the number of peaks on both sides of the disc are equal. The ratios of peak surface densities are independent of and particle semimajor axis, but become more extreme as increases. Finally, whilst the ratios of these overdensities change their order will not; the peak with the highest surface density is always the inner peak at apocentre (), followed by , and finally , the inner peak at pericentre, which has the lowest surface density.
Applying this analysis to discs formed of particles at multiple semimajor axes, we find each overdensity ring has a width
(11) 
Also the two overdense rings will merge if the criterion
(12) 
is satisfied, where was defined in Equation 1. Thus discs with a broad semimajor axis range and/or low eccentricity have a single overdensity ring, whilst a narrow semimajor axis range and/or a high eccentricity yields two rings. Hence observing the number and widths of overdensity rings in an eccentric disc would yield constraints on the perturber, in addition to those from the size and shape of the coherent disc and cleared region.
5.2 Inclined planets
Our simulations show the secular evolution of initially prograde debris follow two distinct modes, depending on the planet’s initial inclination and orientation. Particle orbital planes either precess about that of the planet (outcome Ia, Figure 2), or librate about a plane orthogonal to this (Ib, Figure 5). These two behaviours have been noted in the literature (Verrier & Evans, 2009; Farago & Laskar, 2010; Doolin & Blundell, 2011; Kennedy et al., 2012). Figure 13 shows the secular evolution of test particle planes from an nbody simulation, where the reference (primed) frame is now the planet’s orbital plane with measured from the planet’s pericentre direction. For this simulation , and . Note that this semimajor axis ratio is larger than any from our full simulations, to isolate secular effects from those due to scattering or strong resonances. In spite of this difference the particle behaviours observed in our simulations are very similar to those on this plot, although the former show greater scatter due to the aforementioned effects.
We also plot the initial disc midplanes from our inclined simulations on Figure 13. Because earlier we chose our reference frame to coincide with the disc’s initial midplane with pointing towards the planet’s initial ascending node, the orbital elements of the disc midplane relative to the planet’s orbit (primed frame) are and or for prograde and retrograde planets respectively. The conversion between these frames results in an offset and a rotation in the unprimed frame, so a particle moving along a path from Figure 13 moves about a similar shape on the , plane but offset from the origin and orientated with the planet’s line of apsides. This is exactly the behaviour observed on Figures 2 and 5, and plotting these particles on the , plane shows that they move along similar paths to those on Figure 13.
Figure 13 shows that the separatrix between Ia and Ib behaviours is roughly oval shaped and orientated with its major axis along , so the critical planet inclination separating these behaviours is minimal for and maximal for or . In fact figure 3 from Doolin & Blundell (2011) shows that the critical inclination is when or , so an initially prograde particle always lies in the Ia regime if or (i.e. the planet’s line of apsides lies in the initial disc midplane). Alternatively for () initially prograde particles lie in the Ib regime if the planet is sufficiently inclined; Farago & Laskar (2010) derive this critical inclination in the limit as
(13) 
Doolin & Blundell (2011) find the shape of test particle paths on the , plane to be independent of semimajor axis and binary mass ratio (tested for and ), although not all particles are longterm stable. They also find Equation 13 is valid for all , so this equation should still hold for the regime studied in this paper. Our inclined simulations use , so the above yields a critical inclination of for . The critical planet inclination for all other values of is between that of Equation 13 and . All of these points agree with our results, hence we may generally predict whether outcome Ia or Ib will occur if , or (although intermediate values of require Figure 13 to be generated using planet specific parameters). Our only simulation which did not fall neatly into either outcome Ia, b or c is number 38, which displayed both Ia and Iblike behaviour. This is due to the small initial spread of particle inclinations and longitude of ascending nodes, which caused debris to split between behaviour modes as the initial midplane lay on the boundary between three behaviour regimes.
5.2.1 Planets below the critical inclination
Interactions result in outcomes Ia or Ic if the planet has a small enough initial inclination (as a function of and ) to place the system interior to the separatrix. In this case the predictions of Section 5.1 generally hold. Most of the surviving nonresonant particles still come closest to the planet at the latter’s apocentre, so scattering at this point still sets the innermost stable semimajor axis. This value is smaller than in the coplanar case, as particles must lie closer to the star to pass within the same distance of an inclined planet and the reduced interaction time for mutually inclined bodies lowers the scattering efficiency. Hence the locations of the coherent disc edges are similar to those in the coplanar case, but the innermost semimajor axis predicted by Equation 10 is now an upper limit.
A coherent disc sculpted by a low inclination planet looks very similar to the coplanar case, only now shaped overdensities are observed when viewing the disc side or fronton (Figure 4). These arise from the secular behaviour of the particle elements (Figure 2); particles move in ovals on the , plane, centred on the planet and aligned with . Thus a particle is maximally inclined to the planet when and minimally inclined when . Particle orbits oscillate between these two extrema, spending more time on these orbits than on intermediate ones. Hence the superposition of particles forms an when viewed along the planet’s line of apsides (fronton). This behaviour is similar to that of dust bands either side of the ecliptic in the Solar System (Low et al., 1984; Grogan et al., 2001). The height of the coherent disc is also greatest around planet apocentre, as a particle’s eccentricity is still largest when aligned with the planet; this behaviour gives rise to overdensities when viewed from the side. Thus the overdensity analysis from Section 5.1.3 still applies, although some additional vertical structure will also be present.
5.2.2 Planets above the critical inclination
If the prograde planet’s initial inclination is at least that given by Equation 13 (depending on ), the system may lie on the other side of the separatrix and will thus result in outcome Ib or Ic. In this case the orbital planes of surviving particles librate about a point orthogonal to that of the planet, resulting in a hollow bellshaped debris structure with holes at both ends. This outcome clearly differs from the coplanar case, however we still find several areas in which the previous results apply. Firstly the evolution of this debris remains secular in nature, evolving as predicted by the Gauss averaging code whilst particle energies stay constant. Debris still inhabits a similar region in , space as in the coplanar regime, and the innermost particles are still ejected. However the divide between stable and unstable semimajor axes is now less well defined. An unstable zone still exists around the planet’s orbit; a hollow inner bubble is clear of debris, and this region is roughly axisymmetric about the planet’s line of apsides. The distance between the bubble edge and the planet’s orbit (excluding the holes at both ends) is constant in the latter’s orbital plane, and approximately five apocentre Hill radii wide  similar to the coplanar case. The distance from the star to the bell structure’s outer edge is still of order the initial disc outer radius . Finally, whilst the behaviour of particles on the , plane differs considerably from the coplanar case, they still move in closed loops orientated with respect to the planet. Hence the parameters and are still strongly coupled.
This coupling of and gives rise to the holes at both ends of the bell; Figure 13 shows that no particles have exactly equal to 0 or , hence none have nodes along the planet’s line of apsides. This results in the holes, so we note that these are caused by secular rather than scattering effects. The figure also shows that the size of particle paths on the , plane (and hence the range of and over which particles evolve) is largest when the system only just lies inside the Ib regime. This means that the bell structure is most extended for systems where the mutual inclination between the planet’s orbit and the disc midplane is . For systems farther inside the Ib region (towards and ) the kidneyshaped particle paths shrink, and hence the bell becomes squashed (i.e. the holes at each end become larger) owing to the increasingly narrow and ranges over which particles oscillate. Finally in the extreme case of and all particles have nodes almost orthogonal to the planet’s pericentre, and the bell is squashed to the extent that it becomes a disc orthogonal to the planet’s orbital plane (see Kennedy et al. 2012). Simulations 57 and 58 ( and respectively, ) have such an outcome.
Our simulations highlight the importance of a planet’s argument of pericentre (i.e. the height of its pericentre above the disc midplane) in addition to its inclination in setting the secular behaviour of debris. This angle also affects the scattering of disc particles. If , the planet spends a lot of time near the disc midplane regardless of its inclination. Most importantly its apocentre (where scattering is strongest) lies in this plane, and so the innermost particles can be ejected within a few orbital periods. Alternatively if the planet’s apocentre may lie far out of the plane, and the innermost particles may survive for a few secular periods before coming close to this location. Thus the scattering efficiency is much greater for or than it is for (i.e. scattering is more important for or than for ).
5.3 Evolution of the eccentric planet
Having examined the final state of the debris particles, we now consider the orbital evolution of the eccentric planet. We noted that this evolution was small in our simulations, because we only examined planets at least ten times more massive than the disc. We now make a simple analytical prediction for the planet’s migration distance, to show that it will be small under almost all circumstances with such mass ratios. This will justify our choice of mass ratios to study in this paper, as any larger values could potentially result in significant planetary evolution.
We use scattering arguments to constrain the planet’s migration distance. As described in Section 4.2, scattering is the only mechanism which may significantly change the planet’s semimajor axis. The semimajor axes of surviving particles in the coplanar and highly inclined simulations are unchanged, so these objects are not being scattered. Thus if all unstable particles are eventually ejected, we may calculate the maximum change in the planet’s semimajor axis using simple energy arguments. If a small mass orbits a star of mass , the energy of the system is . Ignoring the additional small amount of energy associated with the potential of the planet on the small bodies, we find the total energy of many disc particles by summing the two body for each object. Hence for a dynamically cold disc with inner and outer radii and respectively, the total energy in disc particles between radial distances of and (where both radii lie within the disc) will be given by
(14) 
if or 2.
If all particles between and are just ejected by the planet, their energy goes to zero and the planet becomes more tightly bound. Therefore equating total energies before and after scattering, we find that the planet’s semimajor axis will have changed by
(15) 
where
(16) 
and is the planet’s initial semimajor axis. Hence if is of order of unity or lower, significant planet migration is impossible if the planet is considerably more massive than the disc. Substituting and equal to the smaller of and , we can predict the migration distances for the coplanar planets. This will be an upper limit; not all unstable particles may have been cleared by the end of the simulations, and some will remain trapped in resonances. This will also serve as an upper limit for the inclined case, where fewer bodies will have been ejected. We plot the planet migration distances against these predictions on Figure 14.
By taking appropriate limits of , we determine the maximum magnitude of the migration in multiple scenarios. If the disc is narrow then , so a narrow debris ring located at planet apocentre may never lead to significant migration if the latter is more massive than the disc. For a very broad disc (), or if the planet has apocentre beyond the outer disc edge, . Therefore may be larger than unity if the inner disc edge is smaller than the planet’s semimajor axis, which requires the planet to be highly eccentric. However for to be orders of magnitude greater than requires a fairly contrived set of circumstances, and we typically expect the two parameters to be of the same order. We conclude that will be of order unity or smaller in the vast majority of plausible system configurations, and hence the planet’s evolution will almost always be negligible if it is at least ten times the mass of the disc.
5.4 Timescales
Finally we will address the important timescales in the interaction. There are several timescales of interest, and we first estimate the time taken for the coherent disc or bell structure to settle into shape. We will explain this timescale in the coplanar case for simplicity, however the conclusion is valid for the inclined planet case too. The final coherent disc shape results from the particles forming well populated circles in the , plane. When the interaction starts the particles begin moving around the circles, at rates which depend on their semimajor axes. After some time the phase of and for each particle will be roughly random, and the disc will assume its final shape. At early times however, particles at similar semimajor axes have similar phases because their secular periods only slightly differ. This means that the circles in the , plane will not be randomly populated in phase (Beust et al., 2014), and the disc will not yet have its final shape; spiral density structures may be present (Wyatt 2005, Faramaz et al. 2014), and the inner and outer edges may not be apsidally aligned with the planet. This behaviour is visible on Figure 15, where we show the time evolution of the coherent disc from Figure 1. The disc will settle into its final state when each particle has undergone a sufficient number of secular oscillations such that it is no longer in phase with particles at similar semimajor axes, so the timescale for this to occur will be some number (between a few and of order 10) of secular periods of the outermost disc particle. An effect of this is that the inner edge of the disc aligns with the planet quickly, whilst the outer edge takes longer to reach a stable shape. This explains why the pericentre of the outer disc occasionally differs from in our simulations; we note that where this is the case the simulation has ended before more than a few secular timescales of the outermost particle have elapsed.
Secondly there is the time taken for the planet to clear particles from the unstable regions. Figure 9 shows that particles are unstable if they come within several Hill radii of the planet’s apocentre. However the secular evolution of a nonplanet crossing particle periodically moves its orbit away from the planet’s apocentre (see Figure 8), effectively shielding it from ejection for long periods. After a secular period the particle comes close to planet apocentre again and may be scattered, possibly repeatedly, before once again moving away through secular evolution. After one or more such interactions the particle will be ejected. Thus we expect the time taken to clear the unstable regions to be set by the longer of the secular and diffusion timescales.
Second order theory predicts a secular timescale for a test mass with semimajor axis (if ) of
(17) 
where is the planet’s orbital period, and is a Laplace coefficient (Murray & Dermott, 1999). For the outermost unstable particle ; if all nonresonant particles have been ejected in the simulation, this is taken to be . The diffusion timescale for this particle will be of order
(18) 
(Tremaine, 1993). On Figure 16 we plot the ratios of (the time taken to clear 95 per cent of nonresonant particles from the unstable regions) to each of these timescales, versus the ratio of the two timescales. We see that the clearing time is roughly equal to ten times the longer of the two timescales for all of our simulations. It should be noted that the majority of unstable particles are removed on timescales much shorter than this, and so the rapid decline in the planet’s eccentricity at the beginning of the simulations may be attributed to the high frequency of scattering events.
Note the debris clump in the first panel of Figure 15; even though particles associated with this structure are planet crossing, their early behaviour is primarily secular in nature. Secular perturbations rapidly pump the eccentricities of these bodies towards unity, whilst their pericentres remain aligned with each other (but not with the planet). The superposition of these particles results in an elliptical overdensity misaligned with the planet’s orbit. These bodies then diffuse in phase space owing to their differing secular periods as previously described, and the clump becomes symmetrical about the planet’s line of apsides. This secular behaviour of planetcrossing bodies is analysed by Beust et al. (2014), and we refer the reader to their figures 6 and 7 for the physical structures formed by such particles and their corresponding phase diagrams. Eventually all nonresonant planetcrossing debris is ejected by scattering, and the inner region of the disc is cleared (final panel of Figure 15).
The third important timescale is that for system precession (see Figure 3). The coherent disc and bell structures exert secular perturbations upon the planet, causing the latter to undergo pericentre precession. However these structures are aligned with the planet. Therefore the debris causes the planet to precess, which in turn causes the debris to precess, so the whole system rotates as one. The precession period is therefore the secular period of the planet. Strictly speaking, this should be calculated using the Gauss averaging technique by summing up the perturbations from the large number of particles at different semimajor axes, eccentricities and orientations. However we find that this may be estimated to order of magnitude accuracy by modelling the debris as a single particle of semimajor axis , where , and calculating the precession rate of the planet using second order theory. For a planet perturbed by an external disc (i.e. ), the precession timescale is
(19) 
where is now defined as and is the fraction of disc particles remaining in the system. The final precession timescale predicted for simulation 21 using Equation 19 is years. The precession rate at the end of simulation 21 is / Myr (Figure 3), which implies an actual precession timescale of years. Thus the above prediction shows reasonable agreement with simulation.
6 Discussion
We have investigated the dynamical interaction between an eccentric planet and a debris disc. In doing so we have gone further than previous studies, having examined a very broad region of parameter space and produced several quantitative results. We now discuss the caveats and potential applications of our work.
6.1 Caveats
Firstly, we did not include any gas dynamics in our simulations. We consider this to be appropriate because the interaction would likely occur after gas has dissipated from the protoplanetary disc, as gas drag would circularise forming bodies. Furthermore the removal of gas could destabilise objects, which could lead to interactions driving up a planet’s eccentricity and thus act as the trigger for this scenario.
Secondly, we only simulated a single eccentric planet. In omitting additional perturbers, we narrowed the parameter space considerably and isolated the effect of the eccentric planet on the disc and vice versa. This also meant that we made no assumptions about how the planet became eccentric. However several scenarios would require additional bodies in the system. The most obvious mechanism to rapidly increase a planet’s eccentricity is scattering by a massive second planet, most likely located at the former’s pericentre. These planets would have to decouple from each other in order to prevent further scattering (also noted by Faramaz et al. 2014), and we see that this is possible from our simulations; the eccentricity of the scattered planet decreases rapidly as it first ploughs into the disc, and an eccentricity change of per cent is possible over tens of orbital timescales (Figure 3). The much smaller semimajor axis change means the eccentric planet’s pericentre rapidly moves outwards (by 2 au in the simulation on Figure 3). This could bring it far out of the scattering object’s Hill radius and which could be sufficient to decouple the two bodies. One could envisage other decoupling mechanisms; for example if the scattering body moved inwards as a result of the interaction, it could then interact with interior bodies and evolve over a much shorter timescale than the more distant eccentric planet. In fact there need not be additional objects at all; if planetary mergers occur, the result may be a single body on a highly eccentric orbit (e.g. Lin & Ida 1997). Thus the results of our investigation should be broadly correct for many cases even if additional planets are present.
However, if there were additional perturbers in the system which remain coupled with the eccentric planet, these could repeatedly scatter that planet. The latter’s orbit would therefore change with each interaction. This would probably prevent the coherent disc or bell structure from forming, as particles would rapidly switch between stability and instability. Many more particles (and possibly the eccentric planet) would be ejected in this regime. Longterm stable particles would also lie farther from the planets, so the resulting debris structure would probably not be as extreme as in the single planet case.
Even in the absence of further planetplanet scattering, secular and resonant interactions in a multiplanet system could result in the outermost eccentric planet (which is presumably that considered here) continually evolving, and hence the system may not have a well defined final state as in our simulations. Even without planet evolution, the debris behaviour could still differ significantly in this regime. For example, secular perturbations would not necessarily cause particles to be apsidally aligned with the outermost planet, since the longitude of forced pericentre (which determines the orientations of the centres of the circular particle tracks on the , plane) would depend on the secular solution for all of the planetary bodies. In the coplanar case this would cause the degree of alignment between the eccentric planet and coherent disc to vary, which would result in a larger innermost stable semimajor axis. Three (or more) body resonances could also trap particles in different locations to the two body resonances in the single planet case, and the reduced secular times could lead to more rapid clearing and disc settling if more than one planet were present.
Finally our simulations only model the largest debris particles, for which gravity is the significant force. Smaller dust grains are also affected by radiation pressure, which is not included here. Optical images only show these small grains, hence the observed debris structures may differ in appearance to those seen in our simulations. As such our results are more comparable to longer wavelength (i.e. ALMA) images of larger particles, and act as tracers of the parent body dynamics.
6.2 Applications
We showed that material surviving an encounter with an eccentric planet forms structures which would not otherwise be expected, and our work could potentially predict the properties and dynamical history of such unseen planets from observed debris. Highly eccentric, coplanar planets force surviving material into a coherently eccentric disc, which may be significantly elliptical. Observing such a structure could be indicative of an unseen eccentric planet in roughly the same plane as the disc and aligned with it. A sharp inner disc edge, possibly combined with resonant structure, would additionally suggest that the planet is strongly scattering disc material; this could be evidence of previous planetary evolution, which would bring the planet close enough to the disc to remove nonresonant material and maintain a sharp disc edge. In the opposite case of a highly inclined eccentric planet, the resulting bellshaped debris structure is even more unusual. Observing this structure would not only yield the orientation of the planet’s current orbit (debris is symmetric above and below the planet’s orbital plane, and is offset from the star in the direction of planet apocentre) but could also point towards a violent dynamical history; this structure only forms if the planet is significantly inclined to the disc at early times, which could hint at some major planetary evolution in the system’s past. Our results could also identify an ongoing planetdisc interaction; a debris structure with a global offset, spiral structure and a central, stationary overdensity could suggest a planet has recently been placed onto an eccentric disc crossing orbit (Figure 15).
Beyond the simple application above, our numerical results could also place more detailed constraints on a sculpting planet if debris has settled into its final state. If the disc and planet lie in the same plane, the forcing eccentricity of a particle forming the inner disc edge may be calculated from the edge shape using Equations 5 and 6 (assuming the preinteraction disc was dynamically cold). Using either second order secular theory (Equation 1) or Gauss averaging would then constrain the planet’s eccentricity as a function of its semimajor axis, and assuming the planet is no longer disc crossing yields further constraints. The planet must be massive enough to have cleared the inner regions within the age of the system, yielding two lower mass limits as functions of semimajor axis (Equations 17 and 18). Finally the width of the instability region (Equations 9 and 10) gives another mass constraint. Thus we may bound the mass, semimajor axis, eccentricity, orientation and orbital plane of an interior perturbing planet using the shape of the coherent disc alone. Resonant structures superimposed on this disc would yield additional semimajor axis constraints, and could even point to the planet’s current location on its orbit. Finally density variations within the disc could produce additional constraints on forcing eccentricity (Figure 11 and Equations 11 and 12). Similar constraints may be found even if the planet is initially inclined with respect to the disc, as many of the numeric results still hold in this regime (see Sections 5.2.1 and 5.2.2). Note however that Gauss averaging must be used in the high inclination regime, as the quantitative predictions of second order secular theory fail here.
Large numbers of debris discs have now been imaged, with many displaying significant asymmetries and/or stellar offsets. We apply the above method to the HD 202628 and Fomalhaut systems as examples (we stress that we are not actually attempting to explain these systems via this scenario, but rather to show how our results could potentially be applied). Both systems host resolved, eccentric debris rings with sharp, elliptical inner edges. For our purposes we assume both rings to be remains of initially broader, circular discs truncated by hypothetical eccentric planets with apocentres interior to their inner edges. The lack of significant vertical structure implies the hypothetical planets are coplanar (whilst such discs could also occur if , coplanar planets seem much more likely). Starting with Fomalhaut, Kalas et al. (2013) finds that the deprojected inner disc edge may be fitted with an ellipse of semimajor axis 138 au and eccentricity 0.12. Assuming a 1.92 star of age 450 Myr (Mamajek, 2012) we plot the constraints derived using the above method on Figure 17, along with the observational upper mass limits for the hypothetical object (Kenworthy et al., 2013). We find that such an object must have a semimajor axis of au, an eccentricity of and a mass of Jupiter masses () to sculpt the observed disc inner edge. These predictions compare favourably with those of Chiang et al. (2009). However problems arise as we consider this system in more detail; we have neither explained the highly eccentric Fomalhautb nor included the additional perturbations it could exert on debris (Tamayo, 2014; Beust et al., 2014). Observed dust also forms a narrow ring whilst our model predicts a broad disc, which we also cannot explain; it is implicit in our model that the proper eccentricities are comparable with the forced eccentricities, which is not the case for Fomalhaut. Therefore we can probably rule out our model as an explanation for the Fomalhaut system, for which it would be more appropriate to study using systemspecific simulations.
Figure 17 also shows our constraints on an eccentric perturber in the HD 202628 system, derived for an elliptical inner disc edge of semimajor axis 158 au and eccentricity 0.18, with a solar mass star of age 2.3 Gyr (Krist et al., 2012). We found no upper mass limits in the literature, but if we assume the perturber is a planet rather than a brown dwarf (i.e. ) then the hypothetical object must have a semimajor axis of au, an eccentricity of and a mass . Thus our model could potentially explain this system. Whilst really more suited to larger disc eccentricities than observed in either of the above systems, these examples show that our results may be used to quickly estimate the major parameters of a hypothetical perturber using just a few observables. These could then form the basis of more detailed system specific studies.
Discs detected in the farinfrared and submillimetre (such as Reticuli: Eiroa et al. 2010) often exhibit what could be interpreted as more extreme asymmetries to those seen in Fomalhaut and HD 202628, and more akin to those predicted here. However, there remains some uncertainty as to whether the offset emission in such systems arises from circumstellar material or from background galaxies (e.g. Panić et al. 2013). Galaxies have a similar temperature to that expected for cold dust at au from a star, and are ubiquitous enough to be commonly found in close proximity to nearby stars. Thus it is often difficult to distinguish between these two interpretations (e.g. Nilsson et al. 2010). Our study shows is that it is not possible to rule out a circumstellar origin for the emission by virtue of its nonaxisymmetry about the star, since we find that highly asymmetric disks are a natural outcome in systems with planets on highly eccentric orbits. Given the large number of highly eccentric extrasolar planets detected, and the significant eccentricity excitation predicted by Solar System evolution models and some Hot Jupiter formation theories (e.g. Boley, Payne & Ford 2012), highly eccentric discs could be more common than one might naively expect. The abundances of both eccentric planets and misshapen debris discs will help determine the frequency of systems with eccentric planets encountering discs, and our work provides testable predictions for such a scenario.
7 Conclusions
We used nbody simulations and analytics to characterise the interaction between an eccentric planet and a less massive debris disc. We identify a chaotic zone of predictable width around the planet’s orbit, and any nonresonant debris here will be ejected within secular or diffusion times (whichever is longer). Surviving nonresonant debris undergoes secular evolution, forming either an eccentric disc, a hollow bellshaped structure or a polar disc, depending on whether the planet’s initial inclination with respect to the disc is low, moderate or high respectively. These structures are apsidally aligned with the planet and centred on its orbital plane, and form within a few secular times. They also contain overdensities at predictable locations, which could be used to constrain parameters of an unseen perturber. The planet’s dynamical evolution will almost always be negligible in this mass regime. Our results may be used to predict the general outcomes of this scenario, or to infer the properties of unseen eccentric planets from observed debris features.
8 Acknowledgements
We thank Dan Tamayo for his insightful review, which was a great help in improving this paper. TDP acknowledges the support of an STFC studentship, and MCW is grateful for support from the European Union through ERC grant number 279973.
References
 Augereau et al. (2001) Augereau J. C., Nelson R. P., Lagrange A. M., Papaloizou J. C. B. & Mouillet D., 2001, A&A, 370, 447
 Backman & Paresce (1993) Backman D. E. & Paresce F., 1993, in Protostars and Planets III, ed. Levy E. H. & Lunine J. I. Univ. Arizona Press, Tucson
 Beust et al. (2014) Beust H. et al., 2014, A&A, 561, A43
 Bitsch et al. (2013) Bitsch B., Crida A., Libert A.S. & Lega E., 2013, A&A, 555, A124
 Boley, Payne & Ford (2012) Boley A. C., Payne M. J. & Ford E. B., 2012, ApJ, 754, 57
 Chambers (1999) Chambers J. E., 1999, MNRAS, 304, 793
 Chambers, Wetherill & Boss (1996) Chambers J. E., Wetherill G. W. & Boss A. P., 1996, Icarus, 119, 261
 Chiang et al. (2009) Chiang E., Kite E., Kalas P., Graham J. R. & Clampin M., 2009, ApJ, 693, 734
 Dent et al. (2014) Dent W. R. F. et al., 2014, Sci, 343, 1490
 Doolin & Blundell (2011) Doolin S. & Blundell K. M., 2011, MNRAS, 418, 2656
 Eiroa et al. (2010) Eiroa C. et al., 2010, A&A, 518, L131
 Farago & Laskar (2010) Farago F. & Laskar J., 2010, MNRAS, 401, 1189
 Faramaz et al. (2014) Faramaz V. et al., 2014, A&A, 563, A72
 Ford & Rasio (2008) Ford E. B. & Rasio F. A., 2008, ApJ, 686, 621
 Golimowski et al. (2011) Golimowski D. A. et al., 2011, AJ, 142, 30
 Greaves et al. (2005) Greaves J. S. et al., 2005, ApJ, 619, L187
 Grogan et al. (2001) Grogan K., Dermott S. F., Durda D. D., 2001, Icarus, 152, 251
 Hamilton & Burns (1991) Hamilton D. P. & Burns J. A., 1991, Icarus, 92, 118
 Hamilton & Burns (1992) Hamilton D. P. & Burns J. A., 1992, Icarus, 96, 43
 Hayashi (1981) Hayashi C., 1981, Progress Theor. Phys. Suppl., 70, 35
 Ida & Lin (2004) Ida S. & Lin D. N. C., 2004, ApJ, 616, 567
 Kalas, Graham & Clampin (2005) Kalas P., Graham J. R. & Clampin M., 2005, Nat, 435, 1067
 Kalas et al. (2013) Kalas P., Graham J. R., Fitzgerald M. P. & Clampin M., 2013, ApJ, 775, 56
 Kennedy et al. (2012) Kennedy G. M. et al., 2012, MNRAS, 421, 2264
 Kenworthy et al. (2013) Kenworthy M. A. et al., 2013, ApJ, 764, 7
 Krist et al. (2012) Krist J. E., Stapelfeldt K. R., Bryden G. & Plavchan P., 2012, ApJ, 144, 45
 Lin & Ida (1997) Lin D. N. C. & Ida S., 1997, ApJ, 477, 781
 Lissauer (1993) Lissauer J. J., 1993, ARA&A, 31, 129
 Low et al. (1984) Low F. J. et al., 1984, ApJ, 278, L19
 Malmberg & Davies (2009) Malmberg D. & Davies M. B., 2009, MNRAS, 394, L26
 Mamajek (2012) Mamajek E. E., 2012, ApJL, 754, L20
 Moór et al. (2013) Moór A. et al., 2013, ApJL, 777, L25
 Murray & Dermott (1999) Murray C. D. & Dermott S. F., 1999, Solar System Dynamics. CUP, Cambridge
 Mustill & Wyatt (2009) Mustill A. J. & Wyatt M. C., 2009, MNRAS, 399, 1403
 Mustill & Wyatt (2012) Mustill A. J. & Wyatt M. C., 2012, MNRAS, 419, 3074
 Najita & Williams (2005) Najita J. & Williams J. P., 2005, ApJ, 635, 625
 Nilsson et al. (2010) Nilsson R. et al., 2010, A&A, 518, A40
 Panić et al. (2013) Panić O. et al., 2013, MNRAS, 435, 1037
 Patel, Metchev & Heinze (2014) Patel R. I., Metchev S. A. & Heinze A., 2014, ApJS, 212, 10
 Schneider et al. (2009) Schneider G., Weinberger A. J., Becklin E. E., Debes J. H. & Smith B. A., 2009, ApJ, 137, 53
 Smith & Lissauer (2009) Smith A. W. & Lissauer J. J., 2009, Icarus, 201, 381
 Stern & Colwell (1997) Stern S. A. & Colwell J. E., 1997, AJ, 114, 841
 Tamayo (2014) Tamayo D., 2014, MNRAS, 438, 3577
 Tremaine & Zakamska (2004) Tremaine S. & Zakamska N. L., 2004, in Holt S. S., Deming D., eds, AIP Conf. Proc. Vol. 713, The Search for Other Worlds. Am. Inst. Phys., New York, p. 243
 Tsiganis et al. (2005) Tsiganis K., Gomes R., Morbidelli A. & Levison H. F., 2005, Nat, 435, 459
 Touma, Tremaine & Kazandjian (2009) Touma J. R., Tremaine S. & Kazandjian M. V., 2009, MNRAS, 394, 1085
 Tremaine (1993) Tremaine S., 1993, ASPC, 36, 335
 Verrier & Evans (2009) Verrier P. E. & Evans N. W., 2009, MNRAS, 394, 1721
 Walsh et al. (2011) Walsh K. J., Morbidelli A., Raymond S. N., O’Brien D. P. & Mandell, A. M., 2011, Nat, 475, 206
 Wisdom (1980) Wisdom J., 1980, AJ, 85, 8
 Wu & Lithwick (2011) Wu Y. & Lithwick Y., 2011, ApJ, 735, 109
 Wyatt et al. (1999) Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., Holmes E. K. & Piña R. K., 1999, ApJ, 527, 918
 Wyatt (2005) Wyatt M. C., 2005, A&A, 440, 937
 Wyatt (2006) Wyatt M. C., 2006, ApJ, 639, 1153
 Wyatt (2008) Wyatt M. C., 2008, ARA&A, 46, 339
 Wyatt & Dent (2002) Wyatt M. C. & Dent W. R. F, 2002, MNRAS, 334, 589
 Wyatt, Dent & Greaves (2003) Wyatt M. C., Dent W. R. F & Greaves J. S., 2003, MNRAS, 342, 876
Appendix A Simulation parameters
Simulation 
(au)  (au)  (au)  ()  ()  Outcome  Notes  

1  25  0.8  20  60  10  1  Ia  
2  25  0.8  20  60  10  0.1  Ia  
3  25  0.8  20  60  100  10  Ia  
4  25  0.8  20  60  1000  10  Ic  
5  25  0.8  20  60  1000  10  Ic  Ended at 8 Myr 
6  25  0.8  20  120  1000  20  Ia  
7  15  0.67  20  60  1000  10  Ia  
8  15  0.67  20  120  1000  10  Ia  
9  25  0.8  20  200  1000  1  Ia  
10  25  0.8  45  120  1000  10  Ia  
11  25  0.8  20  60  1000  10  Ic  
12  25  0.8  20  60  1000  1  Ic  
13  25  0.8  40  60  1000  10  Ic  
14  30  0.83  20  60  1000  10  Ic  
15  17.5  0.71  20  60  1000  10  Ia  
16  25  0.2  20  60  1574  10  Ia  Ended at 5 Myr 
17  29.6  0.83  54.1  139.2  1000  0.1  Ia  
18  25.8  0.81  5  96.1  100  0.1  Ia  
19  14.3  0.65  8.3  170.3  100  1  Ia  
20  18  0.72  19.8  31.6  1000  100  Ic  
21  25  0.8  20  60  100  10  Ia  particles 
22  25  0.8  20  60  1000  10  Ic  particles 
22  40  0.9  65  85  100  10  Ic  Ended at 500 Myr 
23  35  0.8  55  70  100  10  Ic  Ended at 500 Myr 
24  20  0.9  35  39.5  10  1  Ia  Ended at 500 Myr 
25  42.5  0.88  20  60  1000  10  Ic  
26  6.3  0.2  7.5  25  100  1  Ia  
27  6.7  0.25  7.5  25  100  1  Ia  
28  7.1  0.3  7.5  25  100  1  Ia  
29  7.7  0.35  7.5  25  100  1  Ia  
30  8.3  0.4  7.5  25  100  1  Ia  
31  9.1  0.45  10  25  100  1  Ia  
32  10  0.5  10  25  100  1  Ia  

Simulation 
(au)  ()  ()  (au)  (au)  ()  ()  Outcome  

33  25  0.8  5  0  20  60  100  10  Ia 
34  25  0.8  10  0  20  60  100  10  Ia 
35  25  0.8  20  0  20  60  100  10  Ia 
36  25  0.8  30  0  20  60  100  10  Ia 
37  25  0.8  60  0  20  60  100  10  Ia 
38  25  0.8  90  0  20  60  100  10   
39  25  0.8  5  0  20  60  1000  10  Ic 
40  25  0.8  10  0  20  60  1000  10  Ic 
41  25  0.8  19  0  20  60  1000  10  Ic 
42  25  0.8  30  0  20  60  1000  10  Ic 
43  25  0.8  30  0  20  70  1000  10  Ia 
44  25  0.8  5  45  20  60  100  10  Ia 
45  25  0.8  10  45  20  60  100  10  Ia 
46  25  0.8  20  45  20  60  100  10  Ia 
47  25  0.8  30  45  20  60  100  10  Ib 
48  25  0.8  5  45  20  60  1000  10  Ic 
49  25  0.8  10  45  20  60  1000  10  Ic 
50  25  0.8  19  45  20  60  1000  10  Ic 
51  25  0.8  30  45  20  60  1000  10  Ib/Ic 
52  25  0.8  30  45  20  70  1000  10  Ib 
53  25  0.8  5  90  20  60  100  10  Ia 
54  25  0.8  10  90  20  60  100  10  Ia 
55  25  0.8  20  90  20  60  100  10  Ib 
56  25  0.8  30  90  20  60  100  10  Ib 
57  25  0.8  60  90  20  60  100  10  Ib 
58  25  0.8  90  90  20  60  100  10  Ib 
59  25  0.8  5  90  20  60  1000  10  Ic 
60  25  0.8  10  90  20  60  1000  10  Ic 
61  25  0.8  19  90  20  60  1000  10  Ib 
62  25  0.8  30  90  20  60  1000  10  Ib 

Appendix B Derivation of the eccentric Hill radius
For circular orbits the Hill radius may be derived by determining the distance to the first Lagrange point, where the combined gravitational force of the planet and star on a test mass, as well as any fictitious forces arising from system rotation, exactly balance. However the standard concept of the Hill radius does not work in the eccentric case, because the planet’s changing radial distance and velocity means that the width of this region varies over the orbit. Nonetheless we may still calculate the “instantaneous” Hill radius of an eccentric planet at a point on its orbit using the same arguments. A similar approach was taken by Hamilton & Burns (1991, 1992) to find the instantaneous Hill radius at pericentre, only they ommited velocity terms which affect the final result. Here we find the instantaneous Hill radii at pericentre and apocentre without omitting such terms.
Care must be taken with our definitions of the instantaneous Lagrange points for an eccentric body; in the circular case, these were derived by balancing the gravitational and rotational forces. For the eccentric case however we note that these forces should not exactly balance, because the planet is accelerating inwards at apocentre and outwards at pericentre. Therefore the radial acceleration of the test mass at pericentre and apocentre should match that of the planet if the former is at the first Lagrange point. Noting that the centrifugal force is the only nonzero rotational force at these points (see Hamilton & Burns 1991, 1992 for more details), we arrive at
(20) 
where is the gravitational constant, and are the primary (star) and secondary (planet) masses respectively, is the distance between the secondary and the test mass (the distance to the first Lagrange point), is the distance between the test mass and the primary, is the instantaneous angular velocity of the system, is the distance between the test mass and the system centre of mass, and is the instantaneous distance between the primary and secondary. We may calculate by considering the radial acceleration of the secondary, yielding
(21) 
Note that this is positive at pericentre and negative at apocentre. The system angular velocity at apocentre will be
(22) 
where and are the semimajor axis and eccentricity of the  binary. The angular velocity at pericentre may be obtained by substituting for in the above equation. Finally it is clear that and are given by and respectively.
Substituting these into Equation 20, and making the approximation that as , we arrive at the following cubic for the first Lagrange point at apocentre
(23) 
where and is small. Again, the equation for the pericentre case is obtained by substituting for in the above. The term is negligible and may be discarded, hence we find that the Hill radius at apocentre is given by
(24) 
which reduces to the circular case if . The Hill radius at pericentre is obtained by substituting for , and we show the Hill radius at pericentre and apocentre on Figure 18. If we repeat this analysis for the second Lagrange point (which lies along the primary  secondary axis but on the far side of the planet), we find that both of these points are approximately equidistant from the smaller mass, as in the circular case.