Formation of Gaps in Self-Gravitating Debris Disks by Secular Resonance in a Single-planet System I: A Simplified Model

Spatially resolved images of debris disks frequently reveal complex morphologies such as gaps, spirals, and warps. Most existing models for explaining such morphologies focus on the role of massive perturbers (i.e. planets, stellar companions), ignoring the gravitational effects of the disk itself. Here we investigate the secular interaction between an eccentric planet and a massive, external debris disk using a simple analytical model. Our framework accounts for both the gravitational coupling between the disk and the planet, as well as the disk self-gravity -- with the limitation that it ignores the non-axisymmetric component of the disk (self-)gravity. We find generally that even when the disk is less massive than the planet, the system may feature secular resonances within the disk (contrary to what may be naively expected), where planetesimal eccentricities get significantly excited. Given this outcome we propose that double-ringed debris disks, such as those around HD 107146 and HD 92945, could be the result of secular resonances with a yet-undetected planet interior to the disk. We characterize the dependence of the properties of the secular resonances (i.e. locations, timescales, and widths) on the planet and disk parameters, finding that the mechanism is robust provided the disk is massive enough. As an example, we apply our results to HD 107146 and find that this mechanism readily produces $\sim 20$ au wide non-axisymmetric gaps. Our results may be used to set constraints on the total mass of double-ringed debris disks. We demonstrate this for HD 206893, for which we infer a disk mass of $\approx 170$ Earth masses by considering perturbations from the known brown dwarf companion.


INTRODUCTION
Debris disks are ubiquitous around main sequence stars, with current detection rates of ∼ 20% in the Solar neighbourhood (Montesinos et al. 2016;Sibthorpe et al. 2018). They are optically thin, almost devoid of gas, and are believed to be composed of objects ranging from micron-sized dust grains up to kilometre-sized planetesimals. Since the dust grains are short-lived compared to the stellar age (e.g. Dominik & Decin 2003), their sustained presence requires a massive reservoir of large planetesimals continually supplying fresh dust via mutual collisions (Backman & Paresce 1993). Observed disks typically contain 0.01 − 1M ⊕ in mm/cm-sized grains (Wyatt et al. 2003;Holland et al. 2017) which, when extrapolated, yields masses of ∼ 1 − 100M ⊕ for the parent planetesimal population (e.g. Wyatt & Dent 2002;Greaves et al. 2005;Krivov & Wyatt 2020). The spatial distribution of these planetesimals is probed indirectly with observations at millimetre wavelengths, e.g. by ALMA. At such wave-fied, especially in view of observations suggesting that debris disks could contain tens of Earth masses in large planetesimals (Wyatt & Dent 2002;Greaves et al. 2005;Krivov & Wyatt 2020). In this regard, Jalali & Tremaine (2012) have argued that many of observed debris disk features could be ascribed to the slow (m = 1, 2) modes which, if and when excited (e.g. by stellar flybys), could be supported by the disk gravity alone. Despite this fact, gravitational effects of debris disks have not yet been widely appreciated in the literature.
In this paper (the first in a series) we investigate the interaction between an eccentric planet and an external, massive debris disk. The primary aim of this work is to present a novel pathway to sculpting gaps, i.e. depleted regions, in broad debris discs.

