An improved sink particle algorithm for SPH simulations

0.73 MB·English

by  D. A. Hubber

Mon.Not.R.Astron.Soc.000,1–??(2012) Printed22January2013 (MNLATEXstylefilev2.2)
An improved sink particle algorithm for SPH simulations
D. A. Hubber1,2, S. Walch3,4, A. P. Whitworth3
3 1Department of Physics and Astronomy, University of Sheffield, HicksBuilding, Hounsfield Road, Sheffield, S3 7RH, UK
1 2School of Physics and Astronomy, Universityof Leeds, Leeds, LS2 9JT, UK
0 3School of Physics and Astronomy, Cardiff University,Queens Buildings, The Parade, CF24 3AA, UK
2 4Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, Garching, D84758, Germany

n
a
J December4th,2012

8
1
ABSTRACT

] Numerical simulations of star formation frequently rely on the implementation of
M sink particles,(a) to avoidexpending computationalresourceon the detailed internal
physicsofindividualcollapsingprotostars,(b)toderivemassfunctions,binarystatis-
I
. ticsandclusteringkinematics(andhencetomakecomparisonswithobservation),and
h
(c) to model radiative and mechanical feedback; sink particles are also used in other
p
contexts,for example to representaccreting black holes in galactic nuclei. We present
-
o a new algorithm for creating and evolving sink particles in SPH simulations, which
r appearsto representa significantimprovementoverexistingalgorithms– particularly
t
s insituationswheresinksareintroducedafterthegashasbecomeopticallythicktoits
a
own cooling radiation and started to heat up by adiabatic compression. (i) It avoids
[
spurious creation of sinks. (ii) It regulates the accretion of matter onto a sink so as

1 to mitigate non-physical perturbations in the vicinity of the sink. (iii) Sinks accrete
v matter, but the associatedangular momentum is transferredback to the surrounding
0 medium.Withthenewalgorithm–andmodulotheneedtoinvokesufficientresolution
2 to capture the physics preceding sink formation – the properties of sinks formed in
5
simulationsareessentiallyindependentoftheuser-definedparametersofsinkcreation,
4
or the number of SPH particles used.
.
1
0 Key words: stars:formation–galaxies:nuclei–methods:numerical–hydrodynam-
3 ics – gravitation
1
:
v
i
X
1 INTRODUCTION Sinkscanalsobeintroducedintotheinitialconditionsfora
r simulation,torepresentpre-existingcondensations(starsor
a
Starformationisacollectiveprocess:manystarsarebornin black holes). A basic sink is a point mass, and therefore it
multiplesystems(e.g.Duquennoy& Mayor1991),andmost interactswiththesurroundingmattergravitationally. Once
areborninclusters(e.g.Lada & Lada2003).Consequently, created, a sink can grow by accreting nearby matter that
a realistic simulation that captures the important interac- looks likeit would inevitably become part of thesame con-
tionsbetweenneighbouringprotostars,astheyform,isonly densation. Even at this most basic level, there are issues
feasibleifsomeofthedetailedinternalphysicsofindividual with how to identify lumps that should be converted into
protostars is sacrificed. Otherwise all the available compu- sinks,andhowtodeterminewhichmattertheyshouldsub-
tational resource is commandeered by the first protostar to sequentlyaccrete.
form, and events in the rest of the computational domain
At a higher level of sophistication, sink properties
grind to a halt. It is for this reason that Lagrangian sink
(masses,positions,velocities,accretionrates)areusedtoin-
particles were first developed by Bate et al. (1995). Subse-
ferstellarmassfunctions,binarystatistics,clusterkinemat-
quently they have also been used to model single accreting
icsandthestarformationrateinturbulentmolecularclouds
objects,forexampleblackholesingalacticnuclei,wherethe
(e.g. Bonnell et al. 2001; Kitsionas & Whitworth 2002;
focusisnotonthedetailedhydrodynamicsintheimmediate
Bate et al. 2003; Goodwin et al. 2004a,b; Bonnell et al.
vicinity of the black hole, but on feedback from the black
2004;Goodwin et al.2006;Bonnell & Bate2006;Bate2009;
hole on much larger scales (e.g. Springel et al. 2005).
Attwood et al. 2009; Stamatellos et al. 2011; Walch et al.
In the most basic implementation, a sink is created 2012; Federrath & Klessen 2012). They can also be
when and where some lump in the computational domain used to model radiative and mechanical feedback from
lookslikeitwillinevitablycondenseoutgravitationally(into protostars, using simple models or phenomenological
a single star, or a tight multiple system, or a black hole). prescriptions for the physics inside the sink (e.g.
2 Hubber, Walch and Whitworth
Stamatellos et al. 2005; Dale et al. 2005; Krumholz et al. bust, in the sense that there is good correspondence be-
2007; Stamatellos & Whitworth 2009b; Offneret al. 2009; tween the properties of sinks formed in simulations having
Bisbas et al.2009;Bate2010;Peters et al.2010;Wang et al. thesameinitialconditionsbutdifferentℓ .Tocapturethe
MIN
2010; Dale & Bonnell 2011; Peters et al. 2011; Bate 2012; fullrangeoffragmentation,andtheformationoffirstcores,
Cunningham et al. 2011; Dale & Bonnell 2012). requires ℓ ≪ 3AU, for contemporary local star forma-
MIN
In cosmological simulations, sink particles are used tion.WenotethatinsomeAMRsinkimplementations(e.g.
to model the formation and growth of black holes, and Federrath et al.2010)theangularmomentumofmatterac-
henceAGNfeedback (Springelet al. 2005; DiMatteo et al. creted by the sink is not left behind, but is assimilated by
2005, 2008; Johansson et al. 2009; Debuhret al. 2010, thesink.
2011; Choi et al. 2012). For example, in Springel et al. Recently,Federrath et al.(2010)havesimulatedthefor-
(2005) the efficiency of AGN feedback is estimated using mation of a star cluster with SPH and AMR,using compa-
the Bondi-Hoyle-Littleton theory (Hoyle& Lyttleton 1939; rable resolution, and have developed a set of sink-creation
Bondi & Hoyle 1944; Bondi 1952) acting on scales of a few criteria which ensure that the two methods give very simi-
tensofparsecs,whereastheactualSchwarzschildradiusofa lar results, in terms of the masses, locations and accretion
typicalsupermassiveblackholeisoforderafew solarradii. historiesofthesinksformed.However,althoughthesesimu-
Two main algorithms are currently used to model star lationsinvestigatetheconsequencesofdifferentsinkcreation
formation, Smoothed Particle Hydrodynamics (SPH; Lucy criteria,theydonotexploreindetailthedependenceofthe
1977; Gingold & Monaghan 1977) and Adaptive Mesh Re- sinkpropertiesontheresolution, orontheuser-definedpa-
finement (AMR; Berger & Colella 1989; Berger & Oliger rametersofsinkcreation.Thepresentpaperaddressesthese
1984; Dezeeuw& Powell 1993; MacNeice et al. 2000). This issues.
paper is concerned with sinks in SPH. In Section 2 we describe in detail a new algorithm for
In SPH simulations of star formation, two require- the creation and evolution of sinks, explaining the motiva-
ments are critical. First, the mass of a single SPH tion for each of its features; we call these NewSinks. In
particle, m , should be chosen (e.g. Bate & Burkert Sections3 and4,wedescribe thesinkalgorithms that have
SPH
1997; Whitworth 1998; Attwood et al. 2009; Bate 2009; beenusedpreviouslyinSPHsimulationsofstarformation–
Stamatellos & Whitworth 2009a), or adjusted (e.g. hereafterstandard sinks–focussingontheparticularimple-
Kitsionas & Whitworth 2002), so that the Jeans mass mentations that we use for comparison tests. Since all but
is always resolved. Second, sink creation should be consid- one of the new features in NewSinks relate to their evolu-
ered only if the density exceeds a user-defined threshold, tion,ratherthantheircreation,Section3definesOldSinks,
ρ ,andρ shouldinturnexceedthedensityatwhich whicharecreatedinthesamewayasNewSinks,butevolved
SINK SINK
the minimum Jeans mass obtains, so that the code can in the manner of standard sinks, thereby allowing us to
capture properly the full range of fragmentation, and the isolate those aspects of NewSinks that have to do purely
formation of first cores. Thus our main concern here is with sink evolution. On the other hand, since most extant
with sinkscreated from gas that is already being heatedby SPHsimulationsofstarformationthatinvokesinksactually
adiabaticcompression;foradiscussionofhowtotreatsinks use a less robust algorithm to create sinks than that used
created at lower densities from approximately isothermal by NewSinks and OldSinks, Section 4 defines UrSinks,
gas the reader is refered to Federrath et al. (2010). For which are created and evolved in the manner of standard
contemporary local star formation, involving molecular gas sinks.Thebasic propertiesof thedifferenttypesof sink are
at T∼10K,these requirementsreduceto m .10−5M summarisedinTable1.Section5presentstheresultsoftests
SPH ⊙
andρ ≫10−13gcm−3.Oncethedensityexceedsρ , performed with the different types of sink to demonstrate
SINK SINK
the criteria that are applied to determine (i) whether a thatwhenNewSinksareusedtheresultsobtainedareonly
sink actually is created, (ii) what matter it immediately weaklydependentontheresolutionandtheuser-definedpa-
assimilates, and (iii) what matter it subsequently accretes, rametersofsink-creation,whereaswhenOldSinksareused
must be formulated so as to obtain close correspondence the results are strongly dependent and not converged, and
between the properties of the sinks that do form, and the when UrSinks are used the results are divergent. General
stars that should have formed, under the given physical aspectsofthedifferenttypesofsinkarediscussedinSection
circumstances. 6. Ourmain conclusions are summarised in Section 7.
In AMR simulations of star formation, the computa-
tional domain is divided adaptively into finermeshes – and
asappropriatethesemeshesaresubsequentlyde-constructed
– so as to deliver high resolution only when and where it
is needed. There is a minimum mesh spacing, ℓ , below
MIN
which refinement is prohibited. In the original implemen-
tation of sink particles in AMR (Krumholzet al. 2004), a
2 THE NewSink ALGORITHM
sink is introduced wherever the local Jeans length falls be-
low 4ℓ and is therefore no longer adequately resolved. Here we describe and justify the procedures used to define
MIN
There are then rules to determine the rate at which mass theinternalstructureandtranslationofaNewSink(Section
is accreted onto thesink from thesurrounding cells, and to 2.1), to trigger the creation of a NewSink (Section 2.2), to
identify the circumstances under which neighbouring sinks accrete matter onto an existing NewSink (Section 2.3), to
aremerged.Theangularmomentumoftheaccretedmassis update the properties of a NewSink (Section 2.4), and to
leftbehindinthecellfromwhichitisaccreted.Offneret al. pass the angular momentum acquired by a NewSink back
(2008) have shown that this implementation is quite ro- to thesurroundingmatter (Section 2.5).
An improved sink particle algorithm for SPH simulations 3
Type of Sink −→ NewSink OldSink UrSink
Creation criteria: Equation:
density ρi>ρSINK (2) ⋆ ⋆ ⋆
overlap |ri−rs|>XSINKhi+Rs (3) ⋆ ⋆ ⋆
potential minimum φi<min{φj} (4) ⋆ ⋆ —
HillSphere ρi>ρHILL= ... (5) ⋆ ⋆ —
accelerationdivergence (∇·a)i<0 (33) — — ⋆
velocitydivergence (∇·v)i<0 (34) — — ⋆
non-thermalenergy E <0 (35) — — ⋆
NT
Structure: Evolution:
interaction-zone regulatedaccretion (cid:4) — —
exclusion-zone instantaneous accretion — (cid:4) (cid:4)
Table 1.Thecreationcriteria,intrinsicstructures,andevolutionproceduresforthedifferenttypes ofsink.Column1givesthevarious
creationcriteriaandthetwodifferentintrinsicstructures.Column2givesequationsforthevariouscreationcriteria,andthetwodifferent
modes of evolution. Column 3 gives equation numbers. In columns 4 (NewSink), 5 (OldSink) and 6 (UrSink), stars indicate which
creationcriteriaareinvoked,andfilledsquaresindicatewhichstructure/evolutioncombinationisused.RelativetoNewSinks,OldSinks
have the same creation criteria, but different evolution procedures; UrSinks have different creation criteria and different evolution
procedures.MostextantSPHsimulationsofstarformationhavebeenperformedwithstandardsinkslikeUrSinks.
2.1 Internal structure and advection of a NewSink ρ >ρ . (2)
i SINK
A NewSink comprises a central point mass (hereafter, the Weadvocateρ >10−11gcm−3,sothatNewSinksnor-
SINK
point-mass) and a concentric spherical volume of radius Rs mallyform incondensationsthatarewellintotheirKelvin-
(hereafter, the interaction-zone) that moves with, and re- Helmholtz contraction phase.
mains centred on, the point-mass. When we refer to the
mass, M , position, r , velocity, v , or angular momentum,
s s s
L ,ofaNewSink,wemeanthemass,position,velocityand
s
spin angular momentum of its point-mass. Once created, a 2.2.2 Sink-overlap criterion
NewSinkistrackedwiththesameintegrationroutineasan
ThesecondcriterionforSPHparticleitotriggertheforma-
SPHparticle, modulothat it only experiences gravitational tion of a NewSink is that it should not, at the moment of
acceleration. The interaction-zone contains live SPH parti- creation, overlap another NewSink, i.e.
cles that have not yet been assimilated by the point-mass,
andwhichareintendedtorepresentthecontinuationofany |ri−rs′| > XSINKhi + Rs′, (3)
accretiondisc,inflowstream,orothercoherentstructureim-
for all pre-existing NewSinks,s′.
mediatelyoutsidetheinteraction-zone.Theinteraction-zone
isoneofthekeyfeaturesthatdistinguishesaNewSinkfrom
a standard sink.
2.2.3 Gravitational potential minimum criterion
2.2 Criteria for creation of a NewSink The third criterion for SPH particle i to trigger the forma-
tion of a NewSink is that its gravitational potential, φ ,
i
The creation of a NewSink is triggered by an SPH parti- should be lower than that of all it neighbours, j, i.e.
cle. Specifically, if SPH particle i, having smoothing length
hi, satisfies the four criteria below, it is replaced with a φi < min{φj} . (4)
NewSink, s. At the moment of creation, the point-mass of
This criterion, introduced by Federrath et al. (2010), en-
shasthesamemass,position andvelocityasidid,andthe
sures that NewSinks are only created near resolved po-
radius of its interaction-zone is
tential minima (and hence, presumably, resolved density
R = X h . (1) peaks).
s SINK i
X isauserdefinedparameter,chosensothattheneigh-
SINK
bours of i all fall inside the NewSink at the moment of its
creation. Thus, for the M4 kernel, we advocate X = 2. 2.2.4 Hill Sphere criterion
SINK
ByadoptingalargerX oneobtainssmootheraccretion
SINK Thefourthcriterion forSPHparticleitotriggertheforma-
onto a NewSink,but at theexpenseof coarser resolution.
tion of a NewSinkis that its density should satisfy
ρ > ρ ≡ 3XHILL (−∆ris′·∆ais′), (5)
2.2.1 Density criterion i HILL 4πG|∆ris′|2
ThefirstcriterionforSPHparticleitotriggerthecreationof for all pre-existing NewSinks, s′. Here X is a user-
HILL
aNewSinkisthatthedensityatitslocation,ρi≡ρSPH(ri), defined parameter with default value XHILL = 4, ∆ris′ ≡
should exceed a user-definedthreshold, ρSINK, i.e. ri−rs′ and∆ais′ ≡ai−as′ arethepositionandacceleration
4 Hubber, Walch and Whitworth
ofparticleirelativetos′.1 Thiscriteriondealswiththesit- M4 kernel, H =R /2. Consequently SPH particles in the
s s
uation where a condensation undergoing Kelvin-Helmholtz outerregions of theinteraction-zone makea smaller contri-
contraction has a NewSink s′ in its interior, but is much bution than those closer to thepoint-mass. W ensuresthat
moreextendedthanthatNewSink.Eqn.(5)ensuresthata thesum is accurately normalised.
secondNewSink,seededbySPHparticlei,canforminthe
outskirtsof thecondensation only if thereis adensity peak
2.3.3 Timescale for disc accretion
at r that dominates thelocal gravitational field.
i
We model disc accretion using the Shakura& Sunyaev
(1973) prescription, which conflates all possible angular
2.3 Criteria for accretion by a NewSink
momentum transport mechanisms into a single parameter,
2.3.1 The SPH particles inside a NewSink α .Foralow-mass discin approximateKeplerian rotation
SS
aroundastarofmassM ,theaccretion timescale atradius
Under normal circumstances, SPH particles that enter the Risthen ∼α−1(GM R⋆)1/2a−2,whereaisthelocal sound
interaction-zone of a NewSink are not accreted immedi- SS ⋆
speed. We therefore compute a kernel-weighted mean over
atelybyitspoint-mass.Inthefirstinstance,theyaresimply
all the SPH particles in theinteraction-zone,
added to a list of the SPH particles inside the interaction-
zone, its interaction-list. These SPH particles may subse- ht i = (GMs)1/2 |∆rjs|1/2mjW(|∆rjs|,Hs) (.10)
quently be accreted by the point-mass, possibly over sev- DISC αSSW j (cid:26) ρja2j (cid:27)
eral timesteps, but they may leave the interaction-zone be- X
fore this happens. If an SPH particle finds itself inside the Both observational estimates of protostellar accretion disc
interaction-zones of more than one existing NewSink, it is lifetimes (e.g. Hartmann 1998), and theoretical simulations
added to the interaction-list of the NewSink whose point- (e.g. Forgan et al. 2010), suggest that αSS = 0.01 is a suit-
massisclosest.Inthiscontext,itisappropriatetonotethat, able default value, but thereare large uncertainties.
although a NewSink is never created overlapping another
NewSink,itsmotionmaythereafterleadtoanoverlap.We
2.3.4 Net timescale for accretion
do not allow sinksto merge (cf. Krumholzet al. 2004).
To determine whether, and how quickly, the SPH par- The timescale for accretion onto thepoint-mass is given by
ticles in the interaction-zone of an existing NewSink are t = ht i(1−f) ht if , (11)
assimilatedbyitspoint-mass,weconsiderthelimitingcases ACC RAD s DISC s
2E
ofspherically-symmetricradialaccretion,anddiscaccretion. f = min ROT , 1 . (12)
|E |
(cid:26) GRAV (cid:27)
Here E and E are – respectively – the net rota-
2.3.2 Timescale for spherically symmetric radial accretion ROT GRAV
tionalandgravitational energiesoftheSPHparticlesinthe
Insphericalsymmetry,therateatwhichmassflowsinwards interaction-zone, relative tothe point-mass, i.e.
across a spherical surface at radius r is |L |4
E = INT , (13)
M˙(r) = −4πr2ρ(r)vRAD(r). (6) ROT 2 j{mj|∆rjs·LINT|2}
AtthepositionofanSPHparticlej,intheinteraction-zone E = GPMs m φ |∆rjs| +φ |∆rjs|
of NewSink s, this inflow rate becomes GRAV 2 j j(cid:26) (cid:18) Hs (cid:19) (cid:18) hj (cid:19)(cid:27)
X
MT˙hje t=imes−cal4eπf|o∆rrrjasd|∆ialrjasc·c∆revtjisonρji.s obtained by dividi(n7g) +G4Xj jX′6=jmjmj′(cid:26)φ(cid:18)|∆hrjjj′|(cid:19)+φ(cid:18)|∆hrjj′j′|(cid:19)(cid:27).
thenetmassoftheSPHparticlesintheinteraction-zoneby (14)
a weighted sum of theinflow rates:
htRADis = 4π {|∆r |∆r j·∆{mvj}mWW(|∆r |,H )},(8) LINT = Xj {mj∆rjs×∆vjs} (15)
js Pjs js j js s
j is the net angular momentum of the SPH particles in the
W = P{m W(|∆r |,H )/ρ } . (9) interaction-zone,relativetothepoint-mass,andφisafunc-
j js s j
tion that encapsulates the kernel-smoothing of the gravita-
j
X
tional potential (see Eqn.15 in Hubberet al. 2011).
The weighting here uses the kernel-function, W. The
Eqn.11isadopted becauseit givesthecorrect limiting
smoothinglengthofthesink,H ,isadjustedsothattheex-
s behaviour. If the SPH particles inside the sink are in rota-
tentof thekernelfunction equals R ;for example, with the
s tional equilibrium, t →ht i . Conversely, if they are
ACC DISC s
not rotating at all, t →ht i .
ACC RAD s
1 Here,andinthesequel,weadopttheconvention that∆qab≡
qa−qbisthedifferenceinsomequantityqbetweentwoparticlesa 2.3.5 Excising SPH particles
andb,evaluatedatthesametime;thus,forexample,∆rij ≡ri−
rj istheinstantaneouspositionofSPHparticleirelativetoSPH Normally, the mass accreted by the point-mass at the end
particle j. In contrast, we adopt the convention that δqc is the of thecurrent timestep, (t,t+δt ), is
s
changeinsomequantityqc associatedwithparticlec,duetothe
evolution of the system; thus, forexample, δLs isthe increment δM =M 1−exp − δts . (16)
totheangularmomentumofsinks,andδts isthetimestepofs. ACC INT (cid:20) (cid:18) tACC(cid:19)(cid:21)
An improved sink particle algorithm for SPH simulations 5
In the first instance, this mass is removed from the SPH
r′ = M′−1 M r + {δm r } , (18)
particle closest to the point-mass. If the mass of this SPH s s s s j j
!
particle is less than δM , theremainderis removed from Xj
ACC
thesecondclosestSPHparticle,andsoon,untilthewholeof
v′ = M′−1 M v + {δm v } , (19)
δM hasbeenremoved;becausethisonlyaffectsparticles s s s s j j
ACC !
j
close to the point-mass, it is extremely rare that particles X
withreducedmassleavetheinteraction-zone.Therearetwo L′s = Ls + Ms∆rss′×∆vss′
circumstances underwhich the procedure outlined above is + {δmj∆rjs′×∆vjs′} . (20)
superseded.
j
X
(i)IfthetotalmassofSPHparticlesintheinteraction-
Here, the summation is over all SPH particles, j, that lose
zoneofaNewSinkpresentlyexceedsthetotalmassofSPH
mass, δm , to the point-mass, and the updated values are
particles inside the interaction-zone at the time it was cre- j
denoted by primes. Once an SPH particle has zero mass, it
ated, M , then t is decreased artificially by a factor
MAX ACC is removed from thesimulation altogether.
t →t /(M /M )2,inordertoaccretetheexcess
ACC ACC INT MAX
mass more rapidly. This is similar to the procedure used to
transfermassfromnearbygridcellstosinkparticlesinAMR
2.5 Angular momentum feedback from a NewSink
(Krumholzet al. 2004).
(ii) If the timestep, δt , for an SPH particle, j, Ifthepoint-massweretoassimilate alltheangularmomen-
j
in the interaction-zone of NewSink s, satisfies δt < tumoftheSPHparticlesitaccreted,asimplicitinEqn.20,
j
γ (R3/GM )1/2, j is immediately accreted – in its total- it would spin so fast that it could not reach stellar densi-
s s s
ity – by s. γ is a tolerance parameter, with default value ties. This is what happens with all standard sinks in SPH:
s
0.01.Thisdeviceisintendedtomoderatethesituationwhere they are sinks of both mass and angular momentum, and
theSPHparticlesinsidetheinteraction-zonehaveformed a thisishighlyunrealistic.Inreality,ifthematerialinflowing
Toomreunstabledisc.Theassumptionisthat thiswill lead towardsaprotostarhashighspecificangularmomentum,it
toefficienttransportofangularmomentumbygravitational firstfallsontoanaccretiondisc,andthenspiralsinwardsby
torques, and hence veryrapid inspiral ontothe point-mass. transferringthebulkofitsangularmomentumtootherma-
Thispiecemealleachingofmass,fromtheSPHparticles terial further out in the disc or envelope. Therefore, at the
in theinteraction-zone, ontothepoint-massat itscentre,is end of each timestep, we reduce the angular momentum of
termed regulated accretion, and is a critical feature of the the point-mass by transferring some of its angular momen-
NewSink algorithm. It ensures that, as long as the sink is tum to the SPH particles in the surrounding interaction-
accreting, there are SPH particles in the interaction-zone, zone.Sincethis transfer ofangular momentum is presumed
and therefore there are no very steep gradients across the to be effected by viscous torques in an accretion disc, the
boundary of thesink. amount of angular momentum transferred is given by
The idea of transferring mass piecemeal from an SPH
δt
particle to a sink is not new. This device has been im- |δLs| = |Ls| 1−exp −ht s i . (21)
plemented by Anzeret al. (1987), in simulating the accre- (cid:26) (cid:18) DISC (cid:19)(cid:27)
tion of a supersonic wind by a neutron star, and has sub- We achieve this by giving each SPH particle, j, in the
sequently been used in the same context by, for example, interaction-zone of s an impulse of velocity
Boffin & Anzer (1994). In their algorithm, mass is trans-
|δL |L ×∆r
ferred to a pre-existing sink, from all the particles in the δvj = s s js . (22)
computational domain – but at a rate that is strongly j{mj∆rjs×Ls×∆rjs}
weighted towards those presentlynear thesink; theweight- (cid:12) (cid:12)
Thuseach(cid:12)SPPHparticleintheinterac(cid:12)tion-zoneofsreceives
ing function has to be tuned to produce an acceptable rate (cid:12) (cid:12)
an impulse of velocity proportional to its distance from the
of accretion onto the sink. However, here we are consider-
rotation axis of thepoint-mass.
ing dynamically created sinks and subsonic ambient flows.
In order to compensate for these impulses of velocity,
Moreover, in our algorithm, the particles that leach mass
thepoint-massreceivesimpulsesofmomentumandangular
to the point-mass are always the ones near the centre of
momentum given by
theinteraction-zone,whichneverleavetheinteraction-zone;
consequently,weonlyhavetocontendwith asmallnumber δv = −M−1 {m δv } , (23)
s s j j
of particles having time-varying mass, and none outside of
j
X
theinteraction-zone.
δL = − {m ∆r ×δv } . (24)
s j js j
j
X
As a consequence, the net angular momentum (invested in
SPHparticlemotionsandthespinsofpoint-masses)isaccu-
2.4 Updating the properties of a NewSink
rately conserved. At the same time, the amount of angular
momentum invested in the spins of point-masses is small,
Attheendofeachtimestep,δt ,themass,M ,position,r ,
s s s
both because it is continually transferred to the SPH par-
velocity,v , and angular momentum,L ,of thepoint-mass
s s
ticles in theinteraction-zone, and because theclose-in SPH
are updated to takeaccount of accretion,
particles from which it assimilates mass have normally had
M′ = M + {δm } , (17) tolosealotofangularmomentumtogetclose-ininthefirst
s s j
place. This is a critical feature of the NewSink algorithm.
j
X
6 Hubber, Walch and Whitworth
It ensures that NewSinks do not act as sinks for angular First, it must haveentered theexclusion-zone of s, i.e.
momentum,andthatthepoint-masscanonlygrowinmass
|r −r | < R . (25)
at a rate influenced by how fast the close-in SPH particles j s s
are able to transfer their angular momentum to other SPH Second, themutualnon-thermalenergy of j and s (i.e.
particles furtherout in the interaction-zone and beyond. theirmutual bulk-kineticplus gravitational energy),
m M |∆v |2 Gm M |∆r | |∆r |
e = j s js + j s φ js +φ js ,
js 2(m +M ) 2 h H
3 THE OldSink ALGORITHM j s (cid:26) (cid:18) j (cid:19) (cid:18) s (cid:19)(cid:27)
(26)
Many variants of standard sink have been used previously
in SPH simulations of star formation (e.g. Bate et al. 1995; must be negative. If more than one OldSink is minded to
Federrath et al. 2010; Wadsley et al. 2011), and it is not accretej,j isaccreted bytheonewhosepoint-massisclos-
practical to present test results for all of them. Here, and est. When an SPH particle i triggers the formation of an
in Section 4, we describe the particular representative im- OldSink, s, these conditions are normally satisfied by all,
plementationsthatweuseforcomparisontests.Thechoices ormost,oftheneighboursofi,sosimmediatelyassimilates
havebeen madewith aviewtodistinguishing theproblems all, or most, of theneighbours of i.
associatedwithsinkcreationfromthoseassociatedwithsink Bate et al. (1995) use two additional criteria to deter-
evolution.InthisSectionwedefineOldSinks,whicharecre- mine whether an SPH particle, j is assimilated by a stan-
atedinthesamewayasNewSinks,andthereforeavoidthe dard sink, s. First, the angular momentum of j relative to
problemsnormally associated with thecreation ofstandard thepoint-massofsshouldbelessthanitwouldbeifj were
sinks.Thisenablesustoisolatetheproblemsassociatedwith in a circular orbit of radius Rs about thepoint-mass, i.e.
theevolution of standard sinks,viz. that theyall appearto |∆r ×∆v | < (GM R )1/2 . (27)
js js s s
seriously corrupt the hydrodynamics near the sink bound-
ary (by introducing excessively steep gradients), to act as Second, if an SPH particle j passes within a small distance
sinks for angular momentum, and to produce results that (0.1Rs or less) of the point-mass, all other criteria are ig-
are strongly dependent on the user-defined parameters of nored,anditisassimilated bythepoint-mass.Wefindthat
sink creation and evolution. thesecriteria do not have asignificant effect on theresults.
Bate et al. (1995) have also developed a procedure for
extrapolating the density and velocity fields outside a sink,
3.1 Internal structure and advection of an
and using this information to generate correction terms to
OldSink
account for the missing SPH particles inside the exclusion-
An OldSink comprises a central point mass (the point- zone. However, these terms have not been included in sub-
mass)andacomovingconcentricsphericalvolumeofradius sequentimplementationsoftheBate et al.(1995)algorithm
R (hereafter the exclusion-zone). The mass, M , position, (Bate, private communication), and we also do not include
s s
r ,velocity,v ,andangularmomentum,L ,ofanOldSink them in OldSinks(or UrSinks,see Section 4).
s s s
refer to its point-mass. The exclusion-zone, as its name im-
plies, generally contains very few SPH particles – and this
3.4 Updating the properties of an OldSink
isoneofthecriticalfeaturesthatdistinguishesanOldSink
fromaNewSink.Oncecreated,anOldSinkistrackedwith Accretion of an SPH particle by an OldSink, s, is imme-
thesame integration routine as an SPH particle. diate and complete; the SPH particle is removed from the
simulation,andthepropertiesofsareupdatedaccordingto
3.2 Criteria for creation of an OldSink M′ = M + {m } , (28)
s s j
The criteria for creation of an OldSink are the same as Xj
for creation of a NewSink, viz. that an SPH particle, i,
r′ = M′−1 M r + {m r } , (29)
has density exceeding ρ (Eqn. 2), that the resulting s s s s j j
SINK !
sink does not – at its creation – overlap a pre-existing sink Xj
(Eqn. 3), that i has lower gravitational potential than all
v′ = M′−1 M v + {m v } , (30)
its neighbours (Eqn. 4), and that i satisfy the Hill Sphere s s s s j j
!
j
criterion(Eqn.5).IfanSPHparticleisatisfiesthesefourcri- X
teria, an OldSinkiscreated with thesame mass, Ms=mi, L′s = Ls + Ms∆rss′×∆vss′
position, rs = ri, velocity, vs = vi, and angular momen- + {mj∆rjs′×∆vjs′} ; (31)
tum, L =0, as i. The radius of the exclusion-zone is set
s j
X
to R = X h , where X is a user-defined param-
s SINK i SINK theprimed variables represent the new properties of s, and
eter. With the M4 smoothing kernel, the default value is
the sums are over all the SPH particles j accreted during
X = 2, so that the exclusion-zone is the same as the
SINK that timestep. We term this instantaneous accretion.
smoothing volumeof i.
We note that, if an OldSink grows by disc accretion,
most of the SPH particles that it accretes only just satisfy
3.3 Criteria for accretion by an OldSink Eqn.25, and so its specific angular momentum is
An SPH particle j is assimilated by the point-mass of an |L | M 1/2
s ∼ 1020cm2s−1 s
OldSink, s,if it satisfies two criteria. M M
s (cid:18) ⊙(cid:19)
An improved sink particle algorithm for SPH simulations 7
m 1/6 ρ −1/6 4.2.3 ∇·v criterion
× SPH SINK . (32)
10−5M 10−13gcm−3
(cid:18) ⊙(cid:19) (cid:18) (cid:19) ThefourthcriterionforSPHparticleitotriggerthecreation
This is so large that the point-mass can not condense of an UrSink is that the divergence of the velocity at the
to stellar densities (for example, the Sun spinning at position of i should be negative,
break-up speed only has specific angular momentum ∼
(∇·v) < 0, (34)
(GM R )1/2∼3×1018cm2s−1),unlessm isverysmall i
⊙ ⊙ SPH
(impractically high mass-resolution) and/or ρ is very otherwise the density there must be decreasing. This crite-
SINK
high (so high as to negate theadvantages of using sinks). rion is used by Wadsley et al. (2011).
4.2.4 Non-thermal energy criterion
ThefifthandfinalcriterionforSPHparticleitotriggerthe
4 THE UrSink ALGORITHM
creationofanUrSinkisthatthenetnon-thermalenergyof
In this Section we define UrSinks, which are created and i and its neighbours j, in the centre-of-mass frame, should
evolved in the manner of standard sinks, and are represen- benegative,
tative of the sinks that have been used in almost all extant
E = E +E < 0. (35)
SPH simulations of star formation. NT GRAV KIN
Here,
4.1 Internal structure and advection of an UrSink EGRAV = G4Xj jX′6=jmjmj′(cid:26)φ(cid:18)|∆hrjjj′|(cid:19)+φ(cid:18)|∆hrjj′j′|(cid:19)(cid:27)
AnUrSinkhasthesameinternalstructureasanOldSink, (36)
viz.acentralpoint-masssurroundedbyacomovingconcen- is theself-gravitational potential energy,
tric spherical exclusion-zone of radius R . The point-mass
s 1
has mass, Ms, position, rs, velocity, vs and angular mo- EKIN = 2 mj|vj−vj|2 (37)
mentum Ls, and is advected likean SPH particle. Xj (cid:8) (cid:9)
is the net kinetic energy, and v = {m v }/ {m } is
j j j j
thecentre-of-massvelocity.Hereallthesumsovertheneigh-
P P
4.2 Criteria for creation of an UrSink bours j include i itself. This criterion is somewhat weaker
than the energy criteria used by other standard algorithms
Some of thecriteria for creation of an UrSinkare different (see Section 4.2.5 below).
from those used for NewSinks and OldSinks. If an SPH
particleisatisfiesthefivecriteriadetailedbelow,anUrSink
is created with the same mass, M =m , position, r =r , 4.2.5 Alternative creation criteria
s i s i
velocity, v = v , and angular momentum, L = 0, as i.
s i s In this Section, we list – for completeness – some of the
The radius of the exclusion-zone is set to R = X h ,
s SINK i additional criteria used in other standard sink creation al-
and if the M4 smoothing kernel is used, the default value gorithms, butnot for UrSinks.
is X =2, so that the exclusion-zone is the same as the
SINK Bate et al. (1995) and Wadsley et al. (2011) apply ad-
smoothing volumeof i.
ditional energy criteria for creating a standard sink, viz.
1
E + E < 0, (38)
THERM 2 GRAV
4.2.1 Density and overlap criteria E +E +E < 0, (39)
ROT THERM GRAV
E +E +E < 0. (40)
The first two criteria for SPH particle i to trigger the cre- KIN THERM GRAV
ation of an UrSink are the same as for NewSinks and Here
OldSinks,viz.thatitsdensityshouldexceedauser-defined
m a2
threshold,andthattheresultingUrSinkshouldnotoverlap E = j j (41)
a pre-existing UrSink(see Sections 2.2.1 and 2.2.2). THERM j (cid:26)(γj−1)(cid:27)
X
isthenetthermalenergy(a andγ aretheisothermalsound
j j
speed and ratio of specific heats for SPH particle j),
4.2.2 ∇·a criterion |L |4
E = INT (42)
ThethirdcriterionforSPHparticleitotriggerthecreation ROT 2 j{mj|∆rjs.LINT|2}
of an UrSink is that the divergence of the acceleration at
is thenet rotaPtional energy,
theposition of i should benegative,
L = {m (r −r )×(v −v )} (43)
(∇·a) < 0, (33) INT j j j j j
i j
X
otherwise the possibility exists that the dense gas is about isthenetangularmomentum,andr = {m r }/ {m }
j j j j
tobetornaparttidally.ThiscriterionisusedbyBate et al. isthecentreofmass.Allthesumsovertheneighboursj in-
P P
(1995) and Wadsley et al. (2011). clude i itself. Although these energy criteria are somewhat
8 Hubber, Walch and Whitworth
stifferthanEqn.(35),theydonotproducesignificantlydif-
ferent reaults.
Federrath et al. (2010) and Wadsley et al. (2011) also
use the potential-minimum criterion (i.e. Eqn. 4), and in
addition theyadvocatecheckingthattheflowisconvergent
in all directions, by computing the eigenvalues of dv /dx
i j
(here i and j are identifiers for the Cartesian coordinates)
and requiring them to be individually negative. However,
thesecriteriahavenotbeenimplementedinthemajority of
extant SPH simulations of star formation.
4.3 Criteria for accretion by an UrSink
AnSPHparticlejisimmediatelyaccretedbythepoint-mass
ofanUrSinks,ifitfallswithintheexclusion-zoneofs,and
their mutual non-thermal energy is negative, just as for an
OldSink (Eqns. 25 & 26). Under normal circumstances, a Figure 1. The normalised accretion rate (M˙s/M˙BO) as a func-
newly-created UrSink immediately assimilates most of the tion of the normalised sink radius (Rs/RSONIC) for isothermal
neighbours of the SPH particle that triggered its creation. Bondi accretion. Results obtained with NewSinks are shown as
reddiamondsconnected byasolidline,andthoseobtainedwith
OldSinksareshownasblackstarsconnected byadashedline.
4.4 Updating the properties of an UrSink
Each time an UrSink assimilates one or more SPH parti-
neglectingtheself-gravityoftheinflowinggas.Acriticalrole
cles,thepropertiesoftheUrSinkareupdatedaccordingto
is played by thesonic radius at
Eqns. (28) through (31), just as for an OldSink. As with
an OldSink, this often leads to the unphysical acquisition R = GM⋆ . (45)
of large amounts of angular momentum. SONIC 2a2
O
Forr≫R theinwardflowissubsonicandthepressure
SONIC
acceleration(−ρ−1∇P)playsanimportantrole.Conversely,
5 TESTS forr≪R theinwardflowapproachesfreefallandthe
SONIC
pressure acceleration is unimportant.
WehaveusedthefollowingfourteststocompareNewSinks
Thedensityandvelocityprofilesforthetranssonicsolu-
with OldSinks and UrSinks, and to explore the effect of
tionareobtainedbysolvingtheBernouilliequationnumeri-
changing the user-prescribed parameters for sink creation
cally.Theinitialconditionsaresetupbyrelaxingaperiodic
(ρ and X )andtheSPHresolution (N ):isother-
SINK SINK SPH cube to produce a uniform-density glass, cutting a sphere
malBondiaccretion(Section5.1),thecollapseofarotating
containing 5×105 particles from this cube, stretching the
Bonnor-Ebertsphere(5.2),theBoss-Bodenheimertest(5.3),
spheretoproducetherequireddensityfield,andgivingeach
and the fragmentation of a turbulent prestellar core (5.4).
SPHparticletheinwardradialvelocityappropriatetoitspo-
Parametervaluesusedforthedifferenttestsaresummarised
sition. The outer boundary of the sphere is at ∼20R ,
inTable2.Althoughallfourtestshavebeenperformedwith SONIC
andthesimulations are followed totimet =2GM /a3.
allthreetypesofsink,wedonotpresentanddiscussresults END ⋆ O
This test does not involve sink creation, and therefore
obtained with UrSinks for the first three tests; this is be-
wedonotspecifyρ orX .Instead,asinkisplacedat
cause the results obtained with UrSinks are very chaotic SINK SINK
theoriginfromtheoutset.Itspoint-masshasmassM =M
andshownosignofconverging;theturbulentprestellarcore s ⊙
and,althoughitaccretesSPHparticlesandhasanaccretion
test suffices to demonstrate this.
rate,M˙ ,itsmassisheldfixed.Thisisself-consistent,since,
All the tests are performed with the seren SPH code, s
if the accreting gas is sufficiently rarefied to havenegligible
usingitsdefaultoptions(Hubberet al.2011).Wesetη=1.2,
self-gravity (compared with thegravity of thecentralstar),
so that the smoothing-length of SPH particle i is h =
i the increase in the central mass by t is negligible. The
1.2(m /ρ )1/3 and the mean number of neighbours is END
N¯ SP∼H58.i purpose of this test is to evaluate the treatment of radial
NEIB accretion.
ThesystemisevolvedwithaNewSinkoranOldSink,
5.1 Bondi accretion with the sink radius set to different multiples of the sonic
radius, R /R =1/8, 1/4, 1/2, 1, 2, 4, 8. By t , the
s SONIC END
ThefirsttestmodelsBondiaccretion.Theaccretionratefor accretion rateis steady,andtherarefaction wavepropagat-
thetranssonic isothermal solution (Bondi 1952) is ing in from the outer boundary is still in the outer parts
of the computational domain. Fig. 1 shows how the accre-
e3/2πG2M2ρ
M˙ = ⋆ O , (44) tion rates depend on R /R . Since this test does not
BO a3 s SONIC
O involve sink creation, the results obtained with an UrSink
where e is the base of natural logarithms, M is the mass are exactly the same as those obtained with an OldSink.
⋆
of the central star, ρ is the density at infinity, and a As long as R < R , both the NewSink and the
O O s SONIC
is the isothermal sound speed. The solution is obtained by OldSinkmodelBondiaccretionwell,givingaccretionrates
An improved sink particle algorithm for SPH simulations 9
Test−→ Default BondiAccretion RotatingBonnor-Ebert Boss-Bodenheimer TurbulentCore
ρ /(gcm−3) 10−11 — 10−11 10−13,10−12,10−11 10−11,10−10,10−9
SINK
X 2 — 2 2,4 2
SINK
X 4 — 4 4 4
HILL
α 0.01 0.01 0.01 0.01 0.01
SS
Rs/RSONIC — 81, 41, 21,1,2,4,8 — — —
NewSink ⋆ ⋆ ⋆ ⋆
OldSink ⋆ ⋆ ⋆ ⋆
UrSink (⋆) (⋆) (⋆) ⋆
Table 2. Default parameter values and values used for the different tests. Row 1 gives the name of the test. Rows 2 to 6 give the
parametervalues,androws7to9indicatewhichtypesofsinkweretested, withbrackets indicatingthetests thatarenotdiscussed.
∼10% below the analytic value. Under this circumstance, 103AU. The initial conditions are then generated by relax-
the material flowing into the sink is close to free fall and ing a periodic cube of SPH particles to produce a uniform
so the pressure acceleration is unimportant; the fact that glass, cutting a sphere containing 105 particles from this
OldSinksgiverisetoinaccurateevaluationsofthepressure cube, and stretching the particle positions to produce this
and viscous forces nearR is then of little consequence. density profile. The particles start from rest.
s
However,forR >R ,anOldSinkseriously over- Thecloudcollapsestoformasinglesinksurroundedby
s SONIC
estimates the accretion rate. Under this circumstance, the asmallaccretiondisc.Fig.2shows,foralltheSPHparticles
pressure acceleration is an important factor controlling the within∼20AUofthepoint-mass,thetangentialviscousac-
flow of material into the sink, and the steep outward pres- celeration, a ,theradial hydrodynamicacceleration, a ,
TV RH
suregradientattheedgeofanOldSinkartificiallyincreases andthedensity,ρ,allasfunctionsofdistancefromthepoint-
the accretion rate; for R /R =8, M˙ is ten times the mass,shortlyaftersinkcreation.Theupperrowgivesresults
s SONIC s
analyticvalue.Incontrast,withaNewSinkthereisasmall obtained with a NewSink, showing that the SPH particles
increase in the accretion rate with increasing R , but, even near R (=2AU) have smoothly varying a , a and ρ,
s s TV RH
forR /R =8, M˙ iswithin∼1%oftheanalyticvalue. duetothepresenceofSPHparticlesintheinteraction-zone;
s SONIC s
Although this test only deals with one possible accre- inparticular,a ispositive(i.e.outward,thereforeresisting
RH
tion mode, all spherical accretion modes are likely to suffer accretion,asitshould).Thelowerrowgivesresultsobtained
fromthesameproblem,iftheinflowissubsonicnearthesink withanOldSink,showingthatSPHparticlesnearthesink
boundary.Thisisaseriousconcerninsimulationsofstarfor- boundaryexperienceartificiallyreduceda ,a andρ,due
TV RH
mationthatsetρ high,inordertofollowprotostarswell to a lack of neighbours in the exclusion-zone; all three fac-
SINK
into their Kelvin-Helmholtz contraction phase. The flow of torsacttoincreasethenetinwardacceleration,andthereby
gas intothesink isthen subsonic,and an OldSinkwill ex- toartificially increase therateof accretion. Inaddition,be-
perienceartificiallyenhancedaccretion.ANewSinkshould causetheSPHparticlesaccretedbyanOldSinkcarrytheir
capturethe accretion rate muchmore accurately. angular momentum with them, rather than transferring it
to the SPH particles a bit further out, these SPH particles
canmorerapidlymoveinward,tobeaccreted intheirturn.
5.2 Rotating Bonnor-Ebert sphere When this test is performed with UrSinks, the results are
chaotic,andinvolvethecreationofmanydifferentsinks;we
Themainmodeofaccretioninstarformationinvolvesdiscs.
donot discuss these results.
Theaccretionrateisthencontrolledbythetorquesthatre-
Fig. 3 shows the evolution of the sink mass, M , and
distribute angular momentum. There are no analytic solu- s
thesinkaccretionrate,M˙ ,asfunctionsofthetimeelapsed
tions for this mode, so testing is less straightforward. The s
since sink formation, showing that the artificially increased
test we perform involves a rigidly rotating spherical gas
accretionobtainedwithanOldSinkisnotatransienteffect.
cloudhavingmassM ,initialboundaryradius103AU,and
⊙ Bytheendofthesimulation,whichcorrespondsto∼100or-
initialangularspeed3.54×10−12s−1.Thegasisandremains
bital periods at the sink boundary, the OldSink has M ≃
isothermal with sound speed 0.20kms−1 (corresponding to s
0.37M , whereas the NewSink only has M ≃ 0.24M .
molecular gas at T =11K). However, the main purpose of ⊙ s ⊙
ThusanumericalsimulationusingOldSinksoverestimates
thistestisnottoexploresinkcreationfromisothermalgas,
accretion rates. The effect on final protostellar masses is
but ratherto evaluate thetreatment of disc accretion.
hard to call. An enhanced accretion rate would mean en-
The initial density profile is that of a critical Bonnor-
hancedradiativeandmechanicalfeedback (ifthesewerein-
Ebert Sphere, i.e. one with dimensionless boundary radius
cludedinthesimulation),andthismightactuallyterminate
ξ =6.45, but the cloud is not in equilibrium. The initial
B accretion earlier, leading to a lower final mass.
ratios of thermal and rotational energy to gravitational en-
ergy are α=0.09 and β=0.59. The cloud is modelled with
105 SPH particles (hence m = 10−5M ), and we set
SPH ⊙ 5.3 Boss-Bodenheimer test
ρ = 10−11gcm−3 and X = 2. The initial density
SINK SINK
profile is obtained by integrating the Isothermal Equation The Boss-Bodenheimer test (Boss & Bodenheimer 1979) is
(e.g. Chandrasekhar & Wares 1949) from ξ=0 to ξ=6.45, a standard test, applied to many star formation codes.
and scaling the result to give total mass M and radius For example, Bate et al. (1995) used it to test the original
⊙
10 Hubber, Walch and Whitworth
Figure 2.(a,d) Tangential viscous accelerations, (b,e) radialhydrodynamic accelerations, and (c,f)densities, for allthe SPH particles
within∼20AUofthepoint-mass,shortlyaftersinkcreationintheRotatingBonnor-EbertSphereTest.Theresultsonthetoprowhave
beenobtainedwithaNewSink,andthoseonthebottom rowwithanOldSink.
Figure 3.(a)Thesinkmass,and(b)thesinkaccretion rate,intheRotatingBonnor-EbertSphereTest,asafunctionoftimeelapsed
sincetheformationofthesinkatt0,usingaNewSink(redsolidlines)andanOldSink(blackdashedlines).
sink algorithm, although of necessity they could only use a cloudhavingmassM ,radius2×103AU,anddensityfield
⊙
rather small number of SPH particles (N ≃ 8×103);
SPH
Bate & Burkert (1997) repeated the test with ten times as ρ(r,θ,φ) = 1.6×10−17gcm−3 {1 + 0.5cos(2φ)} , (46)
manySPHparticles(N ≃8×104),butdidnotusesinks.
SPH where(r,θ,φ)are spherical polar coordinates. Itis in solid-
The version we consider here starts with a spherical body rotation with angular speed Ω = 1.56×10−12s−1,
O