# The effect of sculpting planets on the steepness of debris-disc inner edges

Tim D. Pearce1,2, Alexander V. Krivov2, Antranik A. Sefilian2,3, Marija R.
Jankovic4, Torsten Löhne2,Tobias Morgner2, Mark C. Wyatt5, Mark Booth6,2,
Sebastian Marino7
1Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4
7AL, UK
2Astrophysikalisches Institut und Universitätssternwarte, Friedrich-Schiller-
Universität Jena, Schillergäßchen 2-3, D-07745 Jena, Germany
3Center for Advanced Mathematical Sciences, American University of Beirut, PO
Box 11-0236, Riad El-Solh, Beirut 11097 2020, Lebanon
4Institute of Physics Belgrade, University of Belgrade, Pregrevica 118, 11080
Belgrade, Serbia
5Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
6UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill,
Edinburgh EH9 3HJ, UK
7Department of Physics and Astronomy, University of Exeter, Stocker Road,
Exeter, EX4 4QL, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Debris discs are our best means to probe the outer regions of planetary
systems. Many studies assume that planets lie at the inner edges of debris
discs, akin to Neptune and the Kuiper Belt, and use the disc morphologies to
constrain those otherwise-undetectable planets. However, this produces a
degeneracy in planet mass and semimajor axis. We investigate the effect of a
sculpting planet on the radial surface-density profile at the disc inner edge,
and show that this degeneracy can be broken by considering the steepness of
the edge profile. Like previous studies, we show that a planet on a circular
orbit ejects unstable debris and excites surviving material through mean-
motion resonances. For a non-migrating, circular-orbit planet, in the case
where collisions are negligible, the steepness of the disc inner edge depends
on the planet-to-star mass ratio and the initial-disc excitation level. We
provide a simple analytic model to infer planet properties from the steepness
of ALMA-resolved disc edges. We also perform a collisional analysis, showing
that a purely planet-sculpted disc would be distinguishable from a purely
collisional disc and that, whilst collisions flatten planet-sculpted edges,
they are unlikely to fully erase a planet’s signature. Finally, we apply our
results to ALMA-resolved debris discs and show that, whilst many inner edges
are too steep to be explained by collisions alone, they are too flat to arise
through completed sculpting by non-migrating, circular-orbit planets. We
discuss implications of this for the architectures, histories and dynamics in
the outer regions of planetary systems.

###### keywords:

planet-disc interactions – planets and satellites: dynamical evolution and
stability – circumstellar matter

††pubyear: 2023††pagerange: The effect of sculpting planets on the steepness
of debris-disc inner edges–E.3

##  1 Introduction

Debris discs, like the Solar System’s Asteroid and Kuiper Belts, are
circumstellar populations of sub-planet-mass objects (Wyatt, 2008). These
objects undergo destructive collisions that release observable dust, and such
dust is detected around 20% of main-sequence stars (Hughes et al., 2018;
Wyatt, 2020). Extrasolar debris discs can be resolved in scattered light and
thermal emission (e.g. Esposito et al. 2020; Lovell et al. 2021; Booth et al.
2023), and are invaluable probes of processes in the outer regions of
planetary systems; modern planet-detection techniques are insensitive to mid-
sized planets orbiting at 10s or 100s of au from stars, but such planets can
be inferred from their influence on observed debris (e.g. Mouillet et al.
1997; Wyatt et al. 1999; Quillen 2006; Pearce & Wyatt 2014; Sefilian et al.
2021; Stuber et al. 2023). The launch of the James Webb Space Telescope (JWST)
has prompted renewed interest in such hypothetical planets, because if they
exist, then many should be detectable by this facility for the first time
(Pearce et al., 2022).

Many theoretical studies have used the shapes and locations of debris discs to
infer the properties of unseen planets. If such planets orbit just interior to
debris discs, and are responsible for sculpting debris-disc inner edges, then
the shapes and locations of these edges can be used to infer the planets’
minimum-possible mass, maximum semimajor axis and minimum eccentricity (e.g.
Pearce & Wyatt 2014; Faramaz et al. 2019; Pearce et al. 2022). However, this
approach suffers from several key unknowns. First, whilst Neptune dominates
the Kuiper Belt’s inner edge (Malhotra, 1993), it is unclear whether planets
actually do reside just interior to extrasolar debris discs. It could be that
the disc shapes and locations are instead set by other, non-planetary
processes, such as debris-debris collisions or system formation. Second, even
if exoplanets do sculpt debris-disc edges, it is unclear whether they do so on
fixed orbits (as often assumed in planet-disc-interaction studies), or whether
they migrate over time. Third, even if planetary sculpting occurs, it is
unclear whether this process has finished in observed discs (as often
assumed), or whether disc edges are still evolving under planet-debris
interactions. These distinctions have significant implications for the
morphology of debris discs, and the masses and locations of any responsible
planets today. Such questions are difficult to answer using only debris-disc
shapes and inner-edge locations.

However, we can gain more insight by also considering the radial surface-
density profile of debris-disc edges. If an edge were sculpted by a non-
migrating planet, then this profile would have a characteristic steepness that
depends on planet mass (Mustill & Wyatt, 2012; Rodigas et al., 2014; Quillen,
2006; Chiang et al., 2009). Conversely, if the edge were set by non-planetary
processes like collisions, then its steepness would be different (Wyatt et
al., 2007; Löhne et al., 2008; Geiler & Krivov, 2017; Imaz Blanco et al.,
2023). Similarly, the profiles of edges sculpted by migrating planets would be
distinct from those due to non-migrating planets (Friebe et al., 2022), and
edges set by combinations of both planetary and non-planetary effects would
also have different profiles (Thebault et al., 2012; Nesvold & Kuchner, 2015).
The profiles of debris-disc edges are therefore powerful tools for probing the
architectures, histories and processes in the outer regions of planetary
systems.

Recently, several studies used ALMA to resolve the radial profiles of debris-
disc edges (e.g. Lovell et al. 2021; Faramaz et al. 2021; Marino 2021). These
observations trace millimetre-sized grains that are unaffected by radiation
forces (Hughes et al., 2018), so offer our best means to probe the
distribution of larger, unseen planetesimals in the discs. Such studies show
that some debris discs have flatter inner edges that are consistent with pure
collisional evolution, whilst others have steeper profiles potentially
indicative of planetary sculpting (Imaz Blanco et al., 2023). Our aim in this
paper is to thoroughly investigate the dynamical effects of sculpting planets
on the steepness of debris-disc inner edges, and establish whether any of the
observed ALMA steepnesses are consistent with planetary sculpting. We will
show that, whilst many inner edges are too steep to be explained through
collisions alone, they are too flat to be purely the result of completed
sculpting by non-migrating, circular-orbit planets. This potentially implies
that different processes, or combinations of processes, are actually operating
in these discs.

Several previous studies also explored the effect of sculpting planets on the
steepness of debris-disc inner edges (e.g. Quillen 2006; Quillen & Faber 2006;
Chiang et al. 2009; Mustill & Wyatt 2012; Rodigas et al. 2014; Nesvold &
Kuchner 2015; Regály et al. 2018). Our study differs from these in several
ways. First, since we do not know debris-disc masses or the sizes of the
largest planetesimals (Krivov & Wyatt, 2021), both of which are critical for
quantifying collisional evolution, we choose to decouple the dynamical and
collisional modelling in this paper. Instead, we initially model the
interaction between a planet and a collisionless disc, to isolate the effect
of a planet on the edge steepness. We then apply a collisional model, to
demonstrate how collisions would change the steepness from the planet-only
regime. Second, we explore a very broad parameter space in both the dynamical
and collisional simulations, to provide general predictive models that are
valid across a wide range of scenarios. Third, we parameterise our results in
a form that directly relates to those fitted to ALMA data, to facilitate the
application of our results to upcoming observations. Finally, we apply our
predictive model to recent ALMA observations of debris discs.

The layout of this paper is as follows. Section 2 describes our $n$-body
simulations of a circular-orbit planet and a collisionless debris disc.
Section 3 provides a predictive model relating the steepness of a disc’s edge
to the parameters of a circular-orbit planet, and specifically Section 3.2.1
gives a step-by-step method demonstrating how to infer an unseen sculpting
planet from a debris disc’s inner edge. Section 4 addresses how collisions
would affect planet-sculpted edges. In Section 5 we apply our results to
observed debris discs, and in Section 6 we discuss our results (including what
the observed-disc profiles may be telling us about planetary systems). We
conclude in Section 7.

##  2 N-body simulations

We use a large suite of $n$-body simulations to explore the dynamical
interaction between a circular-orbit planet and a collisionless debris disc.
The effects of debris-debris collisions will be assessed later (Section 4). We
describe the circular-planet simulation setup in Section 2.1, discuss two
example simulations in Section 2.2, and present the surface-density profiles
of all simulated discs in Section 2.3. We measure the debris-excitation level
at the simulated inner edges in Section 2.4, and show this to arise through
MMR interactions with the circular-orbit planet. In addition, we also extend
the $n$-body analyses to eccentric planets in Appendix B.

###  2.1 Setup of $n$-body simulations

The initial setup of our $n$-body simulations with circular-orbit planets is
shown on Figure 1. Each simulation comprises one star, one planet, and a disc
of 20,000 massless debris particles. We run over 300 simulations, each with a
different combination of star mass, planet mass, planet semimajor axis and
initial-disc excitation level.

Figure 1: Initial setup of our $n$-body simulations with a circular-orbit
planet. A single planet (black circle) is initialised just interior to a
massless debris disc, which spans 1 to 30 Hill radii exterior to the planet
orbit.

Each debris particle is initialised with a semimajor axis between 1 and 30
Hill radii ($r_{\rm H}$) exterior to the planet;

$r_{\rm H}\equiv a_{\rm p}\left(\frac{m_{\rm p}}{3m_{*}}\right)^{1/3},$ (1)

where $a_{\rm p}$ is the planet semimajor axis and $m_{\rm p}$ and $m_{*}$ are
the planet and star masses respectively. The inner edge of 1 Hill radius
ensures that Trojans are omitted, which would otherwise bias the axisymmetric
surface-density profiles that we fit later (we discuss Trojans in Section
6.1.9). The outer edge of 30 Hill radii ensures that the outer edge does not
affect the inner-edge profile; a circular-orbit planet is expected to scatter
all non-resonant material originating within approximately 3 Hill radii of its
orbit (Gladman 1993; Ida et al. 2000; Kirsh et al. 2009; Malhotra et al. 2021;
Friebe et al. 2022), so setting the outer edge at $30r_{\rm H}$ ensures that
this is well beyond the inner region. Debris semimajor axes are drawn such
that the initial surface-density distribution goes as approximately
${r^{-1.5}}$ (where $r$ is stellocentric distance), akin to the Minimum-Mass
Solar Nebula (MMSN; Weidenschilling 1977; Hayashi 1981). Each debris particle
has an initial eccentricity uniformly drawn between 0 and a maximum value
${e_{\rm max,0}}$, and an initial inclination (relative to the planet’s
orbital plane) uniformly drawn between 0 and ${e_{\rm max,0}/2\;{\rm
radians}\;}$. Each particle’s initial longitude of ascending node, argument of
pericentre and mean anomaly are uniformly drawn between 0 and ${360^{\circ}}$.

To ensure that scattering is essentially complete by the end of the
simulations, we set each simulation end time depending on the diffusion
timescale ${t_{\rm diff}}$. The value ${t_{\rm diff}}$ characterises the
scattering timescale, and is roughly the time taken for a planet to
significantly scatter or eject ${70\;{\rm per\;cent}}$ of material originating
within 3 Hill radii of its orbit (Costa, Pearce & Krivov, submitted). The
diffusion time for a planet acting on a body with semimajor axis $a$ is

${t_{\rm diff}}\approx 0.01T_{\rm p}\left(\frac{a_{\rm
p}}{a}\right)^{1/2}\left(\frac{m_{\rm p}}{m_{*}}\right)^{-2},$ (2)

where $T_{\rm p}$ is the planet’s orbital period (Tremaine, 1993).
Equivalently,

${t_{\rm diff}}\approx 1.1\times 10^{4}\;{\rm yr}\;\left(\frac{a_{\rm
p}}{\;{\rm au}}\right)^{2}\left(\frac{a}{\;{\rm
au}}\right)^{-1/2}\left(\frac{m_{*}}{\;{\rm
M_{\odot}}}\right)^{3/2}\left(\frac{m_{\rm p}}{\;{\rm M_{Jup}}}\right)^{-2}.$
(3)

We run each simulation for at least ${10{t_{\rm diff}}}$ calculated at
${a=a_{\rm p}}$, to ensure that scattering is essentially complete by the end.
Simulations are run in rebound (Rein & Liu, 2012) using whfast (a symplectic
Wisdom-Holman integrator; Rein & Tamayo 2015; Wisdom & Holman 1991), with a
timestep of ${1\;{\rm per\;cent}}$ of the planet’s orbital period111rebound
does not define several orbital parameters in the case of a perfectly circular
orbit. To ensure correct behaviour, for our ‘circular-planet’ simulations we
actually implement a planet eccentricity of ${10^{-4}}$.. Simulations are
conducted and analysed in the centre-of-mass frame.

We run simulations with star masses of 1 to ${2\;{\rm M_{\odot}}}$, planet-to-
star mass ratios of ${3\times 10^{-5}}$ to $10^{-1}$, planet semimajor axes of
1 to ${100\;{\rm au}}$, and maximum initial debris eccentricities of ${e_{\rm
max,0}=0.001}$ to 0.3. In addition, we also re-weight debris particles in
post-simulation analyses, so we can test different initial surface-density
profiles; we consider initial surface-density profiles between ${r^{0}}$ and
${r^{-1.5}}$. We assume the planet and debris to be point-like particles; this
lets us apply simple scaling laws to our results, but also means that the
potential removal of debris via accretion rather than ejection is ignored in
our simulations (see Morrison & Malhotra 2015).

###  2.2 Example simulation results

Figures 2 and 3 show two example simulations after 10 diffusion timescales.
Both have a solar-mass star and a ${2\;{\rm M_{Jup}}}$ planet (${m_{\rm
p}/m_{*}=2\times 10^{-3}}$), which is on a circular orbit at ${10\;{\rm au}}$.
The simulation on Figure 2 starts with an initially unexcited debris disc
(${e_{\rm max,0}=0.01}$), and that on Figure 3 has the same setup but an
initially excited disc (${e_{\rm max,0}=0.1}$). The main qualitative features
of these simulations are typical of all our runs.

Figure 2: Simulation of a circular-orbit planet interacting with an exterior
debris disc, for an initially weakly excited disc. The simulation is of a
solar-mass star, a ${2\;{\rm M_{Jup}}}$ planet on a circular orbit at
${10\;{\rm au}}$ (${m_{\rm p}/m_{*}=2\times 10^{-3}}$), and 20,000 massless
debris particles initially spanning 1 to 30 Hill radii exterior to the planet,
with an initial surface-density profile of approximately ${r^{-1.5}}$. Initial
debris eccentricities were each uniformly drawn between 0 and 0.01, and the
plot shows the final state after ${8.7\times 10^{5}\;{\rm yr}}$ (10 diffusion
timescales). The planet clears nearby debris through scattering, and excites
the eccentricities of surviving material through MMRs. Left panel: semimajor
axes and eccentricities of the planet (large circle) and debris (small points,
coloured by initial semimajor axis). Grey points show debris at the start of
the simulation. Solid black lines are the minimum eccentricities required for
debris to come within 3 Hill radii of the planet, and dashed grey lines are
nominal MMR locations. Right panel: radial surface-density profile of the
disc, $\Sigma(r)$. The thin blue line shows the simulation data, and the thick
orange line is the fit (Equation 6). Insets show the inner-edge region in more
detail. The fitted model parameters are shown on the right panel; $\alpha$ was
fixed to 1.5 in the fit. Figure 3: Simulation with the same setup as that on
Figure 2, except that the debris disc is initially more excited (maximum
initial eccentricities of 0.1, rather than the value of 0.01 on Figure 2). The
planet still excites debris through MMRs (left panel), but these populations
are less significant against the background excitation level of the disc. The
resulting inner edge is flatter than in the low-initial-excitation case.

