Mocassin: A fully three-dimensional Monte Carlo photoionization code

18 Pages·1.4 MB·English

by  B. Ercolano

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed2February2008 (MNLATEXstylefilev2.2)
Mocassin: A fully three-dimensional Monte Carlo photoionization
code
B. Ercolano1, M. J. Barlow1, P. J. Storey1, X.-W. Liu1,2
1UniversityCollegeLondon,GowerStreet,LondonWC1E6BT,UK
2Currentaddress:DepartmentofAstronomy,PekingUniversity,Beijing100871,P.R.China
3
0
0
2 Received:

n
a
ABSTRACT
J

7
The study of photoionized environments is fundamental to many astrophysical prob-

lems. Up to the present most photoionization codes have numerically solved the equations
3
v of radiative transfer by making the extreme simplifying assumption of spherical symmetry.
8 Unfortunatelyveryfew realastronomicalnebulaesatisfy this requirement.To remedythese
7 shortcomings,aself-consistent,three-dimensionalradiativetransfercodehasbeendeveloped
3 usingMonteCarlo techniques.The code,Mocassin,is designedto build realistic modelsof
9 photoionizednebulaehavingarbitrarygeometryanddensitydistributions,withboththestel-
0 lar and diffuse radiation fields treated self-consistently. In addition, the code is capable of
2 treatingonesormoreexcitingstarslocatedatnon-centrallocations.
0
ThegaseousregionisapproximatedbyacuboidalCartesiangridcomposedofnumerous
/
h cells. The physical conditions within each grid cell are determined by solving the thermal
p equilibriumandionizationbalanceequations.Thisrequiresaknowledgeofthelocalprimary
- and secondary radiation fields, which are calculated self-consistently by locally simulating
o
theindividualprocessesofionizationandrecombination.Thestructureandthecomputational
r
t methodsusedintheMocassincodearedescribedinthispaper.
s
a Mocassin has been benchmarkedagainst established one-dimensionalspherically sym-
: metriccodesforanumberofstandardcases,asdefinedbytheLexington/Meudonphotoion-
v
ization workshops(Pe´quignot1986; Ferlandetal. 1995; Pe´quignotetal. 2001). The results
i
X obtainedfor the benchmarkcases are satisfactory and are presentedin this paper. A perfor-
manceanalysishasalsobeencarriedoutandisdiscussedhere.
r
a
Keywords: HIIregions–planetarynebulae:general–ISM:abundances–atomicprocesses
1 INTRODUCTION more ions in calculations. Mendoza (1983) presented a compila-
tion of radiative and collisional data for collisionally excited ul-
Amongstthefirstnumericalmodelsforphotoionizedgaseousneb-
traviolet, optical and infrared lines which was widely adopted,
ulae were those calculated by Flower (1968), Harrington (1968)
withsomeofthesedatastillinusetoday,thoughmosthavebeen
and Rubin (1968). These early models included the basic phys-
sepercededbymorerecentcalculationssuchastheR-matrixcalcu-
ical processes of ionization and recombination of hydrogen and
lations of the Iron Project (Hummeretal. 1993) and the Belfast
helium, thermal balance and escape of the emitted photons from
group (e.g. McLaughlin&Bell 1998; Ramsbottometal. 1998).
the nebula. However, the lack of reliable atomic data heavily
Currently, radiative and dielectronic recombination rates are still
limited the success of these models, as well as the fact that
highly uncertain or unavailable for some ions; recent efforts to
a number of important physical processes, such as charge ex-
improve the situation have been reviewed by Nahar&Pradhan
change and dielectronic recombination (Aldrovandi&Pe´quignot
(1999) and Nahar (2000). Most photoionization models include
1973;Pe´quignotetal.1978;Storey1981),werenotaccountedfor
temperature-dependentanalyticalfitstotheserecombinationrates,
atthetime. Theevolution of photoionization modelling hasgone
such as those of Aldrovandi&Pe´quignot (1973) for radiative
handinhandwithadvancesmadeinatomicphysicsandcomputer
and high temperature dielectronic recombination, and those of
technology. Theapplicationofphotoionization modelstoawider
Nussbaumer&Storey(1983)forlowtemperaturedielectronicre-
rangeofionshasbeenaidedbythephotoionizationcross-section
combination.
calculationsbyReilman&Manson(1979),and,morerecently,the
OpacityProject(Hummeretal.1993).Compilationsbasedonthe Available computer power has increased enormously since
latter’s data (e.g. Verner&Yakovlev 1995), have made possible the dawn of photoionization modelling. This has allowed more
theinclusion of accurate photoionization cross-sections for many complexmodelstobebuilt,includingmoreions,morefrequency
2 Ercolanoet al.
points, more linesand moreatomiclevels. Nevertheless, thefun- (Flanneryetal.1980;Boisse´1990),resonance-likescatteringinac-
damental assumption of spherical symmetry has always been re- cretiondiscwinds(Kniggeetal.1995)andpolarizationmapsfor
tained.However,aglanceatanimageofanyGalacticHIIregion the circumstellar envelopes of protostars (Fischeretal. 1994). In
willimmediatelydemonstratethattheseobjectsareneitherspher- theexamples described above theabsorption andscattering coef-
icallysymmetricnorhomogeneous.Inaddition,theyusuallyhave ficientsarenotcoupled totheradiationfieldand, therefore,these
multipleexcitingstarslocatedatnon-centralpositionsintheneb- problemsdonotrequiresolutionbyiteration.
ula.Bycontrast,planetarynebulae(PNe)haveonlyasingle,cen- However, Monte Carlo techniques have also been used for
trallylocated,excitingstar.However,evenforPNe,sphericalsym- dust radiative equilibrium calculations for some time, see e.g.
metry is not a realistic assumption, as demonstrated by observa- Lefevreetal. (1982), Lefevreetal. (1983) and, more recently,
tionswithinstrumentssuchastheHubbleSpaceTelescope,which Wolfetal.(1999).Theseauthorsuseatechniqueinwhichstellar
revealanoverwhelmingvarietyintheshapesofplanetarynebulae. and diffuse photon packets are emittedseparately; thenumber of
Theseobjectsareveryrarelycircularinprojection;arecentstudy diffuse photon packets (i.e.packets emittedby thedust) isdeter-
inferredthatabout50%ofallknownplanetarynebulaearelowec- minedbythedust graintemperature, whichinturnisdetermined
centricityellipticals,whileonlyabout10%arecircularinprojec- bythebalancebetweenthenumberofabsorbedandemittedphoton
tion,withtheremainderhavingmoreextremeellipticalorbipolar packets.Aninitialguessforthedustgraintemperatureisprovided
geometries(Soker1997,2001).Someobjects,forexamplethetwo bythenumberofpacketsabsorbed,andtheiterationcontinuesun-
youngplanetarynebulaeHe2-47andPNM1-37,(alsodubbedthe tilthe grain temperatures converge. Using this method the stellar
starfishtwins;Sahai2000),show evenmorecomplicatedgeome- luminosityisnotautomaticallyconservedduringtheMonteCarlo
tries,withmultiplelobes.OtherPNehaveFLIERs(fast,lowion- simulation;onlyafterthegraintemperatureshavereachedconver-
ization emitting regions; Balicketal. 1993, 1994, 1998), BRETS genceisthestellarluminosityapproximatelyconserved.Thecon-
(bipolar,rotating,episodicjets;e.g.Lopezetal.1993),ansae,jets, vergenceofsuchcodesisoftenveryslowandrequiresalargenum-
knots, filaments, tails or multiple envelopes. (see e.g. Perinotto berofiterationsandsimulationquantainordertoreachtherequired
2000;Corradietal.1999;Garc´ia-Segura1997). accuracy.
To our knowledge, only two three-dimensional photoioniza- Bjorkman&Wood(2001)havedescribedageneralradiative
tioncodeshavebeendevelpedsofar,onebyBaesgenetal.(1990) equilibriumandtemperaturecorrectionprocedureforuseinMonte
and the other by Gruenwald, Viegas&de Broguiere (1997). The Carlo radiative transfer codes having sources of temperature-
firstcodeusedafixednumberofequallysizedcellsandtheon-the- independent opacity, such as dust. Their technique makes use of
spotapproximationforthediffuseradiationfield,withonlythesix informationnaturallygivenbytheMonteCarlomethod,which,by
more abundant chemical elements being taken into account. The tracking every photon/energy packet, makes it easy to determine
work by Gruenwaldetal. (1997) improves on this by allowing a whereinthesimulationgridenergyisbeingabsorbed. Whenen-
more flexible spatial grid and by using an iterative technique for ergyisdepositedatagivenlocation,followingapacket’sabsorp-
thedeterminationofthediffusefieldandalsobyincludingtwelve tion, the local medium is heated. Whenever this occurs the new
chemicalelementsinthesimulations. local temperataure iscalculated and thepacket isthen re-emitted
Since most existing one-dimensional photoionization codes accordingly.Thepacketsarefollowedintheirpaththroughthere-
are based on the numerical solution of the equations of radiative gion,astheyundergo scatteringsandabsorptionsfollowedbyre-
transferassumingsphericalsymmetry,theirexpansiontothreedi- emissions, with the temperatures being updated after each event,
mensions can be either very difficult or impractical, resulting in untilthepacketsreachtheedgeofthenebulaandescapetoinfinity,
verylargecodes.TheMonteCarloapproachtotransferproblems hence contributing tothe emergent spectrum. Once all the stellar
providesageometry-independent techniquewhichcanhandlethe photon packets have escaped, the resulting envelope temperature
radiation transport problem self-consistently. With this in mind, andtheemergentspectrumarecorrectwithouttheneedofanyfur-
the Mocassin code (MOnte CArlo SimulationS of Ionised Nebu- theriterations.
lae)wasdeveloped, inordertoprovideathree-dimensional mod- A great limitation of Bjorkman & Wood’s method is that it
ellingtool capable of dealingwithasymmetric and/or inhomoge- cannotbeappliedtosituationswheretheopacitiesaretemperature-
neous nebulae, as well as, if required, multiple, non-centrally lo- dependent, as is the case in photoionized nebulae. There are two
catedexcitingstars. reasonsforthefailureofthismethodwhentheopacityvarieswith
Section 2contains a description of thegeneral Mocassin ar- the local temperature: firstly, the number of photon packets ab-
chitecture and of some of the main computational methods used sorbedbythecellpriortoatemperatureupdatewouldbeeithertoo
in the code. The code has been benchmarked against established small or too large, and, secondly, achange in temperature would
spherically symmetricone-dimensional photoionization codes for alsoimplyachangeoftheinteractionlocationsofpreviouspack-
a set of standard nebulae and in Section 3 we present the results ets,signifyingthatthepathsofthepreviousphotonpacketsshould
ofthisbenchmarking,togetherwithaperformanceanalysisofthe havebeendifferent.While,itisclearthat,whendealingwithpho-
codes.Insection4wediscusstheresultsofthebenchmarkingand toionised gas, Bjorkman & Wood’s technique is not applicable,
presentsomegeneralguidelinesonhowtorunthecodeefficiently. their work is nevertheless very enlightening and should be taken
intoaccountforfurtherdevelopmentsoftheMocassincode,when
atreatmentfordustgrainswillbeintroduced.
ArecentexampleoftheapplicationoftheMonteCarlotech-
2 DESCRIPTIONOFTHEMONTECARLOCODE nique to problems requiring solution by iteration is the work of
Lucy(1999),whoobtainedthetemperaturestratificationandemer-
2.1 Background
gentspectrumofanon-greysphericallysymmetricextendedstel-
TheMonte Carlomethod hasbeenwidelyappliedtoavarietyof laratmosphereinLTE.Hisresultsshowverygoodagreementwith
astrophysical problems, such as the penetration of ultraviolet ra- thepredictionsofCastor(1974),hencedemonstratingthevalidity
diation into the interiors of uniform or lumpy, interstellar clouds of the Monte Carlotechniques applied, someof which werealso
Mocassin 3
usedinthedevelopmentofMocassin.Thecurrentworkfolowsthe 2.3 Initiation
approach described by Lucy (1999) and also applied in the one-
In our modelling the gaseous region is approximated by a three-
dimensional code developed by Och, Lucy & Rosa (1998). They
dimensional Cartesian grid, where the ionising source can be
employedadifferentMonteCarlotreatmentoftheradiativetrans-
placedatthecentreofthegridoranywhereelseinthegrid.This
ferinordertoiterativelydeterminethetemperatureandionisation
feature is very useful when dealing with axisymmetric nebulae,
stratification for a spherically symmetric photoionised nebula of
since,byplacingthesourceinacornerofthegrid,weneedonly
uniform density. Some of the techniques that they used are also
consideroneeighthofthenebula,whichcanthenbereconstructed
describedindetailbyLucy(1999,2001,2002).Thebasicconcept
infullattheendofthesimulation.Thisallowstherunningofmod-
isthat, when calculating radiative equilibrium temperatures, con-
elswithmuchhigherspatialresolutionthanthosewhichwouldbe
servation of stellar luminosity is more important than the details
possibleifafullnebulahadtobeconsidered,byputtingthesource
ofthespectralenergydistribution.Withthisinmindconservation
inthecentreand,therefore,notmakinguseofanysymmetryprop-
of stellar luminosity is enforced by using energy packets of con-
ertiesof theobject. Switchesbuilt inside the code allow the user
stantnetenergythroughoutthesimulations.Moreover,allabsorbed
tospecifywhetherthenebulahassomedegreeofsymmetryand,if
packets are re-emitted immediately after every absorption event.
so,whetherthesymmetryistobeused.
The frequencies of the re-emitted energy packets are determined
Insideeachgridcellallnebular properties,suchasthemass
bythelocal gasemissivities. Although thefrequency distribution
densityofthegas,ρ;theelectrontemperatureanddensity,T and
ofthere-emittedpacketswillnotbecorrectuntilthenebulartem- e
N ; and the frequency dependent gas opacity and emissivity, κ
peratureshaveconverged,thismethodnaturallyenforcesradiative e ν
andj ,areconstantbydefinition.Thermalbalanceandionisation
equlibrium at each point in the nebula and so naturally provides ν
equlibriumareimposedineachgridcellinordertoobtainthephys-
conservationofenergy.Thisnotonlyresultsinasimplercodebut
icalconditionsinthelocalgas.
also makes the convergence of the gastemperatures easier (Lucy
The energy packets are created at the position of the ionis-
1999, 2001). Energy packets will be discussed in more detail in
ingsource andtheyall carrythesameenergy ε ,asdiscussed in
section 2.2. 0
the previous section. The frequency, ν, of each individual packet
emittedisderived fromtheinput spectrum oftheionisingsource
accordingtotheprobabilitydensityfunction
2.2 EnergyPackets
F dν F dν
p(ν)= ν = ν (3)
Tsihstesmofalioncparlliyncsiipmleuloaftionugrthtreeaintmdievnidtuoaflaprpohcoetsosieosnoizfeidonniezbautiloancaonnd- RννmmianxFν′dν′ L∗/(4πR∗2)
recombination. Theradiationfieldisthereforeexpressed interms whereFν isthestellarfluxandR∗isthestellarradius.Thisisthen
ofenergypackets,ε(ν),whicharethecalculationquanta.ε(ν)isa theprobabilityofanenergypacketbeingemittedwithafrequency
packetconsistingofnphotonsoffrequencyνsuchthat lyingintheinterval(ν,ν +dν).Theupperandlowerintegration
limits,ν andν ,havetobechosenproperly,dependingon
min max
ε(ν)=nhν (1) the input spectrum, in order to ensure that the bulk of the radia-
tionisincludedinthefrequencyrange.Asthesourceemitsenergy
Inaddition,wetakeallpacketstohaveconstantenergyε .There
0 isotropically,thedirectionoftravelofeveryenergypacketemitted
areseveralreasonsforchoosingtoworkwithmonochromatic,in-
is chosen randomly. This is done by choosing two random num-
divisiblepacketsofradiantenergyinsteadofphotons.Firstofall,
bers,αandβ,intheinterval [0,1],andcalculatingthefollowing
energypacketsaremorecomputationallyeconomicand,also,since
quantities:
theyallhavethesameenergy,thenthosepacketsemittedinthein-
fraredwillcontainalargernumberofphotonswhich,asaconse- w = 2α−1
quence,willnothavetobefollowedindividually(Abbott&Lucy t = 1−w2
1985).Notethatallenergypacketsarefolloweduntiltheyescape p
θ = π(2β−1)
the nebula, including infrared energy packets. This is in order to
allow the introduction of dust particles into the radiative transfer u = tcosθ
treatmentofMocassin,whichisplannedforthenearfuture.Also, v = tsinθ (4)
asthetotalstellarluminosity,L∗,isevenlysplitamongstthestellar
TherandomunitvectorinCartesiancoordinatesisthen(u,v,w)
energy packets, the energy carried by a single packet in the time
(Harries&Howarth1997).
interval∆t,whichrepresentsthedurationoftheMonteCarloex-
periment,isgivenby
L∗ = ε0 (2) 2.4 Trajectories
N ∆t Once a stellar packet is created at the source and launched into
where N is thenumber of energy packets used in the simulation thenebula,itstrajectorymustbecomputedasitundergoesabsorp-
(Ochetal. 1998). Most importantly, the use of constant energy tionsfollowedbyre-emissionsduetobound-freeandfree-freepro-
packetsisanaturalwayofimposingstrictenergyconservationat cesses.Thetrajectoryendswhenthepacketreachestheedgeofthe
anypointinthenebula(Lucy1999).So,whenapacketofradiant nebula,whereitescapestoinfinityandcontributestotheemergent
energyε(ν ) = ε isabsorbed,itisimmediatelyre-emittedwitha spectrum.
a 0
frequencyν ,whichisdeterminedaccordingtoafrequencydistri- Therearetwomethodstotrackthepacketsanddeterminethe
e
butionsetbythegasemissivityofthecurrentcell.Thepacketemit- locationsoftheabsorptionevents.Considerapacketoffrequency
ted,ε(ν ),willthenhavethesameenergyastheabsorbedpacket, ν ,emittedinthedirectionuˆ.Thefirstofthesemethodsconsists
e p
ε(ν ), meaning that only the number, n, of photons contained in ofcalculatingtherunofopticaldepth,τ ,attheenergypackets’
a νp
thepacketischanged. frequency ν , fromthelocationat whichthepacket isemittedto
p
4 Ercolanoet al.
theedgeoftheionisedregionalongthedirectionoftravel,uˆ.The energypackets(i.e.thosepacketsre-emittedimmediatelyafteran
probabilityofabsorptionalongthatpathisthengivenby absorptionevent)needstobedetermined.Sinceabsorptionandre-
p(τνp)=e−τνp (5) etemdisissoiotnroapriecatlwlyoainnddeptheenrdeefonrteevtheneitrs,dtihreecdtiifofnusoefptarcakveetlsiasrecheomseitn-
andthenormalisedcumulativeprobabilityfunctionisgivenby randomlyusingequations4
P(l) = R0τν0∞p(le)−eτ−ντpνdpτdντpνp 2.5 TheMeanIntensity
= 1−R e−τνp(l) (6) The success of a Monte Carlo model often relies on the careful
choiceof appropriateestimators. MonteCarloestimatorsprovide
where τ (l) isthe optical depth to theabsorption event and l is
νp the means to relate the quantities we observe during our Monte
the path length. The position at which the energy packet will be
Carloexperimenttothephysicalquantitieswewanttodetermine.
absorbed will then be determined by choosing a random number
In a photoionization model, a measure of the radiation field is
in the interval [0,1] and comparing it against P(l). In reality, it
needed,namelythemeanintensity,J .
ismoreconvenient tousetheinverseapproach, wheretheoptical ν
IntheworkofOchetal.(1998),theMonteCarloestimatorof
depth from the energy packet source to the event can be derived
J isconstructed by usingthe definitionof the specificintensity,
fromtheinverseofequation5 ν
I ,insphericalcoordinates,(r,θ),asastartingpoint:
ν
τ (l)=−ln(1−U ) (7)
νp R ∆E =I (r,θ)∆A|cosθ|∆ν∆ω∆t (9)
ν
whereU isarandom number intheinterval [0,1]. Onceτ (l)
R νp
has been calculated then the path length can be directly derived where∆Aisthereferencesurfaceelement,θistheanglebetween
(Harries&Howarth1997). thedirectionoflightpropagationandthenormaltothesurface∆A
ThesecondmethodwassuggestedbyLucy(1999)anditcon- and∆ωisthesolidangle.Themeanintensitycanthenbeobtained
sists of testing whether an absorption event occurs, on a cell by fromthisbycalculatingthezeroordermomentofIν,whichgives
cell basis. Inother words, assume that, withineach uniform cell,
therandompathofapacketbetweeneventsisgivenbyequation7, 4πJ (r)= I dω= ∆E Nk 1 1 1 (10)
whichcorrespondstoaphysicaldisplacement,l,givenby ν ZΩ ν ∆t Xi=1 cosθi∆A∆ν
τνp =κνρl (8) bycomparisonwithequation9.ThesumisoverallpacketsNkwith
frequencylyingintheinterval(ν,ν+dν),crossing∆Aatanangle
whereκ andρarethefrequencydependentabsorptioncoefficients
ν
θ.Asdiscussedabove,∆E/∆trepresentstheenergycarriedbya
and the density of the current cell respectively. The method then
singlepacketinthetimeinterval∆t,since∆E = ε ,whichis
consists of checking whether the displacement l is large enough 0
givenbyequation2.Equation10thenprovidesarelationbetween
to carry the packet out of its current cell. If this is the case, the
the Monte Carlo observables (i.e. the number of energy packets
packetismovedalongitsdirectionoftravel,uˆ,uptotheboundary
withfrequencylyingintheinterval(ν,ν +dν),crossing∆Aat
withtheadjacentcell,whereanewvalueforU iscast,givinga
R
angleθ andthemeanintensityoftheradiationfield,whichisthe
newτ ,andanyfurthermovementofthepacketinthisnewcell
νp requiredphysicalquantity.
istobe followed. Alternatively, if the displacement l isnot large
Theuseof Ochet al.’sestimatorsfor J ,however, becomes
enough to carry the energy packet across the next boundary, the ν
problematicinthenon-sphericallysymmetriccase,sincetherefer-
packetwillbeabsorbedandthenre-emittedattheend-pointofthe
encesurfaceforthevolumeelementsinanarbitrarytwo-orthree-
displacement.Lucyalsoclarifiesinhispaperthattheselectionofa
dimensionalcoordinatesystemmightnotbeuniqueorasobvious
newvalueofτ atthecrossingofaboundarydoesnotintroduce
νp as in the one-dimensional case. In our work, a more general ex-
abiassinceaphotonalwayshasanexpectedpathlengthtoitsnext
pressionfortheestimatorofJ issought,and,therefore,following
eventcorrespondingtoτ = 1,regardlessofthedistanceitmight ν
ν
Lucy’sargument (Lucy1999), anestimator for J isconstructed
alreadyhavetravelled. ν
startingfromtheresultthattheenergydensityoftheradiationfield
Inthisworkbothmethodswereimplementedinthecode,in
inthefrequencyinterval(ν,ν +dν)is4πJ dν/c.Atanygiven
turn,inordertotesttheirrespectiveperformances.Thefirstmethod ν
time,apacketcontributesenergyε(ν)=ε tothevolumeelement
provedtobemuchmorecomputationally expensive thenthesec- 0
whichcontainsit.Letlbeapacket’spathlengthbetweensucces-
ond. Thisisdue tothefact that, inorder totrackdown theposi-
siveevents,wherethecrossingofcellboundariesisalsoconsidered
tionat whichanenergy packet isabsorbed, usingour knowledge
anevent;thecontributiontothetimeaveragedenergycontentofa
of τ (l), an array searching routine has to be used to locate the
νp volumeelement,duetothelfragmentsoftrajectory,isε δt/∆t,
indexofτ (l)withinthearrayofopticaldepthscalculatedfrom 0
νp whereδt = l/c.Fromthisargumentitfollowsthattheestimator
thepacket’ssourcetotheedgeofthenebula.Althoughthesearch-
forthevolumeelement’senergydensitycanbewrittenas
ingprocedureemploysabisectiontechnique,whichmakesitquite
efficient,thelargenumberofcallstoit,duetothelargenumberof 4πJ dν ε 1 l
ν = 0 (11)
energy packet interactions withinasimulation, means that nearly c ∆tV c
X
60%oftheruntimeisspentcarryingoutthesesearches.Thesec- dν
ondmethoddoesnotrequireanycallstothearraysearchingrou- whereV isthevolumeofthecurrentgridcellandthesummation
tine, asthepackets arefollowed stepby stepthrough thenebula, isoverallthefragmentsoftrajectory,l,inV,forpacketswithfre-
and this results in the run time being considerably reduced. The quencieslyingintheinterval(ν,ν+dν).Again,arelationbetween
currentversionofMocassinthereforeusesLucy’sappoachtotrack MonteCarloobservables(i.e.theflightsegments,l)andthemean
theenergypacketsthroughoutthenebula. intensity of the radiation field, J has been obtained. Moreover,
ν
Finally, the direction of travel of the newly emitted diffuse equation 11 is completely independent of the coordinate system
Mocassin 5
used and, indeed, of the shapes of the volume elements, V. An- estimatedas a function of temperature using the data of Robbins
otherimportantaspectofthisapproachisthatallpacketspassing (1968).Thecontributionsduetotheselinestothetotalenergydis-
throughagivengridcellcontributetothelocalradiationfieldeven tribution,fromwhichtheprobabilitydensityfunctionsarederived,
withoutbeingabsorbed;thismeansthatequation11returnsestima- areaddedintotherespectiveenergybins.Similarly,He II Lyman
torsoftheradiationfieldevenintheextremelyopticallythincase linescanionizebothneutralhydrogenandneutralhelium,aswell
when all packets pass through the nebula without any absorption assomeofthelowionsofheavierelements.Thereforetheemissiv-
events.Fromthisargumentitfollowsthatthistechniqueallowsa itiesofHeIILymanlineswithupperlevelsfromn=2throughn=5
muchbetter samplingand, hence, ingeneral, muchlessnoisy re- (fits to Storey&Hummer 1995) are also estimated as a function
sultscomparedtoothertechniquesbasedonestimatorsforwhich oftemperatureandtheircontributionstothetotalenergydistribu-
onlypacketsabsorbedwithinagivenvolumeelementcount. tionaddedintotherespectivefrequencybin,asfortheHeIlines.
Thismethodisbasedonthefactthatallemissionprofilesarecur-
rentlytreatedasδ functionsandthelineopacityisassumedtobe
2.6 Gasemissivityandthediffusionofenergypackets zero;andtheabsorptionofenergypacketsisonlyduetothecontin-
uumopacity.Finally,theemissivitiesofthecollisionallinesofthe
Aswehavealreadydiscussedinprevioussections,afteranenergy
heavierionsarecalculated.Thisisdonebyusingmatrixinversion
packetisabsorbed,anewpacketisre-emittedfromthesameloca-
procedures inorder tocalculatethelevel populations of theions.
tioninarandomdirection.Thefrequencyofthere-emittedpacket
Appendix1containsreferencesfortheatomicdatausedforeach
iscalculated by sampling the spectral distribution of the total lo-
calemissivity,jtot.Inordertosatisfythethermalbalanceimplied ion.
ν The energy distribution is derived from the total emissivity,
by the Monte Carlo model, all major emission processes have to
summing over all the contributions in a particular frequency in-
betakenintoaccount, includingthecompletenon-ionizing nebu-
terval.Thenon-ionizinglineemissionistreatedseparately,since,
larcontinuumandlineemission,sincetheyarepartoftheenergy
wheneversuchlinepacketsarecreated,theyescapewithoutfurther
budget. The non-ionizing radiation generated in the nebula isas-
interaction1.
sumedtoescapewithoutfurtherinteractionandconstitutestheob-
Once the line and continuum emissivities have been calcu-
servablespectrumwhichcanthenbecomparedwithobservations.
lated, the probability that the absorption of an ionizing energy
Thefollowingparagraphsareconcernedwiththedescriptionofthe
packet willbefollowedbytheemissionof anon-ionizing packet
individualcontributionstothetotalemissivity.
isgivenby:
ThecontinuumemissionduetoHI,HeI,HeIIandtoheavier
ions is included. The H I continuum can be divided into the Ly- jl + νHjcdν
man continuum, which iscapable of ionizing H, andthe Balmer, Pesc = jl + Pjil Xi+ R0jl ν+ νmaxjcdν (13)
Paschen, etc. continua, whicharenotcapable ofionizingH.The i Xi HeI HeII 0 ν
P P P R
emissivity in the Lyman continuum is calculated directly from a whereν isthehigherlimitofthefrequencygrid;thejl arethe
max Xi
combinationoftheSahaandMilnerelations: emissivitiesofthenon-ionizingrecombinationlinesofallspecies
jν = hcν23ωωi+i1(2πmh2kTe)3/2aν(Xi)e−h(ν−ν0)/kTeXi+1Ne(12) cjHolneIsidanerdedj;HljeIνcIiasrethethefrecqounetrnicbyutdioenpsenddueenttocotnhtoisneuurmecoemmbisinsiavtiitoyn;
linesof He I and He II whicharecapableof ionizing neutral hy-
whereω andω arethegroundstatestatisticalweightoftheions
i i+1 drogenandneutralhelium.Thechoicebetweenthere-emissionof
involved, Xi+1 isthe abundance of the recombining ion, aν(Xi) anionizing photon or anon-ionizing oneismade at thispoint in
isthephotoionization crosssection and ν isthephotoionization
0 thecode.
threshold.Theemissivityoftheotherseriescontinuaareobtained
Ifanionizingenergypacketistobere-emitted,thenthenew
by interpolation of published data (Ferland 1980). A similar ap-
frequencywillbecalculatedaccordingtothenormalisedcumula-
proach is used for the He I and the He II continua, where for tive probability density function for the ionizing radiation, given
frequenciesgreaterthan1.8Rydand4.0Ryd,respectively, equa-
by
tion 12 is used, and the emissivities at lower frequencies are ob-
t(ra1iei9ns7.e0Td)hbfeyocrinottnheteripnHoulueamtIiosenemroiiefssstihavenitddyabotayfphFueebarvlliaysnhedelde(m1b9ye8nB0ts)roifswoarnlst&hoecMaHlacetutIhlIaetsweeds- p(ν)= νRνHνmνHaxjνcj′νcd′νd′ν+′+PjHljeHlIe+I+PjHljeHlIeIII (14)
R P P
usingequation12.Inthehydrogeniccase(i.e.HIandHeII),the where,asusual, thecontributions duetotheHe I andHe II lines
two-photoncontinuumiscalculatedusingtheformalismdescribed are added in the corresponding frequency bins. If a non-ionizing
byNussbaumer&Schmutz(1984);thedataofDrakeetal.(1969) energypacketistobere-emitted,thenitsfrequencymustbedeter-
areusedfor He I.Recombination linesbetweenlower levelsn=2 minedfromtheprobabilitydensityfunctionfornon-ionizingradia-
through8andupperlevelsn=3through15,forHI,andlowerlev- tiveenergy,whichisanalogoustoequation14.
elsn=2through16andupperlevelsn=3through30forHeII,are
calculatedasafunctionoftemperatureaccordingtothecaseBdata
2.7 TheIterativeProcedure
published by Storey&Hummer (1995). The He I recombination
linesarecalculatedasafunctionoftemperatureusingthedataof Aninitialguessofthephysicalconditionsinthenebularcells,such
Benjaminetal.(1999).Ingeneral,HeIsingletlinesfollowCaseB astheionizationstructure,electrontemperatureandelectronden-
whereastripletlinesfollowCaseA(asthereisnon=1levelforthe
triplets).Transitionstothe11SgroundstateofHeIproducelines 1 Resonancelineslongwardof912A˚ (e.g.CIVλλ1548,1550)may,infact,
whicharecapableofionizingHandlowionizationstagesofhigher
diffuseoutofthenebulaviaresonantscatteringandmayalsobeabsorbed
elements. In particular, the emissivities of the He I Lyman lines bydustduringsuchdiffusion.Atreatmentofdustgrainswillbeincluded
fromn=2throughn=5(Brocklehurst1972)andtheintercombina- infuturedevelopmentsoftheMocassincode,andsucheffectsmaythenbe
tionlinescorrespondingtothetransitions23S-11Sand23P-11Sare accountedfor.
6 Ercolanoet al.
sity,needstobespecifiedbeforethesimulationcanbegin.Proce- dances given by the final converged model solution to obtain the
duresinMocassinhavebeenconstructedsuchthatonlyaninitial lineemissivitiesforeachgridcell.Theluminosityofthenebulain
guessattheelectrontemperature(whichisinitiallysettoaconstant anygivenlinecanthenbecalculatedeasilybysummingtheemis-
valuethroughoutthenebula)mustbeincludedintheinputfile.Mo- sivityoftherequiredlineoverthevolumeofthenebula.
cassincanthenguessaninitialionizationstructureand,hence,the Acomparison oftheresultsobtainedusingthetwomethods
electrondensity.However,iftheoutputofaonedimensionalmodel described above, provides an indication of the level of accuracy
(oracombinationofmorethanoneofthem)isavailable,thereare achievedduringthesimulation,asthetwomethodswillgivecon-
also procedures built into Mocassin to map these onto the three- sistentresultsonlyifenoughenergypacketshavebeenusedtoyield
dimensionalCartesiangrid,byusingsimpleinterpolationroutines. good statisticsforevery line.Ingeneral, thesecond method (for-
Aone-dimensionalmodeoptionwasimplementedinMocassinfor malsolution)yieldsthemostaccurateresults,particularlyforweak
thispurpose.Severaltestshaveshownthatwhilethechoiceofthe lines,whichmayemit relativelyfewphotons. Forthebenchmark
initialconditionshas,ofcourse,noinfluenceonthefinalresultof casespresentedhere,reasonableaccuracywasdeemedtohavebeen
the simulation, it can, however, have an inpact on the number of achievedwhenthefluxesofthestrongesttransitionsobtainedusing
iterationsrequiredtoreachconvergence. Itishardtoquantifythe thepureMonteCarlomethodwerewithin10%ofthoseobtained
numberofiterationsrequiredforconvergencebyeachmethod,in using the formal solution. Bothmethods can also be used tocal-
particular,itdependsstronglyontheinitialtemperatureinputused culate line of sight results and to simulate long-slit observations.
inthefirstmethod,and,whenapplyingthesecondmethod,onthe However,justasforthecalculationoftheintegratedspectrum,the
deviationof theactual three-dimensional geometry fromthesim- formalsolutionmethodistobepreferred,asityieldsthemostac-
plifiedone-dimensional model used. However, withsufficient en- curateresults,particularlyfortheweakerlines.
ergypackets,thebenchmarkmodelsdescribedhereshouldbefully Inadditiontotheintegratedemergentspectrum,otheruseful
convergedinapproximatelyfifteentotwentyiterations.Astrategy comparisonswiththeobservationscanbecarriedout,e.g.projected
tospeedupthesimulationsisdescribedinSection3.1. imagesofthefinalmodelnebulainagivenlineoratagivencon-
Oncetheinitialconditionsarespecified,thefrequencydepen- tinuum frequency can be produced for arbitrary viewing angles.
dent total emissivities are calculated in each grid cell in order to These can be compared directly with nebular images obtained in
set up the probability density functions for re-emitted radiation, an appropriate filter. Mocassin computes and stores the physical
whichareusedforthedeterminationofthefrequencydistribution properties of the nebula, as well as the emissivities of the gas at
of the re-emitted energy-packets during the Monte Carlo simula- eachgridpoint;thesecanbefedintoIDLplottingroutinesinor-
tion.Theenergy packets arethen firedthrough thegridand their dertoproducemaps(Morisset etal.,2000). Also,byassuminga
trajectoriescomputed.Oncealltheenergypackettrajectorieshave velocityfield,linespectralprofilescanbeproduced,togetherwith
been computed, the Monte Carlo estimators for the mean inten- position-velocitydiagrams.Thesecanbecomparedwithobserva-
sityofthestellarandthediffuseradiationfieldscanbeobtained, tions,ifavailable,todeducespatio-kinematicinformationaboutthe
as described in Section 2.5. The ionization fraction and the elec- objectbeingstudied.MoreinformationabouttheoriginalIDLrou-
tron temperatures and densities must now be updated to be self- tinesisgivenbyMorissetetal.(2000)andMonteiroetal.(2000).
consistent withthecurrentestimatesoftheradiationfieldateach DetailsoftheactualapplicationtoMocassin’sgridfilesareavail-
gridpoint.Thismeanssolvingthelocalionizationbalanceandther- ableinacompanionpaperonthemodellingoftheplanetarynebula
malequilibriumequationssimultaneously.Theentireprocedureis NGC3918(Ercolanoetal.2002,PaperII).
repeateduntilconvergengeisachieved.Theconvergencecriterion At theend of each Monte Carloiterationthephysical quan-
that is used in this work is based on the change of the local hy- titieswhichcharacterisethegridarewrittenouttodiscintothree
drogenionizationstructurebetweensuccessiveiterations.Insome files,namelygrid1.out,grid2.outandgrid3.out.Thefirstfilecon-
cases,however,thisisnotasuitableconvergencecriterion(e.g.in tainsthelocalelectrontemperatureanddensityaswellasthegas
hydrogen-deficientenvironments),forthisreason,othercriteriaare densityateachgridcell,thesecondtheionizationstructureofthe
alsoimplementedinthecode(e.g.basedonthechangeofthelo- nebulaandthethirdanumberofmodelparameters,includingthe
calheliumionizationstructure,orofthelocalelectrontemperature numberofenergypacket tobeusedinthesimulation.Thesefiles
betweensuccessiveiterations),andthesecanbeeasilyselectedby areusedinconjuctionwithawarmstartdriver,whichallowsanin-
usingtheappropriateswitchesintheinputfile. terruptedsimulationtoberesumedfromtheendofthelastMonte
Carlosimulation.Thismeansthatonceasimulationhasbeeninter-
ruptedthenumberofenergypacketsused(andindeedothermodel
2.8 ComparisonoftheModelwithObservations
parameters, if required) can be adjusted, before the simulation is
Whenthemodelhasconvergedtoitsfinalsolution,theoutputspec- restarted,bymodifyingthefilegrid3.out.Thisfeaturecanbeused
trumcanbecomputedandcomparedwiththeresultsobtainedfrom tospeedupthesimulationsbyusingthefollowingapproach.The
other models or with observational data. The total luminosity of firstfewiterationsarerunusingalowernumberofenergypackets
the nebula emitted in various emission lineslongward of the Ly- thanactuallyneeded;so,forexample,iftheoptimumnumberofen-
manlimitcanbeobtainedbyusingtwomethods.Thefirstmethod, ergypacketsforagivenmodelis106,thenthefirstfewiterations
whichisonlyavailabletoMonteCarlocodes,consistsofsumming canbecarriedoutusingonly105 packets,hencereducingtherun
upthenumberofenergypacketsinthegivenline,N ,overthe timeforthesebyafactoroften.Thiswillresultinabout50%-60%
line
gridcells.Hence,thepoweremittedinthelineisgivenby of the grid cells converging; in general, the inner cells converge
morequickly,duetothelargernumberofsamplingunitsavailable
ε imaxjmaxkmax there(due mainlytogeometrical dilutionand tothe reprocessing
L = 0 N (x ,y ,z ) (15)
line ∆t line i j k ofenergypacketstonon-ionizingenergypackets).Atthispointthe
X X X
i=1 j=1 k=1 simulation isinterrupted and then resumed, after having adjusted
where ε0 isgiven byequation 2.Thesecond methodconsistsof thenumberofenergypacketstothefinalrequiredvalue(i.e.106,
∆t
usingthevaluesofthelocalelectrontemperatureandionicabun- inthepreviousexample).Finalconvergencewillthenbeachieved,
Mocassin 7
inmostcases,withinfourorfivefurtheriterations.Theactualnum-
Table1.Lexington2000benchmarkmodelinputparameters.
berofiterationsrequireddependsonthenumberofenergypackets
used: the larger the number of sampling quanta available at each
Parameter HII40 HII20 PN150 PN75
cell,thequickerthecellswillconvergetoasolution.Thenumbers
quotedabove,however,alsodependoneachparticularmodel’sge- L(BB)/1037(esregc) 308.2 600.5 3.607 1.913
ometryandopticalthickness. T(BB)/103K 40 20 150 75
Rin/1017cm 30 30 1 1.5
nH/cm−3 100 100 3000 500
He/H 0.10 0.10 0.10 0.10
2.9 GeneralArchitecture
C/H×105 22. 22. 30. 20.
TheMocassincodewaswrittenusingtheFortran90programming N/H×105 4. 4. 10. 6.
language. The code was developed and run initially on a Com- O/H×105 33. 33. 60. 30.
paq(Dec)XP1000witha500MHzCPUand1Gbofmemoryand Ne/H×105 5. 5. 15. 6.
Mg/H×105 - - 3. 1.
apreliminaryserialversionof thecodestillexists.Afullyparal-
Si/H×105 - - 3. 1.
lel version of the code has since been developed using Multiple
S/H×105 0.9 0.9 1.5 1.
ProcessesInterface (MPI) routines and itcurrently runs onaSil-
iconGraphicsOrigin2000machinewith24processors and6Gb
ofmemoryandaSUNMicrosystemsSunfireV880machinewith ElementalabundancesarebynumberwithrespecttoH.
16processorsand64Gbofmemory.MonteCarlosimulationsare,
bytheirnature,veryparallelizableproblemsand,indeed,Mocassin
canachievealinearspeed-up, i.e.aspeed-up thatisdirectlypro- codeandRubin’sNebula,treatthediffuseradiativetransferexactly.
portionaltothenumberofprocessorsused.Adetaileddescription Theothersusesomeversionsoftheoutward-onlyapproximationof
of all the Mocassin modules, input commands and output filesis varyingsophistication.Inthisapproximationalldiffuseradiationis
givenbyErcolano(2002,PhDThesis).Acopyofthecodeisavail- assumedtobeemittedisotropicallyintotheoutwardhalfofspace.
ablefromtheauthor([email protected])togetherwiththerelevant Thepredictedlinefluxesfromeachcodeforeachbenchmark
thesischapters. casearelistedinTables4to7,togetherwiththevolume-averaged
mean electron temperature, weighted by the proton and electron
densities, N N , <T[N ,N ] >, the electron temperature at the
p e p e
inneredgeofthenebula,T ,andthemeanratiooffractional
inner
3 APPLICATIONTOBENCHMARKCASES He+ to fractional H+, <He+>, which represents the fraction of
<H+>
Numerical simulationsofphotoionized nebulaeareverycomplex helium in the H+ region that is singly ionized. <T[NH+,Ne] >
andanumberoffactors,suchasnumericalapproximationsandas- and <He+> are calculated according to the following equations
sumptions, and thecomplexity of thecalculation itself,introduce <H+>
(Ferlandetal.1995)
a degree of uncertainty into the results. For this reason, it is im-
portant for modelers to have certain standards of comparison, in NeNpTedV
<T[N ,N ]>= (16)
ordertoidentifyproblemsintheircodesandtoreachanadequate p e R N N dV
e p
degree of accuracy in their calculation. A series of meeting have R
and
beenheld,beginninginMeudon,France,in1985(Pe´quignot1986)
aagnadininiLne2x0in0g0to(nP,e´Kqueingtnuoctkey,tfialr.st2i0n011)9,9i5n(oFredrelarntdoedteafil.n1e9a95s)etanodf <<HHe++>> = nn((HHe))R NNeNNHe+dVdV (17)
benchmarkcaseswhichcouldbeusedbyallphotoionizationmod- e p
R
elers to test their codes against. The benchmarks which resulted whereN andN arethelocalelectronandprotondensities,
e p
fromthesemeetingsinclude H II regions,planetary nebulae, nar- respectively,NHe+ isthedensityofHe+,andn(H)andn(He)are
rowlineregions(NLRs)ofAGNsandX-rayslabs.Mocassindoes thetotalhydrogenandheliumdensities.
nothave,atpresent,thecapabilitytotreatNLRsandX-rayslabs, Table1liststheinputparametersforallthebenchmarkmodels
assomerelevantphysicalprocesses,suchasComptonheatingand dicussedhere.AllthebenchmarkcaseslistedinTable1werecal-
innershellionization,arenotyetincluded.Forthisreason,onlythe culatedusingboththethree-dimensionalandtheone-dimensional
HIIregionandplanetarynebulabenchmarksareperformedinthis mode of Mocassin and both sets of results are included here for
work.Theexpansionofthecodetoincludehighenergyprocesses comparison. It is clear from Tables 4 to 7 that the results of the
isplannedinthefuture. three-dimensional and one-dimensional modes of Mocassin are
Resultsfromseveralothercodesareavailableforcomparison; consistentwitheachother.Thesmalldifferencesthatdoexistcan
theseareallone-dimensionalcodesand,apartfromdifferencesin be entirely attributed to the coarseness of the grids used for the
the atomic data used by each of them, their main differences lie three-dimensional calculations. Theaimof thebenchmarking de-
inthetreatment of thediffuseradiationfieldtransfer.Abrief de- scribed in this work is toassess the reliabilityof Mocassin in its
scriptionof each of these codesisgiven by Ferlandetal. (1995). fullythree-dimensionalmode,forthisreasontheone-dimensional
Althoughthemajorityofthesecodeshaveevolvedsomewhatsince mode results will not be included in the following performance
the1995Lexingtonmeeting,mostlyviatheupdatingoftheatomic analysis; moreover theinclusion of twosetsof resultsfromwhat
data sets and the inclusion of more and specialised physical pro- is,essentially,thesamecodewouldintroduceabiasinthemedian
cesses,theirbasicstructureshavestayedthesame.Thesevencodes andisolationfactorscalculationsdescribedbelow.Finally,toavoid
includedforcomparisonareG.Ferland’sCloudy(GF),J.PHarring- anyconfusion,anymentionofMocassinthroughouttherestofthis
ton’s code (PH), D. Pe´quignot’s Nebu (DP), T. Kallman’s XStar paperreferstothefullythree-dimensionalversionofthecode,un-
(TK), H. Netzer’s Ion (HN), R. Sutherland’s Mappings (RS) and lessotherwisestated.
R.Rubin’sNebula(RR).Onlytwoofthesecodes,theHarrington Figures1 and 2 show theelectron temperatures (toppanels)
8 Ercolanoet al.
andthefractionalionicabundancesofoxygen(middlepanels)and
Table2.Summaryofthenumberofenergypacketsneededfor>50%and
carbon(bottompanels)forthefourbenchmarkcasesanalysed.The >95%convergence (seetextforexplanation) foreachofthebenchmark
ionic abundances in every cell in the ionized region are plotted cases
againstradialdistancefromthestar.Theseplotsareinterestingnot
onlybecausetheyprovideaclearpictureoftheoveralltemperature
andionizationstructureofeachmodel,butalsobecausefromthe Case τedge nx×ny×nz Npackets
scatterofthedatapointsonecanestimatetheaccuracyofthefinal H0 He0 He+ >50% >95%
results.(Notethatsuchplotsareonlymeaningfulinthespherically HII40 4.79 1.15 177.8 13×13×13 5·105 5·106
symmetriccase.) HII20 2.95 1.13 91.2 13×13×13 5·106 5·107
Fourbenchmarkmodelnebulaewerecomputed,twoHIIre- PN150 34.0 6.87 57.9 13×13×13 3·105 3·106
gionsandtwoplanetarynebulae.Thesebenchmarksweredesigned PN75 1.16 0.24 31.5 13×13×13 4·105 4·106
tobeuncomplicated yet totest differentaspectsof themodelling
(seeFerlandetal.1995).Thenebulaearehomogeneousindensity
and, for simplicity, blackbodies are used as the ionizing sources
Table3.DeviationoftheMonteCarlomethodfromtheformalsolutionfor
insteadofmodelstellaratmopheres.
thepredictionofsomesignificantlinefluxesinthebenchmarkmodels.
Following the analysis of Pe´quignot (see Pe´quignotetal.
2001), isolation factors, if’s, were computed for each predicted
quantity in each case study. These are defined as the ratio of the
Line HII40 HII20 PN150 PN75
largesttothepenultimatelargestvalueofagiven output quantity
ortheratioofthepenultimatesmallestvaluetothesmallestvalue. Hβ 2.7% 9.5% 5.8% 2.8%
These ratios are computed with the intention to identify aberrant HeI5876A˚ 5.2% 6.3% 0.96% 4.5%
values.Alargeif canbeattributedtoanumberoffactors,butof- [NII]6584A˚ 7.6% 4.9% 8.5% 4.8%
tenthesecanbeattributedtoadifferenceintheatomicdataused [OII]5007A˚ 3.1% 12.0% 4.0% 1.1%
by each modeler. A list of the number of if’s larger than 1.01, [SIII]9532A˚ 5.8% 5.0% 2.0% 2.0%
1.03,1.10,1.30and2.00isgiveninTable8,foreachbenchmark.
Afteranalysing thebenchmark resultsobtained byallthemodel-
ers who participated in the Lexington workshop, Pe´quignotetal. requiredtoachieveagivendegreeofconvergence.Thereasonfor
(2001)suggestedthatanisolationfactorlargerthan1.30isindica- thiseffectisthatinasofterradiationfieldcasethenumberofen-
tiveofasignificantdepartureandapossibleproblem.Alargenum- ergypacketsemittedatwavelengthsshorterthantheLymanlimit
berofoccurencesofif′s > 1.30shouldeitherhaveanacceptable will be less than in the case of a harder radiation field. A larger
explanationorleadtocorrectionstothecode. total number of energy packets then needs to be used in order to
Thenumberofresultsnotpredictedbyanygivencodeislisted obtainanumber ofionizingphotonsadequatetoproperlysample
intheNopredrowofTable8.Pe´quignotetal.(2001)alsonoted,in thenebula.TheaimofTable2ismerelytoprovidesomegeneral
theproceedingsoftheNovember2000Lexingtonmeeting,thatthe guidelinesforselectingtheappropriatenumberofenergypackets
lackofapredictionforaparticularobservablemaysimplyreflecta foraparticularsimulation;however,asstatedbefore,theoptimum
lackofinterestbythemodellerinit;ontheotherhand,afrequent numbershouldbedeterminedforeachgivenmodel,particularlyin
occurenceofNopredmayalsoindicatelimitationsinthepredictive non-sphericallysymmetriccases.
powerofagivencode.
AsarguedbyPe´quignotetal.(2001),alargeerrorcanbein-
troducedwhentheaverageoverasmallsamplecontaininganum- 3.2 BenchmarkResults
berofaberrantvaluesistaken.Inordertominimisethiserror,me-
The Lexington/Meudon Standard H II region model (HII40) was
dianvaluesarecalculatedinsteadofaveragesandthesearegiven
the first benchmark to be run and some very preliminary results
foreachobservablelistedinTables4to7,inthecolumnlabelled
have already been presented, at the November 2000 Lexington
Med.ThemediansarecalculatedtotheprecisionshowninTables4
meeting (Ercolano 2001; Pe´quignotetal. 2001). However, those
to7.Table9liststhenumberofmedianvaluesscoredbyeachcode
resultswereproducedwhenMocassinwasstillunderdevelopment
foreachbenchmark,i.e.thenumberoftimesthecodewastheclos-
andshouldthereforeonlybeconsideredasasnapshotofthecode
esttothemedianvalue.Whenamedianvalueissharedbytwoor
at that particular stage. The code has evolved considerably since
morecodesthescoreisgiventoeachone,thereforethesumofthe
theNovember 2000 Lexington meeting and thenewer resultsare
medianvaluesscoredbyallthecodesishigherthanthenumberof
presentedinthissection(seeTable4).
observables(thecolumnlabelledTotalinTable9).
Table 3 shows the results of a comparison between the line
fluxesobtainedbyMocassinusingtheformalsolutionmethodand
thoseobtainedusingtheMonteCarlomethod(seeSection2.8)for
3.1 SamplingRequirements
someofthemoresignificantlinesinthebenchmarkcases.Itisclear
Table2liststheopticaldepthsattheionizationthresholdfrequen- thattheresultsshownagreewell;however,asexpected,largerdis-
ciesforH0,He0andHe+,attheouteredgeofthegrids,forthefour crepancieswerefoundfortheweakerlines,whoselowernumbers
benchmarkmodelsanalyzedhere.Foreachmodel,thenumberof ofenergypacketsyieldloweraccuracystatistics.
gridpointsisalsogiven(column 5),together withthenumber of
energy packets used, N , according to the two-step stategy
packets
described above, first toachieve convergence in50%-60% of the
3.2.1 TheHII40benchmark
totalnumberofgridcells(> 50%,column6)andthentoachieve
totalconvergence(>95%,column7).Table2showsthatthesofter Mocassinscoredeightif’s>1.01fortheHII40benchmarkmodel
theionizingradiationfield,thelargerthenumberofenergypackets (Table 8); only three of these, however, had values greater than
Mocassin 9
Table4.Lexington2000StandardHIIregion(HII40)benchmarkcasereults.
Line Median GF HN DP TK PH RS RR BE
3-D 1-D
Hβ/1037erg/s 2.05 2.06 2.02 2.02 2.10 2.05 2.07 2.05 2.02 2.10
Hβ4861 - 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
HeI5876 0.116 0.119 0.112 0.113 0.116 0.118 0.116 - 0.114 0.112
CII]2325+ 0.144 0.157 0.141 0.139 0.110 0.166 0.096 0.178 0.148 0.126
CII1335 0.082 0.100 0.078 0.094 0.004 0.085 0.010 - 0.082 0.084
CIII]1907+1909 0.070 0.071 0.076 0.069 0.091 0.060 0.066 0.074 0.041 0.041
[NII]122.µm 0.034 0.027 0.037 0.034 - 0.032 0.035 0.030 0.036 0.034
[NII]6584+6548 0.730 0.669 0.817 0.725 0.69 0.736 0.723 0.807 0.852 0.786
[NII]5755 .0054 .0050 .0054 .0050 - .0064 .0050 .0068 .0061 .0054
[NIII]57.3µm 0.292 0.306 0.261 0.311 - 0.292 0.273 0.301 0.223 0.229
[OI]6300+6363 .0086 .0094 .0086 .0088 .012 .0059 .0070 - .0065 .0080
[OII]7320+7330 0.029 0.029 0.030 0.031 0.023 0.032 0.024 0.036 0.025 0.022
[OII]3726+3729 2.03 1.94 2.17 2.12 1.6 2.19 1.88 2.26 1.92 1.75
[OIII]51.8µm 1.06 1.23 1.04 1.03 0.99 1.09 1.06 1.08 1.06 1.09
[OIII]88.3µm 1.22 1.12 1.06 1.23 1.18 1.25 1.23 1.25 1.22 1.26
[OIII]5007+4959 2.18 2.21 2.38 2.20 3.27 1.93 2.17 2.08 1.64 1.70
[OIII]4363 .0037 .0035 .0046 .0041 .0070 .0032 .0040 .0035 .0022 .0023
[OIV]25.9µm .0010 .0010 .0010 .0010 .0013 .0013 .0010 - .0010 .0010
[NeII]12.8µm 0.195 0.177 0.195 0.192 - 0.181 0.217 0.196 0.212 0.209
[NeIII]15.5µm 0.322 0.294 0.264 0.270 0.35 0.429 0.350 0.417 0.267 0.269
[NeIII]3869+3968 0.085 0.084 0.087 0.071 0.092 0.087 0.083 0.086 0.053 0.055
[SII]6716+6731 0.147 0.137 0.166 0.153 0.315 0.155 0.133 0.130 0.141 0.138
[SII]4068+4076 .0080 .0093 .0090 .0100 .026 .0070 .005 .0060 .0060 .0057
[SIII]18.7µm 0.577 0.627 0.750 0.726 0.535 0.556 0.567 0.580 0.574 0.569
[SIII]33.6µm 0.937 1.24 1.43 1.36 0.86 0.892 0.910 0.936 0.938 0.932
[SIII]9532+9069 1.22 1.13 1.19 1.16 1.25 1.23 1.25 1.28 1.21 1.19
[SIV]10.5µm 0.359 0.176 0.152 0.185 0.56 0.416 0.388 0.330 0.533 0.539
103×∆(BC3645)/A˚ 5.00 4.88 - 4.95 - 5.00 5.70 - 5.47 5.45
Tinner/K 7653 7719 7668 7663 8318 7440 7644 7399 7370 7480
<T[NpNe]>/K 8026 7940 7936 8082 8199 8030 8022 8060 7720 7722
Rout/1019cm 1.46 1.46 1.46 1.46 1.45 1.46 1.47 1.46 1.46 1.49
<He+>/<H+> 0.767 0.787 0.727 0.754 0.77 0.764 0.804 0.829 0.715 0.686
GF:G.Ferland’sCloudy;PH:J.PHarringtoncode;DP:D.Pe´quignot’sNebu;TK:T.Kallman’sXStar;HN:H.Netzer’s
Ion;RS:R.Sutherland’sMappings;RR:R.Rubin’sNebula;BE:B.Ercolano’sMocassin.
1.3. Amongst these, if’s > 1.10 are obtained for the [O III] 3.2.3 ThePN150benchmark
5007+4959 (if = 1.18) and for [O III] 4363 (if = 1.45); the ra-
tio of these lines is often used as a temperature diagnostic (see,
The optically thick high-excitation planetary nebula (PN150) is
forexample,Osterbrock1989,pages119-125).Mocassinpredicts
the most demanding of the benchmark cases in terms of physi-
jλ4959+jλ5007 =745.4,thisvalueishigherthanthevalueobtained
jλ4363 cal processes and atomicdata required. Mocassin’s score for this
by theother codes, infact median value obtained for theratioof
modelwasverygood(Table8),obtainingonlysixif’s,withnone
theselinefluxesbytheothercodesisequalto589.2.Thisisfully
of those being higher than 1.3 and only one slightly higher than
consistent with Mocassin predicting a slightly lower temperature
1.1 (C II λ1335, if = 1.13). As has already been discussed by
(if=1.027)forthisbenchmarkthandotheothercodes.
Pe´quignotetal.(2001),thereseemstobeadichotomybetweenthe
Finally,thenumberofmedianvaluesobtainedforthisbench-
GF,HNandDPmodels(and,now,alsotheBEmodel)ontheone
markcaseisten,whichcomparesverywellwiththeothercodes’
hand, and the TK, PH and RS models on the other. The former
medianscores,rangingfromthreetoten(seeTable9).
groupobtainedveryfewif’slargestthan1.1,indicativeofatighter
agreement.Thiscoherencecanprobablybeattributedtoamorere-
cent updating of atomic dataand amore careful treatment of the
3.2.2 TheHII20benchmark diffuse radiation field transfer. These four codes also obtained a
largerHβ fluxfor thismodel, whichcan probably beascribed to
Mocassin scored seven if’s for the low excitation H II region secondaryphotonsfromheavyions.PHistheonlyclassicalcode
(HII20) benchmark model. None of these, however, has a value herewithafullyiterativesphericallysymmetricradiativetransfer
greater than 1.3. As in the HII40 nechmark case, the mean tem- treatment(sinceRRonly computed H II regions); thiscould also
perature,weightedbyNpNe,predictedbyMocassinforthismodel bethereasonfortherelativelylargernumberofif’sscoredbythe
isalsoslightlylower(if=1.034)thantheothermodels’predictions. PHcodeforthismodel.
FivemedianvalueswereobtainedbyMocassinforthisbench-
markcase,whiletheothercodesscoredbetweenthreeandeleven The score for median values obtained by Mocassin for the
(seeTable9). PN150opticallythickplanetarynebulaisextremelygood,obtain-
10 Ercolanoet al.
Figure1.Electrontemperature(toppanels)andthefractionalionicabundancesofoxygen(middlepanels)andcarbon(bottompanels),asafunctionofnebular
radius,fortheHIIregionbenchmarkcasesHII40(left-handpanels)andHII20(right-handpanels).