<figure><img src="image_1.png" /></figure>

# Debris Disks in Multiplanet Systems: Are Our Inferences Compromised by Unseen Planets?

Jiayin Dong<sup>1,2 Ⓓ</sup>, Rebekah I. Dawson<sup>1,2 Ⓓ</sup>, Andrew Shannon<sup>1,2,3 Ⓓ</sup>, and Sarah Morrison<sup>1,2,4 Ⓓ</sup>

Received 2019 September 20; revised 2019 December 19; accepted 2019 December 23; published 2020 January 24

## Abstract

Resolved debris disk features (e.g., warps, offsets, edges and gaps, azimuthal asymmetries, radially thickened rings, scale heights) contain valuable information about the underlying planetary systems, such as the posited planet's mass, semimajor axis, and other orbital parameters. Most existing models assume a single planet is sculpting the disk feature, but recent observations of mature planetary systems (e.g., by radial velocity surveys or *Kepler*) have revealed that many planets reside in multiplanet systems. Here we investigate if/how planet properties inferred from single-planet models are compromised when multiple planets reside in the system. For each disk feature, we build a two-planet model that includes a planet b with fixed parameters and a planet c with a full range of possible parameters. We investigate these two-planet systems and summarize the configurations for which assuming a single planet (i.e., planet b) leads to significantly flawed inferences of that planet's properties. We find that although disk features are usually primarily dominated by a single planet, when using single-planet models, we are at risk of misinterpreting planet properties by orders of magnitude in extreme cases. Specifically, we are at high risk of misinterpreting planet properties from disk warps; at moderate risk from disk edges and gaps, radially thickened rings, and scale height features; and at low risk from host star-disk center offsets and azimuthal asymmetries. We summarize situations where we can infer the need to use a multiplanet model instead of a single-planet one from disk morphology dissimilarities.

*Unified Astronomy Thesaurus concepts: Debris disks (363); Circumstellar matter (241)*

## 1. Introduction

Thousands of exoplanets have been discovered over the past few decades, including many by the high-sensitivity *Kepler* Mission (Borucki et al. 2010) and high-precision radial velocity (RV) detection (e.g., Cumming et al. 2008; Howard et al. 2010; Mayor et al. 2011). Exoplanets observed with a wide variety of orbital and physical properties (e.g., super-Earths, sub-Neptunes, hot/warm Jupiters; Winn & Fabrycky 2015 and references therein) require us to reassess the classical theory of planetary system formation and evolution built on the solar system prototype. Most known exoplanets reside in systems that are a few gigayears old: they are the end states of a system's early formation and evolution. Planetary system formation theories (e.g., Chiang et al. 2009; Schlichting 2014; Dawson & Johnson 2018, and references therein) have been established based on these end states to explain the origins of planetary systems. One of the most direct ways to test these theories and answer the fundamental questions of planetary system formation and evolution (e.g., where do planets form and how do they migrate?) is to characterize young planetary systems during their early states of formation and evolution (Su et al. 2019).

It is challenging to apply the most prolific exoplanet-detection techniques to young planetary systems. The strong stellar activity noise (e.g., flares, flicker, and jitter; Bastien et al. 2014) of young host stars can easily suppress or mimic exoplanet signals in transit and RV observations. Although direct imaging detection (e.g., Gemini Planet Imager, or GPI; Spectro-Polarimetric High-contrast Exoplanet REsearch instrument, or SPHERE) is most feasible for young planetary

systems, it is challenging as well because the high planet-star flux contrast results in low sensitivity to low-mass planets. However, we can still probe young planetary systems indirectly through their debris disks. A debris disk, a circumstellar disk of dust and planetesimals orbiting its host star, has been detected around roughly 17%-20% solar-type stars (Hughes et al. 2018; Sibthorpe et al. 2018, and references therein). The dust in a debris disk is continuously produced through collisions of planetesimals that are gravitationally redistributed by the posited underlying planets (Wyatt 2008). By investigating debris disk features (here we focus on warps, offsets, edges and gaps, azimuthal asymmetries, radially thickened rings, scale heights) induced by planet-disk interactions, we can characterize the corresponding planetary system properties, such as the planet mass, its orbital parameters, and the timescale expected for a disk to evolve the observed features (see Hughes et al. 2018 and references therein). Debris disk features have been well studied in previous literature using single-planet models. As shown in Figure 1, all of the following disk features may be interpreted as perturbations from a single planet.<fnref n="5" />

1. Warps: a warped disk has an inner disk inclined with respect to the outer disk. The disk needs to almost be edge-on so that the warp will be observed. A warp has often been interpreted as the signpost of a planet on an inclined orbit (e.g., β-Pictoris; Mouillet et al. 1997; Augereau et al. 2001).

<footnote n="1">¹Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA; jdong@psu.edu</footnote>

<footnote n="2">²Center for Exoplanets & Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="3">³LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="4">⁴Department of Physics, Astronomy & Materials Science, Missouri State University, Springfield, MO 65897, USA</footnote>

<footnote n="5">⁵Figure 1 simulation setup can be found in Appendix A.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** A gallery of debris disk features sculpted by a single planet simulated with the *N*-body integrator REBOUND (Rein & Spiegel 2015; Rein & Tamayo 2015) and REBOUNDx (Tamayo et al. 2019). Planets and their orbits are labeled in blue as dots and lines. Planetesimals (i.e., simulated as massless test particles) are labeled in gray. Planet locations are only indicated in systems where disk morphologies depend on their instantaneous locations. The left, middle, and right columns correspond to secular, resonant, and synodic features.</figcaption></figure>

2. Spiral arms: the spiral arm feature refers to a disk with two symmetric spiral arms (e.g., HD 141569;<fnref n="6" /> Biller et al. 2015; Konishi et al. 2016). The feature can be explained as the projected view of a warped disk and thus serves as the signpost of an inclined planet.

3. Offsets: an offset disk has an evident geometric center offset from its central star. This feature has been interpreted as the signpost of a planet on an eccentric orbit (e.g., Fomalhaut, HR 4796A; Telesco et al. 2000; Kalas et al. 2005; MacGregor et al. 2017). Lee & Chiang (2016) demonstrated that an offset ring could present various morphologies (e.g., rings, needles, moths) in scattered-light images caused by different viewing angles.

4. Edges and gaps: the edge-and-gap feature refers to a disk with gaps and sharp edges (e.g., HR 8799; Su et al. 2009; Matthews et al. 2014; Wilner et al. 2018). A void of planetesimals across a certain region of a disk indicates a sculptor with significant mass (i.e., a planet clearing its chaotic zone; Wisdom 1980).

5. Azimuthal asymmetries: clumpy features observed in a disk (e.g., β-Pictoris; Dent et al. 2014) may indicate the existence of a migrating planet that captured and redistributed planetesimals to certain longitudes (Wyatt 2003; Reche et al. 2008; Mustill & Wyatt 2011). Shannon et al. (2015) demonstrated that small dust grains migrating under Poynting–Robertson drag may also get trapped in resonance with planets and produce clumpy structures in scattered-light images.

6. Radially thickened rings: a radially thickened ring (e.g., Fomalhaut, HR 4796A) has often been interpreted as a planet close to the inner edge of a debris ring stirring up nearby material (Chiang et al. 2009; Rodigas et al. 2014).

7. Scale heights: a substantial disk scale height (e.g., AU Microscopii, or AU Mic; Krist et al. 2005; MacGregor et al. 2013; Daley et al. 2019) has been interpreted as a perturbation of a large stirring body in the planetary system (Quillen et al. 2007; Thébault & Augereau 2007).

These resolved debris disk features contain valuable information about the underlying planetary systems, such as the sculpting planet's mass and orbital parameters. Most existing models assume that a single planet is sculpting the disk feature. However, recent observations of mature planetary systems (e.g.,

<footnote n="6">⁶ HD 141569 could be a “hybrid” disk between the stage of protoplanetary and debris disks (Miley et al. 2018). If so, the spiral arms might be a hydrodynamic effect.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The evolution of a warp over time, demonstrated with the planetesimal’s instantaneous inclination vs. semimajor axis (left panel) and the projected image of the warp along the line of nodes of the planet (right panel). The warp is produced by a 10 Jupiter-mass planet with a semimajor axis of 10 au and an inclination of 5°. The warp production time (i.e., the time taken to produce the warp) is labeled in each panel. The planetesimals’ inclinations oscillate between zero and twice the planet’s inclination at different nodal precession rates. The warp at a certain warp production time is located at the semimajor axis where planetesimals are just reaching their peak inclination in their first secular oscillation cycle, labeled with the blue dashed line or arrow.</figcaption></figure>

by RV surveys and *Kepler*) have revealed that many planets reside in multiplanet systems (e.g., see Winn & Fabrycky 2015 and references therein). In this paper, we investigate if/how planet properties inferred from single-planet models are compromised when multiple planets reside in the system. We categorize disk features listed in Figure 1 into three categories: secular, resonant, and synodic, according to their physical processes and timescales. Warps, spiral arms, and offsets are secular features driven by the interchanging angular momentum of planets and planetesimals over a long timescale (Section 2). Edges, gaps, and azimuthal asymmetries are resonant features that result from planets clearing nearby regions due to overlapping resonances or capturing planetesimals into orbital resonances (Section 3). Radially thickened rings and scale heights are synodic features that occur when planets stir up nearby material either radially or vertically over a short timescale (Section 4). For each disk feature, we build two-planet models that include a planet b with fixed parameters and a planet c with a full range of possible parameters. We characterize planet properties from the disk features in these two-planet systems and summarize the configurations for which assuming a single planet (i.e., planet b) leads to significantly flawed inferences of that planet’s properties. Our model considers planetesimals as massless, collisionless, and radiation-free test particles and studies how the gravitational influence of planets affects their distributions. Detailed modeling of other physical processes, such as collisions between planetesimals, is deferred to future work. We discuss secular features (warps, offsets), resonant features (edges and gaps, azimuthal asymmetries), and synodic features (radially thickened rings, scale heights) in Section 2–4, respectively. We summarize our main results and list situations where we can infer the need to use a multiplanet model instead of a single-planet one in Section 5.

## 2. Secular Features in Two-planet Systems

Secular features are driven by the interchanging angular momentum of planets and planetesimals over a long timescale. Through secular interactions, a planet on an inclined or eccentric orbit modifies the disk morphology by changing the planetesimals’ inclinations or eccentricities. Disk feature properties (e.g., a warped disk’s warp location and inclination) can constrain the underlying posited planet’s mass, semimajor axis, and other orbital parameters. In Section 2, we investigate

if/how our inferences of these planet properties are compromised using a single-planet model when multiple planets reside in the system. To feasibly explore a wide range of two-planet configurations each with thousands of planetesimals, we take advantage of the Laplace–Lagrange secular theory, a first-order approximation in *e* and *i* of the secular terms in the disturbing function. The Laplace–Lagrange approximation is appropriate for planets with low inclinations ($i \lesssim 20°$) and eccentricities ($e \lesssim 0.3$; Murray & Dermott 1999, Chapter 7; see Appendix B in a nutshell). The analytical approximation, compared to *N*-body simulations, greatly reduces the computation time and allows us to explore hundreds of two-planet configurations in a reasonable time.

### 2.1. Warps

A warped disk, an inner disk inclined with respect to the outer disk, has been interpreted as evidence for a sculpting planet on an inclined orbit. The vertical warp gradually moves outward in the disk with a rate depending on the sculpting planet’s mass and semimajor axis. Figure 2 demonstrates how a warp moves outward in a disk as a function of the warp production time. The warp production time is the time to produce the observed warp (i.e., calculated from the system’s age minus the planet formation time). We approximate the warp production time by assuming the timescale for the planet to form and get onto its current orbit is either short compared to the warp production time (i.e., the planet instantaneously appears with its present-day properties) and/or planet formation happened during gas-disk stage when planetesimals were protected from the planet by gas. With an estimation of the warp production time $\tau$ and the warp location $a_{\text{warp}}$, the planet’s mass and semimajor axis ($m_b$ and $a_b$) can be constrained, although with degeneracy:

$$ \frac{\pi}{\tau} = n_{\text{warp}} \frac{1}{4} \frac{m_b}{m_*} \alpha_b \bar{\alpha}_b b_{3/2}^{(1)}(\alpha_b), \quad (1) $$

where $n_{\text{warp}}$ is the orbital frequency of a planetesimal at $a_{\text{warp}}$, the semimajor axis of the warp, $m_*$ is the host star mass, $b_{3/2}^{(1)}$ is the Laplace coefficient, and $\alpha_b$ and $\bar{\alpha}_b$ follow the formula:

