Radiation-hydrodynamical simulations of massive star formation using Monte Carlo radiative transfer: I. Algorithms and numerical methods

1.6 MB·English

by  Tim J. Harries

Mon.Not.R.Astron.Soc.000,1–11(2015) Printed26January2015 (MNLATEXstylefilev2.2)
Radiation-hydrodynamical simulations of massive star
formation using Monte Carlo radiative transfer
I. Algorithms and numerical methods
Tim J. Harries(cid:63)
Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, United Kingdom
5
1
0
2 ABSTRACT
We present a set of new numerical methods that are relevant to calculating radiation
n
pressuretermsinhydrodynamicscalculations,withaparticularfocusonmassivestar
a
J formation.TheradiationforceisdeterminedfromaMonteCarloestimatorandenables
3 a complete treatment of the detailed microphysics, including polychromatic radiation
2 andanisotropicscattering,inboththefree-streamingandoptically-thicklimits.Since
the new method is computationally demanding we have developed two new methods

] that speed up the algorithm. The first is a photon packet splitting algorithm that
R
enables efficient treatment of the Monte Carlo process in very optically thick regions.
S The second is a parallelisation method that distributes the Monte Carlo workload
. overmanyinstancesofthehydrodynamicdomain,resultinginexcellentscalingofthe
h
p radiation step. We also describe the implementation of a sink particle method that
- enables us to follow the accretion onto, and the growth of, the protostars. We detail
o the results of extensive testing and benchmarking of the new algorithms.
r
t
s Key words: stars: formation – methods: numerical
a
[

1
v
1 INTRODUCTION tiple dimensions, and the introduction of angular momen-
4
tum, leads to a scenario in which the protostar can accrete
5 Massive stars are hugely important in galactic ecology, en-
via a disc, which intercepts a smaller fraction of the pro-
7 riching them chemically and providing strong feedback ef-
tostellar radiation field, and the radiation escapes preferen-
5 fects via radiation, stellar winds and supernovae. They also
0 tially via low-density cavities along the rotation axis (e.g.
play a pivotal role in measuring star formation rates in dis-
. Yorke&Sonnhalter2002).RecentworkbyKrumholzetal.
1 tant galaxies, and therefore determining the star formation
(2009) used 3-D adaptive mesh refinement (AMR) simula-
0 historyoftheUniverse(Kennicutt1998).However,thepro-
tionscoupledwithaflux-limiteddiffusion(FLD)approxima-
5 cessbywhichmassivestarsformislesswellunderstoodthan
tionfortheradiationtransport(RT)tosimulatetheforma-
1 thatoflower-massstars,bothobservationally,becausemas-
: tionofamassivebinary.Theyfoundaccretionoccurredvia
v siveprotostarsarerareanddifficulttoobserve,andtheoret- aradiatively-drivenRayleigh-Taylorinstabilitythatallowed
i icallybecauseradiationfeedbackhasamuchstrongerinflu-
X materialtopenetratetheradiation-drivencavitiesaboveand
enceonthegasdynamicsthanitdoesforsolar-massobjects.
below the protostars. Their simulation resulted in the for-
r Broadly it appears that observational and theoretical per-
a spectives on massive star formation are converging towards mation of a massive (41M(cid:12)and 29M(cid:12)) binary pair.
ascenarioinwhichamassiveyoungstellarobjectcontainsa Kuiper et al. (2010) implemented a hybrid RT scheme
central (>10M ) protostar surrounded by a stochastically in which ray tracing is used to calculate the radiation field
(cid:12)
accretingKepleriandiscandastellar-ordisc-windoutflow, fromtheprotostellarphotosphereouttowherethematerial
althoughsubstantialconflictsbetweenthepredictionsofnu- becomes optically thick, at which point the ray-traced field
merical models remain (Tan et al. 2014). becomes a source term in the grey diffusion approximation
It has been known for some time that massive proto- used for the rest of the domain. The method shows good
starsundergoingsphericalaccretionwilleventuallyproduce agreementwithmoredetailedtreatmentsoftheRTbuthas
sufficient radiation pressure to halt the inflow of material vastlyreducedcomputationalcost(Kuiper&Klessen2013).
(e.g. Wolfire & Cassinelli 1987). Extending models to mul- Two-dimensionalcalculationsadoptingthismethodresulted
intheformationofmassivestarsviastochasticdiscaccretion
accompanied by strong radiatively driven outflows (Kuiper
(cid:63) E-mail:[email protected] et al. 2010).
(cid:13)c 2015RAS
2 Tim J. Harries
Radiatively-driven Rayleigh-Taylor instabilities have els and observations across a wide range of continuum (e.g.
not been identified in the hybrid RT models (Kuiper near-IRdust,radiofree-free)andline(e.g.forbidden,molec-
et al. 2012), suggesting a fundamental difference with the ularorrecombination)diagnostics,(Rundleetal.2010;Ha-
Krumholzetal.computations.ItappearsthatthegreyFLD worth et al. 2012, 2013).
methodconsiderablyunderestimatesthedrivingforceofthe In this paper we describe the implementation of (i) a
radiation, since the opacity used for the radiation pressure new,highlyparallelisedmethodforcalculatingtheradiation
calculation is the Rosseland opacity at the dust tempera- pressureusingMonteCarloestimators,(ii)apacketsplitting
ture, whereas in the hybrid case the primary momentum method that improves the efficiency of the method at high
deposition occurs where the radiation from the protostar is opticaldepths,and(iii)aLagrangiansink-particlealgorithm
incident on the dust (near the dust sublimation radius). At that enables us to follow the growth of the protostar. We
thispointtheradiationtemperatureismuchhigher,andthe provide the results of extensive benchmarking of the new
dustismuchmoreopaque.ItispossiblethatthegreyFLD methods, with a view to conducting simulations of massive
method underestimates the the driving force by up to two starformationthatincorporatebothradiationpressureand
orders of magnitude (Kuiper et al. 2012). ionisation feedback.
The fundamental point here is that the dynamics of
the simulations can be critically affected by the level of mi-
crophysical detail employed. The lessons of recent numeri- 2 RADIATION PRESSURE
cal simulations are that disc accretion may occur, and fast,
Our Monte Carlo method is closely based on that of Lucy
wide-angle outflows may be driven. It seems that radiation
(1999), in which the luminosity L of the illuminating ob-
pressuredoesnotpresentafundamentalphysicalbarrierto
ject is divided up into N discrete photon packets that have
massivestarformationviaaccretion,althoughthestrength,
a constant energy as they propagate through the adaptive
geometry, and longevity of the outflows is still contentious.
mesh. The energy of each packet is
Furthermore, neither the Krumholz et al. (2009) models
nor those of the Kuiper et al. (2010) include the effects of (cid:15) =w L∆t (1)
i i
photoionization, which will come into play as the protostar
where∆tisthedurationoftheMonteCarloexperimentand
shrinks towards the main sequence and its radiation field
w is a weighting factor, normalized so that
hardens. All these factors point towards the necessity for i
a more thorough treatment of the microphysical processes N
(cid:88)
that underpin that radiation feedback, such as a polychro- wi =1. (2)
matic prescription for the radiation field, consideration of i=1
both dust and gas opacity (absorption and scattering), and Note that the original Lucy method has w =1/N, so that
i
the inclusion of photoionization. each photon packet had the same energy.
Recently significant progress has been made in Asthepacketspropagatethroughthegridtheymayun-
using MC transport in RHD codes. For example dergoabsorptionorscatteringevents.Duringanabsorption
Nayakshin et al. (2009) used MC photon packets event the photon packet is immediately re-emitted, with a
to treat radiation pressure in smoothed-particle hy- newfrequencyselectedatrandomfromaprobabilitydensity
drodynamics. Acreman et al. (2010) presented a function created from the appropriate emissivity spectrum
proof-of-conceptcalculationinvolvingcombining3D at that point in the grid. The radiation field therefore re-
smooth-particle hydrodynamics with MC radiative mains divergence free during the calculation. After the N
equilibrium. In this case the SPH particle distribu- packets have been propagated, the energy density for each
tion was used to construct an AMR grid on which cell can be straightforwardly computed and hence the ab-
the MC radiative-equilibrium calculation was con- sorption rate can be determined.
ducted. We then developed an Eulerian hydrody- Thecalculationoftheradiationpressuremaybeconve-
namics module to perform 3D hydrodynamics on niently conducted in parallel with the radiative equilibrium
the native AMR grid used for the radiation transfer calculation. The simplest method is to assign a momentum
(Haworth & Harries 2012). A similar approach was change to each cell (∆p ), which is zeroed at the start of
j
adopted by Noebauer et al. (2012), who incorpo- each MC loop. The each time a photon packet enters and
rated a radiation-pressure force estimator as well as then leaves a cell the difference in momentum between the
a radiative-equilibrium calculation. Roth & Kasen incoming and outgoing momenta is added to the cell’s mo-
(2014) also coupled MC transport and 1D hydro- mentum change
dynamics, incorporating special relativity and reso- (cid:15)
nance line scattering. ∆pj =∆pj+ ci (uˆin−uˆout) (3)
We have now developed an AMR RHD code (torus)
whereuˆ isthedirectionvectorofthephotonpacketwhenit
that employs a highly-parallelised, time-dependent, Monte- in
entersthecellanduˆ isthedirectionvectorofthephoton
Carlo (MC) method (Harries 2011) to follow the RT at a out
packetwhenitexitsthecell.AttheendoftheMCloopwe
levelofmicrophysicaldetailcomparabletothatofdedicated
can then calculate the radiation-force per unit volume for
RT codes such as Cloudy (Ferland et al. 2013). The conse-
each cell via
quences of this breakthrough are two-fold: firstly, the in-
∆p
creasedsophisticationofthetreatmentoftheradiationfield f = j (4)
rad,j ∆tV
can have a significant impact on the hydrodynamics (Ha- j
worth & Harries 2012), and secondly it becomes possible whereV isthecellvolume.Wewillcallthisthemomentum
j
to make a much more direct comparison between the mod- tracking algorithm.
(cid:13)c 2015RAS,MNRAS000,1–11
Massive star formation simulation algorithms 3
This method has the advantage that it is straightfor- 2.1 Radiation pressure tests
ward to implement, is transparently related to the underly-
We constructed a suite of tests of our radiation pressure
ing physics, is valid for arbitrary scattering phase matrices,
algorithms based on a uniform-density sphere of radius
andisveryfasttocalculate.Itdoeshoweverhavethedisad-
R = 0.1pc and mass 100M illuminated by a solar-
vantagethatthenumberofscatterings/absorptionspercell s (cid:12)
type star. The grey opacity κ was selected in order that
tends to zero in the optically thin limit. This is an analo-
τ = κρR = 0.1,1,10 or 100. The models were run with a
gousproblemtothatassociatedwiththeBjorkman&Wood s
1-Dsphericalgeometryonauniformradialmeshcomprising
(2001) radiative-equilibrium method.
1024 cells. We ran models for pure absorption (α=0) and
An alternative method is to use a MC estimator of the
purescattering(α=1)cases,andforthescatteringmodels
radiation vector flux to calculate the radiation force. We
we further subdivided the models into isotropic scattering
define the energy density u as the energy density per unit
ν and forward scattering (using a Henyey-Greenstein phase
frequency of a beam of radiation. The energy carried in a
function with g = 0.9). We used 105 photon packets for
beam of volume dV is
all the models. For each model we calculated the radiation
dE =u dVdν (5) forceperunitvolume(Fκρ/c)forboththeMonteCarloflux
ν
estimator and the momentum tracking algorithm.
and since dV =cdAdt we get
dE =uνcdAdtdν =IνdtdνdAdΩ (6) 2.1.1 Pure absorption case
where Iν is the specific intensity. Hence we have In this case the radiation force per unit volume (frad) as a
function of radius r in the sphere is expected to be
I
u = νdΩ. (7)
ν c Lκρexp(−τ(r))
f = (13)
rad 4πr2c
If a photon packet traverses length (cid:96) between successive
events(eventsbeingabsorptions/scatteringsorcrossingcell where τ(r) is the radial optical depth from the centre of
boundaries) then the energy density due to this part of the the grid to r. Since the photon packets stop propagating
path of packet (cid:15) is ((cid:15) /V )δt/∆t where δt = (cid:96)/c. We can whentheyareabsorbedthenumberofpacketspassingeach
i i j
therefore obtain an MC estimator of the energy density cell will be monotonically declining as a function of radius,
and we therefore expect the noise on the estimates of the
(cid:15) (cid:96)
u dν = i (8) radiation force to increase radially. Furthermore we expect
ν Vjc∆t that the noise in the momentum tracking algorithm to be
larger than that of the flux estimator method, particularly
This may be related to the specific intensity via equation 7
for the optically thin cases (τ =0.1,1).
and we get
TheresultsofthisbenchmarkaredisplayedinFigure1
I dΩdν = (cid:15)i(cid:96) . (9) andtheMonteCarloestimatorsshowgoodagreementwith
ν ∆tV theanalyticalresult,withtheexpectednoisedependencies.
j
Thesemodelsranrapidly,andsinceeachphotonpacketonly
The radiation flux vector of cell j may then be written
undergoes one absorption event the higher optical depth
(cid:90) 1 1 1 (cid:88) cases ran the most quickly.
F = I (Ω)dΩ= (cid:15) (cid:96)uˆ (10)
j,ν ν ∆tV dν i
j
dν
2.1.2 Pure scattering cases
where uˆ is the unit vector in the direction of the photon
packetasittraverses(cid:96)andonlythepathsforphotonpackets For isotropic scattering the number of photon interactions
withfrequenciesintherange(ν,ν+dν)areconsidered.The before escape from an optically-thick, uniform spherical of
force per unit volume for cell j is then radial optical depth τ will be
max
1(cid:90) 1 1 1 (cid:88) N ≈τ2 /(cid:104)τ2(cid:105)=τ2 /2. (14)
f = κ ρF dν = (cid:15) κ ρ(cid:96)uˆ. (11) scat max max
rad,j c ν ν c∆tV i ν
j Thusfortheτ =100modeleachphotonpacketwillundergo
whereκ istheopacityatthefrequencyν ofpacket(cid:15) asit ∼5000scatteringeventsbeforeescaping.Sincetheradiation
ν i
traverseslength(cid:96),ρisthemassdensity,andthesummation field is divergence free, the outward flux at every radius is
isnowoverallpackets.Notethattheopacityκ mustconsist a constant F = L/(4πr2) and the radiation force per unit
ν
of both the absorption and scattering opacities volume is then
Lκρ
κi =(1−αν)κν+(1−gν)ανκν (12) frad = 4πr2c. (15)
where α is the albedo, and g = (cid:104)cosθ(cid:105) where θ is the Once again we expect the flux estimator method to have
ν ν
scatteringangle(i.e.g =0forisotropicscattering,positive less noise than the momentum tracking method, and that
ν
for preferentially forward scattering and negative for back thelatteralgorithmwillberatherpoorerinthelowoptical
scattering). Although marginally more expensive to calcu- depth models.
late, this flux estimator method should be superior in the Theresultsoftheisotropic,purescatteringbenchmark
optically-thin limit, since even photon packets that do not are displayed in Figure 2. Excellent agreement is seen be-
interact in a cell contribute to the estimator. tweentheMonteCarloestimatorsandtheanalyticalresult,
(cid:13)c 2015RAS,MNRAS000,1–11
4 Tim J. Harries
Figure 1. Benchmark for absorption (α = 0). The solid lines Figure 3. Benchmarks for a scattering medium (α=1) with a
indicate the analytical results (Fκρ/c), while values calculated forward-scattering Henyey-Greenstein phase function (g = 0.9).
from the Monte Carlo flux estimator are shown as plusses (+) SymbolsarethesameasthoseforFigure1.
and those found from following packet momenta are shown as
crosses (×). Benchmarks are plotted for spheres of four optical
depths(τ =0.1,1,10,100). Lκρ(1−g)
f = . (16)
rad 4πr2c
Theresultsoftheforward-scatteringmodelaregiveninFig-
ure 3 and excellent agreement with the expected analytical
radiation force is found.
2.2 Radiation pressure hydrodynamic test
Satisfiedthattheradiationpressureforceswerebeingprop-
erlycapturedbythecode,weprogressedtotestingthetreat-
mentinadynamicmodelusingasimilarbenchmarktothat
presented by Nayakshin et al. (2009). A uniform density
sphere (ρ = 1.67×10−22gcc−1 and R = 1017cm) was
0 s
modelled using 1-d grid comprising 1024 uniformly spaced
radialcells.Thespherewasilluminatedbya1L starplaced
(cid:12)
atthecentre,andweused105 photonpacketsperradiation
step.
For the purposes of testing we assumed isothermal gas
atatemperatureof10Kandconsideredcellsaboveacritical
densitythreshold(1.67×10−25gcc−1)tobecompletelyop-
ticallythickatallwavelengthswithanalbedoofzero,while
cells with densities below the threshold were considered to
betransparenttoallradiation.Photonpacketsenteringop-
tically thick cells were immediately absorbed and were not
Figure 2.Benchmarksforanisotropicscatteringmedium(α= remitted, allowing straightforward comparison with analyt-
1).SymbolsarethesameasthoseforFigure1. ical results.
Inthiscasetheradiationpressuresweepsupathinshell
whose equation of motion (assuming a negligible contribu-
tion from thermal gas pressure) is given by
and as expected the noisiest estimator is the momentum
tracking algorithm for the τ =0.1 model. d (cid:104)(cid:16)4πR3ρ (cid:17)R˙(cid:105)= L. (17)
Finallyweranthepurescatteringmodelforaforward- dt 3 0 c
scattering Henyey-Greenstein phase function (g = 0.9). In Integrating the above expression we get
this model the nett radiation force should be reduced by a
(cid:18) (cid:19)1/4
factor of (1−g) over the corresponding isotropic scattering 3L
R= t1/2. (18)
case, i.e. 2πcρ
0
(cid:13)c 2015RAS,MNRAS000,1–11
Massive star formation simulation algorithms 5
Figure 4. The development of the radiation-pressure driven
sphericalshellasafunctionoftime.Themassdensityasafunc- Figure 5. The evolution of the shell radius with time for the
tion of radius is plotted at 1000kyr intervals from a simulation radiation pressure hydrodynamic test (solid line). The theoreti-
timeof1000kyrto8000kyr.Somebroadeningoftheshellisevi- calshellexpansionbasedonthethin-shellapproximationisalso
dentatlatertimes. shown(dashedline).
Wefollowedthedevelopmentoftheshellfor8kyr,writ-
ingouttheradialdensityprofileat109sintervals(see Fig- ofpacketsthatemergefromthegridisthenN =NhighNlow
Thekeyistocorrectlyidentifythepointatwhichthepacket
ure 4), and we then determined the position of the shell
splitting takes place, in order that (i) high energy packets
foreachradialprofile.Thiswasdonebymakingaparabolic
do not propagate into optically thin regions, since they will
fit to the peak density in the grid and the density at the
increasethenoiseintheMCestimators,(ii)lowenergypack-
adjacent radial grid points. We show the results of the dy-
ets undergo a small (but non-zero) number of interactions
namical calculation in Figure 5. There is good agreement
before leaving the computational domain.
between the model calculation and the theoretical growth
We constructed a one-dimensional test model for the
of the shell predicted by equation 18. Deviations from the-
packetsplittingalgorithm,comprisinga0.1 pcradiussphere
ory at late times may be attributed to departures from the
containing 100M of gas and an r−2 density profile. The
thin shell approximation. (cid:12)
dust opacity was from Draine & Lee (1984) silicates with a
dust-to-gas density ratio of 1 per cent. The dust is heated
byacentralluminousprotostarofT =4000Kandradius
3 PACKET SPLITTING eff
R = 150R . We first calculated the radiative equilibrium
(cid:12)
The splitting of particles in the MC method has a withoutpacketsplittingusingN =105 photonpackets,and
long history, dating back to the first neutron trans- used the final iteration of the radiative equilibrium calcula-
portcodes(Cashwell&Everett1959).Splitting,and tionasourbenchmark.Wesubsequentlyranthesamemodel
its inverse, the so-called Russian Roulette method, with packet splitting and N = 1000 and N = 100,
high low
are variance reduction techniques designed to im- defining the high optical depth region as that for which
prove the efficiency with which the MC sampling the Planck-mean optical depth to the sphere radius was
operates.AlthoughtheoriginalLucy(1999)algorithmhad 20. Note that in two or three-dimensional problems
equal energy photon packets, this is not a fundamental re- the integral of the Planck-mean optical depth is cal-
striction.Inthecaseofstarformationcalculationsthepho- culated in the positive and negative directions of
ton packets may be produced within gas with very high each coordinate axis, defining a typically ellipsoidal
optical depth, and the packets may undergo thousands of boundary for packet splitting.
scattering events prior to emerging from the computational We plot the temperature profile of the sphere in Fig-
domain. The MC estimators of the energy density and ra- ure7forbothcases,andfindgoodagreement.Inparticular
diation pressure will be good quality for the optically thick thereisasmoothchangeinthetemperaturethroughtheop-
cells where packets spend most of their time, and thus it ticaldepthboundary(theregionwherethepacketsplitting
is inefficient to use equal energy packets. Instead we em- occurs).
ploy a packet splitting algorithm, in which a lower number Inthecalculationwithoutsplittingeachphotonpacket
of higher energy packets (N ) are released from the pro- undergoes ∼20000 absorption or scattering events prior to
high
tostar. These propagate through the optically thick region itsemergencefromthecomputationaldomain.Whenpacket
andthenareeachsplitintoanumberoflowerenergypackets splittingisenabledeachhighenergyphotonpacketstillun-
(N )intheopticallythinregion(see6).Thetotalnumber dergoes ∼ 20000 absorptions or scatterings, but the low
low
(cid:13)c 2015RAS,MNRAS000,1–11
6 Tim J. Harries
energy photon packets only undergo ∼ 200. This results in
a speed-up of the radiation step of a factor of ∼50.
4 PARALLELISATION
The primary drawback in such a detailed treatment of the
radiation field is the computational effort involved. Fortu-
natelytheMCmethod,inwhichthephotonpacketsarees-
sentially independent events, may be straightforwardly and
effectively parallelised. The top level of parallelisation in-
volvesdomaindecomposingtheoctreethatstorestheAMR
mesh. Each sub-domain belongs to a separate MPI thread,
enabling the code to run on distributed memory machines
andenablingtheuseofdomainswithalargermemoryfoot-
printthanthatavailableonasinglenode.Wechooseasim-
ple domain decomposition, which although not necessarily
loadbalancing,doesallowastraightforwardimplementation
inthecode.Foreight-waydecompositioneachbranchofthe
octree is stored on an individual MPI process (with an ad-
Figure 6.Acartoonillustratingthepacketsplittingmethod.A ditional thread that performs tasks such as passing photon
highenergypacketisreleasedfromtheprotostarandundergoes packetstothreads).Similarly64-waydecompositionmaybe
manyabsorptionandscatteringevents(blackarrows).Eventually achievedbydomaindecomposingfurtherdownthebranches
thepacketpassesbeyondtheregionidentifiedasbeingoptically oftheoctree,andthisisthedecompositionthatisemployed
thick(dashedline)andatthispointthepacketissplitintoNlow for majority of our runs (although 512-way decomposition
low-energy packets (red arrows) which eventually emerge from
is implemented we do not have access to the resources nec-
thecomputationaldomain(hereN =5).
low essary to regularly run the code in this mode).
The main bottleneck is the communication overhead
whenpassingphotonpacketsbetweenthreads,andweopti-
mize this by passing stacks of photon packet data between
the MPI threads rather than individual packets in order to
reduce latency. Thus when a photon packet reaches
a boundary between domains (naturally this always
corresponds to cell boundary) then the packet posi-
tion, direction, frequency, and energy of the photon
packetisstoredonastack.Oncethestackreachesa
setnumberofpackets(typically200),thenthestack
is passed to the appropriate MPI thread. We note
thatthisalgorithmcloselyresemblestheMILAGRO
algorithm described by Brunner et al. (2006). Fur-
theroptimisationofthestacksizeispossible,includ-
ing varying the stack size across the domain bound-
aries and dynamically altering the stack size as the
calculationprogresses,butwehaveyettoimplement
this.
A further level of parallelisation is achieved by having
many versions of identical computational domains (which
we refer to as hydrodynamic sets), over which the photon
packetloopissplit,withtheresultsderivedfromtheradia-
tion calculation (radiation momenta, cell integrals etc) col-
Figure 7. The results of the packet splitting test. The temper-
latedandreturnedtoallthesetsattheendofeachiteration.
ature of the sphere is plotted as a function of radial distance
Thethermalbalanceandionizationequilibriumcalculations
(in units of the sphere radius Rs = 0.1pc). The model without
(whichcanbetimeconsuming)arethenparallelisedoverthe
packet splitting (+ symbol) shows excellent agreement with the
temperaturefoundusingpacketsplitting(×symbol).Theverti- sets,witheachthreadcorrespondingtoaparticulardomain
callineindicatestheradiusbeyondwhichthehighenergypackets performing the equilibrium over an appropriate fraction of
undergosplitting. its cells before the results are collated and distributed (see
Figure 8).
We performed a scaling test to demonstrate the effi-
cacy of our scheme. We ran the model detailed in section 3
in 2D. The benchmark used one hydrodynamics set with
the quadtree domain-decomposed over 4 threads, giving a
total of 5 MPI threads including the control thread. We
(cid:13)c 2015RAS,MNRAS000,1–11
Massive star formation simulation algorithms 7
Set 1 Set 2
on Domain 1 Domain 2 Domain 1 Domain 2
ati
er
n it
o
ati
di
a
C r Domain 3 Domain 4 Domain 3 Domain 4
M
m
u
bri
uili Set 1 Domain 1 Domain 2 Domain 3 Domain 4
q
e
n
o
ati
niz
o
oi
ot
ph Set 2 Domain 1 Domain 2 Domain 3 Domain 4
al/
m
er
h
T
Figure 9. Results of the scaling test of the parallelisation. The
Figure 8. A schematic showing the parallelisation of the radia-
tionloopintorus.Herethephotonpacketsaresplitacrosstwo same2-Dmodelwasrunwith1,3,6and12hydrodynamicssets,
corresponding to 5, 15, 30 and 60 MPI threads. The wall time
setsoftheAMRmesh(greensquares)andeachmeshisdivided
for one radiation timestep is shown (normalized by the 5-thread
intofourdomains(bluesquares),withinter-domaincommunica-
calculation).Thesolidlineindicatesperfectscaling.
tionofphotonpackets(indicatedbytheblackarrows).Attheend
oftheiterationtheMonte-Carloestimatorsarecombinedforeach
domaininordertocomputerevisedtemperaturesandionization
fractionsofthegas.
recorded the wall time for one radiative transfer step, and
subsequently ran this same model for 3, 6, and 12 hydro- tations described deal with items (i) and (iii) in the same
dynamicssets(correspondingto15,30and60MPIthreads manner,butdifferinthewaythataccretionistreated.The
respectively). We found excellent scaling (see Figure 9) up Federrath et al. method defines an accretion radius (as a
to 60 threads, with the slight deviation from embarrassing small multiple of the smallest cell size in their AMR mesh)
scalability due to the all-to-all MPI communication neces- andsimplyteststhatgasincellswithintheaccretionradius
sary when the results of the MC estimators are gathered isboundtothesinkparticle,andthenaccretesthegasmass
(this is approximately 10 per cent of the wall time for the (and its associated linear and angular momentum) above a
60 thread run) thresholddensityontothesinkparticle.TheKrumholzetal.
method uses a dynamically varying accretion radius, which
ranges from 4-cells to the Bondi-Hoyle radius. The method
thenascribesanaccretionrateontothesinkbasedonBondi-
5 SINK PARTICLE IMPLEMENTATION
Hoyle-Lyttletonaccretionandremovestheappropriatemass
We implemented a Lagrangian sink particle scheme within fromeachcellwithintheaccretionradiusbasedonaweight-
torus in order to follow the formation of stars in our hy- ing function that falls rapidly with distance from the sink
drodynamicscode.Originallyadoptedforsmoothed-particle (cellsoutsidetheaccretionradiushaveaweightingofzero).
hydrodynamics solvers (Bate et al. 1995), sink particles al- Theadvantage ofthismethodis thatanappropriateaccre-
low one to overcome the restrictively small Courant time tionrateismaintainedinsituationswheretheBondi-Hoyle
associated with high densities as the collapse proceeds, by radius is unresolved, whereas the Federrath method breaks
removinggasfromthecomputationandaddingittopoint- down in this regime. However the Krumholz method relies
like particles that interact with the gas via gravity (and on a statistical smoothing of the flow in order to determine
possibly radiation), but not via thermal gas pressure. the appropriate far-field density and velocity to use when
The fundamental principles for creating a reasonable calculating the Bondi-Hoyle accretion rate. In contrast the
sink particle implementation are (i) creating sink particles Federrath method relies primarily on the accretion flow be-
onlyasa‘lastresort’,(ii)accretinggasinamannerthathas ing supersonic outside the accretion radius, and therefore
minimal impact on the dynamics of the gas immediately the algorithm by which material is removed from the grid
around the sink particle, (iii) properly accounting for the has no impact on the upstream dynamics.
gravitational forces between Lagrangian sink particles and In our implementation we use a method akin to Feder-
the gas on the grid. rath’s, since we are always in the regime where the Bondi-
The use of sink particles in grid-based hydrodynamics Hoyle radius is well resolved. Below we describe the details
codes has been investigated by Krumholz et al. (2004) and of our sink particle implementation, and several accretion
Federrathetal.(2010).Broadlyspeakingthetwoimplemen- and dynamics tests of the method.
(cid:13)c 2015RAS,MNRAS000,1–11
8 Tim J. Harries
5.1 Sink particle creation andthesinkparticlerespectively,histhegravitationalsoft-
eninglength,ands(r ,h)isacubicspinesofteningfunction
WeadoptthesamecriteriaforsinkparticlecreationasFed- ij
(see equation A1 of Price & Monaghan 2007).
errath et al. (2010). For the cell under consideration we de-
Forthecellcontainingthesinkparticleweinsteadper-
fine a control volume that contains all cells within a prede-
form a sub-grid calculation by splitting the mass of the cell
fined radius r . Before a sink particle is created a number
acc into 83 subcells and sum the force over the subcells in an
ofchecksonthehydrodynamicalstateofthegasinthecon-
analogous manner equation 20.
trol volume must be passed, briefly:
Thesink-sinkinteractioniscomputedusingasumofthe
• The central cell of the control volume must have the gravitational forces over all sinks (this is computationally
highest level of AMR refinement. tractable since the number of sink particles is small). The
• Thedensityofthecentralcellmustexceedapredefined equation of motion of the sink particles is integrated over a
thresholddensityρ ,thusensuringthat∇·(ρv)<0for timestep by using the Bulirsch-Stoer method (Press et al.
thresh
that cell. 1993), once again using the cubic spline softening of Price
• Flowsincellsalongtheprincipleaxesmustbedirected & Monaghan (2007).
towards the central cell.
• Thegravitationpotentialofthecentralcellmustbethe
5.4 Gas motion
minimum of all the cells in the control volume.
• The control volume must be Jeans unstable i.e. The self-gravity of the gas is solved using a multigrid solu-
|Egrav|>2Eth. tiontoPoisson’sequation,andwiththegradientofthepo-
• The gas must be in a bound state i.e. Egrav +Eth + tential appearing as a source term in hydrodynamics equa-
Ekin <0, where Ekin is the kinetic energy of the gas in the tions.
control volume where the speeds are measured relative to
the velocity of the centre of mass of the control volume.
• Thecontrolvolumemustnotoverlapwiththeaccretion
6 SINK PARTICLE DYNAMICS TESTS
radius of any pre-existing sink particles.
Following Hubber et al. (2011) we adopt the three-body
If these tests are passed then a sink particle is created at
Pythagorean test problem of Burrau (1913) to benchmark
thecentreofthecontrolvolume,andaccretesgasaccording
thesink-sinkgravitationalinteractionsandtesttheintegra-
to the method detailed below.
tor. This problem has masses (taking G = 1) of m = 3,
1
m = 4 and m = 5 starting at rest at Cartesian coordi-
2 3
nates (1,3), (−2,−1) and (1,−1) respectively. The subse-
5.2 Accretion onto sink particle
quentmotionofthemassesinvolvesanumberofcloseinter-
Each sink particle has an associated accretion radius r actions between sinks, and was first numerically integrated
acc
which is defined as 2.5 times the size of the smallest cell by Szebehely & Peters (1967). The results of the test prob-
in the AMR mesh. The accretion radius defines a spherical lem are given in Figure 10 and show excellent agreement
volume, and after each hydrodynamics step each the cells with Figures 2, 3, and 4 of Szebehely & Peters (1967). The
in the accretion volume are checked against a predefined totalenergyofthesystemoverdurationofthecomputation
densitythreshold(ρ ).Cellswithadensityaboveρ is conserved to better than 1 part in 106.
thresh thresh
undergo a further check to see if the gas in them is bound The implementation of the gas-on-sink gravitational
to the sink particle, that is force calculation was tested by using a power-law (ρ ∝
ρ(r)−2) density sphere of 100M and radius 0.1pc, on a
(cid:12)
E +E +E <0. (19)
grav kin th 3D AMR mesh with minimum cell depth 6 (equivalent to a
Note the kinetic energy of the gas is measured relative to fixed-gridresolutionof64×64×64)andmaximumcelldepth
the velocity of the sink particle. If the gas is both bound 10(1024×1024×1024).Thecomputationaldomainwasdi-
and above the density threshold then the mass of the sink videdover64MPIthreads.Threetestsinkparticlesofneg-
particle is increased by (ρ−ρ )V where ρ and V are ligiblemasswereplacedinthegridatradiiof2.5×1017cm,
thresh
the density and volume of the cell. The linear and angular 1×1017cm, and 5×1016cm, at the appropriate Keplerian
momentumoftheaccretedmassisaddedtothesinkparticle orbital speed (∼ 2.074kms−1). The hydrodynamics of the
and subtracted from the cell. gas was neglected for this test, and the sink particles were
allowed to move through the stationary gas under gravi-
tational forces only. The benchmark test was run for 130
5.3 Sink particle motion orbits of the inner-most particle (corresponding to ∼26 or-
bits of the outermost particle). The orbits of the particles
Thegravitationalinfluenceofthegasonagivensinkparticle are overlaid on the density distribution and adaptive mesh
isfoundbysummingupthegravitationalforcesfromallthe inFigure11,indicatingthattheintegrator,domaindecom-
cellsintheAMRgrid.Forallcellsexceptingthatcontaining position,andgas-particleforcesareoperatingsatisfactorily.
the sink particle we use
(cid:88)
F = GM ρ V s(r ,h)ˆr (20)
j j i i ij
7 ACCRETION TESTS
i
whereM isthemassofthesinkparticle,andr andˆr The principle tests of our sink particle implementation are
j ij
arethedistanceanddirectionvectorbetweenthecellcentre those that ensure that the accretion of gas occurs at the
(cid:13)c 2015RAS,MNRAS000,1–11
Massive star formation simulation algorithms 9
Figure 10. The paths of three masses in the Pythagorean three-body benchmark test. The panels shows show the paths of the sink
particlesatt=0−10(leftpanel),t=10−20(middlepanel),andt=20−30(rightpanel).Mass1isshownasagreen,dashedline,
mass2asared,solidline,andmass3asablue,dottedline.Solidcirclesindicatethepositionoftheparticlesatunittimeintervals.
G2M2
M˙ =4πλ ρ (21)
c3 ∞
s
wherec isthesoundspeedandλ=1.12forisothermalgas.
s
The corresponding Bondi radius is given by
2GM
R = . (22)
B c2
s
For this test we adopted M = 1 M and ρ =
(cid:12) ∞
10−25gcc−1,givinganexpectedaccretionrateof5.9×10−11
M yr−1andaBondiradiusof3.75×1017cm.Weranthree
(cid:12)
two-dimensionalmodelswithresolutionsofR /∆xof5,10,
B
and 20. Rapid convergence in the accretion rate with reso-
lution is observed (see Figure 12). The highest resolution
model has an accretion rate of 5.7×10−11 M yr−1, i.e.
(cid:12)
within 2% of the theoretical rate.
7.2 Collapse of a singular isothermal sphere
logdensity Shu (1977) showed that an isothermal sphere with ρ(r) ∝
1/r2 collapses in such a way that there is a constant mass
fluxthroughsphericalshells.Theaddedlevelofcomplexity
over the Bondi test described above is that the self-gravity
of the gas is significant.
Figure11.Theresultsofthegas-on-sinkgravitationalforcetest.
Thepathsofthethreetestsinkparticlesareshownasboldsolid We adopt the same test as Federrath et al. (2010),
lines,whilstthedensityofthegasisshownasalogarithmiccolour with a sphere of radius R = 5×1016cm with ρ = 3.82×
scale.TheAMRmeshisshownbythinsolidlines. 10−18gcm−3withcontains3.02M(cid:12)ofgas.Thesoundspeed
ofthegaswastakentobe0.166kms−1.Theexpectedaccre-
tionrateofthismodelis1.5×10−4M yr−1.Weadoptedan
(cid:12)
correctrate.Wethereforeconductedthreebenchmarktests
adaptive,two-dimensionalcylindricalmeshforfourlevelsof
of increasing complexity.
refinement, with the smallest cells of R/∆x=300
The model demonstrated excellent agreement with the
Shu prediction, starting with an accretion rate of ∼ 1.5×
7.1 Bondi accretion
10−4M yr−1, and only declining when 90% of the origi-
(cid:12)
ApointmassM accretessphericallyfromagascloud,where nal mass had been accreted (see Figure 13). At the end of
far from M the cloud is stationary and has density ρ . the run the local density approached the global floor den-
∞
Bondi (1952) showed that the maximum accretion rate is sityforthesimulation(10−21gcc−1)andtheaccretionrate
given by approaches the Bondi rate for that density.
(cid:13)c 2015RAS,MNRAS000,1–11
10 Tim J. Harries
Figure 12. Results of the Bondi accretion test. The accretion
rateisplottedagainsttheBonditimescale(tB =RB/cs)forthe
RB/∆x = 5, 10, and 20 models (dotted, dashed and solid lines
repsectively). The thick solid line corresponds to the expected
theoreticalrateof5.9×10−11 M(cid:12)yr−1.
7.3 Bondi-Hoyle accretion
Weconstructeda2Dtestcaseinwhicha1 M sinkisplaced Figure 13. Results of the Shu accretion test. The lower panel
(cid:12)
at the origin in initially uniform density gas (10−25gcc−1) showstheaccretionrateasafunctionoftime(crosses),withthe
expected theoretical accretion rates predicted by the Shu model
with molecular weight of 2.33 and a temperature of 10K.
(solidline)andBondiaccretionatthefloordensity(dashedline).
The gas initially had a constant velocity with a mach num-
The upper panel shows the growth of the central object as a
berofM=3paralleltothez-axis,andaninflowcondition
functionoftime(crosses)alongwiththeexpectedtrendassuming
attheupstreamboundaryandandoutflowconditionatthe
aconstantShuaccretionrate(solidline).
downstreamboundary.Thesearethesameinitialconditions
as the Bondi-Hoyle test case in Krumholz et al. (2004), al-
though they ran their simulation in 3D.
The model extent was 0.78pc, and three levels of re- that the new method works in both the pure absorption
finement wereused withthesmallest cellscorresponding to andpurescatteringregimes,andproperlytreatsanisotropic
2.3 × 1014cm. We ran the model until an approximately scattering processes. The method comprises a simple addi-
steady-state of the accretion rate was achieved (Figure 14). tion to the MC estimators required for radiative and pho-
The theoretical Bondi-Hoyle accretion rate for these con- toionisationequilibriumcalculations,anddoesnottherefore
ditions is 1.7 × 10−12M yr−1, but the Krumholz et al. represent a substantial computational overhead to the MC
(cid:12)
(2004) model, and those of Ruffert (1996), found accretion RHD methods described by Haworth & Harries (2012).
rates of close to 2×10−12M yr−1, with considerable tem- However the MC method as a whole is significantly
(cid:12)
poral variation in the accretion rate on the Bondi-Hoyle morecomputationallydemandingthantheFLDandhybrid
timescale. Our simulation reached a steady-state accretion methods,andwehavethereforedevelopedtwonewmethods
rate of 2.4×10−12M yr−1(Figure 15), which is compara- to ameliorate this. The first is the packet splitting method
(cid:12)
ble to the peak accretion rate seen in the Krumholz et al. detailedinsection3,inwhichweextendLucy’soriginalal-
(2004) simulations. The level of variability is substantially gorithmtoincorporatephotonpacketswithvaryingenergy.
lowerinoursimulation,presumablybecausetheinstabilities ThesecondistodistributetheMCphotonpacketloopover
thatbuildupandmodifythedynamicsnearthesinkinthe many instances of the computational domain, allowing an
3D simulations are absent in our 2D models. excellent scalability (see section 4).
The final element needed to compute massive star for-
mation models is a description of the protostar itself, and
to do this we have included a sink particle algorithm. By
8 CONCLUSIONS
interpolating on protostellar evolutionary model grids (e.g.
We have presented a new method for including radiation Hosokawa & Omukai 2009) as a function of mass and ac-
pressureinRHDsimulationsthatincorporatesalevelofmi- cretionratewillwillbeabletodeterminetemperaturesand
crophysical detail that is significantly greater than that of luminosities for our protostar, and assign it a spectral en-
flux-limited diffusion or hybrid techniques. We have shown ergydistributionfromappropriateatmospheremodels.The
(cid:13)c 2015RAS,MNRAS000,1–11