Exoplanets Bouncing Between Binary Stars

0.83 MB·English

by  Nickolas Moeckel

Mon.Not.R.Astron.Soc.000,1–11(XXXX) Printed1February2012 (MNLATEXstylefilev2.2)
Exoplanets Bouncing Between Binary Stars
Nickolas Moeckel1⋆, Dimitri Veras1†
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA
Accepted XXX.ReceivedXXX;inoriginalformXXX
2
1
0
ABSTRACT
2
Exoplanetary systems are found not only among single stars, but also binaries of
n
widely varying parameters. Binaries with separations of 100–1000au are prevalent in
a
theSolarneighborhood;attheseseparationsplanetformationaroundabinarymember
J
maylargelyproceedasifaroundasinglestar.Duringtheearlydynamicalevolutionof
1
aplanetarysystem,planet–planetscatteringcanejectplanetsfromastar’sgrasp.Ina
3
binary,themotionofaplanetejectedfromonestarhaseffectivelyenteredarestricted

three-body system consisting of itself and the two stars, and the equations of motion
]
P of the three body problem will apply as long as the ejected planet remains far from
E the remaining planets. Depending on its energy, escape from the binary as a whole
. may be impossible or delayed until the three-body approximation breaks down, and
h
further close interactions with its planetary siblings boost its energy when it passes
p
close to its parent star. Until then this planet may be able to transition from the
-
o space around one star to the other, and chaotically ‘bounce’ back and forth. In this
r paper we directly simulate scattering planetary systems that are around one member
t
s of a circular binary, and quantify the frequency of bouncing in scattered planets. We
a
find that a great majority (70 to 85 per cent) of ejected planets will pass at least
[
once through the space of it’s host’s binary companion, and depending on the binary

1 parameters about 45 to 75 per cent will begin bouncing. The time spent bouncing is
v roughly log-normally distributed with a peak at about 104 years, with only a small
2 percentage bouncing for more than a Myr. This process may perturb and possibly
8 inciteinstabilityamongexistingplanetsaroundthecompanionstar.Inrarecases,the
5
presence of multiple planets orbiting both stars may cause post-bouncing capture or
6
planetary swapping.
.
1
0 Key words: binaries:general–planet-starinteractions–planetsandsatellites: dynam-
2 ical evolution and stability–scattering
1
:
v
i
X
1 INTRODUCTION tions subsequent to planet–planet scattering and planetary
r ejection (Ford et al. 2005), and/or perturbations from the
a Binary starsystemswhichcontain exoplanetsvarydramat-
stellar companion (Barnes et al. 2011).
ically.From thefirst discoveriesofexoplanetsorbitingstars
Generally, systems with multiple stellar and planetary
inbinarysystems(Butler et al.1997)toapulsarplanetwith
componentspresentrichdynamicallaboratorieswhichcould
a stellar companion (Sigurdsson et al. 2003) to the recent
inform planetary formation and evolution theories. Indeed,
circumbinary“Tatooine”-likeplanetsannouncedbytheKe-
Desidera & Barbieri (2007) estimate that binary compan-
plermissionteam(Doyle et al.2011;Welsh et al.2012),ex-
ions within just 100 − 300 AU are likely to play a role
oplanetshavebeenobservedtoexistinanintriguingvariety
in the formation and evolution of planetary companions.
of multiple-star systems. A few years after the discovery of
In the Solar neighbourhood, such binaries are prevalent;
thefirstmulti-planetsystemaroundamainsequencestar,υ
Raghavan et al.(2010),inthemostcompletesurveyourlo-
And (Butleret al. 1999), Lowrance et al. (2002) discovered
a stellar companion to υ And at ∼ 750 AU. Subsequently, cal volume within 25 pc, find approximately half of solar-
type stars in multiple systems, with a broad peak in the
Curiel et al. (2011) announced a fourth planet orbiting the
(roughlylognormal) separation distributionstretchingfrom
primary (see also McArthuret al. 2010). The orbital archi-
afewtoafewthousandau.Atseparationsaboveafewtens
tecture of this system may be explained by secular interac-
of au, these local multiple systems have the same planet
hosting frequency as single stars, setting the stage for po-
⋆ E-mail:[email protected] tentialinteractionsbetweenstellarandplanetarydynamics.
† E-mail:[email protected] AnexampleisMarzari et al. (2005),who studiedplanetary
(cid:13)c XXXXRAS
2 Moeckel & Veras
systemsundergoingplanet–planet scatteringfrom arounda
2 C =5 C =3.7
starwithabinarycompanion.Whilethefocusofthatwork
was on the orbits of the unejected planets, they found that
1
the binary companion participated directly in the planet
scattering, and a few per cent of the ejectees collided with
0
thebinarypartner;interaction between thebinarystarand
theplanets of themost direct sort. -1
Here, we investigate planetary transfer between stellar
companionssubsequenttoplanet–planetscattering. Incon- -2
trast toMarzari et al. (2005),wearemost interested in the -2 -1 0 1 2
planetsthatareejectedfromtheirhoststarratherthanthe
C =3.5 C =3.1
planetary system left behind. Our motivation is best illus-
tratedbyconsideringthemotionoftestparticlesinthepla-
narcircularrestrictedthreebodyproblem(planarCR3BP).
The planar CR3BP, whose early form was established by
Euler and Lagrange in 1772 (see Barrow-Green 1997), de-
scribes the motion of a massless particle subjected to the
force of two massive bodies moving in circular orbits such
thatallbodiesarecoplanar.TheCR3BP(initsplanarform
or its extension into three dimensions, the 3D CR3BP, e.g.
Szebehely & Peters 1967; Marchal 1990) is well suited to
Figure 1. Curves of zero velocity in the dimensionless synodic
describing the motion of an extrasolar planet orbiting one
frameofacircularbinarysystem.Inthisframethestars(redand
or both stars of a circular binary system; even the most
bluedots)arefixed inspace. Thegrayregions show inaccessible
massive planets, at about 13 MJ (Burrows et al. 1997), are regionsfortestparticleswiththelabeledvalueoftheJacobianin-
<10percentasmassiveasthelowest-knownexoplanethost tegral.Thecontinuouslinesshowexampleorbitsattheseenergies,
star mass (as of 31 January, 2012, 0.13 M⊙ for KOI-691)1. each originating in the space around the more massive star. At
However,inthevastmajorityofcases,theplanet/starmass evenlowervaluesoftheJacobianconstant,theforbiddenregions
ratioismuchlessthan1percent.Hence,inbinarysystems contractaroundtheL4andL5Lagrangepointsbeforevanishing.
with multiple planets, the motion of a planet will be dic-
tated by the CR3BP as long as the perturbations from the
plot, the planet may orbit either the primary or secondary,
other planets are negligible. This situation arises when one
and further, freely travel in between each star’s sphere of
planet is sufficiently far away from the others, which may
influence via the connecting region centered on the stars’
occur subsequentto planet-planet scattering.
The position and velocity of this planet, as well as the L1Lagrangepoint.Thebottomtwoplotsillustratedifferent
waysinwhichtheplanetcanescapethesystem;theC =3.1
masses and separations of the two stars, determine what
plotdemonstrateshowtheplanetcanorbitbothstarsbefore
regimes ofmotion arepossible fortheplanet.Jacobi (1836)
escaping. While this illustration is in the framework of the
showedthatthesequantitiesarerelatedbyaconstantofthe
planar CR3BP, such transiting orbits exist in the 3D case
motion, C, such that
as well (e.g. G´omez et al. 2004).
A planet forming around one star of a widely2 sepa-
C =x2+y2+2 µp + µs −x˙2−y˙2−z˙2, (1) rated binary will be in the limit shown in the upper left
(cid:18)rp rs(cid:19) panel of Figure 1, although the C value of this planet will
constantly change if there are other planets around one or
wheretheposition(x,y,z)andvelocity(x˙,y˙,z˙)oftheplanet
both stars. However, this change is significant only during
arecomputed with respect totheorigin ofa frame rotating
closeapproacheswiththeotherplanets.Planet-planetscat-
with a mean motion of unity about the center of mass of
teringmaycauseasuddendecreaseinC,enoughtoopenup
the stars such that µ = GM , µ = GM , µ +µ = 1
p p s s p s thegap around L1 seen in theC =3.7 plot butnot enough
and where r and r represent thedistance of theplanet to
the primarypand secsondary, respectively. Later, Hill (1878) toenableescapethesystemthroughtheexteriorL2 andL3
Lagrange points, as seen in the C =3.5 and C =3.1 plots.
linked given values of C with regions of space which are
In this case, the planet may ‘bounce’ between the stars. If
allowed or forbidden.
theplanetisscatteredfarenoughawayfromtheotherplan-
Figure 1 shows examples of these regions in the planar
ets,thenC mightmaintainavaluethatallowsthisbehavior
case,i.e.z=z˙ =0inequation1.Supposeamasslessplanet
tocontinueonappreciabletimescales. Here,weexplorethis
is initially orbitingtheprimary.Then,thetrajectory of the
possibility,andcharacterizetheprospectsandconsequences
planet is given by the continuous line, the primary and the
forplanetaryescapeandcapturethatisinternaltoabinary
secondary are represented by the large red and blue dots
stellar system.
on the left and right respectively, and the gray areas show
Because planets are not massless, and we wish to treat
regions that the planet cannot access. The relative masses
planet-planetscatteringaccurately,wecannotrelysolelyon
of the stars here are µ /µ = 0.3. In the C = 5 plot, the
s p
planet is restricted to orbiting the primary. In the C =3.7
2 ‘Widely’forthesepurposesislooselydefinedasasystemwith
semi-majoraxessuchthatabinary ≫aplanets,sothattheplanets
1 http://exoplanets.org/ areinitiallyorbitingasinglestar.
(cid:13)c XXXXRAS,MNRAS000,1–11
Exoplanets Bouncing Between Binary Stars 3
e=0.999binary Pythagoreanproblem Kozaioscillations
-4 0 4 100
a 10-8 4 t=0–10 t=10–20
(cid:144)
Èa 10-1
DÈ 10-9 0 -e
0 250 500 750 1000 -4 1 10-2
t=20–30 t=30–40 10-3
(cid:144)e 10-11 0 5 10 15
Èe
DÈ 10-12 180
0 250 500 750 1000
s
s t=40–50 t=50–80 ee 120
e r
e g
r 0.5 e
eg d 60
d (cid:144)
(cid:144) i
v
3
0 0.0 0
1 0 250 500 750 1000 0 5 10 15
orbits time(cid:144)Myr
Figure 2. Tests of the integrator used inthis paper. Left: Integration of a planet with eccentricity e=0.999 for 1000 orbital periods.
Theerrorsinenergyandangular momentum aresmallandperiodic,withnosystematic drift.Thereis,however, asmalllineardriftin
the longitude of perihelion. Middle: Paths of the stars in the Pythagorean problem through t=80. In each panel the positions of the
bodiesattheendofthetimeintervalareshownwithdots.Right:AtestpresentedinNaozetal.(2011):asystemwithaninnerbinary
consisting of a 1 M⊙ star and a 1 MJ planet on an a = 6 au, e = 0.001, i = 64.7◦ orbit, with a 40 MJ companion on an a = 100
au, e = 0.6, i = 0.3◦ orbit. The longitude of pericenter of both companions are initially zero. The inclination and eccentricity of the
innerplanetareplottedasitundergoes Kozaioscillations,includingswitchesfromprogradetoretrogradeorbits.Theagreementtothe
integrations inNaozetal.(2011)isexcellent.
theanalytical structureof theplanar or3D CR3BP forour 2 THE CODE
explorations. In particular, planets often reach ejection ve-
The N-body code we use is built upon a basic fourth-order
locities via a series of encounters that send them progres-
Hermiteintegrator,writtenbyP.Hut3.Thestraightforward
sively farther from their host star. This rules out introduc-
structure of the code makes it easy to adapt for different
ingtestparticleswithaspectrumofenergiesintheCR3BP.
purposes,althoughthiscomesatthesacrificeof fine-tuning
Ifejections were dominantlytheresult ofsingle encounters,
that could speed up the code if it were tailored to specific
the CR3BP would be more readily and directly applicable,
problems. Ourmodifications and tests are detailed here.
althougheventhenwewouldneedtointroduceastochastic
forcingeffectwithinacertaindistanceofthehoststartoac-
count for possible interactions with the remaining planets.
2.1 Integration Algorthim
We therefore turn to direct numerical simulation. Because
N-bodycodesdesignedtostudyplanetarysystemsoftenrely Inits implicit form, theHermite method istime symmetric
onthepresenceofonemassivecentralobjectthatdominates andexhibitsnodriftinenergywhenintegratingperiodicor-
themotionoftheotherbodies,weconstructedasimpleand bits. In practice, an explicit approximation to the implicit
flexible integrator to investigate planetary bouncing. This methodisusuallyused,following aPEC(predict–evaluate–
codeanditstestsaredescribedinSection2.InSection3,we correct) cycle. This approach is commonly applied to stel-
describe ensembles of simulations with the code, and show lardynamicsproblems, wheresupplementarymethods(e.g.
that bouncing is thedominant behaviour of planets leaving KS regularisation, Kustaanheimo & Stiefel 1965) are used
abinarysystemviaplanet–planetscattering.Wediscussthe to avoid systematic errors in long lived systems with pe-
results and future extensions of this work in Section 4, and riodic motion. In order for the method to be suitable for
then briefly conclude in Section 5. studyingplanetarysystems,theimplicit versionoftheHer-
mite integrator must be used. To achieve this we make two
modifications to the base code. First, we iterate over the
evaluationoftheforceanditstimederivativeandthecorrec-
tiontothepredictedvalues,convertingtoaP(EC)nmethod
(Kokuboet al.1998).Thisiterationprocessconvergestothe
3 Thatcodeisavailableathttp://www.artcompsci.org.
(cid:13)c XXXXRAS,MNRAS000,1–11
4 Moeckel & Veras
implicit Hermite solution, at the expense of multiple force 250aubinary Mh=1.0MŸ,Mc=1.0MŸ
evaluations. In practice, we use a single iteration, n=2. 103
With fixed timesteps, the P(EC)n method would yield u 102
a
excellent energy conservation. When timesteps are allowed (cid:144) 101
tobevariable,thetimestepsmustbechosentobesymmet- a 100 st
1.0 o
rical functions of the system’s state at the beginning and h
endofthetimestepinordertoretaintheenergyconserving e 0.5
0.0
featuresoftheHermitescheme.Whenalltheparticlesshare 0 1 2
thesametimestep,asinourcode,thisisstraightforward to 103
implement, following the method of Hut et al. (1995). This u 102
approachrequiresiterationsoverthetimestepchoice,andwe (cid:144)a 101 on
againuseoneiterationhere.Asingletimesteptheninvolves a 100 ani
oneiteration fortheP(EC)n integrator,wrapped insidean- 1.0 mp
otheriterationforthetimestepdetermination.Theresultis e 0.5 co
arobustandflexibleintegratorforsmall-N systemswithno 0.0
0 1 2
preferred dominant force or geometry.
The timestep requested by each pair of particles i,j time(cid:144)Myr
withrelativeseparationr ,velocityv ,andaccelerationa
ij ij ij
is proportional to the minimum of their estimated free fall Figure3.Semi-majoraxesandeccentricitiesofplanetsina250
and collision times, the latter under unaccelerated motion: au, equal-mass binary. The orbital elements are plotted relative
δt = ηmin[(a /r )1/2,r /v ]. The system timestep for to whichever star’s space the planet is found in at any time, as
ij ij ij ij ij
a given particle state is then the minimum of this quantity definedinthetext.Toppanelsareorbitsaroundthehost,bottom
panelsarearoundthecompanion.
over all particle pairs; this estimate is used in the timestep
symmetry iteration. In the simulations presented here we
set the control parameter η = 0.02, yielding typical rela- lations (Kozai 1962), leading to flipping of the inclination
tive energy errors of the order 10−8, and relative angular from prograde to retrograde. The setup is as follows: a 1
momentum errors of theorder 10−11. M⊙ primary is orbited by a 1 MJ planet with a = 6 au,
◦
e=0.001, and i=64.7 . There is also a 40 M companion
J
with a = 100 au, e = 0.6, inclined 65◦ with respect to the
2.2 Tests
innerplanet.Astheeccentricityandinclinationofthecom-
Inordertoverifytheintegrator’scapabilitiestofollowplan- panions oscillate, the inner planet oscillates between pro-
etary dynamics, we have performed several tests. We show grade and retrograde orbits. This is the same system pre-
the results in Figure 2. In the first, we integrate a binary sentedinfigure3ofNaoz et al.(2011),towhichourresults
system with eccentricity e=0.999 for 1000 orbits. We plot can be compared. Those authors integrated the system us-
thesemi-majoraxis,eccentricity,andlongitudeofperihelion ingaBurlisch-Stoerintegrator,andtheagreementwiththat
overtherunoftheintegration.Thenumericalerrorsarepe- result is excellent.
riodic with minimal systematic drift, with the exception of
a small linear error in thelongitude of perihelion. Similarly
good results are obtained with less extreme eccentricities. 3 THE SIMULATIONS
The example shown has a secondary to primary mass ratio
of0.1,andsimilar resultsareobtainedforratiosfrom unity Weperformedseveralsetsof300simulationseach,exploring
to approximately zero. a limited but illustrative range of binary parameters. Each
ThesecondtestisofthePythagoreanproblem(Burrau simulation is run for 10 Myr.
1913).WithG=1,threebodiesofmass3,4,and5areini-
tially at rest at thecornersof aPythagorean right triangle,
3.1 Initial Conditions
eachbodyatthecorneroppositethesideofthetrianglecor-
respondingtoitsmass.Foradetailedaccountofthecomplex We follow a standard planetary setup for scattering sim-
three body dance that ensues, refer to Szebehely & Peters ulations, following Marzari & Weidenschilling (2002). Each
(1967),thefirstauthorstoresolvethesystemtoits conclu- adjacentpairofplanetsisseparatedbyK mutualHillradii,
sion: after one body is ejected from the system, while the so that the planet i and its outward neighbour j have semi
other two recoil as a binary.Those authors used a regulari- major axes aj =ai+KRHij, with
sation schemetodealwiththedifficultiesofthemanyclose m +m 1/3 a +a
encounters in the system’s evolution. Using the parameter R = i j i j,
Hij (cid:18) 3M (cid:19) 2
for thetimestep control that we use throughout thispaper, ⋆
our code integrates the system through virtually the same and the chosen value of K determines the insta-
trajectory asSzebeheley&Petersfind,with relativeenergy bility timescale (Marzari & Weidenschilling 2002;
error of the order 10−9. Following Szebehely & Peters and Chatterjee et al. 2008). Planetary masses are drawn
Hut et al. (1995), we integrate the system with velocities from a power law mass function with f(m) ∝ m−1.1,
reversed from t = 62 back to t = 0. Achieving similar re- between 0.3 and 5.0 M . We take K =4 and an innermost
J
sultsasthoseauthors,werecovertheinitialconditionswith semi-major axis of 3 au, which together with our mass
errors entering at thethird decimal place. spectrum yields an instability timescale of the order 104
Finally, we integrate a system undergoing Kozai oscil- yr, although the possibility of long-term stable systems
(cid:13)c XXXXRAS,MNRAS000,1–11
Exoplanets Bouncing Between Binary Stars 5
Table1.Parametersandscatteringresultsforthesimulatedbinarysystems.Allbinarieshavezeroeccentricity.Themajorityofscattering
planetarysystemsexperience abouncingplanet.
a/au Mh/M⊙ Mc/M⊙ Nscatter Nswitch ejections collisions
Ntotal nswitch nbounce Ntotal nswitch nbounce
250 1.0 1.0 249 221 236 204(86±7%) 182(77+−65%) 84 14(17+−64%) 0(0+−20.2%)
250 1.0 0.3 231 199 220 192(87+−76%) 138(63+−65%) 67 4(6.0+−42..79%) 0(0+−20.7%)
250 0.3 1.0 209 200 226 185(82±6%) 174(77±6%) 43 10(23+−170%) 3(7.0+−63..88%)
1000 1.0 1.0 235 185 215 180(84+−76%) 133(62+−65%) 69 1(1.4+−11..42%) 0(0+−20.2%)
1000 1.0 0.3 227 157 207 156(75+−76%) 73(35+−54%) 73 0(0+−20.5%) 0(0+−20.5%)
1000 0.3 1.0 165 116 162 112(69±7%) 88(54±6%) 24 2(8.3+−85..33%) 0(0+−70.7%)
MhandMc arethemassesofthehostandthecompanionstars.Nscatter isthenumberofsystemsthatscattered,outof300setsofinitial
conditions. Nswitch isthenumber of planetsthat crossthe L1 thresholdatleastonce. Inthe ‘ejections’and ‘collisions’columns,nswitch
and nbounce are the number that cross L1 at least once and more than once, respectively; the percentages include Poisson confidence
limitscorrespondingtoa1σ Gaussianinterval.
250 au binary
107 107 107
Mh=0.3MŸ,Mc=1.0MŸ Mh=1.0MŸ,Mc=1.0MŸ Mh=1.0MŸ,Mc=0.3MŸ
yr yr yr
(cid:144) 106 (cid:144) 106 (cid:144) 106
n n n
o o o
ni ni ni
a a a
p 105 p 105 p 105
m m m
o o o
c c c
nd 104 nd 104 nd 104
u u u
o o o
ar ar ar
e 103 e 103 e 103
m m m
ti ti ti
102 102 102
102 103 104 105 106 107 102 103 104 105 106 107 102 103 104 105 106 107
timearoundhost(cid:144)yr timearoundhost(cid:144)yr timearoundhost(cid:144)yr
1000 au binary
107 107 107
Mh=0.3MŸ,Mc=1.0MŸ Mh=1.0MŸ,Mc=1.0MŸ Mh=1.0MŸ,Mc=0.3MŸ
yr yr yr
(cid:144) 106 (cid:144) 106 (cid:144) 106
n n n
o o o
ni ni ni
a a a
p 105 p 105 p 105
m m m
o o o
c c c
nd 104 nd 104 nd 104
u u u
o o o
ar ar ar
e 103 e 103 e 103
m m m
ti ti ti
102 102 102
102 103 104 105 106 107 102 103 104 105 106 107 102 103 104 105 106 107
timearoundhost(cid:144)yr timearoundhost(cid:144)yr timearoundhost(cid:144)yr
Figure4.Scatterplotsofthetimespentinthespacearoundeachstarforallplanetsthatbouncedacrosstheboundarybetweenthem.
The zero point for the clock is the moment the planet first crosses over the threshold; all the time spent around host prior to being
scattered across to the companion isnot counted. The slightgridded appearance of points atsmall times is dueto our timeresolution
for these plots, whichis 250 years. Thediagonal graylineshows equal time spent around each star. The errorbars at the bottom and
leftshowthemedianofthedistribution,andthe16thand84thcentiles.
(cid:13)c XXXXRAS,MNRAS000,1–11
6 Moeckel & Veras
exists with this setup. Eccentricities are initially zero, an approximation. With the history of each planet’s stellar
◦ ◦
inclinations are uniformly distributed between 0 and 1 , residencerecorded,wetallythenumberoftimeseachplanet
and all the phase angles are randomized. Changing the crossestheL1threshold,aswellastheamountoftimespent
planets’ eccentricity and mutual inclinations will alter the around each star. When a planet crosses the L1 threshold
instability timescale, but otherwise yield similar outcomes a single time, we say that it has ‘switched’; if it crosses the
(e.g. Juri´c & Tremaine 2008). threshold repeatedly, it has begun to ‘bounce’.
Inthisinitialwork,thebinarysysteminitially haszero The main questions we wish to answer here are: what
eccentricity, and has semi-major axes of 250 or 1000 au. fraction of scattered planets transition to the space sur-
Small changes to these values are possible as the system rounding the companion star, and how much time do they
evolves and planets are ejected or collide with the stars. spend around the companion, or bouncing between the
Thesemi-majoraxisvaluesarechosentoliewithinthepeak stars?
of the local binary separation distribution (Raghavan et al.
2010), but large enough that the planet formation process
3.3 250 au binary results
around a single star is probably not drastically affected by
the presence of its binary companion. The masses of the We turn first to the results for a 250 au binary. Table
planet hosting star and its binary partner (henceforth the 2.2 provides a summary of the outcomes from our sim-
‘host’ and ‘companion’, respectively) are taken to be (1.0, ulations. For the systems with 1 M⊙ hosts the number
1.0), (1.0, 0.3), and (0.3, 1.0) M⊙. We assume coplanarity of systems that go unstable and scatter (N ) dur-
scatter
between the the orbital plane of the binary and the plan- ing the 10 Myr of our simulations is about 230–250, i.e.
etary system; this should have only a secondary effect on 15 to 20 per cent of the runs have not yet scattered4.
theresults, sinceplanetary ejections from scattering donot The systems with a 0.3 M⊙ host seem to be more stable,
occur strictly in theplane of theplanetary system. with about 30 per cent unscattered. Some difference is ex-
Aplanetisconsidered tobeejectedifitsdistancefrom pected; the instability timescale is dependent on the mass
the binary center of mass exceeds 5×104 au, a factor of function (Chambers et al. 1996; Marzari & Weidenschilling
∼2 less than a typical distance at which galactic tides can 2002). Keeping the planetary mass function fixed while re-
causeescapeoveratypicalmainsequencelifetime(Tremaine ducingthehostmassisequivalenttoshiftingtheplanetsto
1993) but well beyond the region in which an appreciable higher masses. The sense of the change in instability times
numberof planetsmaybescattered ontostablewide orbits weseeisthesameasinpreviouswork:relativelymoremas-
(Veraset al. 2009). Collisions between bodies occur when siveplanetsarestableforlongertimeswhenspacedasfixed
their radii are detected to overlap; the planets are all 1.3 multiples of their mutualHill radii5.
Jupiter radii, and the stars are 5 Solar radii. Using this Looking at the number of planets that switch, we see
large stellar radius enables us to avoid the regime where that a large fraction of scattered planets will pass through
tidal interactions overly affect the dynamics. Mergers be- thecompanionstar’sspace;thenumberofswitchingplanets
tween bodies are momentum and mass conserving. N is about 200 for all three binary setups. The over-
switch
whelming contribution to this category comes from plan-
ets that end up ejected from the binary and eventually re-
3.2 Analysis and Characterisation
moved from the calculation, shown as n in the ‘ejec-
switch
To assign a planet at any instant to one of the two stars, tions’ columns. For instance, 204 of the 221 planets that
we use thefollowing simple scheme: first, rotate thesystem switch in the equal-mass binary runs end up ejected, with
about the binary’s center of mass so that the stars are on similar fractionsfor theothertwomass ratios.Crossing the
the x-axis, and the binary’s angular momentum is in the zˆ L1 threshold is symptomatic of a planet that is leaving the
direction.Thelocation ofthebinary’sL1 Lagrangepoint is system. We illustrate one such system in Figure 3, which
determined,andthisiscomparedtothepositionsalongthe plots the semi-major axes and eccentricities of the planets
x-axisoftheplanets;whicheversideoftheL1divideaplanet with respect to whichever star’s airspace they are in at the
liesonisitsinstantaneoushost.Planetswithdistancesfrom time.Thescattering in thiscasetakesplace veryearly, and
thebinary centerof mass greater than twice theseparation the planet shown by the blue line begins to transition back
of the binary components are tagged so that a planet on a and forth. After an uninterrupted period of about 0.5 Myr
longloopingexcursion,oronethathasbeenejectedbutnot around the companion, the final few bounces back to the
yetremovedfromthecalculation,isnotcountedasbeingin hostbringitintofurtherinteractionwiththeotherplanets,
orbit about eitherstar. and its next trip through the companion’s space ends with
While quite simple, this scheme is motivated by the escape from thesystem.
symmetry of the allowed regions of space in the CR3BP Looking further at the ejected planets, note that the
(see Figure 1). After subtracting off the orbital motion of majority (about 85 per cent) cross over to the companion
thebinaryandappropriately nondimensionalising allveloc-
ities,masses,anddistances,wealsotrackanapproximation
to the Jacobian integral of each planet in the 3D CR3BP. 4 Thisfractionisconsistentwiththeresultsofanidenticalplan-
etary setup presented inMoeckel&Armitage (2012), integrated
In doing this we treat the host and its associated planets
withmercury(Chambers 1999).
as a single body at their barycenter. While this serves to 5 Chambersetal.(1996)suggestthattheinstabilitytimescalein
estimate when a planet has sufficient energy to transition
a three planet system should vary independently of the planets’
between topologies of allowed orbits, it is not exact: the mass function when the planets’ spacing is scaled as the one-
non-zero mass of the planets and small eccentricities in the fourthpoweroftheirmasses,ratherthantheone-thirdpowerin
binaryorbitinducedwhenaplanetisejectedrestrictthisto theHillradiusspacingschemeusedhere.
(cid:13)c XXXXRAS,MNRAS000,1–11
Exoplanets Bouncing Between Binary Stars 7
on their way out of the system: ejection from the host star
occurs predominantly through the space around L1 rather 0.2 1000au
than the exterior Lagrange point. Furthermore, a majority 250au
ofejected planets(about 65to75percent)notonly switch al
ot
but bounce, traversing between stars multiple times. Esca- Nt 0.1
pers are not simply passing through the companion star’s (cid:144)
N
space, but rather orbiting there for some time before cross-
ing back to the host, and so onward chaotically until they 0.0
leavethesystemorinteractfurtherwiththerestoftheplan-
102 103 104 105 106 107
ets. Planets that ultimately collide with something are less
interesting from a switching standpoint; very few collision timearoundstar(cid:144)yr
partners are involved in intrabinary excursions.
The answer to the question of what percentage of
ejected planets cross over to orbit the companion is evi- 0.2
dently an interestingly large number:something like 85 per
al
centofescaperswillswitch.InFigure4weaddressthetime ot
theyspendinthespacearoundeachstarbyplottingthethe Nt 0.1
(cid:144)
time around the host versus the time around the compan- N
ion for all stars that make the jump and return, with the
0.0
clock starting from the moment they first pass through the
L1 connection–i.e. we do not count time spent around the 10-2 10-1 100 101 102 103 104
hostpriortoscattering.Thediagonallineintheplotsshows
timearoundstar(cid:144)P
equal time around each star; the points cluster around this binary
lineinaclearcorrelation.Thepointlyingnear1Myraround
Figure5.Thedistributionoftimesthatbouncingplanetsspend
eachstarintheequalmasscaseisthesystemshowninFig-
aroundeach star inyears (top panel) andbinaryorbitalperiods
ure3.Amoretypicalbouncingsystemisshorterlived,with
(bottom panel). Timearoundthehoststarareshownwithsolid
fewerlongtermresidenciesabouteitherstar.Thereisabias
lines,andtimearoundthecompanionasdashedlines.Themedi-
towardsspendingtimearoundthemore massivestarin the
ansofthedistributionsareshownwitharrowsatthetopofeach
unequal mass binary cases, but this is not a strong effect. plot.
Thereareafewoutliersineachcasewithlargetimesaround
the host. These are all cases where a planet made excur-
sionstoorbitaroundthecompanion,andatsomepointfur- 2C=2.56 C=3.24 C=3.64
ther interactions with the planets around the host changed
1
the Jacobian integral of the wandering planet, re-isolating
0
itaroundthehost.TheseadjustmentstotheJacobian inte-
gral are not possible when theplanet is in orbit around the -1
companion. -21000aumedian 250aumedian 250au75-centile
-2 -1 0 1 2
n 1.0
3.4 1000 au binary results ctio
a
Overall, theresults for the1000 au binary are quitesimilar fr
to the 250 au case. This is largely expected; the spacing of ve 0.5
ati
theplanets,setatafixedinnerradiusof3au,togetherwith ul 1000aubinary
the physical radii of the bodies are the physical quantities m 250aubinary
u
thatgeneratetheejectionvelocitiesoftheplanets.Theratio c 0.02.0 2.5 3.0 3.5 4.0
oftheplanets’semi-majoraxesaroundthehosttothebinary
Jacobianintegralatswitch
separationisverysmallforbothofthebinarysizes,andthe
orbitsontowhichtheplanetsarelaunchedwhentheybegin
Figure 6. The bottom panel shows the distribution of the Ja-
bouncingshouldoriginatefromsimilarpointsinspacewhen
cobian integral C at the moment a planet first switches to the
the binary separation is taken as the length scale.The only
companion for the equal mass binaries.The top panel shows, in
remaining explanation for changes in the outcomes is the
grey,therestrictedregionsofspaceinthesynodiccoordinatesys-
relation between the ejection velocities, a physical quantity
temoftheCR3BPforthreevalues:themedianvaluesofthe1000
setbytheplanetarydynamics,andtheJacobianintegral.A and250audistributions,andthe75centileofthe250auruns.
planetwithagivenphysicalvelocityandpositionrelativeto
itshostwillhavealowerJacobian integralinamorewidely
spacedbinary,sincethedimensionlessvelocityistiedtothe tains thesame ejection velocity from around its host in the
binary period. 1000aubinarywillhaveamore‘wideopen’escapethrough
What differences there are in the outcomes compared thecompanion’s exterior Lagrange point; this easier escape
to the 250 au binaries are consistent with this picture. The plausibly reducesthe tendencyto bounceback to the host.
percentageofejected planetsthat switch areessentially the Thetimesspentaroundeachstarafterbouncingbegins
same for both binary separations, but the bouncing per- arelongerforthe1000ausystems.Ifthesystemsweresimple
centage is lower for the 1000 au binaries. A planet that at- rescalingsofanidenticaldimensionlesssetup,thisdifference
(cid:13)c XXXXRAS,MNRAS000,1–11
8 Moeckel & Veras
shouldbeeliminatedbyrescalingthetimetothebinaryor- 1000aubinary Mh=1.0MŸ,Mc=1.0MŸ
bitalperiod. Wetest thiswith theequalmassbinary cases.
In Figure 5 we show the distributions of time spent around 103
u 102
eachstar,whichareroughlylog-normal,forbothbinarysep- (cid:144)a 101
arations. With time measured in the binary orbital period a 100 st
(the lower panel) the histograms become more similar, al- 1.0 o
h
though the relative positions of the histogram peaks swap e 0.5
places.Onceaplanethasbeenlaunchedbytheinternaldy- 0.0
namics of the planets around the host into the three-body 0 1 2 3 4 5 6 7 8 9 10
system consisting of the binary and the ejected planet, the 103
planetsbouncingin the250 au binaryare longer lived rela- u 102
tive to thebinary’s dynamical timescale. This is again con- (cid:144)a 101 on
sistent with the idea that the tighter binary’s planets find a 100 ani
1.0 p
themselveswith ahigherJacobian integral, andless readily m
able to escape through the region surrounding the exterior e 0.5 co
Lagrange points. 0.0
0 1 2 3 4 5 6 7 8 9 10
Wecandirectlytestthisideabyexaminingthedistribu-
tion of theJacobian integral C at themoment of a planet’s time(cid:144)Myr
first switch. In Figure 6 we plot this distribution for the
equal mass cases for both 1000 au and 250 au binaries. We Figure7.Semi-majoraxesandeccentricitiesofplanetsina1000
also showtheintersection with thebinary’sorbitalplaneof au, equal-mass binary. The orbital elements are plotted relative
thezero-velocity surface associated with threevaluesof the to whichever star’s space the planet is found in at any time, as
Jacobian integral: themedian valuefor each binary separa- definedinthetext.Thestarthatisejectedfromthehostsettles
tion, and the 75th centile of the 250 au binary. This latter into a quasi-stable orbit around the companion. Top panels are
orbitsaroundthehost,bottompanelsarearoundthecompanion.
valueisaboutthe95thcentileofthe1000auruns.Thecrit-
ical value of the Jacobian integral when the L1 connection
opens is C1 = 4.0; at C2,3 = 3.46 L2 and L3 open up. As
0.6
noted earlier, thecalculated value of C is not exact, dueto
small errors introduced by the inconsistencies with the 3D 0.5
CR3BP in our actual setups. Despite these minor limita- 0.4
tions, the difference in the distribution of C is clear; more
e 0.3
planets around the 250 au binary have C2,3 < C < C1, D
where access to the outside universe through the exterior 0.2
Lagrangepointsisimpossiblewithoutfurtherperturbations 0.1
from theremaining planets. In contrast, the great majority
0.0
ofswitchingplanetsinthe1000aubinaryarefreetoescape -1.5 -1.0 -0.5 0.0 0.5 1.0
thesystem.
20
Since we observe a dimensional universe, the longer
physical timescale associated with the 1000 au binaries is es 15
of potential interest. In particular there are a handful of e
r
casesthatspendlongperiodsaroundthecompanion,either g 10
e
switching back and forth repeatedly or in rare cases ending d
up in a quasi-stable orbit about the companion. An exam- (cid:144) 5
i
pleofthelatteroutcomeisshowninFigure7.Inthiscasea D
planetbouncesbetweenstarsforseveralMyr,beforesettling 0
aroundthecompanionatabout3Myrandstayingtherefor -1.5 -1.0 -0.5 0.0 0.5 1.0
theremainderofthesimulation.Westressthatthisoutcome
Da (cid:144) au
israre,andisnotreliablystable.Thesemi-majoraxisofthe
planetaround thecompanion isabout 200au,and thecon-
tinuedinfluenceofthehoststarisevidentintheeccentricity,
Figure 8. The change in orbital parameters of a single planet
which oscillates between about 0.5 and 0.95.
around the companion with e0 =0, a0 =10 au, measured at 10
Myr.Theserunswerewitha250aubinaryandM1=M2=1.0
3.5 Response of a planet around the companion M⊙.
The orbits of stars passing through the companion’s space,
even those that are only bouncing for a very short amount of 10 au, and zero eccentricity. This nearly blank slate of a
oftime,arecharacterisedbyhigheccentricities.Ofpotential planetary system records the effect of the wandering plan-
interestistheeffectthismayhaveonplanetsordiscsaround etsonplanetarysystemsnearthelikelyterrestrialandgiant
thecompanion.Weexaminethisverybrieflybyconducting planet forming region around the companion.
simulationsidenticaltothe250au,equalmassbinaryruns, The changes to the semi-major axis, eccentricity, and
but with the addition of a single planet around the com- inclination at 10 Myr of this target planet are plotted in
panion. This planet has a mass of 1 M , a semi-major axis Figure 8. The majority of systems are relatively unaltered
J
(cid:13)c XXXXRAS,MNRAS000,1–11
Exoplanets Bouncing Between Binary Stars 9
planet around companion 250aubinary Mh=1.0MŸ,Mc=1.0MŸ
103
107 u 102
Mp =1.0MŸ,Ms =1.0MŸ (cid:144)a 101
yr a 110.00 ost
(cid:144) 106 h
ry e 0.5
a
0.0
d
n 0.0 0.5 1.0 1.5
co 105 103
e u 102
s (cid:144)a 101 on
nd 104 a 100 ani
u 1.0 p
o m
ar e 0.5 co
me 103 0.00.0 0.5 1.0 1.5
ti
time(cid:144)Myr
102
102 103 104 105 106 107 Figure 10. Switching followed by planet displacement in a 250
au, equal-mass binary with a single planet in orbit around the
time around primary (cid:144) yr
companion.Theorbitalelementsareplottedrelativetowhichever
star’sspace the planet isfoundinatany time,as defined inthe
text.Thestarthatisejectedfromthehostinteractswithaplanet
Figure 9. As in Figure 4, the times spent around each star for
originally around the companion, before ejecting it and settling
bouncing systems in a 250 au binary, but with a single planet
intoastableorbitaroundthecompanion.Theoscillationsinthat
orbiting the companion. Top panels are orbits around the host,
planet’seccentricityareduetointeractions withthehoststar.
bottom panelsarearoundthecompanion.
eccentricbinariesontheseresultstootherwork,theresults
bytheintrudingplanets,butoftheorder10percentofthese
ofoneexamplecase hintatthemagnitudeandcharacterof
planets experience changes in semi-major axis of around 1
thedifferences to expect.
au,increases ineccentricitytotensofpercent,andinclina-
We re-ran the 250 au equal mass case with a bi-
tion from a few to a few tensof degrees.
nary eccentricity of 0.75. The first point to note is that
These changes to the single planet’s orbit imply some
the closer approaches of the stars perturbed the planetary
changetotheorbitoftheintrudingplanet,possiblyaltering
systems so that only 14 cases remained unscattered af-
thetimespent around each star. In Figure9 we showthese
ter the 10 Myr runs. This destabilizing effect is described
times.Thereisverylittlechange,exceptforthepresenceof
in Holman & Wiegert (1999) and noted in Marzari et al.
afewsystemsthatspendmoretimearoundthecompanion,
(2005).Similarlytothelatterauthors,wefindmoreplanets
filling in the space above the diagonal reference line. These
ejected(perscatteredsystem)athighereccentricity;∼1.25
arenotverynotable,withtwoexceptions.Twicetheplanet
persystemcomparedto∼0.95inthezeroeccentricitycase.
ejected from the host displaces the planet around the com-
The fraction of ejected planets that cross L1 is similar, 89
panion, taking its place and ejecting the companion’s orig-
per cent with e=0.75. The number of bouncing planets is
inal planet from the system. These are the only two cases
virtuallythesameat76percentcomparedto77percentin
wheretheplanetoriginallyaroundthecompanionisejected.
thecircularruns.Thetimesspentaroundthehostandcom-
One of these runs is shown in Figure 10. The interactions
panionaresmaller,byaboutafactorof2.Theintroduction
between the incoming planet and the companion’s planet
ofmoderateeccentricity appearstopreservethequalitative
occuroverabout0.5Myr,culminatingintheescapeofthat
andmany quantitativeaspects of thiswork, thoughthede-
planet after briefly crossing over to thehost system.
tailed dependence of these results on the binary eccentric-
ity deserves further attention. We note that as the binary
periastron decreases to becomparable totheplanet’s semi-
3.6 Non-zero binary eccentricity
major axis the companion may take a more direct role in
Forsimplicity,wehavedealtexclusivelywithzeroeccentric- the scattering (Marzari et al. 2005), perhaps leading to a
itybinariesinthisfirstwork.Theeccentricitydistributionof qualitative shift in the outcomes.
long-periodbinariesspansallvalues;Raghavan et al.(2010)
findaflatdistributionuptoe∼0.6,withthedropoffabove
that value probably due at least in part to observational
4 DISCUSSION AND FUTURE DIRECTIONS
bias.Eccentricityinthebinarymayexpectedtochangethe
resultsofthiswork;Hillstabilityofatestparticleisafunc- A natural extension to this work would be to remove var-
tion of the massive bodies’ eccentricity (Marchal & Bozis ious restrictions of the CR3BP. For example, if the two
1982; Holman & Wiegert 1999; Szenkovits& Mak´o 2008), stars are not bound to one another, but are hyperboli-
and this dependence could alter the timescale or qualita- cally scattering off of each other, then one may utilize
tivenatureofthescatteredplanet’strajectoryinthebinary. the parabolic or hyperbolic restricted three-body problem
While will we defer a more complete study of the effect of (Luk’yanov2010;Sari et al.2010)inordertodeterminethe
(cid:13)c XXXXRAS,MNRAS000,1–11
10 Moeckel & Veras
resulting motion of the planet. However, as touched upon panionthatweperformedhereshowsthatasmall butnon-
in section 3.6, a closer extension to our study would in- negligiblepercentageofcompactplanetarysystemscouldbe
volveremovingthecircularorbitrestrictionbutkeepingthe perturbed by the intrusion of a scattered planet. There are
stars bound to each other. If the stellar masses orbit one many possible consequences of this: changes to the eccen-
another on an unchanging, eccentric orbit, then the Jaco- tricity of outer planets can propagate inwards through the
bian integral no longer exists but the 5 Lagrangian points system(Zakamska& Tremaine2004);two-planetsystemsin
still do (Szebehely & Giacaglia 1964; Szebehely 1967). The formally Hillstableconfigurationscouldbepushedintoun-
coordinate system traditionally adopted for this situation stable architectures; mean motion resonances could be bro-
is nonuniformly rotating and pulsating; the zero-velocity ken;andrichersystemsthanthesingleplanetthatwesimu-
surfaces change with time. The consequence is a system lated could beinduced intoscattering themselves, resulting
which is difficult to describe in simple closed analyti- inwholesaletransferofplanetsbetweenthetwostars.These
cal forms (e.g. Contopoulos 1967; Delva& Dvorak 1979; are thetopics of ongoing investigations.
Gawlik et al. 2009;Hou & Liu 2011)but can still admit es-
capeandcapture(e.g.Bailey1972;Huang & Innanen1983;
Llibre & Pinol 1990;Mak´o & Szenkovits2004).
Inallofthesesituations,astarcannotcaptureaplanet 5 CONCLUSIONS
and retain it over the star’s main sequence life without
Via direct integration of a planetary system around one
some additional dissipative force. In the traditional pho-
member of a circular binary system, we have examined the
togravitational restricted three-body problem (Radzievskii
idea that planets can be sent bouncing between the stars
1950), this force arises from the star’s radiation on a mass-
during planet–planet scattering. The internal dynamics of
lessdustparticle. Foraplanetinabinarysystem,however,
theplanetarysystemcanincreasetheenergyofoneormore
a non-negligible force might arise from additional planets
planets,enablingittocrossovertothecompanionstarfrom
(as briefly discussed in Section 3.5) or from mass loss or
its host. Wefind that this is in fact a dominant outcome of
gain from thestarsorplanet (Bekov1988,1991;Luk’yanov
scattering in a binary: about 70 to 85 percent of planets
2009; Letelier & daSilva 2011). This situation is particu-
ejectedfrombinariesintheseparationrange250to1000au
larly relevant after one or both of the stars have left the
willspendtimearoundbothstarsbeforeleavingthesystem
mainsequence,howeverthebouncingtimescaleswefindhere
entirely,and about 45 to 75 percent will bouncerepeatedly
(typically less thanafew Myr) are shorterthan stellar evo-
back and forth. The amount of time spent bouncing is as-
lution timescales. For a discussion of planetary transfer in
tronomicallyshort:afewtotensofbinaryperiods,although
post main sequence binary systems, see Kratter & Perets
withsomesystemssurvivinghundredstothousandsofperi-
(2012, in preparation).
ods.Weproposethattheimpact ofthisbouncingwill most
A more likely dissipative mechanism in the planet likely be seen in the response of the companion star’s own
formation context is interactions with discs, either of planetary system, or with the inclusion of un-modeled dis-
protoplanets, debris, or gas. Scattering can occur when sipative forces, such as discs, around thecompanion.
significant amounts of gas are still around the star
(Moeckel & Armitage 2012), and if the two stars of the bi-
nary are at similar stages in their lives a planet scattered
over to the companion star could encounter a significantly ACKNOWLEDGMENTS
dissipative medium. Ifdissipation is sufficient tostrand the
Our thanks to Dan Fabrycky for a clear and useful ref-
marauding planet around the companion star, the orbit
eree’s report, and to Kaitlin Kratter and Cathie Clarke for
would likely be a wide one, offering a way to produce gi-
comments on the draft. Large portions of this work were
ant planets in the innerregions of a disc and send them off
performed using the Darwin Supercomputer of the Uni-
toawiderorbitaroundabinarycompanion.Thissortofsce-
versityof Cambridge High Performance ComputingService
nario is probably most efficiently explored via a dissipative
(http://www.hpc.cam.ac.uk/), provided by Dell Inc. using
prescription of some sort in the CR3BP rather than direct
Strategic Research Infrastructure Funding from the Higher
integration of a hydrodynamic disc. Finally, while we have
Education FundingCouncil for England.
focused on planetsthemselvesscattering across tothecom-
panion star in this work, other scattering events are likely
in planet formation. Scattered planetesimals and exoplan-
etary analogues to our own outer Solar System structures
REFERENCES
(e.g.theKuiperbelt,orscattereddiscobjects) maybeable
toparticipate in jumping between stars in asimilar fashion Bailey J. M., 1972, AJ, 77, 177
to what we haveshown here. Barnes R., Greenberg R., Quinn T. R., McArthur B. E.,
The results of this work are arguably entertaining, but Benedict G. F., 2011, ApJ,726, 71
are they observable? The short timescales of active bounc- Barrow-Green J.,1997,Poincareandthethreebodyprob-
ing that we observe would suggest not; the consequences of lem
planetsbouncingbackandforth arelikelytobeseenin the Bekov A. A., 1988, Sov. Ast., 32, 106
extensions discussed above, particularly a dissipative envi- Bekov A. A., 1991, Sov. Ast., 35, 103
ronment around the companion, which could lead to per- Burrau C., 1913, Astronomische Nachrichten, 195, 113
manent capture of the scattered planet, or the introduc- BurrowsA.,MarleyM.,HubbardW.B.,LunineJ.I.,Guil-
tion of a planetary system around both stars. The cursory lot T., Saumon D.,Freedman R., SudarskyD., SharpC.,
examination of the response of a planet around the com- 1997, ApJ, 491, 856
(cid:13)c XXXXRAS,MNRAS000,1–11