The Detectability of Exo-Earths and Super-Earths Via Resonant Signatures in Exozodiacal Clouds
aa r X i v : . [ a s t r o - ph ] O c t The Detectability of Exo-Earths and Super-Earths Via ResonantSignatures in Exozodiacal Clouds
Christopher C. Stark and Marc J. Kuchner ABSTRACT
Directly imaging extrasolar terrestrial planets necessarily means contendingwith the astrophysical noise of exozodiacal dust and the resonant structures cre-ated by these planets in exozodiacal clouds. Using a custom tailored hybridsymplectic integrator we have constructed 120 models of resonant structures cre-ated by exo-Earths and super-Earths on circular orbits interacting with colli-sionless steady-state dust clouds around a Sun-like star. Our models includeenough particles to overcome the limitations of previous simulations that wereoften dominated by a handful of long-lived particles, allowing us to quantitativelystudy the contrast of the resulting ring structures. We found that in the caseof a planet on a circular orbit, for a given star and dust source distribution, themorphology and contrast of the resonant structures depend on only two param-eters: planet mass and √ a p /β , where a p is the planet’s semi-major axis and β is the ratio of radiation pressure force to gravitational force on a grain. We con-structed multiple-grain-size models of 25,000 particles each and showed that in acollisionless cloud, a Dohnanyi crushing law yields a resonant ring whose opticaldepth is dominated by the largest grains in the distribution, not the smallest.We used these models to estimate the mass of the lowest-mass planet that canbe detected through observations of a resonant ring for a variety of assumptionsabout the dust cloud and the planet’s orbit. Our simulations suggest that planetswith mass as small as a few times Mar’s mass may produce detectable signaturesin debris disks at a p &
10 AU.
Subject headings: catalogs — circumstellar matter — infrared: stars — inter-planetary medium — methods: N-body simulations — planetary systems Department of Physics, University of Maryland, Box 197, 082 Regents Drive, College Park, MD 20742-4111, USA; [email protected] NASA Goddard Space Flight Center, Exoplanets and Stellar Astrophysics Laboratory, Code 667, Green-belt, MD 20771
1. Introduction
A number of proposed experiments like the Terrestrial Planet Finder (TPF) aim todirectly image the scattered and emitted light from extrasolar planets (Lawson & Traub2006). These experiments will also excel at detecting exozodiacal dust, circumstellar dustanalogous to zodiacal dust in our solar system (e.g. Agol 2007; Beckwith 2007). Zodiacaldust in the solar system consists of ∼ − µ m dust grains released through asteroidalcollisions and the outgassing of comets (e.g. Schramm et al. 1989). This dust forms thezodiacal cloud, extending from the solar corona (e.g. Mann et al. 2000) to beyond Jupiter(e.g. Kr¨uger et al. 1999).Our zodiacal cloud exhibits several structures interpreted as dynamical signatures ofplanets (Dermott et al. 1985, 1994; Reach et al. 1995). Several dusty disks around nearbymain-sequence stars show similar structures (e.g. Greaves et al. 1998; Wilner et al. 2002;Kalas et al. 2005). This trend suggests that exozodiacal clouds may be full of rings, clumpsand other asymmetries caused by planets and other phenomena.This situation raises some important questions. Will the structures in exozodiacal cloudsbe harmful astrophysical noise for direct imaging of extrasolar planets (Beichman 1996;Beichman et al. 1999)? Or can the dynamical signatures of planets in these clouds help usfind otherwise undetectable planets (e.g. Kuchner & Holman 2003)?Several studies have examined the geometry of resonant signatures of planets in debrisdisks (e.g. Kuchner & Holman 2003; Reche et al. 2008). However, most simulations cannotquantitatively study the contrast in these structures: how bright they are relative to thebackground cloud. We need to model the contrast of the structures in exozodiacal clouds tounderstand their roles as astrophysical noise and as signposts of hidden planets. However,accurately simulating the contrast of these structures demands computational resources thathave only recently become available (e.g. Deller & Maddison 2005).In this paper we examine the contrast of resonant structures induced by planets insteady-state exozodiacal clouds and the detectability of these structures via direct imaging.We simulate high-fidelity images of collisionless exozodiacal clouds containing a terrestrial-mass planet—an exo-Earth or super-Earth. By using roughly an order of magnitude moreparticles than most previous simulations, we overcome the Poisson noise associated withconstructing histograms of the column density and populating the external mean motionresonances (MMRs) of planets. We use our simulations to estimate the minimum planetmass that can be indirectly detected via observations of these structures as a function ofthe planet semi-major axis and dominant grain size under the assumption of circular planetorbits. Our models apply to exozodiacal clouds less than a few hundred times the optical 3 –depth of the solar zodiacal cloud, clouds for which the collision time is shorter than thePoynting-Robertson (PR) time for typical grains.Section 2 of this paper describes our numerical techniques. We present a syntheticcatalog of resonant debris disk structures in Section 3. We describe our multiple-particle-size cloud models and discuss their detectability in Section 4. In Section 5, we discuss thelimitations of our simulations; we summarize our conclusions in Section 6.
2. Numerical Method
Dust grains in the inner solar system are primarily released from parent bodies via colli-sions or outgassing. Radiation pressure ejects the smallest particles from the solar system ina dynamical time while the larger particles slowly spiral inward due to PR drag (Robertson1937; Burns et al. 1979). During their spiral toward the Sun, particles may become tem-porarily trapped in the MMRs of planets, extending their lifetimes by a factor of a few toten (Jackson & Zook 1989). This trapping locally enhances the particle density, creatingstructures within the zodiacal cloud, which have been described as circumsolar rings, bands,and clumps (e.g. Kelsall et al. 1998).To model these types of structures in exozodiacal clouds we numerically integrated theequation of motion of dust particles. The equation of motion for a perfectly absorbingparticle orbiting a star of mass m ⋆ is given to first order in v/c by Robertson (1937): d r dt = − Gm ⋆ r (1 − β )ˆ r − (1 + sw) βc Gm ⋆ r [ ˙ r ˆ r + v ] , (1)where r and v are the heliocentric position and velocity of the particle and sw is the ratioof solar wind drag to PR drag. We assume a value for sw of 0.35 (Gustafson 1994). Forperfectly absorbing spherical particles in the vicinity of the Sun, β ≈ . /ρs , where ρ is themass density of the particle in g cm − and s is the radius in µ m. We implemented a customized hybrid symplectic integrator to perform our numericalintegrations. Chambers (1999), hereby referred to as C99, introduced hybrid symplecticintegration as a method for dealing with close encounters in an efficient n-body code. Sym-plectic integrators rely on splitting the Hamiltonian into two easily integrable portions—adominant term, H D , and a smaller perturbative term, H P . However, in the n-body problem, 4 – H P may exceed H D during close encounters. Hybrid symplectic integrators overcome thisproblem by effectively switching from a symplectic integrator to an alternate integrator (e.g.Bulirsch-Stoer).The hybrid method reduces the perturbative term of the Hamiltonian, H P , by a factor K ( r ij ), where r ij is the distance between the two bodies in question, to ensure that theperturbative term remains relatively small. The integrator includes the remaining portion ofthe perturbative term, P i,j H P ,ij [1 − K ( r ij )], in the dominant term which is then integratedusing a method of choice. The “changeover function,” K ( r ij ), is a smooth function thatvaries from 0 for r ij . r crit to unity for r ij & r crit .Using a hybrid integrator requires choosing a changeover function and a value for r crit .We use the same changeover function as C99. We assign a different value of r crit ,i to eachbody, calculated as the larger of 3 R H ,i and τ v i , where R H ,i and v i are the Hill radius andvelocity of the i th body, respectively, and τ is the time step of the integrator. We thencalculate the critical distance for a pair of bodies as r crit ,ij = r crit ,i + r crit ,j .Our integrator also incorporates the effects of radiation pressure, PR drag, and solarwind drag. We implement radiation pressure as a correction to the effective stellar mass(cf. Eq. 1) and treat the drag effects as an additional term in H P , in much the same wayas Moro-Mart´ın & Malhotra (2002), hereby referred to as MMM02. We also use democraticheliocentric (DH) coordinates, composed of the barycentric momenta and heliocentric posi-tions, because of their relative ease of implementation. This choice introduces an additionalperturbative term to the Hamiltonian due to the motion of the star with respect to thebarycenter (Duncan et al. 1998). We checked our integrator using a variety of standard tests. We checked the energy andJacobi constant conservation with the drag terms turned off and examined the evolution ofdust particles’ orbital elements under our implementation of drag effects. We also comparedour hybrid integrator to a Bulirsch-Stoer integrator by examining the path of an individualtest particle during a close encounter and by examining the statistics of a cloud of particlesin a collisionless disk containing a planet.We tested energy conservation in our integration code by integrating the orbits of thefour outer planets and the Sun for 3 × years using a time step of 0.15 years. The energyerror was bounded with a mean value of ∆E / E ≈ × − . 5 –Duncan et al. (1998), which we will refer to as DLL98, tested the relative conservationof energy in their symplectic integrator as a function of planet perihelion distance. Wereplicated their tests using our code. Figure 1 shows the relative energy error in an integrationof the orbit of Jupiter for 3 × years using a time step of 0.15 years. We initially placedJupiter at aphelion. Figure 1 also shows the results of integrating the orbits of Jupiterand Saturn under the same conditions. With the DH method, the perturbative solar termincreases as the perihelion distance of the planet decreases, causing the fractional energyerror to increase similarly. The fractional energy errors shown in Figure 1 agree with thoseobtained by DLL98.We checked the conservation of the Jacobi constant by integrating particles in the Sun-Neptune system. We found results consistent with those of MMM02. Particles that did notundergo close encounters conserved the Jacobi constant at the level of ∼ − to 10 − .To test our implementation of PR drag, we replicated a test performed by MMM02.We integrated the orbit of a particle with β = 0 . .
35 in the presence of theSun. Figure 2 shows the semi-major axis and eccentricity as functions of time. These resultsmatch the results of MMM02 and agree with the analytic solution (Wyatt & Whipple 1950).We tested the performance of our hybrid scheme by integrating the orbit of cometP/Oterma in a close encounter with Jupiter, which has been done previously by Michel & Valsecchi(1996) and C99. The initial conditions for both bodies can be found in Table 3 of Michel & Valsecchi(1996). Figure 3 shows the path of comet P/Oterma for several values of integration timestep τ as seen in the frame co-rotating with Jupiter, which is located at the origin. Theseresults are similar to those obtained by C99. Our code shows a minor improvement over theother codes, most noticeable in the τ = 100 days case, that is likely only due to differencesin the calculation of the changeover distance. C99 explicitly sets r crit = 3 R H for this test;we used our prescription for r crit as described in Section 2.1. We directly compared simulations of resonant structures made with our hybrid integra-tor to simulations made with a Bulirsch-Stoer integrator. During the integrations we recordedthe coordinates of each particle in a 2-D histogram at regular intervals. This histogram mod-els the surface density distribution in a steady-state cloud. Since we only modeled planetson circular orbits, we simply recorded the coordinates in the frame co-rotating with theplanet. This technique has been widely used by dust cloud modelers (Dermott et al. 1994;Liou & Zook 1999; Moro-Mart´ın & Malhotra 2002; Wilner et al. 2002; Deller & Maddison 6 –2005).Figure 4 shows two histograms, one for each integrator, for simulations of 1,000 particleseach in the presence of the Sun/Earth system. We used a histogram bin width of 0 . β = 0 .
02 and initially released the particles with semi-majoraxis, a dust , distributed uniformly between 3 and 5 AU, eccentricity, e , uniformly distributedfrom 0 . .
1, and inclination, i , uniformly distributed between 0 ◦ and 6 ◦ . We used asymplectic time step of 0.02 years and recorded the particle locations every 250 years.Except for a small number of pixels, the middle panel of the figure (simulation using thehybrid integrator) looks qualitatively very similar to the left-most panel (simulation usingthe Bulirsch-Stoer integrator). The right panel of Figure 4 shows the difference of these twoimages divided by the √ n Poisson noise expected for each pixel where n is the number ofparticles in the pixel. This figure demonstrates that the differences between the two modelsare nearly consistent with the Poisson noise of the histograms. The two integrators resultedin histograms with minor structural differences, but the hybrid symplectic integrator runs afew times faster.Besides pixel-to-pixel Poisson noise in the histogram, this method is also sensitive tonoise in the population of MMRs. MMM02 showed that the population of the dominantresonances varied by a factor of ∼
3. Simulations & Results3.1. Cataloging Debris Disk Structure
To explore the range of different types of structures formed by terrestrial-mass planets,we performed 120 simulations of dust interacting with single planets on circular orbits. Thesimulations used 5,000 particles each and covered six values of planet mass, M p (0.1, 0.25,0.5, 1.0, 2.0, and 5.0 M ⊕ ), four values of planet semi-major axis, a p (1, 3, 6, and 10 AU), andfive values of β (0.0023, 0.0073, 0.023, 0.073, and 0.23) corresponding to spherical silicateparticles ranging in radius from ∼ − µ m. We released the particles on orbits withsemi-major axes uniformly distributed between 3.5 and 4.5 times the semi-major axis of theplanet’s orbit—well outside of the strongest MMRs. We used initial eccentricities uniformlydistributed between 0 and 0.2, initial inclinations uniformly distributed between 0 and 20 ◦ ,and the longitude of the ascending node, Ω, and the argument of pericenter, ω , uniformlydistributed between 0 and 2 π . We considered planet semi-major axes of 1 to 10 AU becausetypical designs for TPF can detect an exozodiacal cloud with 10 times the optical depth ofthe solar zodiacal cloud over roughly that range of circumstellar radii (Levine et al. 2006).We chose these initial conditions to model only dynamically-cold dust, i.e. e dust . . i dust . ◦ , since this component of a dust cloud is the dominant contributor to resonantring structure. We neglect dynamically hot dust with the idea that it can always be addedin later as a smooth background (Moran et al. 2004). The asteroid belt probably producesmuch of the solar system’s dynamically cold dust, while comets are thought to contribute amore dynamically hot cloud component (e.g. Liou et al. 1995; Ipatov et al. 2007). We treatonly steady-state dust clouds, assuming dust is continually replenished, and ignore transientcollisional events.Figure 6 shows some examples of the histograms from our simulations, which reveal awide range of trapping behavior. Some histograms show no azimuthal or radial structure,while others show high contrast rings. All of the patterns are Type I structures as identifiedby Kuchner & Holman (2003).Several general trends emerged. The ring contrast increased with increasing planetmass. Reducing β also enhanced trapping, as did increasing the planet’s semi-major axis.These last two trends can be explained by comparing the libration time of a given MMRto the PR time. The PR time scales as a /β , while the libration time for a given resonancescales as a / , where a dust is the semi-major axis of a dust grain’s orbit. The ratio of thesequantities yields √ a dust /β , a parameter that measures the degree to which resonant trappingis adiabatic (e.g. Henrard 1982); the trapping becomes more adiabatic and more efficient at 8 –greater distances from the star, and for larger particles. We discuss this phenomenon furtherin Section 3.3 below.In addition to following these trends, all of the simulated ring structures, like thoseshown in Figure 6, share some salient features:1. For cases in which even a modest amount of trapping occurs (azimuthally averagedcontrasts > ∼ . ≈ . a p . Thisfeature probably appears because the eccentricities of particles trapped in exteriorMMRs are typically pumped up to a limiting value before a close encounter with theplanet ejects them from resonance. For a particular MMR, all particles, regardless of β , tend to approach the same limiting eccentricity and accordingly, a similar pericen-ter distance (Beaug´e & Ferraz-Mello 1994). The limiting eccentricities are such thatthe limiting pericenter distances are nearly equal for the dominant resonances (e.g.the 2:1 and 3:2 resonances have limiting pericenter distances of 0 . a p and 0 . a p ,respectively), creating the ring structure’s sharp inner edge.2. A gap in the ring structure, a local minimum in the surface density, appears aroundthe planet. If we define gap width as the FWHM of the minimum in the azimuthalsurface density profile at r = a p , we find that the gap width is linearly proportionalto the contrast of the ring, as shown in the left panel of Figure 7. A linear fit to thedata shown in this figure gives w gap ≈ ◦ × C AA , IE for C AA , IE > .
6, where w gap is thegap width in degrees and C AA , IE is the azimuthally averaged inner-edge contrast (seeSection 3.2).3. The rings show a leading-trailing asymmetry. The trailing side of the ring structure isnoticeably denser than the leading side, and the structure is rotationally shifted in theprograde direction causing the trailing side to be closer to the planet (Dermott et al.1994). To examine the leading-trailing asymmetry caused by a prograde shift of thering structure, we measured the azimuthal offset of the center of the gap describedabove from the planet. The right panel of Figure 7 shows these measured progradeshifts of our simulations. Kuchner & Holman (2003) showed for a particular first orderexterior MMR, sin φ ∝ β (1 − β ) M p √ a p , (2)where φ is the prograde shift of the pericenter. Therefore, we plotted the sine of eachmeasured prograde shift against β (1 − β ) / ( M p √ a p ) in the right panel of Figure 7. Ourdata reveal the approximate proportionalitysin φ ring ∝ (cid:20) β (1 − β ) M p √ a p (cid:21) . , (3) 9 –where φ ring is the measured prograde shift of the ring structure. While the relationshipin Equation 2 holds for a single MMR, it does not strictly apply to a given ringstructure which consists of several well-populated MMRs. The relative populations ofthese MMRs are also functions of M p , a p , and β . However, this situation seems topreserve a power-law relationship between sin φ ring and β (1 − β ) / (cid:0) M p √ a p (cid:1) , as shownin Equation 3.4. The radial width of the ring increases with the contrast of the ring, ranging from a fewpercent of a p to ∼ . a p in the highest contrast case. As the trapping probabilities ofall the MMRs increase, MMRs farther from the planet’s orbit become populated. Forthis reason, the outer-edge of the ring structure differs significantly among simulations.The outer-edge can be quite blurry or very well-defined, making the radial width of aring structure difficult to quantify.Our catalog of debris disk structures induced by terrestrial-mass planets is publicly avail-able online at http://asd.gsfc.nasa.gov/Christopher.Stark/catalog.php. This online catalogalso contains images synthesized from the density distributions in scattered light and 10 µ mthermal emission assuming blackbody grains. Future studies of resonant ring structures withTPF or other experiments can use our catalog to interpret dust cloud patterns in terms ofplanet and dust parameters, assuming the observed image is dominated by a single grainsize. We envision the following process, inspired by recent papers on disks observed with theHubble Space Telescope (e.g. Clampin et al. 2003; Kalas et al. 2005):1. Deproject the image to remove inclination effects.2. Remove any smooth backgrounds by a power law fit.3. Estimate the dominant grain size in the resonant ring using infrared photometry orother methods.4. Compare the image of the disk to the online catalog to constrain the planet’s mass andlocation. We considered three different metrics for describing the ring contrast in our simulations: C Max : The surface density of the ring at its densest point divided by the surface densityof the background cloud 10 – C AA , Max : The maximum value of the azimuthally averaged surface density divided bythe surface density of the background cloud C AA , IE : The azimuthally averaged surface density at the inner edge of the ring dividedby the surface density of the background cloudWe calculated the above contrast metrics for all 120 simulations. We measured thesurface density of the background cloud at a circumstellar distance r ≈ . a p . The surfacedensity of the background cloud was nearly constant inside and outside of the ring, but didexhibit a small local minimum near r ≈ . a p in a few cases.To calculate C Max , we must search for the densest pixel, which introduces a bias towardpixels that exhibit an extreme amount of Poisson noise. To reduce this noise we averaged thesurface density over nine pixels centered on the densest point. Using C AA , Max or C AA , IE , onthe other hand, automatically averages over the effects of Poisson noise in our simulations.Figure 8 shows two examples of how the contrast, C AA , IE , depends on planet mass and β .Both plots show a similar behavior with three distinct regions: a no-trapping regime (contrast ∼ β , all contrasts converge to the same value for large planetmasses independent of planet semi-major axis, i.e. the contrast becomes “saturated” andincreasing the planet’s semi-major axis has little effect on the contrast. The right panel inFigure 8 illustrates this behavior; all four contrast curves, each of which corresponds to adifferent planet semi-major axis, approach the same value of ∼ M p = 5 M ⊕ . Similarly,for a given planet mass, contrasts converge to the same value for small β independent of a p ,as shown in the left panel in Figure 8. The morphology of the structure can vary, but thecontrast of the ring structure is roughly constant in these saturation regimes. As we mentioned above, dividing the PR time by the libration time of a given MMRyields a parameter, √ a p /β , that indicates the degree to which the resonant trapping isadiabatic. We plot the contrast in our simulations as a function of this parameter in Figure9. This figure demonstrates that for M p . M ⊕ and for a given distribution of parent bodyorbital elements, the ring contrast is a function of only two parameters: planet mass and √ a p /β .The morphology of the resonant rings is also, to good approximation, a function of only 11 –the planet mass and √ a p /β . The models shown in the two right panels of Figure 6 illustratethis phenomenon; for both models √ a p /β ≈
137 AU / . These two models have the samemorphology to a level consistent with pixel-to-pixel Poisson noise. Note that for small β inEquation 3, the prograde shift is approximately a function of (cid:0) √ a p /β (cid:1) − for a given planetmass.For large values of β and M p , the morphology and contrast of the ring structures are notsimple functions of √ a p /β . Simulations with large values of β , but equal values of √ a p /β (e.g. √ / .
073 and √
10 AU / .
23) show morphological differences, including differencesin prograde shift. Our simulations with M p = 5 M ⊕ also show contrast differences amongrings with equal values of √ a p /β .Wyatt (2003) investigated resonant trapping in MMRs for a system of planetesimalsexterior to an outward migrating planet on a circular orbit. Wyatt (2003) plotted thetrapping probability for a single MMR in his model as a function of migration rate and planetmass and found it could be well approximated by a function of the form P = [1+( ˙ a p /p ) p ] − ,where ˙ a p is the migration rate and the parameters p and p are power laws in planet mass.Our trapping scenario assumes dust migrating inward toward the planet, but the concept issimilar. Since contrast is closely related to trapping probability, we decided to fit the datashown in Figure 9 with a function of the form C = 1 + p (cid:18) (cid:18) p √ a p /β (cid:19) p (cid:19) − , (4)inspired by Wyatt (2003). Each of the three parameters, p i , is a power law in planet massof the form p i = p i, M p i, p . We fit all 120 contrast measurements with this six-parameterfunction for each of the three contrast metrics. The best fits, two of which are shown inFigure 9, are: C AA , IE : p ≈ .
38 ( M p /M ⊕ ) . , p ≈
207 ( M p /M ⊕ ) − . , p ≈ .
05 ( M p /M ⊕ ) . C AA , Max : p ≈ .
54 ( M p /M ⊕ ) . , p ≈
205 ( M p /M ⊕ ) − . , p ≈ .
63 ( M p /M ⊕ ) . C Max : p ≈ .
23 ( M p /M ⊕ ) . , p ≈
164 ( M p /M ⊕ ) − . , p ≈ .
72 ( M p /M ⊕ ) . Equation 4, combined with the above values, summarizes our results for all combinationsof planet mass, planet semi-major axis and β we simulated. Figure 9 shows the inner-edge contrast, C AA , IE , deviates significantly from the fits for large M p and large √ a p /β .The increased trapping efficiency for MMRs with these massive planets likely enhances thepopulation of MMRs farther from the planet’s orbit and depletes the inner MMRs that causethe sharp inner edge. 12 –
4. Multi-Particle-Size Models4.1. Composite Simulations
We used our 120 simulations to produce 20 multiple-particle-size dust cloud modelsby forming weighted sums of the histograms assuming a Dohnanyi distribution of particlesizes (Dohnanyi 1969). Each of these composite models effectively utilizes 25,000 particles.Exactly how we apply the ideas in Dohnanyi (1969) has profound effects on our compositemodels, so we present here two different kinds of models.First, we assembled a composite model in which the particles are initially released fromtheir parent bodies according to a crushing law, and do not undergo any further collisionalprocessing as they spiral inward. This scenario models a sparse disk with a belt of dust-producing material, like our own zodiacal cloud. The crushing law for asteroid material atmicron sizes is unknown, so we choose the crushing law used by Dohnanyi (1969): dNds ∝ s − α , (5)where dN is the number of particles with radius s in a bin of width ds , and α = 3 . τ , for our composite models from τ = P i w i A i σ i , where w i , A i , and σ i are the weighting factor, particle cross-section, and surface number densityof the i th single-particle simulation, respectively. We assume that the cross-section of eachparticle is A i ∝ β − . The crushing law in Equation 5 implies a weighting factor for the i th histogram of w i = β α − i ∆ β i , where ∆ β i is the width of the β i bin. For a constant logarithmicspacing in β , like the spacing we used in our simulations, and the Dohnanyi crushing law, w i = β . i .Larger particles have longer PR times, so in the absence of collisional processing, theirdensity is enhanced by a factor of β − under our assumption of a steady-state cloud model.One might expect that this effect must be included in the weighting factor. However, oursimulations include this effect automatically as long as we keep the frequency with whichparticle locations are recorded constant among all of our simulations. We did, in fact, varythe recording frequency with the PR time, but we corrected for the differences in recordingfrequency before summing the histograms.Figure 10 shows the optical depth of one of our 20 composite models ( M p = 2 . M ⊕ , a p = 6 . t PR ∝ β − ) and more likely to be trapped in MMRs. Hence, the upperleft panel of the figure closely resembles the lower right panel.Next, for the purpose of illustration, we ignored the initial size distribution of dustparticles and forced the disk to obey a size distribution of dNds ∝ s − . (6)at a radius of ∼ a p from the star. This scenario probably doesn’t have a physical inter-pretation, but it illustrates an interesting phenomenon: how resonant trapping tends to sortparticles by size. We enforce the size distribution at one location within the disk, but thesize distribution will not follow a Dohnanyi distribution elsewhere in the disk.The top right panel in Figure 10 shows the optical depth of an example of this kind ofcomposite cloud, normalized to a Dohnanyi distribution at ∼
18 AU. Models constructed inthis fashion are more greatly affected by the smallest grains. Hence, the top right panel doesnot greatly resemble the lower right panel in Figure 10.
We can further develop these ideas with a simple semi-analytic treatment. For a givenplanet mass and semi-major axis, the contrast function (see Equation 4) becomes C ( s ),where s is the particle size. We approximate the contrast function in Equation 4 with thepiecewise function C ( s ) = s < s (cid:16) ss (cid:17) m for s < s < s C large for s > s , (7)where C large = 1 + p is the contrast for the largest particles, m is the logarithmic slopeof the contrast in the transition regime, and s and s are the particle sizes that mark thebeginning and end of the transition regime, respectively. We fit our contrast data with thispiecewise function and obtained the following power law estimates assuming silicate grains( ρ ∼ − ): C large;AA , IE ≈ . (cid:18) M p M ⊕ (cid:19) . (8) 14 – m ≈ . (cid:18) M p M ⊕ (cid:19) . (9) (cid:18) s µ m (cid:19) ≈ (cid:18) M p M ⊕ (cid:19) − . (cid:16) a p (cid:17) − . (10) (cid:18) s µ m (cid:19) ≈ (cid:18) M p M ⊕ (cid:19) − . (cid:16) a p (cid:17) − . (11)In the same manner as Section 3.2, we defined the contrast of any ring structure as thesurface density within the ring, σ ring , divided by the background surface density, σ BG . Thecontrast in optical depth of a cloud containing several components of various sized particles,labeled with the index i , is h C τ i = X i C i σ BG ,i A i X i σ BG ,i A i , (12)For a collisionless cloud with a continuous distribution of grain sizes, the contrast in opticaldepth of the composite cloud is given by h C τ i = Z s max s min s − α C ( s ) ds Z s max s min s − α ds , (13)where we have explicitly included the particle cross section ( A ( s ) ∝ s ) and backgroundsurface density ( σ BG ( s ) ∝ s − α ). For a collisionless cloud, the background surface density isenhanced by a factor of s due to the PR time scaling as s (see Section 4.1), and a factor of s − α , which describes the assumed crushing law.Using Equations 7, we can now integrate Equation 13 directly. Assuming s min < s α = 4 + m C large (cid:16) s s min (cid:17) − α − (cid:16) s s min (cid:17) − α (cid:16) (cid:16) s s (cid:17)(cid:17) for α > α = 4 + m. (14) 15 –For cases in which the maximum particle size in a disk is less than s (see Equation 11),simply replace all instances of s in Equations 14 with s max . Similarly, for cases in which theminimum particle size in a disk is greater than s (see Equation 10), replace all instances of s with s min .Equations 14, together with Equations 8–11, give analytic expressions for optical depthcontrast in terms of M p , a p , s min , and s max . Although Equations 14 address all possiblescenarios, the most plausible scenarios have crushing laws with α < h C τ ;AA , IE i ≈ . M ′ . + s ′ α − X h(cid:0) M ′− . a ′− . (cid:1) − α − (cid:0) . M ′ . (cid:1) (cid:0) M ′− . a − . (cid:1) − α i for s max > s Xs ′ α − (cid:0) M ′− . a ′− . (cid:1) − α + (1 − X ) s ′ . M ′ . max (cid:0) M ′− . a ′− . (cid:1) − . M . , for s max < s (15)where M ′ p = (cid:16) M p M ⊕ (cid:17) , a ′ p = (cid:0) a p (cid:1) , s ′ max = (cid:16) s max µ m (cid:17) , and X = . M ′ . − α +0 . M ′ . .If, as in our first composite model in Section 4.1, we assume that the particles arereleased from their parent bodies in accordance with the Dohnanyi (1969) crushing law andthen spiral inward without colliding, α = 3 .
4. In this case, Equations 14 give a contrast inoptical depth of h C τ i ≈ C large in the limit s max ≫ s . This result confirms our numericalresults for our first composite model, shown in the upper left panel of Figure 10; the contrastin optical depth is dominated by the large particles.For our second composite model, we forced the background density to obey a Dohnanyidistribution at ∼ a p , i.e. σ BG ( s ) ∝ s − . , so that α = 4 .
5. This technique essentiallyremoves the factor of s in the background surface density that results from the PR timescaling as s . For the composite cloud shown in Figure 10, M p = 2 . M ⊕ , and a p = 6 . m ≈ . s ≈ . µ m, s ≈ µ m, and C large , AA , IE ≈
6. We let eachsimulated particle size represent a range of particle sizes using the midpoint method, whichgives s min ≈ . µ m. With these values, Equations 14 give a contrast in optical depth of h C τ, AA , IE i ≈ .
4, in agreement with the measured contrast in the top right panel of Figure10. Figure 11 illustrates in general how a distribution of particle sizes affects the contrast ofa ring structure. This figure compares the contrast of a collisionless multi-particle-size cloud(Equations 14) to that of a single-particle-size cloud as a function of √ a p /β min assuming a 16 –Dohnanyi (1969) crushing law. Both kinds of clouds have the same contrast in the adiabaticlimit (large √ a p /β min ), but the contribution of the smaller grains reduces the contrast else-where, effectively broadening the transition between the no-trapping regime and saturationregime. Crushing laws with α < . The detectability of a resonant ring structure depends on many factors specific to thetelescope being used and the observing conditions. We address this complicated issue byimposing one simplifying assumption: a minimum detectable optical depth ring contrast of1.5. This assumption likely underestimates the sensitivity of a TPF-like mission to rings inexozodiacal clouds analogous to the solar zodiacal cloud. In such a cloud, a ring 0.4 AUwide located at 1 AU from the star has ∼
15 times the total flux of an Earth-like planet at 1AU, even for a contrast of unity. Our assumption, conservative on the basis of photon noisealone, allows for the possibility of unknown systematic noise that could hinder the detectionof extended structures.Figure 12 shows the minimum detectable planet mass as a function of semi-major axisand maximum dust particle size based on Equations 15 and a Dohnanyi (1969) crushing law.The masses and semi-major axes of Earth, Mars, and the planet OGLE-2005-BLG-390Lb,detected by the microlensing technique (Beaulieu et al. 2006), are marked for reference. Thisplot shows that an Earth-mass planet at 1 AU might be detectable if the ring contains grainsmore than a few tens of microns in size and a planet with mass equal to a few times thatof Mars might be detectable near 10 AU if the ring contains grains more than one hundred 17 –microns in size.The detectability of a ring structure depends upon the size distribution of dust withinthe ring structure. Dust produced according to a crushing law less steep than the Dohnanyi(1969) crushing law ( α < .
4) will result in more highly contrasted ring structures becauseof the increased relative contribution of the large grains. For crushing laws with α = 3 and α = 2, the curves of constant maximum particle size shown in Figure 12 shift downward bya factor of approximately 1.25 and 1.55, respectively.These values are subject to the assumptions of our simulations, which do not include adynamically hot component in the dust cloud. This component would reduce the contrastin the ring, making planets harder to detect for a given cloud mass. So the detection limitsshown in Figure 12 should be thought of as best-case scenarios.
5. Caveats
Our simulations include a number of simplifying assumptions, which we summarize here.We ignored the effects of dynamically hot dust, like dust that might come from comets.Trapping probability decreases dramatically for particles on highly eccentric and inclinedorbits, so we expect dynamically cold dust to dominate any resonant debris disk structure.As a first approximation, we can treat the contribution from the dynamically hot dust asa constant surface density cloud component, which reduces the contrast of any structureformed from the dynamically cold component. Estimates of the ratio of asteroidal dust tocometary dust in our solar system range from 1:10 to 7:10 (Ipatov et al. 2007). For othersystems, this ratio is also unknown.Our simulations also assumed a single planet on a circular orbit around a Sun-like star.We have performed trial simulations of our solar system and demonstrated that the presenceof Jupiter may reduce the Earth’s ring contrast. Other multiple-planet systems may alsoexhibit a similar effect. Additionally, planets on eccentric orbits give rise to additional MMRswith different capture probabilities and geometries (Kuchner & Holman 2003).The ring contrasts of inclined systems can vary significantly depending on the inclinationand radial extent of the dust cloud. In edge-on systems, resonant features can overlap asseen from the Earth, complicating their interpretation. The contrasts we provide are usefulonly to systems for which projection effects can be taken into account.Finally, our multi-particle-size models demonstrate the subtlety of collisional effects indust clouds. Collisional effects can determine the relative populations of large and small 18 –grains and potentially alter the morphology of the ring structures. Our simulations can notyet handle these effects in detail.
6. Conclusions
We have implemented our own hybrid symplectic integrator for the n -body problem andused it to simulate collisionless debris disks, taking into account solar wind and drag effects.Each simulation contained 5,000 particles. We found that this number of particles sufficesto populate the dominant MMRs of a low-mass planet with an accuracy at the few percentlevel, yielding for the first time models of the surface brightness distributions of exozodiacalclouds that we can use to quantitatively study the contrasts of resonant features—not justtheir geometries.We generated a catalog of resonant structures induced by a single planet on a circular or-bit around a Sun-like star, available online at http://asd.gsfc.nasa.gov/Christopher.Stark/catalog.php.We investigated 120 sets of model parameters, spanning a range of planet masses, planetsemi-major axes, and values for β , assuming dust grains launched from orbits with low e dust and i dust . The resulting ring structures exhibited leading-trailing asymmetries, gaps near thelocations of the planets, and sharp inner edges at ≈ . a p .We performed a detailed analysis of the surface density contrasts of the rings (Figure9). We showed that for a planet on a circular orbit, the contrast and morphology of therings are to good approximation functions of only two parameters, M p and √ a p /β , for agiven stellar mass and distribution of dust sources for M p . M ⊕ and β . .
25. Equation 4summarizes the contrasts of our single-particle-size models as a function of these parameters.Considering only the dynamically cold particles analogous to particles released by asteroidsin the solar system, we find that terrestrial-mass planets are capable of producing resonantring structures with azimuthally averaged contrasts up to ∼ β values, we assembled multi-particle-size models of 25,000 particles each. Releasing the particles according to a Dohnanyi(1969) crushing law without any subsequent collisional processing results in composite cloudswhose optical depths are dominated by large particles; large particles will dominate imagesof these clouds in visible light and throughout the IR. Based on these composite models,we suggested that the best current models for exozodiacal clouds are those with a narrowrange of grain sizes corresponding to grains whose collision time roughly equals their PRtime. Future models should account for processes like grain-grain collisions that destroylarge grains. 19 –Equations 14 and 15 provide semi-analytic predictions for the contrast in optical depthof a multi-particle-size cloud of dynamically cold grains. For ring structures composed ofsilicate grains released according to a Dohnanyi crushing law ( α = 3 . h C τ ;AA , IE i ≈ . (cid:16) M p M ⊕ (cid:17) . − (cid:16) s max µ m (cid:17) − . (cid:0) a p (cid:1) − . (cid:16) M p M ⊕ (cid:17) − . for s max > s “ M p M ⊕ ” . (cid:20) (cid:16) s max µ m (cid:17) − . (cid:16) M p M ⊕ (cid:17) − . (cid:0) a p (cid:1) − . + (cid:18) (cid:16) s max µ m (cid:17) − . (cid:16) M p M ⊕ (cid:17) − . (cid:0) a p (cid:1) − . (cid:19) − “ M p M ⊕ ” . for s max < s (16)where h C τ ;AA , IE i is the ratio of the azimuthally averaged optical depth in the ring structureto the azimuthally averaged background optical depth, s max is the maximum grain size inthe ring structure, and s is given by Equation 11. For the case s max > s , the first twoterms in Equation 16 represent the contrast in the adiabatic limit. The remaining term (andthe terms in the s max < s case) represents deviations from this limit for smaller particles orsmaller semi-major axes.We plotted the mass of the smallest planet that could be detected through observationof a resonant ring structure as a function of planet semi-major axis and particle size inFigure 12. We assumed a cloud composed of a range of particle sizes adhering to a Dohnanyi(1969) crushing law and a minimum detectable optical depth contrast of 1.5:1. We foundthat planets with masses just a fraction of the Earth’s may form detectable ring structuresif the rings harbor grains more than several tens of microns in size.We would like to thank the National Aeronautics and Space Administration and God-dard Space Flight Center for support of this research through funding from the GraduateStudent Researchers Program, and the Harvard-Smithsonian Center for Astrophysics andNASA’s Navigator Program for their financial support via the Keck Interferometer NullerShared Risk Science Program. REFERENCES
Agol, E. 2007, MNRAS, 374, 1271Beaulieu, J.-P., et al. 2006, Nature, 439, 437 20 –Beaug´e, C., & Ferraz-Mello, S. 1994, Icarus, 110, 239Beckwith, S. V. W. 2007, preprint (arXiv:0710.1444v2)Beichman, C. A. 1996, A Road Map for the Exploration of Neighboring Planetary Systems(JPL Publication 96-22)Beichman, C. A., Woolf, N. J., & Lindensmith, C. A. 1999, The Terrestrial Planet Finder:Origins of Stars, Planets, and Life (JPL Publication 99-3)Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1Chambers, J. E. 1999, MNRAS, 304, 793Clampin, M., et al. 2003, AJ, 126, 385Deller, A. T., & Maddison, S. T. 2005, ApJ, 625, 398Dermott, S. F., Jayaraman, S., Xu, Y. L., Gustafson, B. ˚A. S., & Liou, J. C. 1994, Nature,369, 719Dermott, S. F., Nicholson, P. D., Burns, J. A., & Houck, J. R. 1985, in Properties andinteractions of interplanetary dust, ed. R. H. Giese & P. Lamy (Dordrecht: D. ReidelPublishing Co.), 395Dohnanyi, J. S. 1969, J. Geophysical Research, 74, 2531Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067Durda, D. D., Flynn, G. J., Sandel, L. E., & Strait, M. M. 2007, in ESA SP-643, Workshopon Dust in Planetary Systems, ed. H. Krueger & A. Graps, 77Fixsen, D. J., & Dwek, E. 2002, ApJ, 578, 1009Greaves, J. S., et al. 1998, ApJ, 506, 133Gustafson, B. ˚A. S. 1994, Annu Rev. Earth Planet. Sci., 22, 553Henrard, J. 1982, Celest. Mech., 27, 3Ipatov, S. I., Kutyrev, A. S., Madsen, G. J., Mather, J. C., Moseley, S. H., & Reynolds, R.J. 2007, Icarus, in pressJackson, A. A., & Zook, H. A. 1989, Nature, 337, 629 21 –Kalas, P., Graham, J. R., & Clampin, M. 2005, Nature, 435, 1067Kelsall, T., et al. 1998, ApJ, 508, 44Kr¨uger, H., et al. 1999, Planet. Space Sci., 47, 363Kuchner, M. J., & Holman, M. J. 2003, ApJ, 588, 1110Lawson, P. R., & Traub, W. A. 2006, NASA Navigator Science PlanLevine, M., Shaklan, S., & Kasting, J. 2006, Terrestrial Planet Finder Coronagraph: Scienceand Technology Definition Team (STDT) Report (JPL Publication D-34923)Liou, J. C., Dermott, S. F., & Xu, Y. L. 1995, Planet. Space Sci., 43, 717Liou, J. C., & Zook, H. A. 1999, AJ, 118, 580Mann, I., Krivov, A., & Kimura, H. 2000, Icarus, 146, 568Michel, P., & Valsecchi, G. B. 1996, CeMDA, 65, 355Moran, S. M., Kuchner, M. J., & Holman, M. J. 2004, ApJ, 612, 1163Moro-Mart´ın, A., & Malhotra, R. 2002, AJ, 124, 2305Reach, W. T., et al. 1995, Nature, 374, 521Reche, R., Beust, H., Augereau, J.-C., & Absil, O. 2008, A&A, in pressRobertson, H. P. 1937, MNRAS, 97, 423Schramm, L. S., Brownlee, D. E., & Wheelock, M. M. 1989, Meteoritics, 24, 99Wilner, D. J., Holman, M. J., Kuchner, M. J., & Ho, P. T. P. 2002, ApJ, 569, 115Wyatt, M. C. 2003, ApJ, 598, 1321Wyatt, M. C. 2005, A&A, 433, 1007Wyatt, S. P., & Whipple, F. L. 1950, ApJ, 111, 134
This preprint was prepared with the AAS L A TEX macros v5.2.
22 – −8 −7 −6 −5 −4 −3 −2 −1 ∆ E / E Fig. 1.— Maximum fractional error in energy during a 3,000-year integration as a function ofperihelion distance for two scenarios: a two-body system of the Sun & Jupiter (solid line) anda three-body system of the Sun, Jupiter & Saturn (dashed line). For the two-body system,Jupiter’s perihelion distance is plotted. For the three-body system, Saturn’s periheliondistance was altered while Jupiter’s remained fixed. The inclinations and eccentricities ofboth planets remained fixed. cf. DLL98 Fig. 3. 23 – e × × × × Time (yrs)050100150200 a ( AU ) Fig. 2.—
Top: eccentricity as a function of time for a dust particle with β = 0 . . Bottom: semi-major axis as a function of time. Our results match those ofMMM02 (cf. MMM02, Fig. 1) and agree with the analytical solution (Wyatt & Whipple1950). 24 – −2 −1 0 1 2x (AU)−2−1012 y ( AU ) Bulirsch−StoerHybrid 1 dayHybrid 25 daysHybrid 50 daysHybrid 100 days
Fig. 3.— The integrated trajectory of comet P/Oterma during a close encounter with Jupiteras viewed in the frame centered on and rotating with Jupiter with the Sun on the negative x -axis. Shown are the results of a Bulirsch-Stoer integrator and our hybrid symplecticintegrator for four values of integration time step. The hybrid symplectic results overlap theBulirsch-Stoer results for a timestep of 1 day (cf. C99, Fig. 4). 25 – Bulirsch−Stoer Hybrid Symplectic Difference
Surface Density (Linear Scale) Surface Density Surface Density (Linear Scale) Surface Density σ −7.2 −3.2 0.8 4.8 8.8 σ Fig. 4.— Comparison of our hybrid symplectic integrator with a Bulirsch-Stoer integrator.
Left: surface density histogram for 1,000 particles in the Sun-Earth system using a Bulirsch-Stoer integrator.
Middle: surface density histrogram for the same initial conditions using ourhybrid symplectic integrator.
Right:
Bulirsch-Stoer histogram minus the hybrid symplectichistogram (image is in units of σ , the √ n Poisson noise associated with the histograms).Except in a handful of pixels, the difference is roughly consistent with Poisson noise. 26 –
20 25 30 35 40 45 50 55a (AU)02 × × × × × N u m b e r Fig. 5.— Population of Neptune’s MMRs for three independent simulations of 5,000 particleseach (shown in green, red, and black). Populations of the 2:1 and 3:2 resonances differ amongthe three simulations by 6.4% and 4.3%, respectively (c.f. Moro-Mart´ın & Malhotra 2002,Fig. 5). 27 – M p =0.1 M ⊕ , a p =3 AU, β =0.073M p =5 M ⊕ , a p =10 AU, β =0.0023 M p =0.5 M ⊕ , a p =10 AU, β =0.023M p =0.5 M ⊕ , a p =1 AU, β =0.0073 Surface density (Linear Scale) Surface density Fig. 6.— Surface density distributions for four of the 120 simulations (scale is relative).The star is located at the center of the image and the planet is marked with a white dot.The planet orbits counter-clockwise in these images. Integrations were truncated at half theplanet’s semi-major axis. The simulations shown on the right have different values of a p and β , but the same value of √ a p /β . Their surface density distributions are nearly identical;their difference is consistent with Poisson noise (see Section 3.2). 28 – AA,IE )020406080 G a p w i d t h ( F W H M i n d e g r ee s ) β (1− β ) / [(M p /1 M ⊕ )(a p /1 AU) ]0.010.101.00 s i n ( φ r i ng ) Fig. 7.—
Left: the angular size of the “gap” in the ring structure around the location of theplanet in our simulations versus the contrast of the ring structure.
Right: the sine of theprograde shift plotted against the function β (1 − β ) / ( M p a / ). We removed all data with C AA , IE < . β C on t r a s t ( C AA , I E ) a p = 10.0 AUa p = 6.0 AUa p = 3.0 AUa p = 1.0 AU M p = 2.0 M ⊕ p (M ⊕ )110 C on t r a s t ( C AA , I E ) a p = 10.0 AUa p = 6.0 AUa p = 3.0 AUa p = 1.0 AU β = 0.0073 Fig. 8.— The azimuthally-averaged contrast measured at the inner edge of the ring structure(see Section 3.2 for definition of contrast) as a function of β (left panel) and planet mass(right panel). Both figures show a transition from a no-trapping regime to a saturationregime where contrast is independent of semi-major axis. 30 –
10 100 1000(a p /1 AU) / β C on t r a s t ( C AA , I E ) M p = 5.0 M ⊕ M p = 2.0M p = 1.0M p = 0.5M p = 0.25M p = 0.1
10 100 1000(a p /1 AU) / β C on t r a s t ( C M a x ) Fig. 9.— The contrast in surface density of the ring structure compared to the backgroundcloud for all combinations of M p , a p , and β (see Section 3.2 for definitions of contrast). Thecontrast is only a function of two parameters: planet mass and √ a p /β . The solid lines arefits to the data (see Equation 4). 31 – Dohnanyi distribution 1 Dohnanyi distribution 2 β = 0.23 β = 0.0023 Optical depth (Linear Scale) Optical depth Fig. 10.— Comparison of the optical depths for a composite cloud formed by two differentmethods for M p = 2 . M ⊕ and a p = 6 . Top-left:
A composite collisionless cloud where the parti-cles are released with a size distribution equal to the crushing law used by Dohnanyi (1969).
Top-right:
The same composite cloud, but formed by forcing the surface density outside ofthe ring structure to obey a Dohnanyi distribution.
Bottom-left:
The optical depth of thesmallest particles included in the composite clouds.
Bottom-right:
The optical depth of thelargest particles included in the composite clouds. The largest particles dominate the opticaldepth in a cloud of particles released with a Dohnanyi crushing law. 32 – (a p /1 AU) / β min C on t r a s t ( C τ , AA , I E ) multi−particle−size cloudsingle−particle−size cloud Fig. 11.— Contrast in optical depth for multi-particle-size clouds (solid lines) compared tosingle-particle-size clouds (dashed lines) assuming a Dohnanyi (1969) crushing law ( α = 3 . M ⊕ . The contributions of the smallgrains reduce the contrasts of the multi-particle-size clouds compared to single-particle-sizeclouds with the same minimum value of β . 33 – p (AU)0.11.010.0 M i n . D e t ec t a b l e P l a n e t M a ss ( M ⊕ ) s m a x = µ m s m a x = µ m s m a x = µ m s m a x = µ m s m a x = µ m s m a x = µ m E M O
Fig. 12.— Minimum detectable planet mass in a multi-particle-size collisionless cloud as afunction of semi-major axis and maximum grain size, assuming a Sun-like star, a minimumdetectable ring contrast of C τ, AA , IE = 1 .
5, and dust produced according to a Dohnanyi(1969) crushing law ( α = 3 . E and M , respectively. The 5 . ⊕ exoplanet OGLE-2005-BLG-390Lb is denoted with an O (Beaulieu et al. 2006). Listed values for maximum dust size in the ring structure assumeperfectly absorbing spherical grains with mean density ρ = 2 . − and radius s max .The bold line shows the case of the solar zodiacal cloud, for which the observed emission isdominated by 30 µµ