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 computational 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 implementation of sink particles in AMR (Krumholzet al. 2004), a 2 THE NewSink ALGORITHM sink is introduced wherever the local Jeans length falls below 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ρ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 formation 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 boundary (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 energy 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

An improved sink particle algorithm for SPH simulations

Similar An improved sink particle algorithm for SPH simulations