Predictions for Shepherding Planets in Scattered Light Images of Debris Disks
aa r X i v : . [ a s t r o - ph . E P ] N ov Draft version July 16, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
PREDICTIONS FOR SHEPHERDING PLANETS IN SCATTERED LIGHT IMAGES OF DEBRIS DISKS
Timothy J. Rodigas , Renu Malhotra , Philip M. Hinz Draft version July 16, 2018
ABSTRACTPlanets can affect debris disk structure by creating gaps, sharp edges, warps, and other potentiallyobservable signatures. However, there is currently no simple way for observers to deduce a disk-shepherding planet’s properties from the observed features of the disk. Here we present a singleequation that relates a shepherding planet’s maximum mass to the debris ring’s observed width inscattered light, along with a procedure to estimate the planet’s eccentricity and minimum semimajoraxis. We accomplish this by performing dynamical N-body simulations of model systems containinga star, a single planet, and a disk of parent bodies and dust grains to determine the resulting debrisdisk properties over a wide range of input parameters. We find that the relationship between planetmass and debris disk width is linear, with increasing planet mass producing broader debris rings. Weapply our methods to five imaged debris rings to constrain the putative planet masses and orbits ineach system. Observers can use our empirically-derived equation as a guide for future direct imagingsearches for planets in debris disk systems. In the fortuitous case of an imaged planet orbiting interiorto an imaged disk, the planet’s maximum mass can be estimated independent of atmospheric models.
Subject headings: methods: N-body simulations — circumstellar matter — planetary systems INTRODUCTION
Circumstellar dust has been detected aroundseveral hundred main sequence stars, and morethan 30 such systems have now been spatially re-solved at visible, near-, and far-infrared wavelengths(http://circumstellardisks.org). In these so-called debrisdisk systems, the dust is inferred to be generated bycollisions in planetesimal belts that are dynamicallyshepherded by planets, analogous to the solar system’sasteroid and Kuiper belts. In a few systems, notably β Pictoris, HR 8799, and perhaps Fomalhaut, spa-tially resolved debris disk structure as well as one ormore planets have been detected (Lagrange et al. 2010;Marois et al. 2010; Kalas et al. 2008; Currie et al. 2012a;Galicher et al. 2012; Kalas et al. 2013), allowing moredetailed study of planet-planetesimal-disk interactionsand of planetary system formation and evolution.Currently most spatially-resolved debris disks do notalso have associated detections of planets, but many de-bris disks do show signs of shepherding planets, such asclumps, gaps, and sharp edges. Indeed the presence ofa planet in the Fomalhaut system has been predictedbased on the resolved disk’s properties (Quillen 2006),and Chiang et al. (2009) (hereafter C09) have carried outdetailed dynamical modeling to constrain the planet’smass to < J by consideration of the debris disk’s ob-served properties. A faint object scattering star light hasrecently been imaged in the system (Kalas et al. 2008;Currie et al. 2012a; Galicher et al. 2012), but its true na-ture and origin remain ambiguous (Kalas et al. 2013).The Fomalhaut dynamical predictions were very use- Steward Observatory, The University of Arizona, 933N. Cherry Ave., Tucson, AZ 85721, USA; email: [email protected] Carnegie Postdoctoral Fellow; Department of TerrestrialMagnetism, Carnegie Institute of Washington, 5241 BroadBranch Road, NW, Washington, DC 20015, USA Lunar and Planetary Laboratory, University of Arizona,Kuiper Space Sciences Building, Tucson, AZ 85721, USA ful, but they were customized to a specific debris disk.What predictions can be made in the general case? Sev-eral studies have examined how planet mass affects debrisdisk structure (Wyatt et al. 1999a; Kuchner & Holman2003; Quillen 2006; Chiang et al. 2009; Thebault et al.2012; Nesvold et al. 2013), but these are not readilyadapted for the inverse problem, namely estimating theplanet mass and orbit from a given set of debris diskobservations.The goal of the present paper is to provide observerswith a “rough guide” for estimating the mass and orbitof a putative shepherding planet from scattered light ob-servations of a debris disk. The procedure we presentis obtained as follows: we carry out a suite of N-bodynumerical integrations consisting of a single planet inter-acting with an exterior ring of dust-producing planetes-imals; the simulations cover a range of system parame-ters, and we obtain a suite of simulated debris disks. Wecalculate the optical depth profile (a proxy for surfacebrightness) and measure the normalized full-width half-maximum (FWHM) of each simulated disk. We thenexamine the relationships between the FWHM and thesimulation input parameters. Specifically, we provide ob-servers with a simple procedure to estimate the maxi-mum mass of a putative shepherding planet, its orbitalsemimajor axis, and its eccentricity from surface bright-ness profiles of scattered light images of debris disks. InSection 2, we describe our methods. In Section 3, wepresent our results and predictions for planets in five im-aged debris rings. In Section 4, we discuss the implica-tions of our results and outline the observer’s procedurefor estimating a planet’s maximum mass and orbit fromscattered light images of debris disks. In Section 5 wesummarize our main findings. METHODS
We adopt the hypothesis that the scattered light ofthe debris disk arises from stellar radiation scatteredby dust grains produced in an underlying ring of dust- Rodigas, Malhotra, & Hinzproducing planetesimals–“parent bodies”–shepherded bya nearby perturbing planet orbiting interior to its in-ner edge. As in C09, we first numerically integratethe (massless) parent bodies to produce an ensembleof stable orbits. The remaining stable parent bodiesare then assumed to “release” dust grains of a pre-scribed size distribution. The dust-producing parentbodies are large enough that stellar radiation pressurecan be ignored for them, which results in disks thatare intrinsically narrower (C09; Thebault et al. (2012);Nesvold et al. (2013); Boley et al. (2012)). However, ra-diation pressure is not negligible for the much smallerdust grains since it spreads them onto wide orbits (C09).Because radiation pressure causes a radial accelerationthat is inversely proportional to astrocentric distance,it is simple to model it as a fraction, β , of the stellargravitational acceleration, where β is a function of par-ticle size and density (Wyatt & Whipple 1950). In ournumerical simulations, we only include bound particles( β < . β -modified Jacobi constant). To il-lustrate, a simulation of 10,000 particles interacting witha star and a perturbing planet for 1000 orbits of theplanet, and using an integration step size of 5% of theplanet’s orbital period, takes only ∼ . The entire suite of 160independent N-body simulations required to vary all theinput parameters in this study takes only a few hours ofwall clock time to complete. Simulation parameters
The parameters of interest in this study are: the inneredge of the parent body disk, a inner ; the initial widthof parent body disk, w pb ; the initial eccentricity of theparent body disk, e disk,i ; and the longitude of periastronof the parent body disk, ̟ pb . The perturbing planet’sparameters are its mass, m p , its orbital semimajor axis,a p , eccentricity, e p , and longitude of periastron, ̟ p .In spatially resolved scattered light images of debrisdisks, if the disks are not too inclined relative to our lineof sight, we can typically determine the following diskparameters, or observables: the semimajor axis, a peak ,the eccentricity, e disk , and longitude of periastron, ̟ disk ,of the deprojected disk. In this study, we will also uti-lize the normalized FWHM (nFWHM) of the (radiationdilution-corrected) surface brightness profile of the de-projected disk, defined here asnFWHM = a out / − a in / a peak , (1)where a out / and a in / are the outer and inner semimajor Macbook Pro, 4 GB memory, 2.66 GHz Intel Core 2 Duo axis locations where the deprojected disk surface bright-ness is half of the peak surface brightness . Other workshave utilized the sharpness of the inner edge of the diskas the key observable indicating the presence of a disk-shepherding planet. However, the inner edge sharpnesssuffers from a degeneracy in planet mass/semimajor axis(ie, high-mass planets far away from the disk can pro-duce the same sharpness as low-mass planets close to thedisk; see Fig. 3 in C09). The width of the disk, on theother hand, is much less affected by this degeneracy, aswill be demonstrated in Section 3. See Table 1 for a listand description of the relevant parameters utilized in thiswork. Initial conditions
We adopt units whereby the stellar mass M ∗ and theuniversal constant of gravitation, G , are unity, and theinner edge of the initial parent body disk, a inner isadopted as the unit of length. In these units, a par-ticle orbiting at a inner has a period of 2 π . We simu-late systems with planet mass µ in the range (0.3, 1.0,3.0, 10.0) µ J , where µ J = 9 . × − is the mass ratioof Jupiter relative to our Sun. Although smaller planetmasses are likely present in debris disks, we do not sim-ulate these cases for two reasons. First, Thebault et al.(2012) showed that over a wide range of optical depthsand for µ < . µ J , dust grain collisions can wash outobservable effects on debris disks. Second, direct imag-ing instrumentation is currently not capable of detectingplanets less massive than ∼ a few M J .The parent bodies’ semimajor axes are all initializedat a inner (ie, initially infinitesimally-narrow disks with w pb = 0). We also tested disks with non-zero initialwidths but found that these cases introduced a degener-acy with planet mass/semimajor axis combinations. Forexample, consider an imaged debris disk with a FWHM= 15%. This disk could have been broadened from aninitial width of 10% due to interactions with a distant,high-mass planet; or it could have been broadened froman initial width of 14% by a very low-mass, very nearbyplanet. An observer has no way of knowing which planetproduced the observed disk without knowing the initialparent body disk width, which cannot be measured. Ifwe instead consider initially infinitesimally-narrow rings,then for a given observed disk width, we can at least con-strain the maximum mass of the putative planet, even ifnature is unlikely to produce such initial parent bodydisks.The star, planet, and parent bodies all orbit in thesame plane so that their relative inclinations are zero.This assumption is made for simplicity and should notseriously change the results since most debris disks arethought to be flat (aspect ratio, defined as disk heightdivided by disk radius, typically on the order of a fewpercent).We consider five values of the initial eccentricity of theparent body disk, e disk,i : (0, 0.05, 0.1, 0.15, 0.2). Thesevalues adequately span known debris disk eccentricities.The planet and the parent bodies all have the same ini-tial longitude of periastron, which we initialize to zero.The parent bodies have random initial mean anomaliesuniformly spaced between 0 and 2 π . The right-hand side of Eq. 1 is often expressed as ∆
R/R . TABLE 1Simulation Parameters
Parameter Value Description µ/µ J (0.3, 1.0, 3.0, 10.0) planet mass ratio / Jupiter’s mass ratioe p (0.0, 0.05, 0.10, 0.15, 0.20) planet eccentricitya p see Table 2 planet semimajor axis ̟ p disk,i (0.0, 0.05, 0.10, 0.15, 0.20) initial parent body disk eccentricitya inner ̟ pb w pb β (0.0, 0.00625, 0.0125, 0.025, radiation pressure force / gravity0.05, 0.1, 0.2, 0.4)a peak output optical depth profile peak semimajor axisa in / output optical depth profile inner half-peak semimajor axisa out / output optical depth profile outer half-peak semimajor axise disk,f output final parent body disk eccentricitynFWHM output normalized optical depth profile FWHM We set the planet’s eccentricity to be equal to the par-ent body disk’s initial eccentricity. While other studieshave used the forced eccentricity relationship betweenthe planet and the dust (e.g., Quillen (2006) and C09),this relationship is derived from linear secular theory,which is only valid for small mass ratios ( µ . − ) andlow eccentricities (Mustill & Wyatt 2009). Since the par-ent bodies we simulate are effectively indestructible overthe simulation timescale (see the Appendix), assumingthat they can acquire forced eccentricities from nearbyperturbing planets leading up to the start of our simu-lations is reasonable. Moreover the relationship betweenthe planet’s and disk’s eccentricity that we simulate ismerely the initial relationship; the eccentricities of theparent bodies evolve over time. Numerical determination of planet semimajor axis
The initial semimajor axis of the planet is determinedby means of a bootstrap procedure in which we numer-ically determine the inner edge of stable orbits exteriorto the planet’s orbit for the ranges of planet mass andeccentricity of interest. Although we could have madeuse of published formulas for the chaotic zone of a planet(e.g., Wisdom (1980); C09; Mustill & Wyatt (2012)), theresults in the literature do not adequately cover the rangeof planet mass and eccentricity of interest in our study.Therefore we determine the semimajor axis of the planetrelative to the inner stable edge of the parent body diskas follows: we place a planet of a given mass at a inner andplace parent bodies at discrete semimajor axis locationsbeyond a inner , starting far from the planet. The par-ent bodies have the same eccentricity as the planet andhave random mean anomalies between 0-2 π . The sys-tem is then integrated for 1000 orbits of the planet. Inthis simulation and all others described below, we use anintegration step size of 5% of the orbital period of a par-ticle with semimajor axis = a inner . Parent bodies thatapproach within the planet’s Hill radius or cross within0.1 a inner of the star or beyond 100 a inner are discarded.The width of the unstable zone is determined as the dis-tance between the planet and the closest semimajor axisat which at least 90% of the parent bodies remain at theend of the integration. The final planet semimajor axislocations are reported in Table 2.The large separations between the planet and disk ( & µ / from C09, or 1.8 µ / e / p from Mustill & Wyatt (2012))because we are testing larger planet masses and largereccentricities. Additionally, the relationship between theplanet’s and disk’s eccentricity in our study is differentfrom previous works that typically use the forced eccen-tricity relationship. We also used a different stabilitycriterion ( >
90% of particles must survive the integra-tion) to determine the chaotic zone widths, which maybe more stringent than previous works (e.g., C09).
TABLE 2Numerically determined a p /a inner e disk,i µ/µ J = 0.3 µ/µ J = 1 µ/µ J = 3 µ/µ J = 100 0.86 0.80 0.71 0.590.05 0.86 0.80 0.70 0.550.10 0.86 0.80 0.69 0.520.15 0.86 0.79 0.69 0.490.20 0.84 0.79 0.67 0.47 From Table 2, the initial eccentricity of the disk widensthe unstable zone only for the largest mass ratios. How-ever, the initial eccentricities of debris disks cannotbe measured. Therefore we cannot compute a masterchaotic zone equation from the values in Table 2. In-stead, because we are attempting to constrain the max-imum mass of a disk-shepherding planet in this study,we should only constrain the minimum semimajor axisof a putative planet. Therefore we adopt as our unstablezone equation the power-law fit to the e disk,i = 0 .
20 val-ues, since these correspond to the widest chaotic zonesand the planets that are the farthest from the disk’s inneredge: a p = a inner / (cid:0) . µ . (cid:1) . (2)Note that this equation is not used directly in our sim-ulations, as the planet semimajor axes are explicitly de-termined in Table 2. Rather, this equation is to be used Rodigas, Malhotra, & Hinzby observers when conducting the procedure outlined inSection 4.To review, we have two free parameters that will bevaried in the simulations: µ and e disk,i . a p is numeri-cally determined, e p = e disk,i , and ̟ p = ̟ pb = 0. Thesimulations will produce two outputs: nFWHM, calcu-lated using Eq. 1, and e disk,f , the final eccentricity ofthe simulated disks. Procedure
With all parent body and planet input parameters de-termined, we integrate the system of the star, planet, and5000 parent bodies for 1000 orbital periods of a particleat a inner . After 1000 orbits, the parent bodies “release”dust grains that have the same positions and velocitiesas their parents, as in C09. The dust grain orbits arethen numerically integrated, accounting for the effects ofradiation pressure by multiplying the stellar mass M ∗ by1 − β .As in C09, we simulate 8 different values of β ; including β = 0, these are: (0.0, 0.00625, 0.0125, 0.025, 0.05, 0.1,0.2, 0.4). We integrate each system with a given β for 100orbit periods of a particle at a inner . The length of thisintegration is approximately the collisional lifetime of adust particle in the typical debris disks we are simulating(see the Appendix, Eq. A8). Optical depth profiles
To obtain the optical depth profiles of the simulateddebris disks, we follow the procedure outlined in C09.We take each surviving dust particle’s final Cartesianpositions and velocities, construct a grid of concentricellipses with eccentricity given by the final eccentricityof the parent body disk, e disk,f , and ellipse center givenby a disk e disk,f , and count the total number of survivingparticles in a given ellipse. Here, e disk,f and a disk are thelocations of the peak of the eccentricity and semimajordistributions of the surviving parent bodies, respectively.To increase the signal-to-noise (S/N), we “spread”each dust particle out along its orbit (“Gaussian wiremethod”) as in C09. In this method, a surviving particleis cloned and placed at discrete locations along its orbit;the orbital elements are determined from the positionand velocity of the particle at the end of the integrationand account for the effects of radiation pressure. Thenumber of particle clones generated per orbit is chosen sothat the total number of particles equals 10 . For exam-ple, if 5000 particles survive, then each particle would be“spread” along its orbit at 200 locations equally spacedin true anomaly. For the β > inner andthe last at 1.83 a inner . This results in an optical depthresolution of 0.02 a inner per ellipse. This is lower thanthe resolution used in C09 (200 ellipses), but is well-matched to current observations. For example, if a inner = 50-100 AU (most debris disks reside at these sepa-rations), and 1 ellipse = 1 resolution element, then theresolution is 1-2 AU. This is comparable to the typical a disk /a inner no r m a li z ed β p r o f il e β = 0 β = 0.00625 β = 0.0125 β = 0.025 β = 0.05 β = 0.1 β = 0.2 β = 0.4 (a) a disk /a inner no r m a li z ed op t i c a l dep t h µ / µ J = 0.3 µ / µ J = 1 µ / µ J = 3 µ / µ J = 10 (b) Fig. 1.—
Top : β profiles for µ = 0 . µ J and initial e disk,i = 0.10.The profiles are very similar to what C09 observe for µ = 0 . µ J ,e p = 0.12, and parent body disk width = 10%. Bottom : β -summedoptical depth profiles for different µ values, all with e disk,i = 0.10.As the mass ratio increases, the profiles spread out, as is observedby C09 for their Fomalhaut simulations. resolution achieved by HST for debris disks 50-100 pcfrom Earth.For a given planet mass and disk eccentricity, we pro-duce the final optical depth profile τ ⊥ in the same man-ner as C09, by linearly combining the 8 different opticaldepth profiles for each β : τ ⊥ = X β =0 N β max ( N . ) max ( N β ) (cid:18) β . (cid:19) q − + (3) N max ( N . ) max ( N ) (1 + √ , where N β refers to the optical depth profile for a specific β , q is the differential power-law index assuming a col-lisional cascade in the disk and here assumed to be 3.5(Dohnanyi 1969), and the 1 + √ Control simulation
Before running a full suite of N-body simulations, weverified that our simulations were producing results sim-ilar to those obtained by C09. While their input param-eters (mass ratio, planet/disk eccentricity, initial parentbody disk width) differ from ours, the general set up andmethodology are very similar. Therefore we should ex-pect to see similar results for similar inputs.Fig. 1a shows our β profiles for the specific case of µ/µ J = 0.3 and initial disk eccentricity = 0.10. The pro-files are very similar to Fig. 2 in C09. This gives us con-fidence that our simulations will yield accurate results,despite our differences relative to C09.We also verified that the perturbing planet was havingthe expected effect on the dust particles, namely spread-ing them out, resulting in wider optical depth profileswith increasing mass (as was seen by C09). Fig. 1bshows such an example simulation for an initial disk ec-centricity of 0.10. The expected behavior is observed,again validating the methods and parameters chosen forthe simulations. While the disk width increases onlymarginally for small mass ratios, the difference is readilyevident when compared to the 10 µ/µ J model. This issatisfactory for the purposes of our study: discriminat-ing between low-mass ( ∼ µ/µ J ) and high-mass ( ∼ µ/µ J ) planets in a given debris disk system. RESULTS
After running the full suite of simulations, we mea-sured a out / , a in / , and a peak for each β -summed opticaldepth profile (as in Fig. 1b) and computed the normal-ized FWHM of each modeled disk using Eq. 1. We thenexamined the relationships between this output and theplanet mass ratio and final disk eccentricity. We cannotuse the initial disk eccentricity as an independent vari-able because it can increase or decrease from its originalvalue, depending on the mass of the perturbing planet(see Fig. 2). Therefore we first examine the relationshipbetween each disk’s normalized FWHM and its final diskeccentricity (see Fig. 3). From Fig. 3, there is no strongcorrelation between these two variables. This means thatdisk eccentricity may not be an indicator of a nearbymassive planet (neglecting the possible dynamical inter-actions that could have excited the disk into an eccentricstate in the first place).Since the final disk eccentricity has little effect on thedisk’s width, we now examine the relationship betweenthe normalized disk FWHM and planet mass ratio (seeFig. 4). Evidently increasing the mass ratio of the per-turbing planet causes an increase in the range of possibledisk FWHM values.Fitting a line to data yields an equation which canbe used to estimate the maximum mass of a perturbingplanet in a debris disk. To estimate the uncertainties inthe slope and y-intercept, we fit all possible combinationsof points and set the uncertainties as the differences be-tween the minimum and maximum slope and y-interceptvalues. We then propagated these errors to obtain a to-tal error for each predicted mass. Eq. 4 shows our linearfit to the data, along with the uncertainties in the slopeand y-intercept:nFWHM = (0 . ± . µ/µ J + (0 . ± . . (4) initial disk eccentricity f i na l d i sk e cc en t r i c i t y Fig. 2.—
Final parent body disk eccentricity as a function ofinitial parent body disk eccentricity for the various input planetmasses. The points represent the averages of the final disk eccen-tricity values across mass, and the error bars represent the fullrange in values. The dashed line denotes perfect agreement be-tween initial and final eccentricity. The introduction of a massiveperturbing planet into the system can result in both an increaseand a decrease in the disk’s eccentricity, depending on the massand initial eccentricity. final disk eccentricity no r m a li z ed F W H M + c on s t an t µ / µ J = 0.3 µ / µ J = 1 µ / µ J = 3 µ / µ J = 10 Fig. 3.—
Normalized disk FWHM as a function of final disk ec-centricity for the various input planet masses. The different colorsare each associated with a different mass ratio planet, and constantterms have been added to the FWHM values for graphical clarity.The dashed lines are fits to the data for each mass. There is nostrong correlation between a debris disk’s width and its eccentric-ity. This implies that disk eccentricity may not be a good indicatorof a nearby perturbing planet.
Interestingly, our fitted slope (0.019) is close to the slopefrom C09 ( ∼ mass ratio ( µ / µ J ) no r m a li z ed F W H M This workC09
Fig. 4.—
Normalized disk FWHM as a function of planet massratio. The black points represent the averages of nFWHM valuesat each mass for the various final disk eccentricities, and the errorbars represent the full range in nFWHM values for each mass. Thedashed line represents a linear fit to the points. The blue pointsand error bars represent equivalent values taken from Fig. 3 inC09. Increasing the mass ratio of the perturbing planet causes anincrease in the range of possible disk FWHM values. A similarrelationship was seen in C09 for Fomalhaut’s disk. tered light normalized FWHM: m p / M J = (cid:18) nFWHM − (0 . ± . . ± . (cid:19) (cid:18) M ∗ M ⊙ (cid:19) . (5) Predictions for resolved debris rings
We now use our empirically determined relationshipsto constrain the mass and orbit of shepherding planetsin five bright ring-like debris disks.
Fomalhaut
Fomalhaut is a very nearby A star with a debris diskdetected in scattered light (Kalas et al. 2005). Recently apoint source has been imaged orbiting interior to the ring(Kalas et al. 2008; Currie et al. 2012b; Galicher et al.2012). While originally posited as a planet shepherd-ing the ring, it appears to be primarily scattering thestar’s light, its eccentricity is likely to be very large,and its orbit may intersect with the plane of the disk(Kalas et al. 2013)–all of which call into question its na-ture as a disk-shepherding planet. Nonetheless it is usefulto estimate the maximum mass of a planet shepherdingthe disk, since such a planet may still exist in the systemand future observations will seek to detect it. Since thedisk’s eccentricity is not a good predictor of planet mass,we need only the disk’s deprojected normalized FWHM.This is ∼ ⊙ , wefind a maximum planet mass of ∼ ± . J . Theplanet’s eccentricity would be equal to the disk’s (0.11),and its semimajor axis would be &
85 AU from Eq. 2.Taking into account the inclination to the system, theminimum orbit-averaged projected separation would be ∼ . ′′
4. Recent imaging studies have already ruled outplanets more massive than ∼ J at these distances(Janson et al. 2012; Currie et al. 2013). Therefore if aplanet is currently shepherding the debris ring, it mustbe low-mass. HR 4796A
HR 4796A has a bright debris disk that has been re-solved at many wavelengths (e.g., Schneider et al. (2009);Thalmann et al. (2011); Lagrange et al. (2012)). In scat-tered light, the disk appears as a narrow ring, with alarge central gap between the disk and the star. Thegap, the sharp inner and outer edges, and a smalloffset from the star are cited as evidence for a per-turbing planet (Schneider et al. 2009; Thalmann et al.2011; Wyatt et al. 1999a). To date no planet has beendetected, though Lagrange et al. (2012) ruled out 3.5M J planets beyond 0 . ′′ ∼ ⊙ , we estimate the maximum mass of a perturbingplanet would be ∼ ± . J . The planet’s eccentricitywould be equal to the disk’s ( ∼ &
45 AU. Its orbit-averaged minimumprojected separation would be 0 . ′′
14. Given the largeshepherding planet mass this system can tolerate, alongwith the possibility that such a planet may have beenmissed by recent imaging campaigns (Thalmann et al.2011; Lagrange et al. 2012) due to its possible small pro-jected separation, we advocate continued high-contrastimaging of this system in the coming years.
HD 207129
HD 207129 is a Sun-like star that has a large, faintdebris disk, recently resolved in scattered light byKrist et al. (2010) with HST. The disk has a normal-ized FWHM of ∼ ⊙ (Krist et al. 2010), we obtain using Eq. 4 a maximumplanet mass 4.2 ± . J . Since the eccentricity of thedisk is only constrained to be < .
08, the same con-straint applies to the planet’s eccentricity. Its semimajoraxis would be &
92 AU, corresponding to a minimumorbit-averaged projected separation on the sky of ∼ . ′′ ∼ HD 202628
HD 202628 is a Sun-like star with a very wide diskrecently resolved by HST (Krist et al. 2012). The diskhas a normalized FWHM of ∼ ∼ ∼ ± . J . Its eccentricity would be 0.18, and its min-imum semimajor axis would be ∼
71 AU, correspond-ing to a minimum orbit-averaged projected separationof ∼ . ′′
2. Despite the star’s likely old age (2.3 Gyr;Krist et al. (2012)), its broad disk can tolerate a verymassive planet that would likely still be detectable bycurrent direct imaging technology. Therefore we advo-cate high-contrast imaging of this system.
HD 181327
HD 181327 is a Sun-like star with a large, bright debrisdisk resolved by HST (Schneider et al. 2006). The mostcurrent analysis of the resolved images reveals that the
TABLE 3Predicted Masses and Orbits
Star m p /M J a p /AU proj. sep. ( ′′ ) ∗ e p Fomalhaut < . ± . > > . . < . ± . > > .
14 0 . < . ± . > > . < . < . ± . > > . . < . ± . > > .
55 0 ∗ Orbit-averaged disk has a normalized FWHM of 0.32 (Lebreton et al.2012) and an eccentricity consistent with zero (C. Stark,private communication). For a stellar mass of 1.36 M ⊙ (Lebreton et al. 2012), the shepherding planet’s maxi-mum mass would be ∼ ± . J . Its eccentricitywould be ∼ ∼
35 AU, corresponding to a minimum orbit-averagedprojected separation of ∼ . ′′
55. Wahhaj et al. (2013) im-aged this star as part of the Gemini NICI Planet-FindingCampaign and ruled out planets more massive than ∼ J beyond 0 . ′′
36, effectively ruling out a high-mass per-turbing planet. Therefore if this system contains a soli-tary shepherding planet, it must be low-mass.See Table 3 for a summary of the predicted masses andorbits for planets in each system.
Parent body disk widths
The predictions for planet mass and orbit in this studymake use of the dynamical effects on dust grains. Arethere any dynamical signatures on the parent bodies thatproduce the dust? Fig. 5 shows the final normalized par-ent body disk width as a function of planet mass ratio.Clearly there is more scatter such that degeneracies inmass exist for widths . .
10% can-not contain planets with mass ratios & µ/µ J . Whilenot as strong a constraint as can be made with scatteredlight images of a disk’s dust, this is still useful in rulingout very massive planets in narrow parent body disks. DISCUSSION: OBSERVER’S PROCEDURE
Observers can use the results of this study in threeways. (i)
For debris disks that have been resolved inscattered light but in which no planets have been de-tected, they can estimate the maximum mass, minimumsemimajor axis, and eccentricity of a putative solitaryshepherding planet. These can be calculated using thefollowing procedure:1. Calculate the eccentricity of the debris disk, e disk,f ,from the deprojected scattered light image. Theplanet’s eccentricity can be assigned this value.2. Construct the azimuthally-averaged radial profileof the disk’s deprojected scattered light, multiplythe profile by the distance from the star squaredto account for the geometric dilution of star light,and normalize the profile by the peak value.3. Calculate the semimajor axis of the peak in the de-projected azimuthally-averaged radial profile, andthe two locations equal to half the peak (a peak ,a in / , a out / ). Use Eq. 1 to calculate the disk’s nor-malized FWHM. mass ratio ( µ / µ J ) no r m a li z ed pa r en t bod y d i sk F W H M Fig. 5.—
Normalized parent body disk FWHM as a functionof planet mass ratio. The black points represent the averages ofthe values at each mass for the various disk eccentricities, and theerror bars represent the full range in nFWHM values for each mass.While a linear trend is evident, there is much more scatter thanin the plot relating the dust disk FWHM to the planet mass ratio(Fig. 4), making the predictive power less precise. However, wecan infer that parent body disk widths <
10% cannot contain aplanet with mass ratio & µ/µ J .
4. Insert the disk’s normalized FWHM into Eq. 5to solve for the maximum mass of the perturbingplanet.5. Assume a in / = a inner and insert this value, alongwith the calculated maximum mass of the planet,into Eq. 2 to solve for the planet’s minimum semi-major axis.This procedure (and for example, Fig. 1b) requiresthat debris rings be detected at high S/N in scatteredlight. This should be feasible in the coming yearssince HST/STIS is capable of detecting bright ring-likedisks at very high S/N (G. Schneider, private com-munication; C. Stark, private communication). Ad-ditionally, ground-based telescopes with adaptive op-tics (AO) systems are now detecting disks in scatteredlight at high S/N (Esposito et al. 2013; Thalmann et al.2011; Lagrange et al. 2012; Thalmann et al. 2013;Currie et al. 2012c; Boccaletti et al. 2012; Rodigas et al.2012; Buenzli et al. 2010). The visible light camera,VisAO (Kopon et al. 2010), on the Magellan AO system(MagAO; Close et al. (2010)) may be capable of produc-ing the highest resolution images yet on known brightdisks (Rodigas et al. 2014, in prep.), allowing for moreprecise measurements of scattered light disk widths. (ii) While no such system exists yet, for resolved de-bris disks that also contain a directly imaged planet, ob-servers can use the above procedure to place an atmo-spheric model-independent limit on the planet’s mass,semimajor axis, and eccentricity. This is especially use-ful because the masses of directly imaged planets de-pend heavily on atmospheric models (e.g., Baraffe et al.(2003); Burrows et al. (2003)), which themselves also de-pend on the (usually uncertain) age of the host star andon the initial conditions of the planet’s formation (ie,“hot-start” vs. “cold start”; Spiegel & Burrows (2012)). Rodigas, Malhotra, & Hinz (iii)
For debris disks with information only availablefor the parent bodies, rather than the dust, observers canuse Fig. 5 to constrain the maximum disk-shepherdingplanet mass. For example, models of an unresolved diskdetected by its infrared excess can sometimes estimatethe location and width of the parent body disk. This in-formation can be compared with our results to constraina putative shepherding planet’s mass. Additionally theAtacama Large Millimeter Array (ALMA) is now mak-ing it possible to directly image the tracers of parentbodies in debris disks (e.g., Boley et al. (2012)), whichcan provide model-independent measurements of parentbody disk widths for comparison with our results. SUMMARY
Using N-body simulations consisting of a star, a planet,and a disk of parent bodies and dust grains, we haveshown that the width of a debris disk in scattered lightis proportional to shepherding planet mass. This rela-tionship can be used to estimate the maximum planetmass in a debris disk of a given width, as well as theplanet’s eccentricity and minimum semimajor axis.Using the procedure outlined above, we have estimatedthe masses and orbits of putative shepherding planets infive bright debris rings. For Fomalhaut, deep imaging hasalready yielded mass limits below the maximum masseswe report here. For HR 4796A, our results indicate a massive planet might reside close to the star, evading pre-vious imaging detection efforts. For HD 207129, despitethe favorable separation between the star and putativeplanet, the low maximum planet mass and old age makethis system unfavorable for direct imaging. HD 202628can contain a very high-mass planet at a favorable pro-jected separation, making it attractive for direct imagingplanet searches. While HD 181327 can also potentiallyhost a high-mass perturber, Wahhaj et al. (2013) havealready ruled out such planets, implying that if the diskis being shepherded by a planet, the planet must be low-mass.In general, observers searching for planets should pri-oritize systems that contain wide debris disks, as thesecan tolerate more massive interior perturbers. Once aplanet is directly detected orbiting interior to its disk,its mass can be estimated independent of atmosphericmodels, providing a check on this fundamental physicalproperty.We thank the anonymous referee for very helpful com-ments which greatly improved this paper. We thank Eu-gene Chiang for helpful discussions, and for sharing datafrom Chiang et al. (2009). We thank Andy Skemer, JohnDebes, and Chris Stark for helpful discussions. T.J.R.was supported by the NASA Earth and Space ScienceFellowship (NESSF).
APPENDIX
ASSUMPTIONS, FORCES, AND TIMESCALES
Parent bodies
We make several simplifying assumptions regarding the parent bodies to facilitate the many simulations we carryout. First, we assume that all particles have zero mass (as in C09). Doing so allows us to treat each parent bodyas a test particle interacting purely gravitationally with a massive planet and central star (ie, the restricted 3-bodyproblem), which is much less computationally strenuous.The assumption that the parent bodies have zero mass is of course not realistic. The masses of debris disks are notwell constrained, but a few model-dependent measurements of the mass of Fomalhaut’s debris disk estimate a totalmass of ∼ ⊕ (C09; Boley et al. (2012); Acke et al. (2012)), where M ⊕ is the mass of Earth. If we conservativelyassume that our simulated parent body disks have masses at the upper end of this range, what happens to the planet?Bromley & Kenyon (2011) showed that for a single planet orbiting interior to a non-zero mass parent body disk,the planet will migrate away from the disk as long as the separation is & . J planet will have migrated away from a 100 M ⊕ parent body disk of 10% width by only ∼ a few percent of its originalsemimajor axis. More massive planets will migrate even less (Bromley & Kenyon 2011). Furthermore these resultsassume a constant migration rate over time, which is unlikely. Therefore neglecting migration and assuming zero massfor the parent bodies is reasonable for our purposes.Second, we assume (as in C09) that over the course of the numerical integrations, the parent bodies are not destroyedby collisions, and their orbits are unaffected by collisions with smaller bodies. The collisional lifetime of a particle ina debris disk, or time before a particle is catastrophically destroyed via collision with a similarly-sized particle, can bewritten as t collision = t per πτ ⊥ f , (A1)where t per is the orbital period of the particle in the disk, τ ⊥ is the disk’s vertical optical depth, and f is a factorthat depends on the particle size (Wyatt et al. 1999b). For a typical dust grain, f ∼ ∼ years (see the following subsection). From Wyatt et al. (1999b), thecollisional lifetime of a parent body at the top of the collisional cascade can be written as t collision,pb ≈ t collision (cid:18) D pb D dust (cid:19) / , (A2)where D pb and D dust are the diameters of the parent bodies and dust grains, respectively, and a Dohnanyi (1969) sizedistribution has been assumed. For a Fomalhaut-like debris disk, a typical dust grain size is ∼ µ m (C09). Usingthe above dust grain collisional lifetime, Eq. A2 shows that the collisional lifetime of parent bodies as small as just10-100 m is ∼ ∼ the age of Fomalhaut. Such small parent body sizes are rather conservative basedon models of extrasolar debris disks (C09, Wyatt et al. (1999b)), as well as estimates of primordial parent body sizesin the solar system’s asteroid and Kuiper belts (Morbidelli et al. 2009; Schlichting et al. 2013), which typically prefersizes of ∼ Dust grains
Because our study concerns debris disks (where little or no gas is present), the forces acting on dust grains areprimarily: gravity from the host star, gravity from the perturbing planet, radiation pressure from the host star, andPR drag. In addition, the lifetimes of dust grains are limited by grain-grain collisions, for these lead to catastrophicfragmentation. For a geometrically absorbing dust grain of bulk density ρ and radius s , the ratio, β , of the radial forceof stellar radiation pressure to the force of stellar gravity is, β = 3 L ∗ πGM ∗ cρs , (A3)where L ∗ and M ∗ are the stellar luminosity and mass, respectively, G is the universal constant of gravitation, and c isthe speed of light. In response to radiation pressure, dust grains generated by parent bodies moving on circular orbitsacquire larger, more eccentric orbits, with semimajor axis and eccentricity given by a = a pb (1 − β ) / (1 − β ) , e = β/ (1 − β ) . (A4)Only the particles with β < . s b = 3 L ∗ πGM ∗ cρs = 1 . ρ − ) − ( L ∗ L ⊙ )( M ∗ M ⊙ ) − µ m . (A5)Smaller particles, s < s b , acquire hyperbolic orbits upon release from the parent bodies, and their residence time inthe debris disk is t unbound ≈ w Ω( a pb ) = 16( w . a pb
100 AU ) / ( M ∗ M ⊙ ) − / yr , (A6)where Ω( a ) = p GM ∗ /a is the local Keplerian frequency and w is the normalized width of the disk. (For parentbodies on moderately eccentric orbits, the orbit of the dust grain depends on the longitude at which it is released, andso does the minimum blow out size, but the latter is still close to s b given above.) From Wyatt & Whipple (1950), forthe bound dust grains, PR drag causes orbits to circularize on a timescale t P R = (cid:12)(cid:12) ˙ ee (cid:12)(cid:12) − = πρsc a √ − e L ∗ ≈ × − β ) β (1 − β ) / ( M ∗ M ⊙ ) − ( a pb AU ) yr; (A7)the orbit decay timescale is | ˙ a/a | − ≈ − β ) t P R for β . . t collision ∼ (cid:18) τ R Ω( a ) w (cid:19) − ∼ × ( 10 − τ R )( w . / ( a pb
100 AU ) / ( M ∗ M ⊙ ) − / yr , (A8)where τ R is the radial optical depth, and τ R / √ w accounts for the effective optical depth for the path length a pb √ w traversed by a dust grain on an elliptical orbit of eccentricity near unity. In Eqs. A7 and A8 we have also used therelation between the dust grain’s orbital parameters and those of the parent bodies (Eq. A4).Comparing these three timescales for dust grains, we see that for typical debris disks imaged in scattered light (ie,those with large optical depths), t P R ≫ t collision ≫ t unbound (as in C09, Wyatt (2005)). This means that particlesbelow the blowout size do not contribute to the brightness of a debris disk. Therefore we only simulate bound particles.Additionally, we can see that on timescales on the order of the dust grain lifetimes, ∼ t collision , the bound dust grains’orbits do not change significantly from their initial orbits upon “release”. This justifies our neglect of PR drag in thesimulations. LIMITATIONS OF THIS STUDY
Here we describe several limitations of our study that we have not addressed elsewhere. We sought dynamicalstability of all parent bodies for only 10 orbits so that the wall clock time for a given simulation would be short.0 Rodigas, Malhotra, & HinzIntegration times of > orbits would have been more realistic since imaged debris disks can be as old as ∼ ∼ if the observed disk features are being created by a single planet orbiting interior tothe disk , how massive is that planet and what is its orbit like? It is certainly possible that other physical processes (e.g.,dust grain collisions, radiation forces, dynamical interactions with a distant star, the interstellar medium, multiplelow-mass planets (e.g., Raymond et al. (2011, 2012)), or even gas (Lyra & Kuchner 2013)) could produce disk featuresthat we assume are created by a solitary planet., how massive is that planet and what is its orbit like? It is certainly possible that other physical processes (e.g.,dust grain collisions, radiation forces, dynamical interactions with a distant star, the interstellar medium, multiplelow-mass planets (e.g., Raymond et al. (2011, 2012)), or even gas (Lyra & Kuchner 2013)) could produce disk featuresthat we assume are created by a solitary planet.