$$ \alpha_b, \bar{\alpha}_b = \begin{cases} a_b/a_{\text{warp}}, & \text{if } a_b < a_{\text{warp}} \\ a_{\text{warp}}/a_b, & \text{if } a_b > a_{\text{warp}}. \end{cases} \quad (2) $$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Warp location deviation in percentage for two-planet configurations vs. one-dominant-planet configurations, where the dominant planet is defined as the planet that primarily sculpted the warp. The warp location deviation is calculated from Equation (4). Planet b, labeled with a black dot in each panel, has a mass $m_b = 0.1/2/8 M_{Jup}$ and a semimajor axis $a_b = 8/10/20$ au. Planet c has a mass range of 0.001–20 $M_{Jup}$ and a semimajor axis range of 1–100 au. All systems have a warp production time of 10 Myr. The white region in the top-right corner corresponds to systems for which warps move beyond 200 au (i.e., the edge of the disk) in 10 Myr. The hatched region corresponds to configurations that fail to satisfy the analytical dynamical stability criterion given in Equation (5). The black dashed box in the top-right corner represents roughly the parameter space to which GPI and SPHERE are sensitive (for a 30 Myr old star at 30 pc; Bowler 2016). The gray dashed box represents the predicted performance of E-ELT METIS L band (for a 30 Myr old star at 30 pc; Quanz et al. 2015). In each panel, planet b dominates the disk morphology in configurations below the dark blue strip (i.e., the bottom-left corner), whereas planet c dominates beyond the strip (i.e., the upper-right corner).</figcaption></figure>

Moreover, the planet's inclination can be inferred from either the overall disk inclination or the warp inclination:

$$i_b = i_{\text{disk}}. \qquad (3)$$

One of the best-studied warped disks is the debris disk of β-Pictoris (β-Pic). A vertical warp at approximately 85 au from the star is observed in optical and near-infrared wavelengths and has been interpreted as an outcome of a giant planet on an inclined orbit (Mouillet et al. 1997; Augereau et al. 2001). Lagrange et al. (2009, 2010) later discovered β-Pic b, a directly imaged planet. Dawson et al. (2011) further modeled β-Pic b, the posited planet generating the warp, using orbital constraints from the directly imaged planet and showed that planet b, rather than another planet, must be responsible for the warp. Nesvold & Kuchner (2015) modeled the β-Pic system including the effect of collisions between planetesimals and the resulting dust grain dynamics, which allowed them to directly compare models to the thermal and scattered-light observations. We discuss the tentative detection of β-Pic c ($m_c = 9 M_{Jup}$ and $a_c = 2.7$ au; Lagrange et al. 2019) in the context of our findings at the end of Section 2.1.1.

### 2.1.1. Mass Estimation Error Using a Single-planet Model

For systems like β-Pic with the detection of both a warped disk and a directly imaged planet, the planet's mass can be constrained from the warp location and the planet's semimajor axis. This approach can be applied to other warped-disk systems that will be discovered in future missions (e.g., James Webb Space Telescope, or JWST; Beichman et al. 2019; Brande et al. 2020; Chen et al. 2019; The Large UV/Optical/IR Surveyor, or LUVOIR; Debes et al. 2019; Roberge et al. 2019; The European Extremely Large Telescope, or E-ELT; Brandl et al. 2014, 2018). However, our interpretation of the detected planet's mass could be compromised if the system has one or more planets undetected. Here we investigate if/how our inference on planet b's mass is compromised if there is a hidden planet c, assuming planet b is the detected planet and planet c is the hidden planet. In our two-planet model, we fix planet b's parameters and explore a range of

parameters for planet c (i.e., $m_c = 0.001–20 M_{Jup}$ and $a_c = 1–100$ au). Doing so allows us to explore a wide range of two-planet mass ratios and semimajor axis ratios. We first explore a simple case where planet b and c are coplanar ($i_b = i_c$). The forced inclination (i.e., the time-averaged inclination) of the planetesimal from two such planets can be written as $i_{\text{forced}} = i_{b,c}$ assuming all planetesimals begin with zero inclination and longitude of ascending node (i.e., a flat disk; see Appendix B). Planetesimals' inclinations oscillate between 0 and $2i_{b,c}$ at different nodal precession rates. Warps are located at semimajor axes where planetesimals are just reaching their peak inclination in their first secular oscillation cycle. When multiple planets sculpt multiple warps, we characterize the warp location using the semimajor axis of the outermost warp because that warp will be the easiest one to observe.

We first study how far the warp location deviates in a two-planet system from a single-planet system. In other words, we study whether the warp location is mostly contributed by a single planet. For all simulations in this study, we assume a solar-mass host star. We consider a warp production time of 10 Myr, motivated by the earlier studies of the β-Pic warp (e.g., Dawson et al. 2011) that used a timescale of 10 Myr. To calculate the warp location deviation, we have