The left panels of Figures 2 and 3 show the semimajor axes and eccentricities
of the simulated bodies, which demonstrate the two main effects of a circular-
orbit planet on debris. The planet:

  1. 1.

scatters and ejects most debris coming within ${\sim 3}$ Hill radii of its
orbit; this clears the region above the solid black lines on the left panels.

  2. 2.

excites the eccentricities of surviving debris via mean-motion resonances
(MMRs), particularly near the inner edge of the sculpted disc.

The role of MMRs in exciting debris eccentricity is clearest in the initially
unexcited-disc simulation (Figure 2, left panel). Here, populations of debris
in the 2:1 and 3:2 MMRs are clearly visible. For the initially excited disc
(Figure 3), these MMR populations are similar to those in the initially
unexcited simulation, only now they are less distinct due to the higher
intrinsic eccentricity of the disc. The evolution of debris inclination is
much less significant than eccentricity; there is some very small inclination
excitation at specific resonances in these simulations, but the disc still
remains thin across its entire width.

The right panels of Figures 2 and 3 show the azimuthally averaged surface
density profiles of the simulated debris discs (thin blue lines). The planet
imposes a profile on the disc inner edge, whilst the central and outer regions
retain their initial profiles. Individual, strong MMRs can also impose
additional structure, such as the 2:1 MMR at ${15.9\;{\rm au}}$ on Figure 2;
in this example that MMR does not affect the inner-edge fit, but it can do in
other simulations. The thick orange lines are parametric fits to the surface-
density profiles, as detailed in Section 2.3.

Figure 4 shows the final positions of bodies in the initially unexcited-disc
simulation (that on Figure 2). The disc is almost axisymmetric, with a slight
asymmetry between the directions aligned and anti-aligned with the planet due
to MMR structure (an effect noted by Tabeshian & Wiegert 2016, 2017). This
leads to a small difference in the steepness of the inner-edge profile between
the two sides of the disc, which we discuss further in Section 6.2.1. For the
initially excited-disc simulation from Figure 3, the asymmetry is less
pronounced.

Figure 4: Positions of bodies from the simulation with a low-eccentricity disc
(Figure 2). The asterisk, filled black circle and solid line denote the star,
planet and planet orbit respectively, and small points are debris (coloured by
initial semimajor axis). The disc is almost axisymmetric, with a slight
asymmetry between the left and right sides due to MMRs. This causes a small
difference in the inner-edge steepness between the two sides; the right side
(aligned with the planet position) has ${\sigma_{{\rm i}}=0.0286}$, whilst the
value for the opposite side is 0.0354.

###  2.3 Surface-density profile fitting

We aim to quantify how the properties of a sculpting planet affect the
steepness and location of the debris disc’s inner edge. To proceed, we use a
parametric model to fit the radial surface-density profiles of debris in the
$n$-body simulations, and examine how the model parameters change as functions
of system properties. Section 2.3.1 describes the model, Section 2.3.2 the
dependence of the fitted models on system parameters, and Section 2.3.3 the
timescale for a disc edge takes to reach its final state.

####  2.3.1 Surface-density profile model

To quantify the disc surface-density profile, we follow the approach of
Rafikov (2023). They show that, for a low-eccentricity disc with a sharp
semimajor-axis cutoff at the outer edge, the radial surface-density profile
around that edge can be characterised as222Rafikov (2023) derived Equation 4
analytically. Marino (2021) fitted edges numerically using hyperbolic tangents
instead, but the two profiles have similar forms.

$\Sigma(r)\propto 1-{\rm erf}\left(\frac{r-r_{\rm o}}{\sqrt{2}\sigma_{{\rm
o}}r_{\rm o}}\right).$ (4)

Here ${r_{\rm o}}$ characterises the radial location of the outer edge,
$\sigma_{{\rm o}}$ is the ‘flatness’ of the surface-density profile at that
edge, and erf is the Gauss error function:

${\rm erf}(z)\equiv\frac{2}{\sqrt{\pi}}\int^{z}_{0}\exp({-t^{2}}){\rm d}t.$
(5)

Equation 4 is an ‘S-shape’ profile centered on $r_{\rm o}$, which is steeper
for smaller $\sigma_{{\rm o}}$ and flatter for larger $\sigma_{{\rm o}}$. The
value of $\sigma_{{\rm o}}$ is strongly linked to debris eccentricity; if
debris has a sharp cutoff in semimajor axis, and its root-mean-square (rms)
eccentricity is $e_{\rm rms}$, then ${\sigma_{{\rm o}}=e_{\rm rms}/\sqrt{2}}$
(Rafikov, 2023). This means that, for discs with sharp cutoffs in semimajor
axis, lower rms eccentricities mean steeper edges (and smaller $\sigma_{{\rm
o}}$). If the eccentricities are uniformly spread from 0 to $e_{\rm max}$,
then ${e_{\rm rms}=e_{\rm max}/\sqrt{3}}$ and hence ${\sigma_{{\rm o}}=e_{\rm
max}/\sqrt{6}}$. We will later show that a planet-sculpted disc does not have
a sharp cutoff in semimajor axis at the inner edge, so the corresponding
relationship between the edge profile and debris eccentricity is slightly
modified, but the two remain strongly linked.

To parameterise the surface density across a whole disc, we use profiles
similar to Equation 4 at the two edges, combined with an $r^{-\alpha}$ profile
describing the surface density between the edges. By multiplying these three
local profiles together, we arrive at the following model for the overall
surface-density profile:

$\Sigma(r)=\frac{\Sigma_{0}}{2}\left[1-{\rm erf}\left(\frac{r_{\rm
i}-r}{\sqrt{2}\sigma_{{\rm i}}r_{\rm i}}\right)\right]\left[1-{\rm
erf}\left(\frac{r-r_{\rm o}}{\sqrt{2}\sigma_{{\rm o}}r_{\rm
o}}\right)\right]\left(\frac{r}{r_{\rm i}}\right)^{-\alpha},$ (6)

where subscripts ${\rm i}$ and ${\rm o}$ denote terms characterising the inner
and outer edges respectively, and ${\Sigma_{0}}$ is the the surface density at
${r\approx r_{\rm i}}$.

We numerically fit a radial profile of this form to each of our simulated
discs, treating ${\Sigma_{0}}$, $r_{\rm i}$, $r_{\rm o}$, $\sigma_{{\rm i}}$,
$\sigma_{{\rm o}}$ as free parameters. We typically fix $\alpha$ to the
initial surface-density index of the disc (i.e. ${\alpha=1.5}$ for an
initially MMSN profile). The fitting procedure is described in Appendix A, and
the orange lines on the right panels of Figures 2 and 3 show the resulting
fitted profiles for those simulations.

####  2.3.2 Dependence of the final inner-edge profile on system parameters

Having fitted surface-density profiles to each of our $n$-body discs, we now
assess how the profile of a planet-sculpted inner edge depends on system
parameters. The inner-edge profiles are fully characterised by $r_{\rm i}$ and
$\sigma_{{\rm i}}$, and in this section we show how those fitted parameters
vary with simulation setup.

Figure 5 shows how the location of the inner edge of a planet-sculpted disc,
$r_{\rm i}$, scales with the planet-to-star mass ratio, planet semimajor axis,
and the initial disc-excitation level. The location $r_{\rm i}$ is just
exterior to the ‘chaotic zone’ around the planet’s circular orbit, which can
be defined as either ${\sim 3}$ Hill radii or via the Wisdom overlap criterion
(Wisdom, 1980). The value $r_{\rm i}$ is larger for discs with higher initial-
excitation levels, and kinks occur when $r_{\rm i}$ is close to the nominal
location of a strong MMR; this is especially true for discs with low initial
eccentricities. Note that ${r_{\rm i}/a_{\rm p}}$ is independent of the
planet’s semimajor axis, and so depends only on the planet-to-star mass ratio
and the disc’s initial-excitation level.

Figure 5: Location of the sculpted disc’s inner edge $r_{\rm i}$ (in terms of
planet semimajor axis $a_{\rm p}$), from fits to surface-density profiles at
the end of our $n$-body simulations. Coloured points and lines show
simulations with different initial setups, as denoted in the key. Grey solid
lines denote different definitions of the outer edge of the chaotic zone;
after the planet has ejected unstable debris, the disc’s inner edge is
typically just outside this chaotic zone. Dashed horizontal lines denote where
$r_{\rm i}$ coincides with the nominal location of a strong MMR; these
resonances can significantly affect the edge location.

Figure 6 shows how the flatness of the inner-edge profile, $\sigma_{{\rm i}}$,
scales with the planet-to-star mass ratio and the initial disc-excitation
level. It is independent of all other parameters. The inner edges are
generally steeper if sculpted by lower-mass planets, and flatter for higher-
mass planets. The inner edges are also flattest for discs with the highest
initial-excitation levels, although there is less of a dependence for discs
with low initial-excitation levels. The relationships between $\sigma_{{\rm
i}}$ and $m_{\rm p}/m_{*}$ are also not smooth, but show complicated spikes.
These spikes are not numerical effects, because they are replicated across
simulations with the same mass ratios but different planet semimajor axes and
star masses; the spikes actually occur when the inner edge is near the nominal
location of a strong MMR.

Figure 6: Flatness of the sculpted disc’s inner edge, $\sigma_{{\rm i}}$, from
fits to surface-density profiles at the end of our $n$-body simulations. A
smaller value of $\sigma_{{\rm i}}$ corresponds to a sharper edge. Lines and
symbols have the same meanings as on Figure 5. The edge steepness only depends
on the planet-to-star mass ratio and the disc’s initial-excitation level, and
spikes in the plotted relationships occur when the inner edge is near the
nominal location of a strong MMR.

The values of $r_{\rm i}$ and $\sigma_{{\rm i}}$ on Figures 5 and 6 were
fitted for discs with initial surface-density profiles of $r^{-\alpha}$, where
${\alpha=1.5}$, but they are actually independent of $\alpha$ for realistic
setups. This is because the planet-induced edges are typically much steeper
than the overall disc profiles. To check this, we re-scaled the simulated
discs to have initially flat surface-density profiles (i.e. ${\alpha=0}$) and
re-fitted the edges. This resulted in values of $r_{\rm i}$ and $\sigma_{{\rm
i}}$ that are within a few percent of the ${\alpha=1.5}$ values, and hence the
profiles of planet-sculpted edges do not strongly depend on the broader disc
profiles. We further checked this analytically, and show that planet-sculpted
inner edges should be much steeper than the overall disc provided that the
initial-disc profile is flatter than ${r^{-4}}$ (and also ${r^{4}}$); this
calculation is presented in Appendix C. Most measured debris-disc profiles are
shallower than this, so the steepness of a planet-sculpted edge should
typically not depend on the overall disc profile.

####  2.3.3 Edge-sculpting timescale

The timescale for a circular-orbit planet to sculpt the disc inner edge is
expected to scale with the diffusion timescale (Equations 2 and 3), which
quantifies how quickly a planet scatters and ejects debris. We find that this
is reflected in our simulations; after some number of diffusion times, the
inner edges settle into their final configurations. However, the number of
diffusion timescales required appears to have some dependence on the planet-
to-star mass ratio, with different regimes for ratios above and below
$10^{-2}$.

For planet-to-star mass ratios below $10^{-2}$, we find that the inner edge
settles into its final shape within 1 diffusion timescale; therefore, for the
majority of realistic planets, Equations 2 and 3 are reasonable estimates of
the sculpting timescale. However, for mass ratios above $10^{-2}$, more
diffusion timescales are needed; such planets appear to take closer to 10
diffusion times to sculpt the disc. For this reason, we ran all simulations
that had ${m_{\rm p}/m_{*}\geq 10^{-2}}$ until a time of ${100{t_{\rm
diff}}}$, whilst all simulations with ${m_{\rm p}/m_{*}<10^{-2}}$ were run for
${10{t_{\rm diff}}}$. Regardless, despite the difference in the number of
diffusion times required, large planets still sculpt discs much faster than
small planets, because the diffusion timescale strongly decreases with
increasing planet mass.

There may be two reasons for this behaviour change around ${m_{\rm
p}/m_{*}=10^{-2}}$. First, it may mark the transition between a star-planet
interaction, where ${m_{\rm p}/m_{*}\ll 1}$, to a binary interaction, where
the two bodies have comparable mass. In the latter case, dynamical definitions
like the Hill radius break down, so the interaction dynamics may fundamentally
change.

The second reason relates to the ratio of ${t_{\rm diff}}$ to the planet’s
orbital period. For ${m_{\rm p}/m_{*}=10^{-2}}$ the diffusion time is
approximately 100 planet periods, and for ${m_{\rm p}/m_{*}=10^{-1}}$ the
diffusion time is just 1 planet period. For such mass ratios the
approximations used to define the diffusion timescale start to break down. For
example, a planet with ${m_{\rm p}/m_{*}=10^{-1}}$ must take more than one
period (and hence more than one ${t_{\rm diff}}$) to clear unstable material,
because it can only eject material that passes close to it; any material that
is nearly co-orbital with the planet but located on the opposite side of the
star must therefore take several orbital periods (and hence several diffusion
times) to pass close to the planet and get ejected. Conversely, for a small
mass ratio the diffusion time is much longer than the planet period, so such
effects are less important.

###  2.4 Debris eccentricity at the disc inner edge

To understand the profiles of the simulated inner edges, we must understand
the dynamical processes occurring there. Like Marino (2021) and Rafikov
(2023), we expect the edge steepness to be related to debris eccentricity, so
understanding the planet-induced eccentricity is vital for understanding how
planets shape disc inner edges. In this section we quantify the eccentricity
excitation at the simulated inner edges (Section 2.4.1), and show that this
excitation is caused by MMR interactions (Section 2.4.2). Later, in Section 3,
we will use these results to produce a predictive model relating planet
properties to inner-edge profiles.

####  2.4.1 Measuring debris eccentricity at the simulated inner edges

We first directly measure the eccentricities of inner-edge debris at the final
snapshot of each of our $n$-body simulations. To do this, we fit the erf-
powerlaw surface-density model (Equation 6) to each simulated disc, then use
this to define the inner-edge region; the inner edge has a characteristic
width of several times ${r_{\rm i}\sigma_{{\rm i}}}$, so we define the ‘inner-
edge region’ as that centred on $r_{\rm i}$ with an arbitrary full width
${3r_{\rm i}\sigma_{{\rm i}}}$. We then calculate the rms eccentricity of all
debris bodies that are instantaneously located in this radial range, which we
define as the rms eccentricity of the debris-disc inner edge, $e_{\rm i,rms}$.
As an example, for this definition the simulation on Figure 2 has an inner-
edge region spanning from 12.5 to ${13.9\;{\rm au}}$, with ${e_{\rm
i,rms}=0.0433}$.

Figure 7 shows the rms eccentricities at the inner edges of each of our
simulated discs, calculated using the above method. These rms eccentricities
depend only on the planet-to-star mass ratio and the initial-disc excitation
level. For discs with low initial-excitation levels, the eccentricity imparted
by the planet increases with the planet-to-star mass ratio, up to mass ratios
of ${\sim 3.3\times 10^{-3}}$; above this, the inner-edge eccentricity is
roughly independent of mass ratio. Conversely, if the disc’s initial
eccentricity exceeds the level that would be imparted on the edge by the
planet, then the resulting edge eccentricity is essentially independent of the
planet, and remains close to the initial level throughout the simulation.

Figure 7: Root-mean-square eccentricities of surviving debris at the disc
inner edges, at the end of our $n$-body simulations. Lines and symbols are
defined on Figure 5. Dotted horizontal lines are the values expected if the
discs maintain their initial-excitation levels; these should be compared to
simulation lines of the matching colour. The shaded region is a theoretical
prediction for the rms eccentricity acquired by bodies in the planet’s 3:2 or
2:1 MMRs; its width is the spread of these values. The planet excites
surviving debris through MMRs; the black line is a rough fit to the
eccentricity it imparts on debris at the inner edge (Equation 7). If this
planet-induced eccentricity is smaller than the initial-disc excitation level,
then the planet does not significantly increase the excitation level above
that of the pre-interaction disc.

