Collision rates of planetesimals near mean-motion resonances
MMNRAS , 1–16 (2020) Preprint 26 February 2021 Compiled using MNRAS L A TEX style file v3.0
Collision rates of planetesimals near mean-motion resonances
Spencer C. Wallace ★ , Thomas. R. Quinn , Aaron C. Boley Astronomy Department, University of Washington, Seattle, WA 98195 Department of Physics and Astronomy, University of British Columbia, Vancouver BC, Canada
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
In circumstellar discs, collisional grinding of planetesimals produces second-generation dust.While it remains unclear whether this ever becomes a major component of the total dustcontent, the presence of such dust, and potentially the substructure within, it can be used toexplore a disc’s physical conditions. A perturbing planet produces nonaxisymmetric structuresand gaps in the dust, regardless of its origin. The dynamics of planetesimals, however, willbe very different than that of small dust grains due to weaker gas interactions. Therefore,planetesimal collisions could create dusty disc structures that would not exist otherwise. Inthis work, we use N-body simulations to investigate the collision rate profile of planetesimalsnear mean-motion resonances. We find that a distinct bump or dip feature is produced in thecollision profile, the presence of which depends on the libration width of the resonance andthe separation between the peri- and apocenter distances of the edges of the resonance. Thepresence of one of these two features depends on the mass and eccentricity of the planet.Assuming that the radial dust emission traces the planetesimal collision profile, the presenceof a bump or dip feature in the dust emission at the 2:1 mean-motion resonance can constrainthe orbital properties of the perturbing planet. This assumption is valid, so long as radial driftdoes not play a significant role during the collisional cascade process. Under this assumption,these features in the dust emission should be marginally observable in nearby protoplanetarydisks with ALMA.
Key words: planet-disc interactions – planets and satellites: formation – submillimetre:planetary systems – methods: numerical
Recent observations of circumstellar discs by ALMA have revealeda rich variety of substructure. Features such as gaps and asymmetries(ALMA Partnership et al. 2015; Pérez et al. 2016; Isella et al. 2016;Andrews et al. 2016; Cieza et al. 2016) in the emission providediagnostics for the physical processes that drive the evolution of thediscs. In many gas-rich systems, these gap features are argued toindicate the presence of a giant planet, either embedded in the disc(Dipierro et al. 2015) or orbiting adjacent to it, although the detailsof gas-dust interactions can have strong effects on the resultingmorphologies (Dong et al. 2018). In gas-poor/debris systems, agiant planet perturber can also influence the structure of the dustcontinuum emission. For example, a misaligned giant planet canproduce nonaxisymmetric features such as warps (Augereau et al.2001), highly eccentric perturbers can produce structures throughsecular interactions (Pearce & Wyatt 2014; Pearce et al. 2015),and mean motion resonances (MMRs) can open gaps (Nesvold &Kuchner 2015; Tabeshian & Wiegert 2016, 2018).Collisions between planetesimals are thought to be the princi-pal source of dust for debris discs (see Wyatt 2008). Although some ★ E-mail: [email protected] (SCW) amount of primordial dust is likely still present during the earlystages of debris disc evolution, ongoing collisions between smallbodies will augment this and could be used to trace the underlyingdynamical activity of the planetesimals. This process was directlydetected when the putative object Fomalhaut b was confirmed to bean expanding debris cloud, most likely from a planetesimal collision(Gaspar & Rieke 2020).Dobinson et al. (2013) showed that collisional dust that is gen-erated near the gap opened by a giant planet in the inner regions ofa transition disc should produce a distinct observable marker thatcould be used to infer the presence of the planet. In a followupstudy, it was found that morphological differences in the dust emis-sion could be used to determine whether the planet had a circular oreccentric orbit (Dobinson et al. 2016). Moreover, substructure dueto mean-motion resonances (MMRs) with a giant planet may holdadditional clues that could be used to constrain the properties of theplanet. The width of a resonance is set by both the mass of the per-turbing planet and the unperturbed eccentricity of the planetesimals,which is in turn set by the eccentricity of the planet through secularforcing. Therefore, the planetesimal collision profile and potentiallythe structure of second generation dust produced near MMRs mayencode information about both the mass and eccentricity of theplanet. © a r X i v : . [ a s t r o - ph . E P ] F e b S. C. Wallace et al.
The dynamics governing the motion of bodies near MMRsis extremely nonlinear, as is determining what the collision ratesbetween planetesimals should look like in these regions. For a col-lection of bodies massive enough to experience the effects of grav-itational focusing, a large eccentricity dispersion tends to reducethe probability of collision, while enhancements in surface densitytends to increase it. Due to conservation of the Jacobi energy, MMRssimultaneously enhance the local eccentricity dispersion and alsoenhance the surface density adjacent to the resonance (Richardsonet al. 2000; Boley 2017). Unfortunately, collision detection in anN-body simulation is extremely computationally expensive. So far,studies of planetesimal dynamics near MMRs in which the plan-etesimals are directly resolved have involved either collisionless testparticles (Boley 2017; Tabeshian & Wiegert 2016, 2018) or integra-tion times that do not fully capture the dynamics of the resonances(Richardson et al. 2000; Dobinson et al. 2013).To further elucidate this subject, we use the tree based N-bodycode ChaNGa (Jetley et al. 2008; Menon et al. 2015) to follow thecollisional evolution of a planetesimal disc under the gravitationalinfluence of a Jupiter-sized body. Because particle positions aresorted into a tree structure, neighbor finding and collision detectioncan be done quickly and efficiently. This considerably relaxes theconstraints on resolution and integration time. With this toolset,we explore the collision rate structure of a planetesimal disc in thevicinity of mean-motion resonances with a planet. In particular, wewould like to determine (1) what dynamics govern the collision rateprofile near MMRs, (2) whether MMRs leave a detectable signaturein the collisionally-generated dust, and (3) whether these signaturescan be used to determine or constrain the orbital properties of theperturbing planet.This work is organized in the following way: In section 2, weprovide an overview of the relevant dynamics that drive the evolutionof a planetesimal disc under the gravitational influence of an exter-nal perturber. In section 3 we provide an overview of the N-bodycode used and describe the initial conditions chosen for five simula-tions in which a perturbing giant planet is given various masses andeccentricities. Section 4 presents the results of these simulations andwe take an in-depth look at the collision rate profiles of the plan-etesimals near the MMRs. In section 5, we extrapolate the collisionrate down to smaller bodies and motivate a direct correspondencebetween the planetesimal collision profile and the resulting dustprofile. In section 6, we use the resolved collisions to generate syn-thetic dust emission profiles that would be detected with observingfacilities like ALMA. Under our simplifying assumptions, we showthat a characteristic bump or dip feature appears in the dust emissionnear the interior 2:1 MMR, the presence of which depends on themass and eccentricity of the perturbing planet. If the 2:1 MMR canbe identified (presumably, by identifying another prominent MMR,such as the 3:1, and measuring the spacing between the two), wediscuss the potential for using this feature to place constraints onthe mass and eccentricity of the planet. In addition, we highlight thecaveats of connecting the dust and planetesimal collision profiles insuch a simple way. Finally, we conclude in section 7.
We begin by providing a description of the dynamical effects re-sponsible for shaping the orbital distribution of the planetesimals.The purpose of this is twofold: (1) to motivate the initial conditionsused for the simulations described in section 3.2 and (2) to justifythe exclusion of certain physical effects from our simulations. Here, we focus solely on the physics relevant for full-sized ( ∼
100 km)planeteismals and save a discussion of the effects on smaller bodiesgenerated through collisions for section 5.
The most direct and widespread effect that a perturbing giant planetwill have on a planetesimal disc is through secular forcing of theplanetesimals. This will cause the complex eccentricities of theplanetesimals to take on a time-independent forced value, given by(Wyatt et al. 1999) as 𝑧 𝑓 = 𝑏 / ( 𝛼 ) 𝑏 / ( 𝛼 ) 𝑒 𝑔 exp 𝑖𝜛 𝑔 . (1)Here, 𝛼 = 𝑎 𝑔 / 𝑎 where 𝑎 𝑔 and 𝑎 are the semi-major axes of thegiant planet and a planetesimal, respectively. 𝑒 𝑔 and 𝜛 𝑔 are theeccentricity and longitude of pericenter of the giant planet, and 𝑏 𝑗𝑠 ( 𝛼 ) is a Laplace coefficient given by (Murray & Dermott 1999)(ch. 6, pg. 237, eq. 6.67) as 𝑏 𝑗𝑠 ( 𝛼 ) = 𝜋 ∫ 𝜋 cos 𝑗𝜓 𝑑𝜓 (cid:0) − 𝛼 cos 𝜓 + 𝛼 (cid:1) 𝑠 . (2)Without any nearby secular or mean motion-resonances, equa-tion 1 will completely describe the eccentricities and longitude ofpericenter orientations of the planetesimals. Additional forces dueto two-body scattering between planetesimals, along with aerody-namic gas drag will add an additional free component to the com-plex eccentricity, which will be randomly oriented. The magnitudeof the free eccentricity describes how dynamically hot the planetes-imal disc is and sets the random encounter speeds of planetesimals.When the dynamical excitation of the disc is driven by gravitationalstirring, the magnitude of the free eccentricity can be described bya Rayleigh distribution (Ida & Makino 1992).For the case of a Jupiter mass planet at 5.2 au perturbing a testparticle at 3 au, the timescale for secular forcing is approximately12,000 years. The details of this calculation can be found in appendixA. In regions where there are commensurabilities between frequencies,Laplace-Langrange secular theory breaks down and bodies are sub-ject to strong perturbations. For the purposes of this study, we willignore secular resonances, which generally occur on rather largetimescales and will focus on mean-motion resonances. A MMR oc-curs when the orbital period ratio between two bodies is sufficientlyclose to 𝑃𝑃 (cid:48) = 𝑝 + 𝑞𝑝 , (3)where 𝑝 and 𝑞 are integers > 𝑃 / 𝑃 (cid:48) corresponds to a 𝑝 + 𝑞 : 𝑝 resonance. If the perturber is much more massive than theother body and all of the bodies lie in a near-Keplerian potential,the condition for MMR is set by 𝑎𝑎 (cid:48) = (cid:18) 𝑝𝑝 + 𝑞 (cid:19) / . (4) MNRAS , 1–16 (2020) ollision rates of planetesimals If we further assume that the two bodies are orbiting in thesame plane, the motion of the bodies near resonance is determinedby the behavior of a critical angle 𝜙 = ( 𝑝 + 𝑞 ) 𝜆 (cid:48) − 𝑝𝜆 − 𝑞𝜛, (5)where 𝜆 = 𝜛 + 𝑀 is the mean longitude of a body, with 𝑀 be-ing the mean anomaly. For bodies in resonance, the critical anglewill librate around an equlibrium value, while this angle will cir-culate outside of resonance. For small eccentricities, this behavioris analogous to the motion of a pendulum. Furthermore, variationsin the critical angle are coupled to changes in the mean motion andsemimajor axis (Murray & Dermott 1999). An important point tonote, which we will revisit later, is that the variation frequency ofthis angle approaches zero near the edge of a resonance. The widthof a resonance can be defined by determining the largest variationin semimajor axis that permits librational, rather than circulatorymotion of 𝜙 . Calculations for the widths of first and second orderinterior MMRs are shown in appendix B, along with a timescalefor the libration. For the most prominent interior mean-motion res-onances, including the 2:1 and 3:1, the timescale associated withthese oscillations driven by a Jupiter mass planet is ∼ A simple analytic model for the collision rate of a planetesimalpopulation is given by Safronov (1969) as 𝑛𝜎𝑣 = 𝑛𝜋𝑠 (cid:16) + 𝐺𝑚 / 𝑠𝑣 (cid:17) 𝑣, (6)where 𝑛 is the number density of the population, 𝑠 and 𝑚 are theradii and masses of the bodies and 𝑣 is their typical encountervelocity. The encounter velocity is often described in terms of therms eccentricity (cid:10) 𝑒 (cid:11) / and inclination (cid:10) 𝑖 (cid:11) / of the population(e.g. Lissauer (1993)) as 𝑣 = √︃(cid:10) 𝑒 (cid:11) + (cid:10) 𝑖 (cid:11) 𝑣 𝑘 , (7)where 𝑣 𝑘 is the local Keplerian velocity. The second term in equa-tion 6 can be thought of as an additional enhancement to the col-lision cross section due to gravitational focusing. When the typicalencounter velocity is small compared to the mutual escape velocityof the planetesimals, the collision cross section greatly exceeds thegeometric value. An important feature of equation 6 is that for a fixedvalue of 𝑛 , the collision rate exhibits a global minimum as a functionof 𝑣 . For small 𝑣 , gravitational focusing facilitates more collisions,while for large 𝑣 the encounter rate simply becomes so great that thecollision rate again increases, even though gravitational focusing ismostly suppressed. This point will become relevant in section 4.1when we examine the qualitative changes in the collision rate nearmean-motion resonances.Although equation 6 works well to describe the collision ratefor a homogeneous collection of planetesimals, some problems arisewhen regions of commensurabilities are introduced. Namely, theresonances cause the orbits of planetesimals to precess, and an in-terface between secularly aligned and randomly oriented orbits ariseon each side of the resonance. Additionally, an interface betweendynamically cold and hot planetesimals develop. These two effectscause the number density, collision cross section and encounter velocity to rather abruptly vary across the boundaries of the reso-nance. For this reason, an N-body treatment in which collisions aredirectly resolved is necessary to understand how the collision ratevaries near the MMRs. Over the course of many orbits, the residual gas from the primordialnebula can damp the eccentricities and inclinations of bodies. Fora planetesimal-sized body, gas drag operates in the Stokes regimeand the timescale for aerodynamic forces to significantly alter itsrelative velocity is given by (Adachi et al. 1976) 𝑡 𝑠 = 𝑚𝐶 𝐷 𝜋𝑠 𝜌 𝑔 𝑣 𝑔 , (8)where 𝑚 and 𝑠 are the mass and radius of the planetesimal. 𝐶 𝐷 is adrag coefficient which is of order unity, 𝜌 𝑔 is the local density of thegas and 𝑣 𝑔 is the headwind velocity of the gas experienced by theplanetesimal. At 3 au, the gaseous component of the solar nebula hasa density of 3 × − g cm − and the typical headwind experiencedby a body on a Keplerian orbit is ∼ ,
000 cm s − (Hayashi 1981).For a 100 km body with a density of 2 g cm − (assuming a mixtureof ice and rock) the stopping timescale is about 1 Myr. This ismuch longer than the timescales associated with secular forcing andlibration due to mean-motion resonances, as discussed above, butshorter than the typical protoplanetary disk lifetime and potentiallyshorter the planetesimal formation timescale. For this reason, wedo not model the effects of gas drag in the simulations, althoughwe consider its effects when constructing initial conditions, whichis discussed in the next section. To follow the dynamical and collisional evolution of a planetesimaldisc, we use the highly parallel N-body code ChaNGa . This code,which is written in the CHARM++ parallel programming language,was originally designed for cosmology simulations and has beenshown to perform well on up to half a million processors (Menonet al. 2015). Using similar methods to PKDGRAV, which has beenused for numerous studies of planet formation (Richardson et al.2000; Leinhardt & Richardson 2005; Dobinson et al. 2013; Lein-hardt et al. 2015), ChaNGa calculates gravitational forces using amodified Barnes-Hut (Barnes & Hut 1986) tree with hexadecapoleexpansions of the moments and integrates the equations of motionusing a kick-drift-kick leapfrog scheme. All of the simulations weperform use a node opening criterion of 𝜃 𝐵𝐻 = 0.7. More informa-tion about the code can be found in Jetley et al. (2008).We have recently modified ChaNGa to handle solid-body col-lisions between particles, by assigning them a fixed radius, ratherthan treating them as tracers of a fluid with a characteristic soften-ing length. We provide a brief summary of the collision model here,although a full description is provided in Wallace & Quinn (2019).This work is largely based on the solid-body collision implemen-tation in PKDGRAV, which is detailed in Richardson (1994) andRichardson et al. (2000). Imminent collisions are found during the A public version of ChaNGa can be downloaded from
MNRAS000
MNRAS000 , 1–16 (2020)
S. C. Wallace et al.
Figure 1.
The positions of the remaining planetesimals at the end of the e1m2 (left), e2m2 (center) and e3m2 (right) simulations. The red dot indicates theposition of the giant planet, and the dashed line points in the direction of the planet’s longitude of perihelion. Non-axisymmetric gaps are apparent near thelocations of MMRs. At higher eccentricities, more resonances become visible. At the 2:1 MMR, gap features at 𝜃 = 0 and 𝜃 = 𝜋 follow the giant planet in itsorbit. Table 1.
Summary of Simulations RunName Mass of Planet Eccentricity of Planete1m2 1.0 𝑀 𝑗𝑢𝑝 𝑒 𝑗𝑢𝑝 e2m2 1.0 𝑀 𝑗𝑢𝑝 𝑒 𝑗𝑢𝑝 e3m2 1.0 𝑀 𝑗𝑢𝑝 𝑒 𝑗𝑢𝑝 e2m1 0.5 𝑀 𝑗𝑢𝑝 𝑒 𝑗𝑢𝑝 e2m3 2.0 𝑀 𝑗𝑢𝑝 𝑒 𝑗𝑢𝑝 drift phase of each time step by extrapolating the positions of bodiesforward using the velocity calculated during the first kick. For eachbody, the nearest 64 neighbors are tested for an imminent collision.If an imminent collision is detected, the two particles of mass 𝑚 and 𝑚 are merged together together to form one single larger bodywith the same density and a mass of 𝑚 + 𝑚 . The resulting bodyis then imparted with the center of mass position and velocity ofthe two colliders. Because the resolution of a collision can result inanother imminent collision, these events must be handled one at atime, with the soonest collision being resolved first. For this reason,another collision search is run each time a collision is resolved,which continues until there are no more imminent collisions dur-ing the current time step. During the course of the simulation, nodebris is created by collisions. To model the collisionally generateddust distribution, we use the statistics of the resolved planetesimalcollisions to build a dust profile. This process is discussed in detailin section 5. In total, five simulations are run, which are listed in table 1. Thesecond is a “nominal” case, in which the perturbing planet’s massand eccentricity are set to that of Jupiter’s. In the other four cases, themass or eccentricity is altered by a factor of two from the nominalvalue. In all cases, the perturbing giant is placed on a 5.2 au orbitaround a 1 𝑀 (cid:12) star. The planetesimal disc extends from 2.2 to 3.8au, which covers the two most prominent mean-motion resonanceswith the giant planet, the 2:1 at 3.27 au and the 3:1 at 2.5 au. The planetesimal disc loosely follows a minimum-mass solar nebulasurface density profile (Hayashi 1981) Σ = Σ 𝑟 − 𝛼 , (9)with 𝛼 = 3/2 and Σ = 10 g cm − . Planetesimals are given a bulkdensity of 2 g cm − , which corresponds to a mixture of ice androck, and are given a diameter of 300 km. This choice of parametersis meant to mimic the setup used in simulation B from Richardsonet al. (2000). One difference from the aformentioned study is thatwe use a narrower annulus, which reduces the number of particlesrequired from 10 down to roughly 500,000. This allows us to usea finer base timestep size, which, as we will discuss in section 4,appears to be the reason why no features were seen in the collisionprofile near the resonances by Richardson et al. (2000).Another difference from Richardson et al. (2000) is that thedynamical effects of secular forcing by the giant planet, along withthe effects of viscous stirring and gas drag on the planetesimals,are built into the initial conditions. Although we do not model theeffects of gas drag on the planetesimals during the simulation, weconstruct the initial conditions such that the effects of gas drag andviscous stirring are in balance. The viscous stirring timescale of theplanetesimal disc is much longer than our chosen integration time,and so the dynamical excitation of the disc (excluding resonances)stays constant, even without the inclusion of damping forces fromthe gas. This is done by first calculating the equilibrium eccentricity 𝑒 𝑒𝑞 due to viscous stirring and gas drag as a function of semima-jor axis according to equation 12 of Kokubo & Ida (2002). Theeccentricities of the bodies are drawn randomly from a Rayleighdistribution with a mode of 𝑒 𝑒𝑞 , while the inclinations are drawnfrom a similar distribution with a mode of 𝑒 𝑒𝑞 /2 (Ida et al. 1993).The arguments of perihelion 𝜔 , longitude of ascending nodes Ω andthe mean anomalies 𝑀 of the bodies are drawn uniformly ∈ [ , 𝜋 ) .To account for the effects of secular forcing by the planet, theeccentricity vectors of the planetesimals are first decomposed intoreal and imaginary components: 𝑧 = ( 𝑘, 𝑖ℎ ) = 𝑒 𝑒𝑥 𝑝 ( 𝑖𝜛 ) (10) MNRAS , 1–16 (2020) ollision rates of planetesimals and a forced component is added to ℎ according to equation 1 (wherewe have set 𝜛 𝑔 = 0). For the purposes of the integrator, there are two relevant timescalesin this system. The first is the orbital dynamical time √︁ 𝑎 / 𝐺 𝑀 (cid:12) .Through experimentation, we have found that using a base timestepsize with ChaNGa of 3% of an orbital dynamical time at the inneredge of the disc keeps the integration symplectic. Although doing sopreserves orbital frequencies, the errors associated with precessionfrequencies do not average out to zero. This is especially importantgiven the effects of secular forcing by the planet. To mitigate this,we further reduce the base timestep size by a factor of 4 and use Δ 𝑡 = 0.0025 yr. Doing so prevents the longitude of perihelia ofplanetesimals at the inner edge of the disc, where this artificialprecession is most severe, from drifting by more than the intrinsicspread in 𝜛 due to the free eccentricity over the course of ourintegrations.An additional timescale is set by the dynamical time ofthe planetesimals ( ∼ / √︁ 𝐺 𝜌 ), which is about 45 minutes. Thistimescale must be resolved in order to properly follow close gravi-tational encounters between these bodies. To resolve the base timestep and the dynamical timescale of planetesimals simultaneously,we use a two-tiered time stepping scheme, following(Leinhardt et al.2015). To start, all bodies are placed on the base time step. A firstpass of collision detection is then run in which the radii of all bodiesare inflated by a factor of 2.5. Any bodies with imminent collisionspredicted using the inflated radii are placed on a time step that is afactor of 16 smaller than the orbital time step. Although this is still afactor of roughly 7 larger than the dynamical time of a planetesimal,we found no difference in the collision rate when using any smallerof a minimum time step size.The purpose of the two-tiered schemeis to properly resolve the gravitational interactions between any bod-ies that undergo a close encounter. This prevents the coarser basetime step from reducing the effectiveness of gravitational focusing,while minimizing the additional computational expense.
All five simulations are evolved for 5,000 years, which is roughly400 complete orbits of the perturbing planet. Because the effects ofthe mean motion resonances are not built into the initial conditions,the simulations must be run long enough for the distribution oforbital elements to reach equilibrium near the resonances before thecollision rate is measured. To do so, we allow each simulation to runfor 2,000 years before we begin recording any collision statistics.This is comparable to the libration period, which was calculated insection 2.2 for the 3:1 and 2:1 resonance. Although one might worrythat a single libration period may not be long enough for the orbitalstructure due to the resonances to develop, we find that the shapeof the semimajor axis-eccentricity distribution of the planetesimalsreaches a steady state by this time. This is largely due to the factthat there are enough bodies in or near the resonances to allow anensemble average of phase space, rather than a time average. Asan additional confirmation, we find that the time evolution of thecollision rate near the resonances maintains a steady state in all fivesimulations after ∼ Figure 2.
The semimajor axes and eccentricities of the remaining planetesi-mals are shown, with the locations of prominent resonances indicated by thevertical dashed lines. Libration of the critical angle drives large variations ineccentricity, which produce spikes in the a-e plane. Between the resonances,the nonzero eccentricity is due to secular forcing by the planet.
We begin by examining simulations e1m2, e2m2 and e3m2. Thepositions of the planetesimals in the x-y plane after 5,000 years ofintegration are shown in figure 1. In all cases, the coordinate systemis rotated so that the giant planet lies at 𝜃 =
0. The longitude ofperihelion of the planet is shown by the dashed line. Resonanceswith the perturbing planet are visible as nonaxisymmetric gaps inthis figure. Upon close inspection of a series of simulation snap-shots, these gaps appear to follow the planet in its orbit, ratherthan aligning themselves with the longitude of perihelion. A similarsubstructure, which is most noticeable near the 2:1 MMR, revealsitself in Richardson et al. (2000) (see their figures 3c and 3f) andTabeshian & Wiegert (2016) (see their figure 3a). This structure isalso present in the Boley (2017) simulations, although it was notreported at the time (see Appendix C). It is worth noting that bothRichardson et al. (2000) and Boley (2017) started with a completelycold planetesimal disc. Thus, the presence of these features seemsrobust to the choice of initial conditions.The effects of the resonances become much more apparent insemimajor axis-eccentricity space, which is shown in figure 2. In allcases, the 3:1, 2:1 and 5:3 resonances are readily visible as “spikes”
MNRAS000
MNRAS000 , 1–16 (2020)
S. C. Wallace et al.
Figure 3.
Shown here are the longitudes of pericenter and semimajor axes of the remaining planetesimals. The dashed lines again indicate the locations ofprominent MMRs. Close to the resonances, the critical angle librates and drives fast precession of the longitudes of pericenter of the bodies. This overpowersthe secular forcing by the planet and effectively randomizes the orientations of their orbits. in the eccentricity that bend slightly inward (due to the conservationof the Jacobi energy). In the e3m2 simulation, features also appearnear the 5:2, 7:3 and 5:3 resonances. The absence of these finerfeatures from the e1m2 simulation can be explained by the fact thatthe strength of a resonance scales with 𝑒 𝑞 (Malhotra 1994).Another important effect of the resonances is visible in fig-ure 3, which shows the orientation of the longitude of periheliaof planetesimals in the disc. Inside of the resonances, orbits ofplanetesimals quickly precess, and their orientations are effectivelyrandomized. An important point to note is that this strong precessioneffect quickly disappears beyond the boundaries of the resonance.This turns out to be key to explaining the nonaxisymmetric structureseen in figure 1, which will be addressed in more detail below.Next, we examine the statistics of collisions resolved in each ofthe simulations. The 3D positions and velocities of the two collidingbodies are recorded to a table at the moment of impact. We derive theKeplerian orbital elements of a collision using these positions andvelocities. First, we examine the semimajor axis of the first collider,which is shown in figure 4. During a collision, the “first” particleis defined as the more massive of the two. However, nearly all ofthe collisions happen between the initial, equal-mass planetesimals.In this case, the distinction is set by the collision search algorithmand is rather arbitrary. In all of the plots where we show collisionstatistics, we have verified that using the “first” or “second” colliderdoes not qualitatively change any of the features. We find that someof the features present in figure 4 in this and subsequent figures arehighly sensitive to the number of bins and the location of the binedges. For this reason, we construct a probability density function(PDF) of the collisions using a Kernel Density Estimate (KDE). Weuse the neighbors.KernelDensity function from the sklearn(Pedregosa et al. 2011) package to construct our KDEs, using agaussian kernel with a FWHM of 0.02 au. The curves shown infigure 4 are normalized such that the area underneath is equal to 1.Near the stronger resonances, there are noticeable suppressionsor enhancements to the local collision rate. This contrasts with thefindings of Richardson et al. (2000), who simulated a similar setupand found no discernible features near the MMRs. We attributethe differences to a more conservative timestepping criterion inour simulations, along with the inclusion of secular perturbations in our initial conditions. The most prominent features appear asa local maximum in the collision rate near the 2:1 MMR and alocal minimum near the 3:1 MMR. At higher forced eccentricities,features near the 5:2, 8:3 and 5:3 resonances are also visible, dueto the steep sensitivity of higher order resonant perturbations toeccentricity (Malhotra 1994).Although the features in figure 4 near the 3:1 and 2:1 MMRsappear qualitatively different, they can be explained by one singledynamical process. As discussed in section 2.3, the collision ratebetween planetesimals depends on both the encounter rate and thestrength of gravitational focusing. The encounter rate grows as therelative velocity 𝑣 between the planetesimals increases, while grav-itational focusing is most effective for bodies with small relativevelocities. There is therefore an intermediate value of 𝑣 for whichthe average collision rate (cid:104) 𝜎𝑣 (cid:105) is at a minimum, which we will referto as (cid:104) 𝜎𝑣 (cid:105) . The dynamical excitation induced by a mean motionresonance can either increase or decrease the local average collisionrate, depending on the unperturbed value of (cid:104) 𝜎𝑣 (cid:105) relative to (cid:104) 𝜎𝑣 (cid:105) .In figure 5, we show the effect that the 2:1 MMR and 3:1MMR have on the local average collision rate. Each pair of pointsconnected by a line represents the average collision rate before andafter the effects of the mean-motion resonances develop. Here, therelative velocity between bodies is calculated by measuring the ec-centricity and inclination dispersion near each resonance and usingequation 7. We assume that the eccentricity and inclination disper-sions are coupled (Ida et al. 1993) and use the eccentricity dispersionas a free parameter. Because both the strength of secular forcing andthe effects of the resonant perturbations are different near the 3:1and 2:1 MMRs, the planetesimal population at each of these loca-tions in the disk experiences a qualitatively different change in thecollision rate. Near the 3:1 MMR, the unperturbed collision rate iswell above the minimum value (cid:104) 𝜎𝑣 (cid:105) and the dynamical heatingintroduced by the resonance acts to suppress it further. Near the2:1 MMR, the collision rate starts out much closer to the minimumvalue and the additional perturbations instead act to drive it higher.This explains the relative drop in the collision rate near the3:1 MMR and the relative increase near the 2:1 and 5:3 resonancesseen in figure 4. Although this could potentially serve as a usefuldiagnostic of the planetesimal size in a disc (because the size and MNRAS , 1–16 (2020) ollision rates of planetesimals Figure 4.
A PDF of the collision rate in each disc as a function of semimajoraxis, generated using a KDE with a Gaussian kernel with a FWHM of 0.02au. In semimajor axis space, prominent features appear near the 3:1 and 2:1MMRs. Near the 3:1, the collision rate exhibits a local minimum, while anenhancement appears near the 2:1. mass of the planetesimals sets the eccentricity dispersion at which (cid:104) 𝜎𝑣 (cid:105) is minimized), a direct measurement of the semimajor axes ofthe colliding planetesimals is not possible. Furthermore, as we willdemonstrate next, collisions due to bodies in resonance do not havea significant effect on the final shape of the radial dust distributionthat we predict. To construct a radial dust profile from the collision statistics, webegin by making two assumptions: (1) any dust generated by colli-sions is strongly coupled to the gas and (2) any subsequent spatialevolution of the collisional debris is insignificant. So long as bothof these assumptions are true, we can use the radial locations ofplanetesimal collisions in the plane of the disk to generate a dustprofile for each simulation.We first provide a justification for assumption (1), which is asfollows: At 3 au in the protosolar nebula, the mean free path of a gasparticle is ∼
50 cm (assuming a composition of pure hydrogen with a local volume density of 𝜌 𝑔 = × − g cm − ). This places anydust grains in the Epstein drag regime, with a stopping timescalegiven by 𝑡 𝑠 = 𝜌𝑠𝜌 𝑔 𝑣 𝑡ℎ , (11)where 𝜌 is the bulk density of the dust grain, 𝑠 is its size, and 𝑣 𝑡ℎ isthe local thermal velocity of the gas. At 3 au in the protosolar nebula, 𝑣 𝑡ℎ ∼ cm s − . Assuming a 1 mm dust grain with a 𝜌 = 2 g cm − ,the stopping time is ∼ − yr. This is orders of magnitude smallerthan the orbital timescale; therefore, we conclude that collisionallygenerated dust grains will immediately couple to the gas.Assumption (2) is less straightforward to justify and may notbe reasonable in all cases. Although the spatial redistribution of thefull-sized planetesimals and mm-sized dust grains is insignificant,the radial drift forces experienced by intermediate-sized bodies ismuch greater. Collisions between planetesimals do not immediatelyconvert the mass to dust grains. Rather, the debris will undergo acollisional cascade process, which gradually grinds the resultingbodies down to smaller sizes. A proper treatment of the collisionalcascade process is quite complicated and is well beyond the scopeof this work. However, one may be able to ignore this process if thepeak radial drift timescale is longer than the timescale for collisionsbetween these bodies.This occurs for bodies of the size at which 𝑡 𝑠 Ω = 1. In theStokes drag regime at 3 au, this occurs for bodies of 1 meter insize. Following a similar calculation in 2.4 with a 1 m body, wefind that an object of this size will drift across the 2:1 MMR ona timescale of ∼
300 years. Although we do not model collisionsbetween bodies of this size, a lower limit on the collision timescalefor these objects can be obtained by extrapolating the planetesimalcollision rate measured in our simulations down. For the e2m2simulation, we find that 125 collisions occur within the librationwidth of the 2:1 MMR from T=2,000 yr to T=5,000 yr. There are atotal of 15,000 planetesimals within this region, which means thata single planetesimal undergoes a collision every ∼ , while the geometriccollision cross section decreases by a factor of 10 . The collisiontimescale would therefore decrease to ∼ MNRAS000
300 years. Although we do not model collisionsbetween bodies of this size, a lower limit on the collision timescalefor these objects can be obtained by extrapolating the planetesimalcollision rate measured in our simulations down. For the e2m2simulation, we find that 125 collisions occur within the librationwidth of the 2:1 MMR from T=2,000 yr to T=5,000 yr. There are atotal of 15,000 planetesimals within this region, which means thata single planetesimal undergoes a collision every ∼ , while the geometriccollision cross section decreases by a factor of 10 . The collisiontimescale would therefore decrease to ∼ MNRAS000 , 1–16 (2020)
S. C. Wallace et al.
Figure 5.
The collision rate (given by equation 6) of bodies in the vicinity of the 3:1 (left) and 2:1 (right) MMRs, relative to the minimum value, as a functionof the local eccentricity dispersion. The pairs of points connected by lines show the values of the unperturbed (left-hand points in each subplot) and perturbed(right-hand points in each subplot) collision rates in the e1m2, e2m2 and e3m2 simulations. The unperturbed collision rate, along with the amount of dynamicalheating that the planetesimal population experiences determines whether the MMR acts to suppress or enhance the collision rate.
Operating under these assumptions , a map of the relativeconcentration of second-generation dust can be constructed from thecylindrical distance at which the collisions occur. Here, cylindricaldistance is defined as the separation between the central star and thelocation of collision in the ( 𝑟, 𝜃 ) plane. This is shown in figure 6.Similar to figure 4, we use a KDE with a width of 0.02 au to assemblethe collision locations into a radial distribution. Most strikingly,the bump that was present near the 2:1 resonance is no longervisible in cylindrical distance space, although it suddenly appearsagain in the e3m2 simulation. The local minimum that is visiblenear the 3:1 MMR also becomes a local maximum in the highesteccentricity case. To determine how the resonant bodies actuallycontribute to the radial collision profile, we excluded collisions thatfall between 2 . < 𝑎 < .
505 au (near the 3:1 resonance) andbetween 3 . < 𝑎 < .
35 au (near the 2:1 resonance), which isshown by the dashed curve. As in figure 4, the solid curves arenormalized such that the area underneath is equal to 1. For thedashed curves, the normalization factor is scaled according to thenumber of collisions in the entire sample, compared to the numberof collisions in the subsample, which excludes the collisions inresonance. Qualitatively, none of the bump or dip features presentare removed by making this exclusion. This suggests that the featuresnear resonance seen in semimajor axis space become smeared outin cylindrical distance space due to the large spread in eccentricityand orbital orientation.This also suggests that the prominent bumps and dips seen infigure 6 are mainly produced by bodies outside of the libration widthof the resonance, although the resonant population does make someof the features more pronounced. To further understand this, we con-sider the nonaxisymmetric structure seen in figure 1. In figure 7, wecompare the radial collision profiles to the radial structure seen in As discussed in Boley (2017), the local generation of second-generationdust is only one consideration. The dust will also evolve spatially due to drageffects, but can also evolve in size or be re-accreted onto planetesimals. Forsimplicity, we focus on local dust sourcing only while acknowledging thatsubsequent evolution could also affect the morphology. the e1m2, e2m2 and e3m2 simulations. In the lowest forced eccen-tricity case (e1m2), the pileups near the edge of the 2:1 MMR lineup with the boundaries of the dip feature seen in the radial collisionprofile. The locations of the edges of the resonances are calculatedusing equations B5 and B1. The pileups not as pronounced near theother resonances. As discussed previously, the planetesimal densityenhancements and gaps seen in polar coordinates follow the positionof the giant planet in its orbit, which is shown by the vertical line.For the e2m2 case, the density enhancements take up a much largerradial width, which smears out the corresponding collisional dustprofile. In the highest forced eccentricity case (e3m2), the densityenhancements occur well within an annulus containing the resonantregion, completely altering the radial peak-valley dust morphology.Although the nonaxisymmetric pileup of bodies near the edgesof resonances has been seen previously (Richardson et al. 2000;Tabeshian & Wiegert 2016), we could not find a satisfactory ex-planation for why it happens, which we will try to provide here.As mentioned in section 2, the circulation frequency of the criticalangle slows as one approaches the edge of the resonance from theoutside. Following the pendulum analogy, this is equivalent to ap-proaching the point where the pendulum becomes suspended in theat top dead center. When the critical angle stays relatively station-ary during an encounter, the net torque will drive the longitude ofpericenter towards the point of conjunction (Peale 1976). This is thebasic mechanism that causes isolated resonances to be stable. Insideof resonance, the critical angle librates about an equilibrium valueand this configuration is maintained. Just outside of resonance, how-ever, the critical angle slowly drifts away, with the drift directionchanging sign across the location of exact resonance. The net re-sult of this is that bodies just outside of the resonance experiencea torque from the planet that temporarily drives their pericenterstoward the current orbital position of the planet. After the encounterends, differential rotation of the disc slowly undoes the configura-tion. As the planet again passes by, adjacent planetesimals near theresonance edges temporarily align themselves in this configuration.This explains why the concentration of planetesimals appears to‘follow’ the planet. Connecting this back to the radial collision pro-files shown in figure 6, these bodies at the edges of the resonances
MNRAS , 1–16 (2020) ollision rates of planetesimals Figure 6.
Here, collisions are instead ordered by cylindrical distance fromthe star. Features near the 3:1 and 2:1 MMRs are still present, but appearqualitatively different than in semimajor axis space. At low eccentricities, adip appears around the center of the resonance. At the highest eccentricity,a bump is formed instead. The dashed lines shown the collision profile withcollisions between bodies in resonance removed. The bump and dip featuresremain qualitatively the same, which suggests that they are produced bybodies outside of the libration width of the resonance. must be responsible for the bumps and dips morphology seen be-cause excluding collisions between bodies inside of the MMRs doesnot qualitatively alter the profiles.Tabeshian & Wiegert (2016) provide a different explanation forthis phenomenon, suggesting that the nonaxisymmetric structure isthe product of the path that a resonant test particle takes in a frameco-rotating with the planet (see figure 8.4 of Murray & Dermott(1999)). In other words, the gap structure is argued to be due tothe interfaces between the low eccentricity, non-resonant and higheccentricity, resonant regions of the disc. This explanation, however,does not hold for a large collection of planetesimals, especiallywhen there is no forced eccentricity. In such a case, the orbits ofplanetesimals near resonance would be randomly aligned, and theaxisymmetric structure should become washed out. Furthermore,when the perturbing planet is on a circular orbit, this structure is stillpresent (see figure 3 of Tabeshian & Wiegert (2016) and AppendixC). The only connection between the structure of the gaps in thedisc and the path of a resonant particle in the corotating frame is that both phenomena produce the same symmetries, which dependon the particular MMR considered.With no secular forcing, this dynamical phenomena shouldproduce an underdensity in the radial surface density profile near thecenter of the resonance and an overdensity at each edge in cylindricaldistance space. When a forced eccentricity is introduced, the radiallocation of the resonance edges vary over the course of an orbit. Themaximum and minimum radial distance that the resonance edgesoccupy is indicated by the solid and dashed lines in figure 7. If the theaphelion distance of the inner edge becomes close to the periheliondistance of the outer edge, the overdense regions on each side of theresonance meet and we should expect to see the central dip featurein the collision profile disappear. This is exactly what appears to behappening in the middle panel of figure 7. As the forced eccentricityis increased further, a region forms at the center of the resonancewhere planetesimals on both sides of the resonance spend time(although not simultaneously). When this occurs, a bump featureforms in the collision profile near the center of the resonance. Thisis apparent in the bottom panel of figure 7.Although the same phenomenon appears to be happening nearthe 3:1 MMR, the pileups at the edges of the resonance are notas well defined. We attribute this to the more complex symmetry,compared to the 2:1 MMR. In addition, a number of other nearbyresonances, including the 7:2 (at 2.25 au), the 10:3 (at 2.33 au), the8:3 (at 2.71 au) and the 5:2 (at 2.82 au) also likely are contributing tothe collision profile in cylindrical distance space near the 3:1 MMR.The 3:1 MMR is also much narrower than the 2:1, which means thatany dip or bump features produced by it would require much higherresolution observations. For these reasons, we will focus on the 2:1MMR as the main diagnostic indicator.
The dynamical effects of varying the mass of the planet are some-what simpler, in that doing so does not affect the forced eccentricityof the planetesimal disc. Instead, only the width of the resonanceschanges. The width of a first order resonance scales with 𝑚 (seeequation B1), while the leading order terms in the resonant part ofthe disturbing function, which set the strength of the resonance, alsoscale as 𝑚 . For the 2:1 MMR, the dynamics near the resonance areequally sensitive to changes in eccentricity and mass.We show the polar structure of the e2m1, e2m2 and e2m3simulations in figure 8 alongside the radial collision profile. In allthree cases, the eccentricity of the perturbing planet is set to 𝑒 𝑗𝑢 𝑝 .Changes in the apocenter and pericenter distances of the edges ofthe resonances are entirely due to changes in the resonance width.For the e2m1 and e2m2 simulations, the inner apocenter and outerpericenter distances near the 2:1 MMR are quite similar and nostrong features appear in the collision profile near this region. Forthe e2m3 case, the edges of the resonances are sufficiently separatedto allow what appears to be the beginning of gap to form near thecenter of the 2:1 resonance in the collision profile. While the collisional dust profiles show radial amplitude variations,we still need to assess whether those variations could be observable.To proceed, we create a sky model for simulated ALMA observa-tions as follows. First, we use the radially-averaged collision profilefrom each simulation (see figures 7 and 8) as the template for anazimuthally symmetric disc. We then scale the size of each profile
MNRAS000
MNRAS000 , 1–16 (2020) S. C. Wallace et al.
Figure 7.
The collision profiles in figure 6 are shown alongside the positions of the planetesimals, in polar coordinates. The dashed and solid lines show thepericenter and apocenter, respectively, of bodies at each edge of the 3:1 and 2:1 MMR. In these figures, the longitude of pericenter of the planet lies at 𝜃 = , 1–16 (2020) ollision rates of planetesimals Figure 8.
Similar to figure 7, except the eccentricity of the planet is kept constant and the mass is varied. This has the effect of changing the width of theresonances, without altering the relative apo or peri distances of bodies near the edges. Except for the highest mass case, the apocenter and pericenter distancesof the inner and outer edges of the 2:1 resonance are too close together to produce much of a dip or a bump feature.MNRAS000
Similar to figure 7, except the eccentricity of the planet is kept constant and the mass is varied. This has the effect of changing the width of theresonances, without altering the relative apo or peri distances of bodies near the edges. Except for the highest mass case, the apocenter and pericenter distancesof the inner and outer edges of the 2:1 resonance are too close together to produce much of a dip or a bump feature.MNRAS000 , 1–16 (2020) S. C. Wallace et al. by a factor of ten, which places the perturbing planet at 52 au. Thisscaling is permitted by the dynamics (both the resonance widths andthe forced eccentricities are scale-free) and makes the disc compa-rable in size to the many discs that have now been observed byALMA (Huang et al. 2018). The angular size scale is then set byenvisaging a face-on disc at a distance of 100 pc.The sky model intensity (flux per pixel) is produced by in-terpolating the radial collision frequency onto a 2D Cartesian gridwith cell widths of approximately 2 mas. In doing so, we chooseto make the intensity proportional to the collision frequency (i.e.,the dust) because it is the most straightforward case. Many differentparameterizations are possible, depending on optical depth, grainsize, disc temperature, etc., but we will avoid such complicationshere, as they would affect the overall brightness profile, but not thepresence of any gaps or rings. We normalize the intensity by settingthe total disc flux density to be 100 mJy at about 350 GHz.The simulated observations are performed using the simob-serve task in CASA (McMullin et al. 2007). The disc is given anRA = h m s and a 𝛿 = − ◦ (cid:48) (cid:48)(cid:48) (i.e., the J2000 coordi-nates for TW Hydrae) and is “observed” through transit on March20, 2020. We observe the disc using configurations 8, 9, and 10at 350 GHz (cycle 6 antenna file), spending six hours on-sourcein each configuration . The visibilities are corrupted with thermalnoise, setting the precipitable water vapor to 1.5 mm. They are thencombined, imaged, and cleaned using the task tclean . Imaginguses Briggs weighting with a robustness parameter of -1.The results for the nominal simulation (e2m2) are shown infigure 9. The 3:1 and 2:1 resonances can be easily identified in thecleaned image and in the corresponding radial profile. Moreover,the features are qualitatively similar to the bright rings and darkgaps seen in actual disc profiles (ALMA Partnership et al. 2015).In figure 10, we show the azumithally-averaged radial profilesconstructed from the simulated cleaned image from all five simula-tions. In all cases, a bump or a dip feature is clearly visible at thelocation of the 3:1 and 2:1 MMRs, indicated by the dashed verticallines. As mentioned previously, the feature at the 2:1 MMR acts asthe actual diagnostic indicator, while the 3:1 resonance (along withthe gap that would presumably open at the location of the planet)is mainly useful for determining where the 2:1 MMR is actuallylocated. The simple bump vs dip structure that we expect to reveal itself inthe dust emission from colliding planetesimals in near-resonancewith a giant planet could potentially be used to place constraintson the mass and eccentricity of the planet. Near the 2:1 MMR, thecentral dip feature is only produced when bodies at the edges of theresonance stay sufficiently separated in cylindrical distance over thecourse of an orbit. This is achieved when either the resonance widthis large or the forced eccentricity is small. The bump feature, on theother hand, is produced when there is a sufficiently large amountof overlap between the apocenters and pericenters of the inner andouter edges of the resonance, respectively. This is achieved whenthe resonance is narrow and the forced eccentricity is large.In figure 11, we show the constraints that the presence of abump or a dip or dip at the 2:1 MMR places on possible values We have imagined that the TAC really likes us.
Figure 9.
Simulated ALMA observations for the nominal (e2m2) case,operating under the assumption that the dust profile closely traces the plan-etesimal collision profile. Top: The sky model based on the radial collisiondistribution. The size of the disc has been scaled by a factor of 10 and placedat a distance of 100 pc. Bottom: The cleaned image using combined obser-vations in configurations 8, 9, and 10. The circle in the bottom left indicatesthe simulated beam size. Gaps in the dust due to the 3:1 and 2:1 resonancesare visible. of ( 𝑚 𝑔 , 𝑒 𝑔 ) . The dashed lines indicate regions of parameter spacewhere qualitatively different features in the collision profile are ex-pected to form. Above the upper dashed line, the resonance is widewhile the radial excursion distance of planetesimals at the edge ofthe resonance is small, which gives rise to a dip in the collisionprofile. Below the lower dashed line, the resonance is narrow, whilethe radial excursion distance is large, allowing the apocenter andpericenter distances of the edges of the resonance to overlap andproduce a bump in the collision profile. Between the dashed lines,the absolute value of separation between the apocenter and peri-center distances of the inner and outer edges of the resonance isless than the resonance width. In this region of parameter space, nosignificant features in the collision profile near the 2:1 MMR areexpected to form. MNRAS , 1–16 (2020) ollision rates of planetesimals Figure 10.
The azimuthally-averaged radial profile based on the simulated cleaned images from all five simulations. The eccentricity of the perturbing planetincreases from left to right, while the mass of the planet increases from bottom to top. The vertical dashed lines indicate the locations of the 3:1 and 2:1resonances. As in figure 9, the size scale of the disc has been expanded by a factor of 10 and then placed at a distance of 100 pc.
Because the strength of the 2:1 resonance scales with 𝑒 𝑔 and 𝑚 𝑔 , a sufficiently low mass and low eccentricity planet will notcreate enough of a perturbation in the planetesimal disc to producea detectable bump or dip feature. The colored contours in figure 11indicate the relative strength of the resonance as these quantitiesare varied. Although the resonance strength is equally sensitive tochanges in mass and eccentricity, the peri-apocenter overlap of theresonance edges is more sensitive to changes in eccentricity. This isbecause changes to eccentricity affect both the resonance width andthe forced eccentricity. As we showed in figure 8, changing 𝑚 𝑔 bya factor of 4 produced very minimal changes to the resulting dustprofile.In order to actually identify the 2:1 MMR in the dust emission,at least one other resonance must be visible. Assuming the gravita- More specifically, the strength of the resonant term in the disturbingpotential tional field in the disc is near-Keplerian, the distance ratios betweentwo features can be used to determine period ratios, which can beused to confirm whether the features seen are indeed resonances.Although the 3:1 MMR does not appear to follow the simple bumpvs dip dichotomy described above, its presence should be marginallydetectable in all of the cases shown above.We would again like to emphasize that the dust emission pro-files in figure 10 are constructed upon the assumption that the radialdust structure closely traces the collision profile of the planetesi-mals. For this assumption to be valid, the resulting dust must remainlocally confined, which can be achieved if the grains are reaccretedonto planetesimals (Johansen et al. 2015). It is also necessary thatthe collisional cascade from which the dust grains form plays beforeradial drift has a chance to move debris away from the resonances.Although we showed with a back-of-the-envelope calculation (seesection 5) that this is plausible for a typical protoplanetary disk withproperties similar to the protosolar nebula, a more thorough treat-
MNRAS000
MNRAS000 , 1–16 (2020) S. C. Wallace et al.
Figure 11.
The presence of a dip or bump in the collision profile (andtherefore the dust profile) near the 2:1 MMR can be used to constrain themass and eccentricity of the perturbing planet. The color scale indicates therelative strength of the features, while the dashed lines indicate boundariesin parameter space where a dip or a bump will be produced at the resonance.Combinations of mass and eccentricity that fall between the dashed linescorrespond to an inner apocenter and outer pericenter separation at the edgesof the resonance that is smaller than the resonance width and therefore willnot create any significant feature in the collision profile. ment of collisional grinding and the subsequent evolution of thedebris is necessary to be able to interpret the dust emission profilespresented here as anything more than upper limits.
In this work, we have shown that mean-motion resonances with aperturbing planet produce significant local variations in the collisionrate of a planetesimal disc. In contrast to Richardson et al. (2000),we find that the more prominent interior MMRs, including the2:1, 3:1 and 5:3, all produce structure in the collision profile as afunction of semimajor axis. Furthermore, we find that a series ofdistinctly different features appear when collisions are ordered bycylindrical distance from the central star. The morphology of thesefeatures are tied to dynamical behavior of bodies near the edgesof the resonances, where the circulation frequency of the criticalangle approaches infinity. Particularly near the 2:1 MMR, we findthat a bump or dip in the collision profile forms when collisionsare ordered by cylindrical distance. The presence of one of thesetwo features depends on the peri- and apocenter distances of theedges of the resonance, relative to the libration width. Because thesequantities depend on the mass and eccentricity of the perturbingplanet, these properties of the planet are actually encoded in thecollision profile.Near the interior 2:1 MMR, we find that a distinct bump ordip feature is generated in the collision profile, depending on theproperties of the perturbing planet. The presence of one of these two structures can used to constrain the mass and eccentricity of theplanet. For a high mass, low eccentricity planet, a dip will form be-cause the edges of the resonance, where many collisions occur, staysufficiently separated. If the planet has a low enough mass (whichshrinks the size of the resonance) or is sufficiently eccentric (whichdecreases the separation between the apocenter and pericenter dis-tances of the edges of the resonance), a bump will instead form. Wetested this hypothesis for five different combinations of the planet’smass and eccentricity and found that a distinct bump or dip signa-ture is produced as long as the planet properties are sufficiently farfrom the dividing line in parameter space. This diagnostic is moreuseful for massive, eccentric planets, because the strength of theresonant perturbations scale linearly with both of these quantities.Although the planetesimal collision profile would not be di-rectly observable in a planet-forming disc, the dust from the re-sulting collisional cascade could potentially be used to trace it. Solong as the dust grains are well-coupled to the gas and radial driftdoes not have a significant effect on the morphology of the dust, thebump or dip features seen in the collision profile will also appearin the dust profile. Assuming that this is the case, we have gener-ated simulated observations of a protoplanetary disk with ALMAand showed that the bump or dip features seen in the collision pro-files should be detectable through the dust emission. From a dustemission profile with sufficiently strong radial features, one couldconstrain the properties of a perturbing planet on an exterior orbitin the following way: (1) identify the locations of two MMRs inthe disc (2) measure the spacing between them to verify that one ofthe features is indeed associated with the 2:1 MMR (3) determinewhether an under- or overdensity in the radial dust profile existsnear this location (4) use the presence of one of these two featuresto determine the region of the eccentricity-mass plane in which theplanet lies.In line with Boley (2017), the results of this work suggest that,in planet forming discs, the observed dust structure should not beinterpreted as simply the result of primordial dust that is perturbedby the planets and gas. Instead, some morphological features couldbe the product of collisional dynamics between larger bodies. Withinthe confines of the simulations presented here, we have shown thatthis second-generation dust can give rise to significant features inthe dust emission profile. For this reason, we emphasize that planet-forming discs should be thought of as existing along an evolutionarycontinuum, with the youngest discs being dominated by primordialdust and gas and the most evolved discs being comprised of mainlycollisionally generated dust.The results presented here appear to be broadly consistent withthe planetesimal dynamics seen by Tabeshian & Wiegert (2016). Intheir case, a dip feature can be seen in the azimuthally averagedsurface density profile, although they did not test a high enougheccentricity planet in any of their simulations to produce a bump.Another thing to note is that the gap morphology seen by Tabeshian& Wiegert (2016) was markedly different for the 2:1 exterior MMR.Whether this would alter the radial collision profile for the exterior,rather than the interior 2:1 resonance, is not immediately clear.On one hand, the solid surface density diminishes with distance( ∼ 𝑟 − / for the MMSN (Hayashi 1981)), which would weaken thesignal from locally generated collisional dust. On the other hand,exterior MMRs are quite effective at capturing inward drifting bod-ies (Weidenschilling & Davis 1985), which could locally enhancethe surface density well beyond the MMSN value. A full treatmentof the collisional cascade process and the radial drift of debris is asubject that we leave for future work.Compared to using the corotating gap opened by the planet MNRAS , 1–16 (2020) ollision rates of planetesimals as a diagnostic (Dobinson et al. 2013, 2016), measuring variationsin the dust emission near MMRs is much more subtle and requiresmuch higher spatial resolution and sensitivity. As we have shown,the bump vs dip feature near the 2:1 MMR is marginally detectablewith ALMA for the nearest protoplanetary discs, so long as the dustprofile closely matches the collisional profile of the planetesimals.Another complicating factor is that the inner ∼
10 au of most pro-toplanetary discs are optically thick in the sub-mm. The solution tothis problem is to instead observe at longer wavelengths where theinner disc is optically thin. This comes at the expense of much poorerresolution. However, future radio facilities like the NG- VLA are ex-pected to achieve sub-au resolution for nearby planet-forming discs(Ricci et al. 2018). In the more near-term, techniques like Gaussianprocess fitting present a promising way to recover substructure atsub-beam resolution (Jennings et al. 2020).
ACKNOWLEDGEMENTS
This work was facilitated through the use of advanced computa-tional, storage and networking infrastructure provided by the Hyaksupercomputer system at the University of Washington. We ac-knowledge the people of the Dkhw’Duw’Absh, the Duwamish Tribe,the Muckleshoot Tribe, and other tribes on whose traditional landswe have performed this work.
Software:
Astropy (Astropy Collaboration et al. 2013), CASA(McMullin et al. 2007), ChaNGa (Jetley et al. 2008; Menon et al.2015), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011),Pandas (Wes McKinney 2010), PYNBODY (Pontzen et al. 2013),sklearn (Pedregosa et al. 2011)
DATA AVAILABILITY
The data presented in this article are available on our GitHub repos-itory, at https://github.com/spencerw/planetesimal_coll_paper . REFERENCES
ALMA Partnership et al., 2015, ApJ, 808, L3Adachi I., Hayashi C., Nakazawa K., 1976, Progress of Theoretical Physics,56, 1756Andrews S. M., et al., 2016, ApJ, 820, L40Astropy Collaboration et al., 2013, A&A, 558, A33Augereau J. C., Nelson R. P., Lagrange A. M., Papaloizou J. C. B., MouilletD., 2001, A&A, 370, 447Barnes J., Hut P., 1986, Nature, 324, 446Boley A. C., 2017, ApJ, 850, 103Brouwer D., Clemence G. M., 1961, Methods of celestial mechanics. Aca-demic PressCieza L. A., et al., 2016, Nature, 535, 258Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015,MNRAS, 453, L73Dobinson J., Leinhardt Z. M., Dodson-Robinson S. E., Teanby N. A., 2013,ApJ, 777, L31Dobinson J., Leinhardt Z. M., Lines S., Carter P. J., Dodson-Robinson S. E.,Teanby N. A., 2016, ApJ, 820, 29Dong R., Li S., Chiang E., Li H., 2018, ApJ, 866, 110Gaspar A., Rieke G., 2020, Proceedings of the National Academy of Science,117, 9712Hayashi C., 1981, Progress of Theoretical Physics Supplement, 70, 35Huang J., et al., 2018, ApJ, 869, L42Hunter J. D., 2007, Computing in Science and Engineering, 9, 90 Ida S., Makino J., 1992, Icarus, 96, 107Ida S., Kokubo E., Makino J., 1993, MNRAS, 263, 875Isella A., et al., 2016, Phys. Rev. Lett., 117, 251101Jennings J., Booth R. A., Tazzari M., Rosotti G. P., Clarke C. J., 2020, arXive-prints, p. arXiv:2005.07709Jetley P., Gioachin F., Mendes C. Kale L., Quinn T., 2008, Proceedings ofIEEE International Parallel and Distributed Processing SymposiumJohansen A., Mac Low M.-M., Lacerda P., Bizzarro M., 2015, ScienceAdvances, 1, 1500109Kokubo E., Ida S., 2002, ApJ, 581, 666Leinhardt Z. M., Richardson D. C., 2005, ApJ, 625, 427Leinhardt Z. M., Dobinson J., Carter P. J., Lines S., 2015, ApJ, 806, 23Lissauer J. J., 1993, ARA&A, 31, 129Malhotra R., 1994, Physica D Nonlinear Phenomena, 77, 289McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, inShaw R. A., Hill F., Bell D. J., eds, Astronomical Society of the PacificConference Series Vol. 376, Astronomical Data Analysis Software andSystems XVI. p. 127Menon H., Wesolowski L., Zheng G. e. a., 2015, Computational Astrophysicsand Cosmology, 2, 1Murray C. D., Dermott S. F., 1999, Solar system dynamicsNesvold E. R., Kuchner M. J., 2015, ApJ, 798, 83Peale S. J., 1976, ARA&A, 14, 215Pearce T. D., Wyatt M. C., 2014, MNRAS, 443, 2541Pearce T. D., Wyatt M. C., Kennedy G. M., 2015, MNRAS, 448, 3679Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825Pérez L. M., et al., 2016, Science, 353, 1519Pontzen A., Roškar R., Stinson G., Woods R., 2013, pynbody: N-Body/SPHanalysis for python (ascl:1305.002)Ricci L., Isella A., Liu S., Li H., 2018, in Murphy E., ed., AstronomicalSociety of the Pacific Conference Series Vol. 517, Science with a NextGeneration Very Large Array. p. 147Richardson D. C., 1994, MNRAS, 269, 493Richardson D. C., Quinn T., Stadel J., Lake G., 2000, Icarus, 143, 45Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka.Tabeshian M., Wiegert P. A., 2016, ApJ, 818, 159Tabeshian M., Wiegert P. A., 2018, ApJ, 857, 3Wallace S. C., Quinn T. R., 2019, MNRAS, 489, 2159Weidenschilling S. J., Davis D. R., 1985, Icarus, 62, 16Wes McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Pro-ceedings of the 9th Python in Science Conference. pp 56 – 61,doi:10.25080/Majora-92bf1922-00aWyatt M. C., 2008, Annual Review of Astronomy and Astrophysics, 46, 339Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., HolmesE. K., Piña R. K., 1999, ApJ, 527, 918van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Scienceand Engineering, 13, 22
APPENDIX A: TIMESCALE FOR SECULAR FORCING
The timescale for secular forcing is given by the period at whichthe free eccentricity of a body cycles in the complex plane. Thetime-dependent imaginary and real components of the eccentricityof a planetesimal are given by Murray & Dermott (1999) (see ch.7, pg. 282, eq. 7.51) ℎ = 𝑒 𝑠𝑖𝑛 ( 𝑔 𝑡 + 𝛽 ) + 𝑒 𝑠𝑖𝑛 ( 𝑔 𝑡 + 𝛽 ) (A1) 𝑘 = 𝑒 𝑐𝑜𝑠 ( 𝑔 𝑡 + 𝛽 ) + 𝑒 𝑐𝑜𝑠 ( 𝑔 𝑡 + 𝛽 ) , where ( 𝑒 , 𝑒 ) and ( 𝑔 , 𝑔 ) are the corresponding eigenvectorcomponents and eigenvalues of a matrix, which describes the secularinteraction between two bodies, given by Murray & Dermott (1999)(see ch. 7, pg. 276, eq. 7.9 and 7.10, with 𝛼 = 𝛼 and ¯ 𝛼 = MNRAS000
The timescale for secular forcing is given by the period at whichthe free eccentricity of a body cycles in the complex plane. Thetime-dependent imaginary and real components of the eccentricityof a planetesimal are given by Murray & Dermott (1999) (see ch.7, pg. 282, eq. 7.51) ℎ = 𝑒 𝑠𝑖𝑛 ( 𝑔 𝑡 + 𝛽 ) + 𝑒 𝑠𝑖𝑛 ( 𝑔 𝑡 + 𝛽 ) (A1) 𝑘 = 𝑒 𝑐𝑜𝑠 ( 𝑔 𝑡 + 𝛽 ) + 𝑒 𝑐𝑜𝑠 ( 𝑔 𝑡 + 𝛽 ) , where ( 𝑒 , 𝑒 ) and ( 𝑔 , 𝑔 ) are the corresponding eigenvectorcomponents and eigenvalues of a matrix, which describes the secularinteraction between two bodies, given by Murray & Dermott (1999)(see ch. 7, pg. 276, eq. 7.9 and 7.10, with 𝛼 = 𝛼 and ¯ 𝛼 = MNRAS000 , 1–16 (2020) S. C. Wallace et al. 𝐴 𝑗 𝑗 = 𝑛 𝑗 𝑚 𝑘 𝑚 𝑐 + 𝑚 𝑗 𝛼𝑏 ( ) / ( 𝛼 ) (A2) 𝐴 𝑗𝑘 = − 𝑛 𝑗 𝑚 𝑘 𝑚 𝑐 + 𝑚 𝑗 𝛼𝑏 ( ) / ( 𝛼 ) . Here, 𝑛 is the mean motion of a body and 𝑚 is its mass. In the case ofa massless test particle being perturbed by a giant planet, 𝑚 = 𝑚 𝑔 and 𝑚 =
0. The frequency of the secular forcing is simply thedifference of the two eigenvalue frequencies 𝑔 − 𝑔 . APPENDIX B: LIBRATION IN MEAN MOTIONRESONANCE
For second order and higher MMRs, the maximum libration widthis given by Murray & Dermott (1999) (ch. 8, pg. 338, eq. 8.58) as 𝛿𝑎𝑎 = ± (cid:18) | 𝐶 𝑟 | 𝑛 𝑒 | 𝑞 | (cid:19) / , (B1)where 𝑛 , 𝑎 and 𝑒 are the mean motion, semimajor axis and ec-centricity of the unperturbed body. 𝐶 𝑟 is a constant defined by 𝑚 (cid:48) / 𝑚 𝑐 𝑛𝛼 𝑓 𝑑 ( 𝛼 ) , where 𝛼 = 𝑎 / 𝑎 (cid:48) and 𝑓 𝑑 is the resonant part of thedisturbing function. For an interior second-order resonance, 𝑓 𝑑 ( 𝛼 ) = (cid:104) − ( 𝑝 + 𝑞 ) + ( 𝑝 + 𝑞 ) − 𝛼𝐷 (B2) + ( 𝑝 + 𝑞 ) 𝛼𝐷 + 𝛼 𝐷 (cid:105) 𝑏 𝑝 + 𝑞 / , (see Murray & Dermott (1999) table 8.1, third row, with 𝑗 = 𝑝 + 𝑞 )where 𝐷𝑏 𝑗𝑠 is the first derivative with respect to 𝛼 of the Laplace co-efficient defined in equation 2. This is given by Brouwer & Clemence(1961) as 𝑑𝑏 𝑗𝑠 𝑑𝛼 = 𝑠 (cid:16) 𝑏 𝑗 − 𝑠 + − 𝛼𝑏 𝑗𝑠 + + 𝑏 𝑗 + 𝑠 + (cid:17) . (B3)For a first-order resonance, the associated disturbing functionterms are slightly simpler such that 𝑓 𝑑 ( 𝛼 ) = −( 𝑝 + 𝑞 ) 𝑏 𝑝 + 𝑞 / − 𝛼 𝐷𝑏 𝑝 + 𝑞 / , (B4)(see Murray & Dermott (1999) table 8., first row) and there is anadditional contribution to the motion of the critical angle by thepericenter precession rate such that 𝛿𝑎𝑎 = ± (cid:18) | 𝐶 𝑟 | 𝑛 𝑒 (cid:19) / (cid:18) + 𝑝 𝑒 | 𝐶 𝑟 | 𝑛 (cid:19) / − 𝑝𝑒 | 𝐶 𝑟 | 𝑛 (B5)(see Murray & Dermott (1999) ch. 8, pg. 340, eq. 8.76, with 𝑗 = 𝑝 ).In mean-motion resonance, the orbital elements of a perturbedbody, along with the critical angle described by equation 5, willexhibit oscillations. In the circular, restricted three-body case, thesesmall amplitude oscillations exhibit a characteristic frequency, givenby (see Murray & Dermott (1999) ch. 8, pg. 332, eq. 8.47, with 𝑗 = − 𝑝 and 𝑗 = − 𝑞 ) 𝜔 = 𝑝 𝐶 𝑟 𝑛𝑒 |− 𝑞 | . (B6) APPENDIX C: NONAXISYMMETRIC STRUCTURE INBOLEY 2017
As mentioned in section 4.1, the nonaxisymmetric gap structureseen in figure 1 is present in a number of other studies, includingRichardson et al. (2000) and Tabeshian & Wiegert (2016). Onewould expect the same dynamics to be at work in the simulationof Boley (2017), although the structure is not obviously visible infigure 3 of their work.Here, we present a closer look at the positions of the plan-etesimals at the end of the Boley (2017) simulation. In figure C1,we show the planetesimal positions in polar, rather than Cartesiancoordinates as originally presented in the paper. Two planets oncircular orbits are present here. The first has mass of 17 𝑀 ⊕ andlies at 32.3 au. The second planet has a mass of 50 𝑀 ⊕ and lies at68.8 au. A horseshoe shaped collection of planetesimals is visibleto be corotating with each planet. Nonaxisymmetric gaps are alsovisible at the locations of first order MMRs with the more massiveplanet, which are indicated with arrows.Similar to figure 7, a gap is present at the same orbital phaseas the planet at the 2:1 MMR, along with a second gap on theopposite side of the disc. Gaps near the 4:3 and 3:2 MMRs are alsopresent in phase with the planet, although the number and spacingbetween additional gaps appears to be different for each resonance.In addition, the gap structure associated with the exterior MMRsappears markedly different than that of each corresponding interiorMMR, which suggests that the structure is related to the symmetryof the path that a planetesimal in a frame corotating with the planettakes. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–16 (2020) ollision rates of planetesimals Figure C1.
Positions of planetesimals at the end of the simulation presented in Boley (2017) in polar coordinates. The dashed line marks the location of a 50 𝑀 ⊕ planet embedded in the disc. First order MMRs with the planet are indicated with arrows. A second planet is also present at 32.3 au, although it is notmassive enough to induce any significant structure in the disc near associated resonances.MNRAS000