Eccentricity of radiative discs in close binary-star systems

3 MB·English

by  F. Marzari

Astronomy&Astrophysicsmanuscriptno.ms˙final c ESO2012
(cid:13)
January12,2012
Eccentricity of radiative discs in close binary-star systems
F.Marzari1,C.Baruteau2,3,H.Scholl4 andP.Thebault5
1 DipartimentodiFisica,UniversityofPadova,ViaMarzolo8,35131Padova,Italy
2 DAMTP,UniversityofCambridge,WilberforceRoad,CambridgeCB30WA,UnitedKingdom
3 DepartmentofAstronomyandAstrophysics,UniversityofCalifornia,SantaCruz,CA95064,USA
4 Laboratoire Cassiope´e, Universite´ de Nice Sophia Antipolis, CNRS, Observatoire de la Coˆte d’Azur, B.P. 4229, F-06304 Nice
Cedex,France
5 ObservatoiredeParis,SectiondeMeudon,F-92195MeudonPrincipalCedex,France
2 Preprintonlineversion:January12,2012
1
0 ABSTRACT
2
Context.Discsinbinarieshaveacomplex behavior becauseof theperturbationsof thecompanion star.Planetesimalsgrowth and
n
planetformationinbinary-starsystemsbothdependonthecompanionstarparametersandonthepropertiesofthecircumstellardisc.
a
Aneccentricdiscmaysignificantlyincreasetheimpactvelocityofplanetesimalsandthereforejeopardizetheaccumulationprocess.
J
Aims.Wemodeltheevolutionofdiscsinclosebinariesincludingtheeffectsofself-gravityandadoptingdifferentprescriptionsto
1
modelthedisc’sradiativeproperties.Wefocusonthedynamicalpropertiesandevolutionarytracksofthediscs.
1 Methods.WeusethehydrodynamicalcodeFARGOandweincludeintheenergyequationheatingandcoolingeffects.