####  2.4.2 Mean-motion resonances as the edge-excitation mechanism for
circular-orbit planets

We now identify the mechanism that excites inner-edge debris. Despite the
significant eccentricities of this material (Figure 2, left panel), this
excitation is not caused by planet-debris scattering, because the semimajor
axes of this debris remains essentially unchanged. Excitation is also not due
to secular interactions, because the low-eccentricity planet cannot
sufficiently excite debris through secular interactions. This leaves MMRs as
the only possible excitation mechanism, in agreement with previous studies
(Quillen, 2006; Quillen & Faber, 2006). The idea that MMRs are responsible is
also supported by the presence of debris populations with similar excited
eccentricities near nominal MMR locations (e.g. the 2:1 MMR on Figure 2). Note
that a single MMR does not usually dominate all surviving debris at the disc
edge; rather, several nearby MMRs excite debris across a range of semimajor
axes.

Given the resonant nature of debris at the inner edge, we can dynamically
explain the dependence of ${e_{\rm i,rms}}$ on planet-to-star mass ratio. The
specific MMRs at the edge depend on the mass ratio, because as the mass ratio
increases, the nominal MMR locations remain unchanged but the width of the
chaotic zone around the planet increases. However, for low mass ratios the
excitation level appears similar to that expected from the 2:1 or 3:2 MMRs,
even if these are not the specific resonances at the edge. To show this, the
grey region on Figure 7 is a theoretical prediction for the rms debris
eccentricity expected for a population of bodies in the 2:1 or 3:2 MMRs, found
by integrating the equations of motion for those MMRs (Pearce et al., in
prep.). This is similar to the ${(m_{\rm p}/m_{*})^{1/3}}$ scaling from
Petrovich et al. (2013) (their Equations 19 and 34), and agrees with our
simulation results for lower planet-to-star mass ratios. Specifically, if a
circular-orbit, non-migrating planet sculpts a debris-disc inner edge, then
the rms eccentricity at that edge increases with planet mass, provided
${m_{\rm p}/m_{*}\lesssim 3.3\times 10^{-3}}$.

For planet-to-star mass ratios above ${\sim 3.3\times 10^{-3}}$, the inner-
edge excitation is below that expected of 2:1 and 3:2 MMRs; this is because
those MMRs would now lie inside the chaotic zone, and higher-order MMRs at
semimajor axes outside the nominal 2:1 location are weaker and less effective
at exciting debris (Figure 2). Hence for planet-to-star mass ratios above
${\sim 3.3\times 10^{-3}}$, the debris eccentricity at the disc edge does not
further increase with mass ratio.

Given this behaviour, we can roughly quantify the eccentricity that a planet
on a circular orbit induces on initially unexcited debris at the disc inner
edge:

$e_{\rm i,rms,p}\approx\begin{cases}0.4\left(m_{\rm
p}/m_{*}\right)^{1/3},&\text{if $m_{\rm p}/m_{*}\leq 3.3\times 10^{-3}$};\\\
0.06,&\text{else}.\end{cases}$ (7)

This is the black line on Figure 7, which is a good match to the simulations
where the planet-induced excitation dominates over the intrinsic disc
excitation (e.g. the purple lines for mass ratios above $10^{-4}$).
Conversely, since the initial-disc excitation would dominate over the planet-
induced excitation if the former were high enough, we can generally predict
the eccentricity of surviving debris at the inner edge of a planet-sculpted
debris disc as

$e_{\rm i,rms}=\max\left(e_{\rm i,rms,p},e_{\rm i,rms,0}\right).$ (8)

where ${e_{\rm i,rms,0}}$ is the ‘intrinsic’ rms eccentricity at the disc edge
(i.e. the pre-interaction eccentricity), and the planet-induced eccentricity
$e_{\rm i,rms,p}$ is given by Equation 7.

##  3 Analytic model relating the inner-edge profile to the sculpting planet

In Section 2 we showed that, if a circular-orbit planet sculpts a
collisionless debris disc, then the steepness and relative location of the
disc’s inner edge depend only on the planet-to-star mass ratio and the disc’s
initial-excitation level. We also quantified how the planet excites debris
eccentricities at the inner edge. Using these results, we now produce a simple
analytical model to infer the properties of sculpting planets from inner-edge
profiles. Section 3.1 details our analytical model, and Section 3.2
demonstrates how to use it to infer the properties of an unseen sculpting
planet from an observed inner-edge profile.

###  3.1 Model setup and predictions

We assume a simplified model of the inner-edge region, as shown on Figure 8.
We consider a planet on a circular orbit, and a debris disc that initially
extends from the planet’s semimajor axis out to some larger distance. The
debris particles have initial eccentricities uniformly distributed between 0
and some $e_{\rm max}$; for this model, the origin of these eccentricities is
unimportant. To predict the location and steepness of such a planet-sculpted
disc edge, we consider which particles would be scattered by the planet, and
which would survive.

Figure 8: Schematic of our simple model for predicting the disc’s inner-edge
profile (Section 3). The black circle is the planet, and the red region
surviving debris. Debris has eccentricities between 0 and ${e_{\rm max}}$. The
planet eventually ejects all debris originating within its chaotic zone, which
is any debris above the solid line (Equations 9 and 10). The value ${a_{1}}$
is the smallest semimajor axis of any surviving particle (a circular orbit
just exterior to the chaotic zone), and ${a_{2}}$ is the largest semimajor
axis of any unstable particle (an orbit with the highest eccentricity ${e_{\rm
max}}$, whose pericentre is at the outer edge of the chaotic zone).

In our model, any particle whose orbit comes within the planet’s chaotic zone
would be scattered and eventually ejected. This means that, at late times, no
particles should occupy the parameter space above the solid line on Figure 8;
this line is the eccentricity ${e_{\rm q}(a)}$ that results in a particle’s
pericentre coinciding with the outer edge of the chaotic zone, defined by

$e_{\rm q}(a)\equiv 1-\frac{a_{\rm p}+\delta a}{a}.$ (9)

Here $a$ is the semimajor axis, and $\delta a$ is the half-width of the
chaotic zone, taken as 3 Hill radii for a circular-planet orbit:

$\delta a\equiv 3r_{\rm H}=3a_{\rm p}\left(\frac{m_{\rm
p}}{3m_{*}}\right)^{1/3}.$ (10)

Since the planet eccentricity is zero, any non-scattered debris would retain
its initial semimajor axis and eccentricity. Hence any debris with initial
eccentricity above the Equation 9 line would eventually get ejected, whilst
any with eccentricity below the line would remain unperturbed. So in our
model, at late times debris only occupies the shaded region of Figure 8.

In the following sections, we use this model to infer the location and
steepness of the planet-sculpted inner edge. Before doing so, we must first
define two final parameters: the semimajor axes $a_{1}$ and $a_{2}$, as shown
on Figure 8. The span $a_{1}$ to $a_{2}$ roughly defines the inner-edge region
in semimajor-axis space. The values of $a_{1}$ and $a_{2}$ correspond to the
semimajor axes of orbits at the edge of the chaotic zone; $a_{1}$ corresponds
to a circular orbit, and $a_{2}$ to an orbit with eccentricity equal to the
maximum debris eccentricity $e_{\rm max}$. Expressions for $a_{1}$ and $a_{2}$
are given in Appendix D.1.

####  3.1.1 Predicted radial location of the disc inner edge, $r_{\rm i}$

We now use the simple model on Figure 8 to predict $r_{\rm i}$, the
characteristic location of a planet-sculpted inner edge. A full derivation of
the $r_{\rm i}$ prediction is presented in Appendix D.2; our basic method is
to first calculate the time-averaged radial location of all debris at a single
semimajor axis $a$, and then average these for all debris with ${a_{1}\leq
a\leq a_{2}}$. This yields our prediction:

$r_{\rm i}\approx a_{\rm p}\left[1+3^{2/3}\left(\frac{m_{\rm
p}}{m_{*}}\right)^{1/3}\right]\left(1+\frac{2}{\sqrt{3}}e_{\rm i,rms}\right),$
(11)

where we used ${e_{\rm i,rms}=e_{\rm max}/\sqrt{3}}$ for a uniform $e$
distribution to convert $e_{\rm max}$ into $e_{\rm i,rms}$.

Equation 11 predicts that the edge location depends on planet semimajor axis,
planet mass, and debris eccentricity. The origin of this eccentricity is
unimportant, provided that the eccentricity-semimajor axis distribution
resembles that on Figure 8. We can therefore evaluate the eccentricity term in
Equation 11 using Equations 7 and 8. If the debris eccentricity is set by the
intrinsic (i.e. pre-interaction) eccentricity of the disc, which would occur
if the intrinsic eccentricity is higher than that induced by the planet, then
$e_{\rm i,rms}$ in Equation 11 can simply be replaced by the intrinsic disc
eccentricity $e_{\rm i,rms,0}$. Alternatively, for a disc with sufficiently
low pre-interaction eccentricity, we can substitute ${e_{\rm i,rms}=e_{\rm
i,rms,p}}$ from Equation 7 to predict the location of the planet-sculpted
inner edge:

$r_{\rm i}\approx\begin{cases}a_{\rm p}\left[1+2.5\left(\frac{m_{\rm
p}}{m_{*}}\right)^{1/3}\right],&\text{if $m_{\rm p}/m_{*}\leq 3.3\times
10^{-3}$};\vspace{2mm}\\\ a_{\rm p}\left[1.1+2.2\left(\frac{m_{\rm
p}}{m_{*}}\right)^{1/3}\right],&\text{else},\end{cases}$ (12)

where we omit higher-order terms in $\left(m_{\rm p}/m_{*}\right)^{1/3}$.

Having made a prediction for $r_{\rm i}$, we now compare this prediction to
the results of our $n$-body simulations. For each simulation, we predict
$r_{\rm i}$ by evaluating Equation 11; to do this, we take $a_{\rm p}$ and
$m_{\rm p}$ from the simulation setup, and also use the value of $e_{\rm
i,rms}$ measured directly from the simulation. The results are shown on Figure
9. We see that the prediction generally works well; the predicted value of
${r_{\rm i}-a_{\rm p}}$ for simulations with $e_{\rm max,0}<0.3$ is typically
within ${10\;{\rm per\;cent}}$ of that from simulations. However, the
prediction is less accurate for discs with very high intrinsic eccentricities
($e_{\rm max,0}=0.3$) that are interacting with low-mass planets (${m_{\rm
p}/m_{*}\lesssim 10^{-4}}$); for such high eccentricities there is
considerable overlap of debris orbits, so our simple $r_{\rm i}$ approximation
using only the average debris position probably no longer holds. There is also
a noticeable divergence for planet-to-star mass ratios above ${10^{-2}}$,
which could be due to a fundamental shift in the dynamics; this is discussed
in Section 2.3.3. Nonetheless, Equation 11 holds across the large majority of
our explored parameter space, so offers a reasonable means to predict the
location of a planet-sculpted inner edge.

Figure 9: Comparison of our predicted inner-edge location to the results of
$n$-body simulations. The vertical axis is the difference between the planet
orbit and the inner-edge location $r_{\rm i}$ fitted from the simulations,
divided by the equivalent value predicted by Equation 11. The horizontal grey
line denotes unity; all other symbols and lines are defined on Figure 5. Our
predicted edge locations are generally similar to those from the simulations,
with the predictions typically within ${10\;{\rm per\;cent}}$ of the simulated
values for planets below ${1\;{\rm per\;cent}}$ of the stellar mass and
initial debris-excitation levels of $e_{\rm max,0}<0.3$. Predictions are worse
for larger mass ratios, and also for simulations with both high initial-debris
eccentricities and low-mass planets.

####  3.1.2 Predicted steepness of the disc inner edge, $\sigma_{{\rm i}}$

In this section we use our simple model to predict the shape of the planet-
sculpted inner edge, as quantified by $\sigma_{{\rm i}}$. Near the inner edge,
our fitted profile (Equation 6) is essentially an erf function; this has
characteristic width ${xr_{\rm i}\sigma_{{\rm i}}}$, where $x$ is a scalar of
order unity. We hypothesise that this width should be roughly equivalent to
the distance between the pericentre of the innermost stable particle and the
apocentre of the outermost unstable particle, i.e. ${xr_{\rm i}\sigma_{{\rm
i}}\approx a_{2}(1+e_{\rm max})-a_{1}}$ from Figure 8. Substituting
expressions for $r_{\rm i}$, $a_{1}$ and $a_{2}$ yields ${\sigma_{{\rm
i}}\propto e_{\rm max}}$ to first order in $e_{\rm max}$, i.e. the steepness
of a planet-sculpted inner edge depends entirely on the eccentricity of debris
around that edge. Marino (2021) and Rafikov (2023) reached a similar
conclusion for edges with sharp cutoffs in semimajor axis, but we show that
this proportionality is also expected from the smooth distribution arising
from scattering.

Since we expect ${\sigma_{{\rm i}}\propto e_{\rm max}}$, we use our $n$-body
simulations to directly relate $\sigma_{{\rm i}}$ to debris eccentricity.
Figure 10 shows the rms eccentricities measured at the inner edges of our
simulated discs, divided by the $\sigma_{{\rm i}}$ values fitted to those
simulations. This ratio is roughly constant for all simulations, verifying our
prediction that $\sigma_{{\rm i}}$ depends only on $e_{\rm max}$ (and hence
$e_{\rm i,rms}$). Taking the median of this ratio from our simulations, we
find that the flatness of the inner edge of a planet-sculpted disc depends on
the eccentricity of debris at that edge as

$e_{\rm i,rms}\approx(1.2\pm 0.2)\;\sigma_{{\rm i}}.$ (13)

Here the uncertainty on the empirical prefactor is defined from the inter-
quartile range on ${e_{\rm i,rms}/\sigma_{{\rm i}}}$ from our simulations.

Figure 10: Relating the steepness of sculpted-disc inner edges to the debris
eccentricities at those edges, from our $n$-body simulations. Smaller
$\sigma_{{\rm i}}$ values correspond to steeper edges. Our simple model
predicts $\sigma_{{\rm i}}$ and $e_{\rm i,rms}$ to be proportional, in
agreement with the simulations. The horizontal grey line and shaded region
show the value and uncertainty of the prefactor ${1.2\pm 0.2}$ in Equation 13,
which is empirically measured from the simulations. All other lines and
symbols are defined on Figure 5. The breakdown in the fit at planet-to-star
mass ratios between $10^{-3}$ and $10^{-2}$ for low-eccentricity simulations
is due to the 2:1 MMR imposing additional structure on the edge (discussed in
Section 2.3.2).

###  3.2 Inferring the properties of a sculpting planet from the location and
shape of a debris-disc inner edge

Section 3.1 related the disc inner edge to the properties of a sculpting
planet. Since it has historically been easier to resolve a cold debris disc
than to detect a distant planet, it is common to use observed debris discs to
infer the properties of unseen planets (e.g. Pearce et al. 2022). In this
section we provide a method to constrain a sculpting planet on a circular
orbit from the shape and location of a debris-disc inner edge.

To demonstrate the method, we apply it to the simulation on Figure 2 as a
example. We assume that the simulated disc has been observed, but that the
sculpting planet is undetectable. We will constrain the properties of the
unseen planet from the disc alone, then compare these to the known parameters
of the simulated planet to gauge the effectiveness of the method. The various
steps in the calculation are shown on Figure 11, and described below. For this
example we assume that the planet has finished sculpting the disc by the time
the observations are made.

Figure 11: Example application of our results: using a debris disc’s inner-
edge profile to infer the properties of an unseen sculpting planet. The plot
shows the disc from the $n$-body simulation on Figure 2, for which we attempt
to constrain the planet using the step-by-step method from Section 3.2.1. The
disc has a fitted inner-edge radius ${r_{\rm i}=13.2\;{\rm au}}$ and flatness
${\sigma_{{\rm i}}=0.0361}$. The dashed vertical line is $r_{\rm i}$, and the
light brown region around it is the approximate full width of the inner edge
($3r_{\rm i}\sigma_{{\rm i}}$). The large black circle shows the planet
properties we infer, and the red diamond is the actual planet in the $n$-body
simulation. The white square is the planet’s minimum-possible mass and maximum
semimajor axis found using the alternative method of Pearce et al. (2022),
which accounts for the edge location but not its steepness. Line 1 shows the
planet properties required to set the edge at $r_{\rm i}$ (Equation 14), and
Line 2 is the minimum planet mass required to have sculpted the disc edge
within the age of the star (Equation 15); for this example the age is
arbitrarily taken as ${100\;{\rm Myr}}$.

