Modelling Circumbinary Planets: The case of Kepler-38

1.2 MB·

by  Wilhelm Kley

Astronomy&Astrophysicsmanuscriptno.kep38 (cid:13)c ESO2014
January30,2014
Modelling Circumbinary Planets: The case of Kepler-38
WilhelmKley1 andNaderHaghighipour1,2
1 Institutfu¨rAstronomieundAstrophysik,Universita¨tTu¨bingen,AufderMorgenstelle10,D-72076Tu¨bingen,Germany.
e-mail:[email protected]
2 InstituteforAstronomyandNASAAstrobiologyInstitute,UniversityofHawaii-Manoa,Honolulu,HI96825,USA.
e-mail:[email protected]
Received;accepted
ABSTRACT
4 Context.Recently,anumberofplanetsorbitingbinarystarshavebeendiscoveredbytheKeplerspacetelescope.Inafewsystems
1 theplanetsresideclosetothedynamicalstabilitylimit.Duetothedifficultyofformingplanetsinsuchcloseorbits,itisbelievedthat
0 theyhaveformedfurtheroutinthediskandmigratedtotheirpresentlocations.
2 Aims.Ourgoalistoconstructmorerealisticmodelsofplanetmigrationincircumbinarydisks,andtodeterminethefinalpositionof
theseplanetsmoreaccurately.Inourwork,wefocusonthesystemKepler-38wheretheplanetisclosetothestabilitylimit.
n
Methods.Theevolutionofthecircumbinarydiskisstudiedusingtwo-dimensionalhydrodynamicalsimulations.Westudylocally
a
isothermal disks as well as more realistic models that include full viscous heating, radiative cooling from the disk surfaces, and
J radiativediffusioninthediskmidplane.Afterthediskhasbeenbroughtintoaquasi-equilibriumstate,a115Earth-massplanetis

