Table Of ContentMon.Not.R.Astron.Soc.000,000–000 (0000) Printed11January2013 (MNLATEXstylefilev2.2)
Effects of a planetesimal debris disk on stability scenarios
for the extrasolar planetary system HR 8799
Alexander Moore & Alice C. Quillen
3
1 Department of Physics and Astronomy, Universityof Rochester, Rochester, NY 14627, USA
0
2
n
a 11January2013
J
9
ABSTRACT
] HR 8799 is a four planet system that also hosts a debris disk. By numerically inte-
P
gratingbothplanetsandaplanetesimaldisk,wefindinteractionsbetweenanexterior
E
planetesimal disk and the planets can influence the lifetime of the system. We first
.
h consider resonant planetary configurations that remained stable for at least 7 Myrs
p sans debris disk. An exterior debris disk with only ∼ 1% the mass of the outermost
- planet (approximately a Neptune mass) was sufficiently large enough to pull the sys-
o
tem out of resonance after 2 to 6 Myrs. Secondly, we consider configurations which
r
t are unstable in less than a few hundred thousand years. We find that these can be
as stabilized by a debris disk with a mass of more than ∼ 10% that of the outermost
[ planet. Our two sets of simulations suggest that estimates of the long term stability
of a planetary system should take into account the role of the debris disk.
1
v 1 INTRODUCTION reportedvaluesforHR8799b andaddnewobservationsfor
4 HR8799c and d (Soummer et al. 2011).
0 PhotometricsurveysconductedwithKeplerspacetelescope
Dynamical studies of HR 8799 indicate that it is
0 (Borucki & Koch 2011) in addition to radial velocity sur-
likely in a 4:2:1 dual mean motion resonance (MMR).
2 veys such as those conducted by Wright et al. (2009) have
1. indicated that multiple planet systems are common. Nu- In other words, the inner two planets are in a 2:1
mean motion resonance while the outer pair are also in
0 merical integrations can be used to determine if these sys-
a 2:1 mean motion resonance. This architecture is re-
3 tems are stable and estimate the time to first orbit cross-
quired to explain how the system has remained stable
1 ing event or collision (Go´zdziewski & Migaszewski 2009;
over its observed age (Go´zdziewski & Migaszewski 2009;
: Reidemeister et al. 2009; Fabrycky& Murray-Clay 2010;
v Fabrycky& Murray-Clay 2010; Reidemeister et al. 2009;
Xi Mthaartshhaalsl eatsahl.o2rt01l0if;eMtimareoiissectoanls.id20er1e0d).lAesnsylcikoenlfiy.guCroantisoen- Marshall et al. 2010; Marois et al. 2010). Possible orbital
configurations determined by numerical simulations are
r quently, integrations can be used to place constraints on
a summarized by Moro-Mart´ın et al. (2010).
both the orbital elements and masses of the planets. These
HR8799alsohasaninnerdebrisdiskrangingfrom6-15
numerical investigations often neglect planetesimals. How-
ever, extrasolar planetary systems can harbor planetesimal AU, an outer debris disk which is thought to extend from
90to300AUandadustyhalooutto∼1000 AU(Su et al.
debris disks(Su et al. 2009; Moro-Mart´ın et al. 2010).
HR8799isa1.5±0.3M A5Vstarfound39.4±1.0pc 2009). The specific details of this model will be discussed
sol furtherin subsection 1.2.
from Earth (Marois et al. 2008). It has at least four plan-
ets which have been directly imaged (Marois et al. 2010). Planetesimal debris disks can influence the long term
Mass estimates for the planets, even taking into account stability of a planetary system. Within the context of the
HR 8799’s young age of 60+100 Myr determined by var- ‘Nice’ model (Tsiganis et al. 2005), planets migrate due
−30
ious techniques (Marois et al. 2008), or 30+20 as part of to interactions with planetesimals, and instability occurs
−10
the Columba association (Marois et al. 2010), have values when two planets cross a strong mean motion resonance.
of 10±3 M for HR 8799c, d, and e and 7+4 M for Thommes et al.(2008)consideredsystemsputinresonance
Jupiter −2 J
HR8799bforanassumptionofa60Myrage,andmassesof byagasdisk,apossiblescenario explainingHR8799’s cur-
7+3 M for HR 8799c, d, e and 5M for an assumption of rent configuration. However, after the gas disk had been
−2 J J
a 30 Myr age. Projected separations for the planets of HR depleted, they found that planetary interactions with the
8799e, d, c and b from thestar are observed to be14.5, 24, remnant planetesimal disk tended to remove these systems
38 and 68 AU, respectively (Marois et al. 2010). However, from resonances and induce dynamical instability.
thismeasurementlacksalongbaselinehelpfulforconstrain- Itisin thiscontextthatweexaminetheroleof thede-
ingtheplanets’positions.Astrometricmeasurementswitha bris disk in affecting the stability of the HR 8799 system.
much longer baseline taken with Hubble Space Telescope in First, we investigate if it is possible to delay the onset of
1998 for theouterthreeplanetsHR8799b, candd confirm instability for an initially highly unstable orbital configura-
2
tion. Then we consider a configuration with a long lifetime comprise a majority of the total mass of the debris. These
and determine whether adebris disk can cause instability. objects are not bright enough to individually detect in the
wavelengths that they emit. But estimating the total mass
of thedebris is important to determining thedynamics.
1.1 Resonant Structure of HR 8799 The halo and disk dust has been modelled by Su et al.
(2009).Thehalomassisestimatedtobe1.9×10−2M⊕ with
As discussed in the introduction, initial solutions tested
aradiusof upto1000 AU.Estimating thedustmassof the
by Go´zdziewski & Migaszewski (2009); Reidemeister et al.
other disk components is more difficult - particularly for
(2009);Fabrycky& Murray-Clay(2010)hadsuggestedthat
the cold outer disk because the inner and outer edges are
HR8799b, cand dweremost likely in 4:2:1 dualmeanmo-
not well known. In the case of the inner edge, it was first
tion resonances. Nearly all orbital configurations that re-
thought to be at 90 AU from temperature estimates while
main stable during a reverse integration for the estimated
the outer radius had been modelled at 300 AU to account
age of the system which also agree with the observed or-
foralltheobservationalconstraints,givingatotaldustmass
bitalelementsandnominalmassesfortheplanetsminimally of 1.2×10−1M⊕. The inner disk is quite warm at ∼150K,
calledforHR8799canddtobeina2:1MMR.Furthermore, allowingforaverylowdustmassestimateof1.1×10−6M⊕.
thedualmeanmotion 4:2:1 resonancefollowed bytheinner
Abriefsummaryofthemodelledparametersthatwereused
pair of HR 8799c and d in a 2:1 MMR allow for the largest
including assumed surface densities, inner and outer radii,
possible range of masses that still produce simulations sta-
minimumandmaximumgrainsize,aswellastotalmasscan
ble over the lifetime of the system. Alternate solutions in-
befound in Su et al. (2009) in their Table 3.
cluded a 2:1 MMR among the outer pair HR 8799b and
Recentsubmillimeterobservationshaveaddedsomead-
c as well as a few finely tuned solutions. In a later analy-
ditionalconstraints.Hugheset al.(2011)examinedHR8799
sis, Marshall et al.(2010) foundthatwhenthethreeplanet
and its debris disk with the Submillimeter Array at 880µm
configurationisplacedinthe4:2:1MMRwithlowplanetec-
- a wavelength which is suitable for examining larger dust
centricities, HR 8799 could survive the observed age of the
grains.Lowsignal-to-noisepreventedafullmulti-parameter
systemandpotentiallylonger.Re-reductionofHubblespace
modelling of thedustat thiswavelength buta combination
telescope astrometric data of HR 8799 by Soummer et al.
of the SMA data and spectral energy distributions (SEDs)
(2011),originallytakenin1998,robustlysupportstheprevi-
seem to rule out a narrow ring of dust and favor a broad
ous4d:2c:1b dualMMRhypothesiswithresultswhich have outerring starting with an inner edge at ∼150 AU.
only small departures from the exact integer period ratios
Scaling up from these modelled dust masses to a to-
in previously published data.
tal debris disk mass can be difficult. Small dust grains
This previous dynamical work precedes the discovery
are typically evacuated in debris disks by a combination
of the innermost planet HR 8799e by Marois et al. (2010).
of Poynting-Robertson drag and radiation blow-out (de-
However,simulationsbyMarois et al.(2010)haveindicated
pending on grain size) on very short timescales relative to
that the extra planet only places more restrictions on the
the age of an average debris disk. To constantly replenish
possiblesystemarchitectures.Inorderforthosesimulations
the dust in a system that is many millions of years old, a
toremainstableovertheminimumestimatedlifetimeofthe
collisional cascade is required (Safronov & Zvjagina 1969;
HR 8799 system, the planets e, d and c are required to be
Dohnanyi 1969; Williams & Wetherill 1994; Tanaka et al.
in a 4:2:1 MMR along with masses on the minimum end
1996;Kenyon & Luu1999;Kenyon& Bromley2001,2002).
of theestimated values. Notethat the4:2:1 dual mean mo-
It is assumed that the differential number distribution
tion resonance is still the most likely configuration for this
oftheparticlesinmassisapowerlawwhichtakestheform
system to be in if long term stability is desired - the differ-
encebeingwhichplanetsareinthemeanmotionresonance, dn(m)=Am−αdm (1)
theirprojectedmassesduetodynamicalconsiderations,and
or simply
themaximumageofthesystemwhich remainsstablegiven
these orbital elements. dn(a)=Ca2−3αda (2)
in radius. The steady-state solution depends on an α index
of 11. This parameter has been estimated empirically and
1.2 Distribution and Total Mass of HR 8799’s 6
numerically in the previously mentioned literature and, as-
Debris Disk
suming that α is determined only by the mass-dependence
The total mass and distribution of the debris found in the of the collision rate and that the model is self-similar, can
disk affects the dynamics, migration rate, extent of migra- beshowntobethatvalueanalytically(Tanaka et al.1996).
tion, and the smoothness of migration. However, the to- Toestimatethetotalmassofthedisk,weintegratethe
tal mass of the inner and outer debris disk along with the differentialnumberdistribution timesthemass of theseob-
dust halo is not well known for HR 8799. Estimating the jects at a constant density from theminimum tomaximum
total mass in debris within a disk is difficult to compute sizes in thecascade
accurately for both observational and theoretical reasons. amax
Solid(non-gaseous)matter,whichemitscontinuouslyinthe MT = M(a)dn(a). (3)
Z
µm tomm wavelengthsproducesnearly all of theradiation amin
(Beckwith & Sargent 1996). However, while the total cross Theconstantinournumberdistributioncanbesetbyusing
section of the dust particles is many orders of magnitude themodel of theouter debris disk by Su et al. (2009).This
larger than that of larger objects (m to km sized asteroids, modelincludestheminimumandmaximumgrainsizesthat
comets, and planetesimals), it is these later objects which were used in order to create a synthetic SED which was
3
matchedtotheobserveddata.Withthedustmassestimate (Quillen et al. 2007). Thirdly, we make use of a relation
from that model, and a minimum and maximum grain size which describes how the opacity scales with the radius of
of 10µm and 1000µm respectively, the constant C can be theobject,
computed.
a 3−q
To find the total mass of the disk we repeat the previ- τ =τ (6)
(a) d(cid:18)a (cid:19)
ous calculation only with a new a that corresponds to d
max
thelargest radiusobjectsinourcollisional cascade.Theto- whereq=3.5isequivalenttoα= 11 fromaboveanda and
tal mass that this integrand will yield is entirely set by the τ aretheradiusandopacityatasp6ecificradius/waveledngth
d
upper bound because there are many orders of magnitude
(Quillen et al. 2007). Substitutingthe abovethreerelation-
betweentheminimumandmaximum object sizes. Unfortu-
ships into our mass computation yields an integral that is
nately,determining amax is difficult. dependent only on the radius of the dust. Using 1000µm
Assuming the cascade operates for the system age dust particles as the largest grains we find that the opacity
Quillen et al.(2007)estimatedamax,assuminganalphapa- at thespecific dust particle radius is
rameterof11/6,thatusedonlyobservablepropertiesofthe
M 1
debris disk.As pertheir equation 16, τ ∼ d (7)
d A ρa
d
λ M 8/3 r −14/3
atop≈5.4km(cid:18)10µm(cid:19)(cid:18)M⊙⋆(cid:19) (cid:16)100AU(cid:17) awrheear,eρMisdthisetdheenstiottyalofmtahsesdoufstthgeradiunsst,,aAndisathisetthoetarladdiiusks
d
Q⋆ −5/3 t 2 h 10/3 of the largest dust grain in the model. Using dust masses
× D age
(cid:18)2×106erg·g−1(cid:19) (cid:18)107yr(cid:19) (cid:18)0.02(cid:19) from Su et al. (2009) as well as a disk area computed from
the inner and outer edges and a maximum grain size from
τ¯(λ) 2 f −2
× τ the same, along with a density appropriate for dust grains
(cid:20)10−2(cid:21) (cid:18) 4 (cid:19) ρ = 1.5−3.0 g , we compute a value for the opacity of
(4) τ ∼10−4. cm3
d
With these computed disk observables along with the
whereh,thediskaspectratio,andτ¯,thenormaldiskopacity
otherknownparametersforHR87999wefindana value
atwavelengthλ,areourdiskobservables.Otherparameters top
of a ∼ 1km. This a yields a disk mass estimate of
include M , λ, r, Q⋆, t , and f , which correspond to top top
thestellar⋆mass, obseDrvataiogen waveleτngth,theradiiat which Mdisk ∼150M⊕ or about one half of a Jupiter mass. Large
debrisdisksareoftenattributedmassesintheregionof50−
there is a break in the surface brightness profile, specific
energy,age, and an uncertainty factor. 100M⊕.
Whilethisestimatewouldplacethediskmassusablein
However,HR8799hasnoconstraintsonthescaleheight
oursimulations at of ordera Jupitermass or less, thereare
because the disk is not resolved. We therefore adopt the β
anumberofcaveatstoourestimates.Wenotethatouresti-
PictorisestimateforhbyQuillen et al.(2007).β Pictorisis
matefora islowcomparedtootherdisksbyQuillen et al.
alsoanAstarwhichhasasimilarageandmasstoHR8799 top
(2007)whichplacetheradiusatvaluesanywherefromafew
as well as an extended debris disk. Estimates for the dust
toafewhundredtimeslarger. Usinganestimatefora in
massaroundβ Pictorishavebeenfoundtobe7.8M ,or top
Moon line with those for β Pictoris and similar disks would yield
about 0.096M⊕ (Holland et al. 1998).This is similar to the
a significantly more massive disk. Also, we recognize that
0.12M⊕ suggestedbySuet al.(2009)forHR8799.Su et al.
the masses of the planets are extreme, with four planets of
(2009) also notes that theamount of excess emission in the
tenJupitermasseseach.Itisnotinconceivablethatthedisk
HR8799diskissimilartothatofotherdebrisdisksaroundA
maybesignificantly moremassive thanthoseobservedpre-
stars(SeeSu et al. (2006)’sstudyof theevolution ofdebris
viously. We also note that planetesimals with a radius of
disks around A stars.) See table 1 by Quillen et al. (2007)
a arethelargestobjectsthatcontributetothecollisional
and references therein for mass, age, and other parameters top
cascade.Significantmasscouldbefoundinother,moremas-
for β Pictoris.
siveobjectsthatwouldhavecollision timestoolongtocon-
We note that our biggest uncertainty is therefore in h
tributeto dust production. Finally, we note that the major
in the above formula, which goes as 10. This makes even
3 phenomenon discussed in the paper occur at masses of one
small errors in estimates of h have a large impact on a .
top Jupitermassorlessinthedebrisdisk.Itisnotrequiredfor
However,thegoal isnottodeterminetheexactmassofthe
ustoincludemuchofthemoremassivedisks,however,given
disk, only to determine what a reasonable range is.
the difficulty of estimating the mass and the general pecu-
To find a we must also compute the opacity at a
top liarityandsizeofHR8799,wehaveincludedfairlygenerous
specificwavelength.Theopacity can bemeasured bymodi-
debrisdisk masses.
fyingourtotal dust mass integral. Theopacity at aspecific
Asecond moregenericestimateof thelargest sized ob-
wavelength is a measure of the fractional area covered by
jects that can be grown was found by Cuzzi et al. (1993).
particles of radius a (i.e. the opacity depends on the num-
UsinganumericalmodelwhichusestheReynold’saveraged
berofparticlesperunitareatimesthecrosssectionalarea).
Navier-Stokes equations with both turbulence and full vis-
Inthiswaywecanrelatethenumberofparticlesofaspecific
cosity, they found that it was possible to create 10-100km
radius to the opacity and total disk surface area. Secondly,
sizedplanetesimals.Thiswouldplaceourdiskmassestimate
it is possible to equate the differential number distribution
of order ∼M .
to N by J
a Finally,wecan useourown solar system asareference
dN(a) ≡N(a), (5) point. The ‘Nice’ model required approximately 30-50 M⊕
dln(a) toexplaintheoutwardmigrationandeventualconfiguration
4
of our own system (Tsiganis et al. 2005). This corresponds code, it is not possible to run less than one block worth of
to a few tenthsof a Jupiter mass of debris. particles, even in the case where we wish a massless debris
Totaldebrisdiskmassesusedinoursimulationsvaried, disk. Furthermore, due to memory latencies, low particle
but had measurable effects at a Neptune mass, or about countsarenotrunefficientlyontheGPU.Currently,wedo
17 M⊕. This puts those simulated disk masses in a similar nothavetheability todisable theO(N2) kickcomputation
range to that of the ‘Nice’ model. Larger masses than the or a way to offload low particle counts to the CPU. This
‘Nice’ model are justified by way of the discussion above. is something that could be rectified in future revisions of
Giventheuncertaintiesinmeasuringa ,itisverydifficult the code. Additionally, because of the compilation process,
top
to determine if the largest total masses can be used in our our code is restricted to operating on machines which have
simulations. GPU’s attached.
Thetotalparticlenumberusedinoursimulations,1024, These two caveats have an effect on the number and
is also comparable to the 1000-5000 particles used in the size of theparameter space that we are able to explore.
‘Nice’ model simulations. Above we noted that the diame- The inefficiency of the GPU at low particle count and
terofobjectsatthetopendofthecollisional distributionis singleblockrestrictionhastheimpactofreducingthelength
anywherebetweenonetoseveralhundredkm.Ana value for which we can integrate our debris-less test case. We do
top
similar to β Pictoris of 180km is about 6.5 times smaller notwanttorunmultipleintegrators.Twohybrid-symplectic
thanthediameter ofPluto.a valuesnear1km would be integratorswillnotgetidenticalanswersunlessthecollision
top
approximately one thousand times smaller. This suggests integration method is the same. Rather than use another
that our simulated disks which are made up of at least a hybrid-symplectic integrator (like Mercury) to run debris-
Neptune mass in debris (those that are massive enough to less test cases, we continue to use our own code. However,
havemeasurableeffectsonthelifetimeofthesystem)would this does make it difficult to run simulations with as many
have planetesimals that are larger than those that are pre- orbital periodsas thosepossible with integrators which can
dictedbythecollisional cascade. Therefore, whileourplan- berun with very low particle counts.
etesimal masses are close tothose used in the‘Nice’ model, The GPU requirement also makes it more difficult to
they could be more massive than those in HR 8799’s ac- testawideparameterspace.Anodemustbeequippedwith
tual disk. This may result in more stochastic interactions aGPUinordertobeabletorunourcode.TheGPUclusters
between the planetesimals and planets. However, the plan- thatwehaveaccesstoaremuchsmallerthantheCPUclus-
etesimal masses suggested by the collisional cascade are a ters that are available. This limits the number of possible
distribution which couldreasonably involveboth smaller or trials.
larger planetesimals then those which we havesuggested. Last, we note that our integrator is truly O(N2) - all
particles feel the effects of all others. This is unlike many
integratorsavailableinthefieldwhichtypicallydonotcom-
puteplanetesimal-planetesimal effects.Thismakesoursim-
2 INTEGRATOR
ulations more precise, but does increase the arithmetic in-
All simulations were run with the software package QYM- tensity.
SYM1. QYMSYM is a GPU-accelerated hybrid second All simulations were run on either dual or quad core
Intel Core2 Duo architecture CPU’s with either GT200
order symplectic integrator which permits close encoun-
ters similar to the Mercury software package developed by of GF100 architecture NVIDIA GPU’s, both of which are
Chambers (1999). Like Mercury, QYMSYM uses a second CUDAcapabledoubleprecisionarchitectures.Thecodecan
easily be optimized to run on any Linux kernel 2.6+ distri-
order symplectic integrator to advance the positions of all
bution with appropriate GPU hardware.
particles in a simulation in the typical manner described
byDuncan et al.(1998)andLevison & Duncan(1994).Ex-
ploratory work on symplectic maps for N-body integrators
2.1 Accuracy, integration parameters and
was elucidated by Wisdom & Holman (1991). Additional
eccentricity corrections
analytical work on the formulations of the symplectic in-
tegrator can be found by Saha & Tremaine (1992) and by Timestep sizes were set to 0.016 (out of a possible 2π or-
Yoshida (1990, 1993). Unlike these first symplectic integra- bit) corresponding to 42 days in simulations including HR
tors but similar to Mercury, QYMSYM flags any particle 8799eand93daysinsimulationswithout.Weusedsmooth-
from an integration when it is deemed that it has a close ing lengths which correspond to radii at least a few times
approach to another massive particle during a time step. larger than the size of the planets assuming a planet and
Thesearethenintegratedseparatelywithanadaptivestep- planetesimal density of the order of 1g/cm3 and disregard
size conventional integrator. The criterion for closest ap- anysimulationduringwhichtheenergyconservation(given
proach is decided based on the Hill radii of the interact- by∆E/E)dipsbelow1.0×10−4.Additionally,weareonly
ing objects. QYMSYM uses the 4th order Hermite inte- interestedinthetimetotheonsetoffirstinstability,anddo
gration schemedetailed inMakino & Aarseth(1992)rather notneedtoconcernourselveswithlargerenergyerrorsthat
than the Bulirsch-Stoer algorithm used in Mercury. See may occur after a planet-planet interaction or ejection.
(Moore & Quillen 2011) for more details on the QYMSYM In a few simulations, very high eccentricities were de-
integrator. tected for planetesimals which had been ejected from the
We note that due to the nature of our CUDA based system during a close approach. Because our integrator
explicitly conserves angular momentum through the use
of f and g functions when solving Kepler’s equations (see
1 Seeauthor’sWebsiteforsourcecode. Moore & Quillen (2011) for more details), these large ec-
5
centricitiesandtheircorrespondinghighvelocitiesproduced massisconservedtoahighdegreeofaccuracyandwemain-
a drift in the location of the center of momentum. This is tain a satisfactory level of energy conservation throughout
most likely due to the way a hybrid symplectic integrator eithertypeof integration. Conservation of centerof mass is
typicallyupdatestheparticlespositions.Ifaparticlehasits typicallyontheorderofonepartin1012orbetter.Wethere-
position updated such that it is now in close proximity to fore present the results of the corrected and uncorrected
another particle, it would feel a large force. Typically this simulations together - using uncorrected simulations when
force is removed by reverse integration of the Hamiltonian no drift is detected in their outputs and corrected versions
and the particle is then moved to a more accurate integra- elsewhere.
tion regime (the Hermite integrator in our case). However, The above timestep choice and numerical energy
if by chance that particle is updated to be very close to integration error are comparable to those reported by
another, the assumptions made in the perturbation theory Fabrycky& Murray-Clay (2010).
usedtosplittheHamiltonianintocomponents,namelythat
Keplerian motion is dominant over interparticle forces, no
2.2 Simulation configurations of HR 8799
longerapplies.Thismeansthattheparticlewill notbecor-
rectly reverse integrated. This problem is common to many We set up two initial configurations to test in our simula-
symplectic integrators. tionsdetailedineachrespectivesection.Generally,giventhe
Once a planetesimal has been ejected with a very high importance of the 4:2:1 MMR in previous work, our simu-
eccentricity this phenomenon is more likely to repeat. The lations begin with this planetary configuration and havean
high eccentricity means that its angle of impact will be additionaldebrisdiskoralternativelyuseasomewhatmod-
nearly perpendicular to the motion of the other orbiting ified version of this planetary configuration. Minimally, the
objectsexteriortoitscurrentposition.Thisnewlylargeve- simulations havetheinnermost two planets in a 2:1 MMR.
locity relative to the size of the collision detection criteria Three planet (b,c,d) plus debris disk simulations
makes it more likely for a particle to be updated from out- were based off of the stable configurations discussed by
side the collision detection criteria to be in close proximity (Fabrycky& Murray-Clay2010),usingplanetmassesof10,
to thecolliding particle without being previously removed. 10 and 7 M for planets d, c and b.Four planet (b,c,d,e)
Jup
Includingalargercollision detectioncriterion orreduc- plus debris disk simulations are based off the predicted or-
ing the timestep size of the simulation could resolve this bitalelementsfoundby(Marois et al.2010)andusemasses
problem,butat greatly increased simulation time.Decreas- of7,7,7and5M forplanetse,d,candb.Ourdebrisdisk
Jup
ing timestep size goes directly with total simulation time. is simply a uniform distribution of 1024 equal mass parti-
Duetothehighvelocitiesachievedbytheparticles,itwould cles.Dependingontotaldiskmassdesired,theplanetesimal
need to be decreased significantly. Increasing the collision mass was modified accordingly.
detection radius, which is some factor of the Hill radius, Thislatersetofinitialpositionsandmassesreflectsthe
will force more particles to be offloaded to the significantly mostcurrentobservationalresultsbyMarois et al.(2010)by
slower Hermite integration routines. Not only is the Her- including the fourth planet and the corresponding reduced
miteintegratorO(N2)(andthereforesuffersfromnon-linear masses required from the dynamical simulations that were
increases in simulation time based on increasing particle run. The three planet configuration, while no longer repre-
count),butthenumberofparticlesincludedintheHermite sentativeofthosemorerecentobservationsandsimulations,
integrator will increase by a factor that goes with thecolli- is both practical and illuminating for several reasons.
sionalvolume-evensmallincreasesinthecollisiondetection First, we note that the planet which is not included
criteria can lead to significant increases in particle counts. in the three-planet configuration is the most recently dis-
Alternatively,increasingthesmoothinglengthcanpartially covered innermost object. Unless its presence destabilizes
alleviatethisproblembuthasthenegativeeffectofsmooth- the system, we expect its effect on the migration rate of
ing out many of the important planetesimal/planetesimal the outermost planet to be small in those simulations due
dynamics. to the proximity and mass of the debris disk to the out-
To correct for this occasional error without greatly in- ermost planet. Additionally, we only use the three planet
creasingsimulation time,wewroteanadditionalcheckthat configuration for simulations in which the system is arbi-
could remove planetesimals which had eccentricities above trarily placed in an unstable configuration from the begin-
an arbitrary threshold or those that collided with the star. ning. The four planet configuration has much smaller re-
Additionally, particles which are no longer bound and have gions of stability. While the addition of the fourth planet
a high enough semi-major axis are removed. Energy error would make it even easier for us to create an initially un-
checks are common in comparable literature (see appendix stableconfigurationwhichsharesorbitalelementssimilarto
of Raymond et al. 2011 for an example) as is planetesimal those that are observed, it makes it significantly more un-
removal dueto ejection. likely for a system to move from an unstable region to a
This check will have no effect on simulations in which stable region via orbital migration. Last, we recall that the
planetesimals arenotejected orejected withrealistic veloc- recent work by (Soummeret al. 2011) which reduced Hub-
itiesandeccentricitiesthatarestillcloseenoughtointeract ble Space Telescope observations suggests best orbital fits
orcollidewithanotherparticle.Onlysimulationswhichhave for HR 8799b,c,d that coincide with the 1:2:4 mean motion
particles with extremely large eccentricities will be quanti- resonance. This data agrees with previous work that fits
tatively altered. As we would expect, in simulations where only the three planets, but does not agree with dynamical
these particles are removed, it appears that the dynamics simulations which included all four planets (Marois et al.
arequalitativelyidentical,onlywithouttheaforementioned 2010). Marois et al. 2010 found that the innermost three
driftincenterofmomentum.Wealsonotethatthecenterof planets, HR8799c,d,e were the most likely planets to share
6
meanmotionresonanceswithbexcluded.Thereisasofyet islarge,thevalueofKgoestozerowhilewhenthedistance
no consensus on the possible orbital elements or dynamical between two particles is small, the value of K goes to 1 (or
structureofHR8799otherthanthemostlong-livedsystems viceversa).WhenKor(1-K)aremultipliedintotherespec-
having the innermost three planets in a 4:2:1 mean motion tiveseparated Hamiltonians described in (Moore & Quillen
resonance or at least havingtheinnermost two planets in a 2011),itpreventseitherfrombecomingtoolargeandbreak-
2:1 mean motion resonance. Therefore we assume that the ing the perturbation theory used to separate them. K is a
2:1 mean motion is the primary initial requirement for our function of Hill radii, but was created by trial and error.
simulated systems. Thedistributionofplanetesimalsemi-majoraxesisflat
Because we are measuring the effects of planetary mi- with probability independent of a within amin and amax.
gration on systemstability todetermineifallunstablecon- The initial eccentricity and inclination distributions were
figurations remain so, it is useful to run the three planet chosenusingRayleighdistributionswiththemeaneccentric-
system - it has larger regions of stability for the planets to ity e equivalentto twice themean valueof theinclination i
migrateinto.Duetothecomputationalintensityofthesimu- andi=0.01.Theinitialorbitalangles(meananomalies,lon-
lations and their corresponding time requirements, we were gitudes of pericenter and longitudes of theascending node)
only able to run on the order of hundreds of simulations were randomly chosen.
over a few months rather than tens or hundreds of thou- Agroupof150simulationswiththreeplanetswererun,
sands. Given this limitation, it proved difficult to migrate 58 using the eccentricity correction while 92 without. Disk
large numbers of four planet configurations into regions of mass was varied from 1.0×10−30 to 10 MJup, although we
stability. This lower simulation number and corresponding onlypresentresultsfromadiskmassof1.0×10−3MJup and
low number of stabilized configurations would introduce a larger here. 18 simulations with four planets were run, half
fine tuning problem to our analysis so we do not discuss using the eccentricity-correction while the other half with-
those in greater detail, instead focusing on thethreeplanet out.Thesesimulationsdifferedonlyintheinitialconditions
configurations for the stabilization through migration sce- of theplanetesimals.
narios. This fine tuning issue could potentially be resolved
viamorestableinitialconditionswhichwouldrequiresignif-
icantly less planet migration to move from regions of insta- 3 SIMULATIONS OF INITIALLY UNSTABLE
bilitytoregionsofstability.Thisissueisdiscussedinfurther PLANETARY ORBITAL ARCHITECTURES
detail in our results section.
In the three planet configurations the inner disk edge Does the addition of a relatively massive planetesimal disk
hasasemi-major axisof a =2.5 andan outerdiskedge allow for an increase in stability timescale? To answer this
min
of a = 7.0, corresponding to separations of 60 and 170 question we created an unstable initial configuration for
max
AU. In the four planet configurations the inner disk edge both the three and four planet configurations of HR 8799.
hasasemi-majoraxisofa =6.14andanouterdiskedge In either case, this was done by starting with actual ob-
min
of a = 20.69, corresponding to separations of 90 and served orbital elements found by Fabrycky& Murray-Clay
max
300AU.Thissecondsetmatchestheestimatedouterdebris (2010)andMarois et al.(2010)andreducingtheoutermost
disk’s observed inner and outer edges. planet’s semi-major axis in small increments until the sys-
The inner disk edge with reduced semi-major axis for temwasunstableonthetimescaleofthousandsoforbitsor
the three planet configuration was used to encourage rapid less.Inthethreeplanetsimulations,thisgivesaninitialcon-
planet-planetesimalcrossingsandthereforerapidmigration. figurationwithallthreeplanetshavingorbitalelementsthe
Asmentionedpreviously,weforcedourthreeplanetconfigu- same as those reported by Fabrycky& Murray-Clay (2010)
rationtobecomeunstableonveryshorttimescales.Because except that the outermost planet’s semi-major axis was re-
ofthis,werequirerapidmigration inorderforanyoutcome duced by 14%. For the four planet simulations, the plan-
other than planet-planet scattering to be a possibility. The ets’ orbital configuration was identical to those found by
inclusionofadebrisdiskataposition morecoincidentwith Marois et al. (2010) but the outermost planet has a semi-
that which is observed would be possible with either faster major axis reduced by 8%. This arbitrary and largely un-
or more GPUs. This is because configurations that are sta- stable configuration was chosen to both reduce simulation
ble on longer timescales - which allow for reduced amounts time as well as allow for pronounced effects by migration.
of planet migration to be necessary in order to stabilize an Due to the previously mentioned limitation in the number
unstable configuration - could be used. We could addition- of total simulations it is possible to run and the extremely
ally begin to experiment with the inclusion of the fourth limited regions of stability available to the four planet con-
planet as well as keep the outermost planet at its observed figuration,weoptedtorunfarmorethreeplanetsimulations
position rather than moving it in to encourage the system than four planet systems and present only that data.
to become unstable. We found both the removal of the in- Simulations both with and without eccentricity correc-
nermost planet and movement of the inner disk edge were tionsarepresentedsimultaneouslyandseenosignificantdif-
helpfulinfindingthesepost-migrationstableconfigurations. ference between them.
Other simulation parameters include the Hill fac-
tor for encounter detection and the K function (see
3.1 Results
Moore & Quillen (2011)), both of which are set to 2.0. The
K function is an arbitrary function which weights certain Infigure1 weplot thedifferencebetween thestability time
partoftheHamiltoniantobeintegratedinordertoallowit for each simulation and that of a system lacking a debris
to be broken up into evolutions operators which are other- disk. The stability time is the time to first planet-planet
wise not possible. When the distance between two particles encounter. The time difference is plotted as a function of
7
AU in semi-major axis depending on total disk mass. This
80
migration occurred overa few hundredthousand years.
70 If stable regions are small compared to unstable re-
gions,migration isunlikelytoputaninitially unstablecon-
60
figuration into a stable region. Only large migrations could
50 significantly increase the stability time. This is consistent
U] with what we see; rapid migration caused by a massive
A 40
a [ diskis abletosignificantly affect thelifetime of oursimula-
30 tions. Migration is expected to reduce planetary eccentric-
ities and this could increase stability (Lissauer & Stewart
20
1993;Ferna´ndez& Ip1996).Wesearched forsignsthat the
outermost planet’s eccentricity decreased duringmigration.
10
In figure 2 we see oscillations in eccentricity of the outer-
0 mostplanetwerelargethroughoutthesimulationduetothe
0 1 2 3 4 5
proximityandhighmassesoftheinnerplanets.Wetherefore
t [Myr]
attributetheincreasedstabilitytimetothewiderseparation
80 ratherthan any eccentricity damping.
Asaplanetmigratesoutwardsitcancrossmeanmotion
70
resonances with other planets. For simulations with short
60 lifetimes, we searched for evidence that resonant crossing
caused instability. We examined any relevant 1st, 2nd, 3rd
50
and 4th order resonance possible between the outer planet
U]
A 40 and innermost planet as well as the outer planet and the
a [ middle planet. We saw no large semi-major axis or eccen-
30
tricity jumps caused by resonant crossing so they are un-
20 likely to be the cause of instability in the simulations with
shorter lifetimes. Resonances near the outer planet’s posi-
10 tionareprimarily 2nd or3rd orderandcausesmallereccen-
tricityjumpsduetotheweakerdependenceonmass(Quillen
0
0 0.2 0.4 0.6 0.8 1 2006).
t [Myr] As discussed in our section on simulation parameters,
our choice of unstable configurations was done in part to
Figure 2. In this figure we show the evolution in semi-major
reducetheamountofsimulation time.Withaninitialplan-
axisofanexamplethreeplanetconfigurationovertime.a)Shows
etary configuration so far from a region of stability, a large
theevolution intimeofthesemi-majoraxisofthethree planets
amountofmigration isrequiredinordertomovetheplanet
overtheentire5Myrsimulation.b)Showstheevolutionintime
of the semi-major axis of the three planets over the first Myr. toamorestableregionofphasespace.Becauseofthis,only
Inb) wealsoinclude errorbars that indicate the maximum and amassivedebrisdiskcanmigratetheoutermostplanetsuf-
minimum separation (given by a(1+e) and a(1-e) respectively) ficiently far to increase thelifetime.
of the planets given their eccentricity. This second image more Anotherconcernwiththissetofsimulationsisthatwith
clearlyshows themigrationof the outermost planetplanet from aplanetesimalcountof1024andatotaldiskmassofatleast
its initial position at 58.6 AU to approximately 65 AU over the aJupitermass,interactionsbetweenplanetesimalsandplan-
firstfew hundred thousand years. In this example the total disk
etswouldcausestochasticmigration(Zhou et al.2002).We
mass was 1.6MJup and the configuration remained stable over wouldhaveexpectedsuchlargeplanetesimalmassestocause
theentire5Myrsimulation.
instability rather than increase the lifetime of the system.
Therefore, wehavenoreason tosuspect that stochastic be-
haviorisa problem.Because thetypicaldisk mass required
diskmass.Simulationsarerunforamaximumof5Myrs.If to stabilize the system is unrealistically large there may be
noencountershadtakenplaceinthattime,theyareplotted no need to resolve the disk more accurately. As mentioned
as upperlimits in the figure. previously, an approximate upper limit for debris mass is
In figure 2 we plot the evolution in time of the semi- thought to be on the order of 10−3M⊙, or approximately
major axis of the three planets for a sample simulation. In one Jupiter mass. For HR 8799, this is effectively the mini-
this case the total disk mass simulated was 1.6M . mum amount of mass required for any observable effects in
Jup
Figure 1 shows that at low disk masses the difference oursimulations.
in lifetime is near 0. Only massive disks with masses of ten
percent or more of the outermost planet can remain stable
until the end of the simulation. Given enough mass, migra-
4 SIMULATIONS OF INITIALLY STABLE
tion via planetesimal scattering can take a tightly packed
PLANETARY ORBITAL ARCHITECTURES
configuration of HR8799 b, c, d and e, or in the case of fig-
ures1 and 2, simply b, cand d,and migrate theoutermost In this second set of simulations we examine the impact of
planetintoamorestableconfiguration.Thetotalmigration a planetesimal debris disk on the stability time of a con-
for theoutermost planet in system configurations which re- figuration of planets which was stable for about 150k or-
mainedstableoverthe5Myrsimulationrangedfrom5to15 bits (or around 7 Myr) sans debris disk. In these simula-
8
5
4.5
4
3.5
3
]
r
y 2.5
M
[ 2
t
∆
1.5
1
0.5
0
-0.5
0.001 0.01 0.1 1 10
disk mass [M ]
Jup
Figure 1. This figure shows the difference in time until a planet-planet encounter for a group of HR 8799-like simulations with a
debrisdiskfromanidentical simulationwithnodiskmassversustheamountofmassinthatparticulardisk.Anintentionallyunstable
solution was used in order to quickly show the effects of a planetesimal debris disk. All simulations use orbital elements taken from
Fabrycky&Murray-Clay (2010) for the inner two planets (which are in a 2d:1c resonance). The outer planet, HR 8799b, has a semi-
majoraxiswhichisreducedby14%fromapproximately68to58.6AUinordertoachieveaplanet-planetcloseapproachwithinacouple
hundredthousandyears.∆trepresentstheincreaseordecreaseintimebeforethefirstplanet-planetencounter.Simulationswerehalted
afterthefirstplanet-planetencounter orafter5Myr.
tions we begin with all four planets situated in an orbital Matching the exact age found by Marois et al. (2010)
configuration based off the current estimated positions by isdifficult despitethefact that thereisnolonger arequire-
Marois et al.(2010).Thisfourplanetconfigurationisgener- ment of planetesimals with mass. This is due to the way
ally extremely unstableunless situated in a4e:2d:1c MMR. theintegrator operates. Specifically, it is not possible to in-
Asbefore,smoothinglengthsaresettosizesafewtimesthe tegrate less than a certain number of particles, depending
radii of the planet to prevent the unrealistic forces possible on compile time and hardware restrictions. For more de-
if two particles form a very tight binary or collide within tails on theexact natureof thecode, we refer thereader to
whatwouldhavebeenconsideredanapproximateplanetary Moore & Quillen (2011). However, while the control simu-
radius. Energy error is typically even lower in these simu- lation is not as long lived as those found by Marois et al.
lations with a ∆E/E at 1.0×10−5 or 1.0×10−6 due to (2010), it is still stable on long time scales (105 orbits).
reduced planetesimal masses. We stopped simulations after
their first planet-planet encounter regardless of the even-
tual fate of the interacting planets. Stability time in the
Thissimulationisthencomparedto18simulationswith
histogram is thetime to first planet-planet encounter.
Neptune mass disks. These 18 simulations are identical in
We began by simulating the four planets with an ini-
terms of planetesimal mass, total disk mass, initial outer
tialorbitalconfigurationcomparabletothatbyMarois et al.
and inner disk edge, and initial planetary orbital configu-
(2010) with massless debris and plotted the resonant argu-
ration. However, the planetesimals initial orbital configura-
ments. We see constrained resonant angles which indicate
tions havebeen randomized 18 different ways.
the presence of a 2:1 MMR between the inner two plan-
ets, shown in figure 3. While Marois et al. (2010) did not
directly show the resonant argument to illustrate that HR
8799d and e were in resonance, our plot is similar to Fig- Thehistogramshowninfigure4showsthedistribution
ure10byFabrycky& Murray-Clay(2010)althoughforHR inlifetimesoftheeighteensimulationseachwithaNeptune
8799e and d rather than d and c. The integrator confirms massdisk.Figure5showstheevolutionintimeofthesemi-
thatthisplanetaryconfigurationisinresonanceassuggested major axis of all four planets of a sample simulation. The
by Marois et al. (2010). errorbarsindicatetheminimumandmaximumseparations.
9
6
3
2 5
1 4
de 0 3
φ
-1 2
-2 1
-3 0
0 1 2 3 4 5 6 7 0-1 1-2 2-3 3-4 4-5 5+
time [Myr] t [Myr]
3 Figure 4. Here we show a histogram of time to the onset of
instability,whichwedefineasplanetinteractionorejection.The
2 diskmassisaboutoneNeptuneinmass.Theplanetsareinitially
inaconfigurationthatisstableforapproximately7Myrs.Wesee
that the distribution of lifetimes is very broad and many of the
1
simulationshadlifetimeslessthanthis.Thereductioninlifetimes
wasaslargeas∼6Myrs.
de 0
φ
-1
80
-2
70
-3 60
0 0.5 1 1.5 2 2.5 3 3.5
50
time [Myr] U]
Figure3.Hereweshowtheresonantanglefor2:1meanmotion a [A 40
resonance forHR8799d andeas functionof time.a)For asim- 30
ulation without a debris disk. b) For a simulation with a debris
20
diskthathasamassof1/100ththeouterplanet.Theinitialcon-
ditionsfortheplanets werestablefor7Myrs.Inb)theextreme 10
values of φ slowly begin to increase over the length of the sim-
ulation. This effect is not present in a). The resonant island is 0
0 0.5 1 1.5 2 2.5 3
sosmallthat evenadiskmassof1%theouter planetiscapable
of effecting the systems dynamics. Theresonant angle plotted is t [Myr]
givenbyφe=2λd−λe−̟e. Figure 5. In this figure we show the evolution in time of the
semi-majoraxis of an example 4 planet configuration. This par-
ticular simulation became unstable in a little over 3 Myr. We
4.1 Results
includeerrorbarsthatindicatethemaximumandminimumsep-
aration (given by a(1+e) and a(1-e) respectively) of the planets
In figure 4 we see that in 16 of 18 simulations the stability
given their eccentricity. In all four-planet simulations the total
timehasdecreasedbyatleast2Myrsfromtheinitial‘stable’
diskmasswasequivalenttoNeptune.
configuration which had its first planet interactions after 7
Myrs.Weseeabroaddistributionoflifetimes.Eventhough
aNeptunemassdiskisonlyabout1%themassoftheouter Figure3byGo´zdziewski & Migaszewski(2009)demon-
most planet it has an effect on thestability time. strates why a planetesimal disk with a total mass around
WithaNeptunemassdiskeachsimulatedplanetesimal that of Neptune can have such a pronounced effect on the
has a mass of about 8 Plutos. Debris disk models estimate stabilitytime.Inthisplot,darkregionsindicatehighlylikely
the top of the collisional cascade to contain similar (but configurationsofeccentricityandsemi-majoraxisforplanet
slightly smaller) size bodies. For example, Wyatt & Dent d,whereasyellowregionsindicatestrongly chaoticsystems.
(2002) argued for a distribution of planetesimals in the Fo- Hereweseethatsemi-majoraxischangesoflessthan0.1AU
malhautsystemwitharangeofsizesbetween4and1000km. canmoveplanetdfromastableregiontoastronglychaotic
Asdiscussedpreviously,thissecondvalueroughlycoincides one. Similarly, minor changes in eccentricity can also have
with Pluto’s diameter within a factor of a few of what a large effect. Fabrycky& Murray-Clay (2010)’s Figure 7
we used in our simulations. Stochastic perturbations are showsasimilarbehavior.Notethatachangeincurrentsep-
expected in real systems; however they will be somewhat aration of planet d by 0.05 in the figure (corresponding to
smaller than those present in our simulation do to our fac- semi-majoraxischangeof1.22AU)candecreasethestabil-
tor of 8 difference in mass of planetesimals. ity time by as much as four orders of magnitude. A three
10
orderofmagnitudechangein stability timeispossible with disks(10%themassoftheoutermostplanet)couldcausein-
a semi-major axis increase of approximately half an AU.In creasesinstabilitytime.Infourplanetconfigurationswhich
eitherexample theregions of stability havesharp edges be- were stable over long timescales without a debris disk, de-
tween stable and chaotic solutions. bris disks of only a Neptune in mass (1% the outer planet)
Wewouldexpectadiskmassof1%themassofaplanet cancauselargedecreasesinlifetimes.Weattributethissen-
tobeabletomigrateaplanetroughly1%ofthesemi-major sitivity to the small size of HR 8799’s resonant region of
axisifamajority ofitsangularmomentumistransferredto stability.
the planet via scattering. If this were the case than a Nep- Whiletheamountofdiskmassrequiredforsystemsta-
tune mass disk would be able to migrate the outer planet bilizationinoursimulationsisunrealistic,itisonlyrequired
(which has a mass of 5MJup at 68 AU) approximately 0.5 tobethislargeduetotheinitialconditions.Inordertorun
AU.Thisismorethanenoughtovarythelifetime bymany a large number of simulations, a highly unstable configu-
ordersofmagnitude.Thereforewelookedforthenumberof ration that rapidly had planetary encounters was required.
planetesimals which had become orbit crossing by the time Givenan initial configuration which is just outsidearegion
eachsimulation hadbecomeunstable.Inall16simulations, ofstability,itwouldbepossibletocausesufficientmigration
atleast15%ofthediskmasswasorbitcrossing.Thiscorre- to allow an unstable configuration to become stable with a
sponds to a minimum of 2.5 M⊕. By starting a simulation more reasonable disk mass.
neartheedgesofaregionofstability,weseethatevenlower Similarly, while thedisk properties required for system
mass debrisdisks can affect thestability time. disruption in our simulations are reasonable in both total
mass and planetesimal size, the mass distribution for the
planetesimals is not a true mass distribution. In order to
4.2 Pulling HR 8799’s planets out of resonance
keeptotalparticlenumberdownwerestrictallplanetesimals
We continuously monitor the planets’ orbital elements to have an identical 8 Pluto sized mass. We do not believe
throughout all simulations to see if they are in resonance that this had a large effect on the distribution of stability
with one another. We do this by noting whether or not the timesofthesystemwhenadiskisincludedbecausearealis-
resonant angle librates around 0 or π. We looked for 1st, tic mass distribution could encourage even more stochastic
2nd, 3rd and 4th order mean motion resonances that could migration if there are even a few planetesimals of greater
exist near each of the planets’ semi-major axis at varying mass.
times. We did not observe any strong two body resonances These results suggest that HR 8799 may be destined
beyond the 2:1 MMR between e and d, a much weaker 2:1 for eventual planet scattering. The very high masses of the
MMR between d and c, and no clear resonant angles for b. planetscausesregions ofstabilitytoberelatively small and
Weattributethis to theoverwhelmingly large gravitational makes it possible for low mass debris disks to pull the sys-
perturbations caused by the massive planets as well as the tem out of its mean motion resonances and induceinstabil-
factthatmostofthenearbyresonancesfortheplanetswere ity.Wealso see thatmigration ismuchmore likelytomake
ofhigherorder.Thispreventedbehaviorsimilartothe‘Nice’ thesystem unstablethan migrate thesystem to a region of
modelduringwhichmigrationcausesresonancecrossingand stability.Itmaybethecasethatthedebrisdiskhasalready
corresponding eccentricity increases. Finally, we see that in beenresponsibleforremovingthesystemfrom amaximally
all simulations when the planets e and d are pulled out of stableregionandiscurrentlyneartheboundaryofinstabil-
the2:1 MMR thesystem rapidly becomes unstable. ity.Allofthesepossibilitiestendtoinduceinstabilityrather
In figure 3a we see that the resonant angle for the 2:1 thanstability andsuggest thatHR8799 maybeheadedto-
MMR between HR 8799e and d over the 7 Myr simulation wards instability, possibly sooner than would otherwise be
remains largely unaffected until the end of the simulation. predicted.
However, in figure 3b we note that the extreme values of In this studywe haveused initial conditions consistent
φ begin to increase over time, leading to e and d moving with observed positions of the planets when attempting to
out of resonance and therapid disruption of the system. In determine whether a stable orbital configuration could be
this particular example, the system’s lifetime is decreased made unstable by a debris disk. As we have shown here, a
by about 50%. It is possible the debris disk is responsible low mass disk can cause evolution of the planetary config-
for pulling HR 8799 out of resonance and so, reducing its uration. Consequently, in the past the planets could have
lifetime. been in a different configuration. By exploringdifferent ini-
tial conditions future investigations could explore scenarios
for the past evolution of HR 8799 that would be consistent
with its currentconfiguration. Forexample, iftheplanetes-
5 CONCLUSION
imal disk is causing the system to move out of resonance,
Inthispaperwe discusshowaplanetesimal debrisdiskcan in thepast it mayhavebeen in amore stableregion of this
effect the stability of the multiple planet system HR 8799. resonance,insteadofonitsboundary.Increasingsimulation
Twoquestionsareconsidered;whetheritispossibletodesta- length and particle count could also increase the resolution
bilize a stable planetary configuration with a planetesimal of our simulations.
debrisdiskandconversely,whetheritispossibletostabilize Alternatively,it appears that similar simulations could
an unstable configuration with a planetesimal debris disk. beusedtoputupperlimitsonplanetesimaldebrisdiskmass
We examine both three planet (b,c,d) and four planet con- ifitisassumedthatthedebrisdiskhasanegligibledynam-
figurations (b,c,d,e). ical effect. In the case of our simulations of HR 8799, plan-
In three planet configurations which were unstable on etesimal disk masses would need to be much smaller than
short timescales without a debris disk, only massive debris one Neptunein mass.