####  3.2.1 Step-by-step method

The first step is to fit the disc surface-density profile with some function
that quantifies the location and steepness of the inner edge. In this paper we
use an erf-powerlaw function (Equation 6), which yields inner-edge location
$r_{\rm i}$ and flatness $\sigma_{{\rm i}}$; other parameterisations could
alternatively be used, and in Appendix E we provide equations to convert the
outputs of common parametric models into $r_{\rm i}$ and $\sigma_{{\rm i}}$
values. For our example simulation, the fitted inner-edge profile has ${r_{\rm
i}=13.2\;{\rm au}}$ and ${\sigma_{{\rm i}}=0.0361}$.

The second step, now that we have the edge profile, is to infer the sculpting-
planet mass as a function of its semimajor axis. This step makes the implicit
assumption that the planet has a circular orbit and has not migrated, but it
otherwise applies regardless of whether the planet or some other process is
responsible for debris eccentricities. Equation 11 relates planet mass and
semimajor axis to the location and rms eccentricity of the disc inner edge;
using this equation, and noting the relation between $e_{\rm i,rms}$ and
$\sigma_{{\rm i}}$ (Equation 13), we can infer the sculpting-planet mass using

$m_{\rm p}\approx 116\;{\rm M_{Jup}}\left(\frac{m_{*}}{\;{\rm
M_{\odot}}}\right)\left[\frac{r_{\rm i}}{a_{\rm p}}\left(1+1.35\sigma_{{\rm
i}}\right)^{-1}-1\right]^{3}.$ (14)

This is Line 1 on Figure 11 (noting that ${m_{*}=1\;{\rm M_{\odot}}}$ in our
example).

The next step is to put a lower bound on the planet mass, assuming that the
planet has finished sculpting the observed-disc edge. This means that the
sculpting timescale must be shorter than the star age $t_{*}$. In Section
2.3.3 we showed that the sculpting timescale is approximately the diffusion
time if ${m_{\rm p}/m_{*}\lesssim 10^{-2}}$, and ${\sim 10}$ diffusion times
otherwise; we can therefore rearrange Equation 3 and evaluate it at ${a=r_{\rm
i}}$ to show that

$m_{\rm p}\gtrsim 0.105\sqrt{k}\;{\rm M_{Jup}}\left(\frac{t_{*}}{\;{\rm
Myr}}\right)^{-1/2}\left(\frac{a_{\rm p}}{\;{\rm
au}}\right)\left(\frac{m_{*}}{\;{\rm
M_{\odot}}}\right)^{3/4}\left(\frac{r_{\rm i}}{\;{\rm au}}\right)^{-1/4},$
(15)

where

$k\approx\begin{cases}1,&\text{if $m_{\rm p}/m_{*}\lesssim 10^{-2}$};\\\
10,&\text{else}.\end{cases}$ (16)

For our example we ascribe an arbitrary stellar age of ${100\;{\rm Myr}}$,
which results in Line 2 on Figure 11.

The final step is to use the edge shape to break the degeneracy between planet
mass and semimajor axis. Unlike previous steps, this step requires the
implicit assumption that the planet is solely responsible for exciting debris.
It is also only valid if ${\sigma_{{\rm i}}\lesssim 0.05}$ (i.e. ${e_{\rm
i,rms}\lesssim 0.06}$); this is the flattest profile that a non-migrating,
circular-orbit planet could impart on an initially low-eccentricity disc
(Figure 7), so if ${\sigma_{{\rm i}}\gtrsim 0.05}$ then some other process
must have excited debris and this step cannot be applied. If ${\sigma_{{\rm
i}}\lesssim 0.05}$, then rearranging Equation 7 yields the planet mass as

$m_{\rm p}\approx 2.64\times 10^{4}\;{\rm M_{Jup}}\left(\frac{m_{*}}{\;{\rm
M_{\odot}}}\right)\sigma_{{\rm i}}^{3};$ (17)

we can then substitute this mass into Equation 14 to yield the planet’s
semimajor axis as

$a_{\rm p}\approx r_{\rm i}\left(1+1.35\sigma_{{\rm
i}}\right)^{-1}\left[1+0.205\left(\frac{m_{\rm p}}{\;{\rm
M_{Jup}}}\right)^{1/3}\left(\frac{\;{\rm m}_{*}}{\;{\rm
M_{\odot}}}\right)^{-1/3}\right]^{-1}.$ (18)

For our example, this yields a planet mass of ${1.24\;{\rm M_{Jup}}}$ and
semimajor axis ${10.3\;{\rm au}}$; these are in good agreement with the actual
values of ${2\;{\rm M_{Jup}}}$ and ${10\;{\rm au}}$ from the simulation, as
shown on Figure 11. Note that caution must be applied if ${\sigma_{{\rm
i}}\sim 0.05}$ (i.e. ${e_{\rm i,rms}\sim 0.06}$), because such an edge profile
could be imparted by any planet with ${m_{\rm p}/m_{*}\gtrsim 3.3\times
10^{-3}}$ (Figure 7); in this case, the mass estimate from Equation 17 should
be interpreted as a lower bound rather than a single value.

####  3.2.2 Accuracy of planet parameters inferred from disc edges

The above example showed that planet parameters can be well inferred from
inner-edge profiles in at least some setups. We now repeat the above process
for a large fraction of our simulations, to test how well the method applies
in general. We omit simulations where the initial disc-excitation level is
larger than the expected excitation generated by the planet, because for those
simulations, Equation 17 cannot be used to infer planet mass. The results are
shown on Figure 12. The top panel shows that the inferred debris-excitation
level $e_{\rm i,rms}$, calculated from the edge profile using Equation 13,
agrees with the actual simulation values in almost all cases; the largest
discrepancy is for planet-to-star mass ratios around ${5\times 10^{-3}}$,
which is where additional structure from the strong 2:1 MMR coincides with the
disc edge. The bottom panel shows that the planet masses inferred using
Equation 17 agree with the actual planet masses for low mass ratios, but that
the two diverge if the actual planet has ${m_{\rm p}/m_{*}\gtrsim 3.3\times
10^{-3}}$. For the highest mass ratios, Equation 17 can underpredict the
planet mass by over one order of magnitude. The reason for this is that, for
planets with mass ratios above ${3.3\times 10^{-3}}$, the level of
eccentricity excitation imparted on the disc edge is independent of planet
mass (Figure 7); this can cause Equation 17 to significantly underpredict the
planet mass if debris at the disc edge has an excitation level of ${e_{\rm
i,rms}\approx 0.06}$. Note that this degeneracy only affects the use of
Equation 17; regardless of the mass ratio, and whether or not the planet was
responsible for exciting debris, a sculpting planet should still lie close to
the mass-semimajor axis relation from Equation 14 (Line 1 on Figure 11).

Figure 12: Planet properties inferred from simulated-disc inner edges using
the process in Section 3.2, compared to the actual planets. Top panel:
${e_{\rm i,rms}}$ values inferred from our fitted inner-edge profiles
(Equation 13), which agree well with the actual values. Bottom panel: planet
masses inferred from the edge steepnesses (Equation 17), which agree well with
the actual values for planet-to-star mass ratios below ${3.3\times 10^{-3}}$
(equivalent to ${3.5\;{\rm M_{Jup}}}$ around a Sun-like star). For higher
ratios, Equation 17 can significantly underpredict the planet mass, because
the edge excitation level becomes constant with planet mass (Figure 7). Note
that this degeneracy only affects planet masses that were directly inferred
from the edge steepness (black circle on Figure 11); regardless of this, all
planets would still lie along the mass-semimajor axis relation (Line 1 on
Figure 11). Blue circles and red squares are simulations where the inferred
debris excitation ${e_{\rm i,rms}}$ at the disc inner edge is below and above
0.06, respectively (where ${e_{\rm i,rms}\approx 0.06}$ corresponds to an edge
profile with ${\sigma_{{\rm i}}\approx 0.05}$).

##  4 Collisions

Sections 2 and 3 assessed how a circular-orbit planet would sculpt the inner
edge of a debris disc, if the interaction were purely $n$-body dynamics. For
that scenario we showed that the edge steepness could be directly related to
the planet mass. However, in a real disc there are also collisions between
debris bodies, which would alter the edge profile. In this section we consider
the collisional evolution, and answer two basic questions: could a collision-
dominated disc edge be differentiated from a planet-dominated edge, and how
would collisions change the profile of a planet-sculpted edge?

###  4.1 Could a purely collisional disc be differentiated from a purely
planet-sculpted disc?

Destructive collisions cause debris discs to lose mass over time; material
grinds down to dust, which is eventually removed by radiation forces (Wyatt et
al., 1999; Wyatt, 2008; Krivov, 2010). The rate at which this process occurs
depends on several main factors: the number density of debris bodies, their
size distribution, and their relative velocities.

For an initially broad debris disc, these factors result in well-defined
collisional evolution. At early times the surface density is expected to
decrease with distance, for example going as ${r^{-1.5}}$ in the MMSN.
Relative velocities are also higher in the inner regions, and this combination
makes the initial collisional depletion faster closer to the star. However,
this means that the surface density closer to the star drops faster than that
further out, reducing the inner collision rate. The result is that the
initially negative surface-density slope gradually transforms into a positive
one; a ‘collision front’ moves outwards across the disc, where the region
interior to the front tends towards a positive surface-density slope, whilst
the exterior region still has its original profile. In the collisionally
processed region, this profile is well-characterised as a powerlaw $r^{p_{\rm
i}}$, with index ${{p_{\rm i}}=7/3=2.3}$ in Wyatt et al. (2007) and Kennedy &
Wyatt (2010), or ${{p_{\rm i}}=2}$ in more-recent models (Löhne et al., 2008;
Imaz Blanco et al., 2023).

Conversely, if a planet sculpts a collisionless disc, then the inner-edge
steepness is set by the planet mass. To answer whether such a planet-dominated
disc can be differentiated from a collision-dominated disc, we must compare
the inner-edge slopes. Using Equation 53 to relate an erf-fitted slope to a
powerlaw slope, we see that a planet-dominated disc will have a steeper inner
edge than a collision-dominated disc if ${\sigma_{{\rm i}}<0.40}$ (assuming
${{p_{\rm i}}=2}$ for the collisional disc; if ${{p_{\rm i}}=7/3}$ were used
instead, then ${\sigma_{{\rm i}}<0.35}$ for the sculpted edge to be steeper).
Figure 6 shows that these conditions are satisfied by the end of all of our
circular-planet simulations, meaning that the inner edge of a purely planet-
sculpted disc is steeper than that of a purely collisional disc if sculpting
has completed. In some cases it would be much steeper; our steepest $n$-body
discs have ${\sigma_{{\rm i}}\sim 10^{-2}}$, which would correspond to
${{p_{\rm i}}>10}$ in $r^{p_{\rm i}}$ notation. So we should be able to
differentiate purely planet-sculpted discs from purely collision-dominated
discs, with the former having much steeper inner edges.

###  4.2 How would collisions affect a planet-sculpted edge?

The previous section compared a collisionless, planet-sculpted disc to a
collisional disc without a planet. However, in reality a planet-sculpted disc
would also undergo collisions. In this section we consider collisions in more
detail, to assess how they would affect a planet-sculpted disc.

A fully self-consistent model, where a disc undergoes $n$-body interactions
whilst simultaneously collisionally evolving, is difficult to properly
implement. In particular, there are several parameters that are critical for
quantifying the collisional evolution, but that are fundamentally unknown for
real debris discs; these include the debris-disc mass, and the sizes of the
largest planetesimals (Krivov & Wyatt, 2021). For this reason we chose to
decouple the dynamical and collisional modelling in this paper. Instead, we
will take the disc morphologies arising from our $n$-body simulations, and
then input these into a collisional model to estimate how collisions would
modify a planet-sculpted disc edge. This approach is not fully self-
consistent, but would be valid if the planet-sculpting timescale were much
shorter than the collisional timescale; we will later show that this condition
holds for many plausible scenarios.

####  4.2.1 Collisional model

To model collisions, we use a mathematica implementation of the Löhne et al.
(2008) collisional model. The model takes an axisymmetric surface-density
profile, and splits it into radial bins. Each bin is then treated separately,
and the dust mass in each bin is determined after some time. We then use the
masses in each radial bin to calculate the dust surface-density profile. In
this section we briefly describe the model and the parameters we use; we refer
the reader to Löhne et al. (2008) for a much more detailed description.

We assume the eccentricity $e$ and inclination $I$ of the debris population
are related by $I\approx e/2$, and treat $e$ as a free parameter. Following
Löhne et al. (2012), Schüppler et al. (2016) and Krivov et al. (2018), we
assume a critical fragmentation energy of

$Q_{\rm D}^{*}=\left(\frac{v_{\rm col}}{v_{0}}\right)^{1/2}\left[A_{\rm
s}\left(\frac{s}{1\;{\rm m}}\right)^{3{b_{\rm s}}}+A_{\rm
g}\left(\frac{s}{1\;{\rm km}}\right)^{3{b_{\rm g}}}\right]$ (19)

with ${v_{0}=3\;{\rm km\;s^{-1}}}$, ${A_{\rm s}=A_{\rm g}=5\times 10^{6}\;{\rm
erg\;g^{-1}}}$, ${b_{\rm s}=-0.12}$ and ${b_{\rm g}=0.46}$. The planetesimal-
collision speed $v_{\rm col}$ is given by

$v_{\rm col}=\sqrt{\frac{\mathcal{G}m_{*}}{r}}f(e,I),$ (20)

where $\mathcal{G}$ is the gravitational constant and

$f(e,I)=\sqrt{\frac{5}{4}e^{2}+I^{2}}.$ (21)

The value of $Q_{\rm D}^{*}$ is dominated by the material strength for smaller
bodies, and gravity for larger bodies. The transition between these can be
defined either as the size at which the strength and gravity terms in Equation
19 are equal, or as the size at which $Q_{\rm D}^{*}$ is minimised. Following
Löhne et al. (2008), here we use the former definition. For our parameters,
the strength-gravity transition then occurs at the “breaking” radius ${{s_{\rm
b}}=232\;{\rm m}}$.

We assume that debris in each bin follows a three-slope size distribution as
in Löhne et al. (2008), with slope ${{q_{\rm p}}=1.87}$ for primordial bodies,
${{q_{\rm g}}=1.68}$ in the gravity-dominated quasi-steady state and ${{q_{\rm
s}}=1.89}$ in the strength-dominated quasi-steady state. We choose minimum and
maximum dust-grain radii of ${s_{\rm min}=2\;{\rm\mu m}}$ and ${{s_{\rm
dust}}=1\;{\rm mm}}$, and treat the largest-planetesimal radius $s_{\rm max}$
as a free parameter. We use a bulk density of solids ${\rho=3.5\;{\rm
g\;cm^{-3}}}$.

According to the model, there are two key timescales that determine the
collisional evolution. The shorter, $\tau_{\rm b}$, is the time when the
weakest bodies begin to collide. At this time, the collisional decay of the
dust density sets in (if ${q_{\rm p}}\leq{q_{\rm s}}$) or at least speeds up
(if ${q_{\rm p}}>{q_{\rm s}}$); see Equation 43 and Figure 9 in Löhne et al.
(2008). The value of $\tau_{\rm b}$ is given by Equation 31 in Löhne et al.
(2008):

$\tau_{\rm b}=\frac{16\pi\rho}{3M_{0}}\left(\frac{{s_{\rm b}}}{s_{\rm
max}}\right)^{3{q_{\rm p}}-5}\frac{s_{\rm
max}r^{7/2}}{\sqrt{\mathcal{G}m_{*}}}\frac{\delta r}{r}\\\ \times\frac{{q_{\rm
p}}-5/3}{2-{q_{\rm p}}}\left[1-\left(\frac{s_{\rm min}}{s_{\rm
max}}\right)^{6-3{q_{\rm p}}}\right]\frac{I}{f(e,I)G({q_{\rm p}},{s_{\rm
b}},r)},$ (22)

