Resolution requirements for Smoothed Particle Hydrodynamics simulations of self-gravitating accretion discs

0.16 MB·English

by  Giuseppe Lodato

Mon.Not.R.Astron.Soc.000,1–7(2010) Printed14January2011 (MNLATEXstylefilev2.2)
Resolution requirements for Smoothed Particle Hydrodynamics
simulations of self-gravitating accretion discs
1 1⋆ 2
Giuseppe Lodato and C. J. Clarke
1
0 1DipartimentodiFisica,Universita` DegliStudidiMilano,ViaCeloria,16,Milano,20133,Italy
2 2InstituteofAstronomy,MadingleyRoad,Cambridge,CB30HA

n
a
J
Submitted:Revised:Accepted:
2
1

ABSTRACT
]
R Stimulated by recent results by Meru&Bate (2010a,b), we revisit the issue of resolution
S requirementsforsimulatingself-gravitatingaccretiondiscswithSmoothedParticleHydrody-
. namics (SPH). We show that the results by Meru&Bate (2010a) are consistent with those
h
ofMeru&Bate(2010b)iftheyarebothinterpretedasdrivenbyresolutioneffects,therefore
p
implyingthatthe resolutioncriterionfor coolinggaseousdiscsis a functionofthe imposed
-
o coolingrate.Wediscusstwopossiblenumericaloriginsofsuchdependence,whichareboth
r consistentwiththelimitednumberofavailabledata.Ourresultstentativelyindicatethatcon-
t
s vergenceforcurrentsimulationsisbeingreachedforanumberofSPHparticlesapproaching
a 10 millions (for a disc mass of order 10 per cent of the central object mass), which would
[ set the critical cooling time for fragmentation at about 15Ω−1, roughly a factor two larger

1 thanpreviouslythought.Moreingeneral,wediscusstheextenttowhichthelargenumberof
v recentnumericalresultsarereliableornot.Wearguethatthoseresultsthatpertaintothedy-
8 namicsassociatedwithgravitationalinstabilities(suchasthelocalityofangularmomentum
4 transport, and the relationship between density perturbation and induced stress) are robust,
4 whilethosepertainingtothethermodynamicsofthesystem(suchasthedeterminationofthe
2
criticalcoolingtimeforfragmentation)canbeaffectedbypoorresolution.
.
1
Keywords: accretion,accretiondiscs–gravitation–instabilities–stars:formation–galax-
0
1 ies:active
1
:
v
i
X 1 INTRODUCTION fragmentintoboundobjects.Thisissueisextremelyimportantin
r thecontextofplanetformationtheories,asthismechanismforms
a Gravitationalinstabilitiesinaccretiondiscshavebeen intensively
thebasisoftheso-calledgravitationalinstabilitymodelforplanet
studied in the last decade. Several analytical models of grav-
formation.Moreingeneral,suchmechanismmightberesponsible
itationally unstable discs have been proposed over the years
for the formation of stellar clusters in galactic centers, including
(Bertin&Lodato 1999; Balbus&Papaloizou 1999; Clarke 2009;
ourown(Nayakshin2006;Nayakshinetal.2007).Themechanism
Rafikov2009;Riceetal.2010;Krumholz&Burkert2010).How-
leading to fragmentation is fairly well understood. In a cooling
ever, a large number of important properties of self-gravitating
disc, thermal saturation of the instability occurs when the ampli-
discshavebeen determinedthrough theuseofnumerical simula-
tudeofthedensityperturbationislargeenoughthattheassociated
tions(Gammie2001;Lodato&Rice2004,2005;Riceetal.2005;
shockheatingprovidedbytheinstabilitybalancesthecoolingrate
Mejiaetal. 2005; Boleyetal. 2006; Caietal. 2006; Boss 2004,
(Cossinsetal.2009).Iftherequireddensityperturbationexceedsa
2006;Cossinsetal.2009,2010,justtomentionafew).Suchnu-
criticalvalue,thediscfragments.Determiningthecriticalvaluefor
mericalstudieshaveinvestigatedseveral differentaspectsofself-
thedensityperturbationisthusequivalenttodeterminingacritical
gravitatingdiscdynamics.Some(e.g.Lodato&Rice2004,2005;
coolingtimebelowwhichfragmentationoccurs.Suchcriticalcool-
Boleyetal. 2006; Cossinsetal. 2009) have considered the angu-
ing time has been variously determined over the years: Gammie
larmomentumtransportinnon-fragmentingdiscs,withtheaimof (2001) sets it to 3Ω−1 (where Ω is the local angular velocity of
determining whether the long termdisc evolution induced by the
thedisc),basedonhistwodimensionalshearingsheetsimulations,
instability can be treated within a local framework. This is par- Riceetal.(2005)setitto6Ω−1foraspecificheatratioofγ =5/3,
ticularlyrelevant for constructing local analytical models of self-
basedontheirthreedimensionalglobalSPHsimulations,andnote
gravitatingdiscs(Clarke2009;Rafikov2009).Insomeothercases
adependencyofthecriticalcoolingtimeonγ.Otherstudieshave
(Gammie 2001; Riceetal. 2005; Boss 2004; Mayeretal. 2004,
found other dependencies eitheron thespecificformof thecool-
2005; Boss 2006; Mayeretal. 2007; Cossinsetal. 2010), the fo-
ingrate(Johnson&Gammie2003;Cossinsetal.2010)oronthe
cus was on the conditions required for a self-gravitating disc to
2 Lodato& Clarke
thermalhistoryofthedisc(Clarkeetal.2007).Itisconvenientto uatingthisquantityattheonsetoftheinstability,thatisinacon-
expressthecoolingtimeinunitsofthelocaldynamicaltimeinthe ditionwherethediscismarginallystable.Letusthenestimatethe
discandintroducea“coolingparameter”β=Ωt . ratioh/Hforamarginallystabledisc.Wefirstrecallthanmarginal
cool
Recently,somenewresultshavecastdoubtsonthereliability stabilityforaKepleriandiscimpliesthat
ofthedeterminationofsuchcriticalcoolingtime,andinparticular
onthose obtained by SPHsimulations. Inparticular,weconsider Q= csΩ ≈1, (2)
πGΣ
heretwonewresults:first,Meru&Bate(2010a)haverunanum-
berofsimulationswithdifferentslopesofthediscsurfacedensity wherecsisthesoundspeed.Thisisequivalentto
profile,differentdiscmasses,andhaveconsideredthelocationat H ΣR2
whichthefragmentappearedintheirfragmentingruns.Theircon- =πQ =πQm(R), (3)
R M
⋆
clusionisthatthecriticalvalueofβ atagivenradiusRisafunc-
tionofthe“localdiscmass”,m(R) = Σ(R)R2/M⋆,whereΣis wherem(R)=ΣR2/M⋆.
the surface density and M is the mass of the central object. In TheSPHsmoothinglengthissetbythecondition
⋆
ptaatritoincuilnacrr,etahseeysfirnodugthhalyttahsemcr1it/i2c.alSceocoolnindgly,paMraemruet&erBfoartefr(a2g0m10ebn-) ρh3=η3mp, (4)
havenotedthat,foragivendiscsetup(intermsoftotaldiscmass, whereρisthelocalSPHdensity,ηisanumericalparametergener-
surfacedensity profile,etc.)thecriticalβ for fragmentation isan allysetto1.2andmp =Mdisc/N istheSPHparticlemass,where
increasingfunctionofthetotalnumberofSPHparticlesused,and Mdiscisthetotaldiscmass,andN isthetotalnumberofSPHpar-
—moreimportantly—theyfindnohintofconvergence. ticles.Withthisprescription,asmoothingspherecontainsroughly
In this paper, we reconsider the issue of resolution require- 60SPH“neighbours”.Weevaluatethesmoothinglengthatthedisc
mentsinordertosimulatefragmentationinacoolinggaseousdiscs. midplane.Themidplanedensityρcisrelatedtothesurfacedensity
Weintroduceanewresolutionrequirementthatrelatesthenumber byρc = Σ/2H,which,insertedinEq.(4),andusingEq.(3)and
of SPHparticlesused to theexternally imposed cooling rate. We thedefinitionofm(R),gives:
thenshowthatthewholesetofresultsbyMeru&Bate(2010a,b)
h η 2q 1/3
can be naturally explained as a consequence of failing to satisfy =
H m(R)(cid:18)π2Q2N(cid:19)
this condition. Since the Meru&Bate (2010a,b) simulations are
amongthehighestresolutionSPHsimulationsoffragmentingdiscs 5×10−3 q 1/3 N −1/3
≈ , (5)
todate,thisfailureissharedalsobyallpreviousSPHsimulations. m(R) (cid:16)0.1(cid:17) (cid:18)2.5×105(cid:19)
Wewillthendiscusswhichofthepreviousresultsarereliableand
whereq=M /M .Actually,ifthediscisverypoorlyresolved,
whichshouldbeconsideredwithmorecare. disc ⋆
thatisifh/H > 1,theSPHmethodunderestimatesthemidplane
densitybyafactorH/h,andwehaveρc ≈ ρ(h/H).Thiscorre-
spondstoanequivalent2Dversionofequation(4):
2 RESOLUTIONREQUIREMENTSFOR Σh2 ≈2η3mp. (6)
SELF-GRAVITATINGCOOLINGDISCS
Inthispoorlyresolvedcase,onereadilygetsthath/H isgivenby
Several papers in the past have investigated the resolution re-
quirements for resolving gas fragmentation in numerical simula- h η 3/2 2q 1/2
≈
tions. In particular, Bate&Burkert (1997) require that the mass H (cid:18)m(R)(cid:19) (cid:18)π2Q2N(cid:19)
ofSPHparticleswithinasmoothingspherebelessthattheJeans
3.4×10−4 q 1/2 N −1/2
mass.Asimilarcondition,applicablealsotogrid-basedsimulations ≈ . (7)
(Trueloveetal.1997),isthattheminimumresolvablelength(i.e. (m(R))3/2 (cid:16)0.1(cid:17) (cid:18)2.5×105(cid:19)
the mesh sizeor the smoothing length h for grid-based and SPH It is also useful to evaluate h/H at the outer edge of the
simulations, respectively) be less than the Jeans length, although disc, Rout. If the surface density profile is a power-law with re-
notethatfailingtosatisfytheTrueloveconditioningrid-basedsim- spect to R with index p < 2, we have that q = M /M =
disc ⋆
ulations generally induces spurious fragmentation, while in SPH 2πm(Rout)/(2−p),fromwhichweget(forh/H <1)
poor resolution may well lead to suppression of fragmentation.
Nelson(2006)liststhreedifferentconditions:thefirstisthatadap- h 2πη 2 1/3
=
tivegravitationalsofteningisused(thisisactuallyusuallydonein H(cid:12) 2−p(cid:18)q2π2Q2N(cid:19)
(cid:12)out
mostmodernSPHsimulations,includingthosediscussedhere),the (cid:12)(cid:12) q −2/3 N −1/3
secondisavariantoftheBate&Burkert(1997)andTrueloveetal. ≈ 0.3 , (8)
(1997)conditions,andthethirdisthatthesmoothinglengthhbe (cid:16)0.1(cid:17) (cid:18)2.5×105(cid:19)
lessthanthediscthicknessH.However,alsothislastconditionis wherethelastequalityholdsforp=1.Theequivalentof(8)when
equivalenttotheBate&Burkert(1997)andTrueloveetal.(1997) h/H >1is
conditionforamarginallystabledisc,forwhichtheJeanslengthis
actuallyequaltothediscthickness.Alltheseconditionscanthere- h = 2πη 3/2 2 1/2
forebecastinthesimplerequirement: H(cid:12)(cid:12)out (cid:18)2−p(cid:19) (cid:18)q2π2Q2N(cid:19)
Hh .1. (1) (cid:12)(cid:12) ≈ 0.16(cid:16)0q.1(cid:17)−1(cid:18)2.5×N105(cid:19)−1/2. (9)
The natural measure of resolution for an SPH simulation of Note however that equation (9) is of little practical use, since in
anaccretiondiscisthusthequantityh/H.Inordertodiscussthe ordertohaveh/H >1atthediscouteredge,thenumberofSPH
outcomeofgravitationallyunstablediscsweareinterestedineval- particlesneedstobeverysmall,N .6000.
ResolvingfragmentationinSPH simulations 3
WecannowreconsidertheresultsofMeru&Bate(2010a,b)
in the light of Equations (5) and (8) (or (7) and (9)). In fact, on 10
the one hand Meru&Bate (2010a) find that, for a given number
ofparticles,thereisarelationbetweentheimposedvalueofβand
thequantitym(R),evaluated atthelocationofthefirstfragment
8
appearing in their simulations. On the other hand, Meru&Bate
(2010b)demonstratethat,foragivendiscsetup,thecriticalvalue
of β for fragmentation is an increasing function of N. Both re-
sultscan be thus explained if (a) fragmentation issuppressed for 6
lowresolutionand(b)theresolutionrequirement,asmeasuredby
theparameterh/H,isafunctionoftheimposedcoolingtime.In
this picture, Meru&Bate (2010a) thus determine the location at
which,foragivenN,theresolutionrequirementissatisfied,while 4
Meru&Bate (2010b) determine the minimum N at which, for a
givensetup,therequirementissatisfied.
In fact, it would be puzzling to interpret the results of
2
Meru&Bate(2010a)otherwise.Byconstruction,suchsimulations
are self-similar and, as long as global effects do not play a role
(whichisthecasewhenm(R).0.1),thereshouldbenopreferred
locationinthedisc.Theonlypossiblemechanismthatbreaksthe 0
0 1 2 3 4 5
self-similarityisinfactresolution.
Statedotherwise,wecansaythatinordertodetectfragmen-
tationinacoolingdisc,theresolutionrequirementmustdependon
Figure1.Correlationbetweenlinearresolutionatthelocationoffragmen-
β,insuchawaythatdiscswithlongercoolingtimesrequireabet-
tation,asmeasuredbyh/H,andtheimposedcoolingtime,asmeasuredby
terresolution(thatis,asmallerh/H).Inthiscase,fragmentation
theparameterβ.Datapointsaretakenfromtable2ofMeru&Bate(2010a).
wouldthenfirstoccurattheminimumradiusforwhichsuchresolu- Aleftwardpointingerrorbarindicatescaseswherefragmentationoccurred
tionrequirementissatisfied,sincethedynamicalandcoolingtimes nearthediscinnerboundary.Thesolidlineshowsasimpleleastsquares
areshorterthanatlargerradius.Wecanexpresssuchaconditionin fittothedata,assumingalinearrelation(Eq.(11)).Thebestfitmodelin
theform thiscaseisβ=βres=2(H/h).Thedashedlineshowsaleastsquaresfit
tothedataassumingEq.(12).Thebestfitparametersareβ0 = 14.7and
β<βres ≡βres(h/H) (10) a=1.77.
Whatistheactualfunctionalrelationbetweenβresandh/H?
Weexplore thisby plotting in Fig. 1the values of β inthe frag-
expect that at high H/h the data would start deviating from the
menting simulations of Meru&Bate (2010a) versus H/h at the
pointoffragmentation, andassumingthatthecorrelationinthese relationβ=βresandwillsaturateatafinitevalueofβ.
Alternatively, the dashed line shows a fit to the data of the
twoquantitiesispurelydrivenbyresolution,thatisthatβ =βres. form
ThequantityH/hiscomputed usingequation (5)or(7),andus-
ianlsgothineddicaatatepwroitvhidaeldefbtywaMrderpuo&intBinagteer(r2o0r1b0aar)tihnetchaesirestawblheer2e.oWuer βres = (1+aβh0/H)2 (12)
estimateofm(R)(whichisbasedontheinitialdiscsurfaceden- (seeSection3foradiscussionofthenumericaleffectsthatcould
sityprofile)isanover-estimate:inthesecasesfragmentationoccurs giverisetosucharelation).Theimportantdifferencebetweenthis
closetothediscinneredgeandthesurfacedensityhasbeensignif- functional form and that given by equation (11) isthat it implies
icantlydepletedbyboundaryeffectsbythestageoffragmentation. thatβres convergestoafinitevalue(β0)athighH/h.Thebestfit
Thereareevidentlyavarietyoffunctionalformsthatwillpass values of theparameters areβ0 = 14.7 and a = 1.77. Wehave
through the data- wearehampered both by thefact that there is excludedfromthefitthedatapointsthatcorrespondtocompletely
considerable scatter and by the fact that we have no simulation unresolvedsituations,thatish/H >1.
pointswithH/h> 3.Onepossibilityisthatthereisalinearrela- Evidently, either of the lines shown in Figure 1 provide an
tionshipbetweenβandH/h(seeSection3foradiscussionofthis adequate fit to the data over the limited dynamic range of H/h
formintermsoftheroleofartificialviscosity).Wehavefittedthe available.1
datawithaleastsquaresfit(shownwiththesolidlineinFig.1)and We now turn our attention to the results of Meru&Bate
haveobtained: (2010b),whohadfoundanincreaseinthecriticalβ forfragmen-
h −1 tionatanyradiuswithN.Theseauthorsconsiderthecasep = 1
βres ≈2 . (11) andq = 0.1. Inthecasewherewearesimplyinterestedinfrag-
(cid:18)H(cid:19)
Itisworthstressingthatifthetrendβ =βres,givenbyequa- 1 Wealsonotethatthedataisadequatelyfittedbythesuggestedparameter-
tion (11), extended to indefinitely high H/h then it would mean
that therewould be no physical limittothecooling raterequired isationofMeru&Bate(2010a),inwhichβres(inscalingwithm(R)0.5)
effectivelyscaleswith(H/h)0.5.Itisalsoworthnotingtheslightlydiffer-
forfragmentation.Insteaditwouldimplythat-howeverslowthe
entdependenceofourcoordinateH/honstellarmasscomparedwiththat
cooling ratein thedisc- it would always fragment if modeled at ofMeru&Bate:wenotethatalldiscquantitiesdependonthecentralobject
high enough resolution. In contrast, if there is indeed a physical massonlyviatheKeplerianangularvelocity,Ωandthatconsequentlyβres
limittothecoolingraterequiredforfragmentationthenonewould mustdependonm(R)andqinthecombinationm(R)/q1/3.
4 Lodato& Clarke
mentationandnotintheexactlocationwithinthediscwherethis
happens,weonlyneedtomakesurethattheresolutioncriterionis 100
satisfiedatthediscouteredge,thatis
q 2/3 N 1/3
β<βres(h/H)out ≈6.6(cid:16)0.1(cid:17) (cid:18)2.5×105(cid:19) , (13)
wherethelastapproximateequalityholdsinthecasethatweadopt
thelinearfit(equation11),andwherewehaveusedequations(8)
and(10).
ThedatapointsofMeru&Bate(2010b)areplottedinFig.2,
alongwithourpredictionforβresbasedonEq.(13)(solidline)and 10
alsothecorrespondingrelationshipbasedonequation(12)(dashed
line).Weobtainthethresholdforfragmentationβ fromthedata
frag
of Meru&Bate (2010b) by averaging the highest value of their
fragmentingrunswiththelowestnonfragmentingrun.Forthispur-
pose,weconsiderasnonfragmentingthoserunsthatwereclassi-
fiedas‘borderline’byMeru&Bate(2010b).
WeimmediatelyseethattheresultsofMeru&Bate(2010a)
andMeru&Bate(2010b)aremutuallyconsistentifbothsetsofre-
sultsareexplainedintermsofresolutioneffects(thatis,ifβ = 1
frag
βres), regardless of the exact parameterisation that is used to de-
scribe how the cooling rate for fragmentation depends on H/h.
Wealsonotethat,forthehighest number of particlesadopted by
Figure2.Correlationbetweenthecritical coolingtimeforfragmentation
Meru&Bate(2010b),theresultsappeartodeviatefromtheN1/3 and the number of SPH particles used, for the simulations described in
(solid)lineimpliedbythelinear fit(equation 11).Inthesecases, Meru&Bate (2010b).Thethreshold isobtained averaging thehighest β
we thus have non-fragmenting simulations that are well resolved for fragmenting simulations and the lowest β for non-fragmenting ones.
accordingtothiscriterion. Thesolidanddashedlinesshowequations(11)and(12),respectively,eval-
ThesimulationresultsinFigure2arealsobroadlyconsistent uatedatthediscouteredge,andwherethefreeparametersaresetbythefit
with the dashed line (constructed using equation (12)) which is showninFig.1.
aparameterisation that impliesa converged value of therequired
coolingrateforsimulationsthataresufficientlywellresolved.
Figure2,andespeciallythedeviationfromtheN1/3 depen-
denceathighN,thereforecontainstentativeevidencethatconver- surethatsuchexternalheatingsources,andinparticularnumerical
gencemaybebeingreachedatthelargestnumberofparticles,and heatingsources,arenegligible.
thatthe’true’valueofthethresholdcoolingtimeforfragmentation Itiswellknown(Pringle1981)thattheconditionofthermal
wouldthereforeliebetweenβ ≈ 10−15,thatisroughlyafactor equilibriumforaviscousdisccanbeeasilyexpressedintermsof
oftwolargerthanpreviouslythought. arelationbetweentheviscosityparameterα(Shakura&Sunyaev
1973)andthecoolingparameterβ:
4 1
α= . (14)
3 PHYSICALORIGINOFTHERESOLUTION 9γ(γ−1)β
CONDITION
For SPHsimulations, the most significant source of numeri-
Intheprevioussection,wehaveexplainedtheavailablesimulation calheatingisprovidedbyartificialviscosity.Inparticular,itcanbe
dataintermsoftwoalternativeparameterisationsofhowthecool- shown(Artymowicz&Lubow1996;Murray1996;Lodato&Price
ingraterequiredforfragmentationmaydependontheratioofthe 2010)thatthelinearterminthestandardartificialviscosityimple-
smoothinglengthtothediscscaleheight.Wenowturnouratten- mentationofSPHleadstoanequivalentShakura-Sunyaevparam-
tiontopossiblephysicalexplanationsofsuchrequirements. eterαartgivenby
The first possibility is related to artificial viscosity. Indeed,
1 h
theartificialviscosityterminSPHislinearlyproportionaltoh/H αart = αSPH , (15)
10 H
(Monaghan1992),anditseffectmightinprinciplealterthethermal
balanceofthedisc. where αSPH is an SPH input parameter (usually set to 1, but re-
As discussed above, for cooling discs the saturation of the ducedto0.1inthesimulationsdiscussedhere).Actually,theabove
gravitational instability is due to the balance of internal heating relationholdsonlyinthecontinuumlimitandifartificialviscosity
due to shocks induced by the instability itself and the externally isappliedbothinregionsof converginganddivergingflow.Usu-
imposedcooling. Clearly,ifanyother formofexternal heatingis ally,aswitchenablingartificialviscosityonlyforconvergentflows
present,thesaturationamplitudeoftheinstabilitywillbedecreased isused,andsometimesadditionalswitches(Balsara1995)areused,
andfragmentationmightthereforebeinhibitedevenforcoolingpa- sothattheaboverelationisactuallyanoverestimateoftheresulting
rametersthatwouldleadtofragmentationintheabsenceofexter- viscosity.
nalheating.Foractualastrophysicallyrelevantdiscs,suchexternal Whenαart becomescomparabletothevalueimpliedbyEq.
heatingsourcesmightwellbeanimportantphysicalingredientto (14),artificialheatingoffsetscompletelytheimposedcooling,and
betakenintoaccount.Ontheotherhand,innumericalsimulations, no gravitational perturbations will be able to grow. This occurs
whereonewishestodeterminethecriticalcoolingrate,onehasbe when
ResolvingfragmentationinSPH simulations 5
40 1 h −1 h −1 4 DISCUSSIONANDCONCLUSIONS
β≈ =40 , (16)
9γ(γ−1)αSPH (cid:18)H(cid:19) (cid:18)H(cid:19)
In this paper we have shown that the results of Meru&Bate
(2010a) (concerning the location at which fragmentation occurs
where the last equality refers to the case of γ = 5/3 and in SPH simulations of cooling self-gravitating discs) are quanti-
αSPH = 0.1. The resolution criterion proposed in the previous tativelyconsistent(ifinterpretedasbeingdrivenbyresolutionef-
section(equations10and11)thuscorresponds tothecasewhere fects)withtheresultsofMeru&Bate(2010b)inwhichitisshown
artificialviscosity provides 5 per cent of theheating required for thatthecriticalcoolingrateforfragmentationdependsonthenum-
thermalequilibrium. berofparticlesN.Indeedtheresultsofeachpaperimplytheresults
oftheotherandbothimplythatthemeasuredcoolingrateforfrag-
Wedrawattentiontotwoaspectsofthisresult.Firstly,itwould
imply-rathersurprisingly-thatartificialviscosityeffectscanpre- mentation is a function of h/H (where h is the SPH smoothing
ventfragmentationevenwhentheycontributeatinyfraction(5%) lengthandHisthediscscaleheight).
ofthetotalthermalenergybalance.Secondly,itshouldbestressed It isdifficult toassignaprecise functional formfor the way
thatthisconditionisanecessary(butnotnecessarilyasufficient) thatthecriticalcoolingratedependsonresolution,giventhescat-
criterion for fragmentation. As noted in Section 2, there istenta- ter in the simulation data and the relatively small range in H/h
tive evidence from Figure 2 that the points at the highest N fall for which simulations are currently available. We have explored
belowthesolidline-inotherwordswehavesimulationsthatare twopossibilitiesforeffectsthatmayexplainthesimulationresults.
notfragmentingeventhoughtheyareintheregimewherethearti- In one case we explore the possibility that a necessary criterion
ficialviscositycontributionislessthan5%.Weinterpretthisresult forfragmentationtobedetectedisthattheartificialviscositycon-
asevidence thatthereisindeedaphysicalmechanism preventing tributeslessthanacertainfractionofthethermalenergyinput(and
fragmentationforsmallcoolingrates. associatedangular momentumtransport) inthedisc.Wefindthat
the results at low resolution can be explained in these terms but
In addition to the thermal effects described above, artificial
that the requirement on the contribution from artifical viscosity
viscosity might also dynamically stabilize the disc. The specific
issuprisingly stringent(i.e.fragmentation occursonlywherethis
effect of viscosity inthisregard depends onhow viscosity scales
with surface density and might also lead to a secular instability contributeslessthan5%oftheenergyinputtothedisc).Alterna-
tively,weconsiderthepossibilitythatfiniteresolutionjustsmooths
(Schmit&Tscharnuter 1995). While this is beyond the scope of
outdensitypeakssothatfluctuationsthatwouldcollapseifwellre-
thepresentpaper,sucheffectsshouldbefurtherinvestigated.
solveddonotachieveacriticalamplitudeatfiniteh.
Sincewefinditrathersurprisingthatartificalviscosityshould
Irrespective of how one explains how resolution affects the
suppress fragmentation when it contributes such a minor compo-
cooling requirements for fragmentation, an important issue - as
nenttothedisc’sthermalbalance,wealsoexplorethehypothesis
stressed by Meru&Bate (2010b) - is whether there is indeed a
that theeffect of under resolution isoneof artificiallysmoothing
convergenceintherequiredcoolingrateathighresolution.Inother
densityenhancementsovertheSPHsmoothinglengthh.Inpoorly
words,istherealevelofresolutionabovewhichonerecoversthe
resolvedsimulations,wherehisnotmuchsmallerthanthelength
resultthathasbeenassumedhitherto-i.e.thatthereisaphysical
scaleofdensitypeaks(λ),thenthiswillsuppressthepeakampli-
conditiononthecoolingraterequiredforfragmentation?Weevi-
tudeoftheresultingdensityfluctuationsbyafactorof1+h/λ(as-
dentlycannotanswerthisquestionwithcalculationsthathavenot
sumingthat gravitational instabilitiesgeneratepredominantly lin-
earstructures).AccordingtoCossinsetal.(2009),thegravitational attained convergence (see Figure 2 where βfrag rises mildlywith
heating rate associated with modes of r.m.s. fractional amplitude N evenatthehighestN valuesstudied).However,Figure2con-
∆Σ/Σisproportionalto(∆Σ/Σ)2andthereforeinthermalequi- tainshintsofapproaching convergence. Onewaytoseethisisto
libriumthisquantity scaleswiththe cooling rate(∝ β−1).If we comparethesimulationdatawiththesolidlinewhichcorresponds
assumethatinaperfectlyresolvedcalculationthepeakamplitude torequirement of a5% contributionfromaritificalviscosity. The
wouldscalewithther.m.s.amplitude(butbedegradedbyafactor factthattherighthandpointsliebelowthissolidline(whereβfrag
1+h/λ in the case of finite h) and if we furthermore associate scaleswithN1/3)isevidencethatthisconditionmaybenecessary
butnotsufficientforfragmentation.Similarly,thesimulationdata
fragmentationwiththepeakfluctuationsachievingacriticalvalue
of∆Σ/Σ(oforderunity)thenitfollowsthatthecoolingratere- isconsistentwiththedashedlineinwhichtheactualβrequiredfor
quiredforfragmentationisincreasedinthecaseofsimulationsof fragmentationconvergestoavalueofβ0 ≈14.7athighresolution.
finiteh.Specificallywehave Clearlythequestionofultimateconvergencewillonlybeset-
tledbyfurthersimulationsathigherN.Ifconvergenceisnotulti-
matelyattainedthenitwould involvearadical re-thinkof allour
βres =β0(1+h/λ)−2 (17) understanding ofgravitationallyunstablediscssinceitwouldim-
plythatanyself-gravitatingdisc(howeverslowlycooling)should
inrealityfragment.Thiswouldbesurprisingbecause, inthecase
where β0 isthevalue required for fragmentation in thecase of a ofveryslowcooling,ther.m.s.modeamplituderequiredtoachieve
wellresolvedsimulation(i.e.smallh).Inestimatingλwenotethat thermal balance should be very low and it would be unexpected
theJeanslengthinadiscwithToomreQparameterclosetounity (thoughpossibleinprinciple)thatsuchadiscneverthelessexhib-
isjusttheverticalscaleheightH;itthereforeseemslogicalthatλ itedlocallynon-lineardensityfluctuationsgivingrisetocollapsing
shouldscalewithH andthismotivatestheformofequation(12). fragments.
Wenotethat both our resolutioncriteriadonot abandon the The more conservative conclusion (if convergence is ulti-
notionthat(inawellresolvedsimulation)thereisacriticalcooling matelyattained)isthatwehavesimplyunder-estimatedthecritical
rate associated with fragmentation but take into account that the βvaluehitherto.Thiswouldinvolvesomequantitativeadjustment
factthatmorevigorouscooling(lowerβ)isrequiredinthecaseof of previous conclusions. For example, if the critical cooling time
poorlyresolvedsimulations. forfragmentationisroughlyafactor2higherthanhadbeenprevi-
6 Lodato& Clarke
ouslythought thenthisaffectsthelocationatwhichgiant planets codes.Artificialviscosityisalsopresenttovariousdegreesisgrid-
canformthroughgravitationalinstability(Meru&Bate2010b):in basedcodes,andonemightthusconcludethatananalogousreso-
theregimeoficecoolingappropriatetotheouterpartsofprotostel- lutionrequirementshouldbeformulatedalsointhuscase.Whether
lardiscs,thischangewouldbringintheminimumradiusforfrag- alsogrid-basedsimulationsdonotconvergeinthedeterminationof
mentation by afactor ∼ 22/9 (Clarke 2009; Cossinsetal.2010). thecriticaltimescaleforfragmentationandtheexactformofthe
Although thischange isquite modest, it may beimportant tothe resolutionrequirement inthiscasearequestions thatstillneedto
discussion of whether recently imaged planets (e.g., around HR beworkedout.
8799, Maroisetal. 2008, around β Pic, Lagrangeetal. 2009 or
around Fomahault, Kalasetal. 2008) could have been formed by
gravitationalinstability.Note,however,that—foragivencritical
ACKNOWLEDGEMENTS
cooling time — the uncertainty in the determination of the frag-
mentation radius due to uncertainties in the relevant opacities, in We thank Phil Armitage, Matthew Bate, Farzana Meru and Ken
thedetailedverticalstructureofthediscandduetotheeffectsof Ricefor interesting discussions. We also thank Charles Gammie,
magnetic fields are probably larger than the modest change dis- foraninsightfulreferee’sreport.
cussedabove.
More in general, we might ask which ones of the properties
ofnonfragmentingself-gravitatingaccretiondiscsaregoingtobe
REFERENCES
significantlyaffectedbypoorresolution,inthesenseexpressedby
Eq. (10). Probably the most important property, apart from frag- ArtymowiczP.,LubowS.H.,1996,ApJ,467,77L
mentation,isrelatedtothestressinducedbytheinstability,andin BalbusS.A.,PapaloizouJ.C.B.,1999,ApJ,521,650
particularitsmagnitudeandthelocalityoftheassociateddissipa- BalsaraD.S.,1995,JournalofComputationalPhysics,121,357
tion(whichareinturnrelatedtothepowerspectrumofthemodes BateM.R.,BurkertA.,1997,MNRAS,288,1060
excitedbytheinstability). BertinG.,LodatoG.,1999,A&A,350,694
As mentioned in the Introduction, several studies BoleyA.C.,Mej´ıaA.C.,DurisenR.H.,CaiK.,PickettM.K.,
(Lodato&Rice 2004, 2005; Cossinsetal. 2009) have con- D’AlessioP.,2006,ApJ,651,517
sidered the issue of locality of the transport induced by the BossA.P.,2004,ApJ,610,456
instability.Herethequestioniswhetherthetypicalwavelengthof BossA.P.,2006,ApJ,641,1148
thedisturbancesisasignificantfractionoftheradiusR.Theabove CaiK.,DurisenR.H.,MichaelS.,BoleyA.C.,Mej´ıaA.C.,Pick-
studies(andinparticularCossinsetal.2009)haveshownthatthe ettM.K.,D’AlessioP.,2006,ApJ,636,L149
spectrum of the disturbances peaks at a wavelength λ ≈ H and ClarkeC.,Harper-ClarkE.,LodatoG.,2007,MNRAS,381,1543
becomesnegligibleforλ & 3H.Inthisrespect,thus,therelevant ClarkeC.J.,2009,MNRAS,396,1066
resolution requirement is the generally less restrictive condition CossinsP.,LodatoG.,ClarkeC.,2009,MNRAS,393,1157
h/H . 1 already discussed in Nelson (2006). Furthermore, not CossinsP.,LodatoG.,ClarkeC.,2010,MNRAS,401,2587
resolving the smallest scales would only prevent the disc from GammieC.F.,2001,ApJ,553,174
developing small scale structures, which would not contribute to JohnsonB.M.,GammieC.F.,2003,ApJ,597,131
globaltransportanyway. KalasP.,GrahamJ.R.,ChiangE.,FitzgeraldM.P.,ClampinM.,
Letusnowdiscusswhethertheevaluationofthestressfrom KiteE.S.,StapelfeldtK.,MaroisC.,KristJ.,2008,Science,322,
suchsimulations isaffectedby Eq.(10).Clearly,foragivendis- 1345
turbance,regardlessofthewayithasbeenexcited,acorrecteval- KrumholzM.R.,BurkertA.,2010,ArXive-prints
uation of the resulting stress only requires that the typical wave- LagrangeA.,GratadourD.,ChauvinG.,FuscoT.,EhrenreichD.,
lengthsofthedisturbanceareresolved.However,herethequestion MouilletD.,RoussetG.,RouanD.,AllardF.,GendronE´.,Char-
iswhetherthediscresponsetoanexternalboundarycondition(in tonJ.,MugnierL.,RabouP.,MontriJ.,LacombeF.,2009,A&A,
this case, to an externally imposed cooling rate) depends on res- 493,L21
olution. The condition of thermal equilibrium relates in a unique LodatoG.,PriceD.J.,2010,MNRAS,405,1212
waytheexternalcoolingtothemagnitudeofthestressinducedby LodatoG.,RiceW.K.M.,2004,MNRAS,351,630
theinstability(eq.(14)),underthehypothesisthattheonlyheating LodatoG.,RiceW.K.M.,2005,MNRAS,358,1489
sourceinthediscisprovidedbytheinstabilityitselfandisdissi- MaroisC.,MacintoshB.,BarmanT.,ZuckermanB.,SongI.,Pa-
patedlocally.Here,wehaveshownthatthenewempiricallydeter- tienceJ.,LafreniereD.,DoyonR.,2008,ArXive-prints
minedresolutionconditioncorrespondstothecasewhereartificial MayerL.,LufkinG.,QuinnT.,WadsleyJ.,2007,ApJ,661,L77
viscosityonlyprovides5percentoftheheatingnecessaryforther- MayerL.,QuinnT.,WadsleyJ.,StadelJ.,2004,ApJ,609,1045
malbalance.Therefore,unlesstheeffectsofartificialviscosityare Mayer L.,Wadsley J., Quinn T., Stadel J., 2005, MNRAS,363,
muchlargerthanexpected, onewouldnotexpect achange inthe 641
stressabovethe5percentlevel.Inthisrespect,theeffectofpoor MejiaA.C.,DurisenR.H.,PickettM.K.,CaiK.,2005,ApJ,619,
resolutionisjusttosuppressthedevelopmentofsmallscaledistur- 1098
bances,andthusreducethepeakvalueofthedensityperturbation, MeruF.,BateM.R.,2010a,MNRAS,410,559
whilekeepingtheoverallaveragevalueofthestressessentiallyun- MeruF.,BateM.R.,2010b,MNRAS,inpress
changed. MonaghanJ.J.,1992,ARA&A,30,543
Finally, it is worth noting that the discussion above is fo- MurrayJ.R.,1996,MNRAS,279,402
cussedonSPHsimulations.Manypreviousresults(Gammie2001; NayakshinS.,2006,MNRAS,372,143
Mejiaetal.2005;Boleyetal.2006),includingseveraldetermina- NayakshinS.,CuadraJ.,SpringelV.,2007,MNRAS,379,21
tions of the fragmentation boundary (Gammie 2001; Mejiaetal. NelsonA.F.,2006,MNRAS,373,1039
2005; Boss 2004, 2006), have been obtained using grid-based PringleJ.E.,1981,ARA&A,19,137
ResolvingfragmentationinSPH simulations 7
RafikovR.R.,2009,ApJ,704,281
RiceW.K.M., LodatoG.,ArmitageP.J.,2005, MNRAS,364,
L56
RiceW.K.M.,MayoJ.H.,ArmitageP.J.,2010,MNRAS,402,
1740
Schmit,U.,Tscharnuter,W.M.1995,Icarus,115,304
ShakuraN.I.,SunyaevR.A.,1973,A&A,24,337
TrueloveJ.K.,KleinR.I.,McKeeC.F.,HollimanIIJ.H.,Howell
L.H.,GreenoughJ.A.,1997,ApJ,489,L179