Existing mechanisms and this work
To date, four debris disks are known to exhibit double-belt structures that are separated by depleted gaps in their dust distribution as traced by ALMA: HD 107146 (Ricci et al. 2015;Marino et al. 2018), HD 92945 (Marino et al. 2019), HD 15115 (MacGregor et al. 2019, and HD 206893 (Marino et al. 2020). These systems (except HD 206893) have no known companions or planets to date, and the disks are gaspoor. In this work we focus on the nearly face-on disk of HD 107146, a nearby ∼80-200 Myr old G2V star (Williams et al. 2004). This disk, extending from ∼30 au to ∼150 au, features a circular ∼40 au wide gap centred at around 70−80 au in which the continuum emission drops by ∼ 50% (Ricci et al. 2015;Marino et al. 2018).
Various mechanisms have been explored for explaining the origin of gaps in debris disks. In analogy with the asteroid and Kuiper belts, the most popular scenario involves the presence of single or multiple planets orbiting within the depleted region (Schüppler et al. 2016;Shannon et al. 2016;Zheng et al. 2017). For instance, it has been suggested that a planet of few tens of Earth masses on a near-circular orbit at ∼ 70 − 80 au could reproduce HD 107146's gap (Ricci et al. 2015;Marino et al. 2018). However, forming (in situ) planets at such large distances represents a challenge to the current paradigm of planet formation (e.g. Kenyon & Bromley 2015).
Other scenarios involving planets interior to the disk have also been considered, avoiding the aferomentioned issue. For instance, Tabeshian & Wiegert (2016) showed that a loweccentricity planet can carve a gap at its external 2:1 mean motion resonance. On the other hand, Pearce & Wyatt (2015) demonstrated that HD 107146-like disks could be produced as a result of secular interactions and scattering events between a massive (∼ 10 − 100M ⊕ ) planetesimal disk and an initially high-eccentricity (∼ 0.5) planet of comparable mass to the disk. In the course of evolution, the planetary orbit is then circularized due to scattering events. However, Pearce & Wyatt (2015) consider only the back reaction of the disk on the planet (and vice versa) in their simulations, neglecting the disk self-gravity.
Finally, Yelverton & Kennedy (2018) considered a scenario whereby two coplanar planets carve a gap through their secular resonances within an external debris disk, which was assumed to be massless. In their model, the secular resonances occur at sites where the precession rates of the planets (i.e. system's eigenfrequencies) match that of the planetesimals in the disk (due to planetary perturbations). They find that at and around one of the two resonant sites planetesimal eccentricities are excited, triggering a depletion in the disk surface density of the kind seen in HD 107146.
The model proposed by Yelverton & Kennedy (2018) requires (at least) two planets to ensure that their orbits are precessing due to planet-planet interactions, a condition necessary for establishing secular resonances. However, another mechanism which may drive planetary precession is the secular perturbation due to the disk, which was ignored by Yelverton & . This motivates our investigation into whether gaps could be carved in self-gravitating debris disks via secular resonances when perturbed by single rather than multiple inner planets. A related scenario was studied by Zheng et al. (2017) which showed that a single planet embedded within a decaying gaseous disk (i.e. transitional disk) could carve a wide gap around its orbit via sweeping secular resonances assisted by the waning disk gravity.
In this paper we propose that double-ringed structuresakin to that of HD 107146 -could be explained as the aftermath of secular resonances in systems hosting a single eccentric planet and an external self-gravitating debris disk. The mechanism we invoke here is different from those of Pearce & Wyatt (2015) and Yelverton & Kennedy (2018). It is realized through a secular resonance between the apsidal precession rate of planetesimals due to both the disk and planet, and that of the planet due to disk gravity (c.f. Yelverton & Kennedy 2018). Additionally, our mechanism does not require scattering events between the planet and disk particles (c.f. Pearce & Wyatt 2015). As we show below, this mechanism is robust over a wide range of parameters; particularly when the disk is less massive than the planet.
Our work is organized as follows. In Section 2 we describe our model system and present the equations governing planetesimal dynamics. In Section 3 we characterize the features of the secular resonances over a wide range of parameter space. In Section 4 we apply these considerations to HD 107146, and identify the planet-disk parameters which could reproduce the observed gap. In Section 5, using some of these parameters, we investigate the evolution of disk-planet systems and present our main results. We discuss our results along with their implications in Section 6, where we also consider the application of our results to other systems (HD 92945 and HD 206893). In Section 7 we critically assess the limitations of our model, discuss the implications of relaxing some of them, and propose future work. Our findings are summarized in Section 8.

ANALYTICAL MODEL
We describe a simple model to analyze the long-term dynamical evolution of planetesimals embedded within a massive debris disk in a single-planet system. In our notation, a planetesimal orbit is characterized by its semimajor axis a, eccentricity e, and longitude of pericenter . Orbital elements subscripted with 'p' and 'd' refer to the planet and the disk, respectively.

Model system
Our model system consists of a broad debris disk of mass M d orbiting the host star M c exterior to, and co-planar with, a planet of mass m p (M d , m p M c ). We assume the planet is initially on a low eccentricity orbit (e p ≤ 0.1) and that it does not intersect the disk along its orbit. We consider the debris disk to be razor-thin and initially axisymmetric. The disk surface density is characterized with a truncated powerlaw profile given by for a in ≤ a ≤ a out , and Σ d (a) = 0 elsewhere. Here, a in and a out are the semimajor axes of the inner and outer disk edges, respectively. Defining δ ≡ a out /a in > 1, the total mass M d of such a disk can be written as which allows us to express Σ d in terms of M d . This setup is very similar to that explored in Rafikov (2013) and Silsbee & Rafikov (2015a) in the context of planetesimal dynamics in circumbinary disks. In this work, unless otherwise stated, we adopt a fiducial disk model with p = 1, a in = 30 au and a out = 150 au (i.e. δ = 5). This choice of p corresponds to a disk with constant amount of mass per unit semimajor axis.

Secular gravitational effects
We are primarily interested in the long-term dynamics of large (∼km-sized) planetesimals. Since the latter are effectively insensitive to radiative non-gravitational forces, we focus purely on gravitational effects accounting for perturbations due to both (1) the debris disk and (2) the planet. For simplicity, the non-axisymmetric component of the disk gravity is ignored in this work, although, as we will see later, the disk naturally develops non-axisymmetry (a discussion of the implications of this omission is provided in §7.1.2). We perform calculations within the framework of secular (orbitaveraged) perturbation theory to second order in eccentricities (Murray & Dermott 1999).

Effects of the disk and planet on planetesimals
The secular dynamics of planetesimals is described by the disturbing function R which consists of contributions due to the planet R p and due to the disk R d . An analytic expression for the disturbing function R d due to an axisymmetric disk with surface density of the form (1) has been previously derived in Silsbee & Rafikov (2015b) (see also Heppenheimer 1980;Ward 1981;Sefilian & Touma 2019). Combining R d with the contribution R p due to the planet (e.g. Murray & Dermott 1999, equation 7.7), the total disturbing function R = R d + R p to second order in eccentricities reads as: where n = GM c /a 3 is the planetesimal mean motion and the meaning of different constants is explained below. In Equation (3), A = A d + A p is the precession rate of the free eccentricity vector of a planetesimal. It has contributions from both the gravity of the disk (A d ) and the planet (A p ). The contribution of the planet is (Murray & Dermott 1999) where a out,150 ≡ a out /(150 au), and the numerical estimate is for p = 1 and δ 1 such that ψ 1 ≈ −0.5. In general, the coefficient ψ 1 in Eq. (6) depends on the power-law index p as well as the planetesimal semimajor axis with respect to the disk edges (Silsbee & Rafikov 2015b, equation A33). As the sharp edges of the disk are approached, ψ 1 formally diverges. However, when the planetesimal is well separated from the edges (i.e. a in a a out ), ψ 1 is effectively a constant of order unity (depending on p) which can be well approximated by equation (A37) in Silsbee & Rafikov (2015b). It is very important to note that the disk and the planet drive planetesimal precession in opposite directions, A p > 0 and A d < 0, with A p (a) falling off more rapidly with a than |A d (a)|.
The term B p in Eq. (3) represents the excitation of planetesimal eccentricity due to the non-axisymmetric component of the planetary potential. It is given by (Murray & Dermott 1999) 3/2 (a p /a)e p .
Note that the analogous term due to the disk is absent in Eq.
(3), since we have neglected the non-axisymmetric component of the disk self-gravity.

Effect of the disk on planet
Next we consider the effect of the disk on the planet. Since the disk is taken to be axisymmetric, it simply causes the planetary apsidal angle to advance linearly in time such that p (t) = A d,p t+ p (0), i.e.˙ p = A d,p , without exchanging its angular momentum with the planet. In this work, without loss of generality, we set p (0) = 0. In Appendix A we show that the planetary precession rate A d,p due to the disk with surface density (1) is given by (see also, Petrovich et al. 2019): where n p = GM c /a 3 p is the planetary mean motion, a in,30 ≡ a in /(30 au), and the numerical estimate is for p = 1 and a p = 20 au such that φ c 1 ≈ 1.8. Here φ c 1 = φ c 1 (a p /a in , p, δ) is a factor of order unity accounting for contributions of the disk annuli close to the planet (Eq. A7). Its behavior as a function of a p /a in and for various disk models (i.e. p, δ) is shown in Fig. 13. For a p /a in 1, we have φ c 1 ≈ 1 regardless of (p, δ).

Combined planet-disk effects
The fact that the planet is precessing renders the forcing term in R (Eq. 3) time-dependent. This time dependence could be eliminated upon transferring to a frame precessing with the planetary orbit: i.e. by subtracting ΦA d,p from Eq.
(3) where Φ = na 2 1 − √ 1 − e 2 ≈ na 2 e 2 /2 is the action conjugate to the angle ∆ ≡ − p . As a result, we obtain the following expression: This completes our development of the disturbing function. Note that for the particular set of parameters in equations (4), (6), (8), the planetesimal free precession rate A at a = 70 au is comparable to that of the planetary orbit, A d,p . In Figure 1 we show the radial behavior of A = A d + A p , together with the curve for A d,p . The fact that A(a) = A d,p at certain semimajor axes has very important implications for planetesimal dynamics; see Section 2.4.

Evolution equations and their solution
The secular evolution of a planetesimal orbit in the combined potential of the planet and the disk can be determined by Lagrange's planetary equations (Murray & Dermott 1999). Introducing the eccentricity vector e = (K, H) = e(cos ∆ , sin ∆ ), convenient for describing the dynamics in the frame corotating with the planet (e.g. Heppenheimer . Planetesimal free precession rate A = A d + Ap due to both the planet and the disk as a function of semimajor axis (red curve). Dotted and dashed curves represent Ap(a) and A d (a), respectively. The blue line represents the rate of planetary precession A d,p due to the disk. Calculations assume a 20M⊕ disk with p = 1 extending from ain = 30 au to aout = 150 au, and a 0.6MJ planet at ap = 20 au around a 1.09M star (Model A, Table 1). Note that A(a) = A d,p at two locations: at 70 au and at ain. 1980), we find that: Note that in the case of a massless disk (A d,p = 0, A = A p ), one recovers the evolution equations due to a non-precessing perturbing planet (e.g. Murray & Dermott 1999). The system of equations (10) admits a general solution given by the superposition of the 'free' and 'forced' eccentricity vectors, e(t) = e free (t) + e forced (t) (Murray & Dermott 1999). In particular, when planetesimals are initiated on circular orbits, K(0) = H(0) = 0, we have e free = e forced and the evolution of planetesimal orbits is described by: where the forced eccentricity is given by (13) Equations (11)-(13) represent the key solutions needed for our work. We remark that this framework has been previously verified against direct orbit integrations of test particles in disks (e.g. Silsbee & Rafikov 2015b;Fontana & Marzari 2016;Davydenkova & Rafikov 2018).
For illustrative purposes, in Figure 2 we show the radial profiles of instantaneous eccentricities (left panels) and longitudes of pericenter (relative to the planet, right panels) of planetesimals computed using Eqs. (11) and (12) (i.e. for e(0) = 0) at different times, as indicated in each panel. The calculations assume the same disk-planet parameters as in Fig. 1 and we have taken e p = 0.05 -the parameters of the fiducial disk-planet model (Model A, Table 1) which we consider in details later in this work (Section 5). Furthermore, here we have sampled secular evolution using N = 5000 planetesimals with semimajor axes distributed logarithmically between a in and a out , i.e. with a ratio of spacing β = (a out /a in ) 1/N ≈ 1.0003, each of which is represented by a blue dot in Fig. 2. We note that, as is typical for secular evolution, the eccentricity oscillation at a given semimajor axis is bounded between the initial value of 0 and e m (a) = 2|e forced (a)| (the red lines in left panels of Fig. 2). Moreover, as expected, the period of each eccentricity oscillation in the frame corotating with the planet is given by

Planetesimal eccentricity behavior and secular resonances
We now describe the essential features of planetesimal dynamics in the combined disk-planet potential 1 . In general, planetesimal orbits evolve differently depending on their free precession rate A(a) relative to that of the planet A d,p , i.e. for A(a) > A d,p or A(a) < A d,p -see Eqs. (11), (12).
For the particular set of parameters in Figs. 1 and 2, we see that the regime A(a) > A d,p is realized at small separations from the planet, where the precession rate of planetesimals is dominated by the planet so that A ≈ A p (except near a in where A d diverges due to disk edge effects, Silsbee & Rafikov (2015b)); see also Eqs. (4), (6). In this planetdominated regime planetesimal orbits precess in the same direction as the planet (i.e. prograde, see Eq. 12 and right panels of Fig. 2), and we have e forced > 0 (Eq. 13). Thus, as planetesimal orbits evolve, the apsidal angles ∆ remain constrained within [−π/2, π/2] at all times. Moreover, planetesimals attain their maximum eccentricity when their orbits are aligned with that of the planet, i.e. when ∆ = 0; see Eq. (11) and Fig. 2. Assuming A p A d,p , the maximum planetesimal eccentricity in this regime is e m,p ≈ |2e forced,p | with (e.g. Murray & Dermott 1999) ≈ 1.8 × 10 −2 a p,20 a 70 e p 0.05 , 1 For detailed summary of the dynamics in an analogous setup (in application to planetesimal dynamics in circumbinary disks), see Rafikov (2013) and Silsbee & Rafikov (2015a). Planetesimal eccentricities (left panels) and apsidal angles (right panels, measured relative to that of the precessing planet) as a function of semimajor axis after t = 1, 10, 30, 50, 80 and 100 Myr of evolution (top to bottom). The time is also indicated relative to τ , Eq. (16). The planetesimals were initiated on circular orbits in the fiducial disk-planet model (Model A, Table 1). The maximum of eccentricity oscillations em = 2|e forced | (Eq. 13) is shown by the red lines. For reference, the solid black lines show the maximum planetesimal eccentricities driven by the planet in the absence of the disk (em,p, Eq. 14). The dashed vertical lines show the secular resonance location (ares = 70 au), where eccentricities diverge in the course of evolution. One can clearly see that at the resonance ∆ = −π/2 at all times. Note also the resonance near the disk inner edge. This figure is available as an animation in the electronic edition of the journal. see Eq. (13), where we have used the approximations b 3/2 (α) ≈ 3α and b 3/2 (α) ≈ (15/4)α 2 valid for small α. This is the limit of a massless disk, a configuration most often adopted in studies of debris disks. In the course of evolution, planetesimals in this regime will form an eccentric structure largely aligned with the planetary orbit (e.g. Wyatt et al. 1999).
In the opposite disk-dominated limit, far from the planet (and for a ≈ a in , which we discuss later), Figure 1 shows that the precession rate of planetesimals is dominated by the disk so that A ≈ −|A d | A d,p . In this regime planetesimal orbits undergo retrograde free precession (see Eq. 12 and right panels of Fig. 2), and we have e forced < 0. Thus, the apsidal angles ∆ are confined within the range ±[π/2, π] at all times. Moreover, planetesimals attain their maximum eccentricity when their orbits are anti-aligned with the planetary orbit, i.e. when |∆ | = π; see Eq. (11). Assuming A d,p → 0 for simplicity, the maximum eccentricity in this regime is e m,d ≈ |2e forced,d | with where the numerical estimate assumes p = 1 and a in a a out so that ψ 1 ≈ −0.5. Equation (15) shows that planetesimal eccentricities in the disk-dominated regime decline more rapidly with a than in the planet-dominated regime, and their magnitude is suppressed -an effect pointed out in Rafikov (2013). In the course of evolution, planetesimals in this regime will form an eccentric structure anti-aligned with the planetary orbit.

Main secular resonance
More importantly, one can clearly see that the transition between planet-and disk-dominated regimes occurs via a secular eccentricity resonance where A(a) = A d,p ; see Fig. 1 (see also Rafikov 2013;Silsbee & Rafikov 2015a). This resonance emerges because the relative precession between the planetesimal orbits and the planetary orbit vanishes, while the torque exerted by the non-axisymmetric component of the planet is non-zero. At and around the locations of secular resonances, a = a res , planetesimal eccentricities are forced to arbitrarily large values (in linear approximation), see left panels of Fig. 2. This is because the denominator in Eq. (13) becomes small, introducing a singularity into the secular solution 2 (Rafikov 2013). By taking a limit A(a res ) → A d,p in Eq. (13) we find that the growth of eccentricity at the resonance occurs linearly in time, e(t) = t/τ , with a characteristic timescale given by where the approximation is valid for a p a res . Eq. (16) also explains why the eccentricities at the resonance near the disk inner edge are pumped up more quickly than at the resonance at 70 au, see left panels of Fig. 2.
Moreover, we can see from the right panels of Fig. 2 that at the resonance ∆ remains fixed at −π/2, as expected from  Eq. (12). In Section 3.1 we will show that such secular resonances are generic: they occur for a large range of disk-toplanet mass ratios, 10 −4 M d /m p 2, for all a p a in . To further illustrate the analysis above, Figure 3 shows the radial profiles of planetesimal forced eccentricities computed for different values of disk mass. The calculations are done for the same planetary parameters as in Figs. 1, 2. The most pronounced feature in Fig. 3 is the occurrence of a secular resonance within the disk (apart from the one very close to a in , see below) for 10 −3 ≤ M d /m p ≤ 1, where e forced diverges. At the same time, e forced asymptotically approaches e forced,p inward of the resonance, i.e. where A A d,p , whereas e forced → e forced,d external to it, i.e. where A A d,p (which is, of course, possible only if a in a res a out ). At the highest disk mass, M d /m p = 2, there are no secular resonances as the disk dominates planetesimal precession throughout the whole disk.
We note that in the region where the dynamics is dominated by the disk e forced (a) does not follow the simple power law profile ∝ a −4 given by Eq. (15). By and large, this is because the disk edge effects neglected in computing Eq. (15) render ψ 1 = ψ 1 (a) in a non-trivial manner, even when a in a a out (Silsbee & Rafikov 2015b). For instance, it is evident in Fig. 1 that A d (a) behaves more like a constant for a in a a out rather than as A d ∝ a −1/2 (Eq. 6), implying that |ψ 1 | ∝ a 1/2 for the employed disk model. This will be important in §3.1. As a matter of fact, ψ 1 becomes independent of semimajor axis only in disks of infinite radial extent (Silsbee & Rafikov 2015b), whereas the radial range of our adopted disk is finite with δ = a out /a in = 5 ( §2.1) .

Secular resonance at ain
Finally, we clarify that the origin of the resonance at ≈ a in (apart from the one at a in ) lies in the fact that A d ∝ −|ψ 1 | diverges as the sharp edges of a razor-thin disk are approached, see black dashed lines in Fig. 1. This makes |A d (a)| ∼ A p (a) as a → a in , even for a modest value of disk mass. However, it is also known that disks with Σ d dropping continuously near the edges rather than discontinuously, or disks with small but non-zero thickness, should exhibit finite A d near the edges (Davydenkova & Rafikov 2018;Sefilian & Rafikov 2019); different from our disk model. Thus, in such more realistic disks, only a single resonance -rather than two -will occur. This is portrayed in Fig. 3 for M d /m p = 1 by artificially stipulating ψ 1 (a) = −0.5, i.e. by ignoring the edge effects (Silsbee & Rafikov 2015b).

Secular resonances and gaps in debris disks
To summarize, the analysis presented here elucidates that the disk gravity can have a considerable impact on the secular evolution of planetesimals. In the remainder of this paper, we exploit the feasibility of the discussed secular resonance as the basis of a mechanism for sculpting depleted regions, i.e. gaps, in debris disks.
The emergence of a gap could be understood as follows. Planetesimals on eccentric orbits spend most of their time near their apocenter, further away from their orbital semimajor axes. Thus, provided that a secular resonance occurs within the disk, we expect the surface density of planetesimals to be depleted around the resonance location where planetesimal eccentricities grow without bound. This reasoning, in essence, is similar to that presented by Yelverton & Kennedy (2018) where the authors show that two planets could carve a gap in an external massless debris disk through their secular resonances. Additionally, given that generally planetesimals in the inner disk parts tend to apsidally align with the planet while those in the outer parts tend to anti-align, we expect the depleted region to have a nonaxisymmetric shape. This effect has been previously pointed out by Pearce & Wyatt (2015) in the context of secular interaction between a debris disk and an interior, precessing planet.

CHARACTERIZATION OF SECULAR RESONANCES
We now investigate how the characteristics of the secular resonances -i.e. their locations, their associated timescales for exciting eccentricities, and their widths -depend on the properties of the disk and the planet. This will guide us in putting constraints on the possible disk-planet parameters that could reproduce the structure of an observed debris disk featuring a gap ( §4).

Location of secular resonances
As mentioned in Section 2.4, secular resonances occur at semimajor axes a = a res where the apsidal precession rates of both the planet and planetesimals are commensurate, Using Equations (4), (6) and (8), we can express the resonance condition (17) in terms of the disk-to-planet mass ratio M d /m p and the relevant semimajor axes, i.e. a res , a p , and a out , scaled by a in : Here C 1 = (2−p)/(δ 2−p −1) and C 2 = C 1 (1−δ −p−1 )/(p+ 1) are constants depending on the disk model. It follows from Eq. (18) that the locations of secular resonances can be computed relative to the disk inner edge as functions of a p /a in and M d /m p . This is illustrated in Figure 4, where we plot the contours of M d /m p in the (a p /a in , a res /a in ) plane computed using our fiducial disk model i.e. p = 1 and δ = 5 ( §2.1). Figure 4 shows that for any given planet, two or no secular resonances occur within the disk provided that 10 −4 M d /m p 2. Additionally, we can see that for any a p /a in one of the resonances always occurs in the vicinity of the disk inner edge as described in §2.4.2, i.e. a res,1 a in , and its location varies weakly with M d /m p . On the other hand, the second resonance occurs at semimajor axis a res,2 a res,1 whose location changes significantly with varying M d /m p . Indeed, with increasing M d /m p (at fixed a p /a in ) this resonance is pushed inwards from a out towards the inner resonance at a in until both resonances 'merge', i.e. the distance between them approaches zero. Figure 3 provides a complementary view of this behavior. Looking at Fig. 4 we also see that, for planets closer to the disk, larger M d /m p is necessary to maintain the resonance at a given semimajor axis.
We recall that the existence of the inner resonance is mainly due to the disk edge effects. That is, the divergence of A d (a) ∝ −|ψ 1 (a)| as a → a in allows the resonance condition (17) to be satisfied around ≈ a in , even for relatively small values of M d (Section 2.4). This explains why for a given a p /a in the resonance at a res,1 is constrained to be very close to a in irrespective of M d /m p . In the absence of edge effects this inner resonance will not exist, resulting in a single resonance for fixed system parameters rather than two. This is illustrated in Fig. 4 for M d /m p = 1 by setting ψ 1 (a) = −0.5 (white full line).
The behavior of the resonance locations can be explained analytically. Consider the approximate form of the resonance condition, Eq. (17), in the limit of a p /a in → 0 so that A d,p is negligible and one can use the asymptotic limit of b (1) 3/2 , and the two terms on the left hand side of Eq. (18) balance each other (recall that ψ 1 < 0). It is then easy to demonstrate that for a resonance to occur at a in a res a out , the disk mass must be given by where the numerical estimate is obtained for our fiducial disk model (p = 1, δ = 5), for which |ψ 1 (a)| ∝ a 1/2 when a in a a out , see Section 2.4 3 . Fixing M d /m p in Eq. (19) then approximates the slopes of the contours in Fig. 4 reasonably well -see the white dashed line. As expected, the numerical results deviate from the scaling in Eq. (19) both as a res → a in or a out , where ψ 1 diverges, and as a p → a in , since A d,p becomes non-negligible.

Timescale for eccentricity excitation
We now consider how the eccentricity excitation timescale varies as a function of model parameters. To this end, we make use of the definition of τ given by Eq. (16), which quantifies the time it takes for initially circular orbits to reach e = 1 at the resonance. We note that τ is a strong function of the resonance location, and it explicitly depends on the parameters of the planet but not the disk. This is because the disk, assumed to be axisymmetric in our model (Section 2), does not contribute to eccentricity excitation. In Figure 5 we plot the contours of τ in the (a p /a in , a res /a in ) plane for a particular choice of planetary mass and eccentricity, m p = 100M ⊕ and e p = 0.1, assuming a solar-mass star. It is evident that the timescales are shorter when the planet and the resonance location are closer together, i.e. in the lower-right corner of parameter space where a res /a p → 1. Note that for the adopted planetary parameters, over a broad range of parameter space the timescales range from ∼ 10 Myr to few Gyr; this is comparable to the ages of observed debris disks. Moreover, the slopes of the contours in Fig. 5 can be explained by setting τ to a constant in Eq. (16): this yields the scaling a res ∝ a 2/3 p illustrated by the white dashed-line in Fig. 5.
Finally, Equation (16) shows that τ is inversely proportional to both the planetary mass and eccentricity. Thus, more massive or eccentric planets exert larger torque and excite planetesimal eccentricities more quickly, shortening the timescale τ when a p /a in and a res /a in are kept fixed. This means that in Figure 5 the contours of τ will be shifted to the left (right) when the product of m p and e p is increased (decreased).

Resonance width
We now quantify the range of semimajor axes w over which resonances act to significantly excite planetesimal ec-centricities. To this end, we follow 4 Yelverton & Kennedy (2018) and calculate the distance over which the forced planetesimal eccentricities e forced (a) exceed a constant threshold valueẽ. That is, we define w as the difference (in absolute values) between the two values of semimajor axis a i (i = 1, 2) satisfying in the vicinity of a given resonance. Here, we clarify that this definition serves as a proxy for the significance of a given resonance, and it does not necessarily correspond to the actual widths of gaps that we expect to observe 5 . In Equation (20), the planetary and disk masses appear only through their ratio M d /m p , and the two relevant semimajor axes -a i and a p -could be expressed relative to a in ; see Eqs. (4) -(8). Furthermore, the ratio M d /m p could be related to a p /a in and a res /a in by using the condition for secular resonance, Eqs. (17), (18). Thus, we can compute the resonance width w relative to a in as functions of a p /a in and a res /a in only, onceẽ and e p are specified (recall that B p ∝ e p , Eq. 7).
The threshold eccentricityẽ in Eq. (20) represents an ad hoc parameter, necessitating a physical justification for a particular choice of its value. To this end, we note that the presence of a physical gap within the disk is subject to the condition that planetesimal eccentricities are larger around the resonances than elsewhere. Away from the resonances, the forced planetesimal eccentricity is maximized near the disk inner edge where, approximately, e forced (a in ) → e forced,p (a in ) which can not exceed e p ; see Eq. (14), Fig.  3. Based on this reasoning we adoptẽ = e p in what follows, unless stated otherwise.
In Figure 6 we plot the contours of w/a in in the (a p /a in , a res /a in ) plane for our fiducial disk model with p = 1 and δ = 5 (see §2.1), assumingẽ = e p . Looking at Figure 6, we see that increasing the planetary semimajor axis for a fixed a in tends to generally broaden the width of a given resonance. This is, though, less obvious in the range 1.1 a res /a in 1.5 as the width there is a weaker function of a p /a in . Secondly (and relatedly), we see that for a given planetary semimajor axis, resonances occurring closer to the disk inner edge generally have larger widths compared to resonances further away; see also To understand this behavior, we recall that for a given a p /a in , our disk model with sharp edges has two resonance sites: one always at a res,1 a in and another further away at a res,1 a res,2 a out ; see Section 3.1. In terms of a res,2 , the inner resonance will be much narrower than the other. This behavior could be understood for instance by looking at the curves in Fig. 3 for M d /m p = 10 −3 , 10 −2 or 10 −1 , which show that the inner resonance width is insignificant.
On the other hand, for fixed (a p /a in , M d /m p ), if the resonances are close to each other such that a res,2 /a in 1.5 and a res,1 a in (see Fig. 4), the resonances 'merge' together yielding relatively large values of w/a in . What we mean by 'merging' here is that e forced (a) in-between the resonances stays larger thanẽ, and our definition of w does not disentangle the two resonances 6 . This could be understood, for instance, by looking at the curve for M d /m p = 1 in Fig. 3. These considerations explain why the contours of constant w/a in in Fig. 6 behave differently for a res /a in 1.5 compared to a res /a in 1.5.
To better understand the behavior of w/a in , in Appendix B we derive an analytic expression for the resonance widths showing that, to a good approximation, where the scaling holds for p = 1 in the limits of a p /a in → 0 and a in a res a out . First, Equation (21) shows that the width is inversely proportional to the gradient of A at a res . This explains why resonances in proximity of the disk edges are relatively narrow: in the limit of a res → a in , a out we have A → A d which diverges due to edge effects (Fig. 1), and dA/da is very large. Second, we see from Eq. (21) that the width is directly proportional to B p ∝ e p : this makes intuitive sense since e p controls the amplitude of planetesimal eccentricities (Eq. 13). It follows that more eccentric planets tend to produce wider resonances, provided thatẽ can be chosen independently from e p (though this is not clear a priori). Third, and more importantly, the scaling of Eq. (21) adequately explains the slopes of the w/a in contours: setting w/a in to a constant in Eq. (21) yields the scaling a res ∝ a 2 p , obvious in Fig. 6. Indeed, by fitting the numerical results in Fig. 6 with the functional form of Eq. (21), we find that the following expression provides an acceptable approximation of the resonance widths for our fiducial disk model (Section 2.1).
4. EXAMPLE: APPLICATION TO HD 107146 For a given debris disk exhibiting a depletion in its surface density, we can hypothesize that this depletion is due to eccentricity excitation by secular resonances mediated by the gravity of the disk and an unseen planet. We can then employ the characteristics of the secular resonances analyzed in §3 to constrain the disk-planet parameters that could configure the secular resonances appropriately and produce a depletion similar to the observations. In this section, as an exemplary case, we apply these considerations to the HD 107146 disk and identify the "allowed" parameter space subject to observational constraints. The detailed investigation of the dynamical evolution in models chosen from the allowed parameter space is carried out in the next section.

Constraints from gap location
As noted in Section 1, ALMA observations show that the HD 107146 disk, spanning from a in ∼ 30 au to a out ∼ 150 au, features a gap centered at a g ∼ 70 − 80 au (Ricci et al. 2015;Marino et al. 2018). Thus, we must choose the diskplanet parameters such that a secular resonance occurs within the depleted region. Here we opt to fix the resonance location at a res = 70 au. The analysis in Section 3.1 then allows us to uniquely determine the ratio M d /m p as a function of a p /a in , see also Eq. 19. In other words, for a given disk mass, we can deduce the planetary mass and semimajor axis that configure the resonance location appropriately (or vice versa). This is displayed by the black solid lines in Figure 7 for various values of disk mass (in M ⊕ ).
However, the disk mass can not be arbitrarily large and must be constrained. To this end, we note that observations of HD 107146 have detected around 0.25M ⊕ of dust at millimeter wavelengths (Ricci et al. 2015;Marino et al. 2018). By extrapolating this up to planetesimals of ∼ 100 km in diameter the estimated total disk mass is M d ∼ 100 − 300M ⊕ Planet semimajor axis, a [au] The grey region is ruled out as the disk would be too massive. The green region shows the excluded region where the eccentricity excitation timescales are much longer than the stellar age. The blue region is ruled out as the resulting resonance width would be much narrower than the observed gap. A planet close to the disk inner edge is ruled out (yellow region) by considerations of overlapping mean motion resonances. The red region is ruled out by direct imaging. The remaining white area represents the region where the disk-planet parameters meet all the above conditions. The lettered points represent the model parameters discussed in Sections 5.1, 5.2.1 and listed in Table 1. See the text ( §4) for details.
(assuming a size distribution with an exponent of −3.5, Ricci et al. (2015); Marino et al. (2018)). Here we choose to take 100M ⊕ as the upper limit of the disk mass. Based on this, we exclude regions in the (a p , m p ) parameter space that require more massive disks -see the gray shaded area in the upper part of Figure 7.

Constraints from stellar age and disk asymmetry
We can further constrain the parameter space by considering the age of HD 107146, which is estimated to be t age ∼ 80 − 200 Myr (Williams et al. 2004). Specifically, we require the timescale for eccentricity excitation at the resonance τ to be less than around the age of the system, i.e. τ t age . From Section 3.2, however, we know that τ depends not only on the planet's mass and semimajor axis but also on its eccentricity, see Eq. (16). To this end, we note that ALMA observations have found that the HD 107146 disk is roughly axisymmetric, with a 2σ upper limit of ∼ 0.03 for the global disk eccentricity (Marino et al. 2018). This suggests that the invoked planet must be of relatively low eccentricity. Thus, in what follows, we limit ourselves to e p ≤ 0.1.
The green curves in Figure 7 show contours along which the excitation timescale τ is 20, 200, and 2000 Myr (dashed, solid, and dotted lines, respectively) at a res = 70 au. The calculations assume e p = 0.1 -the maximum value of e p that we consider in our subsequent calculations -and use the stellar mass of HD 107146, namely M c = 1.09M (Watson et al. 2011). We first note that by definition τ ∝ 1/e p (Eq. 16): thus, for less eccentric planets the contours shown in Fig. 7 will correspond to longer timescales. Second, recall that τ is a measure of the time within which initially circular planetesimal orbits become radial, e → 1 ( §3.2). Thus, even if τ t age for a given planet (such that e(t age ) 1), we might still expect sufficient eccentricity excitation for depletion to be apparent at the resonance within the stellar lifetime. Given these considerations and the uncertainty on the age of the system, we exclude the region in (a p , m p ) parameter space corresponding to τ > 200 (0.1/e p ) Myr. This is illustrated by the green shaded region in Figure 7.

Constraints from gap width
As noted in Section 1, the gap width in the HD 107146 disk is estimated to be w obs ≈ 40 au (Marino et al. 2018). Given this, the planet's semimajor axis could, in principle, be constrained by using the analysis of resonance widths w in Section 3.3 (recall that w ∝ a p , Eq. 22). However, we recall that the resonance widths as defined in Section 3.3 do not necessarily correspond to the physical width of gaps that we expect to form. Nevertheless, we could still use the definition of w to rule out the range of planetary semimajor axes for which the resonance widths would be negligible, i.e w/w obs 1. Here we consider resonance widths to be negligible if w/w obs ≤ 0.1 (this choice is somewhat arbitrary). The blue solid line in Fig. 7 corresponds to w/w obs = 0.1; planetary semimajor axes to the left of this line are ruled out (blue shaded region).

Considerations of mean-motion resonances
Finally, we note that the planet can not be arbitrarily close to the disk. This is because the planetary orbit is surrounded by an annular 'chaotic zone' wherein particles will be quickly ejected from the system due to overlapping first-order mean motion resonances (MMR). Moreover, the secular approximation of Section 2 would break down within this zone. The half-width of the chaotic zone on either side of the planetary orbit depends on the planet's mass (Wisdom 1980;Duncan et al. 1989) such that, to lowest order 7 : 7 Strictly speaking, Eq. (23) is valid for circular orbits in the absence of collisions. The chaotic zone is known to broaden with both increasing eccentricity (Mustill & Wyatt 2012) and due to collisional effects (Nesvold & Kuchner 2015). For simplicity, we have ignored these effects.
We thereby can rule out the region in the (a p , m p ) parameter space wherein the planet's chaotic zone would lie within the disk, i.e. a p + δa p > a in . This is illustrated by the yellow shaded region near the right boundary of Fig. 7. Planetary parameters lying along the yellow solid line correspond to a p + δa p = a in ; thus, they could be responsible for setting the inner disk edge (e.g. Quillen 2006) at a in = 30 au (orange line).
We have now identified the 'allowed' range of disk-planet parameters that can produce an HD 107146-like disk structure. This is represented by the white (unshaded) region in Fig. 7, and roughly defined by a p in the range ∼ 5 − 27 au, m p between ∼ 0.1 and 25M J , and 3 M d /M ⊕ 100. Note that the allowed combinations of m p and a p are consistent with the limits placed by direct imaging of HD 107146 (Apai et al. 2008), see the dashed red curve in Fig.  7. For reference, the combinations of m p , a p and M d which we consider later in this work are labelled as models A-C in Fig. 7, see also Table 1. Note that each of these configurations correspond to τ ≈ 135 × (0.05/e p )Myr, and model A represents the fiducial configuration considered next in Section 5.1.
We remark that in the above discussion we have implicitly ignored the occurrence of an inner secular resonance at a in ; apart from the one already fixed at a res = 70 au in Fig. 7, see §3.1. This can be justified on the grounds that the inner resonance is of very narrow width except if the two resonances are close to each other, which is not the case here ( §3.3). As a result, and as we will see next, the inner resonance is irrelevant and does not have any observable effect.
Finally, we point out that equations (19), (16) and (22), combined with Eq. (23), can be applied to generate an approximate version of Figure 7 for any other observed debris disk with a gap.

EVOLUTION OF THE DISK MORPHOLOGY
In the previous section, we identified the combinations of the 'allowed' disk-planet parameters that could reproduce the observed depletion in the HD 107146 disk, see Fig. 7. We now investigate the dynamical evolution of disk-planet systems using some of these parameters. Our specific aims here are two-fold: to illustrate how secular resonances sculpt depleted regions, and to analyze more fully the disk and gap morphology in the course of secular evolution.

A Fiducial Configuration
We begin by presenting results showing the evolution of the disk surface density in the fiducial configuration, i.e. model A (see Table 1). We recall that model A is the configuration that was considered in Section 2.4, where we discussed the temporal evolution of planetesimal eccentricities and apsidal angles as a function of semimajor axis -see Fig.  2. To this end, we convert the orbital element distributions of planetesimals shown in Fig. 2 into surface density distributions. Technical details about this procedure can be found in Appendix C. The resulting maps of the (normalized) disk  Fig. 7. Column 5 presents the disk-to-planet mass ratio. Columns 6-7 present the planet's eccentricity and initial apsidal angle, whose precession period is given in column 8. (a) Model A is the fiducial configuration adopted in this work. (b) Each of the considered models have τ ≈ 135 × (0.05/ep) Myr.
surface density Σ at times corresponding to those in Fig. 2 are shown in Figure 8. For reference, in this figure we also show the planet's orbit and its pericenter position, which precesses with a period of τ sec ≡ 2π/A d,p ≈ 33 Myr (Eq. 8).
To facilitate the interpretation of our results, in Figure 9 we also show the profiles of the azimuthally-averaged disk surface density Σ as a function of radial distance r at the same times as in Fig. 8. Below we provide a detailed description of the different evolutionary stages that we identified. Stage 1 (0 ≤ t τ sec ): At early times, the disk quickly evolves away from its initial axisymmetric state by developing a trailing spiral structure (see Figs. 8a, b). This spiral structure initially starts off at the inner disk edge and propagates radially outwards with time as it wraps around the star; see also the animated version of Fig. 8. For instance, by 1 Myr at least two windings are noticeable (Fig. 8a), with the outermost prominent spiral arm occurring at ∼ 40 au. This arm moves out to ∼ 60 au by 10 Myr (Fig. 8b). A complementary view of this behavior is provided by Figs. 9(a),(b).
We note that the outermost portion of the spiral is associated with planetesimal orbits that have attained their maximum eccentricity, i.e. have completed half a precession period -see Fig. 2. Interior to this, the spirals become difficult to discern since planetesimals in this region have completed more than one precession period and their orbits are phasemixed, i.e. ∆ (a) spans the range [−π/2, π/2] -see Figs. 2(a), (b). As a result, the surface density distribution interior to the outermost spiral looks roughly axisymmetric; see e.g. panel (b) of Fig. 8. We also note that the spiral propagates outwards at a slower rate as it extends to larger radii; see panels (a)-(c) of Fig. 8 and its animated version. This follows from the fact that the planetesimal precession rate is a decreasing function of the semimajor axis (Fig. 1).
We remark that the behavior described thus far shows some parallels with the findings of Wyatt (2005), which showed that an eccentric planet launches a spiral wave which propagates throughout a massless disk. The main difference is that, in our setup, the spiral wave extends out to only about a radius of 70 au and not to the outer disk edge (as would happen in a massless disk), see Fig. 8. This is to be expected, since in our model planetesimal dynamics is dominated by the planet only within ≈ 70 au, beyond which the disk gravity becomes important -see Fig. 1 and §2.4.
Stage 2 (t ∼ τ sec ): By the time the planet has nearly completed its first precession cycle, the disk develops a clear depletion in its surface density, which effectively splits the disk into an internal and an external part (Figs. 8c, 9c). The depletion occurs around the location of the secular resonance, i.e. at a res = 70 au, where the system was designed to emplace one -see §4. The appearance of the gap is evidently correlated with the excitation of planetesimal eccentricities at and around a res , where e = t/τ ≈ 0.22 by 30 Myr (Fig. 2c).
An interesting feature of the gap is that it is of a crescent shape which points in the direction of the planet's pericenter (Fig. 8c). In other words, the gap is asymmetric in the azimuthal direction such that it is wider and deeper towards the planetary pericenter. This asymmetry is associated with the inner and outer disk components being offset relative to the star in opposite directions (Fig. 8c). Indeed, the inner part forms an eccentric structure which is apsidally aligned with the planet while the outer part is anti-aligned (see also Section 2.4) -the latter though is difficult to discern in Fig.  8 due to the smaller eccentricities in the outer parts (Fig. 2). Nevertheless, by simply looking at the azimuthally-averaged density profile we find that the gap has a radial width of ∼ 20 au (measured relative to the initial density profile, Fig. 9c). Looking at Fig. 9(c), it is also clear that this region is not depleted fully but only partially -by about a factor of two relative to the initial density distribution.
Finally, we note that the gap is surrounded by narrow overdense regions, with the one just exterior to the gap being sharper than that interior to it (see Figs. 8c, 9c). These overdensities correspond to the apocentric positions of planetesimals with semimajor axes in the depleted region. The contrast between the sharpness of the overdensities is mainly due to the apsidal angles of planetesimals at a a res being more phase-mixed than at a a res (Fig. 2c). This also justifies why these sharp overdensities are transients: they taper with time as planesimal orbits around the resonance are perturbed further (see panels d-e in Figs. 2, 8, and 9).
Stage 3 (τ sec t τ ): Further into the evolution, the structure of the gap practically remains invariant without being significantly affected by the continued growth of eccentricity around a res = 70 au (see panels d-e in Figs. 2, 8). Indeed, the gap maintains its crescent shape along with its   Table 1). The snapshots correspond to the same moments of time as in Fig. 2, and are indicated in each panel for reference. All panels have 400 × 400 pixels and share the same surface density scale (and normalization constant) as shown in the colour bar. In each panel the stellar position is marked by the yellow star, while the planet's orbit and its pericenter position are shown by the white solid line and green circle, respectively. To enhance the resolution of the images, the orbit of each planetesimal (N = 5000 in number) has been populated with 10 4 particles with the same orbital elements but with randomly distributed mean anomalies (see Appendix C). At early times (panels a, b), the planet launches a trailing spiral wave at the inner disk edge which is quickly wrapped around the star. By the time the planet has completed around one precession cycle (panel c), a crescent-shaped gap forms around the secular resonance at ares = 70 au, which is both wider and deeper in the direction of planet's pericenter. Beyond this time (panels d-e), the shape of the gap practically remains the same as it precesses while maintaining its coherence with the planet's pericenter. Note that the disk part interior to the gap is offset relative to the exterior part, where a wound spiral pattern is visible at late times (panels d-e). It is also clear that no gap forms around the secular resonance at ∼ ain. See the text (Section 5.1) for more details. This figure is available as an animation in the electronic edition of the journal. alignment with the planet as it co-precesses with the planet's apsidal line. At the same time, since the inner component of the disk precesses much faster than the outer component (Fig. 1), the degree of offset between them varies as the system evolves. This causes the gap width w g to fluctuate in time, see e.g. Figs. 9(d)-(e), with a time-averaged value of w g ≈ 18.13 ± 1.04 au. Looking at Figs. 9(d)-(e), it is also clear that the gap depth remains roughly constant such that, in a time-averaged sense, about 50 ± 3% of the initial density is depleted at the resonance.
Note that, at this stage, i.e. at t τ sec , at least one secular period has elapsed for planetesimals interior to the depletion, causing them to settle into a lopsided, precessing coherent structure (Figs. 8d-e). It is also noticeable that this structure reveals little or no evidence for surface density asymmetry between its apocenter and pericenter directions, as would have otherwise been the case if the disk were massless (i.e. pericenter or apocenter glow; see Wyatt et al. 1999;Wyatt 2005;Pan et al. 2016). This can be understood by noting that in this region, although planetesimal dynamics is dominated by the planet, the disk gravity renders the forced eccentricity  Figure 9. The azimuthally-averaged surface density of the disk Σ as a function of radial distance r from the star (solid blue lines). Each panel corresponds to each of the snapshots of the fiducial configuration (Model A, Table 1) shown in Fig. 8. The results are obtained by splitting the disk into 200 annular bins (Appendix C), and are all normalized with respect to the initial analytic surface density Σ d (a) (Eq. 1 with p = 1) at the inner disk edge, a = ain. For reference, the normalized profile of the initial Σ d (a) is shown in each panel with the solid black lines. At early times (panels a, b), the overall shape of Σ is similar to the initial profile, but with some peak features around ∼ 40 au at 1 Myr and ∼ 60 au at 10 Myr, respectively. At all times after 30 Myr (panels c-e), a clear depletion in the surface density is evident around the location of the secular resonance (ares = 70 au, dashed vertical lines). One can see that the width and the depth of the depletion are effectively constant in time (panels c-e). Note also the peak structure in the density just exterior to the depletion in panels (c)-(e). See the text (Section 5.1) for more details. This figure is available as an animation in the electronic edition of the journal.
to be more of a constant with semimajor axis rather than scaling as 1/a (see Figs. 1, 2). This hinders the occurrence of a pericenter or apocenter glow (for a more detailed discussion, see section 2.4 in Wyatt (2005)).
On the other hand, planetesimal orbits exterior to the depletion have not yet had the time to be randomly populated in phase (Fig. 2). Hence, a spiral pattern develops in this region as planetesimals undergo eccentricity oscillations. The spirals appear to wrap almost entirely around the star, and these are more noticeable closer to the depletion than to the outer disk edge (Figs. 8d-e). This can also be seen in Figs. 9(d)-(e) as a series of narrow peaks in the radial profile of Σ . This behavior can be understood by noting that planetesimals closer to the outer disk edge have smaller eccentricities (e.g. Fig. 2) and that their orbits are quickly phase-mixed as a result of their rapid orbital precession due to disk edge effects, particularly at a 130 au (e.g. Fig. 1,  §2.2). Relatedly, if we were to evolve the system for longer, planetesimals exterior to the depletion would become phase-mixed and the spiral structure would fade away. We note that, depending on the resolution of observations, the spirals in this region may or may not be visible.
Before moving on, we note that already by 1 Myr into the evolution, planetesimal eccentricities around the inner secular resonance (i.e. a res ≈ a in ) are excited to ≈ 1; see e.g. Fig. 2(a). Evidently, however, this occurs over a narrow radial range that it does not lead to the emergence of a gap (see Figs. 8,9), in agreement with our expectations from Section 3.3. This also justifies our assertion in Section 4 about ignoring the occurrence of an inner secular resonance for the purposes of Fig. 7.

Parameter variation
We now analyze the variation of the disk morphology associated with varying the disk-planet parameters relative to the fiducial values (Model A).

Variation of the planetary semimajor axis ap
We first consider the effects of varying the planetary semimajor axis a p which, we remind, all else being kept the same, is equivalent to changing the ratio M d /m p ( §3.1, §4). For ease of comparison, we choose the combinations of a p , m p and M d from Fig. 7 such that they yield the same eccentricity excitation timescale at the secular resonance τ as in model A. The parameters of the chosen models, which we label as B and C, are listed in Table 1 and are marked on Fig.  7. Note that the planet in Model C could be responsible for truncating the disk at a in = 30 au; see §4.4.
Generally, we find that the evolution of the disk morphology in each of models B and C proceeds in a similar manner as in the fiducial model (i.e. stages 1-3 in §5.1). Indeed, we observe the same qualitative behaviour: the launching of a spiral arm at a in and its outward propagation in time, the sculpting of a crescent-shaped gap around a res = 70 au by ∼ τ sec , the development of a spiral pattern exterior to the depletion at t τ sec and its subsequent potential disappearance at late times (depending on the period of secular precession at a a res ). Figure 10 summarizes the snapshots of models B and C at 100 Myr (i.e. t/τ ≈ 0.74) into their evolution. A comparison of the results shown in this figure with those of Model A (Figs. 8f, 9f) indicate that the only obvious difference is in terms of the radial width of the gaps w g . Indeed, the gap is radially narrower when the planet is closer to the star than to the inner disk edge: for a p = 7 au (i.e. Model B), on  Table 1. The results are shown after 100 Myr of evolution, corresponding to t/τ ≈ 0.74 for both models. Rows (a) and (b) show the planetesimal eccentricities and apsidal angles (relative to that of the planet) as a function of semimajor axis, respectively. The corresponding snapshots of the disk surface density and radial profiles of the azimuthally-averaged surface density are shown in the rows (c) and (d), respectively. All other notations are the same as in Figs. 2, 8 and 9. One can see that wider gaps are carved around the secular resonance at ares = 70 au when the planet is closer to the disk inner edge than to the star. It is also evident that the resultant gaps are asymmetric and of approximately the same depth in both models. See the text ( §5.2.1) for more details.
time-average, w g ≈ 11.32 ± 0.05 au, while for a p = 26.93 au (i.e. Model C) we have w g ≈ 20 ± 2 au. This dependence will be investigated in the future (Rafikov & Sefilian, in preparation), though for now we note that it is in qualitative agreement with our expectation from Section 3.3 regarding the resonance widths. Finally, we note that the gap depth is not affected by variations in planetary semimajor axis: on average, about a half of the initial density is depleted around the secular resonance regardless of a p .

Variation of the planetary eccentricity ep
The models presented thus far assumed the same planetary eccentricity of e p = 0.05. To examine its effect on the disk morphology, we considered the evolution in otherwise identical setups but differing in the value of e p by a factor of two from model A. These are referred to as models A-Loep (with e p = 0.025) and A-Hiep (with e p = 0.1) in Table 1.
Once again, we found that the evolution of the disk morphology qualitatively follows the same stages outlined in Section 5.1, but on a shorter timescale when the planet is more eccentric (recall that τ ∝ 1/e p , Eq. 16). Additionally, we identified subtle differences in the structure of the spiral arms with increasing e p . First, the spiral initially launched at a in by the planet became more open for larger e p -in agreement with the results of Wyatt (2005). Second, and relatedly, the spirals beyond the gap became more prominent with increasing e p due to the higher forced eccentricities in that region.
More importantly, however, we found that more eccentric planets give rise to wider gaps 8 -in qualitative agreement with our expectations from Section 3.3, see Eq. (21). Indeed, on time-average, we find that w g ≈ 12.8 ± 0.2 au when e p = 0.025, and w g ≈ 24.6 ± 2.8 au when e p = 0.10. This can be seen in Figure 11, where we summarize the results for models A-Loep and A-Hiep. Note that, for ease of comparison, the results are shown at different times such that t/τ (e p ) ≈ 0.74 for both models -the results must be compared with those of model A at 100 Myr (Figs. 8,9). Looking at Fig. 11, it is also evident that variations in e p do not significantly affect the fractional depth of the gap. Note also that, while planets with lower e p reduce the offset of the inner disk component, the gap retains its non-axisymmetric feature. This is largely related to the fact that for narrower gaps a smaller offset suffices for the inner component to occupy about the same fraction of the gap.

Variations with disk and planet masses
We now discuss the effects of varying the disk and planet masses while keeping other parameters unchanged. To begin with, we first recall that this requires varying both M d and m p simultaneously, i.e. while keeping M d /m p constant, to ensure that the secular resonance location where a gap is expected to form remains the same (i.e. a res = 70 au); see §3.1 and §4.1. In Figure 7, this is equivalent to moving vertically up or down relative to any of the simulation setups we have considered thus far.
As we know from Section 2, the secular precession rates scale linearly with masses (Eqs. 4 -8), whereas the forced eccentricities depend only on the ratio M d /m p (Eq. 13). Thus, varying the disk and planet masses (while M d /m p = cte) should only change the secular evolution timescale, but not the details of the secular dynamics. This simply is a restatement of the fact that scaling both M d and m p does not  Figure 11. Similar to Fig. 10, but for models A-Loep (left panels) and A-Hiep (right panels); see Table 1. Models A-Loep and A-Hiep are identical to the fiducial model A, except that they are initiated with planets with eccentricities that are lower and higher by a factor of two than in model A (i.e. ep of 0.025 and 0.10), respectively. For ease of comparison, results for each model are shown at different times (as indicated in the top panels) such that they both correspond to t/τ ≈ 0.74. One can see that increasing ep leads to a wider gap around the secular resonance at ares = 70 au, without significantly affecting the asymmetric shape of the gap and its depth. See the text ( §5.2.2) for more details.
affect the relative strength of perturbations due to the disk and the planet. Consequently, if we decrease both the disk and planet masses in any of our simulations, then the very same dynamical end-states -hence, disk morphology -will be achieved within shorter timescales, and vice versa. We note that, in principle, this scaling rule applies as long as M d , m p M c , since otherwise the Laplace-Lagrange description in Section 2 becomes unreliable (Murray & Dermott 1999). However, looking at Figure 7 we see that this limitation is not a concern in our case: the most massive allowed planet has m p ∼ 10 −2 M c .

Variations with the mass distribution in the disk
Our calculations so far have assumed a disk with density profile Σ d ∝ 1/a, i.e. with a power-law index of p = 1 in Eq. (1). We now discuss how our results would change for different values of p, when all else is kept the same. Since the slope of the surface density p effectively controls the precession rate of both the planetesimals and the planet (Eqs. 6, 8), it is natural to expect that the location of the secular resonance will shift as the mass distribution in the disk is varied; see also Eq. (18). We found that this is indeed the case, and we further confirmed that it does not qualitatively affect the evolutionary stages presented in Section 5.1.
We generally find that when a in a res a out , the resonance location shifts at most by only about 10 per cent as p is varied between 0.5 and 1.5. However, the direction in which the resonance shifts in a given setup is rather subtle to characterize for the following reasons. First, larger values of p lead to larger A d,p (and vice versa) as now more mass will be concentrated in the inner disk parts than in the outer regions, causing the planet to precess at a faster rate. Secondand relatedly -the disk induced precession rate of planetesimals A d at a a in decreases in absolute magnitude, since it is proportional to the local surface density of the disk (Eq. 6) 9 . To summarize, varying p has opposite effects on A d,p and |A d |, and it is the detailed balance between these two effects that determines whether the resonance shifts outwards or inwards in a given setup, see Eq. (17). For the parameters of HD 107146 in Figure 7, we find that the resonance shifts inwards from its nominal location, i.e. a res = 70 au, when a larger value for p is adopted (and vice versa). Thus if we were to generate a version of Figure 7 with e.g. p = 1.5 rather than p = 1, the values of M d required to reinstate the resonance at a res = 70 au would be a factor of ∼ 1.1 lower.

DISCUSSION
The results of previous sections show that the secular interaction between a low-eccentricity planet and an external, co-planar debris disk can lead to the formation of a gap in the disk. This occurs through the excitation of planetesimal eccentricities at around one of the two secular resonances arising due to the combined gravitational influence of the disk 10 and the planet. The novelty of this mechanism is that it requires the presence of only a single planet interior to a less massive disk, and is also robust, in the sense that it operates over a wide range of parameters.
As an example, we applied our model to the HD 107146 disk and investigated the general features of the disk and gap morphology in the course of secular evolution. In the following, we first discuss (in a general context) how the results of our model compare with the observed features in HD 107146 ( §6.1). We also discuss the application of our model to other systems ( §6.2). Finally, we discuss the implications of our study for determining the masses of debris disks ( §6.3), and for their dynamical modeling in general ( §6.4),

Comparison with observed structure in HD 107146
By applying our model to HD 107146, we have shown that a gap can be readily sculpted at the observed location, i.e. around 70 au (Marino et al. 2018), for a wide range of planet-disk parameters; see e.g. Fig. 7, Section 5. Additionally, our results show that the produced gaps invariably have a fractional depth of about 0.5 (Section 5), which is consistent with that observed in HD 107146 (Marino et al. 2018). While these results are encouraging, there are some issues with our model that need to be highlighted when it comes to comparing with the observational data of HD 107146 (Marino et al. 2018).
First, as already mentioned in Section 4.2, ALMA observations of HD 107146 indicate that its disk is axisymmetric and characterized by a circular gap (Marino et al. 2018). Our model, however, produces gaps that are asymmetric in the azimuthal direction (Section 5), with the disk surface density being depleted to a greater extent and over a wider region in the direction of planet's pericenter. We further found that the gap asymmetry can not be mitigated, as one might naively expect, by adopting lower values for the planetary eccentricity -see Section 5.2.2.
Second, as already stated in Section 4.3, the observed gap in HD 107146 is ∼ 40 au wide. This is larger by about a factor of two compared to the gap in our fiducial configuration (Section 5.1). In principle, our model can yield such wide gaps with a combination of high-eccentricity and large semimajor axis for the planetary orbit; see Sections 5.2.1 and 5.2.2. However, this would also impose more notable nonaxisymmetric structure on the disk which, given the discussion above, is problematic for HD 107146. Thus the conclusion is that, within the limitations of our model (for a detailed discussion, see Section 7), it is difficult to sculpt a gap as wide and as axisymmetric as that in HD 107146 without invoking additional processes. We discuss a way in which a wider gap could form as a result of disk mass depletion and secular resonance sweeping in Section 7.2.
Third, observations of HD 107146 indicate that the surface brightness of the outer and inner rings are comparable (see fig. 2 in Marino et al. 2018). Since sub-mm dust emission at a distance r scales as T (r) ∝ r −1/2 (assuming black body emission in the Rayleigh-Jeans limit), this observation suggests an increasing surface density with radius, which may seem unnatural in the context of protoplanetary disks. As a result, this has been taken as evidence for collisional depletion of planetesimals in the inner disk regions (Ricci et al. 2015;Yelverton & Kennedy 2018). Thus, if our collisionless model were applied to any physically realistic profile (i.e. with p > 0, Eq. 1), it is unlikely that we would reproduce the observed brightness peaks. However, it is possible that a shallower density slope than p = 1 could generate comparable brightness peaks at times t ∼ τ sec , when our model pro-duces an overdensity just exterior to the depletion (see Stage 2 in §5.1).
The above discussion suggests that although our mechanism acting alone can produce a structure qualitatively similar to that observed in HD 107146, it does not provide a quantitative interpretation of the observations. However, we re-emphasize that our aim in this work was not to provide a complete description of the HD 107146 disk, but rather to provide a proof-of-concept for our mechanism and its feasibility. We also stress that the limitations of our simple model need to be assessed before making any definitive conclusions (see §7 for a detailed discussion). Our results serve as a starting point to guide future, more comprehensive studies which aim to match the observations of the HD 107146 disk, or any other disk with an observed gap.
Given the potential ubiquity of gaps in debris disks (e.g. Kennedy & Wyatt 2014;Marino et al. 2020), it is also possible that future surveys will reveal a sample of disks with asymmetric gaps. Two potential candidates for such systems are HD 92945 (Marino et al. 2019) and HD 206893 (Marino et al. 2020), which we discuss next.

HD 92945
We first consider the system HD 92945 (Golimowski et al. 2011) which is often viewed as a sibling to HD 107146 in many ways. Both systems not only have stars with similar masses and ages (1M and 100 − 300 Myr, Plavchan et al. (2009)), but also their disks show some similarities in terms of their radial structure. Indeed, ALMA observations of Marino et al. (2019) show that the HD 92945 disk, extending from ∼ 50 to 140 au, is double-peaked with a gap centered at about ∼ 73 au, roughly coincident with that in HD 107146. However, and in contrast to HD 107146, the gap in HD 92945 appears to be asymmetric and is relatively narrow with an estimated width of 20 +10 −8 au (Marino et al. 2019).
These features speak in favor of our model, so we could use our results ( §3) to determine the properties of the planet and disk such that the gap is sculpted by secular resonances. Figure 12 summarizes the results of our analysis (following a similar reasoning as for HD 107146 in Section 4). We find that a companion with a semimajor axis a p in the range ∼ 3 − 50 au and mass m p between ∼ 10 −2 and 10 2 M J can produce a wide enough gap at the observed location within the stellar age, provided that 1 M d /M ⊕ 100 -see the white region in Fig. 12. These limits are in agreement with (i) direct imaging constraints (Biller et al. 2013, red curve in Fig. 12), and (ii) disk mass estimates of ∼ 100 − 200M ⊕ derived from collisional models (Marino et al. 2019).
Finally, we note that since the inner disk edge in HD 92945 is located at ∼ 50 au, i.e. further out than in HD 107146, it is possible for the planet to be on a more distant orbit than in HD 107146 (Fig. 12). However, we confirmed that this is only necessary if the true gap width is towards the upper end of its estimated range (recall that increasing a p /a in in our model leads to wider gaps). For instance, we find that invoking a planet similar to that in Model A (but with a disk of mass M d ≈ 16.4M ⊕ ) produces a ∼ 16 au wide gap, which is comparable to that observed. Future observations of this system could help to put better constraints on the disk mass and planetary properties.

HD 206893
We next consider HD 206893, a 50 − 700 Myr old F5V star, which hosts a debris disk (Marino et al. 2020) as well as one brown dwarf companion, HD 206893 B, detected using direct imaging (Milli et al. 2017). ALMA observations of Marino et al. (2020) show that this disk, extending from ∼ 30 to 180 au, features an asymmetric ∼ 27 au wide gap centered at ∼ 75 au. Given that HD 206893 B orbits interior to the disk with a p ∼ 11 au (Delorme et al. 2017), this system is ideally suited to test whether our model can reproduce the observed gap.
To assess this, we adopt the minimum possible mass of HD 206893 B (∼ 12M J , Delorme et al. 2017) and calculate, using Eq. (17), the disk mass that would place a secular resonance at the observed gap location, i.e. a res = 75 au. Assuming a surface density profile with p = 1 (Eq. 1), we find that the required disk mass is M d ≈ 170M ⊕ ; see also Eq. (19). This is roughly consistent with the disk mass estimates of Marino et al. (2020) based on collisional models. Moreover, we also confirmed that the gap width w g obtained from our model agrees well with that observed: adopting the best-fitting eccentricity of HD 206893 B, e p ∼ 0.15 (Marino et al. 2020), we find that w g ≈ 26 au after ∼ 20 Myr of evolution. If future observations with better resolution confirm that the gap in the HD 206893 disk is indeed wider towards the companion's pericenter position, this will then provide a strong support to our model.
Finally, we note that recent analyses of HD 206893 have indicated that it is likely that this system harbors a second inner companion at ∼ 2 au (Grandjean et al. 2019;Marino et al. 2020). While in this work we only considered single-planet systems, our model may easily be extended to two-planet systems (or more). In this case, depending on the strength of perturbations from the companion(s), our results both in general (e.g. Section 2) and for HD 206893 may or may not be affected significantly. However, such an analysis is beyond our scope here.

Implications for disk mass estimates
Our results may be used to infer the presence of a yetundetected planet in any system harboring a double-ringed debris disk. The inferences are, of course, degenerate with the assumed system parameters but, more importantly, they are subject to the condition that there be sufficient mass in the disk (Sections 3, 4). Thus, the detection of planets with the inferred properties will not only provide strong support to our model, but also -and more importantly -provide a unique way to indirectly measure the total mass of the debris disk M d (see e.g. Section 6.2.2). This is particularly appealing, considering the fact that M d can not be accessed using other techniques -not least without invoking theoretical collisional models to extrapolate observed dust masses to the unobservable larger planetesimals that carry most of the disk mass (see Krivov & Wyatt 2020, for a detailed discussion). This represents a promising avenue to consider in the future, in particular with the advent of new generation instruments such as JWST which could detect planets with m p 10M J at a p ∼ 10 au separations. Conversely, the results of Section 3 may be used to investigate whether or not the debris disk of a known planet-hosting system should have a gap. Future observations of such systems e.g. with ALMA looking for evidence -or lack thereof -of a gap could help in constraining the total disk mass.
6.4. The importance of disk self-gravity in dynamical modelling of debris disks The study presented here has further consequences beyond an explanation of gap formation in debris disks. Particularly, our findings strongly emphasize the need to account for the (self-)gravitational effects of disks in studies of planet-debris disk interaction. As we showed in this study, the end-state of secular interactions between a single planet and a disk having only a modest amount of mass can be radically different from the naive expectations based on a massless disk. Indeed, if it were not for the disk gravity in our model, secular resonances would have not been established and so no gap would have formed in the disk -at least not without invoking two or more planets (e.g. as done by Yelverton & Kennedy 2018), or a single but precessing planet (Pearce & Wyatt 2015).
This also highlights an important caveat related to the dynamical modelling of debris disks in general. While studies treating debris disks as a collection of massless particles seem to successfully reproduce a large variety of observed disk features by invoking unseen planets (e.g. see reviews by Krivov 2010;Wyatt 2018), their inferences about the underlying planetary system architecture may be compromised. The inclusion of disk gravity would -at least -impose modifications on the masses and orbital properties, if not numbers, of invoked planets. Thus, caution must be exercised in the interpretation of observed disk structures when the disk mass is ignored.
Recently, Dong et al. (2020) raised a similar point when it comes to ascribing observed morphologies of disks (assumed to be massless) to single planets in situations where the potential presence of a second planet is ignored. We urge a similar analysis to be performed by considering a natural hypothesis of having non-zero disk mass in contrast to the potential presence of additional planets. Although this is beyond the scope of our current work, the formalism outlined in Section 2 could provide a useful starting point for such an analysis. To summarize, the inclusion of disk self-gravity in studies of planet-disk interactions should be considered in dynamical modelling of debris disks.

LIMITATIONS AND FUTURE WORK
We now review some of our model assumptions and limitations, and discuss how relaxing them would affect our results. We plan to address these issues in future papers of this series. In this work we treated planetesimals as massless testparticles, and analyzed their secular evolution under the influence of gravity from both the planet and the debris disk. To this end, we modelled the debris disk as being passive: that is, as a rigid slab that provides fixed axisymmetric gravitational potential (see Eq. 3 and §2, disk non-axisymmetry is discussed next in §7.1.2). Thus, at first glance, it appears that instead of the planetesimals to be contributing to the collective potential of the disk, they are enslaved by the fixed disk potential given in Eq. (3). In reality, though, these two approaches are subtly similar. This is because the orbitaveraged disturbing function for a planetesimal of mass m j due to all other N massive planetesimals in a disk -in the continuum limit (i.e. N → ∞, m j ∼ N −1 ) -is equivalent to that in Eq. (3). This can be verified by somewhat tedious but straightforward calculation which requires softening the gravitational interaction between massive planetesimals, integrating radially over all planetesimals, and taking the limit of zero softening (Hahn 2003;Sefilian & Rafikov 2019).
To further justify this equivalence, we simulated the secular dynamics of disk-planet systems by modelling the disk as a swarm of N massive planetesimals, each represented as a ring 11 , that interact via softened gravity (e.g. Hahn 2003;Touma et al. 2009;Batygin 2012). We found that simulations carried out with negligible softening parameter accurately reproduce the analytical solutions presented in Section 2.3 (which is, of course, possible only when the nonaxisymmetric perturbations due to simulated disk particles are neglected, i.e. as in §2). We will present further details about this softened 'N -ring' method in an upcoming work (Paper II).

Non-axisymmetric component of disk gravity
A major limitation of this work is that we only accounted for the axisymmetric contribution of the disk gravity, ignoring its non-axisymmetric component (Section 2). That is to say, our model does not account for the non-axisymmetric perturbations that disk particles can exert both among themselves and onto the planet (see §7.1.1), even though we find that the disk naturally develops non-axisymmetry (Sections 5.1, 5.2). This omission allowed us to elucidate the key effects of disk gravity (semi-)analytically. This comes at the expense of reduced coupling within the system that inhibits the exchange of angular momentum between the disk and planet. Thus, the outlined theory serves as a first step towards a comprehensive understanding of the role played by disk gravity and its observational implications.
Previous studies of gravitating disk-planet systems (which include the full gravitational effects of disk particles) have shown that an eccentric planet could launch a long, onearmed, spiral density wave at a secular resonance in the disk (Ward & Hahn 1998;Hahn 2003Hahn , 2008. Such spiral waves propagate away from the resonance location as trailing waves with pattern speed equal to the planetary precession rate. These waves also transfer angular momentum from the disk to the planet in a way that damps the planet's eccentricity, without affecting its semimajor axis 12 (Goldreich & Tremaine 1980;Ward & Hahn 1998;Tremaine 1998;Ward & Hahn 2000).
Our idealized model is not designed to capture the full richness of such dynamical phenomena. Thus, a more sophisticated analysis is crucial, and will be the subject of future work (Paper II, in preparation). For now, we note that the non-axisymmetric component of disk gravity is not going to qualitatively affect the gap-forming picture. This is because the divergence of eccentricities at the resonance ensues from the commensurability between planetesimal and planetary precession rates, while the torques due to the planet's and disk's non-axisymmetric potentials are non-zero (e.g. Silsbee & Rafikov 2015b,a;Davydenkova & Rafikov 2018). Nevertheless, the generation of long spiral waves exterior to the depleted region may affect the disk structure and its evolution; this could be of observational relevance. Additionally, the damping of planetary eccentricity could reduce the gap asymmetry observed in our simulations via lowering e forced over time, especially in the inner disk parts. Preliminary simulations carried out with the softened 'N -ring' model confirm these expectations (Paper II).

Collisional depletion of planetesimals
We modelled the debris disk as an ensemble of collisionless planetesimals. In practice, once the disk is sufficiently stirred, planetesimals collide and break up into smaller fragments, initiating a collisional cascade (e.g. Wyatt 2008). In this process, colliding planetesimals are gradually ground to dust until they are removed from the system by radiation effects; causing the disk mass to collisionally deplete over time.
We expect collisions to preferentially deplete the disk density around the secular resonance (where e → 1 and relative velocities between planetesimals are high), in addition to the purely dynamical depletion illustrated in Sections 5.1, 5.2. This may enhance the gap depths arising from our collisionless model. Collisional evolution may also contribute to widening the gaps resulting from our model. This can be understood as follows: as the total disk mass is depleted over time, the system's precession frequencies get altered, modifying the location of the secular resonances in a timedependent way (e.g. Heppenheimer 1980;Ward 1981;Nagasawa & Ida 2000). Looking at Fig. 4, we can infer that the resonance would sweep through the disk outwards as M d decreases, potentially producing a wider gap than in our model (as then eccentricities could be excited over a larger range in semimajor axis). This could be important e.g. for the HD 107146 disk, for which our fiducial model produces gaps that are narrower than observed (see §6.1). Furthermore, we expect the shape of the resulting gap to provide information on the initial and final disk masses along with the history of mass loss. We defer detailed investigation of collisional effects to future work.

Coplanarity of the disk-planet system
Another assumption of our model is the coplanarity of the debris disk and the planetary orbit, which can be easily relaxed in future studies. Generally, however, we believe that a small but non-zero relative inclination (e.g. 5 • ) would not affect our results for eccentricity dynamics (e.g. Pearce & Wyatt 2014), as the evolution of eccentricities e and inclinations I are decoupled from each other when e, I 1 (Murray & Dermott 1999). Nevertheless, it is possible for planetesimal inclinations -similar to eccentricities -to be excited significantly at inclination resonances (e.g. Hahn 2003Hahn , 2007, where the precession rates of both planet's and planetesimal's longitudes of ascending node are commensurate. Future studies should investigate this intriguing phenomenon.

Secular approximation
We limited the expansion of the secular disturbing function to second order in eccentricities ( §2). Hence, our results are only approximate at high eccentricities, e.g. in the vicinity of the secular resonances, where it is necessary to include higher-order terms in the disturbing function (e.g. see Sefilian & Touma 2019). Such an exercise would, primarily, limit the eccentricity amplitude at the resonance (Malhotra 1998). Nevertheless, it seems unlikely that this would affect the gap formation. For instance, from Figs. 8, 9, we can see that the gap is already well-developed when eccentricities at the resonance are still rather modest, i.e. e ∼ 0.2. Higher-order terms, however, could give rise to mild quantitative differences in terms of the dynamical timescales, e.g. period of eccentricity oscillations.
We also ignored mean-motion resonances between the planet and the planetesimals. Previously, Tabeshian & Wiegert (2016) found in simulations of synthetic debris disks that gaps can be carved at the 2:1 MMR with an internal lowe planet (e p 0.1, see also Regály et al. 2018). In our simulated systems, this can occur around a in . However, as the authors explain, MMR gaps will be blurred or even washed out by high-eccentricity planetesimal orbits further out in the disk. In our case, this could be easily achieved by planetesimals in the vicinity of the secular resonance.

SUMMARY
In this work we explored the secular interaction between an eccentric planet and an external self-gravitating debris disk, using a simplified analytic model. The model is simplified in the sense that it only accounts for the axisymmetric component of the disk (self)-gravity, ignoring its non-axisymmetric contribution. Despite this limitation, however, this is the first time (to our knowledge) that the effects of disk gravity have been considered analytically in such detail in the context of debris disks. We used the analytic model to assess the possibility of forming gaps in debris disks through excitation of planetesimal eccentricities by the secular apsidal resonances of the system. We summarize our key results below.
(i) When the debris disk is less massive than the planet, 10 −4 M d /m p 1, the combined gravity of the disk and the planet can mediate the establishment of two secular apsidal resonances in the disk.
(ii) We map out the behavior of the characteristics of the secular resonances -i.e. locations, time-scales, and widths -as a function of the disk and planet parameters. In particular, we find that one of the secular resonances can lead to the formation of an observable gap over a broad region of parameter space.
(iii) As an example we applied our results to HD 107146 and HD 92945, and showed how the properties of a yet-undetected planet, together with the mass of the debris disk, can be constrained to produce a gap at the observed location. In the case of HD 206893, we find that the directly imaged companion can sculpt the observed gap if the debris disk is ≈ 170M ⊕ in mass.
(iv) By investigating the secular evolution in such systems, we identified three distinct evolutionary stages which occur on timescales measured relative to the planetary precession period. We find that the gap forms by the time the planet has completed around one precesional cycle, on a timescale of tens of Myr.
(v) Independent of the system parameters, the gap carved around the secular resonance is asymmetric: it is both wider and deeper in the direction of the planetary pericenter. Additionally, its fractional depth is always about 0.5. The gap width, however, increases with increasing planetary semimajor axis and/or eccentricity.
(vi) More generally, our results suggest that the gravitational potential of debris disks can have a notable effect on the secular evolution of debris particles. We advocate the inclusion of disk gravity in studies of planetdebris disk interactions.
The mechanism presented here represents what is arguably the simplest pathway to forming gaps in debris disks, akin to those observed in HD 107146, HD 92945 and HD 206893. It may indeed obviate the need for invoking more complicated scenarios, e.g. multiple planets interior to or within the disk.
Finally, we remark that the present work should be envisaged as a first step towards an in-depth exploration of the effects of disk gravity in planet-debris disk interactions. In a forthcoming paper (Paper II), we will extend our current calculations using numerical techniques to properly account for the full gravitational effects of the disk. In the future, we also plan to investigate the role of disk gravity in shaping debris disk morphologies other than gaps. (2) 3/2 1 u a p a in du. (A8) Figure 13 shows the behavior of φ c 1 and φ c 2 as a function of a p /a in , computed for different values of p, q and δ. For clarity, we have plotted the curves of φ c 1 and φ c 2 in separate panels. We see that φ c i (i = 1, 2) mainly depend on a p /a in , showing weak dependence on the disk model. Indeed, regardless of (p, q, δ), we have φ c i → 1 for a p /a in → 0, while in the limit a p /a in → 1 we see that φ c i diverge. This divergence follows from the fact that b (m) 3/2 (α) → (1 − α) −2 when α → 1.

B. ANALYTIC EXPRESSION FOR RESONANCE WIDTHS
The width w of a given resonance at a = a res can be approximated by using the fact that Additionally, Equation (20)

allows us to write
A (a res ± w/2) ≈ A d,p ∓ẽ −1 B p (a res ) × sgn [dA/da] ares , (B10) where sgn(x) = x/|x| is the sign function introduced to account for the fact that resonances occurring at a in have dA/da > 0, while those further away have dA/da < 0; see Fig. 1. Substituting Eq. (B10) into Eq. (B9), we thus arrive at w a in ≈ 2 a in B p (a)ẽ −1 dA/da ares .
The above expression can be further simplified by considering the approximate forms of A p and A d in the limits of a p /a res → 0 and a in a res a out , respectively. In this case, we can approximate the derivative of A = A p + A d in the following fashion and expression (B11) reduces to w a in ≈ 4 a res a in B p (a res )ẽ −1 7A p (a res ) + (2p − 1)A d (a res ) .
Inserting the condition for secular resonances, i.e. Eq. (17) or Eq. (19), into the above expression for p = 1, and taking the limits a p /a in → 0 (so we can use the asymptotic behavior of b (m) s (α)) and a in a res a out , we arrive at the scaling relationship given by Eq. (21).

C. CONSTRUCTING MAPS OF DISK SURFACE DENSITY
Here, we provide some technical details about how we convert the eccentricity-apsidal angle distribution of planetesimals into maps of disk surface density.
We first begin by assigning a mass m i to each considered planetesimal in a given annulus of the disk (which in this work are N = 5000 in number, §2.3). Given that in our calculations the planetesimals are initiated on circular orbits, the planetesimal masses can be determined from their initial semimajor axis distribution -which remains constant in the secular approximation. This can be done by using the relationship dm(a) = 2πaΣ d (a)da (Statler 2001;Davydenkova & Rafikov 2018) relating the mass distribution per unit semimajor axis to the density distribution (which in our case is given by Eq. 1 with p = 1, §2.1). The self-consistency of this initial mass assignment to planetesimals -which are essentially treated as massless particles in our analytical model (see §2) -is discussed in Section 7.1.1.
At a given time of the evolution, we then populate every planetesimal's orbit with N np = 10 4 new particles: each with mass m i /N np , orbital elements similar to the parent planetesimal, but with randomly distributed mean anomalies l between 0 and 2π. This procedure is motivated by the orbit-averaging principle (Murray & Dermott 1999). We also note that this procedure effectively increases the number of evolved planetesimals (from N to N × N np ), enhancing the quality of the resultant maps of disk surface density. Next, we numerically solve for each new particle's eccentric anomaly using Kepler's equation (Murray & Dermott 1999), and compute the position of each particle along its orbit via (Sridhar & Touma 1999;Binney & Tremaine 2008): Finally, we bin the positions of all N × N sp particles in the Cartesian system centered at the host star (with a resolution of 400 × 400 pixels in this work), compute the total mass per bin and divide by its area to arrive at the disk surface density distribution, Σ, at a given time. Note that this also allows us to trivially obtain the azimuthally-averaged surface density profile Σ as a function of radial distance r, where r = √ X 2 + Y 2 = a(1 − e cos ), by splitting the disk into annular bins.