Results. Radiative discs have a lower disc eccentricity compared to locally isothermal discs with same temperature profile. Its
]
averaged value isabout 0.05, and it isalmost independent of the eccentricity of thebinary orbit, in contrast tolocally isothermal
P
discmodels.Asaconsequence,wedonotobservetheformationofaninternalellipticallowdensityregionasinlocallyisothermal
E
discmodels.However,thedisceccentricitydependsonthediscmassthroughtheopacities.Akintolocallyisothermaldiscmodels,
h. self-gravityforcesthedisc’slongitudeofpericentertolibrateaboutafixedorientationwithrespecttothebinaryapsidalline(π).
p Conclusions. Thedisc’sradiativepropertiesplayanimportantroleintheevolutionofdiscsinbinaries.Aradiativedischasanoverall
- shapeandinternalstructurethataresignificantlydifferentcomparedtoalocallyisothermaldiscwithsametemperatureprofile.This
o isanimportantfindingbothfordescribingtheevolutionarytrackofthediscduringitsprogressivemassloss,andforplanetformation
r sincetheinternal structureofthediscisrelevant forplanetesimalsgrowthinbinarysystems. Thenon-symmetrical distributionof
t
s massinthesediscscauseslargeeccentricitiesforplanetesimalsthatmayaffecttheirgrowth.
a
[ Keywords.Protoplanetarydisks—Methods:numerical—Planetsandsatellites:formation

1
v
1. Introduction Of particular importance are the initial stages of planetesi-
3
9 mal evolution when the mutual impact velocities must remain
2 Protoplanetary discs are known to exist in pre-main-sequence low to allow mass growth despite the perturbations by the
2 binary star systems through direct imaging (Koerneretal. companionstar. Previousstudies (e.g., Marzari&Scholl 2000;
. The´baultetal. 2006, 2008, 2009; Paardekooperetal. 2008;
1 1993; Stapelfeldtetal. 1998; Rodr´ıguezetal. 1998) and spec-
Xie&Zhou2009;Paardekooper&Leinhardt2010)haveshown
0 tral energy distributions (Ghezetal. 1993; Pratoetal. 2003;
2 Moninetal. 2007). Petr-Gotzensetal. (2010) have recently thatthecombinationofgasdragforceandsecularperturbations
1 claimed that as much as 80% of all binary star systems in by the secondarystar has strongeffects on planetesimalorbits,
: the 1 Myr old Orion Nebula Cluster might have an ac- leading to pericenter alignmentand to an equilibrium distribu-
v
∼ tionfortheeccentricityofsmallplanetesimals.Sincethemagni-
i tive accretion disc. Trillingetal. (2007) showed that the inci-
X dence of debris discs in binary systems is not that different tudeofthegasdragforcedependsontheplanetesimalssize,the
variation of eccentricity and pericenter longitude with particle
r than that for single stars. These observations seem to suggest
a sizemaywellinhibittheformationoflargerbodiesbyexciting
that the ground for planet formation in binary star systems is
largemutualimpactvelocities.
presentevenif,atsubsequentstages,thegravitationalperturba-
tions by the companionstar may negatively impact the growth One question that needs to be addressed is to what extent
process. Bonavita&Desidera (2007); Mugrauer&Neuha¨user the disc’s perturbation due to the companionstar (e.g., Lubow
(2009) have shown that the frequency of planets in binaries is 1991a,b;Kleyetal. 2008)affectsthedistributionsofeccentric-
not very different from that around single stars. However, this ityandlongitudeofpericenterforplanetesimalsofdifferentsize.
frequencycriticallydependsonthe binaryseparationandclose Models in which the planetesimals orbital evolution was com-
binariesappeartobelessfavorableenvironmentsforplanetfor- putedalongwiththetimeevolutionofthedisc(Kley&Nelson
mation.Onlyafewplanetsarepresentlyknownintightsystems 2008;Paardekooperetal.2008)haveshownthatthedisceccen-
(GI86,HD41004andγCephei)anditisexpectedthatwhenthe tricity may strongly influence the dynamicsof small planetesi-
binary separation is in the range 20-100 AU the gravitational malsbyintroducingradialgasdragforcesandnon-axisymmetric
influenceofthecompanionstarmusthaveinfluencedtheplanet componentsinthegravityfieldofthedisc.Thisadditionalper-
formationprocess,affectingeitherthediscevolutionortheplan- turbingforcemayacteitherinfavororagainstplanetesimalac-
etesimalaccumulationprocess(Desidera&Barbieri2007). cretion by affecting the relative impact velocities between the
1
F.Marzarietal.:Radiativediscsinclosebinaries
bodies. An in depth analysis is needed to explore the strength 2. Modeldescription
of these perturbations. This first requires investigating the im-
2.1.Code
pact of the companion star on the disc’s eccentricity. We fo-
cus on relatively close binary star systems with large mass Two-dimensionalhydrodynamicalsimulationswere carriedout
ratios and potentially large eccentricities which, according to with the ADSG version1 of the codeFARGO. The codesolves
Duquennoy&Mayor(1991),populatethepeakofthefrequency the hydrodynamical equations on a polar grid, and it uses an
distribution in our neighborhood. Modelling the response of a upwind transportscheme along with a harmonic, second-order
protoplanetary disc to the gravitational perturbationsby a sec- slopelimiter(vanLeer1977).TheADSGversionoftheFARGO
ondary (companion)star requires by necessity a numericalap- code includes an adiabatic energy equation and a self-gravity
proach,sinceitisdifficulttopredictanalyticallythediscshape module based on fast Fourier transforms (Baruteau&Masset
andevolutionwithtime. 2008).Heatingandcoolingsourcetermshavebeenimplemented
Inapreviouspaper(Marzarietal.(2009),hereafterPaperI) intheenergyequationasdescribedinSect.2.2.Thespecificity
weusedthenumericalcodeFARGOtomodelthetimeevolution oftheFARGOalgorithmistouseachangeofrotatingframeon
ofatwo-dimensional(2D)circumstellardiscinclosebinarystar eachringofthegrid,whichincreasesthetimestepsignificantly.
systems,includingtheeffectsofdiscself-gravity.Wefocusedon Results of simulations are expressed in the following code
massivediscs,withaninitialmassequalto0.04M ,aswellasa units: the mass unitis the mass of the primarystar, Mp, which
largebinary’smassratio(µ = M /M = 0.4).Thi⊙svalueissta- is taken equal to 1M . The length unit is set to 1 AU, and the
tistically themostfrequentamonsgthpebinarysystemsobserved orbitalperiodat1AU⊙is2πtimesthecode’stimeunit.
so far (Duquennoy&Mayor 1991). We found that self-gravity
significantlyaffectstheincreaseinthedisceccentricityinduced
2.2.Physicalmodelandnumericalsetup
bythecompanionperturbations,andalsoinfluencestheorienta-
tionofthediscrelativetothebinaryreferenceframebycausing We adopta 2D disc modelin whichself-gravityandan energy
librationinsteadofcirculation.Wealsosampleddifferentvalues equation are included (unless otherwise stated). The hydrody-
ofthebinary’seccentricity,e ,rangingfrom0to0.6.Themain namicalequationsaresolvedinacylindricalcoordinatesystem
b
findingsofpaperIwerethefollowing: r,ϕ centeredontotheprimarystar,withr [0.5AU 15AU]
{ } ∈ −
andϕ [0,2π].ThegridusedinourcalculationshasN = 256
r
∈
radial zones and N = 512 azimuthal zones, and a logarithmic
s
– self-gravityplaysa significantrolein shapingthedisc.The spacing is used along the radial direction. The frame rotates
dynamical eccentricity e (defined in Sect. 2.3) of the disc with the Keplerian angular velocity at the binary’s semi-major
d
typically ranges from 0.05 to 0.15, depending on the axis, and the indirect terms accounting for the acceleration of
∼ ∼
binary eccentricity, e . It is smaller than in models without the primarydueto the gravityof the secondaryand of the disc
b
selfgravity, areincluded.
– e isinverselyproportionaltoe ,withthecaseofacircular
d b
binary(e =0)beingthemostperturbingconfiguration, Binary parameters— Throughout this study, we adopt as
b
– thediscorientationω (definedinSect.2.3)libratesaround standardmodela binarysystem wherethesecondarystarhasa
d
π,whileitwascirculatingintheabsenceofselfgravity, mass M = 0.4M .Thebinaryisheldonafixedeccentricorbit
s
– an eccentric low-density region develops in the inner disc with semi-major⊙axis a = 30 AU and eccentricity e = 0.4,
b b
parts because of the large eccentricity and aligned pericen- correspondingtoanorbitalperiod 134yr.
≈
tersofthegasstreamlinesthere.
Energyequation—Sincethepurposeofthisworkistoexamine
the impact of an energy equation on the disc’s response to the
TheresultsofpaperIhavebeenobtainedassumingalocally periodicpassagesofaclosecompanion,wecarriedsimulations
isothermal equation of state for the gas, wherein the initial ra- includingeitheranenergyequation,oralocallyisothermalequa-
dialprofileofthetemperatureremainsconstantintime,itsvalue tionsuchthattheinitialtemperatureprofileremainsconstantin
being set by the choice for the disc aspect ratio h = H/r, with time.Inallcases,thediscverifiestheidealgaslaw,
H the pressure scale height. This approximation is well suited
in disc regions where radiation cooling is efficient and the gas p= ΣT/µ, (1)
R
is optically thin. However, in particular in the initial stages of
wherepandT denotethevertically-integratedpressureandtem-
theirevolution,discsmaybeopticallythickandamoredetailed
perature,respectively,Σisthemasssurfacedensity, istheideal
treatmentoftheenergybalanceisrequired.Inaddition,whenit
R
gasconstant,andµisthemeanmolecularweight,takenequalto
passesatitspericenter,thesecondarystartriggersspiralshocks
2.35.Whenincluded,theenergyequationtakestheform:
that may generate local strong shock and compressional heat-
ing,whichmayviolatethelocalisothermalapproximation.The
∂e
propagationof these shock waves may also be significantly al- ∂t +∇·(ev)=−p∇·v+Q+visc−Q−cool+λe∇2log(p/Σγ) (2)
teredinradiativediscs.Inthispaper,wefocusonhowthedisc
eccentricity and orientation depend on the disc radiative prop- wheree= p/(γ 1)isthethermalenergydensity,γistheadi-
−
erties.Ourapproachisthefollowing.Wefirstexaminehowour abatic index, and v denotes the gas velocity. We take γ = 1.4
previousresultsonlocallyisothermaldiscsdependonthechoice throughoutthisstudy.InEq.(2),Q+ denotestheviscousheat-
visc
forthe(fixed)temperatureprofileofthedisc.Wethenconsider ing.Weusebothaconstantshearkinematicviscosity,ν = 10 5
−
discmodelswitharadiativeenergyequation,forwhichwefind incodeunits,andavonNeumann-Richtmyerartificialbulkvis-
thatthe averageddisc’s eccentricityinits innerpartsissmaller cosity, as described in Stone&Norman (1992), where the co-
comparedtolocallyisothermaldiscmodelswithsimilartemper- efficient C is taken equal to 1.4 (C measures the number of
2 2
atureprofile.Wefinallydiscusstheimpactofourresultsinterms
ofplanetesimalsdynamics. 1 See:http://fargo.in2p3.fr
2
F.Marzarietal.:Radiativediscsinclosebinaries
zonesoverwhichtheartificialviscosityspreadsoutshocks).The and
coolingsourceterminEq.(2),Q ,istakenequalto2σ T4 ,
−cool SB eff
wfehcetirveeσteSmBpisertahteurSete(Hfaunb-Benoylt1zm99a0n)n, constant,andTeff istheef- ̟d = Md−1×" ̟(r,ϕ)Σ(r,ϕ)rdrdϕ (6)
T4 =T4/τ , (3) where Md denotes the disc’s mass, and e and ̟ are the eccen-
eff eff tricity and pericenter longitude of each grid cell, respectively,
witheffectiveopticaldepth assuming that the local position and velocity vectors uniquely
definea2-bodyKeplerianorbit.
3τ √3 1
τ = + + . (4)
eff 8 4 4τ
3. Locallyisothermaldiscmodels
The vertical optical depth, τ, is approximated as τ = κΣ/2,
Beforeexaminingtheimpactofanenergyequationonthedisc’s
where for the Rosseland mean opacity, κ, the formulae in
response to the tidal perturbations by the secondary star, we
Bell&Lin (1994) are adopted. Following Paardekooperetal.
adopt in this section the simpler case for a locally isothermal
(2011), we also model thermal diffusion as diffusion of the
gas entropy, s, defined as s = (γ 1) 1log(p/Σγ). This equationofstate,wheretheinitial(axisymmetric)profileofthe
−
R − disc temperature remains fixed. This first step will help us an-
corresponds to the last term in the right-hand side of Eq. (2),
alyze the more complex situation of a disc whose temperature
whereλis a constantthermaldiffusioncoefficient.Throughout
thisstudy,weadoptλ=10 6incodeunits. evolvesin time due to radiativecooling and varioussourcesof
−
heating,includingthatarisingfromtheshockwavesinducedby
thesecondary.
Initial conditions— The disc is initially axisymmetric and the
angular frequency Ω(r) about the primary star is computed
taking into account the radial acceleration due to the pressure 3.1.Disc’seccentricityandsurfacedensityprofiles
gradient, the gravitational acceleration due to the primary
star only and the self–gravity of the disk. The initial gas Wecarriedoutaseriesofsimulationsusingarangeofvaluesfor
surface density, Σ , is set as Σ (r) = Σ (1AU) (r/1AU) 1/2 thedisc’saspectratioat1AU, h(1AU).Thiscomestovarying
0 0 0 −
with Σ (1AU) = 2.5 10 4 in code units. ×Assuming the thedisc’stemperatureatthesamelocation.Otherdiscandbinary
0 −
mass of the primary sta×r is 1M , this corresponds to setting parametersareasdescribedinSect.2.2.Weconsidertwodiffer-
the disc surface density at 1 A⊙U to 2.2 103 g cm 2. entvaluesforthesecondary-to-primarymassratio:q=0.4(our
−
Beyond 11 AU, Σ (r) is smoothly red≈uced to×a floor value, fiducialvalue),andq=0.1.Resultsofsimulationswithq=0.4
0
Σ = 10 9 = 4 10 6Σ (1AU), using a Gaussian function areshownintheupperpanelsinFig.1,andthosewithq = 0.1
floor − − 0
with standard devia×tion equalto 0.4 AU. To preventnumerical inthelowerpanels.Azimuthally-andtime-averagedprofilesof
instabilities caused by low density values or by steep density the disc’s eccentricity and surface density are displayed in the
gradients near the grid’s outer edge, the gas density in each left and right columns of Fig. 1, respectively. Profiles are dis-
grid cell is reset to Σ whenever it becomes smaller than playedafter 60orbitsofthesecondary,andtimeaveragingis
floor ∼
this floor value (Kleyetal. 2008, Paper I). Similarly, we adopt doneover5orbits.Bycheckingthetimeevolutionofthedisc’s
a floor value for the thermal energy density, e = 10 18. eccentricity profile, we find that a steady state is reached after
floor −
Our choice for this parameter is conservative, and we checked typically20binaryrevolutionswithq = 0.4forallvaluesofh,
that with a larger e , the modelled disc behaved similarly. and after 40 revolutions for q = 0.1. The value of h(1AU) is
floor
The initial gas temperature, T , is taken proportional to r 1: indicatedineachpanel(andsimplydenotedbyh).
0 −
T = T (1AU) (r/1AU) 1, where the value for T (1AU), We are primarily interested in the disc’s averaged eccen-
0 0 − 0
whichistakenas×afreeparameterinourstudy,willbespecified tricity in its inner parts, below the truncation radius located at
below. Writing the pressure scale height H as H = c Ω 1, r 5 6 AU for q = 0.4, where planet formation is more
S −K ∼ −
where c denotes the sound speed and Ω the Keplerian likely to occur. As shown in Marzarietal. (2009), the trunca-
S K
angular frequency, the disc temperature can be conveniently tion radius of the disc closely matches the critical limit for or-
related to the disc’s aspect ratio, h = H/r. In our code units, bitalstabilityduetothesecularperturbationsofthecompanion
where /µ=1,thisyieldsT (1AU)=630K h(1AU)/0.05 2. star (Holman&Wiegert1999). From Fig. 1, it is clear that the
0
R ×{ } disceccentricityincreasesfromr 4AUdownwards,peaksat
Boundary conditions— By default, an outflow zero-gradient r.1AU,thendecreasestowardst∼helocationofthegrid’sinner
boundaryconditionisadoptedatboththegrid’sinnerandouter edge, where the boundary condition imposes zero eccentricity.
edges. The azimuthalvelocityis set to its initial, axisymmetric Interestingly, we see that the averaged peak eccentricity of the
value. No mass therefore flows back into the system, and the disc(reachedatr . 1AU)increaseswithhuptoh = 0.05and
diskmassdeclineswithtime.Theimpactoftheinnerboundary decreasesbeyondthisvalue.Thesamebehaviorisobtainedwith
conditiononourresultsisexaminedinSect.4.3. both values of q. This behavior may be interpreted as follows.
The disc’s density perturbationdue to the secondary decreases
with increasing disc aspect ratio (h or, equivalently, increas-
2.3.Notations ing sound speed). In the limit of large aspect ratios, the disc’s
eccentricity thus decreases with increasing h. As h decreases,
We will make use of the following quantities. We denote
thedisc’sperturbeddensityincreases,butshockwavesbecome
by e and ̟ the disc’s eccentricity and perihelion longi-
d d
moretightlywound,andhavetotravelalongerdistancebefore
tude, respectively. They are defined as in Kleyetal. (2008);
reachingthedisc’sinnerparts.Beingthenmorepronetoviscous
Pierens&Nelson(2007);Marzarietal.(2009):
damping,shockwavesbecomelessandlessefficientatdeposit-
ing angular momentum in the disc’s inner parts, where the ec-
e = M 1 e(r,ϕ)Σ(r,ϕ)rdrdϕ (5)
d d− ×" centricityremainssmall.Thisexplainswhywhendecreasingh,
3
F.Marzarietal.:Radiativediscsinclosebinaries
theaveragedeccentricityofthedisc’sinnerpartsfirstreachesa 4. Radiativediscsmodels
maximumandthendecreases.
In this section, we describe the results of our hydrodynamical
We cansee in Fig. 1 thatthe radialdependenceof the den-
simulationsthatincludeanenergyequation.Ouraimistoassess
sity and eccentricity profiles reasonably match each other, as
theimpactoftheenergyequationonthedisc’saveragedeccen-
expected. The general trend is that the larger the disc’s eccen-
tricity.
tricity profile, the smaller the surface density profile. Note that
this is more visible for q = 0.4, where the density’s perturba-
tionbythe secondarytakeslargervaluesthanforq = 0.1.The 4.1.Disceccentricity
significantdecreaseinthedisc’sdensityprofileinthediscinner
parts(typicallybelow 2AU)takestheformofanellipticinner
∼
hole,whichhasbeenobservedinanumberofpreviousnumeri-
0.12
calstudies(e.g.,Kleyetal.2008,PaperI).
Isothermal
Our results concerningthe dependenceof the disc’s eccen- Radiative
tricityanddensityprofileswithhdifferfromthoseofKleyetal. 0.1
(2008), but this dependson intrinsic differencesin the models.
Inparticular,weconsidermoremassivediscsandweincludethe 0.08
effectsofself-gravitywhichprovedtobeefficientinaffectingthe city
discevolution. ntri 0.06
e
c
c
E
0.04
3.2.Fourieranalysisofthedensitydistribution
Toprovidesomeinsightintothedependenceofthedisc’seccen- 0.02
tricityanddensityprofileswith varyingtemperature,weexam-
ine in this paragraphthe Fourier componentsof the disc’s sur-
0
face density.Fig. 3 displaysthe instantaneousamplitudeof the 0 20 40 60 80 100
Fouriermodecoefficientswithazimuthalwavenumber1 m Binary revolutions
≤ ≤
5. Resultsare shownat2550yr,thatis after 20orbitsof the
≈
secondary,whenthediscistruncatedatabout7-8AU.Wecom- 6 Isothermal
parethe profilesobtainedforpreviousseries oflocallyisother- Radiative
mal disc models for q = 0.4, with h = 0.02 (left panel) and 5
h = 0.06(rightpanel).Thesecondarystar ishalf-waybetween
thepericenterandapocenteratthisparticularpointintime. d) 4
a
efficFieronmtsdFeigc.re3a,sietsiswclietharinthcaretathsienagmmp,liatunddethoaftFtohuermier=m1odmeocdoe- nter (r 3
prevails.Therunwithh=0.06displayslargermodeamplitudes ce
throughoutthe disc’s inner parts, whichaccountsfor the larger eri
P 2
disceccentricityobtainedinthiscase.Also,thelargeamplitude
ofthem = 1modebelowr 1forh = 0.06isdirectlyassoci-
≈ 1
atedtothepresenceofanellipticinnerhole,whosepresencehas
alreadybeenpointedoutintheupperpanelsinFig.1.
0
0 20 40 60 80 100
Binary revolutions
350
Fig.5. Timeevolutionoftheaverageddisceccentricitye (up-
d
perpanel)andpericenterlongitude̟ (lowerpanel)forthelo-
300 d
callyisothermalandradiativerunsofSect.4.1.
250
g)
e
d Asafirststep,wecomparearadiativeandalocallyisother-
er ( 200 malmodelwithsameinitialaspectratio,h=0.05.Allotherdisc
nt
ce 150 parametersareasdescribedinSect.2.Fig.4showscontoursof
eri the disc’ssurfacedensityat 90binaryorbitsobtainedwith two
P
100 valuesofthe binary’seccentricity:eb = 0.4(upperpanels)and
e = 0(lowerpanels).AsalreadypointedoutinSect.3,eccen-
b
50 tric streamlines (m = 1 density mode) have different values of
thepericenterdependingontheradialdistance,and,asaconse-
0 quence,theycombineintoapatternofspiralstructure.However,
0 1 2 3 4 5 6 7
the discs computedwith the radiative model (left-handedplots
R (AU)
inFig.4) appearsmootherandmoresymmetricthanthecorre-
Fig.2. Azimuthally-averagedprofile of the disc’s pericenter in spondinglocallyisothermaldiscs,and,inparticular,theydonot
thelocallyisothermaldiscmodelwithh = 0.05andq = 0.4.It feature any elliptic hole near the inner edge, independently of
iscomputedafter100orbitsofthesecondarystar. e .Thespiraldensitywavesarelessstrongintheradiativecase
b
andthisispossiblyrelatedtotheabsenceofalowdensityregion
closetothestar.Thedisceccentricitye andperihelionlongitude
d
4
F.Marzarietal.:Radiativediscsinclosebinaries
0.6 -3
h = 0.02
h = 0.03
0.5 h = 0.04
h = 0.05 -4
h = 0.06
0.4 h = 0.08
city Σ) -5
ntri 0.3 10(
e g
cc Lo -6 h = 0.02
E h = 0.03
0.2 h = 0.04
h = 0.05
-7 h = 0.06
0.1 h = 0.07
h = 0.08
Unperturbed
0 -8
0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8
R (AU) R (AU)
0.6 -3
h = 0.02
h = 0.03
0.5 h = 0.04
h = 0.05 -4
h = 0.06
0.4 h = 0.08
city Σ) -5
ntri 0.3 10(
e g
cc Lo -6 h = 0.02
E h = 0.03
0.2 h = 0.04
h = 0.05
-7 h = 0.06
0.1 h = 0.07
h = 0.08
Unperturbed
0 -8
0 2 4 6 8 10 0 2 4 6 8 10
R (AU) R (AU)
Fig.1.Azimuthally-andtime-averagedprofilesofthedisc’seccentricity(leftcolumn)andsurfacedensity(Σ,rightcolumn)obtained
withtheseriesoflocallyisothermaldiscmodelsdescribedinSect.3.Inthedensityprofileplotstheinitial,unperturbedprofileis
shownasreference.Resultsareshownat60orbitsofthesecondary,andtime-averagingisdoneover5orbits.Severalvaluesofthe
disc’s aspectratio at 1 AU are consideredandtwo secondary-to-primarymass ratios: q = 0.4 (upperpanels) and q = 0.1 (lower
panels).
areshownasa functionoftime in Fig.5 andtheyconfirmthat profileoftheaboveradiativerun.Inthe secondmodel,theini-
theradiativecasehas,onaverage,alowereccentricityevenifit tial temperatureis chosento givethe same soundspeed profile
takesmoretimetoreachasteadystate.Theperihelionlibration, as the time-averaged one in the radiative run. In Sect. 3.1, we
observed also for the radiative case, strongly suggests that this pointedoutthatthedisceccentricityinlocallyisothermalmod-
behaviouris solely due to the disc self-gravity and that it does elsstronglydependsontheaspectratiohand,asaconsequence,
notdependontheenergyequation.Byinspectingtheradialperi- on the temperature profile (see Fig. 1). In Fig. 6, bottom plot,
centerprofile,evenintheradiativecasetheazimuthallyaveraged we compare the temperature profile of the locally isothermal
pericenterofthegasstreamlineschangeswithradialdistance(as model with h = 0.05 to the time-averaged temperature profile
isillustratedinFig2foralocallyisothermalmodel).Thisvari- of the radiative model. The temperature in the radiative run is
ation,inadditiontocausingspiralwaves,leadstoanasymmet- significantlylarger,implyingthatweshouldcomparetheradia-
ricdistributionofmassandthentoanon-symmetricdiscgravity tive modelto a locally isothermalmodelwith h 0.07.In the
∼
field.Thisiscriticalforplanetesimalsembeddedinthediscsince toppanelofFig.6,wedisplaytheeccentricityanddensitypro-
thenon-homogeneousdiscforcesnon-radialcomponentsonthe files in all four models: (i) the radiative model, (ii) the locally
gravity field felt by planetesimals significantly perturbingtheir isothermalmodelwithsameinitialaspectratioasintheradiative
orbits. These perturbationsare comparablein magnitudeto the model,(iii)thelocallyisothermalmodelwithsametemperature
secular effects of the companionstar but are irregularandmay profile as in the radiative case, and (iv) the locally isothermal
thencauselargerchangesintheplanetesimalsorbitalelements. modelwithsamesoundspeedprofileasintheradiativeone.The
Theyareindeedanindirecteffectofthecompaniongravitybut disceccentricityinmodels(iii)and(iv)issmallerthaninmodel
for planetesimals they represent an independentsource of per- (ii),asexpectedfromtheresultsshowninFig.1(top-leftpanel).
turbation. Consequently,thediscsurfacedensityremainslargerinmodels
To getfurtherinsightintothedifferentdensitiesandeccen- (iii) and (iv) compared to in model (ii). Still, models (iii) and
tricities with and without an energy equation, we carried out (iv)do notmatchthe low eccentricityof the radiativecase and
twoadditionallocallyisothermaldiscmodels.Inthefirstmodel, itssmoothdensitydistributionclosetothestar.Thestrengthof
the initial temperature is set to the time-averaged temperature the m = 1 Fourier mode in the radiative run is smaller in the
5
F.Marzarietal.:Radiativediscsinclosebinaries
0.7 0.7
m=1 m=1
m=2 m=2
m=3 m=3
0.6 mm==45 0.6 mm==45
0.5 0.5
0.4 0.4
m m
Σ Σ
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8
r (AU) r (AU)
Fig.3.InstantaneousFouriercomponentsofthedisc’ssurfacedensityforthelocallyisothermaldiscmodelsofSect.3withh=0.02
(leftpanel)andh=0.06(rightpanel),respectively.Resultsareshownatt=2550yr(thatis,afterabout20orbitsofthesecondary),
andareobtainedwithq = 0.4.TheazimuthalwavenumberforeachFouriercomponentisindicatedinthetop-rightcornerofeach
panel.
disc’sinnerparts,andthisdifferencedependsontheformofthe the disc to be no less than that at 6 AU from the star (that is,
energyequation.InFig.7,wedisplaydensitycontoursobtained about150yr). To keep the temperatureprofile as close as pos-
withtheradiativerun(leftpanel),andwiththelocallyisothermal sibletothatofthestandardradiativemodel,andthustoprevent
modelwithsametemperatureprofile(rightpanel).Intheradia- discrepanciesrelated to differenttemperatureprofiles, we limit
tivemodel,thewaveperturbationsarelessstrongandthewaves the viscousheatingtimescale by a similar amount.In a second
appeartobemoredamped.Beingweakerintheinnerdiscparts run, we increase the cooling rate throughoutthe entire disc by
in the radiative case, waves should be less efficient at deposit- afactorof10,andtheviscousheatingrateisincreasedaccord-
ingangularmomentumthere.Thiswouldexplainwhyradiative ingly.Bothrestartsimulationswererunfor50additionalbinary
discs are less eccentric in the regions close to the central star, orbits,overwhichthedisc’stemperatureprofiledoesnotevolve
andwhythesurfacedensityismorehomogeneous. significantly,exceptwithin1AUfromthecentralstar.
Fig.8comparesthedisctime-averagedeccentricityandtem-
peratureprofilesofthestandardradiativemodelwiththoseofthe
4.2.Radiativedampingofwaves
additionalmodelswithlongerandshortercoolingtimes.Thein-
Themainquestionarisingfrompreviousresultsiswhatcauses nerdisc’seccentricityissignificantlyincreasedwithlongercool-
astrongerdampingofdensitywavesintheradiativemodelthan ing timescale, exceeding 0.1 almost uniformly from R = 2
∼
in the locally isothermal one, while the same (time-averaged) AU to R = 6 AU. Thisvalueis in goodagreementwith thatof
sound speed profile is adopted in both models. In Fig. 6, the the locally isothermal run with same sound speed profile (see
radiativemodelshowsindeeda significantlylower eccentricity Fig.6).However,closetotheinneredgethedisceccentricityis
profile.Assumingthedisceccentricityisrelatedtothestrength stillsmallbutthisispossiblyduetothefactthatwecannotmain-
of the density waves induced by the binary gravitational per- tainconstantthetemperatureprofilewithin1AUofthecentral
turbations,we have shownin Fig. 7 thatspiralwaves are more starinthemodelwithreducedcooling.Intherestartsimulation
dampedintheradiativecase. withacoolingrate10timeslarger,thedisceccentricityisabout
There are two potentialsources of wave dampingthat may halfthatofthestandardradiativemodel,whilethecorrespond-
help understandwhy the disc eccentricityremainslower in the ingtemperatureprofileshardlydiffer.
radiativecase:shockdamping(Goodman&Rafikov2001)and TheresultsshowninFig.8suggestthatradiativedampingis
radiativedamping(Cassen&Woolum1996).Toexploretheef- responsibleforthelimitedgrowthinthedisceccentricitycom-
ficiency of the shock damping mechanism, we examined the paredtotosimilarlocallyisothermaldiscmodels.Thismecha-
vortensity distribution in the disc, since it experiences a jump nism,likeself-gravity,helpsmaintainingthedisceccentricityto
at shocks, the magnitude of which depends on the strength of amildvalue.
theshock.Wedidnotobservesignificantdifferencesinthedisc
vortensitydistributionbetweenthetwomodelsandthuswebe-
4.3.Dependenceontheboundaryconditions
lievethatshockdampingdoesnotsignificantlycontributetothe
smootherbehaviouroftheradiativedisc. We examineinthisparagraphthedependenceofourresultson
Theenergylossbyradiationisanadditionalpossiblesource the choice for the boundary condition at the grid’s inner edge.
ofwavedamping(Cassen&Woolum1996).Wavepropagation For this purpose, we compare the results of simulations using
throughadiabaticcompressionsandexpansionsmaybedamped our fiducial boundary condition (zero-gradient outflow bound-
by radiative losses which, in our model, are controlled by the arycondition,see Sect.4.1)withthoseoftwoadditionalsimu-
coolingtermQ−cool.Totestthishypothesis,werestartedourstan- lations:
dard radiative simulation (e = 0.4, M = 0.4M and a = 30
b s b
AU)after165binaryrevolutionsadoptingdiffere⊙nttwocooling - One using a viscous outflow boundary condition
prescriptions.Inafirstrun,welimitthecoolingtimethroughout (Kley&Nelson 2008). It is very similar to our zero-
6
F.Marzarietal.:Radiativediscsinclosebinaries
-012345.5 -5.0 -4.5 -4.0 -3.5 -3.0 -012345.5 -5.0 -4.5 -4.0 -3.5 -3.0
-5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0
Log (Surface density) at t = 90Torb Log (Surface density) at t = 90Torb
5 5
U] U]
A 0 A 0
y [ y [
-5 -5
-5 0 5 -5 0 5
x [AU] x [AU]
-012345.5 -5.0 -4.5 -4.0 -3.5 -3.0 -012345.5 -5.0 -4.5 -4.0 -3.5 -3.0
-5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0
Log (Surface density) at t = 90Torb Log (Surface density) at t = 90Torb
10 10
5 5
U] U]
A 0 A 0
y [ y [
-5 -5
-10 -10
-10 -5 0 5 10 -10 -5 0 5 10
x [AU] x [AU]
Fig.4. Contoursof the gas density after 90 binary revolutions.Results with an energyequation are shown in the left part of this
figure,andthose with a locallyisothermalequationofstate with fixedtemperatureprofile are shownin the rightpart. Theupper
plotsarefore =0.4andthelowerplotsfore =0.0.Thesameinitialaspectratio(h=0.05)isusedinthesesimulations.
b b
gradient outflow boundary condition, except that the disc at the grid’s inner edge is no longer forced to remain
radial velocity at the inner boundary is set to the local circular.
(azimuthally-averaged) viscous inflow velocity of a disc
The results of this comparison for locally isothermal runs
in equilibrium with locally uniform surface density profile
with same temperature profile as in the radiative model are
( 3ν/2R ). Theazimuthalvelocityisalso setto its initial
− min showninFig.9.Theprofilesareingoodagreementforallthree
axisymmetricvalue,
different boundary conditions. Note that the surface density in
- Another simulation also using our standardoutflow bound-
the innerdisc fortheviscousboundaryconditionis largerthan
arycondition,butwheretheazimuthalvelocityatthe inner
with the two other boundaries. This is expected, since the im-
boundaryisextrapolatedfromthatinthefirstactiveringwith
posedviscousdriftvelocityissmallerthattheradialvelocityset
ar 1/2 law.Incontrasttopreviousboundaryconditions,the
− bythepropagationofthespiralwavesinducedbythesecondary.
Asimilargoodagreementisobtainedwitharadiativemodel,as
7
F.Marzarietal.:Radiativediscsinclosebinaries
0.5 -3
Radiative Radiative
Isoth. with h=0.05 Isoth. with h=0.05
Isoth. with same T profile -3.5 Isoth. with same T profile
0.4 Isoth. with same c profile Isoth. with same c profile
s s
-4
city 0.3 Σ)
ntri 10( -4.5
e g
c o
c 0.2 L
E
-5
0.1
-5.5
0 -6
0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
R (AU) R (AU)
1400
Radiative
Isoth. with h=0.05
1200
1000
K)
ure ( 800
at
per 600
m
e
T
400
200
0
0 1 2 3 4 5 6 7
R (AU)
Fig.6.Time-andazimuthally-averageddisceccentricity,densityandtemperatureprofilesforourradiativediscmodelandvarious
locallyisothermaldiscmodels:(i)onewiththesameinitialaspectratioasintheradiativemodel(h=0.05),(ii)onewiththesame
temperatureprofileasintheradiativemodel,and(iii)anotherwithsamesoundspeedprofileasintheradiativemodel.
-012345.0 -4.5 -4.0 -3.5 -012345.0 -4.5 -4.0 -3.5
-5.0 -4.5 -4.0 -3.5 -5.0 -4.5 -4.0 -3.5
Log (Surface density) at t = 35.86Torb Log (Surface density) at t = 35.86Torb
6 6
4 4
2 2
y 0 y 0
-2 -2
-4 -4
-6 -6
-6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6
x x
Fig.7.Contoursofthediscsurfacedensityfortheradiativemodel(leftpanel),andthelocallyisothermalmodelwithsameimposed
temperatureprofileasintheradiativecase(rightpanel).Intheradiativecase,spiralwavesappearsmootherandmoredamped.
8
F.Marzarietal.:Radiativediscsinclosebinaries
0.25 10-3
0.20
10-4
y
t
ntricity 0.15 e densi 10-5
e c
c 0.10 a
c f
E ur
S 10-6
0.05 Open BC Open BC
Viscous BC Viscous BC
Open BC+vθ extrapolated Open BC+vθ extrapolated
0.00 10-7
0 2 4 6 8 0 2 4 6 8
R [AU] R [AU]
Fig.9. Eccentricityandsurfacedensityprofiles,time-averagedbetween35and40binaryorbits,obtainedwiththelocallyisothermal
disc model with imposed temperature profile, using three different boundary conditions at the inner edge: (i) our standard zero-
gradientoutflowboundary(labelledasopen,solidcurve),(ii)aviscousoutflowboundary(labelledasviscous,dash-dottedcurve,
seetext),and(iii)ourstandardzero-gradientoutflowboundarywiththeazimuthalvelocityattheinneredgeextrapolatedfromthat
inthefirstactivering(dashedcurve).
0.08 10-3
0.06 10-4
y
t
ntricity 0.04 e densi 10-5
e c
c a
c f
E ur
0.02 S 10-6
Open BC Open BC
Viscous BC Viscous BC
Open BC+vθ extrapolated Open BC+vθ extrapolated
0.00 10-7
0 2 4 6 8 0 2 4 6 8
R [AU] R [AU]
Fig.10. SameasinFig.9,butforourstandardradiativediscmodel.
illustratedinFig.10.Theseresultsconfirmtherobustnessofour ing from 0 to 0.6. We note that the eccentricity profiles signif-
resultsagainstthechoicefortheinnerboundarycondition. icantlydiffer evenif the positionof the star with respectto the
initial reference frame is the same. All radiative runs show an
increasing eccentricity profile towards large radii, whereas lo-
4.4.Dependenceonthebinaryeccentricity cally isothermal runs feature a peak in the eccentricity in the
inner disc parts for most values of h. It should be pointed out
InpaperI,weshowedthatthedisceccentricityeddecreaseswith thatthe abovecomparisonis doneby adoptingthe same initial
increasing binary eccentricity, the circular case (eb = 0) being aspectratio forradiativeand locallyisothermaldisc models.A
the mostperturbingconfiguration.We interpretedthis resultas differentapproachwouldbetocomparediscmodelswithsame
being due to the larger size of the disc for lower values of eb, temperatureorsoundspeedprofiles,asisdoneinFigs.7and6.
and to the consequent larger number of resonant perturbations But, when comparing the temperature profiles of the radiative
whichmayaffectit.Theeccentricityofradiativediscs,however, discs, we find that they correspond to locally isothermal discs
seems to be rather insensitive to eb. The outcome of the simu- withh 0.06 0.07,whichanywayfeaturealargereccentricity
lationsatdifferenteb isshownintheleftplotofFig.11,where inthed∼iscinn−erparts.
wecomparethevaluesofe obtainedforlocallyisothermaland
d
radiativediscs.Incontrasttolocallyisothermaldiscs,theaver- In the rightpanelof Fig. 11 we also comparethe valuesof
agede ofradiativediscsisalmostconstantaround0.05,andit theoutersemi-majoraxisa oftheellipsebestfittingtheouter
d d
does not show a significant dependence with e . However, the edge of the disc (at a density level of Σ = 10 6, code units).
b −
azimuthally-averagedprofile of the disc’s eccentricity over the The radiative discs extend slightly farther out compared to the
radiativediscdoesdependone .InFig.12,theradialprofileof corresponding locally isothermal discs, except for our fiducial
b
the disc’s eccentricity is shown for different values of e rang- binary’seccentricity(e = 0.4),wherethedisc’ssizeisapprox-
b b
9
F.Marzarietal.:Radiativediscsinclosebinaries
0.14 12
e (isoth) e (isoth)
d d
0.12 ed (rad) ed (rad)
10
y 0.1
ntricit 0.08 AU) 8
ce a (
ec 0.06
6
0.04
0.02 4
0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6
e e
B B
Fig.11.Disceccentricity(leftplot)fordifferentvaluesofthebinaryeccentricitye .Theaverageofe isdisplayedafter70binary
b d
revolutions, and time-averaging is done over 20 revolutions. For comparison, we also depict the corresponding values of e for
d
isothermaldiscswithsamebinaryconfigurations(Marzarietal.2009).Intherightplot,weshowthedisc’ssizefordifferentvalues
ofe ,morespecificallythesemi-majoraxisoftheellipsethatbestfitsthediscdensityatalevelΣ=10 6(codeunits).
b −
imately equal to its isothermal counterpart. The mass loss rate whylessmassiveradiativediscshavedifferentdiskeccentricity.
forradiativediscs,fordifferentvaluesofe ,isshowninFig.13. Even self–gravity will play a minor role for less massive discs
b
Itiscomputedwithalinearfitofthetimeevolutionofthedisc and in Paper I we showed that indeed self–gravity is effective
masswhenasteadystateisreached,soitisindependentfromthe in reducing the disc eccentricity. This, however, is not the full
initialdisctruncationprocess.Itshowsanalmostlineardepen- story and the different energy equation also plays a significant
dencewithe anditisrelatedtothestrongperturbationeffects role. The disc eccentricity radial profile of radiative discs, also
b
onthediscasthesecondarystarpassesatpericenter.Theinspec- the less massive ones (see Fig. 14), does not peak close to the
tionoftheamountofmassstrippedvs.timerevealsthatalmost star like for the isothermal discs (Fig. 1). The value of e for
d
allthemasslossoccursduringandafterthepericenterpassageof radiativediscs growsfor largervaluesof R as predictedby the
thecompanionwhilea loweramountisconstantlylostthrough analytical theory of Paardekooperetal. (2008). This behaviour
the inner borderdue to the disc eccentricity and viscosity. The isfurtherillustratedbyexaminingtheFouriercomponentsofthe
absolutevalueofthemasslossislarge,anditpredictsareduc- normalizedsurface density. Fig. 16 comparesthe Fouriercom-
tioninthediscmassbyafactoroftwoinabout3 104 yrsfor ponents for both our standard disc with Σ = Σ , and the
0 MMSN
×
our standard case (e = 0.4). However, this value depends on disc model with Σ = 0.1Σ . The m = 1 and m = 2 com-
b 0 MMSN
thediscmass,andasimulationwiththesamebinaryparameters ponents are much stronger for the less massive disc, and they
butan initialdisc mass 10 times smaller givesa mass loss rate accountfor the overallhigher disc eccentricity. In the standard
3.7 10 9 M /yr. case, the high disc density is able to damp efficiently the two
−
≈ × ⊙ Fouriercomponentsofthebinaryperturbationsastheymoveto-
wardthediscinnerparts.Ontheotherhand,thesecomponents
4.5.Dependenceonthediscmass
propagate further in for less massive discs, leading to a larger
InPaperI,wefoundthatforlocallyisothermaldiscs,adecrease valuefored.However,alsointhelessmassivedisctheydonot
in the disc mass had no significant impact on the disc eccen- reachtheinnerpartofthedisc,likeinisothermaldiscs,prevent-
tricity.Belowweshowthatradiativediscsbehavedifferently.In ingtheformationofaneccentriclow-densityregionclosetothe
Fig. 14, we compare the disc eccentricity in the nominal case, centralstar.
wheretheinitialdensityat1AUisthatoftheMMSN,tothatof Accordingtotheresultspresentedsofar,isothermalandra-
adisc initially2and10timesless massive.Thedisceccentric- diativediscsbehavedifferently.Massiveisothermaldiscsinlow
itye inasteadystateismuchlargerforthelessmassivediscs. eccentricitybinariesareexpectedtobeeccentricforreasonable
d
Thisisfurtherillustratedbythecontoursofthediscdensityfor values of h while radiative discs always have low eccentricity.
Σ = 0.1Σ in Fig. 15. The disc appears smaller and very For isothermal discs there is a weak dependence of the disc
0 MMSN
eccentric,particularlyinitsouterregions. eccentricity on the disc mass (see Paper I) while, for radiative
It is difficult to interpret this outcome on the basis of the discs, this dependence is strong with less massive discs being
moreeccentric.Hotisothermaldiscsdevelopaninnereccentric
isothermal simulations. Less massive radiative discs, like that
holewhileradiativediscsaresmoothalsoattheinneredge.
shown in Fig. 15, have lower temperatureprofiles. This is pre-
sumablyduetoafastercoolingrate.Accordingtoourmodel,a
lowersurfacedensityofthediscimpliesasmalleropticalthick- 5. Whyisdisceccentricitysoimportant?
nessthatleadstoashortercoolingtimescale.Asaconsequence,
Planetesimals!
ourstandardmodelwithΣ =1Σ hasastationarytempera-
0 MMSN
tureprofilewhichisequivalenttothatofalocallyisothermalrun The shape and profile of the gas disc may have a crucial role
withh 7 8%whilelessmassivediscshavealowertemper- intheearlystagesofplanetaryformation,whenkilometer-sized
∼ −
atureprofileinasteady-statewithatemperaturetypically3 4 planetesimalsarecollidingwitheachother.Evenwhenneglect-
−
times less high. This translates into equivalent locally isother- inggasdiscgravity,severalstudieshaveshownthatthecoupling
mal models with h down to 3 4%. By inspecting Fig. 1 (top betweensecularperturbationsandgasdrag,whichincreasesim-
−
leftplot)wenoticethatinthisrangeofaspectratios,thedisc’s pactvelocitiesbetween non-equalsized bodiesand can lead to
eccentricityactuallychangesquicklywithh.Thismightexplain accretion-hostileenvironments(e.g.The´baultetal.2006, 2008,
10