where $M_{0}$ is the total initial mass of solids in the bin, ${\delta r}$ is
the radial width of the bin,

$G(q,s,r)=\left[X_{\rm c}(s,r)^{5-3q}-\left(\frac{s_{\rm
max}}{s}\right)^{5-3q}\right]\\\ +2\frac{q-5/3}{q-4/3}\left[X_{\rm
c}(s,r)^{4-3q}-\left(\frac{s_{\rm max}}{s}\right)^{4-3q}\right]\\\
+\frac{q-5/3}{q-1}\left[X_{\rm c}(s,r)^{3-3q}-\left(\frac{s_{\rm
max}}{s}\right)^{3-3q}\right],$ (23)

and

$X_{\rm c}(s,r)=\left[\frac{2Q_{\rm
D}^{*}(s,r)r}{f^{2}(e,I)\mathcal{G}m_{*}}\right]^{1/3}.$ (24)

We arbitrarily use ${\delta r/r=0.1}$ in our code, but this later cancels out
so its value does not affect any results.

The longer timescale is $\tau_{\rm max}$: the collisional lifetime of the
largest planetesimals. This is Equation 42 in Löhne et al. (2008)333Equation
42 in Löhne et al. (2008) has an erroneous index of -1 around the square
bracket, which we omit here; this does not affect our results because ${s_{\rm
min}\ll s_{\rm max}}$.:

$\tau_{\rm max}=\frac{16\pi\rho}{3M_{0}}\frac{s_{\rm
max}r^{7/2}}{\sqrt{\mathcal{G}m_{*}}}\frac{\delta r}{r}\\\ \times\frac{{q_{\rm
g}}-5/3}{2-{q_{\rm p}}}\left[1-\left(\frac{s_{\rm min}}{s_{\rm
max}}\right)^{6-3{q_{\rm p}}}\right]\frac{I}{f(e,I)G({q_{\rm g}},s_{\rm
max},r)}.$ (25)

We will later show that $\tau_{\rm max}$ is the most important factor for
determining the collisional behaviour. We provide an open-access python code
to calculate $\tau_{\rm max}$ given the assumptions in this
paper444http://www.tdpearce.uk/public-code, and in Section 4.2.3 we write
$\tau_{\rm max}$ in an alternative form to make its dependencies more
explicit.

Given the above setup, the model yields the dust mass $m_{\rm dust}$ in each
radial bin after some time $t$. This is Equation 43 in Löhne et al. (2008):

$m_{\rm dust}(t)=\frac{M_{0}}{1+t/\tau_{\rm max}}\mathcal{T}^{({q_{\rm
g}}-{q_{\rm p}})/[{q_{\rm p}}-5/3+({q_{\rm p}}-1){b_{\rm g}}]}\frac{2-{q_{\rm
p}}}{2-{q_{\rm s}}}\\\ \times\left(\frac{{s_{\rm b}}}{s_{\rm
max}}\right)^{6-3{q_{\rm p}}}\left[\left(\frac{{s_{\rm dust}}}{{s_{\rm
b}}}\right)^{6-3{q_{\rm s}}}-\left(\frac{s_{\rm min}}{{s_{\rm
b}}}\right)^{6-3{q_{\rm s}}}\right]^{-1}$ (26)

for ${\tau_{\rm b}<t<\tau_{\rm max}}$, where ${\mathcal{T}=t/\tau_{\rm b}}$ if
${t<\tau_{\rm max}}$ or ${\tau_{\rm max}/\tau_{\rm b}}$ otherwise555Equation
43 in Löhne et al. (2008) has erroneous indices of ${2-q}$ instead of
${6-3q}$, for ${q_{\rm p}}$ and ${q_{\rm s}}$.. For ${t<\tau_{\rm b}}$,
${q_{\rm g}}$ and ${b_{\rm g}}$ should be replaced by ${q_{\rm s}}$ and
${b_{\rm s}}$ respectively.

For each collisional simulation, we initialise the disc with a surface-density
profile similar to the outcomes of our $n$-body simulations. This is a
powerlaw with a truncated inner edge:

$\Sigma(r)={x_{\rm m}}\left(\frac{r}{\;{\rm au}}\right)^{-\alpha}\frac{\;{\rm
M_{\oplus}}}{\;{\rm au}^{2}}\times\frac{1}{2}\left[1-{\rm
erf}\left(\frac{r_{\rm i}-r}{\sqrt{2}\sigma_{{\rm i}}r_{\rm
i}}\right)\right],$ (27)

where a disc with ${{x_{\rm m}}=1}$ and ${\alpha=3/2}$ would have the MMSN
profile beyond the inner edge. This initial profile sets $M_{0}$ in the above
equations. The initial values of ${{x_{\rm m}}}$ and ${\alpha}$ are free
parameters that we vary between models. After setting up the disc, we use
Equation 26 to determine the dust mass in each radial bin after some time.
Using these bin masses, we compute the dust surface-density profile. Finally,
we fit this profile with an erf-powerlaw function similar to Equation 27, to
measure $r_{\rm i}$, $\sigma_{{\rm i}}$ and $\alpha$ at that time.

####  4.2.2 Collision results

The previous section described the collisional model. In this section we apply
the model to assess the impact of collisions on planet-sculpted disc edges.

We first apply it to the $n$-body disc from Figure 2 as an example. We
initialise the disc with ${{x_{\rm m}}=1}$ and ${\alpha=3/2}$, with an inner
edge defined by ${r_{\rm i}=13.2\;{\rm au}}$ and ${\sigma_{{\rm i}}=0.0361}$,
and set the largest-planetesimal radius to ${s_{\rm max}=50\;{\rm km}}$. We
also set the debris eccentricity to 0.043, the rms eccentricity at the disc
inner edge in the $n$-body simulation. This results in ${\tau_{\rm
max}=1.1\;{\rm Gyr}}$ and ${\tau_{\rm b}=0.14\;{\rm Myr}}$ at ${13.2\;{\rm
au}}$ (Equations 22 and 25). We then use the collisional model to calculate
the edge profile at various times up to ${10\;{\rm Gyr}}$, the lifetime of the
solar-type star. The first snapshot is made shortly after $\tau_{\rm b}$, to
ensure there is sufficient dust in the system. Figure 13 shows the results;
the inner edge is largely unchanged for the first ${100\;{\rm Myr}}$, but from
${1\;{\rm Gyr}}$ collisions begin to make the edge flatter. We quantify this
flattening using the erf-powerlaw function; at 0.2, 10, 100, $10^{3}$ and
${10^{4}\;{\rm Myr}}$, the respective fitted values are ${r_{\rm i}=13.1}$,
13.0, 13.0, 13.0 and ${13.2\;{\rm au}}$, and ${\sigma_{{\rm i}}=0.034}$,
0.042, 0.043, 0.064 and 0.11. Note that the erf function is an increasingly
poor fit to the edge shape as $\tau_{\rm max}$ is approached, due to the
complex edge profile.

Figure 13: Collisional evolution of the planet-sculpted disc from Figure 2,
according to the Löhne et al. (2008) model (Section 4.2). We initialise the
disc with ${r_{\rm i}=13.2\;{\rm au}}$ and ${\sigma_{{\rm i}}=0.0361}$, and
assume ${m_{*}=1\;{\rm M_{\odot}}}$, ${x_{\rm m}}=1$ and $s_{\rm max}=50\;{\rm
km}$. This yields ${\tau_{\rm max}=1.1\;{\rm Gyr}}$ and ${\tau_{\rm
b}=0.14\;{\rm Myr}}$ at ${13.2\;{\rm au}}$. Solid lines show the edge at
various times; collisions flatten the edge, tending towards an $r^{2}$ profile
(dashed line). Note that our simulated collisions only deplete material; the
lines are scaled to the same value at ${13.2\;{\rm au}}$ for gradient
comparison, giving the impression that the surface density increases with time
at some radii, but in reality it decreases at each successive snapshot and
unscaled lines would not cross (the lines are scaled because the surface
density drops considerably with time, and a logarithmic plot would skew the
gradient comparison). The edge undergoes little evolution until the time
becomes comparable to $\tau_{\rm max}$; after that, the edge starts to
flatten.

Figure 13 demonstrates the importance of $\tau_{\rm max}$ in setting the
collisional evolution. Before $\tau_{\rm max}$, collisions have little effect
on the inner-edge profile. However, once $\tau_{\rm max}$ is reached,
collisions start to have a significant effect, making the edge flatter. This
collisional profile is expected to eventually tend towards $r^{2}$, the
profile for a broad collisional disc, as shown by the dashed line on Figure
13. However, in this example the profile is never reached within the stellar
lifetime; in fact, the inner edge assumes two slopes, with an erf-like profile
interior to $r_{\rm i}$ and a flatter, positive profile beyond this. This
means that, in this example, the effect of planet sculpting would still be
visible on the disc’s inner-edge even at very late times666Since dust mass,
and hence brightness, decreases with time, a real observation of a
collisionally evolved disc might require a high signal-to-noise ratio to
detect the deviation from an $r^{2}$ profile at $r_{\rm i}$..

Figure 13 was for one example setup. Figure 14 shows how the inner-edge
steepness collisionally evolves in a number of different setups, demonstrating
that this evolution is both qualitatively and quantitatively similar across a
broad parameter space. Starting with the setup from Figure 13 as a reference,
we vary the input parameters and re-run the collisional simulation multiple
times. We see that all of the setups have the same general collisional
evolution. First, between $\tau_{\rm b}$ and $\tau_{\rm max}$ the inner edge
becomes slightly flatter, with $\sigma_{{\rm i}}$ increasing by a factor of
${\sim 1.2}$. This remains constant, until the time approaches ${\tau_{\rm
max}}$; around this time the edge becomes much flatter, and tends towards the
$r^{2}$ collisional profile. The plot shows the importance of $\tau_{\rm max}$
in setting the collisional evolution of the sculpted-disc edge.

Figure 14: General collisional evolution of the inner edge of a planet-
sculpted debris disc. The plot shows the evolution of the inner-edge profile
over time, where time is expressed relative to $\tau_{\rm max}(r_{\rm i,0})$,
the collisional lifetime of the largest planetesimals at the initial inner
edge. The simulations start at time $\tau_{\rm b}$. The black line is the
reference simulation from Figure 13, and other lines show simulations that
differ from the reference setup by the parameters listed in the key. At early
times collisions make the edge slightly flatter (i.e. increase $\sigma_{{\rm
i}}$), before making it significantly flatter as the time approaches
$\tau_{\rm max}(r_{\rm i,0})$. The local peak around $\tau_{\rm max}$ should
be ignored; it arises because the erf-powerlaw function used to fit the
profile struggles with the complex edge shape at this time.

####  4.2.3 Implications of collisions for planet-sculpted debris discs

The previous sections showed that collisions have a very small effect on the
inner-edge profile, until the time $\tau_{\rm max}$. After this, a planet-
sculpted edge would become significantly flatter. Consideration of $\tau_{\rm
max}$, the stellar lifetime and the planet-sculpting timescale ${t_{\rm
diff}}$ are therefore vital for assessing the effect of collisions on the edge
profile.

We can use the above timescales to identify three possibilities for the
evolution of the debris-disc inner edge, depending on ${t_{\rm diff}}$ and
$\tau_{\rm max}$:

  * •

${t_{\rm diff}}\ll\tau_{\rm max}$: the planet sculpts the disc, producing a
sharp inner edge like our $n$-body simulations. This edge then gradually gets
flatter due to collisions, but does not change significantly until the time
approaches $\tau_{\rm max}$. Then, if left long enough, it may eventually
assume an $r^{2}$ profile.

  * •

${t_{\rm diff}}\sim\tau_{\rm max}$: unknown outcome. In this regime planetary
sculpting and collisions are equally important at all times, and our
simulations are insufficient to model this.

  * •

${t_{\rm diff}}\gg\tau_{\rm max}$: the disc collisionally evolves almost as if
no planet were present. The disc tends toward an $r^{2}$ profile.

So if ${t_{\rm diff}}$ and $\tau_{\rm max}$ are known, then the relative
importance of collisions to planetary sculpting can be ascertained, and the
evolution of the disc edge predicted. However, whilst ${t_{\rm diff}}$ is
straightforward, Equations 23 to 25 show that $\tau_{\rm max}$ is a
complicated function depending on many parameters, several of which are
interdependent and unknown.

We can gain deeper insight by rewriting $\tau_{\rm max}$, to make its
dependencies more explicit. Equations 23 to 25 give $\tau_{\rm max}$, $G$ and
$X_{\rm c}$ in their general forms; we can simplify these by substituting our
assumed form for the initial $\Sigma(r)$, as well as $Q_{\rm D}^{*}$, $f$,
$I/e$ etc. We also assume ${s_{\rm max}\gg s_{\rm min}}$. For the region
beyond the inner edge, this yields

$\tau_{\rm max}=\frac{649\;{\rm yr}}{{x_{\rm m}}}\;\frac{{q_{\rm
g}}-5/3}{2-{q_{\rm p}}}\left(\frac{r}{\;{\rm
au}}\right)^{3/2+\alpha}\left(\frac{s_{\rm max}}{\;{\rm
km}}\right)\left(\frac{m_{*}}{\;{\rm M_{\odot}}}\right)^{-1/2}\\\
\times\left(\frac{\rho}{\;{\rm g\;cm^{-3}}}\right)G^{-1}({q_{\rm g}},s_{\rm
max},r),$ (28)

where

$G({q_{\rm g}},s_{\rm max},r)=\left[X_{\rm c}(s_{\rm max},r)^{5-3{q_{\rm
g}}}-1\right]\\\ +2\frac{{q_{\rm g}}-5/3}{{q_{\rm g}}-4/3}\left[X_{\rm
c}(s_{\rm max},r)^{4-3{q_{\rm g}}}-1\right]\\\ +\frac{{q_{\rm g}}-5/3}{{q_{\rm
g}}-1}\left[X_{\rm c}(s_{\rm max},r)^{3-3{q_{\rm g}}}-1\right]$ (29)

and

$X_{\rm c}(s_{\rm max},r)=0.0138\;e^{-1/2}\left(\frac{s_{\rm max}}{\;{\rm
km}}\right)^{{b_{\rm g}}}\left(\frac{r}{\;{\rm
au}}\right)^{1/4}\left(\frac{m_{*}}{\;{\rm M_{\odot}}}\right)^{-1/4}.$ (30)

Expressing the equations in this form demonstrates that $\tau_{\rm max}$ is
inversely proportional to the unknown ${x_{\rm m}}$, which is related to the
initial mass of the debris disc in MMSN units. We use this to better visualise
$\tau_{\rm max}$.

Figure 15 shows $\tau_{\rm max}$ multiplied by ${x_{\rm m}}$, for a Solar-type
star and a disc with an initial surface-density profile of $r^{3/2}$. The
figure shows that, for reasonable debris parameters, $\tau_{\rm max}$ is
typically extremely long, unless $s_{\rm max}$ is very small or ${x_{\rm m}}$
very large. However, $s_{\rm max}$ is unlikely to be smaller than a few
kilometers, otherwise it would violate planetesimal-formation models and the
statistics of discs of various ages (Krivov & Wyatt, 2021, and refs therein).
A ${100\;{\rm au}}$ disc with debris eccentricity 0.01 would therefore have
${{x_{\rm m}}\tau_{\rm max}\sim 10^{10}\;{\rm yr}}$. Furthermore, ${{x_{\rm
m}}>10}$ is unlikely because such debris discs would have masses exceeding
those of solids in the preceding, protoplanetary stage (Krivov & Wyatt, 2021).
As a result, for many observed debris discs $\tau_{\rm max}$ may be longer
than the ${\sim 10^{7}}$ to ${10^{9}\;{\rm yr}}$ system age. In fact,
${\tau_{\rm max}}$ can even be infinite, which means that collisions never
disrupt the largest bodies; this occurs in the large regions of parameter
space where ${X_{\rm c}(s_{\rm max})\geq 1}$ (for which ${G\leq 0}$).
Rearranging Equation 24, $\tau_{\rm max}$ becomes infinite if

