Stellar irradiated discs and implications on migration of embedded planets II: accreting-discs

6 MB·

by  Bertram Bitsch

Astronomy&Astrophysicsmanuscriptno.Mdot2D (cid:13)c ESO2014
January14,2014
Stellar irradiated discs and implications on migration of embedded
planets II: accreting discs
BertramBitsch1,Alessandro Morbidelli1,ElenaLega1,Aure´lienCrida1
UniversityofNice-SophiaAntipolis,CNRS,ObservatoiredelaCoˆted’Azur,LaboratoireLAGRANGE,BP4229,06304NICEcedex
4,FRANCE
Preprintonlineversion:January14,2014
4 ABSTRACT
1
0 Context. The strength and direction of the migration of embedded low mass planets depends on the disc’s structure. It has been
2 shownthat,indiscswheretheviscousheatingisbalancedbyradiativetransport,themigrationcanbedirectedoutwards,aprocess
thatextendsthelifetimeofgrowingplanetaryembryos.
n
Aims.InthispaperweinvestigatetheinfluenceofaconstantM˙-fluxthroughthedisc,aswellastheinfluenceofthedisc’smetallicity
a
onthedisc’sthermodynamics.WefocusonM˙ discs,whichhaveanetmassfluxthroughthem.Utilizingtheresultingdiscstructure,
J
wedeterminetheregionsofoutwardmigrationinthedisc.
3 Methods. Weperformednumerical hydrosimulations of M˙ discswithviscousheating, radiativecooling, and stellarirradiationin
1 2Dinther-zplane.Weusedtheexplicit/implicithydrodynamical codeFARGOCAthatincludesafulltensorviscosityandstellar
irradiation, as well as a two-temperature solver that includes radiation transport in the flux-limited diffusion approximation. The
]
migrationofembeddedplanetsisstudiedbyusingtorqueformulae.
P
Results. ForadiscofgassurfacedensityΣ andviscosityν,wefindthatthedisc’sthermalstructuredependsontheproductΣ νand
E G G
theamountofheavyelements,whilethemigrationofplanets,besidesthementionedquantities,dependsontheamountofviscosityν
h. itself.Asaresultofthis,thediscstructurecannotbeapproximatedbysimplepowerlaws.Duringthelifetimeofthedisc,thestructure
p ofthediscchangessignificantlyinanon-linearwayintheinnerparts.Inthelatestagesofthedisc’sevolution(characterizedbylow
- M˙),outwardmigrationisonlypossibleifthemetallicityofthediscishigh.Forlowmetallicity,planetswouldmigrateinwardsand
o couldpotentiallybelosttothestar.
r Conclusions. Thepresenteddiscstructuresandmigrationmapshaveimportantconsequencesontheformationofplanets,sincethey
t
s cangivehintsonthedifferentformationmechanismsfordifferenttypesofplanetsasafunctionofmetallicity.
a
[ Keywords.accretiondiscs–planetformation–hydrodynamics–radiativetransport–planetdiscinteractions–stellarirradiation

2
v
1. Introduction embeddedbodies.Roughlyspeaking,alocalincreaseintheas-
4
3 pect ratio H/r of the disc results in inward migration, while a
In recent years, it has been shown that the migration
3 local decrease in H/r leads to outward migration. As a result,
of low mass planets (≈ 10 − 20M ) can be changed
1 E planetscanmigrateoutwardsinshadowedregionsofthedisc.
from inwards to outwards in discs with heating and cool-
1. ing (Paardekooper&Mellema 2006, 2008; Baruteau&Masset InPaperI,wepointedoutthattheopacityin thediscplays
0 2008; Paardekooper&Papaloizou 2008; Kley&Crida 2008; animportantroledeterminingthediscstructure,sincetheopac-
4 ityisrelevantfortheabsorptionofstellarirradiationandforthe
Kleyetal.2009;Ayliffe&Bate2010).Theheatinginthesepre-
1 cooling of the disc. We only considered equilibrium discs that
viousworksisprovidedbyviscousfriction,whileforthecool-
: havenonet-massfluxthroughthedisc.Inthiswork,wewantto
v ing,alocalcoolingrate(in2Dsimulations)oraradiativediffu-
expandto discs with a constantmass flow M˙ throughthe disc.
i sion(in3Dsimulations)isutilized.Thischangeofmigrationhas
X Inprinciplea viscousaccretiondiscisonlyaccretinginthein-
importantconsequencesfortheformationofplanets,asitwould
r nerpartsofthedisc,whileitisviscouslyspreadingintheouter
provide for a zero-migration radius in the discs, where plan-
a partsofthedisc(Lynden-Bell&Pringle1974)andactuallybe-
ets can survive (Lyraetal. 2010; Bitsch&Kley 2011). Along
havesasadecretiondiscthere.If M˙ stronglydependsonr,the
with that, these points of convergentmigration in the disc can
surfacedensityprofilevaries,sothatthehighM˙ regionsempty,
trap planetary cores in resonances (Cossouetal. 2013). But
andthelowM˙ regionsfill,whichsmoothsthevariationsinM˙ on
these cores could then break free from the resonance (due to
aquickertimescaleforsharpervariations.Thus,itisreasonable
turbulence), then collide and form the cores of giant planets
toassumethat M˙ isalmostindependentofratagivenmoment
(Pierensetal.2013).
in theinnerpartsofthe disc,onwhichwe focushere.Thisac-
Only recently have Bitschetal. (2013) (hereafter Paper I)
cretionrateM˙ canactuallybemeasuredinobservablediscs.
showntheimportanceofstellarirradiationonthediscstructure
Bystudyingdifferent M˙ ratesfordiscs,weexplorethedisc
andthemigrationoflowmassplanetsinnumericalsimulations.
structures as a function of disc evolution, where M˙ decreases
InPaperI,wedidshowthatthechangeofthediscstructurefrom
with time. Additionally, we study the influence of metallicity
ashadowedtoaflaringdiscgreatlyinfluencesthemigrationof
on the disc structure, because the dust-to-gas ratio can change
Sendoffprintrequeststo:B.Bitsch, in time as the disc evolves. The disc’s dust-to-gas ratio could
e-mail:[email protected] decreaseiftheradialdriftgetsridofthedustparticlesfasterthan
1
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
the gas falls onto the star or if planetesimals form efficiently,
1
binding the dust. But it could increase as the gas is removed,
and dustis left behindor is regenerated(e.g.due to collisional 0.5
grinding);intheextremecase, thisleadstodebrisdiscs, where
0
thedust-to-gasratioisinfinite.
Additionally,we investigatethe influence of the disc struc- g) -0.5
tureonthemigrationofembeddedplanets.Since3Dsimulations 2m/
c
ofplanetsarecomputationallyveryexpensive,werelatebackto n -1
usingtorqueformulaeforprescribingtheexpectedmigrationof κg ( i -1.5
planetsin the disc. We thereforeare able to lay outa complete lo
historyofmigrationduringtheevolutionofprotoplanetarydiscs, -2
whichwouldnotbepossiblein3D. Transition
Thepaperisstructuredasfollows.First,wegiveanoverview -2.5 κ inκ BR2 =0 1κ3P
overthe numericalmethodsused in this work in section 2. We * κ*
-3
thendiscussaboutthediscstructureandplanetarymigrationin
10 100 1000
discs that have a metallicity of 0.01 in section 3. The implica-
T in K
tions of different rates of metallicity on the disc structure and
migrationisdescribedinsection4.Wethenpointoutwhatroles Fig.1. Opacity profiles of the opacities used in the code. The
viscosityandgassurfacedensityplayforthemigrationofplan- different bumps in the profile correspond to opacity transition,
etsindiscswithconstant M˙ ratesinsection5.Insection6, we where the ice line is located at ≈ 170K. The grey area marks
discusstheimplicationsofourresultsontheformationandmi- the region where the opacity changes due to the melting of
gration of planets in accretion discs. Finally, we summarize in ice grains. The maximumtemperature during the simulation is
section7. reachedin the inner parts of the disc, where viscousheating is
dominant and can reach up to ≈ 900K for the 1×10−7M /yr
⊙
case.Opacitiesareplottedforametallicityof0.01.
2. Methods
The protoplanetary disc is treated in this study as a three- 2.1.Opacity
dimensional(3D),non-self-gravitatinggaswhosemotionisde-
Theopacityis a crucialparameterwhenusing anenergyequa-
scribed by the Navier-Stokes equations. Without any perturber
tionwithstellarheating.ThePlanckmeanopacityκ isusedin
like planets, the disc is an axisymmetric structure. We there- P
thecoupledtwo-energyequations(radiativeenergyandthermal
fore can use only one grid cell in azimuthal direction, making
energy),theRosselandmeanopacityκ isusedfortheradiative
the computational problem de-facto 2D in the radial and ver- R
cooling,and the stellar opacityκ isused forthe absorptionof
tical (r − z plane) directions, where we utilize spherical coor- ⋆
stellarphotonsintheupperlayersofthedisc.Formoreinforma-
dinates (r-θ). The colatitude θ is measured in such a way that
θ=90◦isthemidplaneofthedisc.Aclearpictureofthesource tion on how the differentopacities are connectedto the energy
equation,seePaperI.
of turbulenceinside accretion discs is understrong debate (see
In Fig. 1, the three different opacities used in the code are
e.g. Turneretal. (2014)). Nevertheless, some sort of viscosity
displayed.FortheRosselandmeanopacityκ andforthePlanck
is needed to drive the accretion disc, so we treat the viscosity R
meanopacityκ ,weusethesameopacitylaw,sincetheseopac-
withanα-prescription(Shakura&Sunyaev1973).Thedissipa- P
ities do not differ that much in the temperature region we are
tiveeffectscanthenbedescribedviathestandardviscousstress-
interestedin(PaperI).However,comparedtoourpreviouswork
tensorapproach(e.g.Mihalas&WeibelMihalas1984).Wealso
inPaperI,wechangedthestellaropacityκ .InPaperI,wetook
includetheirradiationfromthecentralstar,whichwasdescribed ⋆
the stellar opacity to be 0.1 of the Rosseland mean opacity, as
in detail in PaperI. Forthat purposewe modifiedand substan-
indicatedinFig.1.However,thestellaropacitydependsonthe
tially extended an existing multi-dimensional hydrodynamical
temperatureofthestarT andnotonthetemperatureofthegas,
codeFARGOCA,aspresentedinLegaetal.(2013). ⋆
sincedustgrainsabsorbphotonsregardlessoftheirtemperature.
The radiative energy associated with viscous heating and
Asaconsequence,κ isshownasaconstantlineinFig.1.The
stellar irradiation is then diffused through the disc and emitted ⋆
valueofκ =3.5cm2/gcorrespondstothetypicalvalueinadisc
from its surfaces. To describe this process we utilize the flux- ⋆
with a mixture of ice and silicate grains for T = 5600K.The
limiteddiffusionapproximation(FLD,Levermore&Pomraning ⋆
changein κ comparedto PaperI will leadto a higherlevelof
1981),anapproximationthatallowsthetransitionfromtheopti- ⋆
absorptionofstellarphotonsintheupperlayersofthedisc,be-
callythickmid-planetothethinregionsnearthedisc’ssurface.
cause the optical depth τ = κ ρ∆r is much larger. Eventually,
The hydrodynamicalequations solved in the code have al- ⋆
since the opacitywas verylow in the upperregionsof the disc
ready been described in detail (Kleyetal. 2009), and the two-
inthePaperIstudy,thiswillleadtoflareddiscsformuchlower
temperatureapproachforthestellarirradiationwasdescribedin
gasdensitiesthaninPaperI.Butatthesametime,themainar-
detailinPaperI,sowerefrainfromquotingithereagain.Asin
gumentsaboutthe disc structure (flared vs. non-flared)and for
Paper I, we take R = 3R and T = 5600K, which gives the
⋆ ⊙ ⋆ themigrationstillholdinthesenseofthePaperI.
fluxfromthestar:
Inthefollowing,wedefinemetallicityastheratioofheavy
R2σT4 elementstogas.Inthecondensedform,asiceorsilicategrains,
F⋆ = ⋆r2 ⋆ , (1) theseheavyelementsare,inourwork,onlyµminsize.We de-
fine Σ as the surface density of heavy elements in vapour or
Z
where σ is the StefanBoltzmannconstantand r the distance to micrometergrainformandΣ asthegassurfacedensity.Thus,
G
thestar.Stellarheatingisresponsibleforkeepingthediscflared themetallicityZ is simplytheratioZ = Σ /Σ , assumedtobe
Z G
intheouterpartsofthedisc(PaperI). independent of r in the disc. This means that if grain growth
2
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
occursandthetotalamountofheavyelements(independentof 2.4.Discevolution
size)staysthesame,thenthemetallicityinoursenseisstillre-
Theinitialconditionsprovidedin section2.2giveaninitialas-
duced.
pectratioofthedisc,(H/r) ,whichhastobesettoahighervalue
0
astheequilibriumH/rinordertogetaflareddisc,becauseanon
2.2.Surfacedensity flareddisc can notbecomeflared anymore(Dullemond2002),
see Fig. 2. In time the viscosity would then drop (as H drops
Inasteadystateaccretiondisc,theinwardvelocityofthegasis andν = αH2Ω ),untilanequilibriumstateisreached.Starting
K
asimulationlikethatwouldimplythatH,thereforeΣ,changea
3ν lotintimetoguaranteeaconstant M˙.If H dropsbyafactorof
v =− , (2)
r 2r 2−3asobservedinPaperI,thesurfacedensitythereforehasto
rise bya factorof 4−9. The time for the disc to adjust to that
whereν=αH2Ω istheα-viscosity(Shakura&Sunyaev1973) wouldbe extremelylong,since the disc has to be refilled from
K
withHtheheightofthedisc,Ω theKeplerianorbitalfrequency, the outerboundaryona viscoustimescale,whichis essentially
K
and r denotes the orbital distance. With the radial velocity we thelifetimeofthedisc.
candefineanaccretionrateM˙ as Weimposeseveralstepsinthesimulationstosavecomputa-
tiontime.Wefirstkeeptheviscosityconstantintimeandequal
M˙ =−2πrΣGvr =3πνΣG , (3) toitsinitialvalueν0 =αh20r15/14,byusingtheinitialh0 =(H/r)0
configuration to determine the viscosity during the whole first
where Σ is the gas surface density. The slope of the surface step. As the disc evolves, it compresses towards midplane and
G
densitycanbecalculatedusingtheviscosity we reach an equilibrium state with a flared disc profile in the
outerpartsofthedisc.Wethenfitthisprofile,whichislabelled
(H/r) inFig.2,withanewH/rvalue,whereh ∝ r2/7 (shown
ν=αH2Ω∝αh2r2ar−3/2 (4) 1 1
0 as (H/r) in Fig. 2). With this new estimate of H/r, we can
1,fit
compute a new α value for the viscosity by comparing it to
with h = H/r and where a = 9/7 as H/r ∝ r2/7 (H/r) .ThechangeofH intheviscosityiscompensatedbythe
0
(Chiang&Goldreich (1997), Paper I). The proportionality, ∝
changeofα, so thatthe disc will havethe same viscosity as in
r2/7 describes the flaring index of the disc, which gives for a thefirststepandthereforetheM˙ rateisconservedaswell.
steady state accretion disc, where M˙ is constant for all radii
Thesimulationisthenrestartedwiththenewαvalue(which
((d/dr)M˙ =0), is given by α = α h2/h2), and the viscosity ν is now given
end 0 0 1
by eq. 4, with the local value of H to account for bumps and
M˙ ∝h20r2ar−3/2ΣG,0r−sM˙ , (5) dipsthatmighteventuallyformin thedisc.Forthe outerparts,
these changesof ν are notimportantbecause the disc structure
wheresM˙ denotesthepowerlawindexofthegassurfacedensity is determinedby stellar irradiation and not by viscousheating.
inthe M˙ disc.Here,ΣG,0 isthevalueofthegassurfacedensity In the inner parts, where viscous heating is the dominant heat
atr =1AU.Thisleadsto source,theviscositychangescomparedtothe(H/r) case,lead-
0 0
ingto changein thedisc structure.However,these changesare
2a−3/2−s =0 ⇔ 18/7−3/2= s ⇔ s =15/14.(6) minimal. The advantage of this procedure is that the M˙ rate is
M˙ M˙ M˙
thesameasintheinitialsimulationwithnoneedtorefillΣ as
G
Weusethisslopeofsurfacedensityfortheinitialconditionsin Hdrops.
all our simulations. Vertically, we initially set the density ρ of Wethenreachanewequilibriumstate(H/r)2,wherethein-
thedisctobeinhydrostaticalequilibrium nerpartsofthedischavechangedcomparedto(H/r)1,because
theviscositychanged.But,sincethediscstillneedstore-arrange
thesurfacedensityonaviscoustimescale,werelatebacktoalo-
(rcosθ)2
ρ(r,θ)=ρ (r) exp − , (7) callyisothermalsimulation,wherewepluginthediscstructure
0 " 2H2 #
and sound speed as shown in the (H/r) disc. The disc is then
2
evolved until a time independentmass flux through the disc is
whereρ (r)istheinitialdensity.
0 achieved. After that, we go back to the fully radiative energy
Forallsimulationsanadiabaticindexofγ=1.4andamean equation with stellar irradiation to account for changes in the
molecularweightofµ=2.3isset. heating due to the re-arranged gas surface density Σ and ar-
G
rive at (H/r) , which gives the final state of the disc. Please
final
note that this whole process can only work as the outer disc is
2.3.Boundaryconditions
supportedbystellarirradiation,keepingtheinputofmassatthe
To get a disc with a given M˙ rate, we impose an M˙-rate at the outer boundary constant, not only in time, but also in vertical
outer boundary,while we leave the inner boundary open. This spacingbecauseHstaysconstant.
way,thedensitystructureinsidethedisccanchangefromtheini- In this way, much computation time is saved, because the
tialconditionsandformbumpsinthesurfacedensitythatcould computationtimeismainlydominatedbythestellar-irradiation
notbecreatedifanM˙ boundarywasappliedattheinnerbound- algorithm that can take up to 90% of the whole computation
aryaswell,becausethetotalmassinsidetheM˙ discwouldstay time.Especiallyinthestatewherethesurface-densityhastobe
constantandnorefillingcouldtakeplace.Theimportantquanti- re-arranged,duetochangesinviscositythismakesabigdiffer-
tiesthathavetobetreatedinaspecialwayattheboundariesare ence.
theradialvelocityv andthegassurfacedensityΣ intheghost Inprinciplewithtimethemetallicityofsmalldustgrainsin
r G
rings. The details for the boundary conditions can be found in the disc can change(e.g. due to grow of planetesimals),which
AppendixA. wouldchangethediscstructure(seesection4).Thiswouldim-
3
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
0.16 haviourofplanetarycoresinsidetheM˙ discsbyusingthetorque
((HH//rr))10 formulaprovidedbyPaardekooperetal.(2011).
0.14 (H/(rH)1/r )fi2t
(H/r)final
3.1.Discstructure
0.12
InFig.3,the H/r,temperature,andgassurfacedensityprofiles
0.1 of discs with different M˙ are displayed.The H/r profilesshow
H/r a flaring partin the outerdisc and somebumpsanddipsin the
0.08 innerpartofthedisc.Thesebumpsanddipsarecausedbytran-
sitions in the opacity of the disc, for example, at the ice line.
0.06 This kind of structure was also observed for equilibrium discs
withconstantviscosityinPaperI.Theheightofthebumpsand
0.04 thedepthofthedipsinthediscisdecreasingwithdecreasingM˙
rate(whichmeansdecreasingΣ ).As M˙ decreases,sodoesthe
G
0.02 viscousheating,hencethetemperatureandH/rintheinnerpart
10 20 30 40 50
ofthedisc.Thisalsomeansthattheshadowedregionofthedisc
r [AU]
isbecomingsmallerasthediscsdensityreduces.Likewise,the
Fig.2.EvolutionofH/rinthediscforthedifferentstages.The bumpsofthedisc structuremovetowardsthecentralstar as M˙
discfeaturesM˙ =2×10−7M /yr.Atthebeginning,inthe(H/r)
⊙ 0 decreases.Thebumpsinthedisccomefromtheopacitytransi-
stage,α =0.001,whichwechangetoα =0.003afterthestage
0 1 tion,whichalwayscorrespondstothesametemperatureregion.
(H/r) isreached.Thestage(H/r) showstheaspectratioafter
1 2 Thereforelessheatingintheinnerpartsofthediscmovesthese
theviscosityhaschangeand(H/r) showsthefinalstageafter
final bumpsfurtherinside.
thesurfacedensityhasbeenre-arranged.
Because the gassurface densityis reducedmoreand more,
the bumps in the disc become smaller and smaller (as the disc
heatsless)untiltheyfinallydisappearcompletelyforlowM˙.For
ply thatdisc needsto re-arrangethe gas inside, whichhappens the 1×10−8M /yr disc there is only a verysmall bump in the
⊙
on a viscous timescale. However, as we show below, this re- H/r profilevisible. Thisbumpvanishesforthe 5×10−9M /yr
⊙
arrangingofmassismuchlesspronouncedintheouterpartsof disc,whichonlyshowsanincreaseinH/r.Intheouterpartsof
thedisc,wheretheviscoustimescaleismuchlongercompared the disc, the structure is very similar for all M˙. The heating in
totheinnerparts(τ = r2/ν = r2/(αH2Ω )).Thisisbecause theouterdiscisdominatedbystellarirradiation,whichdepends
visc K
the change in viscosity in the outer parts is much smaller than ontheopticaldepth(∆τ = ρ κ ∆r) ofthedisc. Aslongasthe
G ⋆
thechangeintheinnerpartsofthedisc,becausetheH/rprofile disc is opticallythick (∆τ > 1), it can efficientlyabsorbstellar
staysconstantasitisdominatedbystellarirradiation.Thisindi- irradiation and be heated. The disc is then flared with H/r ∝
catesthatthesteadystateofconstant M˙ canbeachievedinthe r2/7 asfoundfortheoreticalcalculationsinChiang&Goldreich
innerpartsonatimescalethatisshorterthanthelifetimeofthe (1997).
disc,whichisdeterminedbytheviscoustimescaleofthe outer Inhydrostaticequilibrium,thetemperatureisrelatedinmid-
partsofthedisc. planetoH/rthrough
H 2 GM µ
T = ⋆ , (8)
2.5.Numericalsetup (cid:18) r (cid:19) r R
We performed several simulations for different M˙ rates corre- where µ is the mean molecular weight and R the gas constant.
spondingtodifferentstagesofthediscs’lifetimethatarelisted The temperature profile (middle in Fig. 3) therefore reflects
in Table 1. In the first sets of simulations, we assumed thatthe thefluctuationsinthe H/r profilebyshowingsteepergradients
dust is bound perfectly to the gas and as the gas accretes onto whenH/rdropsandflattergradientswhenH/rincreases.These
thestar,thedustdecreasesatthesamerate,keepingthedust-to- changes are caused by the transition of opacity. The region of
gasratioconstant.Inthesecondsetofsimulations,weassumed change in the opacityis marked in grey in the middle panelof
thatthedustisnotperfectlyboundtothegasatthelaterstages Fig.3.
of the disc’s evolution, resulting in an increasing and decreas- Thesurfacedensityisreducedbytheappropriatefactorfor
ingdust-to-gasratio,dependingonthefavouredscenarioforthe changing the M˙ rate. In the outer parts of the disc, this can be
evolutionofdustparticles.
Allsimulationswereperformedonagridwith386×66active M˙-rate[M /yr] Metallicity Σ α
⊙ G,0 end
cells in r−θ direction,where the openingangleof the numer- 2×10−7 0.01 5.385×10−3 0.003
ical simulationis 20◦, meaningthat70◦ < θ < 90◦. The radial 1×10−7 0.01 2.693×10−3 0.003
domainextendsfrom1AUto50AU. 5×10−8 0.01 1.346×10−3 0.003
1×10−8 0.01 2.693×10−4 0.003
5×10−9 0.01 1.346×10−4 0.003
3. M˙-discswithconstantgas-to-dustratio
5×10−8 0.02 1.346×10−3 0.003
Inthissectionweanalysethediscstructurefordiscswithdiffer- 1×10−8 0.03 2.693×10−4 0.0026
entM˙ rates,butwithaconstantmetallicityof0.01.Thedifferent 5×10−9 0.05 1.346×10−4 0.0026
M˙ ratescorrespondtodifferentevolutionarystatesofthedisc’s 5×10−8 0.001 1.346×10−3 0.0026
1×10−8 0.001 2.693×10−4 0.0026
lifetime. We reduce the surface density of the gas, Σ , to real-
G
izedifferentstatesoftheevolution,butkeeptheviscosityatthe Table1.Simulationparametersfortheusedmodels.
samevalueforallsimulations.Wethenanalysethemigrationbe-
4
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
0.08 M˙ =1×10−7 byasmaller H,resultsinhighersurfacedensitybecause M˙ has
M˙ =5×10−8 tobeconstant,thusexplainingthesurfacedensityprofiles.
M˙ =3×10−8 Additionally, there are wiggles in the Σ profile at ≈ 5AU
0.07 M˙ =1×10−8 G
M˙ =5×10−9 for all M˙ discs. This also can be explained by changes in the
r2/7fit viscosityduetochangesin H (seetopofFig.3),and M˙ hasto
0.06 beconstantthroughoutthedisc.Thechangesofthegradientsin
/r thesurfacedensityΣG andinthetemperatureT haveimportant
H
consequencesforthemigrationinsidethesediscs.
0.05
3.2.Migrationmaps
0.04
To estimate the torqueactingon planetsembeddedin thediscs
with different M˙, we use the formula by Paardekooperetal.
0.03
(2011), which captures the effects of torque saturation in con-
2 5 10 20
trast to Paardekooperetal. (2010), where the torques are fully
r[AU]
1000 unsaturated.However,theformulabyPaardekooperetal.(2011)
mightnot be accurate for low-mass planets. In fact, Legaetal.
500 (2013) have found a new additionalnegative torque, which di-
minishes the total torque acting on embedded low-mass plan-
ets. But these differences do not seem to be that great, so in
200
the absence of a formula that would capture this effect and to
K avoid large scale numericalsimulations, we use the formulaof
n 100
i Paardekooperetal. (2011), which givesa good approximation.
T
TheformulacapturesthetorquecausedbyLindbladresonances
50 Transition
M˙ =1×10−7 and horseshoe drag on low-mass planets embedded in gaseous
M˙ =5×10−8 discs in the presence of viscous heating and thermal diffusion.
M˙ =3×10−8
M˙ =1×10−8 Theformuladoesnotincludestellarheating,butstellarheating
M˙ =5×10−9 onlychangesthediscstructureandnotthemechanismresponsi-
r−1fit
10 bleforoutwardmigration(PaperI).Thisformulawasalsotested
2 5 10 20 against 3D simulations in Bitsch&Kley (2011), who found a
r[AU] goodagreementfor20M planets.
Earth
TheformulaofPaardekooperetal.(2011)isverycomplex,
1000
so we do not explain the whole formula here. The total torque
500
acting on an embeddedplanet is a compositionof its Lindblad
200 torqueanditscorotationtorque:
2m 100 Γtot =ΓL+ΓC . (9)
/c 50
g
n The Lindblad torque depends on the gradients of temperature
i
Σ T ∝ r−β and gas surface density Σ ∝ r−s. It is given in
M˙ =1×10−7 G
10 M˙ =5×10−8 Paardekooperetal.(2011)by
M˙ =3×10−8
M˙ =1×10−8 q 2
M˙ =r−515×/1410fi−t9s γΓL/Γ0 =−2.5−1.7β+0.1s and Γ0 =(cid:18)h(cid:19) ΣPr4pΩ2P, (10)
1
where q is the mass ratio between planet and star, Σ the gas
2 5 10 20 P
surface density of the disc at the planet’s location, and r the
r[AU] P
distance of the planet to the host star. One can clearly see that
Fig.3.H/r-profile(top),temperature(middle),andsurfaceden- achangeinthegradientoftemperatureinfluencestheLindblad
sityprofile(bottom)ofdiscswithdifferentM˙ rates.Theprofile torque.Thesameappliestothecorotationtorque,whichstrongly
is taken when the discs have reached their radially constant M˙ dependsonthegradientofentropy,S ∝ r−ξ,withξ = β−(γ−
rate.Theblackdottedlinesrepresentpower-lawfitstoguidethe 1.0)s. The largest contribution of the corotation torque arises
eyes.Thegreyareain the temperatureplotmarksthetempera- fromtheentropyrelatedhorseshoedrag,whichisgivenby
turerangewheretheopacityprofilechangesduetothemelting
ofice grains(see Fig. 1). Please notethatwe cutthedisplayed ξ
γΓ /Γ =7.9 . (11)
discat30AUtoenhancethedetailsintheinnerpartsofthedisc. hs,ent 0 γ
This indicates that we expect a change in the migration rate,
when the temperaturechangessignificantly. This is the case in
observedwell(bottominFig.3).However,intheinnerpartsof the inner parts of the discs, where the aspect ratio H/r shows
thediscthisdoesnotseemtobethecaseexactly.Forexample, fluctuations (as H/r is proportional to T). In fact, in Paper I
between the 1×10−8M /yr disc and the 5×10−8M /yr disc, we found that outward migration only seems possible in re-
⊙ ⊙
thedifferenceinsurfacedensityattheinnerboundaryisonlya gionsofthediscwhere H/r drops.Thetorqueformulathatac-
factorof≈2insteadof5.Thiscanbeexplainedbythedifferent counts for saturation effects (Paardekooperetal. 2011) will re-
viscosity in the inner partsof the disc, since ν = αH2Ω and H sultin a smaller regionof outwardmigrationthan with the un-
clearlyvariesforthedifferentM˙ values.Lowerviscositycaused saturatedtorqueformula(Paardekooperetal. 2010), whichcan
5
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
clearly be seen in Fig. 4. Additionally, the formula for the un- 50 4
saturated torques is independentof the planetary mass, so that 45 3
outwardmigrationwillbefoundforallplanetarymassesforthe arth 40
Paardekooperetal.(2010)formula,incontrasttothetorquefor- ME 35 2
ppmeruensldCaesnothotmeandpttaahirnceecdpPolautaponnetetsrhtafIeorytrehqmsaautatiuhlsisrabad(rtPiiuosaemnaq,rdbd=eeiksc0coas.ou5psw,eewirtthoeertncqaooul.nwes2st0haa1antuv1tre)va.tiaisocmnousdicteyh- anet mass in 12235050 - 011 γΓΓ / tot0
steepersurfacedensitygradientwith sM˙ ≈ 15/14.Thisreduces Pl 10 -2
not only the entropy gradient in the disc and therefore the en- 5
tropyrelatedhorseshoedrag,butalsothebarotropicpartofthe -3
2 4 6 8 10 12 14 16 18
horseshoedrag,whichisgivenby
r [AU]
γΓ /Γ =1.1(1.5−s) , (12) 50 4
hs,baro 0
45
by ≈ 50% compared to equilibrium discs with s = 0.5. The 3
effectontheLindbladtorque(eq.10),however,isemq inimal.The arth 40
expectedtotaleffectwouldbe areducedregionofoutwardmi- ME 35 2
gratIinonFiing.th4ewaeccprreetisnegntdtihsecsmciogmraptaiorendmtoapthsefoerqduiislicbsriwumithdMi˙sc=s. mass in 223050 01 ΓΓ / tot0
w5×h1ic0h−8thMe⊙d/iysrc,Ms˙tru=ct1u×re1s0a−r8eMd⊙i/syprla,yaenddiMn˙ F=ig5.×31.0F−o9rMt⊙h/eyhri,gfohr- anet 15 -1 γ
est M˙ disc model, we have two separated regions of outward Pl 10 -2
5
migration,whichareencircledbyblacklinesinthefigure.The
-3
firstregionisfromthe inneredgeofthedisc upto≈ 3AU, the
2 4 6 8 10 12 14 16 18
second is from ≈ 5AU< r < 7AU. These regions compare
P r [AU]
nicely to regions in the disc where H/r drops. A drop in H/r
50 4
corresponds to an increase in the temperature gradient result-
45
3
rinelgatiendacnoirnoctaretiaosnedtoernqtureoplyeagdriandgietnot,ohuetwncaerdamstirgornagtieorne.nTtrhoepsye MEarth 3450 2
cfsohhroaewnqgnueisilniibnFriiHugm./r3da(irmsecisdciadnulesP)ea,dpthebreyIg.orIpenayfcaaictryte,atirnmatnhasreikttiseomtnhspe,earresagtwuioraenspwsrtohafiteerldee mass in 223050 01 ΓΓ / tot0
dthieenot.paFcoirtythcehaM˙ng=es,5w×h1ic0h−8sMhow/syradsihsacl,lothwisergtreemypreegraiotunre(fgrorma- anet 15 -1 γ
⊙ Pl 10
≈ 3AU to ≈ 5AU) corresponds nicely to the region of inward 5 -2
migrationshowninFig.4. -3
Compared to the studies of equilibriumdiscs with constant 2 4 6 8 10 12 14 16 18
viscosity (as in Paper I) where the minimummass for outward r [AU]
migrationwas≈ 5M ,itisnow≈ 8M .Thisdifferencehasits
origininthedifferentEdiscstructures.SEincethe M˙ discfeatures Fig.4. Torque acting on discs with different M˙, with 5 ×
a muchsteepersurfacedensitygradient sM˙ ≈ 15/14compared 10−8M⊙/yr (top), 1×10−8M⊙/yr (middle)and 5×10−9M⊙/yr
to the equilibriumdisc s = 0.5, it reduces the differentposi- (bottom)fortheappliedPaardekooperetal.(2011)formula.The
eq
tivecontributionsofthetorque,sothatverylow-massplanetsno black lines encircle the regions of outward migration for the
longer migrate outwards. For low planetary masses, the horse- Paardekooperetal. (2011) formula. The vertical solid yellow
shoedragtendstowardsthelinearcorotationtorque,resultingin linesmarktheouteredgeofthezero-migrationregionfortheun-
inwardmigration.Infact,forlow-massplanets(<5M ),new saturated torques(Paardekooperetal. 2010), while the dashed-
Earth
studies have shown an additional negative torque (Legaetal. yellow line marks the inner edge of the zero migration region
2013). The maximum mass of outward migration, ≈ 40M , oftheunsaturatedtorques.Theverticalredlinesindicatetheice
E
seemstobethesameasinPaperI.However,oneshouldbescep- lineat170K.
ticalabouttheresultsforhigh-massplanets,becausetheformula
of Paardekooperetal. (2011) was derived in the linear regime
forlow-massplanetsandnotforlargeplanetsthatcanevenstart ThemigrationmapforM˙ =5×10−9M /yr(bottominFig.4)
⊙
toopenuppartialgaps. lookscompletelydifferentthanforthehigherM˙ models.Infact,
Wedonotshowthemigrationmapfor M˙ = 3×10−8M /yr there are no more regions of outward migration left. As stated
⊙
because it shows just a transition state between M˙ = 5 × before,outwardmigrationis strongerif H/r drops,but forthis
10−8M /yr and M˙ = 1×10−8M /yr, since they have no new case, no drop in H/r can be observed (see Fig. 3), leading to
⊙ ⊙
features, which can also be guessed from the structures shown inwardmigrationforallplanetarymassesandatallorbitaldis-
Fig.3.However,inthiscasetheinnerandouterregionsofout- tances.However,intheregionwhereweobservedoutwardmi-
wardmigrationshrinkandareshiftedinwards.Infact,whenar- gration for higher M˙ discs, the inward migration is slowest. In
rivingatM˙ =1×10−8M /yronlyoneregionofoutwardmigra- theregionbetween1.5and2AU,theinwardmigrationisactually
⊙
tion remains in our computational domain. This region is then increasedcomparedtotherestofthedisc.Thisissupportedby
muchnarrowerinsizeandvalidforasmallerrangeofplanetary thetemperaturegradientofthedisc,whichactuallyseemstobe
masses. But, it also implies that outward migration is possible partlypositive(seeFig.3),whichcanresultinanegativecorota-
forplanetswithslightlylowermass(≈5M ). tiontorque.But,keepinmindherethattheactualmigrationrate
Earth
6
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
is lower for reduced disc masses, since Γ scales linearly with 0.08
0 M˙ =5×10−8,met=0.001
ΣG,seeeq.10. 0.075 M˙ =1×10−8,met=0.001
Interestingly,onlythetorqueformulathataccountsforsatu- M˙ =5×10−8,met=0.02
rationgivesanegativetorqueforthewholedisc,whilethefully 0.07 M˙ =1×10−8,met=0.03
M˙ =5×10−9,met=0.05
unsaturatedtorqueformula(Paardekooperetal. 2010)still pre- 0.065
dictsoutwardmigration.Infact,inanM˙ disc,onecanrelatethe
0.06
total unsaturatedtorque, which consists of the Lindbladtorque
(eq.10),andthebarotropic-(eq.12)andentropy-related(eq.11) /r 0.055
H
parts of the corotation torque, to the flaring index of the disc.
0.05
Assuming that Σ ν is independantof r, as well as of α (which
G
givess=3/2−β),andthatγ=1.4,onefindsbyusingeq.8 0.045
Γ /Γ =1−10b, (13) 0.04
tot 0
wherebistheflaringindexofthedisc.Thisindicatesthateven 0.035
the unsaturated torque is negative in discs with flaring indices 0.03
greater than 0.1. By interpolating the trend in the H/r profile 2 5 10 20
towards even lower M˙ discs, we speculate that the unsaturated r[AU]
torque formula will predict inward migration as M˙ becomes
Fig.5.H/r-profileofdiscswithdifferentM˙ ratesandmetallici-
lower and lower. However,the torquessaturate in time, so that
ties.Theprofileistakenwhenthediscshavereachedtheirradi-
the formulaforthe unsaturatedtorquesshould notbe used, es-
peciallyforlow-M˙ discs(seediscussioninsection6). ally constant M˙ rate. The colourcodingof the high-metallicity
discs matches the M˙ rates displayed in Fig. 3. Please note that
wecutthedisplayeddiscat30AUtoenhancethevisibilityinthe
4. M˙-discswithdifferentmetallicity innerpartsofthedisc.
Inthissectionweexploretheinfluenceofdifferentmetallicities
on the disc structure and the correspondingmigration maps. A
reducedmetallicityin smallgrainscouldhappenin time asthe
coolingturnstohighertemperatures,hencetohigherH.Thisef-
firstplanetesimalsandembryosforminthediscandreducethe fect,ofcourse,diminishesastheM˙ rateisreducedandtherefore
amountofheavyelementsin thedisc. Additionally,radialdrift
viscousheatingaswell.
couldwashoutthedustgrains,reducingthemetallicity.Ahigher
Interestingly, not only is the aspect ratio greater for higher
metallicityin smallµm size dustgrainsinthe late stagescould
metallicitydiscs,butthebumpsintheinnerpartsofthediscare
beachievedbycollisionalgrindingofplanetesimalsthatproduce
also more pronounced. For the M˙ = 5× 10−9M /yr disc, the
smallparticles. ⊙
bumpsareactuallyvisibleforthefirsttime.Theincreaseinthe
The change in metallicity has many implications, also for
bumps inside the disc structure is related to the reduced cool-
planet formation. For example, a higher metallicity helps the
ing rate in the high-metallicity discs. This can have important
streaming instability operate (Johansen&Youdin 2007) and
implicationsonthemigrationofgiantplanetcores.
planetesimals to form (Johansenetal. 2007), while a lower
In Fig. 6, we present the migration maps for the same M˙
metallicityhindersthestreaminginstabilitytowork.
discs as in Fig. 4, but with different metallicities. The metal-
licities are 0.02 (top), 0.03 (middle), and 0.05 (bottom). The
4.1.Highermetallicity cleardifferenceforthemigrationmapsdisplayedherewith the
onesinFig.4isthattheregionsofoutwardmigrationaremuch
Achangeinthemetallicityofthedisctranslatestoachangein
larger and more extended. Additionally, the migration map for
the opacity of the disc. If the metallicity of a disc is doubled,
the 5 × 10−9M /yr disc shows outward migration, which was
so is the opacity. If the metallicity increases, the cooling rate ⊙
notvisibleforthediscwithametallicityof0.01.Thisisclearly
ofthediscdecreasesascanbeseenfromeq.15.InFig.5,itis
relatedtothechangeinthediscstructureintheinnerparts(see
shownthattheaspectratioH/rofdiscsincreaseswithincreasing
Fig.5),whereabumpnowappears.
metallicities,comparedtothelow-metallicitycases(Fig.3).
Intheouterpartsofthedisc,thestructureisquitesimilarto
theoneshownforthelow-metallicitydiscs(Fig.3).Thisisnot 4.2.Lowermetallicity
obvious,sincethecoolingratereduceswithincreasingmetallic-
itywhiletheabsorptionofstellarirradiationincreases.However, Owingtotheformationofplanetesimalsandplanetaryembryos,
owing to more grains in the disc, the absorptionof stellar irra- whichrequiredustandicegrains,theamountofsmallsolidscan
diation(whichis∝ ρκ )willhappenforlowerdensities,which decrease,if there is noadditionalsourceof small grains.Discs
⋆
impliesthatithappensatahigheraltitudeinthedisc.Therefore, with reduced metallicity have reduced opacity inside the disc,
theheatthattransfersdowntomidplaneisactuallyreduced.But which works in the same way as described in section 4.1. The
thiseffectiscompensatedbythereducedcooling,sotheresult- aspectratioofdiscswithametallicityof0.001arealsodisplayed
ingdiscstructureintheouterpartsisthesame. inFig.5.
Clearlytheaspectratioishigherintheinnerpartsofthedisc Clearly, lower metallicity leads to a higher cooling rate
fortheM˙ =5×10−8M /yrdiscwithametallicityof0.02com- (eq. 15), which explains the drop in H/r in the inner parts
⊙
paredtothediscwithametallicityof0.01showninFig.3.We of the disc. The outer parts of the disc remain flared and at
recallthattheinnerpartofthediscisdominatedbyviscousheat- about a constant level, as explained in section 4.1. For the
ing.Theviscousheatingdoesnotdependonthemetallicity,but M˙ = 1×10−8M /yrdisc,thereisadramaticchangeinthedisc
⊙
onlyonthegasdensity.Thehighermetallicityreducesthecool- structure because no bumps in the discs profile are visible any
ingrate andthereforethe balancebetweenviscousheatingand more,whicharestillvisibleinthecaseofametallicityof0.01.
7
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
50 4 50 4
45 45
3 3
arth 40 arth 40
ME 35 2 ME 35 2
mass in 223050 01 ΓΓ / tot0 mass in 223050 01 ΓΓ / tot0
net 15 -1 γ net 15 -1 γ
a a
Pl 10 Pl 10
-2 -2
5 5
-3 -3
2 4 6 8 10 12 14 2 4 6 8 10 12 14
r [AU] r [AU]
50 4 50 -1.5
45 3 45 -2
MEarth 3450 2 MEarth 3450 -2.5
net mass in 12235050 - 011 γΓΓ / tot0 anet mass in 12235050 ----4433..55 γΓΓ / tot0
Pla 10 -2 Pl 1 50 -5
5
-5.5
-3
2 4 6 8 10 12 14
2 4 6 8 10 12 14
r [AU]
r [AU]
50 4 Fig.7. Torque acting on discs with different M˙ and metallic-
45 3 ity, with 5 × 10−8M⊙/yr and 0.001 metallicity (top) and 1 ×
arth 40 10−8M⊙/yr and0.001metallicity (bottom).The blacklinesen-
ME 35 2 circle the regions of outward migration. The vertical red (top)
net mass in 12235050 - 011 γΓΓ / tot0 itasinoadndbmilffuaeeprs(ebcnoottrctrooemlsop)uolrinnsdecstaolienthdfoeicrdatithseecthbpeorotitcfioelmelsipnslehooat.wt1n7i0nKF.iTg.h5e.mTihgerrae-
a
Pl 10
-2
5
-3 ture,whichinturnleadstoanegativeentropy-relatedcorotation
2 4 6 8 10 12 14
torque, resulting in an even stronger negativetotal torque. The
r [AU] reason for that is probably related to the transition of viscous
Fig.6. Torqueactingon discs with different M˙ and metallicity, heatingto stellar heating.Since viscousheating is operatingin
with5×10−8M /yrand0.02metallicity(top),1×10−8M /yrand the inner (midplane) parts of the disc, heat is produced there.
⊙ ⊙
0.03metallicity(middle),and5×10−9M /yrand0.05metallic- Stellarirradiationheatstheupperlayersofthedisc,whichthen
⊙
transportheatdowntowardsmidplane.At≈2.5AU,stellarheat-
ity(bottom).Theblacklinesencircletheregionsofoutwardmi-
ingisgettingabsorbedeffectively,whichheatsthemidplane,but
gration.Theverticalredlinesindicatetheicelineat170K.The
viscousheatingalreadyseemsto diminishata shorterdistance
migrationmapscorrespondtothediscprofilesshowninFig.5.
to the star. Therefore a higher temperature can be observed at
≈2.5AUcomparedtotheinnerpartsofthedisc,whichleadsto
apositivetemperaturegradientinthedisc.
In Fig. 7 the migration maps for the low-metallicity discs
are displayed. In the M˙ = 5 × 10−8M /yr case (Fig. 7, top),
⊙
the regions of outward migration are much smaller than in the
0.01metallicitycase(topofFig.4).Infact,theinnerpartregion 5. InfluenceoftheinterchangebetweenΣ andν
G
ofoutwardmigrationhasnowcompletelydisappeared,andthe
outerregionisonlyvalidformuchlowermasses.Additionally, BecausethemassfluxM˙ throughadiscisproportionaltothegas
theouterregionofoutwardmigrationisshiftedinwardsfrom≈ surfacedensityΣG andtheviscosityν,achangeinΣG couldbe
5AUto≈ 3.8AU,whichcanbeexplainedbythehighercooling compensatedbyachangeinνtostillhavethesameM˙.Wenow
rateinthedisc,whichshiftstheicelinefurtherinwards. wanttoconstructascenariowhere M˙ andthethermalstructure
In the case of M˙ = 1 × 10−8M⊙/yr (Fig. 7, bottom), the are the same, but where ΣG and ν are not. For changing ν, we
changes are dramatic because no region of outward migration changeα.
existsanymore.Infact,theregioninsidethe iceline(r > r ) Wenowlookatthedifferentheatingandcoolingsourcesof
ice
nowfeaturesaregionofenhancedinwardmigration.Inallother thedisc.Theviscousheating,whichisthedominantheatsource
displayedmigrationmaps,we observea regionofoutwardmi- intheinnerpartsofthedisc,isproportionaltoΣ ν(M˙):
G
gration(oratleastwithslowerinwardmigration).Atthispoint
in the disc, the aspect-ratio profile undergoes a steep gradient 9
(Fig.5),whichactuallyreferstoapositivegradientintempera- Q+ = 8ΣGνΩ2K2πrδr, (14)
8
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
50 4 ofoutwardmigrationisvisible.However,outwardmigrationin
45 3 thiscasestartsatmuchhigherplanetarymasses(12MEarth)than
arth 40 in the normalviscosity case presentedin fig. 4, whereoutward
ME 35 2 migrationstartsfor≈5M .
Earth
net mass in 12235050 - 011 γΓΓ / tot0 icsoneffioTthcleiiennreetaaχrsioinnnρtfhGoeκrRf.othIritmsrueilasadtoshfeaPssacaarldinekgooofpethreettahle.r(m20a1l1d)i,ffwuhsiiochn
Pla 10 16γ(γ−1)σT4
-2 χ= , (16)
5 3κ ρ2H2Ω2
-3 R G P
2 4 6 8 10 12 14 whereσ is theBoltzmannconstant.A factorof 4is missing in
r [AU] Paardekooperetal. (2011) according to Bitsch&Kley (2011).
50 4 Additionally,thesaturationparametersforthecorotationtorque
45 3 and horseshoe drag also depends on the viscosity of the disc.
arth 40 These saturation parameters are responsible for the transition
ME 35 2 fromhorseshoedragtowardsthelinearcorotationtorque,which
anet mass in 12235050 - 011 γΓΓ / tot0 iwdsiasmrcdTsuhmcwihsiitglhneroasdtsniiff-otlhneina,reencnahtrtahiνnteygthoicnoaΣgrusGtsehesreshaotmdieoiiffsgd,errareaetvgineo.tnnTmimhfiigasthprcaesaty.inofnethabetuenrhelaevtahidoeutsoramfinoe-r
Pl 10 M˙ value and the same thermal structure. In total, one could
-2
5 roughly say that lower viscosity allows outward migration for
-3 lower planetary masses (unless it gets too low and no outward
2 4 6 8 10 12 14
migrationispossibleanymore),whileahigherviscosityallows
r [AU] outwardmigrationforhigherplanetarymasses, wherethemin-
Fig.8. Torque acting on discs with M˙ = 1×10−8M /yr for a imum mass required for outward migration is also increasing.
⊙
rescaling of Σ ν by a factor of s = 0.5 (top) and a factor of This is crucial for N-Body simulations of planetary embryos
s = 2(bottomG).Thismeansthatfthedisc displayedontophas in evolving discs, where just stating an M˙ value would not be
f
enough.
ametallicityof0.005,whilethedisconthebottomhasametal-
licity of 0.02. The black lines encircle the regions of outward
migration.Theverticalredlinesindicatetheicelineat170K.
6. Consequencesforplanetformation
Thepresentedmigrationmapshaveimportantconsequencesfor
whichimpliesthataninterchangebetweenΣ andνwouldnot
G the resulting structures of planetary systems. Lyraetal. (2010)
changetheviscousheating.Theradiativecoolingisgivenby
statedthatplanetsformedin adisc, wouldmigrateandthensit
λc at the zero-migrationradius, which is independentof the plan-
F=−ρ κ ∇ER, (15) etary mass in their case, because they use the torque formula
G R of Paardekooperetal. (2010) that only accounts for the unsat-
wherecisthespeedoflight,λtheflux-limiter(Kley1989),and urated torques. This picture is contradicted by the results pre-
E theradiationenergy.Togetthesamecoolingrate,achange sented here, since we use the torque formula with saturation
R
in ρ has to be balanced by a change in κ . The same applies (Paardekooperetal.2011).
G R
to the absorptionof stellar photons,whichscales with ρ κ . If TheplanetsinLyraetal.(2010)consequentlymoveinwards
G ⋆
thegasdensityisincreasedbyafactors ,theopacity,hencethe withthezero-torqueradiusintimeasthediscdisappears,which
f
metallicity, has to decrease by a factor s . The ratio of Σ /Σ can also be observed from the migration maps presented here.
f G Z
will change, but this will not change the thermal structure of However,theyuseda simple 1D modelfortheirdisc evolution
the disc. This actually means that the surface density of heavy and did not account for stellar heating, which clearly changes
elementsΣ staysconstant. the disc structure and the migration pattern compared to fully
Z
Asaresult,wecanusethedataofthediscspresentedinthe radiativediscs(PaperI).Inparticular,eq.13showsthattheflar-
previoussectionandscalethembyscalingfactors toexchange ing index should be determined accurately for a good estimate
f
Σ and ν. When increasingν by s , the gassurfacedensity Σ of the unsaturated torque, but non-stellar irradiated discs can-
G f G
hasto decreaseby s . Atthesametime,Σ hasto beconstant, not be flared in the outer parts (Paper I). The main difference
f Z
but this means that the metallicity of the disc increases, so in between our model and the Lyraetal. (2010) model originate
scenariolikethat,theopacitiesκ ,κ ,andκ increaseby s as inthemorecomplexdiscstructurepresentedhereandfromthe
R P ⋆ f
well. differenttorqueformulathataccountsforsaturationeffects.
In Fig. 8 we present the migration maps of an M˙ = 1 × Lyraetal. (2010) also stats that in the late evolutionary
10−8M /yr disc, fortwodifferentcasesforexchangingΣ and stages of the disc, starting at ≈ 4.0Myr, the surface density is
⊙ G
ν. Atthetop,theviscosityisdecreasedbyafactorof2,andin solowthatitcannottransferenoughangularmomentumtothe
thebottomplotitisincreasedbyafactorof2.Thismeansthat planet for its orbitalradius to evolve as fast as the equilibrium
thetopdischas0.005metallicityandthebottomone0.02. radius.Thismeansthattheplanetsdetachfromtheevolutionof
These two migration maps significantly differ from each thediscandareleftbehind.Thebeginningoftheirlatestageof
other and from the one presented in fig. 4. In the case of low ≈ 4.0Myrmatchesthesurfacedensityat1AUquitewellinour
viscosity (top in Fig. 8), no outward migration is detected any M˙ = 5 × 10−9M /yr model with a metallicity of 0.01, where
⊙
more. But in the case of high viscosity, a much larger region we no longer observe a zero-torque radius any more for the
9
Bitschetal.:Stellarirradiateddiscsandimplicationsonmigrationofembeddedplanets
torques that are prone to saturation (Paardekooperetal. 2011). once formed,migratesfaster than Jupiter and catchedit in res-
Thiswouldindicatethatplanetarycorescouldstartinwardmi- onance inside their common gap, which causes the migration
gration while the disc is disappearing and cannot safely wait reversal (Masset&Snellgrove 2001; Morbidelli&Crida 2007;
and detachthemselvesfrom the disc as proposedby Lyraetal. Pierens&Nelson 2008). A potential problem that is often in-
(2010).Thereforethiscouldmakethesurvivalofplanetsinthe vokedwiththisscenarioisthataprioriSaturnandJupitershould
outerpartsofthediscinthelatestagesmuchharder.Thisprob- have followed exactly the same growth-migration histories, so
lembecomesevenmoresevereifthemetallicityofthediscdrops that Saturn could never catch up with Jupiter. Our work now
in time (e.g. due to the formation of planetesimals out of dust showsthatthis assumptionis actually nottrue. IfSaturn forms
grains)ascanbeseeninFig.7,whereonlyinwardmigrationis afterJupiter, M˙ willhavedecreasedin thedisc, andthe migra-
observedforM˙ =1×10−8M /yr. tion map will differ. In fact, as times goes by, the maximum
⊙
TheresultsofLyraetal.(2010)couldbeinterpretedinsuch mass for being still trapped in a zero-torque zone decreases,
awaythattheevolutionofthediscnaturallywouldaccountfor so Saturn should start migration inwards at a lower mass than
theoccurrenceofSuper-EarthandMini-Neptuneplanetsinthe Jupiterdid.Possibly,Jupiterneverhadinwardtype-I-migration,
innerpartsofthesolarsystem.However,ifthemetallicityofthe whileSaturnhad,andmayevenhaveexperiencedfasttype-III-
discishighinthelatestagesoftheevolution,thezero-migration migration,whichleadsSaturntocatchupwithJupiter.
radiusfortheseplanetsisat≈ 4AU(seeFig.6),whichismuch Of course this depends on the ratio between the migration
larger than in the model of Lyraetal. (2010), even though we speedofthegiantplanetsandtheactualaccretionalevolutionof
used the torque formula when accounting for saturation here. thedisc.Inaddition,openingofagapbyJupiterwouldchange
Additionally,inthelow-metallicitycase,thezero-torqueradius the thermalstructureof the disc, hencechangingthe migration
for M˙ = 5 × 10−9M /yr using the fully unsaturated torque is map for Saturn. A detailed study of this process is beyond the
⊙
much further out in our simulations (≈ 5AU) compared to ≈ scopeofthispaper,butitisworthnotingthatourresultsquali-
2AUinLyraetal.(2010),whichisclearlyaconsequenceofthe tativelysupporttheconceptoftheGrandTackscenario.
differentdiscmodels. Italso seemsthatit is easier to keep the first coresin discs
Fromoursimulations,weproposetwodifferentwaystoform withhighermetallicity,sincetheminimummassneededforout-
hotNeptunes.Eithertheyformintheouterdisc,wheretheyget wardmigrationisreducedcomparedtolowermetallicitydiscs.
trapped in the outer region of outward migration and are then Thiswouldpotentiallymakeiteasiertokeepthefirstcoresthat
releasedata latestagewheretheymigrateinwardstotheinner formgiantplanets.Because morecorescanpotentiallybekept
edgeofthedisc(low-metallicityscenario).Ortheyassemblein insidethedisc,itismorelikelythatthese corescanformgiant
theinnerpartofthediscfrommigratingembryosthatweretoo planets.Therefore,high-metallicitydiscswouldimplyahigher
small to be trapped in the regions of outward migration. The probabilityofforminggasgiantplanets,whichissupportedby
secondscenariocanhappenforlowandhighmetallicities,since observations.
bothscenariosrequirea minimalplanetarymasstohaltinward Hasegawa&Pudritz(2011)statethatintheicelineanddead
migration. zone,theheattransitionbetweenviscousandstellarheatingcan
Themigrationmapsforametallicityof0.01(Fig.4)canalso actasa planettrap.Ina disc,stellar heatingis dominantin the
be used as a guidelineto discuss the evolution of giantplanets outer parts of the disc and therefore responsible for the flaring
during their formation. The core of a giant planet most likely part of the disc (see Chiang&Goldreich (1997), Paper I). In
forms inside the ice line (r > r ) at the early stages of the the inner parts, viscousheating is dominant.The differentheat
ice
disc. The reason for the core to actually grow to become a gi- sources result in a change of the gradient of temperature be-
antplanetandnotstayahotNeptunemightbefoundinthedif- tween these two different regions, which can lead to a region
ferent formation time of the core itself. A faster formation of ofoutwardmigration,accordingtoHasegawa&Pudritz(2011).
the coreprovidesmoretime to accretegasto becomea gasgi- Additionally, they stated that the ice-line (located at r ) can
ice
ant than a later formationof the core where the planetis stuck functionasaplanettrap.
at Neptune size. We assume that a core is trapped at the outer However,inthesimulationspresentedhere,weseeadiffer-
zero-migrationradiuswhereitcangrowuntilitreachesamass ent picture. The ice line functionsas a region of divergentmi-
of ≈ 40M . Thus, at that stage the planet is in a region be- gration,meaningthatplanetsoutside the ice line (r < r ) mi-
Earth ice
tween≈ 4 and≈ 5AU, dependingonwhichevolutionarystage grateinwards,whileplanetslocatedinsidetheiceline(r > r )
ice
thediscisin(thelower M˙,theclosertothestartheequilibrium migrate outwards. Additionally, one could argue that the disc
radius).Whenthe≈40M massisexceeded,theplanetisre- is locally radialisothermalat the ice line, because melting and
Earth
leasedandstartsmigratingtowardsthecentralstar.Aplanetof resublimating of ice grains keep the temperature constant for
thatsizegrowsrapidlyduetogasaccretion(Pollacketal.1996) a certain distance around the actual ice line (Lyraetal. 2010;
and it starts to open a significant gap inside the disc changing Hasegawa&Pudritz 2011), which is not explicitly taken into
its migrationregimeto thatoftype-II-migration.Owing tothis account in our simulations, but it is slightly reflected by the
migration phase, the giant planet, once fully formed, is signif- smoothed transition between the different branches of opacity
icantly closer to the star than the original location of its core, (Fig.1).Sincethediscislocallyradialisothermal,notempera-
i.e. closerthan∼ 4AU. Theformationofan innercavityin the ture gradientexists aroundthe ice line, so that outward migra-
discbyphoto-evaporationcaneventuallystoptheinwardmigra- tionisnotpossiblethere,andthereforetheicelinecannotactas
tionoftheplanet,assuggestedbyAlexander&Pascucci(2012), apointofconvergentmigration.
preventingitfrombecomingahotJupiter. Theoutwardmigrationinsidetheicelinethenstopsinare-
However,inthecaseofourownSolarSystem,thisisnotsuf- gionofthediscwhereH/rreachesalocalminimum,whichisa
ficient.Amechanismisneededtomovethegiantplanet(Jupiter) convergentmigrationradiusinthedisc.AminimumofH/rcan
outwards and bring it to 5AU. This consideration supports the onlyexist,ifH/rincreasesagain,which,inthiscase,isachieved
GrandTackscenario(Walshetal.2011),inwhichtheformation bystellarirradiation.Inthatsense,achangeforminwardtoout-
of Saturn eventually forces Jupiter to reverse its migration di- wardmigrationiscausedbystellarirradiation,becauseitflares
rectionandmoveoutwards.IntheGrandTackscenario,Saturn, thedisc.Butatthesametimethechangeoftheopacityattheice
10