$$\frac{a'_{\text{warp}} - \max(a_{\text{warp,b}}, a_{\text{warp,c}})}{\max(a_{\text{warp,b}}, a_{\text{warp,c}})} 100\%, \qquad (4)$$

where we compare the warp location sculpted by two planets, $a'_{\text{warp}}$, to the warp location sculpted by either planet b or planet c depending on which one has a greater value, $\max(a_{\text{warp,b}}, a_{\text{warp,c}})$. We summarize our results in Figure 3. The left, middle, and right panels of Figure 3 present warp location deviations for three different planet b setups ($m_b = 0.1/2/8 M_{Jup}$ and $a_b = 8/10/20$ au). In the middle and right panels, it would be feasible to detect planet b by direct imaging using GPI or SPHERE (Bowler 2016). Planet c has a mass range of 0.001–20 $M_{Jup}$ and a semimajor axis range of 1–100 au. Dynamically unstable two-planet configurations are shown by
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mass estimation error in percentage for the detected planet b for two-planet configurations vs. the planet b-only configuration, calculated from Equation (6). Planet b's mass and semimajor axis are indicated as the black dot in each panel. The white region in the upper-right corner represents systems whose warps move beyond 200 au in 10 Myr. The hatched region is labeled for unstable configurations given in Equation (5). The black dashed box in the top-right corner represents roughly the parameter space to which GPI and SPHERE are sensitive (for a 30 Myr old star at 30 pc; Bowler 2016). The gray dashed box represents the predicted performance of the E-ELT METIS L band (for a 30 Myr old star at 30 pc; Quanz et al. 2015). Color scales indicate levels of mass estimation error. The darker the color, the more significant the error is caused by the undetected planet c. Configurations under the orange curve and the white curve correspond to the disk median inclination being reduced by 85% and 65%, respectively, compared to single-planet configurations. We categorize these configurations into two types of disk morphology dissimilarities between a single-planet system and a two-planet system, demonstrated in Figures 5 and 6 and explained in Section 2.1.2.</figcaption></figure>

the hatched region using the criterion given in Petrovich (2015):

$$ \frac{a_{\text{out}}(1 - e_{\text{out}})}{a_{\text{in}}(1 + e_{\text{in}})} > 2.4[\max(\mu_{\text{in}}, \mu_{\text{out}})]^{1/3} (a_{\text{out}}/a_{\text{in}})^{1/2} + 1.15, \quad (5) $$

where $a_{\text{out}}$, $e_{\text{out}}$, and $\mu_{\text{out}}$ are the semimajor axis, eccentricity, and planet-star mass ratio for the outer planet, respectively, and $a_{\text{in}}$, $e_{\text{in}}$, and $\mu_{\text{in}}$ are for the inner planet. As shown in Figure 3, the warp location is dominated by a single planet in most two-planet configurations (i.e., configurations in light blue). Planet b dominates the disk morphology in configurations below the dark blue strip (i.e., the bottom-left corner), whereas planet c dominates beyond the strip (i.e., the upper-right corner). Which planet dominates the disk purely depends on the planets' contribution to the nodal precession rate. When planet b and planet c have a similar contribution, the warp location deviation is the highest and can be at most 22%.<fnref n="7" />

We also notice that in systems with one or two massive planets with large semimajor axes, the warp sculpted by it/ them can easily move beyond the observable outer edge of the debris disk (i.e., ~200 au) in a short timescale (i.e., <10 Myr). This region of parameter space is colored white in Figure 3 and subsequent figures.

On the one hand, Figure 3 demonstrates in most configurations one planet dominating the warp location, but on the other hand, it could be problematic when the dominating planet is not the one directly imaged. The mass estimation for the detected planet b could be significantly flawed, if the warp is majorly sculpted by a hidden planet c. To characterize the mass estimation error in situations like this, we estimate planet b's mass from the warp location sculpted by two planets and compare the estimated mass to planet b's true mass. The mass estimation error is

$$ \frac{m'_{b} - m_{b}}{m_{b}} 100\%, \quad (6) $$

where $m'_{b}$ is the estimated mass of planet b, and $m_{b}$ is planet b's true mass. As shown in Figure 4, when planet c gradually dominates the disk morphology, the mass estimation error for planet b grows because we credit the disk feature to the wrong planet. The upper-right black box in each panel indicates the parameter space to which GPI and SPHERE are sensitive (assuming a 30 Myr old star at 30 pc; Bowler 2016); within that box, planet c is detectable around some stars. The gray box represents the predicted parameter space to which the E-ELT METIS L band is sensitive.<fnref n="8" /> The mass estimation error for planet b could be huge if planet c is either massive or close to the warp but not detected. The left panel of Figure 4 ($m_{b} = 0.1 M_{\text{Jup}}$, $a_{b} = 8$ au) emphasizes two-planet configurations for which planet b's estimated mass could deviate significantly from its true mass. As shown in dark blue pixels, a large semimajor axis but low-mass planet c may cause a mass estimation error as high as 10,000%-100,000% for planet b. The reason is that the observed warp is predominantly contributed by planet c with a large semimajor axis (i.e., >100 au). Due to planet c's low mass, the warp has not moved beyond the outer edge of a debris disk yet. Attributing a warp so far from the star to a small semimajor axis planet b results in a huge mass estimation error for the planet. In contrast, if a massive planet b with a large semimajor axis is directly imaged (i.e., the right panel in Figure 4), the mass estimation for the planet will be close to its true mass. The small estimation error in this case is primarily because a massive, large semimajor axis planet c is incompatible with a warp location inside the disk over 10 Myr.

We consider the tentative detection of β-Pic c ($m_{c} = 9 M_{\text{Jup}}$ and $a_{c} = 2.7$ au; Lagrange et al. 2019) in the context of our findings. How would β-Pic c contribute to the observed warp? We can answer the question from Figure 3, the warp location deviation plot. β-Pic b's mass and semimajor axis ($m_{b} = 12 M_{\text{Jup}}$ and $a_{b} = 8.9$ au) are similar to those in the planet b setup in the

<footnote n="7">⁷ The maximum warp location deviation can be derived from setting the same nodal precession rate in a single-planet system and a two-planet system and solving for the warp location ratio.</footnote>

<footnote n="8">⁸ We estimate the E-ELT performance with an inner working angle of 0.0"038 and a planet-star flux contrast of $10^{-8}$ for a 30 Myr old star at 30 pc (Quanz et al. 2015). The planet luminosity is estimated using Equation (1) in Bowler (2016).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Snapshot (10 Myr) of a disk's instantaneous inclination (upper panel) and morphology (lower panel) sculpted by planets b and c ($m_b = 0.1 M_{Jup}$, $a_b = 8$ au, $m_c = 0.01 M_{Jup}$, $a_c = 50$ au; in black) vs. by planet b only ($m_b' = 2.35 M_{Jup}$, $a_b = 8$ au; in blue). To emphasize the morphology differences, the axis aspect ratio is not unity (i.e., y:x is 3:2).</figcaption></figure>

right panel of Figure 3 (Lagrange et al. 2012; Bonnefoy et al. 2014). Figure 3 shows that the parameter space planet c lies in has minimal influence on the warp location. Therefore, β-Pic b remains as the dominant planet sculpting the observed warp. Furthermore, Dawson et al. (2011)<fnref n="9" /> demonstrated β-Pic b must be responsible for the warp, given the planet’s inclination constraint from direct imaging and the observed disk inclination. A second planet is unlikely responsible for the warp because it would be hard for that planet to produce the warped disk morphology without exciting planet b’s inclination. Our finding that a second planet could affect the warp location in a system similar β-Pic (i.e., Figure 3) is compatible because we consider a more general case where the planet’s inclination is unknown.

## 2.1.2. Disk Morphology Dissimilarities between a One-planet System and a Two-planet System

It is no surprise that the disk morphology sculpted by two planets will on some level differ from the disk morphology sculpted by a single planet. If we are able to distinguish these differences, we may infer the need to use a multiplanet model instead of a single-planet one and thus avoid potential mass estimation errors. Here we summarize two types of disk morphology dissimilarities between a one-planet system and a two-planet system. The first is the thickness of the outer disk (i.e., how steeply the outer disk flattens from the warp), labeled as configurations within the orange dashed curve in Figure 4. When the outer disk beyond the warp is primarily sculpted by a low-mass planet c, the warp flattens steeply. Compared to a scenario where planet b is responsible for the warp, the outer disk is expected to be thinner, illustrated in Figure 5. As shown in Figure 5, the disk instantaneous inclination beyond the warp is much lower in a planet b and c-sculpting system (in black) than a planet b-sculpting system (in blue), even though both systems have the same warp location. Consequently, we observe a much sharper slope beyond the warp. Arrows in Figure 4 indicate the direction in which this feature becomes more obvious and observable (i.e., lower masses and larger semimajor axes for planet c). The second morphology difference between warps sculpted by one versus two planets

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Snapshot (10 Myr) of the disk's instantaneous inclination (upper panel) and morphology (lower panel) sculpted by planets b and c ($m_b = 0.1 M_{Jup}$, $a_b = 8$ au, $m_c = 0.001 M_{Jup}$, $a_c = 100$ au) vs. by planet b only ($m_b' = 13.2 M_{Jup}$, $a_b = 8$ au; in blue).</figcaption></figure>

corresponds to configurations under the white curve in Figure 4. In these configurations, we expect to observe two warps sculpted by planet b and c individually, illustrated in Figure 6. Given planet c’s large semimajor axis and low mass, it sculpts an outer warp at around 105 au after 10 Myr evolution. Meanwhile, planet b sculpts a warp at around 28 au. The region between planet b and c remains flat and thus we observe two separate warps. For a given evolution timescale, this double-warp feature is only present when planet b’s mass and/or semimajor axis is sufficiently small, so we do not see this feature in the middle and right panels of Figure 4 (i.e., more massive planet b). It has a strong dependence on the warp production time, because after the inner warp and the outer warp meet, we will no longer observe multiple warps.

### 2.1.3. The Effect of Secular Resonances

The forced inclination of planetesimals sets the warp height and the overall disk inclination. We can understand the forced inclination as a time-averaged inclination on which planetesimals have their inclination oscillations centered. In systems where planet c has the same inclination as planet b, the forced inclination is a constant with semimajor axis and is the same as the planets’ inclination. However, in systems where planet c and planet b have different inclinations, the forced inclination is not a constant but depends on the planets’ masses, semimajor axes, and inclinations. First, both planets will force nearby planetesimals to their inclined orbits, and therefore, the observed disk inclination will differ from a single-planet case. Second, secular resonances modify the disk’s forced inclination. When the nodal precession rates of planetesimals reach one of the characteristic frequencies of the planetary system, the planetesimals get excited to highly inclined orbits.<fnref n="10" /> Usually, the secular resonance has little effect on the morphology of a warped disk. Because we observe the disk edge on, gaps cleared by resonances would not be easily observable. At certain warp production timescales, however, a warp could move near one of the secular resonance locations, and the warp location could be modified by the resonance. In Figure 7, we present a case where the secular resonance reduces the semimajor axis of a warp but increases the warp height. In this situation, we are at risk of moderately underestimating the planet’s mass but overestimating the planet’s inclination. Fortunately, secular

<footnote n="9">⁹ We found that Dawson et al. (2011) mistakenly used 1 solar mass instead of 1.75 solar mass in calculating the warp production time in their Figure 1. This error does not affect the conclusions of Dawson et al. (2011).</footnote>

<footnote n="10">¹⁰ In Equation (18), the forced inclination or eccentricity approaches infinity as $B \to f_i$ or $A \to g_i$ in the first-order secular approximation.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Illustration of an inclination-type secular resonance modifying the location and height of a warp. Both simulations (i.e., blue and black) have an evolution timescale of 5.14 Myr (i.e., one secular cycle) and the same planet masses and semimajor axes: $m_b = 2 M_{\text{Jup}}$, $a_b = 10$ au, $m_c = 0.1 M_{\text{Jup}}$, and $a_c = 50$ au. In the blue system, $i_b = i_c = 10^\circ$ (i.e., the planets are coplanar), whereas in the black system, $i_b = 10^\circ$ and $i_c = 0^\circ$ (i.e., the planets have a mutual inclination of $10^\circ$). The upper panel presents the instantaneous inclinations of planetesimals at different semimajor axes. The gray curve is the forced inclination of the black system, with a secular resonance around 75 au that modifies the warp location. The lower panel is a projected view of disk morphology showing differences caused by the secular resonance.</figcaption></figure>

resonance can only modify the warp at a particular warp production time. Nevertheless, it would be helpful to take this possibility into consideration when characterizing planet properties.

### 2.1.4. Inclination Estimation Error Using a Single-planet Model

In the last part of our exploration of the warp feature, we investigate how our interpretation of the detected planet's inclination is affected by the presence of an additional planet. We mentioned earlier that the disk inclination could be different from the single-planet case if planet c were to have a different inclination from planet b. To study the strength of this effect, we build a two-planet model in which planet b has an inclination of $10^\circ$ ($i_b = 10^\circ$), whereas planet c has an inclination of zero ($i_c = 0^\circ$). We calculate the median forced inclinations of 2000 evenly distributed planetesimals from 20 to 200 au and use it as a disk inclination indicator. The following equation is used to compare forced inclinations of the disk:

$$ \frac{|i'_{\text{disk}} - i_{\text{disk}}|}{i_{\text{disk}}} 100\%, \qquad (7) $$

where $i'_{\text{disk}}$ is the median disk forced inclination in systems with planet b and planet c, and $i_{\text{disk}}$ is in systems with planet b only with $i_b = 10^\circ$. As shown in Figure 8, the disk forced inclination is reduced significantly when planet c is more massive than planet b. In other words, the median forced inclination of the disk is dominated by the inclination of the more massive planet. To illustrate this feature, we plot the disk forced inclination for configurations with different planet mass ratios in Figure 9. As planet c's mass increases, the disk forced inclination decreases from $10^\circ$ to $0^\circ$ (dark blue curve). It is noticeable that the singularity in the forced inclination shifts to greater semimajor axes as planet c's mass increases. This shift is mostly caused by a higher nodal precession rate of nearby particles as a more massive planet c shifts the secular resonance to a greater distance. The triangle-shaped feature in the contour in Figure 8 when planet c has a semimajor axis $a_c > 70$ au and low mass is caused by our definition of the outer edge of a

debris disk as 200 au. Figure 10 illustrates that as planet c's semimajor axis increases, the secular resonance occurs at a larger semimajor axis and its width increases. Because we truncate the disk at 200 au (i.e., shown by the vertical gray dashed line in Figure 10), our calculation of the median forced inclination is affected by the high forced inclinations of planetesimals near resonance. In addition, as planet c's semimajor axis increases, the semimajor range of planetesimals with zero forced inclination increases. If we extend the disk outer edge to 500 au, we find that the forced inclination of the planetesimals drops back to $10^\circ$, which still follows the inclination of the more massive planet. In interpreting observations, we should make use of the observed disk outer edge, on which the disk median inclination depends.

In summary, for the warp feature, our inference of planet b's mass could be significantly flawed using a single-planet model if an undetected planet c is either massive or close to the warp. Disk morphology dissimilarities may let us distinguish between one-planet and two-planet models. The overall warped disk inclination follows the inclination of the more massive planet. The inference of the detected planet's inclination could be significantly different from the true value if a hidden planet is more massive than the detected planet. We also find secular resonances could affect the warp location and height and require detailed modeling for specific systems.

### 2.2. Offset

In some exoplanetary systems (e.g., Fomalhaut, HR 4796A; Telesco et al. 2000; Kalas et al. 2005; MacGregor et al. 2017), the geometric center of the debris disk or the debris ring has been observed to be offset from its central star. Wyatt et al. (1999) explained the observed offset as a result of gravitational sculpting of a nearby planet on an eccentric orbit. The forced eccentricity ($e_{\text{forced}}$), a time-averaged eccentricity of a planetesimal at the semimajor axis $a$ from planet b, can be expressed as

$$ e_{\text{forced}} = \frac{b_{3/2}^{(2)}(\alpha)}{b_{3/2}^{(1)}(\alpha)} e_b, \qquad (8) $$

where $e_b$ is the posited planet b's eccentricity that generates the offset, $\alpha$ is the planet-planetesimal semimajor axis ratio in Equation (17), and $b_{3/2}^{(1)}$ and $b_{3/2}^{(2)}$ are the Laplace coefficients in Equation (21). Unlike the forced inclination, the forced eccentricity depends on the semimajor axis ratio of the planetesimal to the planet (i.e., the forced eccentricity gets smaller when the planetesimal is more widely separated from the planet). Although the forced eccentricity itself is independent of the sculpting planet's mass,<fnref n="11" /> the semimajor axis and eccentricity of the posited planet can be inferred from the offset feature using the forced eccentricity equation.

For the offset feature, we focus on the question: how does a second planet with an eccentricity different from the dominating planet affect the observed offset of the debris ring and the debris disk? To study planet c's effect on a debris ring, we consider two-planet configurations for which planet b has $e_b = 0.2$ and planet c has $e_c = 0$. We study the forced eccentricity of a ring near planet b and see how planet c's low eccentricity reduces the offset of that ring. In Figure 11, we

<footnote n="11">¹¹ The secular precession timescale still depends on the mass of planet b.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Median disk forced inclination deviation in percentage for two-planet configurations vs. planet b-only configuration, calculated from Equation (7). Planet b, labeled with a black dot in each panel, has a mass of $m_b = 0.1/2/8 M_{Jup}$ and a semimajor axis of $a_b = 8/10/20$ au. Planet c has a mass range of 0.001–20 $M_{Jup}$ and a semimajor axis range of 1–300 au. Planet b has an inclination of 10°, whereas planet c has zero inclination. The hatched region corresponds to dynamically unstable configurations using the criterion given in Equation (5). The black dashed box in the top-right corner represents roughly the parameter space to which GPI and SPHERE are sensitive (for a 30 Myr old star at 30 pc; Bowler 2016). The gray dashed box represents the predicted performance of the E-ELT METIS L band (for a 30 Myr old star at 30 pc; Quanz et al. 2015). The darkened region in each panel corresponds to planet c causing significant disk forced inclination deviations (i.e., $i'_{disk} \to 0^{\circ}$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Disk forced inclination as a function of the semimajor axis. In each system, planet b is a 2 Jupiter-mass planet with a semimajor axis of 10 au and a 10° inclined orbit, the same as the middle panel of Figure 8. Planet c has a semimajor axis of 30 au and zero inclination. Different curves represent disk forced inclination profiles in systems with different masses for planet c ($m_c = 0.1, 1, 10 M_{Jup}$). Increasing planet c's mass forces the disk to lower forced inclinations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Disk forced inclination as a function of the semimajor axis. In each system, planet b is a 2 Jupiter-mass planet with a semimajor axis of 10 au and a 10° inclined orbit, the same as the middle panel of Figure 8. Planet c has a mass of 0.1 $M_{Jup}$ and zero inclination. Different curves represent disk forced inclination profiles in systems with planet c having different semimajor axes ($a_c = 100, 150, 200$ au). Systems in which planet c has a large semimajor axis are affected by a secular resonance (i.e., the spike around 180 au in $a_c = 100$ au case).</figcaption></figure>

present the forced eccentricity reduction for an offset ring 2.5 Hill radii away from planet b caused by planet c. We place the debris ring 2.5 Hill radii outside planet b's orbit, assuming it is at the boundary of planet b's chaotic zone (Wisdom 1980; Duncan et al. 1989; Murray & Holman 1997).<fnref n="12" /> As shown in

Figure 11, in most parameter space, planet c barely affects the forced eccentricity of planetesimals near planet b. Planet c only significantly affects the offset ring when planet b has a low mass (i.e., left panel) and planet c is both massive and close to the ring. The forced eccentricity of the ring in these configurations is reduced by at most 35%. In other words, eccentricities of planetesimals are primarily dominated by the planet nearby.

We also generalize our study from a debris ring to a full debris disk. We investigate two planets with different eccentricities and their combined effect on the median forced eccentricity of 2000 planetesimals evenly distributed from 20 to 200 au. Because we study the forced eccentricity of the whole disk instead of a ring close to a planet, secular resonances become important and complicate the characterization of eccentricities. Yelverton & Kennedy (2018) recently investigated how planets' masses, semimajor axes, and eccentricities determine the locations, timescales, and widths of two exterior secular resonances in two-planet systems. In our study, two interior resonances and two exterior resonances all affect disk eccentricity. We find that the median forced eccentricity of the disk can be strongly affected by secular resonances, making it challenging to characterize the planet eccentricity from the disk eccentricity. In general, the planet with greater semimajor axis dominates the disk eccentricity, assuming the system is old enough for the disk to have experienced one secular oscillation cycle.

In summary, for the offset feature, we find the forced eccentricity of a debris ring is usually primarily sculpted by the nearest planet. A disk's forced eccentricity is more complicated to model due to secular resonances, but the overall disk eccentricity follows the eccentricity of the planet with the greater semimajor axis inside or within the disk.

In our model, we assumed the eccentricity of the debris ring is secularly excited by a nearby planet on an eccentric orbit. Nesvold et al. (2016, 2017) demonstrated an exterior companion outside the disk may also excite planetesimal eccentricities through either secular perturbations (if the companion and disk are coplanar) or Kozai–Lidov mechanism (if the companion and disk are mutually inclined). Our conclusion still stands for the coplanar case: the forced eccentricity of a debris ring is usually primarily sculpted by

<footnote n="12">¹² $\Delta a_{chaotic} \sim 1.3\mu^{2/7} \approx 2.5 R_{Hill}$ for Jupiter-mass planets.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Forced eccentricity reduction on a ring 2.5 Hill radii away from an eccentric orbit planet b ($a_{\text{ring}} = a_b + 2.5 R_{\text{Hill,b}}$, $e_b = 0.2$) caused by a zero-eccentricity planet c ($e_c = 0$). Planet b, labeled with a black dot in each panel, has a mass of $m_b = 0.1/2/8 M_{\text{Jup}}$ and a semimajor axis of $a_b = 8/10/20$ au. Planet c has a mass range of 0.001–20 $M_{\text{Jup}}$ and a semimajor axis range of 1–300 au. The hatched regions are for dynamically unstable configurations using the criterion given in Equation (5). The black dashed box in the top-right corner represents roughly the parameter space to which GPI and SPHERE are sensitive (for a 30 Myr old star at 30 pc; Bowler 2016). The gray dashed box represents the predicted performance of the E-ELT METIS L band (for a 30 Myr old star at 30 pc; Quanz et al. 2015). In most parameter space, planet c barely affects the forced eccentricity of the debris ring near to planet b.</figcaption></figure>

the nearest perturber. How a misaligned, massive perturber affects the eccentricity and inclination of the debris ring when there is an inner planet is deferred to future work.

### 3. Resonant Features in Two-planet Systems

Debris disk features are sculpted by various types of planet-disk resonant interactions. A planet can remove nearby planetesimals, creating a gap, and sharpen the disk's edge because of the planet's overlapping resonances. In another type of resonant interaction, a planet captures planetesimals into its orbital resonances while migrating and redistributes the planetesimals to certain longitudes, creating an azimuthal asymmetry in the disk. Properties of these resonant disk features (e.g., the gap width or the shape of clumpy features) can be used to constrain the sculpting planet's mass, semimajor axis, and migration history (Wisdom 1980; Quillen 2006; Wyatt 2018). Here we investigate if/how our inferences of these planet properties are compromised using a single-planet model when multiple planets reside in the system for edges and gaps (Section 3.1) and azimuthal asymmetries (Section 3.2).

#### 3.1. Edges and Gaps

Sharp edges and gaps have been observed in many debris disk systems (e.g., Fomalhaut, HR 8799, HD 107146; Kalas et al. 2005; Ricci et al. 2015; Wilner et al. 2018). The width of a gap and the sharpness of an edge can constrain the underlying planet's mass. Wisdom (1980) derived a simple relation of the width of the chaotic zone (i.e., a region where planetesimals are removed because of the planet's overlapping resonances) and the posited planet's mass for a planet on a circular orbit:

$$ \Delta a_{\text{chaotic}}/a_b \propto \mu^{2/7}, \qquad (9) $$

where $\Delta a_{\text{chaotic}}$ is the width of the chaotic zone, $a_b$ is the sculpting planet's semimajor axis, and $\mu$ is the planet-star mass ratio. This relation is also known as the “2/7” law. Follow-up studies took into account the planet's eccentricity (e.g., Gladman 1993; Chiang et al. 2009) and the planetesimal removal timescale (e.g., Morrison & Malhotra 2015). One of the most well-studied disks featuring edges and gaps is Fomalhaut's debris disk. Quillen (2006) applied the “2/7” law to Fomalhaut and constrained the mass of Fomalhaut b (Fom b), the posited

planet that would clear the chaotic zone and sharpen the inner edge of the belt. Chiang et al. (2009) later constrained Fom b's mass, semimajor axis, and eccentricity using orbital constraints from the directly imaged planet (Kalas et al. 2008). Unlike the β-Pic system, later studies showed that the detection by Kalas et al. (2008) might not be the planet responsible for the feature. On the one hand, Chiang et al. (2009) still provide a good example of how to characterize the sculpting planet's mass from the width of the chaotic zone and the sharpness of the inner edge; on the other hand, Fom b is a cautionary tale about the dangers of assuming we know what planet is causing the features.

#### 3.1.1. Mass Estimation Error Using a Single-planet Model

If an edge or gap is sculpted by multiple planets, we expect the width of the chaotic zone to expand, and the sharpness of the inner edge may also be affected. In this section, we aim to answer the question: how would an undetected planet c impact our inference of planet b's mass from an edge or gap? To do so, we build a two-planet model including a Jupiter-mass planet b with a semimajor axis of 16 au ($m_b = 1 M_{\text{Jup}}$ and $a_b = 16$ au) and a planet c with a mass range of 0.0005–1 $M_{\text{Jup}}$ and a semimajor axis range of 1–15 au. This parameter space avoids planet c's chaotic zone exceeding planet b's chaotic zone and thus becoming the dominant planet. For simplicity, we also assume planet b and planet c have zero eccentricity and are coplanar with the disk. Resonant features, unlike secular features, do not have a known analytical solution in the four-body case (i.e., a host star, planet b, planet c, and a planetesimal). To investigate these features, we perform N-body simulations using REBOUND (Rein & Liu 2012; Rein & Spiegel 2015; Rein & Tamayo 2015) and REBOUNDx (Tamayo et al. 2019). For each two-planet configuration, we include 2000 test particles evenly distributed from 10 to 28 au with zero eccentricity and inclination, and random mean anomalies. Each system is integrated for 5 Myr with a time step of 1% of planet c's orbital period using the WHFast integrator (Rein & Tamayo 2015). The width of the chaotic zone is computed from the difference in semimajor axis between planet b and the median of the 10 innermost planetesimals beyond planet b's orbit. We choose this metric because of its insensitivity to transient unstable test particles. We first apply the metric in the single-planet case and characterize how well it works. We consider a planet mass range of 1–5 $M_{\text{Jup}}$, characterize the widths
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Mass estimation error for planet b ($m_b = 1 M_{Jup}$ and $a_b = 16$ au) calculated from Equation (11) according to the width of the chaotic zone sculpted by both planet b and planet c. The hatched region represents dynamically unstable two-planet configurations over 10 Myr (i.e., $e_c > 1$, collisions or ejections). The color scale demonstrates levels of mass estimation error labeled by the color bar. A planet c with an orbit close to the inner disk edge and with a comparable mass to planet b leads to the greatest mass estimation error for planet b.</figcaption></figure>

of the chaotic zones after 5 Myr evolution, and fit the widths with Equation (9). We find that the planet's mass scales with the width of the chaotic zone as

$$ \Delta a_{\text{chaotic}} = 1.75 \left( \frac{m_b}{m_*} \right)^{2/7} a_b, \qquad (10) $$

where the coefficient 1.75 is from our fitting. The relation has some scatter and this variation can cause a small mass estimation error for the planet (i.e., $\lesssim 25\%$). For two-planet configurations, the mass estimation error is calculated as follows:

$$ \frac{|m'_b - m_b|}{m_b} 100\%, \qquad (11) $$

where $m'_b$ is the estimated mass for planet b using Equation (10) and $m_b$ is planet b's true mass ($m_b = 1 M_{Jup}$). In Figure 12, we present the mass estimation error for planet b in such two-planet configurations. The hatched region represents dynamically unstable two-planet configurations over 10 Myr (i.e., if $e_c > 1$ or any ejections or collisions occur). As shown in Figure 12, planet c causes negligible mass estimation errors in most configurations, due to its modest expansion of the chaotic zone. Only a planet c with a large semimajor axis and comparable mass to planet b leads to a relatively large estimation error (i.e., 100%–250%). In other words, our inference of planet b's mass may be flawed when an undetected planet c has a comparable mass and orbit to planet b.

### 3.1.2. Further Planet Mass Constraints from the Edge Sharpness

The sharpness of the disk edge may place a further constraint on the planet mass. Quillen (2006) has shown that the eccentricity dispersion of planetesimals at the boundary of the chaotic zone is correlated with the sculpting planet's mass as $u_e \sim \mu^{3/7}$. For example, a low-mass planet disperses planetesimals' eccentricities by a relatively small amount and thus sculpts a sharp edge, and vice versa for a massive planet. Given the observed sharpness of the disk edge, we may constrain the sculpting planet's mass. In a single-planet system, we find that the observed edge sharpness is affected by planet

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Kernel density estimations for the planetesimals' semimajor axis distribution (upper panel) and stellar radial distance distribution (lower panel) sculpted by planet b only (in darkest blue) or by both planet b and planet c with different masses for planet c (in lighter blues). We use Gaussian kernels with a bandwidth of 0.25 au. Planet b is a Jupiter-mass planet with a semimajor axis of 16 au. Planet c has $a_c = 12$ au and $m_c = 0.05/0.5/1 M_{Jup}$. Gray dashed lines indicate the location of mean-motion resonances of planetesimals with planet b.</figcaption></figure>

b's mass via its effect on edge location and its perturbation of planetesimal eccentricities, in agreement with Quillen (2006). Chiang et al. (2009) demonstrated that the semimajor axis distribution of planetesimals could also affect the observed edge sharpness, because planetesimals on elliptical orbits at greater semimajor axes could cross the inner edge boundary and smooth the disk edge. To constrain the sculpting planet's mass from the observed sharpness of the disk edge, we have to include the contribution from planetesimals at greater semimajor axes in the model.

We conduct a series of case studies for two-planet configurations, aiming to explore an additional planet c's contribution to the edge sharpness. If planet c modifies the edge sharpness from that sculpted by planet b only, we can infer the need to use a multiplanet model instead of a single-planet one. For each selected configuration, we include 4000 test particles near the disk inner edge (i.e., 18–23 au) and integrate the system for 5 Myr using the IAS15 integrator (which treats close encounters more accurately than other integrators, e.g., WHFast; Rein & Spiegel 2015). For two-planet configurations, we set $m_b = 1 M_{Jup}$, $a_b = 16$ au, $m_c = 0.05/0.5/1 M_{Jup}$, and $a_c = 12$ au. Planet c is placed close to planet b to have an observable effect on the disk edge. Kernel density estimation is applied to the semimajor axis and radial distance distribution of planetesimals. A Gaussian kernel with a bandwidth of 0.25 au is used.

As shown in the upper panel of Figure 13, increasing planet c's mass moves the disk edge to a greater semimajor axis. Each peak of semimajor axis distribution is slightly offset from locations of planet b's mean-motion resonances (i.e., labeled in gray dashed lines), because of the planets' perturbation of the eccentricity of planetesimals (the higher the planet mass, the stronger the perturbation). We observe a slightly steeper inner edge slope for the planetesimal semimajor axis distribution when planet c is 0.5 or 1 Jupiter mass, compared to planet c with 0.05 Jupiter mass. This slight sharpness change could be a
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Kernel density estimations for the planetesimals’ semimajor axis distribution (upper panel) and stellar radial distance distribution (lower panel) sculpted by both planet b and planet c ($m_b = 1 M_{Jup}$, $a_b = 16$ au, $m_c = 0.5 M_{Jup}$, $a_c = 12$ au; in black) or by planet b only ($m_b = 2 M_{Jup}$, $a_b = 16$ au; in blue). We use Gaussian kernels with a bandwidth of 0.25 au. Gray dashed lines indicate the location of mean-motion resonances of planetesimals with planet b. Edge sharpness dissimilarity is ambiguous in distinguishing between one-planet and two-planet systems.</figcaption></figure>

result of a 0.5 or 1 Jupiter-mass planet c expanding the chaotic zone so the disk edge is close to planet b’s 3:2 resonance; therefore, this subtle steepness change is caused by the change in edge location rather than the eccentricity dispersion. However, in the radial distance distribution presented in the lower panel of Figure 13, the steepness change is less observable. The flatter slope in the radial distance distribution is because, as described by Chiang et al. (2009), planetesimals on elliptical orbits from greater semimajor axes cross the inner boundary and smooth the disk edge. It is difficult to distinguish among various two-planet configurations using the edge sharpness in the radial distribution of planetesimals, which is closer to what we observe than the semimajor axis distribution.

We also compare the edge sharpness for two-planet configurations to single-planet configurations. We show one example in Figure 14. In the single-planet system, planet b is a 2 Jupiter-mass planet with a semimajor axis of 16 au ($m_b = 2 M_{Jup}$, $a_b = 16$ au). In the two-planet system, planet b is a Jupiter-mass planet with the same semimajor axis and planet c is a 0.5 Jupiter-mass planet with a semimajor axis of 12 au ($m_b = 1 M_{Jup}$, $a_b = 16$ au, $m_c = 0.5 M_{Jup}$, $a_c = 12$ au). The two configurations produce slightly different chaotic zone widths. However, as shown in Figure 14, we observed no significant difference in edge sharpness in either the planetesimal semimajor axis distribution or stellar radial distance distribution. We conclude edge sharpness dissimilarity cannot unambiguously distinguish between one-planet and two-planet systems.

If a planet has its mass estimated from the width of the chaotic zone, we could in principle model the disk edge using a single-planet model and compare the modeled edge sharpness to the observed one to validate the model. For example, imagine we directly imaged a planet with a semimajor axis of 16 au in a gap disk. Using the mass-chaotic zone width relation in Equation (10), we would infer the planet to be 3 Jupiter masses. If the planet is truly a Jupiter-mass planet with an interior companion, we would

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Eccentricity vs. semimajor axis of planetesimals in a system with planet b only (gray dots), planets b and c in 2:1 mean-motion resonance (blue triangle up), and planets b and c with a similar configuration but not in 2:1 mean-motion resonance (orange triangle down). Planet b is a Jupiter-mass planet with a semimajor axis of 16 au in all three systems. Planet c is a Jupiter-mass planet with a semimajor axis of 10 au in both systems, having its resonant argument librating in the blue system but not in the orange system.</figcaption></figure>

have overestimated that planet’s mass by 200%. However, a 3 Jupiter-mass planet would produce a greater planetesimal eccentricity dispersion compared to a Jupiter-mass planet with an interior companion, which would affect the edge sharpness. Nevertheless, we would still face the problem (i.e., Chiang et al. 2009) that the semimajor axis distribution of planetesimals also affects the observed edge sharpness.

### 3.1.3. The Effects of Two Planets in Mean-motion Resonance

We also investigate planet b and c’s joint effect on the disk when they are in mean-motion orbital resonance (MMR). When planet c is in resonance with planet b and at the same time planet b is in resonance with outer planetesimals, a resonant chain is formed, which could affect planetesimals’ orbits and stability. In Figure 15, gray dots show the eccentricity versus semimajor axis of planetesimals in a planet b-only system ($m_b = 1 M_{Jup}$ and $a_b = 16$ au). Adding in a Jupiter-mass planet c with a semimajor axis of 10 au ($m_c = 1 M_{Jup}$ and $a_c = 10$ au) expands the chaotic zone from 19.5 au to 21.5 au, illustrated by blue and orange triangles. In the two-planet-in-resonance configuration (in blue), planets b and c are trapped in 2:1 MMR with the libration of the resonant argument $\phi$ written as

$$ \phi = 2\lambda_b - \lambda_c - \omega_{b,c}, \quad (12) $$

where $\lambda_{b,c}$ are the mean longitudes of planet b and c, and $\omega_{b,c}$ are the longitudes of pericenter of planets b and c. For the not-in-resonance configuration (in orange), planets b and c have similar masses and semimajor axes but are not librating in resonance. We observe no significant difference in the width of the chaotic zone or the sharpness of the edge between two types of configurations. Planet b and planet c in MMR can help stabilize planetesimals at certain semimajor axes, e.g., planet b’s Trojans, compared to the planet b and planet c not-in-resonance case.

### 3.1.4. Summary of Edge and Gap Features

In summary, we explored how our inference of the posited planet’s mass from the width of the chaotic zone and the edge sharpness can be affected by a hidden planet. For all two-planet configurations, the posited planet remains the dominant planet by construction. We find that the posited planet’s mass estimation may be flawed only for a hidden planet with a comparable mass and orbit to planet b, given the mass estimation from the width of
the chaotic zone. (An additional planet has little effect on the sharpness of the planetesimals' radial distribution but may subtly affect the semimajor axis distribution.) In observations, two planets with similar masses and semimajor axes are likely to be both detected. In the unfortunate case where only one planet is detected, the mass estimation error for the detected planet could be as high as 250%. Furthermore, planet b and planet c in MMR do not enhance the chaotic behavior, compared to planet b and c near MMR.

## 3.2. Azimuthal Asymmetries

Clumpy structures observed in debris disks (e.g., β-Pic; Matrà et al. 2017, 2019) may result from resonant interactions between planets and planetesimals (e.g., Wyatt 2003). As a planet migrates inward/outward, it captures and redistributes planetesimals inside/outside its orbit to certain longitudes, resulting in disk asymmetries. Wyatt (2003) showed that the probability of planetesimals being captured by a planet depends on the planet's mass, semimajor axis, and migration rate. Follow-up studies took into account additional factors such as a planet's eccentricity (e.g., Reche et al. 2008). Consequently, if an azimuthally asymmetric distribution of planetesimals is observed, we can place constraints on the underlying planet's mass, semimajor axis, and migration history.

A second planet, migrating along with the planet primarily sculpting the disk, may modify the resonant behavior of planetesimals and thus affect the asymmetries. To understand planet c's effect on the disk morphology, we investigate a few two-planet configurations including a Neptune-mass planet b and a Neptune-mass planet c. Both planets migrate outward through the disk by exchanging angular momentum with planetesimals (Fernandez & Ip 1984). We consider a semimajor axis evolution of $a_0e^{t/20 \text{ Myr}}$ for both planets, where $a_0$ is the planet's initial semimajor axis and $t$ is time since migration began. With an initial semimajor axis of 40 au for planet b ($a_{0,b} = 40$ au), this evolution means planet b migrates from 40 au to 66 au over 10 Myr. The migration rate we use here is motivated by Wyatt (2003) to optimize the capture probabilities into the 3:2 resonance for planetesimals. Three different initial semimajor axes of planet c are investigated: $a_{0,c} = 36$ au, $a_{0,c} = 30.5$ au, and $a_{0,c} = 25$ au, which correspond to final semimajor axes of 59.4 au, 50.3 au, and 41.2 au, respectively. We investigate how a planet c with different semimajor axes affects the resonant behavior of planetesimals. We simulate planet migration using the *modify_orbits_forces* routine in the `REBOUNDx` library<fnref n="13" /> (Tamayo et al. 2019). We integrate 5000 planetesimals (i.e., massless test particles) evenly distributed from 40 to 140 au for 10 Myr and use the IAS15 integrator to accurately handle close encounters. Our simulation outputs for the planet b-only configuration and three two-planet configurations are presented in Figure 16.

Planetesimal positions (left panels) and eccentricity versus semimajor axis distributions of planetesimals (right panels) are shown. The blue and orange dots represent planet b and planet c, respectively. In the planet b-only configuration (a1, a2 in Figure 16), we find planetesimals outside planet b's orbits are captured in several resonances (e.g., 2:1, 5:3, 3:2, 4:3, 5:4, 6:5). A significant amount of planetesimals are captured into the 3:2 resonance with a wide range of eccentricities. Two symmetric clumps observed in Figure 16 (a1) are specifically these with

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Case studies of the azimuthal asymmetry feature for a single-planet configuration (panel (a)) and three two-planet configurations (panel (b)–(d)). Left: planetesimal positions; right: eccentricity vs. semimajor axis distributions after 10 Myr of evolution. The blue and orange dots represent the final properties of planet b and planet c, respectively. Planets in all simulations are Neptune mass and migrate outward through the disk, following $a_0e^{t/20 \text{ Myr}}$. Gray dashed lines indicate locations of resonances with planet b. The orange dashed line demonstrates where periapses of planetesimals cross planet c's orbit (i.e., $e = 1 - a_c/a$). Panel (a) presents a disk sculpted by the migration of planet b from 40 to 66 au. Panel (b) presents the same planet b along with the migration of planet c from 36 to 59.4 au. Panel (c) presents the same planet b along with the migration of planet c from 30.5 to 50.3 au. Panel (d) presents the same planet b along with the migration of planet c from 25 to 41.2 au.</figcaption></figure>

high eccentricities (i.e., $e > 0.3$). When planet c has a large semimajor axis (e.g., $a_{0,c} = 36$ au or $a_{0,c} = 30.5$ au; b1, b2 or c1, c2 in Figure 16), planetesimals captured in 6:5, 5:4, and 4:3 resonances are greatly reduced because the planetesimals cross the orbit of planet c. Although a large amount of planetesimals are still trapped in 3:2 resonance with planet b, highly eccentric planetesimals ($e > 0.3$) are missing because the inner planet c destabilizes the outer eccentric planetesimals. Orange dashed lines indicate where periapses of planetesimals cross planet c's orbit (i.e., $e = 1 - a_c/a$). In contrast, when planet c has a small semimajor axis ($a_{0,c} = 25$ au; c1, c2 in Figure 16), planet c's orbit is too small to affect planetesimals trapped in resonance with planet b, even those on highly eccentric orbits. Meanwhile, planet c clears the inner region of the disk and makes the clumpy structures even easier to observe, compared to the planet b-only case.

In summary, the azimuthal asymmetry feature is primarily dominated by a single planet. However, a second planet may affect the asymmetries by destabilizing planetesimals crossing its orbit. Disks with no observed clumpy structures could be the result of multiple migrating planets.

We only explored a few case studies for a second planet with different semimajor axes. Although we placed two planets

<footnote n="13">¹³ REBOUNDx (2.19.2), https://github.com/dtamayo/reboundx.</footnote>
close to resonance for some configurations (e.g., Figure 16 panel (c) near to 3:2 and panel (d) near to 2:1), in none of our case studies were the two planets captured into resonance with each other during migration, and this scenario's effect on the feature would be interesting to explore further. A full parameter space study of two-planet configurations (e.g., mass ratios, eccentricities, in MMR, etc.) is beyond the scope of this paper, but worthwhile to explore to fully understand a second planet's effect on azimuthal asymmetries.

## 4. Synodic Features in Two-planet Systems

A planetesimal's orbit can be modified by a nearby planet through close encounters. At each conjunction, the planetesimal's eccentricity and inclination are altered by the sculpting planet, along with a minor change in its semimajor axis. If a planet sculpts a ring of planetesimals, these planetesimals will show a distribution of eccentricities and inclinations, resulting in a radially thickened ring with a vertical scale height. We categorize this type of planet-disk interaction as a synodic feature. The underlying posited planet's mass and semimajor axis can be constrained from the normalized thickened ring width (e.g., Rodigas et al. 2014) and the vertical scale height (e.g., Quillen et al. 2007). Both disk features are particularly useful for characterizing planet properties in systems with no detected planet. In this section, we investigate if/how using a single-planet model compromises our inferences of these planet properties when multiple planets reside in the system.

### 4.1. Radially Thickened Rings

As a planet clears nearby material to sculpt the disk edge (Section 3.1), it simultaneously stirs up planetesimals outside its chaotic zone; if the planet is near a ring of planetesimals (e.g., Fomalhaut or HR 4796A), we will observe a radially thickened ring. How much a ring gets thickened is directly related to the sculpting planet's mass and how close the planet is to the ring. With the assumption that the observed ring lies just beyond the sculpting planet's chaotic zone, we can infer the planet mass from the thickened ring width (Chiang et al. 2009; Rodigas et al. 2014). Rodigas et al. (2014) confirmed a simple linear relation between the planet mass and the thickened ring width:

$$ \mu \propto \text{nFWHM}, \qquad (13) $$

where $\mu$ is the planet-star mass ratio and nFWHM is the normalized FWHM of the scattered-light debris ring. Because we do not know the initial ring width for an observed system, this relationship places an upper limit on the planet's mass.

We are interested in how a second planet in the system affects the ring width and consequently compromises our estimation of the first planet's mass. To begin with, we confirm the mass-normalized ring width relation given in Rodigas et al. (2014) for a single-planet system. This step is to affirm that the linear relation found by Rodigas et al. (2014) for scattered light also applies to planetesimal distributions. We consider a planet mass range of 0.3–10 $M_{Jup}$ and place a zero-width ring of 2000 test particles outside the planet's orbit. Both the planet and planetesimals have zero eccentricity. The ring is placed as close to the planet as possible to have the strongest gravitational effect. The ring's semimajor axis is determined by the following criterion: we require that more than 99% of planetesimals must have a final semimajor axis between 0.1 and 10 of the initial semimajor axis but outside the sculpting planet's Hill sphere after 2000 orbit periods of the sculpting planet. This criterion, compared to Rodigas et al. (2014), is more stringent to make sure the debris ring is stable for at least a Myr timescale (i.e., ~50% of planetesimals remain after 1 Myr). The ring width is computed as the range of the radial distances of the inner 99% of planetesimals and further normalized by dividing by the semimajor axis of the ring center. We find the normalized ring width varies in time due to planet–planetesimal synodic and resonant interactions.<fnref n="14" /> To take into account the ring width variation, we take 80 evenly spaced snapshots of the system from 1000 to 1100 orbital periods of planet b and use the median of the 80 normalized ring widths as the ring width for the modeling. We start at 1000 orbital periods of planet b, adopted from Rodigas et al. (2014), and choose a time range of 100 orbital periods of period b because it is comparable to or longer than the synodic or resonance libration timescale for the planet mass range we study. The WHFast integrator is used with a time step of 1% of the planet's orbital period. For the planet mass range we study (0.3–10 $M_{Jup}$), a linear mass-normalized ring width relation is found:

$$ w = 0.019\mu/\mu_{Jup} + 0.082, \qquad (14) $$

For two-planet configurations, we adapt a similar setup to the edge and gap feature (Section 3.1): planet b is a Jupiter-mass planet with a semimajor axis of 16 au ($m_b = 1 M_{Jup}$ and $a_b = 16$ au); planet c is interior to planet b with a mass range of 0.0005–1 $M_{Jup}$ and a semimajor axis range of 1–15 au. We determine the stable outer ring locations and characterize the normalized ring widths following the procedures described for the single-planet case. The ring locations we find in these two-planet configurations are similar to the inner edge locations we find in Section 3.1, as expected.

We draw several conclusions about the normalized ring widths for two-planet configurations, which are summarized in Figure 17. First, if planet c is a low-mass planet with a small semimajor axis (i.e., the bottom-left corner), the normalized ring width is similar to the single-planet case. Second, the normalized ring width could be reduced by 8%–23% in systems with planet c with either a comparable mass or a comparable orbit to planet b but not both (i.e., configurations colored in white/light blue in Figure 17). However, for a Jupiter-mass planet c with a semimajor axis of 1 au, we find the ring width is further thickened compared to a single-planet model. The reason could be that the ring is placed closer to planet b for the configuration. Lastly, a planet c with mass and orbit comparable to planet b (i.e., configurations colored in dark blue in Figure 17) significantly thickens the normalized ring width. However, we notice one exceptional two-planet configuration that causes less spread of the ring compared to other

<footnote n="14">¹⁴ The variation follows the synodic timescale (i.e., the time taken for a conjunction), the repeated conjunction timescale (i.e., the time taken to repeat conjunctions at the same location), and the resonance libration timescale.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Normalized ring width sculpted by both planet b ($m_b$ = 1 $M_{Jup}$ and $a_b$ = 16 au) and planet c. The hatched region represents dynamically unstable two-planet configurations over 10 Myr from *N*-body simulations (i.e., $e_c$ > 1) and is color coded as the normalized ring width of 0.13, which corresponds to the width sculpted by planet b only. The ring width is reduced in some two-planet configurations (i.e., colored in white/light blue) compared to a single-planet system and thickened for a planet c with a comparable mass and orbit to planet b.</figcaption></figure>

two-planet systems with similar configurations ($m_c$ = 1 $M_{Jup}$ and $a_c$ = 12 au; the upper-right pixel in white). Future investigation is required to understand why the normalized ring width does not spread out as much in this configuration. The corresponding mass estimation errors for planet b in these two-planet configurations are presented in Figure 18. We estimate planet b's mass ($m'_b$) from Equation (14) and calculate the mass estimation error from Equation (11). Because of the large single-planet relationship intrinsic variation, we set the minimum mass estimation error percentage as 200%.<fnref n="15" /> We reach a conclusion similar to that for edge and gap features: an overestimate of the mass upper limit for planet b occurs when a planet c with mass and orbit comparable to planet b is present but undetected.

We also investigate how a planet c exterior to the ring affects the normalized ring width (e.g., Fomalhaut's ring; Boley et al. 2012). We consider the same planet b ($m_b$ = 1 $M_{Jup}$ and $a_b$ = 16 au) and a Jupiter-mass or Neptune-mass planet c. We tune the semimajor axis of planet c and slightly increase the semimajor axis of the ring if necessary to attain stable system configurations. For a Jupiter-mass planet c, we find that the normalized ring width could be thickened by ~50%, which corresponds to a mass estimation error of ~500% for planet b. In such configurations, planet c could even be tens of astronomical units apart from planet b but still have a strong gravitational effect on the ring width. For a Neptune-mass planet c, the normalized ring width could be either thickened or reduced by roughly 10%–20%, which may also induce a small mass estimation error for planet b. However, this error is comparable to the uncertainties of the single-planet relationship and may be considered negligible.

In summary, the thickened ring width is usually primarily dominated by a single planet. However, the mass estimation for the posited planet (i.e., an upper limit on the planet's mass) could be severely flawed only for (1) an inner, hidden planet with mass and orbit comparable to the posited planet or (2) an outer, wide-separation planet with mass comparable to the posited planet. Moreover, the intrinsic variations of the single-planet relationship already cause a large mass estimation error

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Mass upper-limit estimation error for planet b ($m_b$ = 1 $M_{Jup}$ and $a_b$ = 16 au) using the normalized ring width sculpted by both planet b and planet c. The hatched region represents dynamically unstable two-planet configurations over 10 Myr from *N*-body simulations (i.e., $e_c$ > 1). The color scale demonstrates the levels of mass estimation error labeled by the color bar. Scatter in the single-planet relationship can cause a large mass estimation error for planet b (~150%). The mass estimation error for planet b is severely flawed only for a planet c with mass and orbit comparable to planet b.</figcaption></figure>

for the posited planet and usually dominates over the error caused by an undetected planet.

## 4.2. Scale Heights

In Section 4.1, we investigated the radially thickened ring feature, in which planets radially stir nearby planetesimals to a wider radial range. Similarly, embedded planets can also vertically stir nearby material, resulting in substantial scale heights (e.g., AU Mic, β-Pic; Krist et al. 2005; Matrà et al. 2017; Daley et al. 2019). Previous studies on scale height focused on how small stirring bodies (e.g., a Pluto-sized planetary embryo) imparts random kinetic energy to smaller planetesimals inside or near its Hill sphere (e.g., Quillen et al. 2007; Thébault & Augereau 2007). The scale height-to-radius aspect ratio (i.e., $H/r$) can be used as a mass or size indicator of the largest stirring bodies. In this work, we focus on scale heights sculpted by much larger stirring bodies—planets. We study orbital changes of planetesimals (i.e., inclination changes) due to the close encounter with the sculpting planet and the resulting disk scale height.

Thébault (2009) raised concerns about degeneracy in the scale height feature caused by the radiation pressure on small dust grains (i.e., with a size of 10 µm or less). Fortunately, this degeneracy only affects observations at short wavelengths (i.e., comparable to the small dust grain sizes). For longer wavelength debris disk observations (i.e., $\lambda$ > 50 µm), the scale height will still be dominated by the stirring bodies.

We first investigate the scale height feature in a single-planet case. We are not aware of any analytical expressions that describe the scale height of planetesimals generated by a single, massive planet separated by several Hill radii from a ring. To understand how the planet's mass affects the ring's scale height, we use *N*-body simulations with a similar setup to the radially thickened ring feature. We consider a planet mass range of 0.3–10 $M_{Jup}$ and put a ring of test particles outside the sculpting planet's orbit with initial inclinations of 0.01 radian longitudes of the ascending nodes, and random mean anomalies. We study how a planet with zero inclination increases the disk aspect ratio from the initial value of 0.01. Although we find in general the scale height increases as the planet mass increases, there is not a monotonic relation

<footnote n="15">¹⁵ The mass estimation error for a Jupiter-mass planet b with the normalized ring width of 0.13 is around 150% using the mass-normalized ring width single-planet relation.</footnote>
between planet mass and scale height. The reason could be the planet's resonances affect the ring location, and the scale height is sensitive to the proximity of the planet. For example, for a 5 $M_{\text{Jup}}$ planet with a semimajor axis of 16 au, the scale height barely increases above the initial value of 0.01 because the ring needs to be 3.65 $R_{\text{Hill}}$ apart from the planet to avoid the planet's 5:3 resonance, which destabilizes the ring. The closest location for a stable ring is still too far from the planet to have its scale height effectively sculpted. For a 6 $M_{\text{Jup}}$ planet, however, the scale height is greatly increased because a stable ring can be placed at a semimajor axis similar to the 5 $M_{\text{Jup}}$ planet (also avoiding the 5:3 resonance) but only 3.35 $R_{\text{Hill}}$ apart from the planet. Due to the planet's resonances and the sensitivity of the scale height to the proximity of the planet, we cannot derive a simple relation between the disk scale height and the planet mass. When interpreting observations, if both the planet and the ring location are known, we could characterize the sculpting planet mass from the width of the chaotic zone (e.g., Equation (10)) and use that mass in numerical simulations to compare the simulated scale height to the observed one.

We further investigate how a second planet (i.e., planet c), in addition to the first planet (i.e., planet b), affects the scale height. If the scale height is increased by planet c, we are at risk of overestimating planet b's mass. In the single-planet configuration, planet b is a Jupiter-mass planet with a semimajor axis of 16 au ($m_b = 1 M_{\text{Jup}}$, $a_b = 16$ au), which is the same setup we used for the thickened ring feature (Section 4.1) and the edge-and-gap feature (Section 3.1). We place a zero-width ring of 5000 test particles as close to planet b as possible to have the strongest gravitational effect on the ring while keeping it stable ($a_{\text{ring}} = 20$ au). We integrate the system for 1000 orbit periods of planet b with a time step of 1% of planet b's orbital period using the WHFast integrator. In the left two panels of Figure 19, we plot the initial ring setup (i.e., a zero-width ring with an aspect ratio of 0.01) and the ring sculpted by the Jupiter-mass planet b only. In the planet b-only case, both the ring width and the scale height are moderately increased and the scale height is elevated from 0.01 to ~0.015.

We then study how a planet c in the same system as planet b further increases the scale height. Two types of two-planet configurations are investigated: (1) planet c is interior to planet b's orbit (i.e., two planets on the same side of the ring) and (2) planet c is exterior to the ring (i.e., two planets on the opposite sides of the ring). In both setups, we find configurations for which the scale height is significantly increased by adding in planet c, demonstrated in Figure 19. In the planet c-interior-to-planet b case, we explore planet c with a mass range of 0.0005–1 $M_{\text{Jup}}$ and a semimajor axis range of 1–15 au. Similar to the radially thickened ring feature, we take 80 evenly spaced snapshots of the system from 1000 to 1100 orbital periods of planet b and use the median of the 80 aspect ratios as the aspect ratio for the modeling. As shown in Figure 20, for a planet c more massive than 0.03 Jupiter mass ($m_c > 0.03 M_{\text{Jup}}$ and $a_c = 1–15$ au), the scale height could be moderately increased from 0.015 in the single-planet case to roughly 0.02–0.025 in the two-planet case. For a few configurations for which planet c is both massive and at a large semimajor axis, the scale height could be significantly increased to ~0.05.

We demonstrate one example in Figure 19 as “plt b + c (same side),” where $m_c = 0.46 M_{\text{Jup}}$ and $a_c = 8.8$ au. Planet b and the ring have the same setup as the single-planet configuration ($m_b = 1 M_{\text{Jup}}$, $a_b = 16$ au, and $a_{\text{ring}} = 20$ au).

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Comparison of a debris ring's aspect ratio ($H/r$) and normalized ring width sculpted by planet b only (“plt b-only”; $m_b = 1 M_{\text{Jup}}$, $a_b = 16$ au), by planet b with a planet c interior to planet b's orbit (“plt b+c (same side)”; $m_c = 0.46 M_{\text{Jup}}$, $a_c = 8.8$ au), and by planet b with a planet c exterior to the debris ring (“plt b+c (opposite sides)”; $m_c = 0.054 M_{\text{Jup}}$, $a_c = 32$ au) after 1000 orbit periods of planet b. “Initial” refers to the initial orbital parameters of the debris ring: 5000 test particles with a semimajor axis of 20 au and an inclination of 0.01 rad. The two blue dashed lines indicate an aspect ratio of ±0.01 for reference. An additional planet can significantly increase the scale height, compared to the scale height sculpted by planet b only.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Scale height-to-radius aspect ratio ($H/r$) sculpted by both planet b ($m_b = 1 M_{\text{Jup}}$ and $a_b = 16$ au) and planet c. The hatched filled region represents dynamically unstable two-planet configurations over 10 Myr from N-body simulations (i.e., $e_c > 1$) and is color coded relative to the aspect ratio of 0.015, which corresponds to the scale height sculpted by planet b only. The aspect ratio could be moderately increased from 0.015 to roughly 0.02–0.025 for a planet c more massive than 0.3 $M_{\text{Jup}}$.</figcaption></figure>

For this configuration, planet c is massive and close to planet b but does not increase the width of the chaotic zone so that the ring is still close to planet b. The scale height is significantly increased to nearly 0.05 due to the gravitational perturbation from both planets. However, if the ring is sculpted by either planet b or planet c only, the scale height is barely enhanced.

In the planet c-exterior-to-the-ring case, we consider a planet c of Jupiter mass or Neptune mass and a wide range of semimajor axes. In both cases, we find that planet c could increase the scale height from that sculpted by planet b. Specifically, a Neptune-mass planet c close to the ring could significantly increase the scale height. We demonstrate this effect in Figure 19, “plt b + c (opposite sides),” where $m_c = 0.054 M_{\text{Jup}}$ and $a_c = 32$ au. Planet b and the ring have the same setup as the single-planet configuration ($m_b = 1 M_{\text{Jup}}$, $a_b = 16$ au, and $a_{\text{ring}} = 20$ au). Compared to the scale height of ~0.015, if sculpted by planet b only, the scale height increases to ~0.04 by adding a Neptune-mass planet c. For a Jupiter-mass planet c, the scale height may also be increased but usually only by a moderate amount because planet c could not be placed close to the ring without destabilizing the system.

For both types of two-planet configurations we explored, the proximity of the planets to the disk is a key factor in increasing
disk scale height. The sculpting planets' masses play a
secondary role. The dependence on proximity explains why a
Neptune-mass planet c has a stronger effect on the disk scale
height than a Jupiter-mass planet c when two planets are on the
opposite sides of the debris ring. The Neptune-mass planet can
be closer to the ring without destabilizing it and therefore can
increase the scale height more. Similarly, a planet c interior to
planet b's orbit that does not increase the width of the chaotic
zone has a stronger effect on the scale height, compared to
those with higher masses and greater semimajor axes that will
increase the width of the chaotic zone.

In summary, the disk scale height can be used as a mass
indicator of the stirring planet. Due to the complex interplay
between Hill and resonant stability, there is no simple
relationship between planet mass and disk scale height.
Therefore, deriving a mass requires detailed modeling of the
individual system. Although the aspect ratio is usually
primarily dominated by the nearby planet, we find a second
stirring planet could affect the aspect ratio by a significant
amount if it is close to the disk. For these configurations, our
inference of the assumed planet's mass may be overestimated.

5. Conclusion

Resolved debris disk features contain valuable information about the underlying planetary systems, and can reveal the mass, semimajor axis, and other orbital parameters of the planets by which they are created. Most existing planet–disk models assumed a single planet sculpts the disk feature. We investigated if/how planet properties inferred from single-planet models are compromised when a second, unseen planet resides in the system. On the bright side, we found debris disk features (warps, offsets, edges and gaps, azimuthal asymmetries, thickened rings, and scale heights) are usually primarily dominated by a single planet. However, the inferred planet properties can be severely flawed if the wrong planet is identified as the source of the morphology. We summarize our findings for each feature and highlight situations where we can infer the need to use a multiplanet model instead of a single-planet one from the disk morphology and thus avoid potential planet characterization errors.

5.1. Warps

The warp feature (Section 2.1) allows for constraints on the sculpting planet's mass, semimajor axis, and inclination. If a planet is detected in the system and its semimajor axis is known, the planet's mass can be inferred from the warp location and the warp production time from Equation (1). We studied how our inference of the planet's mass could be compromised by an undetected planet present in the system. We found that mass estimations for the detected planet could be significantly flawed when the detected planet is not the dominant planet in the system, demonstrated in Figure 4. In extreme cases, the detected planet's mass could be misestimated by orders of magnitude. To avoid such estimation errors, we identified two types of disk morphology dissimilarities between a single-planet system and a two-planet system (Section 2.1.2): the thickness of the outer disk (Figure 5) and multiple warps (Figure 6). If the disk is well resolved, comparing the single-planet model to the observed disk image allows one to infer whether or not a multiplanet model is needed.

We also explored how an undetected planet could compromise our inference of the detected planet's inclination. We found that a warped disk's overall inclination follows the inclination of the more massive planet (Figures 8 and 9). The inference of the detected planets inclination could be significantly different from the true value if a hidden planet is more massive than the detected planet. If an observed system has no detected planet, we can infer the inclination of the most massive planet from disk inclination. If the system has any detected planets, ruling out the presence of any more massive planets is necessary before characterizing the detected planet's inclination from the disk.

5.2. Offsets

The offset feature (Section 2.2) is useful for inferring the sculpting planet's eccentricity. In the case where the observed ring lies just beyond the sculpting planet's chaotic zone, we can approximate the sculpting planet's eccentricity by the observed ring eccentricity (Equation (8)). We investigated two types of disk configurations: a debris ring and a debris disk. We found that the forced eccentricity of a debris ring is primarily sculpted by the nearest planet, demonstrated in Figure 11. If we can confirm the observed planet is the one sculpting the ring, we can infer the planet's eccentricity from the ring eccentricity.

We found it challenging to characterize the planet's
eccentricity from the overall disk eccentricity. First, the disk
eccentricity can be much smaller than the sculpting planet's
eccentricity (Equation (8)). Second, the forced eccentricity of
an entire disk in a multiplanet system can be complicated by
secular resonances. A two-planet system would introduce two
interior resonances and two exterior resonances that all affect
disk eccentricity. Fortunately, we may be able to identify a
multiplanet system from disk features sculpted by secular
resonances (e.g., multiple gaps; see Yelverton & Kennedy 2018
for a detailed discussion).

5.3. Edges and Gaps

Width and edge features (Section 3.1) place constraints on the sculpting planet's mass and semimajor axis. We confirm the relation between planet mass and the width of the chaotic zone, also known as the “2/7” law (Equation (10)). For two-planet configurations, we studied how a second, hidden planet would further increase the width of the chaotic zone and thus impact our inference of the detected planet's mass. We found the chaotic zone width is usually primarily dominated by a single planet. The mass estimation for the detected planet is flawed only for a hidden planet with mass and orbit comparable to the detected planet, as shown in Figure 12. The mass estimation error could at most be ~250% in such cases. Fortunately, in observed systems, we would expect two planets with similar masses and semimajor axes to be both detected (via either direct imaging or RV) if one is detected. We can be less concerned about the planet's mass estimation if we can rule out another, more massive planet closer to the ring.

The gap's edge sharpness is affected by the planet's mass (Quillen 2006), but it is difficult to constrain the planet's mass uniquely even in the single-planet case (Chiang et al. 2009). We further explored how a second planet modifies the edge sharpness and found edge sharpness dissimilarity cannot unambiguously distinguish between one-planet and two-planet systems, demonstrated in Figures 13 and 14.
We also investigated how would two planets in MMR, compared to two similarly spaced, similar-mass planets not in MMR, affect the disk morphology. We found that two planets in MMR do not significantly enhance the chaotic behavior of the planetesimals or modify the edge sharpness, as demonstrated in Figure 15.

### 5.4. Azimuthal Asymmetries

The azimuthal asymmetry feature (Section 3.2) places constraints on the underlying planet's mass, semimajor axis, and migration history. We studied how a second planet, migrating along with the planet primarily sculpting the disk, modifies the resonant behavior of planetesimals and thus affect the asymmetries. We found a second planet could affect the azimuthal asymmetries by destabilizing captured planetesimals crossing its orbit, as shown in Figure 16, panels (b) and (c). If no clumpy structure is observed in a system, we cannot rule out migration: the lack of resonant objects could be the result of multiple migrating planets. However, if clumpy structure is observed, it is likely to be sculpted by one dominant planet (e.g., Figure 16 panel (d)). Using a single-planet model to characterize that planet's properties and migration history is generally acceptable.

### 5.5. Thickened Rings

The thickened ring feature (Section 4.1) places an upper limit on the sculpting planet's mass, which is particularly useful for characterizing planet properties in systems with no detected planet. We confirmed that the planet mass-normalized ring width relation derived in Rodigas et al. (2014) for scattered-light images also works for planetesimal distributions for single-planet systems, described in Equation (14). We investigated two types of two-planet configurations: two planets on the same side of the ring and two planets on the opposite sides of the ring. We found that the thickened ring width is usually primarily dominated by the more massive planet closer to the ring. For configurations with both planets interior to the ring, the ring width could be thickened by 40%-70% if two planets have comparable mass and orbits (Figure 17). Using a single-planet model for such configurations, we will greatly overestimate the upper limit on the sculpting planet's mass by 400%-650% (Figure 18). For configurations with both planets on the opposite side of the ring, if two planets have comparable mass, the ring could also get significantly thickened even if one planet is far from the ring. Ignoring such hidden planets could cause a mass estimation error of ≲500%. For both types of two-planet configurations, a second planet could further thicken the ring width and lead to an incorrect inference on the assumed planet's maximum mass. When interpreting observations, we must be cautious using the planet mass-normalized ring width relation when we do not know how many planets may contribute to the observed ring feature. However, if the ring system we study has one or more detected planets likely sculpting the ring, we can be less worried because another planet can only compromise our inference if it has comparable mass and comparable/greater semimajor axis. If one is detected

via either direct imaging or RV, the other one is likely to be detected, too.

### 5.6. Scale Heights

The scale height feature (Section 4.2) is useful to characterize the mass of the sculpting planet. We examined how a second planet would further increase the disk scale height and thus lead to an overestimation for the assumed planet's mass. We explored two types of two-planet configurations: two planets on the same side of the debris ring (Figure 20) and two planets on opposite sides of the ring. We found the scale height is usually primarily dominated by the stirring planet with greater mass and closer to the ring. However, in both types of configurations, we found systems for which the scale height is significantly increased by the second planet, demonstrated in Figure 19. The proximity of the planets to the ring is the key factor that increases scale height.

### 5.7. Summary

Using our findings, we rank the risk levels of using a single-planet model to characterize planet properties for the disk features we studied. We are at high risk of misinterpreting planet properties from disk warps; at moderate risk from disk edges and gaps, radially thickened rings, and scale height features; and at low risk from host star-disk center offsets and azimuthal asymmetries.

We thank the referee for a helpful report that improved the clarity of the paper. We are grateful to Meredith Hughes for helpful discussion and feedback. We thank Ruobing Dong, Tom Esposito, Brad Foley, Eric Ford, Kevin Luhman, and Jason Wright for helpful comments. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. We gratefully acknowledge support from the Alfred P. Sloan Foundation's Sloan Research Fellowship and the Oak Ridge Associated Universities Ralph E. Powe Junior Faculty Enhancement Award. A.S. is supported by funding from the European Research Council under the European Community's H2020 (2014-2020/ERC Grant Agreement No. 669416 "LUCKY STAR"). Computations for this research were performed on the Pennsylvania State University's Institute for CyberScience Advanced CyberInfrastructure (ICS-ACI). This content is solely the responsibility of the authors and does not necessarily represent the views of the Institute for CyberScience. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound.

*Software: REBOUND (Rein & Liu 2012; Rein & Spiegel 2015; Rein & Tamayo 2015), REBOUNDx (Tamayo et al. 2019), Matplotlib (Hunter 2007; Droettboom et al. 2016), Numpy (van der Walt et al. 2011), Jupyter (Kluyver et al. 2016).*

## Appendix A
### Simulation Setup for Figure 1

Table 1 shows setup parameters for N-body simulations in Figure 1.
<table><caption>**Table 1**
Setup Parameters for N-body Simulations in Figure 1</caption><thead><tr><th>Debris Disk Feature</th><th>m<sub>b</sub><br>(M<sub>Jup</sub>)</th><th>a<sub>b</sub><br>(au)</th><th>i<sub>b</sub><br>(°)</th><th>e<sub>b</sub></th><th>a<sub>par</sub><br>(au)</th><th>N<sub>tot</sub></th><th>Integration Time<br>(Myr)</th><th>Integrator</th><th>Notes</th></tr></thead><tbody><tr><td>Warp</td><td>10</td><td>10</td><td>10</td><td>0</td><td>20–200</td><td>1e4</td><td>10</td><td>whfast</td><td>edge-on view</td></tr><tr><td>Spiral Arm</td><td>10</td><td>10</td><td>10</td><td>0</td><td>20–200</td><td>1e4</td><td>10</td><td>whfast</td><td>view at 30°-inc, 150°-pos Angle</td></tr><tr><td>Offset</td><td>1</td><td>10</td><td>0</td><td>0.4</td><td>20–200</td><td>1e4</td><td>10</td><td>whfast</td><td>...</td></tr><tr><td>Azimuthal Asymmetry</td><td>1</td><td>100</td><td>0</td><td>0</td><td>20–200</td><td>7.5e3</td><td>10</td><td>ias15</td><td>a<sub>b</sub> decaying as e<sup>-t/(20 Myr)</sup></td></tr><tr><td>Edge and Gap</td><td>1</td><td>16</td><td>0</td><td>0</td><td>16–26</td><td>4e3</td><td>10</td><td>ias15</td><td>...</td></tr><tr><td>Radially Thickened Ring</td><td>1</td><td>10</td><td>0</td><td>0</td><td>13</td><td>1e4</td><td>0.3</td><td>whfast</td><td>...</td></tr><tr><td>Scale Height</td><td>1</td><td>10</td><td>10</td><td>0</td><td>9–11</td><td>1e4</td><td>0.3</td><td>whfast</td><td>...</td></tr></tbody></table>

**Note.** Columns 2–5 illustrate planet b’s mass, semimajor axis, inclination, and eccentricity, respectively. Columns 6–7 present the semimajor axis range and the total number of test particles. Integrator types for the simulations are indicated in column 8. (a) $\Omega_b$ and $\bar{\omega}_b$ are fixed to be zero in all simulations. (b) All test particles are massless. (c) We assume a solar-mass host star. (d) N-body simulations run on REBOUND (Rein & Spiegel 2015; Rein & Tamayo 2015) and REBOUNDx (Tamayo et al. 2019).

## Appendix B
### Secular Perturbation Code

The secular perturbation code is based on the Laplace-Lagrange approximation, a first-order approximation in $e$ and $i$ of the secular terms in the disturbing function, which works appropriately for planets with low inclinations ($i \lesssim 20°$) and eccentricities ($e \lesssim 0.3$; e.g., Murray & Dermott 1999, Chapter 7). We describe a planetesimal’s inclination and longitude of ascending node using its inclination vector ($p, q$) as $i = \sqrt{p^2 + q^2}$, $\Omega = \tan^{-1}(p/q)$ and a planetesimal’s eccentricity and longitude of pericenter using its eccentricity vector ($h, k$) as $e = \sqrt{h^2 + k^2}$, $\bar{\omega} = \tan^{-1}(h/k)$. The inclination vector ($p, q$) and the eccentricity vector ($h, k$) have the following form:

$$ 
\begin{aligned}
p &= i_{\text{free}} \sin(Bt + \gamma) + p_{\text{forced}}(t), \\
q &= i_{\text{free}} \cos(Bt + \gamma) + q_{\text{forced}}(t), \\
h &= e_{\text{free}} \sin(At + \beta) + h_{\text{forced}}(t), \\
k &= e_{\text{free}} \cos(At + \beta) + k_{\text{forced}}(t).
\end{aligned} 
\quad (15) $$

The first term in both the inclination and eccentricity vector is the free element. The free inclination $i_{\text{free}}$, the phase offset $\gamma$, the free eccentricity $e_{\text{free}}$, and the phase offset $\beta$ depend on the initial condition of the planetesimal. Unless stated otherwise, we set initial orbital parameters of a planetesimal to 0 except its semimajor axis. The nodal regression rate $B$ and the pericenter precession rate $A$ for a planetesimal located at the semimajor axis $a$ have a similar form:

$$ B = -A = -n \frac{1}{4} \sum_{j=1}^{N_{\text{plt}}} \frac{m_j}{m_*} \alpha_j \bar{\alpha}_j b_{3/2}^{(1)}(\alpha_j), \quad (16) $$

where $n$ is the orbital frequency of the planetesimal, $N_{\text{plt}}$ is the total number of planets in the system, $m_*$ is the stellar mass, $\alpha_j$ and $\bar{\alpha}_j$ for the $j$th planet at the semimajor axis of $a_j$ follow the formula:

$$ \alpha_j, \bar{\alpha}_j = \begin{cases} a_j/a, & 1 & \text{if } a_j < a \\ a/a_j, & 1 & \text{if } a_j > a \end{cases} \quad (17) $$

and $b_{3/2}^{(1)}$ is the Laplace coefficient. The second term in the inclination/eccentricity vector is the forced element, ($p_{\text{forced}}, q_{\text{forced}}$) and ($h_{\text{forced}}, k_{\text{forced}}$), reflecting the secular perturbation

of a planetesimal from planets:

$$ 
\begin{aligned}
p_{\text{forced}} &= \sum_{i=1}^{N_{\text{plt}}} \frac{-\mu_i}{B - f_i} \sin(f_i t + \gamma_i), \\
q_{\text{forced}} &= \sum_{i=1}^{N_{\text{plt}}} \frac{-\mu_i}{B - f_i} \cos(f_i t + \gamma_i), \\
h_{\text{forced}} &= \sum_{i=1}^{N_{\text{plt}}} \frac{-\nu_i}{A - g_i} \sin(g_i t + \beta_i), \\
k_{\text{forced}} &= \sum_{i=1}^{N_{\text{plt}}} \frac{-\nu_i}{A - g_i} \cos(g_i t + \beta_i),
\end{aligned} 
\quad (18) $$

where $\mu_i, f_i, \gamma_i, \nu_i, g_i$, and $\beta_i$ are the $i$th mode of an $N_{\text{plt}}$-planet system. $\gamma_i$ and $\beta_i$ can be determined from the initial condition of the planets. $f_i$ and $g_i$ are the $i$th eigenvalue of matrices $B$ and $A$, written as

$$ 
\begin{aligned}
B_{jk} &= +n_j \frac{1}{4} \frac{m_k}{m_* + m_j} \alpha_{jk} \bar{\alpha}_{jk} b_{3/2}^{(1)}(\alpha_{jk}), & B_{jj} &= -B_{jk}; \\
A_{jk} &= -n_j \frac{1}{4} \frac{m_k}{m_* + m_j} \alpha_{jk} \bar{\alpha}_{jk} b_{3/2}^{(2)}(\alpha_{jk}), & A_{jj} &= -B_{jj},
\end{aligned} 
\quad (19) $$

where $j \neq k$, $n_j$ is the orbital frequency of the $j$th planet, $m_k$ and $m_j$ are the masses of the $k$th and $j$th planets. $\alpha_{jk}$ and $\bar{\alpha}_{jk}$ in the equation are unitless numbers written as

$$ \alpha_{jk}, \bar{\alpha}_{jk} = \begin{cases} a_k/a_j, & 1 & \text{if } a_j > a_k \\ a_j/a_k, & 1 & \text{if } a_j < a_k, \end{cases} \quad (20) $$

and $b_{3/2}^{(1)}$ and $b_{3/2}^{(2)}$ are the Laplace coefficients written as

$$ 
\begin{aligned}
b_{3/2}^{(1)}(\alpha) &= \frac{1}{\pi} \int_{0}^{2\pi} \frac{\cos \psi d\psi}{(1 - 2\alpha \cos \psi + \alpha^2)^{3/2}} \\
b_{3/2}^{(2)}(\alpha) &= \frac{1}{\pi} \int_{0}^{2\pi} \frac{\cos 2\psi d\psi}{(1 - 2\alpha \cos \psi + \alpha^2)^{3/2}}.
\end{aligned} 
\quad (21) $$

$\mu_i$ and $\nu_i$ in Equation (18) can be expressed as

$$ \mu_i = \sum_{j=1}^{N_{\text{plt}}} B_j I_{ji}, \quad \nu_i = \sum_{j=1}^{N_{\text{plt}}} A_j e_{ji}, \quad (22) $$

where $I_{ji}$ and $e_{ji}$ are the $i$th eigenvector of matrices described in Equation (19) normalized by the initial conditions of
planets, and

$$ 
\begin{aligned}
B_j &= +n \frac{m_j}{4 m_*} \alpha_j \bar{\alpha}_j b_{3/2}^{(1)}(\alpha_j), \\
A_j &= -n \frac{m_j}{4 m_*} \alpha_j \bar{\alpha}_j b_{3/2}^{(2)}(\alpha_j).
\end{aligned} 
\quad (23) $$

$\alpha_j$ here has the similar expression to Equation (17). Given the secular code written from the analytical solution above, we are able to explore a full range of configurations of a multiplanet model and study the planets' effects on the disk.

## ORCID iDs

Jiayin Dong https://orcid.org/0000-0002-3610-6953

Rebekah I. Dawson https://orcid.org/0000-0001-9677-1296

Andrew Shannon https://orcid.org/0000-0002-0711-4516

Sarah Morrison https://orcid.org/0000-0002-2432-833X

## References

Augereau, J. C., Nelson, R. P., Lagrange, A. M., Papaloizou, J. C. B., & Mouillet, D. 2001, A&A, 370, 447

Bastien, F. A., Stassun, K. G., & Pepper, J. 2014, ApJL, 788, L9

Beichman, C., Barrado, D., Belikov, R., et al. 2019, BAAS, 51, 58

Biller, B. A., Liu, M. C., Rice, K., et al. 2015, MNRAS, 450, 4446

Boley, A. C., Payne, M. J., Corder, S., et al. 2012, ApJL, 750, L21

Bonnefoy, M., Marleau, G. D., Galicher, R., et al. 2014, A&A, 567, L9

Borucki, W. J., Koch, D., Basri, G., et al. 2010, Sci, 327, 977

Bowler, B. P. 2016, PASP, 128, 102001

Brande, J., Barclay, T., Schlieder, J. E., Lopez, E. D., & Quintana, E. V. 2020, AJ, 159, 18

Brandl, B. R., Absil, O., Agócs, T., et al. 2018, Proc. SPIE, 10702, 107021U

Brandl, B. R., Feldt, M., Glasse, A., et al. 2014, Proc. SPIE, 9147, 914721

Chen, C., Ballering, N., Duchene, G., et al. 2019, BAAS, 51, 342

Chiang, E., Kite, E., Kalas, P., Graham, J. R., & Clampin, M. 2009, ApJ, 693, 734

Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531

Daley, C., Hughes, A. M., Carter, E. S., et al. 2019, ApJ, 875, 87

Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175

Dawson, R. I., Murray-Clay, R. A., & Fabrycky, D. C. 2011, ApJL, 743, L17

Debes, J., Choquet, E., Faramaz, V. C., et al. 2019, BAAS, 51, 566

Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Sci, 343, 1490

Droettboom, M., Hunter, J., Caswell, T. A., et al. 2016, Matplotlib: Matplotlib V1.5.1, vv1.5.1, Zenodo, doi:10.5281/zenodo.44579

Duncan, M., Quinn, T., & Tremaine, S. 1989, Icar, 82, 402

Fernandez, J. A., & Ip, W. H. 1984, Icar, 58, 109

Gladman, B. 1993, Icar, 106, 247

Howard, A. W., Marcy, G. W., Johnson, J. A., et al. 2010, Sci, 330, 653

Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

Hunter, J. D. 2007, CSE, 9, 90

Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Sci, 322, 1345

Kalas, P., Graham, J. R., & Clampin, M. 2005, Natur, 435, 1067

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in 20th Int. Conf. Electronic Publishing, Jupyter Notebooks—A Publishing Format for Reproducible Computational Workflows (Southampton: CED), 87

Konishi, M., Grady, C. A., Schneider, G., et al. 2016, ApJL, 818, L23

Krist, J. E., Ardila, D. R., Golimowski, D. A., et al. 2005, AJ, 129, 1008

Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Sci, 329, 57

Lagrange, A. M., De Bondt, K., Meunier, N., et al. 2012, A&A, 542, A18

Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21

Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019, NatAs, 3, 1135

Lee, E. J., & Chiang, E. 2016, ApJ, 827, 125

MacGregor, M. A., Matrà, L., Kalas, P., et al. 2017, ApJ, 842, 8

MacGregor, M. A., Wilner, D. J., Rosenfeld, K. A., et al. 2013, ApJL, 762, L21

Matrà, L., Dent, W. R. F., Wyatt, M. C., et al. 2017, MNRAS, 464, 1415

Matrà, L., Wyatt, M. C., Wilner, D. J., et al. 2019, AJ, 157, 135

Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780, 97

Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv:1109.2497

Miley, J. M., Panić, O., Wyatt, M., & Kennedy, G. M. 2018, A&A, 615, L10

Morrison, S., & Malhotra, R. 2015, ApJ, 799, 41

Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ. Press)

Murray, N., & Holman, M. 1997, AJ, 114, 1246

Mustill, A. J., & Wyatt, M. C. 2011, MNRAS, 413, 554

Nesvold, E. R., & Kuchner, M. J. 2015, ApJ, 815, 61

Nesvold, E. R., Naoz, S., & Fitzgerald, M. P. 2017, ApJL, 837, L6

Nesvold, E. R., Naoz, S., Vican, L., & Farr, W. M. 2016, ApJ, 826, 19

Petrovich, C. 2015, ApJ, 808, 120

Quanz, S. P., Crossfield, I., Meyer, M. R., Schmalzl, E., & Held, J. 2015, IJASB, 14, 279

Quillen, A. C. 2006, MNRAS, 372, L14

Quillen, A. C., Morbidelli, A., & Moore, A. 2007, MNRAS, 380, 1642

Reche, R., Beust, H., Augereau, J. C., & Absil, O. 2008, A&A, 480, 551

Rein, H., & Liu, S.-F. 2012, A&A, 537, A128

Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424

Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376

Ricci, L., Carpenter, J. M., Fu, B., et al. 2015, ApJ, 798, 124

Roberge, A., Fischer, D., & Peterson, B. 2019, BAAS, 51, 199

Rodigas, T. J., Malhotra, R., & Hinz, P. M. 2014, ApJ, 780, 65

Schlichting, H. E. 2014, ApJL, 795, L15

Shannon, A., Mustill, A. J., & Wyatt, M. 2015, MNRAS, 448, 684

Sibthorpe, B., Kennedy, G. M., Wyatt, M. C., et al. 2018, MNRAS, 475, 3046

Su, K., Ballering, N., Ertel, S., et al. 2019, arXiv:1903.10616

Su, K. Y. L., Rieke, G. H., Stapelfeldt, K. R., et al. 2009, ApJ, 705, 314

Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2019, MNRAS, in press

Telesco, C. M., Fisher, R. S., Piña, R. K., et al. 2000, ApJ, 530, 329

Thébault, P. 2009, A&A, 505, 1269

Thébault, P., & Augereau, J. C. 2007, A&A, 472, 169

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Wilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ, 855, 56

Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409

Wisdom, J. 1980, AJ, 85, 1122

Wyatt, M. C. 2003, ApJ, 598, 1321

Wyatt, M. C. 2008, ARA&A, 46, 339

Wyatt, M. C. 2018, Debris Disks: Probing Planet Formation (Berlin: Springer), 146

Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918

Yelverton, B., & Kennedy, G. M. 2018, MNRAS, 479, 2673