Figure 15: Plots of $\tau_{\rm max}$ (collisional lifetime of the largest
debris) multiplied by ${x_{\rm m}}$ (initial surface density relative to the
MMSN), for a Solar-type star and a disc with an initial surface-density
profile of $r^{3/2}$. The plots show that $\tau_{\rm max}$ strongly depends on
distance, $s_{\rm max}$ and debris eccentricity, and can sometimes be very
long. It is unlikely that ${{x_{\rm m}}>10}$, so $\tau_{\rm max}$ can be much
longer than the ${\sim 10^{7}}$ to ${10^{9}\;{\rm yr}}$ ages of debris-disc
systems. If the largest bodies lie above the solid line (Equation 31), then
they would never be destroyed and $\tau_{\rm max}$ becomes infinite. Note that
${x_{\rm m}}$ and $s_{\rm max}$ cannot be arbitrarily varied for an observed
debris disc; only certain combinations are allowed, since they set the
observed-dust mass at the system age (Equation 26). Increasing or decreasing
the star mass respectively raises or lowers the contours slightly, but does
not have a large effect.

$Q_{\rm D}^{*}(s_{\rm max})\geq\frac{f^{2}Gm_{*}}{2r},$ (31)

which yields the solid lines on Figure 15.

Conversely, calculating ${t_{\rm diff}}$ (Equation 2) shows that a Jupiter-
mass planet located within ${100\;{\rm au}}$ of a Solar-mass star would finish
sculpting a disc within ${<10\;{\rm Myr}}$. So if such planets were sculpting
debris discs, then they would likely do so long before collisions had any
effect on the inner edges. This would make the inner-edge profiles long lived,
and since even ${10\tau_{\rm max}}$ is insufficient to fully erase the
planet’s signature in the example on Figure 13, it is likely that a disc
sculpted by a massive planet would maintain a distinctive shape even if
collisions were fully accounted for. For such cases, the planet can be
constrained using the process in Section 3.2. For smaller, Neptune-mass
planets the sculpting timescale is 100 times longer, in which case both
collisions and planets could simultaneously sculpt the discs. Such systems may
not have reached their final configurations by the time they are observed,
meaning that the disc edges would be somewhere between the very steep values
expected for such planets, and the $r^{2}$ profile expected from collisions.
This is discussed further in Section 6.1.

##  5 Application to observed debris discs

Finally, we now compare the inner-edge slopes we predict from planetary
sculpting to the inner edges of observed debris discs. We consider seven ALMA-
resolved, extrasolar debris discs, each of which has a literature measurement
of its inner-edge steepness. The systems are listed in Table 1, and were
analysed by Lovell et al. (2021), Faramaz et al. (2021) and Imaz Blanco et al.
(2023).

Table 1: Constraints on the inner edges of several axisymmetric, extrasolar debris discs. The parameters $r_{\rm t}$, ${p_{\rm i}}$, ${p_{\rm o}}$ and $\eta$ are from literature fits to the surface-density profiles in the inner-edge regions (Equation 32), and $r_{\rm i}$ and $\sigma_{{\rm i}}$ are the equivalent parameters in our erf-powerlaw model for comparison with our simulations (Equation 33). Note that ${p_{\rm i}}$ and ${p_{\rm o}}$ are often 0.5 larger than the values provided in disc-imaging papers, because we consider the surface-density slope whilst those papers often give the intensity slope (e.g. Imaz Blanco et al. 2023). System | Name | Star mass /$\;{\rm M_{\odot}}$ | Age /${\;{\rm Myr}}$ | $r_{\rm t}\;/\;{\rm au}$ | ${p_{\rm i}}$ | ${p_{\rm o}}$ | $\eta$ | $r_{\rm i}\;/\;{\rm au}$ | $\sigma_{{\rm i}}$ | Refs.
---|---|---|---|---|---|---|---|---|---|---
HD 9672 | 49 Ceti | $1.98\pm 0.01$ | $45\pm 5$ | $131^{+13}_{-12}$ | $1.3^{+0.4}_{-0.3}$ | $-3.0^{+0.4}_{-0.5}$ | 2 | - | $0.6^{+0.1}_{-0.2}$ | 1, 3
HD 10647 | q1 Eri | $1.13^{+0.03}_{-0.04}$ | $1700\pm 600$ | $76\pm 1$ | $>5.6$ | $-2.1\pm 0.2$ | 2 | 70 | $<0.14$ | 1, 4
HD 92945 | V419 Hya | $0.87\pm 0.01$ | $250\pm 100$ | $54\pm 2$ | $8\pm 2$ | $-0.8^{+0.4}_{-0.6}$ | 2 | 52 | $0.11\pm 0.03$ | 1, 3
HD 107146 | - | $1.03^{+0.02}_{-0.04}$ | $150^{+100}_{-50}$ | $44\pm 2$ | $7.2^{+0.9}_{-0.7}$ | $-0.2^{+0.1}_{-0.2}$ | $2.8^{+1.2}_{-0.7}$ | 41 | ${0.11\pm 0.01}$ | 1, 3
HD 197481 | AU Mic | $0.59\pm 0.03$ | $24\pm 3$ | $36.4\pm 0.7$ | $1.4\pm 0.4$ | $-9\pm 1$ | 2 | - | $0.6\pm 0.3$ | 1, 3
HD 206893 | - | $1.32^{+0.07}_{-0.05}$ | $160\pm 20$ | $35^{+7}_{-10}$ | $>1.1$ | $-2.2\pm 0.2$ | 2 | - | $<0.76$ | 2, 3
HD 218396 | HR 8799 | $1.59\pm 0.02$ | $42^{+6}_{-4}$ | $180^{+10}_{-20}$ | $3.0^{+0.9}_{-0.5}$ | $-0.6^{+0.5}_{-0.3}$ | 2 | 160 | $0.27^{+0.04}_{-0.08}$ | 1, 5
References: Star masses and ages from (1) Pearce et al. (2022), (2) Hinkley et
al. (2023). ALMA surface-density fits from (3) Imaz Blanco et al. (2023),
(4) Lovell et al. (2021), (5) Model 2 in Faramaz et al. (2021).

The literature works used various parameterisations to quantify disc profiles,
which we must first convert to erf-powerlaw functions for comparison with our
simulations. The literature works generally fitted the inner edges using a
radial powerlaw, with some turnover occurring further out; this turnover could
be the disc outer edge, or a gap within the disc. For our analyses the nature
of the turnover is unimportant, but its location is needed to define the
‘inner’ region of the disc. Following Lovell et al. (2021), we quantify the
surface densities around the inner edges of the literature discs using a
double powerlaw:

$\Sigma(r)\propto\left[\left(\frac{r}{r_{\rm t}}\right)^{-\eta{p_{\rm
i}}}+\left(\frac{r}{r_{\rm t}}\right)^{-\eta{p_{\rm o}}}\right]^{-1/\eta},$
(32)

where ${p_{\rm i}}$ is the slope of the inner edge, $r_{\rm t}$ is the radial
location of the turnover, ${p_{\rm o}}$ is the slope exterior to the turnover,
and $\eta$ sets how sharp the turnover is. The assumed values for each disc
are listed in Table 1; for fits where $\eta$ was undefined or unconstrained,
we assumed a value of 2 as in Lovell et al. (2021).

We next convert these powerlaw fits into erf-powerlaw profiles for comparison
with our simulations. This is not strictly valid as we should really fit our
erf-powerlaw profile to the underlying data, but it should be sufficient for
this simple analysis. We consider the profile

$\Sigma(r)\propto\left[1-{\rm erf}\left(\frac{r_{\rm
i}-r}{\sqrt{2}\sigma_{{\rm i}}r_{\rm i}}\right)\right]\left(\frac{r}{r_{\rm
i}}\right)^{{p_{\rm o}}},$ (33)

which is equivalent our Equation 6 near the inner edge. We set ${p_{\rm o}}$
to the literature values describing the region just beyond the edge, noting
that ${p_{\rm o}}$ is typically negative. We then use Equation 41 to convert
the inner-edge powerlaw index ${p_{\rm i}}$ to an equivalent flatness
$\sigma_{{\rm i}}$, which approximates how the inner edge would appear if
fitted with an erf function instead. Finally, we deduce $r_{\rm i}$, which
characterises the location of the inner edge in the erf model, such that
Equations 32 and 33 peak at the same radial location. A comparison of the erf-
powerlaw and double-powerlaw profiles for an example system is shown on Figure
18, and the inferred values of $\sigma_{{\rm i}}$ and $r_{\rm i}$ are listed
for all systems in Table 1. In some cases $r_{\rm i}$ could not be fitted
because the region beyond the turnover was so steep that Equation 33 does not
turn over; however, this does not affect the edge flatness $\sigma_{{\rm i}}$,
which is the parameter of interest in this section.

Figure 16 shows the inferred $\sigma_{{\rm i}}$ values for the observed discs,
which exhibit a range of inner-edge slopes. At least two have relatively flat
slopes consistent with pure collisional evolution, which would yield powerlaw
indices of ${{p_{\rm i}}\approx 2}$ as described in Section 4 (equivalent to
${\sigma_{{\rm i}}\approx 0.4}$, shown as the dashed line on Figure 16).
Several have much steeper inner edges, which Imaz Blanco et al. (2023) argued
are indicative of planetary sculpting. However, whilst those edges are indeed
steeper than expected from collisions alone, they are still flatter than our
simulations predict for in-situ sculpting of low-eccentricity discs by planets
on circular orbits. The shaded region on Figure 16 shows the maximum
$\sigma_{{\rm i}}$ expected if a ${m_{\rm p}/m_{*}\leq 10^{-3}}$ planet
sculpts a disc with ${e_{\rm max,0}\leq 0.1}$; all of the observed inner edges
are flatter than this777Whilst a planet with mass ${m_{\rm p}/m_{*}\approx
0.07}$ could drive $\sigma_{{\rm i}}$ to required values of ${0.11}$ (Figure
6), we disfavour this possibility because it seems unlikely that planets with
such specific masses exist in two of our seven systems.. In Section 6.1 we
discuss the potential implications of these flatter edges, which could teach
us about architectures, dynamical processes and histories in the outer regions
of planetary systems.

Figure 16: Flatness of the inner edges of observed discs (Table 1). Blue
circles and orange triangles are $\sigma_{{\rm i}}$ values and upper limits
respectively. Several discs (49 Ceti, AU Mic and HD 206893) have flat inner-
edge slopes consistent with pure collisional evolution, as shown by the dashed
horizontal line. Other discs (q1 Eri, HD 92945 and HD 107146) have much
steeper profiles, but these are still flatter than expected from in-situ
sculpting of low-eccentricity discs by circular-orbit planets.

##  6 Discussion

In Section 6.1 we discuss the possible implications of observed discs having
flatter-than-expected edge profiles, and in Section 6.2 we compare the
analyses in this paper to previous studies.

###  6.1 Why are observed inner edges flatter than expected from sculpting
planets?

In Section 5 we showed that, whilst the edges of several ALMA-resolved discs
are steeper than expected from collisions alone, they are flatter than
expected from planetary sculpting of an initially low-eccentricity disc. Here
we discuss several possible reasons for this, and what these could teach us
about planetary systems.

####  6.1.1 Debris discs have high intrinsic-excitation levels

One possibility is that planetesimals in debris discs have higher intrinsic
excitation levels than is often assumed. For example, if debris at a disc
inner edge had ${e_{\rm i,rms}\approx 0.13}$, then that edge would have
${\sigma_{{\rm i}}\approx 0.11}$, which is similar to the steepest discs in
Table 1. This value of $e_{\rm i,rms}$ is higher than the maximum of 0.06 that
could be imparted by a non-migrating, circular-orbit planet (Figure 7),
raising the possibility that these eccentricities arise through non-planetary
processes. Examples of possible excitation mechanisms include self stirring by
the debris disc’s self gravity, or processes that occurred during system
formation. Marino (2021) used the outer edges of several observed discs to
infer the debris-excitation levels, and found these to be consistent with
${e_{\rm i,rms}\gtrsim 0.1}$; that sample included ${{\rm HD}\;92945}$ and
${{\rm HD}\;107146}$, which have the steepest edges in Table 1. It is
therefore plausible that in-situ planets sculpted the inner edges of some
observed discs, but that this debris has separately been excited by some non-
planetary processes.

One means to test this would be to measure scale heights at the inner edges of
debris discs. We find that MMRs of a coplanar planet excite debris
eccentricity but not inclination (Section 2.2), so if in situ, coplanar
planets sculpt debris-disc inner edges, then these edges would be more excited
radially than vertically.

####  6.1.2 Planets are present at the inner edges, but sculpting has not yet
finished

Planets take time to eject debris and fully sculpt a debris-disc inner edge.
The required time is characterised by 1 to 10 diffusion timescales (Section
2.3.3). Until this time, the disc’s inner edge could have almost any profile,
which would depend on its initial configuration. It is therefore possible that
planets reside just interior to debris discs, but the planets have not yet
finished sculpting the edges into steep profiles. This is particularly
plausible for young systems, or those with low-mass planets.

The observed discs in Table 1 with the flattest inner-edge profiles are also
the youngest, with ages less than ${100\;{\rm Myr}}$. These are 49 Ceti, AU
Mic and ${{\rm HR}\;8799}$. Equation 15 shows that, if planets with masses
below 0.1 to ${1\;{\rm M_{Jup}}}$ lie at the inner edges of those discs, then
these systems would be younger than one diffusion time and hence planetary
sculpting would not yet have finished. It is therefore plausible that inner-
edge profiles are flatter than expected not because sculpting planets are
absent, but because sculpting has not yet finished. If true, then care must be
taken when inferring planet properties from debris discs, because the
hypothesised planet-disc interactions may not have reached their final state.
This may be a prevalent problem because the discs most favourable for
observation are generally the youngest, for which the dust mass (and hence
brightness) is highest.

####  6.1.3 Discs are sculpted by migrating planets

In this paper we modelled planets on non-evolving orbits, but a migrating
planet would impose a very different morphology on the disc inner edge.
Specifically, a migrating planet causes MMR sweeping, where the nominal
locations of MMRs move across the debris disc and trap large numbers of bodies
in resonance. During this process, resonant debris is excited to increasingly
high eccentricities (e.g. Wyatt 2003; Reche et al. 2008; Friebe et al. 2022),
which would result in a flatter edge profile.

An example of this occurred in the Solar System. Neptune is located at the
inner edge of the Kuiper Belt, but the orbits of Kuiper-Belt objects (KBOs)
show that the planet migrated outwards in the past (Malhotra, 1993). As a
result, KBOs at the Kuiper Belt’s inner edge have relatively high
eccentricities, which make the Belt’s inner edge flatter than would be
expected from perturbations by a non-migrating Neptune. In an upcoming paper
(Morgner et al., in prep.), we take the observed KBOs and, using a method
similar to Vitense et al. (2010), de-bias these observations to estimate the
true KBO population distribution. We then fit the surface density of this de-
biased population with Equation 33, to get the steepness of the Kuiper-Belt
inner edge. This yields ${\sigma_{{\rm i}}=0.076}$, which is much flatter than
would be expected for a non-migrating Neptune; Neptune has ${m_{\rm
p}/m_{*}=5\times 10^{-5}}$, so the inner edge should have ${\sigma_{{\rm
i}}<0.02}$ if Neptune’s orbit never changed (Figure 6). The reason for this
difference is that Neptune excited KBOs through MMR sweeping as it migrated,
resulting in eccentricities at the Kuiper Belt’s inner edge of ${e_{\rm
i,rms}\sim 0.2}$, compared to just 0.01 expected from excitation by a non-
migrating Neptune (Equation 7).