9 embeddedanditsevolutionisfollowed.
2 Results.Inallcasestheplanetsstopinwardmigrationneartheinneredgeofthedisk.Inisothermaldiskswithatypicaldiskscale
heightofH/r = 0.05,thefinaloutcomeagreesverywellwiththeobservedlocationofplanetKepler-38b.Fortheradiativemodels,

] thediskthicknessandlocationoftheinneredgeisdeterminedbythemassinthesystem.Forsurfacedensitiesintheorderof3000
P g/cm2at1AU,theinnergapliesclosetothebinaryandplanetsstopintheregionbetweenthe5:1and4:1mean-motionresonances
E
withthebinary.Amodelwithadiskwithapproximatelyaquarterofthemassyieldsafinalpositionveryclosetotheobservedone.
. Conclusions. For planets migrating in circumbinary disks, the final position is dictated by the structure of the disk. Knowing the
h
observedorbitsofcircumbinaryplanets,radiativedisksimulationswithembeddedplanetscanprovideimportantinformationonthe
p
physicalstateofthesystemduringthefinalstagesofitsevolution.
-
o
Keywords.circumbinarydisks–hydrodynamics–planetformation
r
t
s
a
[ 1. Introduction quence of sticking collisions with subsequent gas accretion for

themoremassiveobjects(Armitage2010).Inabinarystarsys-
1
Aninterestingdevelopmentinexoplanetarysciencehasbeenthe tem,thepresenceofthebinaryinthemiddleoftheprotoplane-
v
8 recent discovery of circumbinary planets (CBPs) by the Kepler tarydiskwillalterthisprocessandmakeplanetformationmore
4 space telescope. In contrast to systems where a planet revolves complicated. The perturbation of the binary and its (eccentric)
6 around onestar of abinary, herethe planet orbitsthe entire bi- diskwilldynamicallyexcitetheorbitsofplanetesimalsandhin-
7 narysystem.Currentlyknownmainsequencebinarieswithcir- dertheirgrowthtolargerobjects(Scholletal.2007;Meschiari
. cumbinary planets are: Kepler-16 (Doyle et al. 2011), Kepler- 2012;Marzarietal.2013).Becauseinacircumbinarydisk,the
1
0 34 and 35 (Welsh et al. 2012), Kepler-38 (Orosz et al. 2012a), outer edge of the central cavity, that is created due to the ef-
4 Kepler-47 (Orosz et al. 2012b) and Kepler-64 (Schwamb et al. fect of tidal forces from the binary on the disk material, coin-
1 2013). In all these systems, the binary is close with an orbital cides with the boundary of the planetary stability, the previous
: periodof7to41days.Theorbitalperiodsoftheirplanetsrange statementimpliesthattheinsituformationofCBPsclosetothe
v
from50to300days. stability limit may not be possible (Paardekooper et al. 2012).
i
X Since the discovery of the first circumbinary planet in the However, because of the small separations of CBP-hosting bi-
r Kepler 16 system, it has been noted that in some of these sys- naries (each of the above-mentioned systems fits into the orbit
a tems, the innermost planet is very close to the boundary of the of Mercury), their effects on the formation of planets at large
dynamical stability (Dvorak 1986; Holman & Wiegert 1999). distances will be negligibly small suggesting that these planets
Given that, as indicated by the observations, the orbital planes couldhaveformedfartheroutandmigratedtotheircurrentorbits
of these binaries are perfectly aligned with those of their plan- (Dunhill&Alexander2013;Marzarietal.2013).
ets,thisimpliesthattheplanetsareformedinflat,circumbinary Planetmigrationisanaturalconsequenceofplanet-diskin-
disks.TheproximityoftheorbitsoftheseCBPstothestability teraction (Nelson et al. 2000; Kley & Nelson 2012). To study
limitthenraisesthequestionastowhethertheseplanetsformed the migration of planets in circumbinary disks, the structure of
intheircurrentorbits,oratfartherdistancesfromthebinaryand the disk has to be analyzed and compared to those around sin-
migratedtotheirpresentlocations. glestars.Themostimportantdynamicaldifferencebetweenthe
Althoughnotallaspectsofplanetformationarefullyunder- formerandlatterdisksistheexistenceofacentralcavityinthe
stood,itiswidelyacceptedthatplanetsformaroundsinglestars circumbinarydisks.AsshownbyArtymowicz&Lubow(1994),
throughabottom-upprocesswheregrowthisachievedviaase- the radial extent of this cavity is a function of the binary semi-
1
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
major axis, eccentricity, and mass-ratio as well as the disk vis- 2. Thehydrodynamicmodel
cosity. As indicated by these authors, for typical values of the
As mentioned above, we consider Kepler-38 to be our test bi-
diskviscosity,anddependingonthebinaryeccentricity,thelo-
nary star and carry out simulations for that system. To model
cationoftheinneredgeofthediskwillbeabout2to3timesthe
the evolution of a disk around a binary star, we consider a cir-
separationofthebinary.
cumbinary disk around Kepler-38, and assume that the disk is
The migration of planets in circumbinary disks was first vertically thin and the system is co-planar. The assumption of
studied by Nelson (2003) for massive, Jupiter-type planets. He co-planarityiswelljustifiedasthemutualinclinationofKepler-
found that for small values of the binary eccentricity, the mi- 38 and its planet is smaller than 0.2 degrees. We then perform
gratingplanetwillbecapturedina4:1mean-motionresonance two-dimensional(2D)hydrodynamicalsimulationsintheplane
(MMR) with the binary whereas in more eccentric systems ofthebinary.
(ebin (cid:38)0.2),theplanetwouldbecapturedinastableorbitfurther To simulate the evolution of the disk, the viscous hydrody-
out (see also Pierens & Nelson 2008b). Subsequent studies by namic equations are solved in a polar coordinate system (r,φ)
Pierens&Nelson(2007)extendedtheseanalysestoplanetswith with its origin at the center of mass of the binary. In our stan-
massesassmallas20 MEarth showingthattheseplanetsusually dard model, which we use as a basis for comparing the results
stop near the edge of the cavity, and they suggested that plan- of our subsequent simulations, we consider the radial extent of
etsincircumbinarydisksshouldbepredominantlyfoundinthat thedisk(r)tobefrom0.25to2.0AU,andφtovaryinanentire
area. As the inner edge of the disk roughly coincides with the annulus of [0,360◦]. This domain is covered by an equidistant
boundaryofplanetarystability,thepredictionsoftheseauthors grid of 256×512 gridcells. For testing purposes, we have also
turned out to be in a very good agreement with the orbital ar- usedagridtwiceasfine.Inallmodelsweevolvethevertically
chitectureofseveralKeplerCBPs.Moregeneralcasesinwhich integrated equations for the surface density Σ, and the velocity
accretion and multiple planets were also considered were later components(v ,v ).
r φ
studiedbythesameauthors(Pierens&Nelson2008a,b). Whensimulatinglocallyisothermaldisks,wedonotevolve
the energy equation, and instead use an isothermal equation of
In a recent paper Pierens & Nelson (2013) revisited planet
stateforthepressure.Whenconsideringradiativemodels,aver-
migration in circumbinary disks and developed models to ex-
ticallyaveragedenergyequationisused,whichevolvesthetem-
plaintheorbitsoftheplanetsaroundbinariesKepler-16,Kepler-
perature of midplane (Mu¨ller & Kley 2012). Radiative effects
34 and Kepler-35. Similar to the majority of the simulations of
are include in two ways. First, a cooling term is considered to
thiskind,theseauthorsconsideredalocallyisothermalapproxi-
account for the radiative loss from the disk surface (D’Angelo
mationwherethediskthermodynamicsismodeledbyprescrib-
et al. 2003). Second, we include diffusive radiative transport in
ing a given temperature profile. They also considered a closed
the midplane of the disk using flux-limited diffusion (Kley &
boundaryconditionattheedgeofthecavityassumingzeromass
Crida 2008). The use of a flux-limiter is required because near
accretionontothecentralbinary.AlthoughtheworkbyPierens
itsinneredge,thediskisveryopticallythin,whichwouldleadto
&Nelson(2013)representssignificantresults,theirapplicability
anunphysicalhighenergyfluxinthediffusivepartofradiative
maybelimitedduetothesimplifyingassumptionsofanisother-
transport equations. For more details on the implementation of
maldiskwithaclosedboundaryconditionsforthecentralcav-
thisdiffusivetransport,seeMu¨ller&Kley(2013).Inourradia-
ity.NumericalsimulationsbyArtymowicz&Lubow(1996)and
tivesimulations,weconsiderthefull,verticallyaverageddissi-
Gu¨nther&Kley(2002)haveshownthatdespitetheappearance
pation(Mu¨ller&Kley2012).However,stellarirradiationisnot
of the above-mentioned cavity in the center of a circumbinary
takenintoaccount.
disk, material can still flow inside and onto the central binary.
TocalculatethenecessaryheightofthediskH ataposition
Inotherwords,amorerealisticboundaryconditionwouldallow
r,wefirstnotethatinacircumstellardisk(aroundasinglestar
forthein-flowofthematerialthroughtheinneredgeofthedisk
atthepositionr ),theverticalheightofthediskisgivenby
intothediskcavity. ∗
Here,wepresentanimprovedandextendeddiskmodelthat (cid:32)|r−r |3(cid:33)−21
allowsformuchmorerealisticsimulations.Weconsiderthesys- H(r)=c ∗ . (1)
s GM
tem of Kepler-38 as this system represents one of the binaries ∗
in which the circumbinary planet is in close proximity to the
Inthisequation,M isthemassofthecentralstar,c isthespeed
stabilitylimit.Wehavedevelopedanimprovedapproachwhich ∗ s
ofsoundinthediskatthelocation r,andG isthegravitational
includesdetailedbalanceofviscousheatingandradiativecool-
constant. In a binary star system, H(r) has to be calculated by
ingfromthesurfaceofthedisk(D’Angeloetal.2003),aswell
as additional radiative diffusion in the plane of the disk (Kley taking the contributions of both stars into account (Gu¨nther &
Kley2002).Thatis,
& Crida 2008). For a planet embedded in disks, this improved
thermodynamicscanhavedramaticeffectsontheplanetorbital
 −1  −1
ddyirneacmtioicnsosfucmhigthraattifoonrl(oKwle-yma&ssCprliadnaet2s0,0it8c;aKnleevyenetreavl.er2s0e0t9h)e. H(r)=(cid:88) c2√G|rM−ir|3 2 =(cid:88)Hi−2(r) 2 , (2)
Also,unlikePierens&Nelson(2013),weallowfreein-flowof i=1,2 s i i=1,2
materialfromthediskintothecentralcavityandconstructmod-
whereM = M andM arethemassesofthestarsofthebinary.
elswithnetmassin-flowthroughthedisk. i 1 2
Equation2indicatesthatinacircumbinarydisk,thetotalheight
Thispaperisorganizedasfollows.InSection2,wepresent HisalwayssmallerthantheindividualheightsH.
i
thephysicalsetupofourdiskmodels.InSection3,westudylo- The equation of state of the gas in the disk is given by the
cally isothermal disks and in Section 4, we present the results ideal gas law using a mean molecular weight, µ = 2.35 (in
of our full radiative models and compare them to the results atomic mass units), and an adiabatic exponent of γ = 1.4. For
from standard isothermal cases. Section 5 concludes this study theshearviscosityweusetheα-parametrizationwithα = 0.01
bysummarizingtheresultsanddiscussingtheirimplications. forourstandardmodel,andwesetthebulkviscositytozero.For
2
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
Table 1. The binary parameter and the observed planetary pa- density at 1 AU. Such a radial profile for Σ corresponds to that
rameter of the Kepler-38 system used in the simulations. The usuallyadoptedfortheminimumsolarmassnebula.Asthetotal
massoftheprimarystaris0.949M .Thevalueshavebeentaken stellarmassintheKepler-38systemisapproximately1.2solar-
(cid:12)
fromOroszetal.(2012a). masses, one can assume that, to a first order of approximation,
the mass of the circumbinary disk would lie in the same range
a)BinaryParameter as that of the disk around the protosun. Assuming that the ob-
Massratio Period abin ebin servedcircumbinaryplanetshavereachedtheircurrentpositions
q=M2/M1 [days] [AU] in the late phase of their evolution, such a disk surface density
0.2626 18.6 0.1469 0.1032
may be on the high side. However, together with the relatively
large value of the viscosity parameter (α = 0.01), it allows for
b)PlanetParameter asufficientlyfastevolutionofthesysteminourstandardmodel
Mass Period a e
p p such that it can be numerically simulated. We will vary these
M [days] [AU]
Jup parametersfortheradiativemodelsbelow.
0.34 105.6 0.46 <0.03
We choose the initial temperature of the disk to vary with
r as T(r) ∝ r−1 such that, assuming a central star of M , the
bin
verticalthicknessofthediskwillalwaysmaintainthecondition
the Rosseland opacity, we use analytic formula as provided by
H/r =0.05.Theinitialangularvelocityofthediskatadistance
Lin&Papaloizou(1985),andtheflux-limiterasinKley(1989).
rischosentobeequaltotheKeplerianvelocityatthatdistance,
To calculate the gravitational potential of the binary-planet
andtheradialvelocityissettozero.
system at a position r in the disk, we use, for all 3 bodies, a
Theparametersofthecentralbinaryhavebeenadoptedfrom
smoothedpotentialfunctionoftheform
Oroszetal.(2012a)andcanbefoundinTable1,togetherwith
Ψ (r)=− GMk , (3) theplanetarydata.Duetotheinteractionofthebinarywiththe
k [(r−r )2+((cid:15)H)2]1/2 disk (and planet), these parameters will slowly vary during a
k
simulation. For the models with an embedded planet, we con-
where M denotes the masses of the objects (stars and planet),
k sider the planet to have a mass of m = 3.63×10−4M which
r−r isthevectorfromapointinthedisktoeachobject,andHis p 1
k isequivalenttom = 0.34M orabout115Earth-masses(see
calculatedusingeq.(2).Thequantity(cid:15) ineq.(3)isasmoothing p Jup
Table1).Atthebeginningofeachsimulation,westartthebinary
parameterthataccountsfortheeffectofthefinitethicknessofthe
atitsperiastron.
disk.Inoursimulations,weassume(cid:15) =0.6or0.7(Mu¨lleretal.
Theboundaryconditionsofthesimulationsareconstructed
2012).Wenotethatfortheplanet,thesmoothinglengthcannot
such that at the outer boundary of the disk, r = 2.0 AU, the
besmallerthan(cid:15)r ,wherer =a (m /3M )1/3istheradiusof max
H H p p bin surface density remains constant. This is achieved by using a
theHillspherearoundtheplanet,a isthesemi-majoraxisofthe
p dampingboundaryconditionwherethedensityisrelaxedtoward
planetsorbit,and M = M +M .Theindirecttermsaretaken
bin 1 2 itsinitialvalue,Σ(r ) = 2121g/cm2,andtheradialvelocityis
intoaccountthroughashiftofthepositionsoftheobjectssuch max
damped toward zero. For this, we use the procedure specified
thatthebinary’sbarycenteralwayscoincideswiththeoriginof
bydeVal-Borroetal.(2006).Theangularvelocityattheouter
thecoordinatesystem.
boundary is also kept at the initial Keplerian value, and for the
Thehydrodynamicalequationsareintegratedusingamixed,
temperatureweuseareflectingconditionsothattherewillbeno
explicit/implicit, scheme. The standard hydrodynamical terms
artificialradiativefluxthroughtheouterboundary.Thesecondi-
areintegratedexplicitlyusingastandardCourantconditionwith
tionsattheouterboundaryleadtoadiskwithzeroeccentricityat
Courantnumber f = 0.66.Theviscousandthediffusiveradia-
C r .Hence,r hastochosenlargeenoughsuchthattheinner
tivepartsareintegratedimplicitlytoavoidinstabilities.Forthe max max
regionsarenotinfluenced.
time-integrationofthehydrodynamicalequations,weutilizethe
At the inner boundary (r ), we consider a boundary con-
RH2D-codewiththeFARGO(Masset2000)upgrade. min
dition such that the in-flow of disk material onto the binary is
Themotionsofthestellarbinaryandtheplanetareintegrated
allowed. This means, for the radial boundary grid cells at r ,
usinga4thorderRunge-Kuttaintegrator.Allobjects,including min
we choose a zero-gradient mass out-flow condition, where the
thetwostars,feelthegravityofthediskaswellastheirmutual
material can freely leave the grid and flow onto the binary. No
gravitation.Theself-gravityofthediskisnotconsidered.
mass in-flow into a grid is allowed at r . The zero-gradient
Tocalculatetheforcefromthediskactingontheplanet,we min
condition is also applied to the angular velocity of the material
exclude parts of the Hill sphere of the planet using a tapering
since due to the effect of the binary, no well-defined Keplerian
function
(cid:34) (cid:32) d/r −p(cid:33) (cid:35)−1 velocitycanbefoundwhichcouldbeusedotherwise.Thiszero-
f(d)= exp − H +1 . (4) gradient condition for the angular velocity implies a physically
p/10
morerealisticzero-torqueboundary.
Inthisequation,disthedistancefromthegridcelltotheplanet, Withtheseboundaryconditions,thediskcanreachaquasi-
and p is a dimensionless parameter that is set to 0.6 for the stationary state in which there will be a constant mass-flow
isothermal models and 0.8 for the radiative runs. We calculate throughthedisk.
theorbitalparametersoftheplanetusingJacobiancoordinates,
assumingtheplanetorbitsastarwithmass M = M +M ,at
bin 1 2 3. Thelocallyisothermalcase
thebinarybarycenter.
Inthissectionwedonotevolvetheenergyequationbutinstead
leavethetemperatureofthediskatitsinitialvalue.Thisproce-
2.1. Initialsetupandboundaryconditions
dure has the advantage of making the simulations much faster
In all models, the disk initial surface density is chosen to have asnoheatingandcoolingofthediskhastobeconsidered.This
a Σ(r) = Σ r−1/2 profile where r is the distance from the cen- isanapproximationthathasoftenbeenusedinplanetdisksim-
0
ter of mass of the binary, and Σ = 3000g/cm2 is the surface ulations and gives the first order results, but it also has its lim-
0
3
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
Table2.Propertiesforthelocallyisothermaldiskmodelswith- 1.5e-05
outaplanet.Model1isourreferencemodelthatisdescribedin r = 0.31
r = 0.68
detail in Sect. 2.1. In the other two models certain (numerical) 1e-05 Mean
parametervariationsthathavebeenappliedinordertocompare
tTohepyreavrieoudsesmcroibdeedls,ininthpeaqrtuiocuteldarsbecytiPoinesr.ens & Nelson (2013). M/yr]sun 5e-06
e [
at 0
Name Characterization Section n R
o
MMooddeell12 sctlaonsdedaridnnmeordbeolundary i3n.22.1 ccreti -5e-06
Model3 differentinnerradii 3.3 A
-1e-05
0.0038
-1.5e-05
3950 3950.5 3951 3951.5 3952
0.0037
Time [Years]
0.0036
ass [M]sol 0.0035 Fm−i8og.d.6e28l.×1Gr1ma0pe−ha7s,ouifrsetdahlesaotratptwelooottfdeidmffefaorsersnratcd=cirset0tai.no6cn8estAh.UrTo.huAeghsmttehhaeendcviesanklturiaenl,
M 0.0034
k binary accretes material, the mass flow is directed inward and
s
al Di 0.0033 becomesnegative.However,inthemaintext,wereporttheac-
ot cretionrateinformofpositiveabsolutevalues.
T 0.0032
0.0031
6000
0.003
tt == 00
0 500 1000 1500 2000 2500 3000 3500 4000 4500
tt == 330000
Time [Years] 5000 tt == 664400
Fig.1. Graph of the total mass of the disk (without the planet) 2cm] 4000 tttt ==== 1111282850500000
inourlocallyisothermalstandardmodel(model1).Afterafew y [g/ tt == 33775500
thousand years, the disk reaches equilibrium and a nearly con- nsit 3000
e
stantmassisestablished. e d
c
a 2000
urf
S
itations when planetary migration is concerned (Kley & Crida
1000
2008).Wewillfirstshowtheresultsofourstandardcase(model
1), and then compare these results to those recently obtained
0
byPierens&Nelson(2013)whoalsousedalocallyisothermal
model.Table2givesaverybriefoverviewofthemodelsandcan
0.25
serveasareference.
We begin by discussing models without a planet where we
y 0.2
bring the disk into equilibrium and analyse its structure. This cit
equilibration process is required because the disk structure is ntri
ce 0.15
determinedbytheactionofthebinaryforwhichnosimpleana- Ec
lyticmodelcanbeprescribed.Additionally,theevolutionofthe sk
Di 0.1
planetinthediskoccurssoslowlythatthediskalwaysremains
nearequilibrium.
0.05
3.1. Thestructureofthecircumbinarydisk 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Toanalysethedisk’sowndynamics(i.e.,withoutaplanet)due Radius [AU]
totheeffectofthecentralbinary,wesimulatedthedynamicsof
thediskinourstandardmodelwithazero-massplanet.Dueto Fig.3.Azimuthallyaveragedsurfacedensityanddiskeccentric-
the fact that the disk inner boundary is open, and at the outer ityformodel1.Displayedaretheprofilesatdifferentevolution-
boundary, its surface density is considered to be constant, the arytimes(inyrs),wheretheredcurvedenotestheinitialsetup.
disksettlesintoafinalquasi-stationarystateinwhichthemass
ofthediskisconstantandthereisaconstantmassaccretionrate
ontothecentralbinary.ThishasbeenshowninFig.1wherethe Thesmallremainingoscillationsarecausedbytheeccentricmo-
totalmassinthecomputationaldomainisplottedvs.timeforthe tion of the inner binary which induces variations of the disk
locallyisothermalmodelwithnoplanet.Afterabout2000-3000 structure. In this simulation where we have been interested in
yrs a quasi-stationary state has been reached at which state the the structure of the disk without the planet, the disk does not
outflow through the inner boundary is exactly balanced by the affect the motion of the binary. Nevertheless, the motion of the
’inflow’fromoutsideasestablishedbykeepingthevalueofthe binaryhasbeennumericallyintegratedtoexaminetheaccuracy
densityconstantatr .Duringthisprocess,thediskmasshas ofthenumericalmethod.Wefindthatduringthetotalevolution
max
decreasedbyabout20%from3.8×10−3M to3.04×10−3M . timeofabout4000years,whichcorrespondstoover80,000bi-
(cid:12) (cid:12)
4
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
naryorbitsand4.5million(!)timesteps,thesemi-majoraxisof
thebinaryhasonlyshrunkby0.014%,avaluesufficientlysmall
forourpurposes.
At equilibrium, the rate of the mass-flow through the disk,
M˙ , becomes on average nearly constant. Fig. 2 shows this
disk
where M˙ ,measuredattwodifferentdistances,isplottednear
disk
the end of the simulation. The resulting accretion rate equals
approximately M˙ = 8.7 × 10−7M /yr. The oscillations are
disk (cid:12)
much larger in the inner parts of the disk due to the perturba-
tion from the binary, and they occur at the binary’s orbital pe-
riod.Thisvaluemaybecomparedtothetypicaltheoreticalequi-
librium disk accretion rate, M˙ = 3πΣν. At the outer bound-
th
ary, where the mass density has been kept constant, we find
M˙ =5.4×10−7M /yr.
th (cid:12)
During time, the surface density slowly evolves away from
its initial profile until it settles in a new equilibrium state. This
isshownintheupperpanelofFig.3wheretheazimuthallyav-
eragedradialsurfacedensityprofileisshownatdifferentevolu-
Fig.4.Twodimensionaldensitystructureofanisothermaldisk
tionary times. In agreement with Fig. 1, the equilibration time
(model1)aroundthecentralbinarystar.Thegraphshowsalocal
takes about 2000 yrs. At this state, the surface density profiles
view around the binary. The computational grid extends from
will no longer change with time. Because of the tidal action of
r = 0.25 AU to r = 2.0 AU. The white inner region lies inside
the binary on the disk, a central gap is formed with a surface
the computational domain and is not covered by the grid. The
density many orders of magnitude smaller than inside the disk.
positionsoftheprimaryandsecondaryareindicatedbythegray
Theinneredgeofthedisk,whichwecandefine,veryapproxi-
dots.ThetheRoche-lobeofthesecondaryisalsoshown.
mately,asthatradiusatwhichthesurfacedensityisabouthalf
the maximum value, lies here at around r ≈ 0.45AU. This is
slightlylargerthan3a andinagoodagreementwiththefind-
bin
ings of Artymowicz & Lubow (1994). We note that for the bi- simulationsusingtheirboundaryconditions.InFig.5,weshow
narysystemconsideredhere,thestabilitylimitforplanetaryor- thesurfacedensityanddiskeccentricityfortwosimulations;our
bits is about 0.4 AU (Dvorak 1986; Holman & Wiegert 1999). standardmodelwithanopeninnerboundary,anddiskmodel2
Outsideofthegapbeyond r = 0.6AU,thesurfacedensitypro- (seeTable1)withreflectingconditionsatrminasusedbyPierens
file becomes relatively flat. One can notice that the mass-flow &Nelson(2013).Becauseinthelattermodel,materialisnotal-
ratesshowninFig.2correspondtoaregiondeepinsidethegap lowedtoleavethediskthroughtheinnerboundary,itwillaccu-
(r=0.31≈2a ),andimmediatelyoutsideofit(r=0.68). mulateneartheouteredgeofthegap,asseenbythespikeinthe
bin
As was shown by Pierens & Nelson (2013), circumbinary densityatr=0.8AUinthetoppanelofFig.5.Atthesametime,
diskscanattainsignificanteccentricities.Wecalculatetheeccen- theeccentricityofthediskinthismodelbecomessubstantially
tricityofthediskbytreatingeachgridcellasanindividualpar- largerthanthatinthemodelwiththeopenboundarycondition.
ticlewithamassandvelocityequaltothemassandvelocityof Our results of the simulations with the closed inner boundary
thecell(Kley&Dirksen2006).Tocalculatearadialdependence are very similar to those reported by Pierens & Nelson (2013),
forthediskeccentricity,e (r),weaverageovertheangulardi- although these authors did not model the system of Kepler-38,
disk
rection.Inoursimulations,thediskeccentricityremainssmall, butsomerelatedones.
asshowninthelowerpanelofFig.3.Intheregionsoutsideof
thecentralgap(beyondr = 0.6),thediskeccentricityisalways
3.3. Varyingthelocationoftheinnerradius
belowe < 0.07.Atradialdistancesr > 1.0,thiseccentricity
disk
becomes smaller than about 0.01. These values of the disk ec- In this section, we examine the effect of the location of the in-
centricityareincontrasttothosereportedbyPierens&Nelson nerboundaryonthediskstructure(model3).Pierens&Nelson
(2013), who found, on average, larger values. We attribute this (2013) found that, at least in the case of a closed inner bound-
differencetotheassumptionofanopeninnerboundary,andwe ary,theamplitudeofthediskeccentricitydependsontheloca-
willanalysethisinmoredetailinthenextsection. tion of r . In particular, if r is chosen so large that some
min min
In Fig. 4, we show a the two-dimensional surface density MMRs between the binary and the disk fall outside the com-
distributionforourisothermaldiskmodels.Notethatthefigure putationaldomain,theeccentricityofthediskwillremainvery
shows only the inner part of the computational domain around small.Toexaminewhetherthediskwillbehavesimilarlyforan
thecentralbinary.TheRochelobeofthesecondarystarisshown openinnerboundary,wecarriedouttwoadditionalsimulations
aswell.Asshownhere,aneccentriccentralbinarystronglyper- with slightly smaller inner radii, r = 0.22 AU and 0.20 AU.
min
turbsthediskandproducestimevaryingpatterns.Asindicated Giventhebinary’ssemi-major,a =0.1469AU,theseradiicor-
bin
bytheflatsurfacedensityprofile,outsideofthecentralgap,the respondto1.5and1.36a ,respectively.
bin
diskshowslessstructureandisrelativelyhomogeneous. AsshowninFig.6(toppanel),forsmallvaluesofr ,the
min
densityofthediskincreasesoutsideofthegap,andatthesame
time, the flow of the mass onto the binary is slightly reduced
3.2. Aclosedinnerboundary
(to about 8 × 10−7M /yr). The profile of the edge of the gap,
(cid:12)
Inthissection,westudytheeffectofdifferentboundarycondi- where the slope of the density is positive, is hardly affected by
tionsattheinneredgeofthecomputationaldomainonthefinal the choice of r . However, the gap is slightly wider. The ec-
min
results. Inspired by the work of Pierens & Nelson (2013), who centricityinthemainpartofthedisk(outsider = 0.7)doesnot
focusedonmodelswithclosedinnerboundaries,wecarriedout vary with the location of r (lower panel of Fig. 6). At r = 1
min
5
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
3000
6000
OOppeenn
CClloosseedd
2500
5000
2m] 2m]
g/c 4000 g/c 2000
y [ y [
ensit 3000 ensit 1500
d d
e e
c c
urfa 2000 urfa 1000
S S 00..2255 AAUU
1000 500 00..2222 AAUU
00..2200 AAUU
0 0.4 05
0.3
0.4
0.25 0.35
y y 0.3
cit 0.2 cit
ntri ntri 0.25
e e
Ecc 0.15 Ecc 0.2
k k
s s
Di 0.1 Di 0.15
0.1
0.05
0.05
0 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Radius [AU] Radius [AU]
Fig.5.Azimuthallyaveragedsurfacedensityanddiskeccentric- Fig.6. Azimuthally averaged surface density and eccentricity
ityprofilesfortwomodelswithdifferentinnerboundarycondi- profilesforthreediskmodelswithdifferentlocationsofthein-
tions.Theredcurvecorrespondstothestandardmodelwithan ner boundary, r . The red curve shows the standard model
rmin
openinnerboundaryandthebluecurverepresentsasystemwith (model1)withr =0.25AU,thebluecurveisforadiskwith
min
areflectinginnerboundary(model2).Inbothmodels,thedisk r = 0.22 AU, and the green curve represents the disk model
min
isatequilibrium. withr =0.20AU(model3).Inallthreemodels,thediskisat
min
equilibrium.
thediskeccentricityisaboute =0.1andbeyondr =1.0itis
disk
alwaysbelow0.01.Withinthegapregion,wherethedensityis
(A.Pierens,privatecommunication).Foroursystem,Kepler-38,
small,e risestoamaximumof0.4attheinnerboundaryfor
themoddeislkswiththesmallerr .Thediskswithr = 0.2AU wefoundtheoppositebehaviour;areduceddiskeccentricityfor
min min
andr =0.22AUshowverysimilarbehaviour.Forthesakeof anopeninnerboundary.Thisfollowsfromthefactthatanopen
min
innerboundaryreducesthestrengthofaglobaleccentricmode.
comparison,wealsoransimulationswithhigherspatialresolu-
Becauseinthispaperwearemoreinterestedintheevolutionof
tion.Theresultswereverysimilartothosedepictedbythetwo
planets in circumbinary disks, we do not study this interesting
panelsoffigure6.
behaviourofthediskanyfurtheranddeferthattoafuturetime.
Overall, results of our simulations indicate that in a model
We note that among the three models shown in figure 6,
withanopeninnerboundary,thediskeccentricitye remains
disk
the model with the smaller inner radii can have a considerably
small in a large portion of the disk and independent of the
smaller numerical time step. In addition to the smaller size of
location of the inner boundary. This is unlike the findings of
thegridcells,thiseffectiscausedbythelargerdeviationfroma
Pierens & Nelson (2013) who showed that in cases with close
uniformkeplerianvelocityasinducedbythebinarystar,which
boundaries, the disk eccentricity will rise to much larger val-
makes the FARGO-algorithm less effective. Therefore, to save
ues(Pierens&Nelson2013).Thelattercanbeattributedtothe
computationaltimeandresources,intheremainderofthisstudy,
factthatareflectinginnerboundarycreatesacavitywhichsus-
wewillusefortheinnerradiusofthedisk,astandardvalueof
tainstheeccentricglobalmodeofthedisk.Ourfindingsarealso
r =0.25AU,althoughwewillpresentsomecomparisonsim-
inagreementwithKleyetal.(2008)whofoundaneccentricity min
ulationsusingr =0.22AUaswell.
reduction for circumstellar disks in close binary systems when min
using outflow boundary conditions. In our opinion, open inner
boundaries are more realistic because they allow for the accre-
3.4. Planetsinisothermaldisks
tion of material onto the binary stars. However, we would like
tocautionthatinthecaseofmoreeccentricbinaries,ithasbeen To study the evolution of the system with an embedded planet,
foundthatanopeninnerboundarymaycausethecentralcavity we start our simulations with a planet initially at different dis-
to be very large since disk material cannot re-enter the compu- tances(semi-majoraxes,a )fromthecenterofmass(barycen-
0
tational domain once it has been lost through r (see Marzari ter)ofthebinaryandinacircularorbit.Theplanet’sevolutionis
min
et al. (2009)). This situation occurs primarily in highly eccen- determined by the gravitational action of the disk and the cen-
tricbinariessuchasKepler-34whichhasaneccentricityof0.52 tral binary. As mentioned earlier, during the evolution of the
6
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
1.8 the disk and the planet evolve simultaneously. In the other two
initial disk: a0=1.75 simulations, the planet is started in a disk that is already in an
1.6 rreellaaxxeedd ddiisskk:: aa00==11..7050 equilibrium. The blue line corresponds to the planet starting at
0.47 a = 1.75 AU and the green line is when the planet starts at
0
AU] 1.4 0.46 a=1.0AU.Wehaveshiftedthelastmodelintimetoensurean
major axis [ 1.2 0.45 overAlaspswhoitwhnthienfitrhsettFwiog.. 7, the semi-major axis of the planet
Planet Semi- 0 .18 000...444234 cbToihnneatirsnytu)roounungstliydlrsiothprrieniankcshth(eies.edt.h,isethkeinspnulerafrnaeccteavmditieygnrscairtteeysatwteodiwthbairyndttthhheeebcceeinnnattrrraayll.
6000 6500 7000 cavityandtheresultingpositivedensitygradientactasaplanet
0.6 trap(Massetetal.2006).HerethenegativeLindbladtorquesare
balanced by positive corotation torques. As a result, the planet
0.4 stopsitsinwardmigrationinthisregion.Asshowninthefigure,
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Time [Years] in all runs, the planet stops at approximately the same distance
from the binary. The semi-major axis of the planet in the three
Fig.7. The evolution of the semi-major axis of the planet in runs has an average value of a = 0.436 AU (blue line in inset
an isothermal disk (model 1). In two simulations, the planet is p
ofFig.7).Theresultsindicatethatinallsimulations,thesystem
started in a circular orbit at a distance of a = 1.75 AU from
0 evolvestothesamefinalstate.Atthisstate,theplanetresidesin
thecenterofmassofthecentralbinary.Inthesimulationshown
anorbitwithaperiodofabout100dayswhichisslightlyoutside
in red, the planet is included in the disk from the very begin-
ofitsobservedorbit(seeTab.1).Thisimpliesthattheoutcome
ningofthesimulationanditevolvesasthediskevolvesstarting
of the migration process does not depend on the history of the
fromitsinitialdensityprofile.Inthesimulationshowninblue,
system and is determined solely by the physical parameters of
thediskisfirstrelaxedtotheequilibriumandthentheplanetis
the disk. In an isothermal system and for a given binary, these
embedded in it. In the simulation shown in green, the planet is
parametersaredeterminedbythevaluesofH/randviscosity.
started at a = 1.0 AU in a relaxed disk. The inset shows the
0 Fig. 8 shows the evolution of the eccentricity of the planet
results near the end of the simulation for the first model. The
correspondingtotheredmodelinFig.7.Asshownhere,during
dashedblacklineintheinsetcorrespondstothemeanvalueof
theplanetsinwardmigration,itseccentricityincreasescontinu-
theplanetsemi-majoraxisata =0.436AU.
p ously until at t = 4500 when it reaches a constant value close
to0.05andoscillatesaroundthatvalueforapproximately1000
0.25 years. At time (t ∼ 5700), the planet is temporarily captured
inamean-motionresonancewiththebinaryanditseccentricity
undergoesasuddenjump(Murray&Dermott1999)whereitos-
0.2
cillatesbetween0.15and0.20.Wewilldiscussedthisresonant
city captureinmoredetailinthenextsection.
ntri 0.15 Fig. 9 shows four two-dimensional snapshots of the evolu-
e
c tion of the disk surface density during the inward migration of
c
e
net 0.1 an embedded planet in one of our models. One can clearly see
Pla thespiraldisturbancescreatedbytheplanetduringitsmigration.
Thelastpanel,correspondingtot = 7417,showsthefinalposi-
0.05 tionoftheplanetinanorbitwithasemi-majoraxisofa =0.436
p
AU(seeFig.7).InFig.10,weplottheazimuthallyaveragedsur-
0 facedensityofthediskcorrespondingtothesametimesasinthe
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 snapshots of Fig. 9. The colored circles in this figure represent
Time [Years] theplanetinitsactualradialdistancefromthecenterofmassof
thebinary.Aplanetofthismassdoesnotproduceafullgap,but
Fig.8.Graphoftheevolutionoftheeccentricityoftheplanetin
onlyclearsoutaboutathirdofthesurfacedensity.Astheplanet
anisothermaldiskshowninredinFig.7(model1).Theplanet
moves inward, the surface density of the inner part of the disk
wasembeddedinthediskpriortothestartofthediskevolution
atasemi-majorofa =1.75AU. isstronglydistorted.However,itregainsitsunperturbedprofile
0
aftertheplanethasreacheditsfinalposition,(seeFig.3).
planet, its orbital elements are calculated using Jacobian coor-
3.4.1. CaptureinResonance
dinateswithrespecttothebarycenterofthecentralbinary.The
planetmassisfixedto115Earthmasses,andthereisnoaccre- Theevolutionofthesemi-majoraxisoftheplanetandinpartic-
tion onto the planet. As a reference, we present in Tab. 1, the ular,thesuddenincreaseinitseccentricityattimet≈5700seem
orbital parameters of the planet of Kepler-38, as inferred from topointtoatemporarycaptureinamean-motionresonancewith
theobservations. thecentralbinary.Tofurtheranalyzethis,weplottedthevaria-
Fig. 7 shows the evolution of planet through the disk. We tions of the period-ratio of the planet and binary, P /P ,
planet bin
present here the results of three simulations that were carried duringtheinwardmigrationoftheplanet.Fig.11showsthere-
outfordifferentinitialsetupstoexaminethedependenceofthe sults.Asshownhere,astheplanetmigratesinward,thisperiod-
outcome on the initial state of the disk and planets orbital ele- ratio continuously decreases until it reaches a constant value at
ments. In the first model (red line) the planet is started directly t≈5000.TheinsetinFig.11showstheevolutionofthesystem
intheinitialdiskwhichhasnotbeenbroughtintoanequilibrium neartheendofthesimulationforonethousandyears.Theblack
(i.e. using the initial disk setup of model 1). In this simulation, dashedlineinthisinsetdenotesthearithmeticaveragedvalueof
7
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
45
40
5.7
35 5.6
5.5
30
P/Pplanetbin 2205 5555....1234
5
15 7000 7200 7400 7600 7800 8000
10
5
0 1000 2000 3000 4000 5000 6000 7000 8000 9000
Time [Years]
Fig.11.Theevolutionoftheperiodratiooftheplanetandbinary
forthemodelwheretheplanetisembeddedintheinitialdiskat
a =1.75AU.Theinsetshowstheperiodratioforthefinalstage
0
ofitsevolution.Thedashedlineintheinsetcorrespondstothe
meanvalueofP /P =5.22.
planet bin
Fig.9. Two dimensional snapshots of the evolution of the disk
surface density with an embedded planet. The snapshots corre-
spondtothetimes449,729,2807,and7414years.Theplanet,
shownasthegraydotwiththeRocheradius,wasembeddedin
thediskataninitialdistanceofr = 1.75AU.Thedashedgray 150
linearoundthecentralbinaryshowstheapproximateboundary
100
ofstabilityforplanetaryorbits.Amovieofthesimulationshown
herecanbefoundinthesupplementaryonlinematerial. 50
hi1 0
P
3000
-50
2500
-100
2m]
c 2000 -150
g/
y [
nsit 1500
e
d 150
e
c
urfa 1000 100
S
t = 449
500 t = 729 50
t = 2807
0 t = 7414 Phi2 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Radius [AU] -50
Fig.10. Graph of the averaged radial surface density of a disk -100
withanembeddedplanet.Thecolorscorrespondtothetimesof
thesnapshotsshowninFig.9.Thebigcircleineachgraphrep- -150
resentstheradialpositionoftheplanetatthetimecorresponding
5620 5640 5660 5680 5700 5720 5740 5760 5780 5800
tothegraph.Forillustrativepurposes,wehavemovedthecircles
Time [Years]
closetotheircorrespondingcurves.
the period ratio, P /P = 5.22, suggesting a possible cap- Fig.12.GraphsoftheresonantanglesΦ1andΦ2(seeEqs.6and
planet bin 7)duringthecapturephaseina5:1mean-motionresonance.
tureintoa5:1MMRwiththebinary.
To determine whether the planet is truly captured in a 5:1
MMR, its corresponding resonant angles have to be analyzed. anglesaredefinedas
Thecaptureinaresonancebetweentwobodiesoccurswhenat
leastoneoftheresonantanglesofthesystemlibratesarounda Φ = pλ −qλ −p(cid:36) +q(cid:36) +k((cid:36) −(cid:36) ), (5)
k p b p b p b
certain value and does not cover the full range of 0 to 2π (see
Murray & Dermott (1999) for details). For a general p:q com- whereλ ,λ ,(cid:36) and(cid:36) denotethemeanlongitudesandlongi-
b p b p
mensurability with p > q and in a planar system, the resonant tudeofperiapsefortheinnerbinary(b)andouterplanet(p).The
8
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
1 1
Model 1: a0=1.00 rmin = 0.25, closed
U] 0.9 00..5205 DDiisskk MMaassss U] 0.9 rrmmiinn == 00..2252,, ooppeenn
A A
xis [ 0.8 xis [ 0.8
a A
major 0.7 major 0.7
mi- mi- 0.6
net Se 0.6 net se 0.5
a a
Pl Pl
0.5
0.4
0.4 0.3
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 1000 2000 3000 4000 5000 6000 7000
Time [Years] Time [Years]
Fig.13.Theevolutionofthesemi-majoraxisofaplanetembed- Fig.14.Theevolutionofthesemi-majoraxisofaplanetembed-
ded in an isothermal disk for different values of the disk mass. dedinanisothermaldiskfordifferentlocationsofthediskinner
Thegraphinredcorrespondsthesimulationofthestandarddisk boundary. The graph in red corresponds to the standard model
modelwiththeplanetinitiallystartingata =1.0AU.Thegraph withanopeninnerboundaryat0.25AU.Thegraphinbluerep-
0
in blue shows a similar simulation but with half the disk mass. resents a model with a closed inner boundary at 0.25 AU, and
The graph in green corresponds to a simulation with a quarter thegraphingreenis foradiskwithanopeninnerboundary at
oftheinitialdiskmassandcontinuedfromthemodelinblueat 0.22AU.
t=4300.
integer k in this equation satisfies the condition q ≤ k ≤ p and thatdiskswithreducedsurfacedensityshowsimilardensityevo-
has p−q+1possiblevalues.Ofthe p−q+1resonantangles,at lutionasthoseshownbythetoppanelofFig.3,butonlyscaled
mosttwoarelinearlyindependent.Thatmeans,inaactualreso- down by an appropriate factor. This identical density evolution
nantconfiguration,atleastoneoftheseanglewilllibrate(Nelson can be attributed to the fact that in an isothermal disk with an
&Papaloizou2002).Asournumericalsimulationsindicatethe α-type viscosity and no self-gravity, the surface density of the
possibility of a capture in the 5:1 MMR, we studied the time- disk is scaled out of the normalized equations of state making
variationsofallresonantanglesofoursystem.Notethatinthis thediskevolutionindependentoftheactualdensity.
case p = 5, q = 1, and k runs from 1 to 5. Fig. 12 shows the Inadiskwithasmallermassandsurfacedensity,theinward
evolutionofthefollowingtworesonantanglesfor180yrsfrom migration of the planet will be slower. Fig. 13 shows the evo-
t=5620tot=5800, lution of a planet in an isothermal disk with different masses.
The red line in this figure represents the inward migration of a
Φ1 =5λp−λb−4(cid:36)p, (6) planet started at a0 = 1.0 AU in our standard disk model. The
blue line corresponds to the migration of the same planet in a
and
disk with half the mass of our standard model. As shown here
Φ =5λ −λ −3(cid:36) −(cid:36) . (7)
2 p b p b andasexpected,theinwardmigrationisabouttwotimesslower
As shown here, before t = 5670, the evolution of these angles thaninthestandardcase.Westartedourthirdsimulations(green
is rather irregular confirming that the planet was not in a reso- line)withaquarterofthediskmassandtosavecomputational
nance.Afterthistime,Φ showsaregularcirculatingbehaviour, time, continued it from the 2nd mode (blue line) at time 4300.
1
whereastheangleΦ beginstolibratearoundzerowithanam- Theplanetmigrationrateisagainreducedbyafactoroftwo.As
2
plitude of less than 40 degrees. This librating evolution of Φ pointedout,forinstance,byPierens&Nelson(2013)whenan-
2
indicates that the planet is indeed captured in a 5:1 MMR with alyzingthesystemofKepler-16,inasystemwithaslowerrate
the binary star. Our analysis of all other resonant angles of the ofmigration,acaptureinhigherorderMMRsismoreexpected.
system, Φ , Φ , and Φ showed that these all have circulating These authors found that Kepler-16b could have undergone a
3 4 5
behavioursaswell,sometimesprograde(asΦ )andsometimes temporarycaptureina6:1MMRwiththecentralbinary.Inour
1
retrograde.ForallcasesshowninFig.7,theplanetremainscap- simulations,theplanetisalwayscapturedina5:1resonanceand
turedinthe5:1resonanceforthedurationofthesimulations.We inthesameconfiguration,whereonlyΦ2librates.Theeccentric-
also studied the behaviour of the longitude of periastron of the ityofthediskalsoreachesthesamevalueinallthreecases.
binarystar,(cid:36) ,duringtheslowinwardmigrationoftheplanet.
bin
Ouranalysisindicatedthatpriortothecaptureoftheplanetina
3.4.3. Planetsindiskswithdifferentinnerboundary
5:1 MMR, (cid:36) was precessing in a prograde sense. After cap-
bin conditionsatr
tureintothe5:1resonance,thisprecessionbecameretrograde. min
Totestthesensitivityofoursimulationstothechoiceofthedisk
boundary conditions at r , we ran three different models. In
3.4.2. VariationofDiskmass min
twomodelsweconsideredr =0.25AUandonceweassumed
min
Up to this point, the presented results and analyses were based the inner boundary to be closed to mass in-flow, and once we
onsimulationsinwhich,asshowninFig.1,thediskhadagiven considered it to be open (see model 2 in section 3.2). We also
mass.Toexaminewhethertheresultswilldependonthemassof ran a third model for which we moved r to a different inner
min
thedisk,wecarriedoutsimulationswherethedisksurfaceden- radiusat0.22AU(seemodel3insection 3.3).Fig.14showsthe
sity was reduced by a factor of two and four. Results indicated evolution of the planets semi-major axis for these three cases.
9
WilhelmKleyandNaderHaghighipour:ModellingCircumbinaryPlanets:ThecaseofKepler-38
0.147
Model 1: a0=1.00 6000
0.50 Disk Mass
0.1468 0.25 Disk Mass 5000 iinniittiiaall
axis [AU] 0.1466 2y [g/cm] 4000 iiss oorraatthhddeeiiaarrmmttiivvaaeell
major 0.1464 densit 3000
mi- ce 2000
ary Se 0.1462 Surfa 1000
n
Bi
0.146 300 00
2500
0.1458
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Time [Years] e [K] 2000
ur
Fig.15.Theevolutionofthesemi-majoraxisofthebinaryina erat 1500
system with an isothermal disk for different values of the disk mp
e 1000
mass.ThegraphsusesthesamemodelsandcolorsasinFig.13 T
fortheevolutionoftheplanets. 500
0
0.3
The planet was started at a = 1.0 AU in an already relaxed
0 0.25
disk in all three cases. See Figs. 5 and 6 for the density and y
eccentricitydistributionofmodels2and3. ntricit 0.2
Due to the higher density in the simulation with a closed cce 0.15
inner boundary (shown in blue), the planet migrates inward at k E
amuchfasterspeed.Similartotheprevioussimulations,inthis Dis 0.1
casetheplanetiscapturedina5:1MMRwiththebinary.Atthis
0.05
state,theplanetsresonanceangleΦ begantolibratearoundzero
2
and its other resonance angles start circulating. In contrast to 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
theprevioussimulations,thefastinwardmigrationoftheplanet
Radius [AU]
causes its orbital eccentricity to be slightly higher (e = 0.25)
p
thaninthecasewhenthediskhadanopeninnerboundary.This Fig.16.Graphsofthesurfacedensity(ing/cm2,toppanel),tem-
largervalueofeccentricityleadstoorbitalinstabilityasaresult perature (in K, middle panel) and eccentricity (bottom panel)
ofwhichtheplanetleavestheresonanceattimet ≈2300andis ofaradiativediskintheinitialmodel(red),anisothermalcase
scatteredoutintoanorbitwithmuchlargersemi-majoraxisand with H/r = 0.05(model1,blue),andaquasi-equilibriumstate
veryhigheccentricity. (green).
In the simulations where the open inner boundary has been
shifted to r = 0.22AU, the planet migrates inward at a re-
min
ducedspeedandreachesa =0.47AUattheendofthesimula- diskandtheplanet.InFig.15,weplottheevolutionofthesemi-
p
tion.Thisisinanexcellentagreementwiththeobservedorbitof majoraxisofthecentralbinaryinthreesystemswithdiskswith
Kepler-38b.Theeccentricityoftheplanetatthisstateoscillates masses as in Fig. 13. One can see from this figure that during
between 0.01 and 0.06 with an average around 0.03. These re- theevolutionofthesystem,thesemi-majoraxisofthebinaryre-
sultsconfirmtheexpectationthatthefinalpositionoftheplanet duceswithtime.Thiscanbeexplainednotingthatwhilethedisk
is determined by the location of the outer edge of the central interactswiththebinary,angularmomentumistransferredfrom
cavity. thebinarytothediskmaterialcausingtheorbitofthebinaryto
Insummary,inthemodelinwhichtheouteredgeofthecen- shrink.Asexpected,diskswithsmallermassesleadtoaslower
tralcavityisatafartherdistanceandhastheunrealisticbound- declinein theorbit ofthe binary.At thepoint wherethe planet
aryconditionofbeingclosedtomassin-flow,theplanetcrosses iscapturedintoresonance,thebinarylosesangularmomentum
the5:1MMRandeventuallybecomesunstable.Whentheouter morerapidly.
edgeofthecentralcavityisatacloserdistance(e.g.,0.22AU)
and has the more realistic boundary condition of being open to
4. Theradiativecase
massin-flow,theplanetstopsitsinwardmigrationoutsideofthe
5:1 MMR and maintain its orbital stability for a long time. We As mentioned earlier, although the study of a locally isother-
willreturntothispointwhenwediscussradiativemodelsinthe mal disk can portray a general picture of the evolution of the
nextsection. system and its planet, it does not present a realistic model. For
instance,inanisothermaldisk,despitethestrongheatingaction
ofthespiralwavesthatareinducedbythebinarystar,thetem-
3.4.4. Theorbitalelementsofthebinary
perature does not vary and is held constant. However, because
Asmentionedearlierinthedescriptionofourmodelsandtheir the migration of the planet and the evolution of its eccentricity
setups, the motions of all objects (including the stars of the bi- depend on the disk temperature (Ward 1988), a more realistic
nary)areaffectedbythegravityofthedisk.Asaresult,theorbit modelshouldallowthedisktemperaturetovary,aswell.Inthis
ofthecentralbinarywillchangeintimeasitinteractswiththe section, we consider the latter. Because the evolved quantity in
10