It is possible that some extrasolar discs were also sculpted by migrating
planets, which made their edge profiles flatter than our in-situ model
predicts. Dust clumps that may be indicative of planetary migration are
inferred in debris discs (Lovell et al., 2021; Booth et al., 2023), and Pearce
et al. (2022) argue that migration may be required to relate the location of
debris-disc edges to system-formation theory. These results, combined with
debris-disc edges being flatter than expected for sculpting by non-migrating
planets, could imply that planetary migration is common in the outer regions
of debris-disc systems. Such migration could be caused by planetesimal
scattering (e.g. Friebe et al. 2022), planet-gas interactions in the
protoplanetary disc phase, or planet-planet scattering.

It may be possible to test exoplanet-migration scenarios in the near future.
JWST will search for planets in the outer regions of debris-disc systems, and
should detect many of the inferred planets if they exist (Pearce et al.,
2022). If a planet were found just interior to a debris disc, and that disc
had a flatter edge profile than would be expected from in-situ planetary
sculpting, then it could be evidence that the planet migrated outwards in the
past, exciting debris as it did so.

####  6.1.4 Collisions have flattened the planet-sculpted edges

Debris-disc edges could initially have been sculpted by planets, and these
edges could since have undergone significant collisional evolution. The
steepest edges of the Table 1 discs have ${\sigma_{{\rm i}}\approx 0.11}$,
which is at least 1.3 times flatter than expected from sculpting by non-
migrating planets (Figure 6). However, Figure 14 shows that collisions can
quickly increase $\sigma_{{\rm i}}$ by a factor of 1.2, long before the
largest bodies start colliding at $\tau_{\rm max}$. This could explain the
steepest inner edges in Table 1; in these cases, the system ages could be
longer than $\tau_{\rm b}$ but shorter than $\tau_{\rm max}$. After $\tau_{\rm
max}$, collisions would flatten the edge further, increasing $\sigma_{{\rm
i}}$ by a factor of 3 to 6 within $10\tau_{\rm max}$; this could potentially
reproduce some of the flatter edges in Table 1, if those systems are already
older than $\tau_{\rm max}$. Since some observed edges are too steep to be
explained by collisions alone, yet too flat to be explained by non-migrating
planets alone, this could suggest that both processes play a significant role
in some systems.

####  6.1.5 Planets are present at debris-disc inner edges, but the edges are
set by planet formation rather than scattering

A circular-orbit planet would eject debris that initially lies within 3 Hill
radii of its orbit (Figure 5). Therefore, if a newly formed system had a
planetesimal disc extending down to less than 3 Hill radii exterior to a
planet’s orbit, then that planet would eventually impart a sharp inner edge on
the disc. This is the scenario we modelled in our simulations. However, a real
planet would also have accreted material from around its orbit as it formed.
If the forming planet’s feeding zone were wider than 3 Hill radii, then much
of the debris within the scattering radius would already have been accreted by
the time the planet formed. Also, other processes in the protoplanetary disc
phase could further affect planetesimal profiles, such as the forming and
trapping of material in planet-induced pressure bumps near a planet.
Therefore, it is possible that a debris disc’s edge profile is set by
processes occurring during planet formation, rather than by pure planetary
scattering; in this case, its edge profile could differ from those that we
model (e.g. Eriksson et al. 2020).

####  6.1.6 Debris-disc edges are set by planetesimal formation alone

Planets are often invoked to explain the shapes and locations of debris-disc
inner edges, motivated in part by Neptune dominating the inner edge of the
Kuiper Belt. However, an alternative possibility is that planetesimals
naturally form in distinct radial zones; in this case debris-disc edges may be
completely unrelated to planets, but instead mark locations were planetesimal
formation transitioned from being inefficient to efficient. Such formation
zones are predicted by some streaming-instability models (Carrera et al.,
2017), as well as models involving snowlines (Ida & Guillot, 2016;
Schoonenberg & Ormel, 2017; Drążkowska & Alibert, 2017; Schoonenberg et al.,
2018; Izidoro et al., 2022; Morbidelli et al., 2022). This possibility can
soon be tested because, if sculpting planets exist, then many should be
detectable by JWST (Pearce et al., 2022); the absence of such detections would
potentially imply that other processes set disc-edge locations.

####  6.1.7 Debris discs do not have sharp cutoffs in semimajor axis

Debris-disc inner edges could be broad in semimajor-axes space, meaning that
the number density of debris gradually increases with semimajor axis. The
radial profile of such a disc would mimic that of a higher-eccentricity disc
with a sharper semimajor-axis distribution. Shallow semimajor-axis
distributions could be a natural outcome of system formation (Section 6.1.6),
arise through MMR sweeping during planet migration (Friebe et al., 2022), or
result from debris diffusion during self stirring (Ida & Makino, 1993). Since
the average eccentricities and inclinations would be related in a relaxed
debris disc, the disc vertical thickness would test this; Marino (2021) argued
that the outer region of the AU Mic disc is more extended radially than
vertically, which could imply low eccentricities and hence a slowly changing
semimajor-axis profile.

####  6.1.8 Debris discs are sculpted by planets with low but non-zero
eccentricities

Eccentric planets can excite debris to higher eccentricities than circular-
orbit planets, due to additional secular perturbations (eccentric planets are
considered in Appendix B). This would result in flatter disc inner edges.
Combining Equations 13 and 38, we can show that a Neptune-mass planet orbiting
a solar-type star with eccentricity 0.10 could induce an inner-edge steepness
of ${\sigma_{{\rm i}}\approx 0.11}$, like the sharpest observed edges. Such a
planet would make the inner edge elliptical with global eccentricity 0.052
(Equation 40), which is high enough to detect with modern observations (e.g.
Lovell et al. 2021). Therefore, it should be possible to rule out an eccentric
planet as being responsible for a flat disc edge if the disc’s global
eccentricity were robustly measured as negligible.

However, an intriguing possibility arises if the star position is poorly
constrained. In this scenario, a disc with a low global eccentricity could
masquerade as axisymmetric, owing to a poorly constrained stellar offset. Were
this the case, then an eccentric planet could flatten the disc inner edge,
whilst no global disc eccentricity would be detected. One way to rule out this
scenario would be to measure the azimuthal variation in the disc profile,
because a planet sculpting a truly eccentric disc would make the inner edge
flatter at pericentre than apocentre (Section B.2). Depending on the disc
width and the intrinsic debris excitation, the eccentric planet may also
flatten the outer edge, which could be observable (compare Figures 2 and 17).

####  6.1.9 Unresolved Trojans are present interior to the inner edges

Unresolved Trojans could be present at the disc inner edges. In our
simulations, we purposefully truncated the initial-disc edge at one Hill
radius exterior to the planet, to omit any Trojans and ensure that our fitted
${\sigma_{{\rm i}}}$ values describe the ‘true’ inner-edge slopes. However, if
sculpting planets are present at the inner edges of real discs, but there are
also unresolved Trojans co-orbiting with the planets, then the Trojans would
flatten the azimuthally averaged inner-edge profiles. The degree of flattening
would depend on the surface density of Trojans relative to the disc.

####  6.1.10 Dust transport is important

Dust-transport processes could be operating in addition to planetary
sculpting, which would make the disc edges more radially extended and their
profiles flatter. Poynting-Robertson (PR) drag can be significant for small
grains, causing discs to be radially extended in scattered light; whilst the
larger, ALMA-imaged grains we consider should be unaffected by PR drag, they
could be affected in a similar way by stellar winds. Winds would be
particularly important for discs around late-type stars, and would drag
material inwards from the planetesimal belt (e.g. Plavchan et al. 2005;
Reidemeister et al. 2011; Schüppler et al. 2015). This would flatten the
inner-edge profile. Similarly, CO gas is detected in many debris discs around
early-type stars, which could also cause radial drift through gas drag (e.g.
Krivov et al. 2009; Marino et al. 2020; Pearce et al. 2020). Again, this would
flatten the inner edges.

####  6.1.11 More complex disc-planet interactions occur in some debris discs

The inner edges of some specific systems in Table 1 could be excited by more
complex planet-disc interactions. There are gaps in the discs of both ${{\rm
HD}\;92945}$ and ${{\rm HD}\;107146}$, which are the discs with the steepest
inner edges. Friebe et al. (2022) showed that the simplest and most self-
consistent way to explain the morphology of ${{\rm HD}\;107146}$ is if a
planet has migrated across the gap; such migration would cause sweeping by
mean-motion resonances, which would excite debris at the disc inner edge
(their Figure 9, lower-right panel). In that model the planet would now lie
just exterior to the inner edge of the gap, and may have swept up a Trojan
population which could resemble the additional gap features in Imaz Blanco et
al. (2023). It is therefore possible that planets lying further out in more-
complex discs could excite debris at the inner edge, leading to flatter edge
profiles.

Alternatively, planets could excite inner-edge debris through secular
resonances rather than mean-motion resonances, which could lead to higher
debris eccentricities and thus flatter edges. A secular resonance occurs where
the apsidal precession rate of debris (due to the planet and disc self-
gravity) matches that of the planet (due to the disc or other planets), and
such resonances can drive up debris eccentricities and even open remote gaps
in broad discs (Pearce & Wyatt, 2015; Yelverton & Kennedy, 2018; Sefilian et
al., 2021, 2023). If a secular resonance were located near the disc inner
edge, then it could dominate over scattering and MMRs, and hence produce a
different edge location and profile to those in our simulations (e.g.
Smallwood 2023). Given the very large uncertainties on both debris-disc masses
and orbital architectures in the outer regions of systems, we cannot state for
certain whether this secular-resonance excitation occurs, but we can estimate
the disc masses that could lead to this scenario.

In a coplanar system with a single planet and an external, self-gravitating
disc, a secular resonance occurs at semimajor axis ${a_{\rm res}}$. The
location of ${a_{\rm res}}$ depends on the total disc mass $M_{\rm d}$
relative to that of the planet, as well as the planetary semimajor axis. A
general expression describing this relationship is derived in Sefilian et al.
(2021) (their Equation 19) which, assuming a power-law disc with surface
density $\propto a^{-1.5}$ spanning semimajor axes $a_{\rm in}$ to $a_{\rm
out}$, with ${a_{\rm p}\approx a_{\rm in}}$, reads as

$\frac{M_{\rm d}}{m_{\rm p}}\approx\frac{1.5}{|\psi_{1}({a_{\rm
res}})|}\left(\frac{a_{\rm out}}{a_{\rm in}}\right)^{1/2}\left(\frac{{a_{\rm
res}}}{a_{\rm in}}\right)^{-2.5},$ (34)

where $\psi_{1}(a)$ is a factor of order unity (see also Sefilian et al.
2023). A secular resonance will hence lie near the disc inner edge if ${M_{\rm
d}/m_{\rm p}\gtrsim 0.89\left(a_{\rm out}/a_{\rm in}\right)^{1/2}}$; here we
arbitrarily define ‘near’ as ${{a_{\rm res}}<1.2a_{\rm in}}$, where
${\psi_{1}(1.2a_{\rm in})\approx-1.07}$ (computed using Equation 18 in
Sefilian & Rafikov 2019, assuming ${H=0.05}$ for the disc scale height888For
comparison, defining the inner-edge region as ${<1.5a_{\rm in}}$ instead
yields ${M_{\rm d}/m_{\rm p}\gtrsim 0.75\left(a_{\rm out}/a_{\rm
in}\right)^{1/2}}$, with ${\psi_{1}(1.5a_{\rm in})\approx-0.73}$.). This
implies that, for disc-to-planet mass ratios of roughly unity or higher,
secular resonances could play a role in shaping the disc inner edge. Since
debris discs could have masses up to ${1000\;{\rm M_{\oplus}}}$ (Krivov &
Wyatt, 2021), this effect could be significant even for Jupiter-mass planets.
We will further investigate the effect of secular resonances on debris-disc
edges in a future work (Sefilian et al., in prep.).

###  6.2 Comparison to literature studies

Several literature studies also assessed the effect of a sculpting planet on
debris-disc inner edges. These used various different approaches, and yielded
various different results. In this Section we compare our paper to literature
works, specifically focussing on the $n$-body results (Section 6.2.1),
collisional modelling (Section 6.2.2) and the planet-inferring technique of
Pearce et al. (2022), which uses the locations of debris-disc inner edges but
not their steepnesses (Section 6.2.3).

####  6.2.1 Gravitational effects

The fact that MMRs of planets on circular orbits excite debris at disc inner
edges, and that more-massive planets cause higher excitation and hence flatter
edge profiles, was also identified in previous works. Mustill & Wyatt (2012)
used encounter maps to demonstrate this effect on planetesimals (their Figure
4), and showed that a similar effect arises in the prescription of Wisdom
(1980). Chiang et al. (2009) and Rodigas et al. (2014) consider small dust
released from a population of planet-sculpted planetesimals, and show that
this decrease in edge steepness with increasing planet mass is also expected
in scattered-light observations (their Figures 3 and 1b respectively).

Quillen (2006) and Quillen & Faber (2006) used MMR theory to predict that a
low-eccentricity planet should impose an eccentricity dispersion at the disc
inner edge that is proportional to $(m_{\rm p}/m_{*})^{3/7}$, which is
slightly steeper than our empirical $(m_{\rm p}/m_{*})^{1/3}$ for a circular-
orbit planet (Equation 7). However, we note that the zero-planet-eccentricity
simulations of Quillen & Faber (2006) appear more consistent with our flatter
index of 1/3 (filled pentagons on their Figure 2), and this index is also
predicted theoretically by Petrovich et al. (2013) (their Equations 19 and
34). We therefore argue that $(m_{\rm p}/m_{*})^{1/3}$ is a better estimate of
debris eccentricities at the inner edge of a disc sculpted by a circular, non-
migrating planet, provided that ${m_{\rm p}/m_{*}\lesssim 3.3\times 10^{-3}}$.
Above this mass ratio, we observe a turnover where increasing the planet mass
no longer excites debris to higher eccentricities; this is due to the inner
edge of the planet’s chaotic zone extending beyond the location of strong
MMRs, and was not observed by the above authors because their simulations did
not extend to such high mass ratios.

Tabeshian & Wiegert (2016, 2017) modelled the interaction between a planet and
a planetesimal disc, and showed that the MMRs of a circular-orbit planet can
impose asymmetric structure on an external disc (Figure 5 in Tabeshian &
Wiegert 2016). In particular, they show that MMRs can carve crescent-shaped
gaps in the disc, and that the inner-edge profile can vary with azimuth. We
find the same result; in the simulation on Figure 4, a crescent-shaped gap can
be seen aligned with the planet, and the inner-edge steepness differs slightly
on the left and right sides of the disc. However, we do not investigate this
asymmetry in more detail, because the difference is small enough that it is
unlikely to significantly affect any inferred-planet parameters. Also,
Tabeshian & Wiegert (2016) show that the degree of asymmetry scales with the
planet-to-star mass ratio, and this asymmetry is already small for the
comparatively high mass ratio on Figure 4; for smaller mass ratios, the
asymmetry would be even less pronounced.

####  6.2.2 Collisional effects

We find that collisions make the disc edge flatter, and move its
characteristic radius outwards (Figures 13 and 14). This is in qualitative
agreement with Figure 3 of Nesvold & Kuchner (2015), which shows that the same
effects manifest in the smack collision prescription; that prescription
simultaneously models the collisional and n-body dynamical evolution of debris
(Nesvold et al., 2013). Nesvold & Kuchner (2015) also find that it takes 10 to
100 collisional timescales for this process to occur (their Figure 4).
However, our collisional models, and our definitions of collisional
timescales, are fundamentally different.

We employ the analytical model of Löhne et al. (2008), where the most relevant
collisional timescale is $\tau_{\rm max}$, the collisional lifetime of the
largest bodies. This predominantly depends on disc location, debris
eccentricity, the disc’s initial mass and the size of the largest bodies. The
latter two are unknown, but for physically plausible values, a disc with inner
edge at ${50\;{\rm au}}$ would have $\tau_{\rm max}$ of at least ${\sim
10^{7}\;{\rm yr}}$, and potentially much longer (Figure 15). This means that
$\tau_{\rm max}$ should be at least comparable to the ages of debris-disc
stars (Table 1), and may be much longer; based on this, we argue that
collisions should not significantly reduce the steepness of planet-sculpted
edges in many observed discs.

This conclusion differs from Nesvold & Kuchner (2015), who suggest that
collisions typically would have a significant effect on the observed edges of
planet-sculpted discs. They model the collisional cascade in detail and self-
consistently with the dynamical evolution, unlike our simulations. However,
the difference between our conclusions is unlikely to be due to our separation
of dynamics and collisions, but rather differences between our collisional
assumptions and those of Nesvold & Kuchner (2015).

There are two main reasons for this difference. First, the maximum grain
radius in Nesvold & Kuchner (2015) is ${s_{\rm max}=10\;{\rm cm}}$ (their
Table 1); these pebbles are depleted much faster than the ${\gtrsim 1\;{\rm
km}}$ planetesimals in our model, so their collisional erosion occurs on much
shorter timescales. Second, the amount of dust in the Nesvold & Kuchner (2015)
models is higher because they assume a greater minimum grain size;
extrapolating from their ${s_{\rm min}=1\;{\rm mm}}$ to the typical ${s_{\rm
min}=1}$ to ${10\;{\rm\mu m}}$ expected from radiation pressure blowout
corresponds to vertical optical depths that are higher than ours by factors of
10 to 30, assuming a size distribution index of 3.5.

Accounting for the respective differences in $s_{\rm min}$ and $s_{\rm max}$
each increase the ratio between the collision timescales of smallest dust
(their Eq. 3) and that of the biggest objects (which govern the long-term
evolution). Nesvold & Kuchner (2015) define smallest-grain collisional
timescales of $10^{3}$ to ${10^{6}\;{\rm yr}}$ for a sample of observed discs
(their Section 4); instead of the factors of 10 to 100 in their Figure 4,
planetesimals would evolve collisionally on timescales that are longer by
several orders of magnitude, in line with our results. So this difference in
timescales appears to be the cause of our different conclusions regarding the
relative effects of planets and collisions on the disc inner edge.

Nonetheless, the Nesvold & Kuchner (2015) model captures additional physics
that is omitted in our simpler prescription. In particular, our model would
not capture the expected collisional enhancement in resonant populations
(Stark & Kuchner, 2009); this means we probably underestimate collision rates
at the MMR-dominated inner edges, and therefore overestimate collisional
timescales. This effect could be significant in some cases, because the inner-
edge eccentricities in our low-eccentricity planet simulations are determined
by MMRs. Detailed, self-consistent modelling is therefore required in future,
to properly assess the balance between collisions and planet-debris
interactions at the inner edges of debris discs.

####  6.2.3 Comparison to Pearce et al. (2022) planet predictions

Pearce et al. (2022) give a general model to infer the minimum mass, maximum
semimajor axis and minimum eccentricity of a perturbing planet based on the
shape and location of a debris disc’s inner edge. Those predictions are based
on Pearce & Wyatt (2014), and do not consider the steepness of the disc edge;
as a result they have a degeneracy in planet mass and semimajor axis, because
they cannot distinguish a low-mass planet near the disc from a high-mass
planet away from the disc. However, in this paper we make predictions that
also include information about the edge steepness; this can break the
degeneracy, because edge steepness can be directly related to planet mass (at
least in the collisionless regime).

Section 3 describes our simple analytic model to predict planet properties
from an observed-disc edge. These predictions are comparable to those of
Pearce et al. (2022), with three important differences. First, we have an
extra dependence on the edge profile, and hence the debris eccentricity
$e_{\rm i,rms}$; our predicted inner-edge location (Equation 11) is increased
by a factor of ${(1+2e_{\rm i,rms}/\sqrt{3})}$ relative to that in Pearce et
al. (2022).

Second, Pearce & Wyatt (2014) use a stricter definition of the sculpting
timescale; they assume sculpting takes 10 diffusion times, which is the time
it takes an eccentric planet to eject ${95\;{\rm per\;cent}}$ of unstable
material. However, we find that an edge sculpted by a circular-orbit planet
with mass ${m_{\rm p}/m_{*}<10^{-2}}$ takes just one diffusion timescale to
assume roughly its final configuration (Equation 15). This is because it takes
a planet exponentially longer to clear increasingly large fractions of
unstable debris; a circular-orbit planet removes ${\sim 70\;{\rm per\;cent}}$
of unstable debris within 1 diffusion time, but takes 10 diffusion times to
remove ${\sim 90\;{\rm per\;cent}}$ (Costa, Pearce & Krivov, submitted). Hence
our minimum-allowed planet mass can be a factor of ${\sqrt{10}\approx 3}$
times smaller than Pearce et al. (2022).

Finally, we take the width of the chaotic zone around the planet’s orbit to be
3 Hill radii, whilst Pearce et al. (2022) used 5. The reason for this is that
Pearce & Wyatt (2014) studied eccentric planets, for which the innermost
stable semimajor axis is 5 eccentric Hill radii exterior to planet apocentre
(their Figure 9), whilst for circular-orbit planets a value of 3 provides a
better fit (Gladman, 1993; Ida et al., 2000; Kirsh et al., 2009; Malhotra et
al., 2021; Friebe et al., 2022). It is unclear at what planet eccentricities
the transition from 3 to 5 Hill radii occurs, though Figure 17 shows that a
planet eccentricity of 0.2 is already high enough for 5 Hill radii to provide
a better estimate. This is consistent with Regály et al. (2018), who find the
transition occurs at planet eccentricities somewhere between 0 and 0.2 (their
Figure 5c).

The result of these three differences is that any planets predicted using our
model will differ slightly from the Pearce et al. (2022) model. However, the
mass predictions should typically be within an order of magnitude or so of
each other. For the example on Figure 11, the lower bound on planet mass
predicted by the Pearce et al. (2022) model is 10 times smaller than the mass
predicted from our model; the two predictions still lie in a similar region of
parameters space, along Line 1 on that figure.

###  6.3 Application to narrow discs

We only considered broad discs in our simulations, to isolate the inner edge
from the outer edge. This ensured that our inner-edge profiles were dominated
by planet-disc interactions, rather than the initial profile of the disc.
However, this means that our quantitative results may not hold for very narrow
discs, where the outer- and inner-edge profiles could overlap.

In Section 2.4 we argued that the characteristic inner-edge width is several
times ${r_{\rm i}\sigma_{{\rm i}}}$. Hence our results should hold if the disc
is sufficiently wide for the edges to be well separated, i.e.

$r_{\rm o}-r_{\rm i}\gtrsim r_{\rm i}\sigma_{{\rm i}}+r_{\rm o}\sigma_{{\rm
o}}.$ (35)

Equation 13 shows that $\sigma_{{\rm i}}$ is proportional to the eccentricity
of inner-edge debris, and a similar relation holds for the outer edge.
Therefore, if the eccentricity distribution were constant across the disc,
then Equation 35 implies that our results should hold provided the disc
fractional width is larger than the rms debris eccentricity; this is in
agreement with Marino (2021). Our full criterion (Equation 35) is slightly
more complicated, because the debris-eccentricity level would vary across a
planet-sculpted disc.

###  6.4 Range of explored mass ratios

We explored planet-to-star-mass ratios ranging from ${3\times 10^{-5}}$ to
$10^{-1}$. We did not consider lower-mass planets, because in many cases the
dynamical timescales would then become comparable to the stellar lifetimes.
For example, an Earth-mass planet at ${10\;{\rm au}}$ from a Solar-type star
(${m_{\rm p}/m_{*}=3\times 10^{-6}}$) would have a diffusion timescale of
${40\;{\rm Gyr}}$, which would exceed the stellar lifetime of ${\sim 10\;{\rm
Gyr}}$. Hence our results would not hold for very low mass ratios, because the
interaction we study would not have time to occur.

Our upper limit of ${m_{\rm p}/m_{*}=10^{-1}}$ represents the upper end of
what could realistically be called a ‘planet’. For example, this could be a
${10\;{\rm M_{Jup}}}$ super-Jupiter orbiting a ${0.1\;{\rm M_{\odot}}}$ M6
dwarf. In many cases, such high-mass companions would be detectable via
imaging, radial velocity or astrometry, particularly in younger systems (e.g.
Carter et al. 2023).

Were a companion detected, then our results could be used to infer its
evolutionary history. For example, if a debris disc in the system had an edge
steepness consistent with sculpting by the observed companion, but the
companion were detected far interior to the disc, then this could be evidence
that the detected companion sculpted the disc historically but has since
migrated inward. This scenario has been suggested as a way to explain non-
detections of planets at the inner edges of debris discs (Pearce et al.,
2022).

##  7 Conclusions

We perform a dynamical investigation into the effect of planets on the
profiles of debris-disc inner edges. We consider both planet-disc interactions
and debris collisions, and explore the interaction across a broad parameter
space. We quantify our simulated surface-density profiles using an erf
function, as in Rafikov (2023), for direct comparison with ALMA-resolved inner
edges. Our main conclusions are as follows:

  1. 1.

For a non-migrating, circular-orbit planet, in the case where collisions are
negligible, the steepness of the disc inner edge is set by the planet-to-star
mass ratio and the initial-disc excitation level. Lower-mass planets lead to
steeper inner edges, with the edge width proportional to ${(m_{\rm
p}/m_{*})^{1/3}}$ (Equations 7, 8 and 13).

  2. 2.

There is a maximum eccentricity that a planet on a circular, non-migrating
orbit can generally impart on the population of planetesimals at a debris-disc
inner edge, which is a root-mean-squared eccentricity of 0.06 (Equation 7).

  3. 3.

Considering the steepness of a debris disc’s inner edge when inferring the
properties of unseen sculpting planets can break the degeneracy between planet
mass and semimajor axis. We provide a step-by-step method for inferring such
planets from disc inner edges in Section 3.2.1.

  4. 4.

Eccentric planets make inner-edge profiles flatter, and introduce azimuthal
asymmetries (Appendix B).

  5. 5.

The inner edge of a purely planet-sculpted debris disc is much steeper than
that of a purely collisional debris disc.

  6. 6.

Collisions flatten the profile of a planet-sculpted inner edge. The effect of
collisions is small before the largest bodies start to collide; after this
time, collisions start to significantly flatten the inner edge.

  7. 7.

In most cases, collisions would never fully erase the signature of a sculpting
planet on a debris-disc inner edge.

  8. 8.

Whilst the inner edges of many ALMA-resolved debris discs are too steep to be
caused by collisions alone, they are also too flat to arise through pure
sculpting by non-migrating, circular-orbit planets, unless sculpting is still
ongoing. Possible implications of this for the outer regions of planetary
systems are discussed in Section 6.1.

## Acknowledgements

We thank the anonymous referee, whose comments improved the paper. We are
grateful for exploratory work by Zuzanna Jonczyk as a summer intern at the
Institute of Astronomy, University of Cambridge; her work helped identify some
of the challenges in connecting inner-edge profiles to sculpting planets. We
also thank Brenda Matthews for useful discussions. TDP, AVK, TL and MB were
supported by Deutsche Forschungsgemeinschaft (DFG) grants Kr 2164/13-2, Kr
2164/14-2, Kr 2164/15-2 and Lo 1715/2-2. TDP is also supported by a Warwick
Prize Fellowship, made possible by a generous philanthropic donation. MB also
received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No. 951815 (AtLAST). AAS
acknowledges support by the Alexander von Humboldt Foundation through a
Humboldt Research Fellowship for postdoctoral researchers. MRJ acknowledges
support from the European Union’s Horizon Europe programme underthe Marie
Sklodowska-Curie grant agreement No. 101064124, and funding provided by the
Institute of Physics Belgrade, through the grant by the Ministry of Science,
Technological Development, and Innovations of the Republic of Serbia. SM is
supported by a Royal Society University Research Fellowship (URF-R1-221669).

## Data Availability

The data underlying this article will be shared upon reasonable request to the
corresponding author.

## References

  * Booth et al. (2023) Booth M., et al., 2023, MNRAS, 521, 6180
  * Carrera et al. (2017) Carrera D., Gorti U., Johansen A., Davies M. B., 2017, ApJ, 839, 16
  * Carter et al. (2023) Carter A. L., et al., 2023, ApJ, 951, L20
  * Chiang et al. (2009) Chiang E., Kite E., Kalas P., Graham J. R., Clampin M., 2009, ApJ, 693, 734
  * Drążkowska & Alibert (2017) Drążkowska J., Alibert Y., 2017, A&A, 608, A92
  * Eriksson et al. (2020) Eriksson L. E. J., Johansen A., Liu B., 2020, A&A, 635, A110
  * Esposito et al. (2020) Esposito T. M., et al., 2020, AJ, 160, 24
  * Faramaz et al. (2014) Faramaz V., et al., 2014, A&A, 563, A72
  * Faramaz et al. (2019) Faramaz V., et al., 2019, AJ, 158, 162
  * Faramaz et al. (2021) Faramaz V., et al., 2021, AJ, 161, 271
  * Friebe et al. (2022) Friebe M. F., Pearce T. D., Löhne T., 2022, MNRAS, 512, 4441
  * Geiler & Krivov (2017) Geiler F., Krivov A. V., 2017, MNRAS, 468, 959
  * Gladman (1993) Gladman B., 1993, Icarus, 106, 247
  * Hayashi (1981) Hayashi C., 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Hinkley et al. (2023) Hinkley S., et al., 2023, A&A, 671, L5
  * Hughes et al. (2018) Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541
  * Ida & Guillot (2016) Ida S., Guillot T., 2016, A&A, 596, L3
  * Ida & Makino (1993) Ida S., Makino J., 1993, Icarus, 106, 210
  * Ida et al. (2000) Ida S., Bryden G., Lin D. N. C., Tanaka H., 2000, ApJ, 534, 428
  * Imaz Blanco et al. (2023) Imaz Blanco A., et al., 2023, MNRAS, 522, 6150
  * Izidoro et al. (2022) Izidoro A., Dasgupta R., Raymond S. N., Deienno R., Bitsch B., Isella A., 2022, Nature Astronomy, 6, 357
  * Kennedy & Wyatt (2010) Kennedy G. M., Wyatt M. C., 2010, MNRAS, 405, 1253
  * Kirsh et al. (2009) Kirsh D. R., Duncan M., Brasser R., Levison H. F., 2009, Icarus, 199, 197
  * Krivov (2010) Krivov A. V., 2010, RAA, 10, 383
  * Krivov & Wyatt (2021) Krivov A. V., Wyatt M. C., 2021, MNRAS, 500, 718
  * Krivov et al. (2009) Krivov A. V., Herrmann F., Brandeker A., Thébault P., 2009, A&A, 507, 1503
  * Krivov et al. (2018) Krivov A. V., Ide A., Löhne T., Johansen A., Blum J., 2018, MNRAS, 474, 2564
  * Löhne et al. (2008) Löhne T., Krivov A. V., Rodmann J., 2008, ApJ, 673, 1123
  * Löhne et al. (2012) Löhne T., et al., 2012, A&A, 537, A110
  * Lovell et al. (2021) Lovell J. B., et al., 2021, MNRAS, 506, 1978
  * Malhotra (1993) Malhotra R., 1993, Nature, 365, 819
  * Malhotra et al. (2021) Malhotra R., Castro-Cisneros J. D., Fitzgerald J., Rosengren A. J., Ross S., Todorovic N., Wu D., 2021, in AAS/Division of Dynamical Astronomy Meeting. p. 305.04
  * Marino (2021) Marino S., 2021, MNRAS, 503, 5100
  * Marino et al. (2020) Marino S., Flock M., Henning T., Kral Q., Matrà L., Wyatt M. C., 2020, MNRAS, 492, 4409
  * Morbidelli et al. (2022) Morbidelli A., Baillié K., Batygin K., Charnoz S., Guillot T., Rubie D. C., Kleine T., 2022, Nature Astronomy, 6, 72
  * Morrison & Malhotra (2015) Morrison S., Malhotra R., 2015, ApJ, 799, 41
  * Mouillet et al. (1997) Mouillet D., Larwood J. D., Papaloizou J. C. B., Lagrange A. M., 1997, MNRAS, 292, 896
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge University Press, Cambridge, UK
  * Mustill & Wyatt (2012) Mustill A. J., Wyatt M. C., 2012, MNRAS, 419, 3074