DART-Ray: a 3D ray-tracing radiative transfer code for calculating the propagation of light in dusty galaxies

1 MB·

by  G. Natale

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed8January2014 (MNLATEXstylefilev2.2)
DART-RAY: a 3D ray-tracing radiative transfer code for calculating
the propagation of light in dusty galaxies
4 G. Natale1, C. C. Popescu1,2, R. J. Tuffs2, D. Semionov1
1
0
2 1JeremiahHorrocksInstitute,UniversityofCentralLancashire,Preston,PR12HE,UK
2MaxPlanckInstitutefu¨rKernphysik,Saupfercheckweg1,D-69117Heidelberg,Germany
n
a
J

7 8January2014

]
M
I ABSTRACT
h. We present DART-RAY, a new ray-tracing 3D dust radiative transfer (RT) code designed
p specificallytocalculateradiationfieldenergydensity(RFED)distributionswithindustygalaxy
- modelswitharbitrarygeometries.Inthispaperweintroducethebasicalgorithmimplemented
o
r inDART-RAYwhichisbasedonapre-calculationofalowerlimitfortheRFEDdistribution.
t Thispre-calculationallowsustoestimatethe extentofregionsaroundtheradiationsources
s
a withinwhichthesesourcescontributesignificantlytotheRFED.Inthisway,ray-tracingcal-
[
culationscanberestrictedtotakeplaceonlywithintheseregions,thussubstantiallyreducing

1 thecomputationaltimecomparedtoacompleteray-tracingRTcalculation.Anisotropicscat-
v teringisincludedinthecodeandhandledinasimilarfashion.Furthermore,thecodeutilizes
7 a Cartesian adaptivespatial gridand an iterative methodhas beenimplementedto optimize
3
theangulardensitiesoftheraysoriginatedfromeachemittingcell.Inordertoverifytheac-
4
1 curacyoftheRTcalculationsperformedby DART-RAY, wepresentresultsofcomparisons
. withsolutionsobtainedusingtheDUSTY1DRTcodeforadustshellilluminatedbyacen-
1
tral pointsource and existing 2D RT calculationsof disc galaxieswith diffusely distributed
0
4 stellar emission and dust opacity. Finally, we show the application of the code on a spiral
1 galaxymodelwithlogarithmicspiralarmsinordertomeasuretheeffectofthespiralpattern
: ontheattenuationandRFED.
v
i
X
r
a 1 INTRODUCTION
Interstellardustisofprimaryimportanceindeterminingthespectralenergydistribution(SED)oftheradiationescapingfromgalaxiesat
wavelengths rangingfromtheultraviolet(UV) tothesubmillimetre(submm) and radio. Dustattenuates andredistributesthelight,origi-
natingmainlyfromstarsand, ifpresent, fromanactivegalacticnucleus (AGN),by eitherabsorbing or scatteringphotons. Theabsorbed
luminosityisthenre-emittedintheinfraredregime.Thedustmaybesituatedincomplexgeometrieswithrespecttothesesources,affecting
theobservedstructureofthegalaxyateachwavelengthaswellasitsintegratedSED.Modelingthepropagationoflightwithinrealgalaxies
isthusachallenging task. Nevertheless, itisessential todo suchmodeling, ifphysical quantitiesof interest,such asthedistributionand
propertiesofthestellarpopulationsandtheinterstellarmedium(ISM)astracedbydustandtheinterstellarradiationfields,aretobederived
frommultiwavelengthimagesandSEDs.
Takingadvantageoftheapproximatecylindricalsymmetryofgalaxies,2Ddustradiativetransfer(RT)models,suchastheonepresented
byPopescuetal.(2011),alreadycontainthemainingredientsneededtopredictintegratedgalaxySEDs,averageprofiles,dustemissionand
attenuationforthecaseofnormalstar-formingdiscgalaxies.However,thereareanumberofreasonswhy3DdustRTcodesaredesirable.
First,spiralgalaxies,althoughwellmodeledwith2Dcodes,showthepresenceofmultipleandirregularfeaturessuchasspiralstructures,
bars,warpsandlocalclumpinessoftheISM.Also,galaxiesmayhostacentralAGNwhosepolaraxismaynotbealignedwiththatofthe
galaxy.Formergersorpost-mergergalaxiesthereisclearlynofundamentalsymmetryofthedistributionofstarsanddust.Finally,solutions
forthedistributionofstarsandISMprovidedbynumericalsimulationsofformingandevolvinggalaxiesgenerallyrequireprocessingwith
a3DRTcodeinordertopredicttheappearanceindifferentbands.
2 Nataleetal.
Themainchallengeinrealizing3DsolutionsofthedustRTproblemisthecomputationalexpense.Thestationary3DdustRTequationis
anon-localnon-linearequation:non-localinspace(photonspropagatewithintheentiredomain),direction(duetoscattering,absorption/re-
emission)andwavelength(absorption/re-emission).Evenusingarelativelycoarseresolutionineachofthesixfundamentalvariables,namely
thethreespatialcoordinates,thetwoanglesspecifyingtheradiationdirectionandthewavelength,solvingthe3DdustRTproblemrequire
animpressiveamountofbothmemoryandcomputationalspeed,atthelimitsofthecapabilitiesofcurrentcomputers.
PossiblythequickestwaytocalculateanimageofagalaxyindirectandscatteredlightinaparticulardirectionisbyusingMonteCarlo
(MC)methods(includingmodernaccelerationtechniques,seeSteinackeretal.2013).ThereisarichhistoryofapplicationsofMCcodes
todustRTproblems,startingwiththepioneeringworksofe.g.Mattila(1970),Roark,Roark&Collins(1974),Witt&Stephens(1974)and
Witt(1977).InthefollowingdecadestheMCRTtechniquewasfurtherdevelopedbymanyauthorssuchase.g.Witt,Thronson&Capuano
(1992),Fischer,Henning&Yorke(1994),Bianchi,Ferrara&Giovanardi(1996),Witt&Gordon(1996)andDullemond&Turolla(2000).
Nowadays, thismethodcanbeconsideredasthemainstreamapproach to3D dustRTcalculations(seee.g.Gordon etal.2001, Ercolano
etal.2005,Jonsson2006,Bianchi2008,Chakrabarti&Whitney2009,Baesetal.2011,Robitaille2011,butalsoseetable1ofSteinacker
etal. 2013 forarecent listof published 3D dust RTcodes). TheMCapproach todust RTconsistsof asimulationof thepropagation of
photonswithinadiscretizedspatialdomain,basedonaprobabilisticdeterminationofthelocationofemissionofthephotons,theirinitial
propagation direction, theposition where an interactionevent (absorption or scattering) occurs and the new propagation direction after a
scatteringevent.Thus,theMCtechniquemimicscloselytheactualprocessesoccurringinnaturewhichshapetheappearanceofgalaxiesin
UV/opticallight.However,sinceitisbasedonaprobabilisticapproachtodeterminethephotonpropagationdirections,anRTMCcalculation
doesnotnecessarilydeterminetheradiationfieldenergydensity(RFED)accuratelyintheentirevolumeofthecalculation.Thereasonis
thatregionswhichhavealowprobabilityofbeingilluminatedarereachedbyonlyfewphotonsunlessthetotalnumberofphotonsinthe
RTrunissubstantiallyincreased.Nonetheless,inthecaseofdiscgalaxies,accuratecalculationofradiationfieldintensitiesthroughoutthe
entirevolumeisneeded,inparticularforthecalculationofdustemission.Indeed,far-infrared/submmobservationsofspiralgalaxiesshow
thatmostofthedustemissionluminosityisemittedlongwardsof100µm(seee.g.Sodroskietal.1997,Odenwaldetal.1998,Popescuet
al.2002,Popescu&Tuffs2002,Daleetal.2007,2012,Bendoetal.2012)throughgrainssituatedinthediffuseISMwhicharegenerally
locatedatveryconsiderabledistancesfromthestarsheatingthedust.
AnothermethodtosolvetheRTproblemingalaxies,alternativetothemainstreamMCapproach,isbyusingaray-tracingalgorithm.
Thismethodconsistsinthecalculationofthevariationoftheradiationspecificintensityalongafinitesetofdirections,usuallyreferredto
as”rays”.Ray-tracingalgorithmscanbespecificallydesignedtocalculateradiationfieldintensitiesthroughouttheentirevolumeconsidered
intheRTcalculation.Also,itshouldbepointedoutthatMCcodesalreadymakelargeuseofray-tracingoperations(seeSteinackeretal.
2013).Itisthusinterestingtopursueinthedevelopingofpureray-tracing3DRTcodes,whichcanbesufficientlyefficientforthemodelling
ofgalaxieswith3Darbitrarygeometries,ifappropriateaccelerationtechniquesareimplemented.SimilartoMCcodes,ray-tracingdustRT
codes have had arichhistory inastrophysics (seee.g. Hummer &Rybicki 1971, Rowan-Robinson 1980, Efstathiou& Rowan-Robinson
1990,Siebenmorgenetal.1992,Semionov&Vansevic˘ius2005,2006).Applicationtoanalysisofgalaxiesstartedwiththe2DcodeofKy-
lafis&Bahcall(1987).Althoughoriginallyimplementedonlyforthecalculationofopticalimages(seealsoXilourisetal.1997,1998,1999),
thisalgorithmwaslateradaptedbyPopescuetal.(2000)forthecalculationofradiationfieldsandwascoupledwithadustemissionmodel
(includingstochasticheatingofgrains)topredictthefullmid-infrared(MIR)/FIR/submmSEDofspiralgalaxies(seealsoMisiriotisetal.
2001,Popescuetal.2011).Thusfar,extensionsoftheray-tracingtechniqueto3Dhavebeenimplementedbutarespecificallydesignedfor
solvingtheRTproblemforstarformingclouds(e.g.Steinackeretal.2003,Kuiperetal.2010),heatedbyfewdominantdiscretesources,
ratherthanforveryextendeddistributionsofemissionanddustasencounteredingalaxies.
Inthispaper,wepresentDART-RAY1,anewray-tracingalgorithmwhichisoptimizedforthesolutionofthe3DdustRTproblemfor
galaxieswitharbitrarygeometriesandmoderateopticaldepthatoptical/UVwavelengths2.Themainchallengefacedbythismodelisthe
constructionofanefficientalgorithmfortheplacingofraysthroughoutthevolumeofthegalaxy.Infact,acompleteray-tracingcalculation
betweenallthecells,usedtodiscretiseamodel,isnotaviableoption,sinceitisbyfartoocomputationallyexpensiveevenforrelatively
coarsespatialresolution.Ouralgorithmcircumventstheproblembyperforminganappropriatepre-calculation,whosegoalistoprovidea
lowerlimittotheRFEDdistributionthroughoutthemodel.Inthisway,therayangulardensityneededintheactualRTcalculationcanbe
dynamicallyadjustedsuchthattheraycontributionstothelocalRFEDarecalculatedonlywithinthefractionofthevolumewherethese
contributionsarenotgoingtobenegligible.
Furthermore, the code we developed can be coupled with any dust emission model. Applications of the 3D code for calculation of
infraredemissionfromstochasticallyheateddustgrainsofvarioussizesandcomposition,includingheatingofPolycyclicAromaticHydro-
carbonsmolecules,utilizesthedustemissionmodelfromPopescuetal.(2011),andwillbegiveninafuturepaper.
Thepaper isstructured asfollows.In§2weprovide some background informationand motivation behind our particular ray-tracing
1 Thenameofthecodecanbeseenastheacronymfor“DustAdaptiveRadiativeTransferRay-tracing”
2 Theactualversionofthecodedoesnotconsiderthedustabsorption/scatteringoflightemittedbydustatotherpositions(socalled“dustself-heating”).This
effectcanbeneglectedincasethegalaxiesareopticallythinatinfraredwavelengths
DART-RAY 3
solutionstrategy.In§3wegiveatechnicaldescriptionofourcode.In§4weprovidesomenotesonimplementationandperformanceofthe
code.In§5wecomparesolutionsprovidedbyourcodewiththosecalculatedbythe1DcodeDUSTYandthe2DRTcalculationsperformed
byPopescuetal.(2011).In§6weshowtheapplicationofthecodeonagalaxymodelincludinglogarithmicspiralarms.Asummarycloses
thepaper.AlistofdefinitionsforthetermsandexpressionsusedthroughoutthepapercanbefoundinTable1.
2 BACKGROUNDANDMOTIVATIONS
Inthissectionwedescribethebasiccharacteristicsofthetime-independent3DdustRTequation,brieflyintroducetheray-tracingapproach
usedinourcodeandprovidethemainmotivationsbehindoursolutionstrategy.Finally,wedescribethemainstepsofthenewalgorithmwe
presentinthiswork.Amuchmoretechnicaldescriptionofourcodecanbefoundin§3.
2.1 Dustcontinuum3DRadiativeTransfer:Ray-tracingandSolutionStrategy
Givenaninputdistributionofstellarluminosityanddustmass,solvingtheRTproblemrequiresinprincipletheresolutionofthefollowing
equationforthespecificintensityI(λ,x,n),whichrepresentstheluminosityperunitarea,solidangleandwavelengthintervalpropagating
atpointxintothedirectionn:
n∇xIλ(x,n)=−kλ(x)ρ(x) Iλ(x,n)−ωλ Φλ(n,n′)Iλ(x,n′)dΩ′ +jλ(x) (1)
(cid:20) ZΩ (cid:21)
wherek (x)isthetotalextinctioncoefficientperunitmassofdust(includingbothabsorptionandscattering),ρ(x)isthedustmassdensity,
λ
ω isthealbedo,definedsuchthatω ×k (x)givesthefractionofextinctionduetoscattering,Φ (n,n′)isthescatteringphasefunction,
λ λ λ λ
whichgivestheprobabilityforradiationcomingfromdirectionn′ tobescatteredintodirectionn,andj (x)isthedistributionofstellar
λ
volumeemissivity3.Thefirsttermontheright-handsideofEq.1actstoreducetheradiationspecificintensityI(λ,x,n)byaquantitythat
isproportionaltotheradiationintensityitself.Thesecondterminsteadactsasasourcetermandgivesapositivecontributiontotheradiation
intensitybyaddingthelightcomingfromalldirectionstopointxandthenscatteredintothedirectionn.Thethirdtermgivesapositive
contributiontothepropagatingradiationspecificintensity,whichisduetostellaremissionanditisassumedtobeisotropicateachposition
x.Sincein3DRTtherearesixindependentvariables,namelywavelength,threespatialcoordinatesandtwoangulardirections,thesolution
vector for I (x,n) can be extremely large, making the problem very challenging also from the point of view of memory requirements,
λ
apartfromthedifficultyofsolvingtheintegro-differentialequationitselfinthreedimensions.Notealsothatwedidnotincludeatermthat
representsthere-emissionbydust,importantatinfraredwavelengths.Theresolutionalgorithmpresentedinthisworkisdesignedonlyto
handlethepropagationofdirectandscatteredlightfromstellarpopulationsanddoesnotconsidertheself-heatingofdust.
Instead of seeking to obtain a solution for I (x,n), the main aimof this work is toconstruct an algorithm optimized toderive the
λ
radiationfieldenergydensityU (hereafteralsorefereedtoasRFED)ateachpositionx,inordertobeabletocalculatesuccessivelythedust
λ
emissionspectraassuminglocalenergybalancebetweendustradiativeheatingandemission4
IntermsofI (x,n),onecanexpressU as:
λ λ
I (x,n)dΩ
U (x)= λ (2)
λ c
R
InordertocalculateU (x)ataspecificpointx,onecansimplysumupthecontributionsδU providedbytheradiationcomingfrom
λ λ
alltheemittingsourcestothevalueofU (x)atthatposition.AnumericalmethodtocalculateU (x)atanypositioncanbeimplementedby
λ λ
considering“rays”originatingfromeachemittingsourceandpropagatingthroughoutthewholevolumeconsideredinthecalculation.Along
eachraypathonefollowsthevariationoftheradiationintensityandonecanthuscalculatetheδU contributionsatafinitesetofpositions
λ
(thissolutiontechniqueisamongthoseknownas“ray-tracing”methods).Morespecifically,asolutionalgorithmonecouldusetoderivethe
distributionofU (x)forasinglewavelengthλ,givenaninput3Ddistributionofstellarluminosityanddustmass,isthefollowing.
λ
First,theentiremodelissubdividedinanadaptivegridofcubiccellsandtoeachcelloneassignstheaveragevaluesforboththedust
densityandstellarvolumeemissivitywithinthecellvolume.Acellforwhichtheaveragevalueofstellarvolumeemissivityishigherthan
3 Throughoutthetextby“volumeemissivity”wewillalwaysmeantheluminosityperunitvolumeperunitsolidangleandperunitwavelengthintervalof
thestellarradiationateachposition.Tonotbeconfusedwiththeemissivitycoefficientusedtocharacterisetheemissionpropertiesofe.g.gasordust.
4 Forexample,thelatterconditionforasinglegrainstochasticallyheatedcanbeexpressedas:
Z Qλ,absZ Bλ(T)P(T)dTdλ=(c/4π)Z Qλ,absUλdλ
whereQλ,absisthegrainabsorptionefficiencyatwavelengthλ,Bλ(T)isthePlanckfunctioncalculatedfordusttemperatureT andP(T)istheprobability
forthedustgraintohaveatemperatureequaltoT.Numericalmethods,suchastheonepresentedinGuhathakurta&Draine(1989),allowustoderiveP(T)
andthereforethedustemissionspectra,oncetheabsorptionefficiencyQλ,absandtheRFEDUλareknown..
4 Nataleetal.
Table1.Tablesoftermsanddefinition.Thesubscriptλdenotesadependencefromthewavelengthoftheradiation.
Term Definition
AEM Projectedareaofanemittingcell
AINT Projectedareaofanintersectedcell
fL Inputparameterneededtosetthethresholdvaluefor∆Lλbelowwhichscatteringitera-
tionsarestopped
fU InputparameterneededtosetthethresholdvalueforδUλbelowwhichraysarestopped
gλ Henyey-Greensteinscatteringphasefunctionparameter
Iλ Radiationspecificintensity
<Iλ> AverageofIλoverthepathcrossedbyaraywithinanintersectedcell
Iλ,i Radiationspecificintensityofaraybeforecrossingacelli
Iλ,esc Radiationspecificintensityofarayonceithasreachedthemodelborder
Iλ,ext Amountofradiationspecificintensityofarayextinctedwithinanintersectedcell
δIλ,sca(θ,φ) Raycontributiontothespecificintensityoftheradiationscatteredintodirection(θ,φ),
seeFig.6
jλ Volumeemissivity,thatis,theluminosityemittedatacertainpositionperunitvolume,
unit solid angle and unit wavelength interval (also jλ,c when refereed to the average
volumeemissivitywithinacell)
κλ Dustextinctioncoefficientperunitdustmass(alsoκλ,absorκλ,scawhenreferredtothe
extinctioncoefficientsduetodustabsorptionorscattering)
Lλ Totalstellarluminositydensityofthemodel
∆Lλ Amountofluminositystilltobeprocessedduringscatteringiterations
Lλ,ray Luminosityassociatedwitharaybeam
Nrays Inputparameterspecifyingtheminimumnumberofrayscrossingacellwithinthefully
sampledregion
Φλ Scatteringphasefunction
ρ Dustmassdensity(alsoρcwhenreferredtotheaveragedustmassdensitywithinacell)
∆r Raypathwithinanintersectedcell
SCATTEN(θ,φ) Arraystoringtheluminosityoftheradiationscatteredbydustycellsintoafinitesetof
solidangles
τλ Opticaldepth
Uλ Radiationfieldenergydensity(alsoreferredtoasRFED)
δUλ ContributiontothelocalRFEDcarriedbyaray(alsoUλ,INTwhenreferredtothecon-
tributiontoanintersectedcellRFED)
Uλ,LL LowerlimittotheRFED
UTEMP TemporaryarrayusedtostoreRFEDcontributionsfromanemittingcellthroughoutthe
model
Uλ,FINAL ArraystoringtheRFEDdistributionwhichisoutputbythecode
VINT Intersectedcellvolume
ωλ Albedo
ΩHP,EM SolidangleassociatedwiththeHEALPixsphericalpixelsusedtodefinethedirectionsof
raysfromanemittingcell
ΩHP,INT SolidangleassociatedwiththeHEALPixsphericalpixelsusedtodefinethedirectionsof
theradiationscatteredbyanintersectedcell
ΩHP,MS SolidangleassociatedwithanHEALPixmainsector(seeFig.4)
ΩINT SolidanglesubtendedbyprojectedareaoftheintersectedcellAINT(seeFig.6
“Dustycell” Acellwheretheaveragevalueofdustdensityishigherthanzero
“Dustself-heating” Thedustabsorptionofradiationemittedbydustatotherpositions(notincludedinthis
versionofthecode)
“Emittingcell” Acellwheretheaveragevalueofthestellarlightorscatteredlightvolumeemissivityis
higherthanzero
“Escapingradiation” Thedirect orscattered stellar radiation propagating outside the borders ofthe volume
consideredintheRTcalculation
“Fullsampling”(ofaregion) Theprocessoflaunchingenoughraysfromasource,suchthatallthecellswithinaregion
areintersectedbymultiplerays
“Intersectedcell” Acellintersectedbyaray
“Leafcell” Acellofthe3Dadaptivegridwhichisnotfurthersubdivided
“Lostluminosity” Amount of stellar luminosity not considered in the RT calculation (to be kept low to
guaranteeapproximateenergybalance,see§4)
MC MonteCarlo
RFED RadiationFieldEnergyDensity
RT RadiativeTransfer
“Sourceinfluencevolume” ThefractionofmodelvolumewithinwhichasourcecontributessignificantlytotheRFED
“Volumeemissivity” seejλ
DART-RAY 5
zerocanbetreatedapproximately asadiscreteradiationsource. Inthefollowing, wewillrefertothiskindof cellasan“emittingcell”.
Similarly,cellswithaveragedustdensityhigherthanzerowillbereferredtoas“dustycells”.
Then, one performs ray-tracing for a large set of directions originating from the centres of the emitting cells. That is, one follows
thevariationof thespecificintensityI (x,n)fromthecellcentresalongrayscorresponding toeach directionn.Whilefollowingaray,
λ
one considers the increase of radiation intensity I (x,n) due to the stellar volume emissivity in the cell originating the ray but not in
λ
theintersectedcells,whereonly thedecrease of intensityduetodust absorption and scatteringisconsidered. Thisallowsustocalculate
separatelythecontributionsδU providedonlybytheemittingcelloriginatingtheraytothefinalvalueofU (x)inallthecellsintersected
λ λ
bythesameray5.Whenarayintersectsacellidifferentfromtheoriginalcell,thenewvalueofthespecificintensitywillthenbe:
I =I e−kρc,i∆r (3)
λ,i+1 λ,i
whereρ isthecelldustdensityand∆risthecrossingpath.Foreachray-cellintersectiononecalculatesanappropriateaverageofthevalue
c
ofI withintheraycrossingpathand,thus,thecontributionδU tothelocalvalueofU (x)byusingadiscreteversionofEq.2.Whenall
λ λ λ
theraysfromoneemittingcellhavebeenprocessed,theray-tracingisperformedfromanotheremittingcellandsoonuntilalltheemitting
cellshavebeenconsidered.Scatteredradiation,whoseintensityinafinitenumberofdirectionshasalsobeenstoredlocallyforeachray-cell
intersection,isthenprocessedinasimilarfashion.
Ifonecreatesagridsufficientlyfineinresolutionandlaunchesasufficientlylargenumberofrays,thecalculatedcellRFEDvaluesare
veryclosetotheexactvaluesofU (x)atthecellcentresandcanbeusedtocalculatethedustheating.Furthermore,thedescribedprocedure
λ
is extremely flexible and capable to handle completely arbitrary 3D distributions of stellar luminosity and dust mass. Unfortunately, the
implementationoftheprocedureinthesimpleformdescribedaboveistoocomputationallyexpensive(scalingapproximatelyasN5/3,with
N beingthetotalnumberofcells,inthecaseofuniformspatialsampling),consideringalsothatthecalculationhastobeperformedinan
iterativewayforthescatteredlightandfordifferentwavelengths.Nonetheless,forwellmixedemitters–absorbersgeometries,suchasinthe
caseofgalaxydust–stellardistributions,thefullray-tracingcalculationforrayspropagatingthroughouttheentiremodelfromeachemitting
celldoesnothavetobenecessarilyperformedtoobtainareasonablyaccuratesolutionforU (x).
λ
Infact,eachradiationsourceinageneralRTproblemnotnecessarilycontributessignificantlytoU ateachpositionwithinthevolume
λ
considered for the RT calculation but often only within a fraction of it, which we call the “source influence volume”. In the ray-tracing
methoddescribedabove,ifoneknewinadvancetheextentoftheinfluencevolumeforeachemittingcell,itwouldbepossibletoreducethe
numberofcalculationsbysimplyperformingray-tracingonlywithinthisvolume.Unfortunately,itisnotpossibletoknowthisinformation
apriori,sincethefinalvaluesofU ateachpositionareavailableonlyafterperformingtheRTcalculation.However,onecanestimatein
λ
aconservativewaytheextentoftheinfluencevolumeforagivenemittingcell.Thebasicideaistoidentifyapointatsomedistancefrom
theemittingcellwhereitcanbeprovedthatthecontributionδU carriedbyarayisnegligiblecomparedtothelocalfinalvalueofU .If
λ λ
thisisthecase,thiscanimplythattheemittingcellwillnotcontributesignificantlytoU beyondthatposition,whichisalreadyoutsidethe
λ
influencevolumeoftheemittingcell[seeFig.1forsimpleexampleswherethiscriteriawillwork(left-handpanel)orcanfail(middleand
right-handpanels)ifnotproperlyapplied].
In practice, one can implement thismethod in the following way. First,one calculates alower limit U for U withinthe entire
λ,LL λ
volumebyperformingray-tracingfromeachemittingcellonlyuntilacertainarbitrarydistance(seeFig.2).OncethelowerlimitforU is
λ
calculated,onecanstarttheRTcalculationagainfromthebeginningbutthistimeisabletocheckifacontributionδU ,carriedbyarayto
λ
aparticularcell,isgoingtobesignificantornot.Infact,ifforacertainintersectedcellδU isnegligiblecomparedtothelocalvalueofthe
λ
lowerlimitU ,thenitwillbenegligiblealsocomparedtothefinalvalueofU atthatposition.Theactualimplementationofthismethod
λ,LL λ
consistsincheckingateachcellintersectionifδU <f ×U withf equaltoaverysmallnumbertobechosenappropriately(see§4
λ U λ,LL U
fordiscussiononthispoint).Whenthisconditionisrealizedandifthechosenvalueoff issufficientlysmall,thenthecontributionsδU
U λ
willalsobenegligibleforallthecellsbeyondthatpositioninthesameraydirection.
Thisimpliesthat,onceδU < f ×U atacertaindistancefromanemittingcellalongaray-path,onecanstoptheray-tracing
λ U λ,LL
calculationforthatparticularrayatthatposition.Infact,inthatcase,theintersectedcellisalreadyoutsidetheinfluencevolumeoftheemit-
tingcelloriginatingtheray(seeFig.3).Inthisway,thetotalamountofcalculationstobeperformedcanbesubstantiallyreduced,makingit
feasibletouseamodifiedversionoftheaboveray-tracingalgorithmtoinfertheU distributionwithincomplexdust/starsstructuresasthose
λ
observedingalaxies.Inthefollowingsubsection,weprovideasimplifieddescriptionoftheRTalgorithmwehavedevelopedbasedonthis
approach.
2.2 BasicDescriptionoftheDART-RayAlgorithm
Inthefollowingweprovideasimpledescriptionofthethreestepsperformedbythealgorithmwhichimplementsthesolutionstrategyout-
linedabovetocalculatetheRFEDdistributionU atasinglewavelengthλ.Thisdescriptionassumesthatagridofcellssubdividing the
λ
5 Thisapproach canberedundant because the rays originating from different cells canoften gothrough almostthesamepaths. Ontheother side, itis
convenientbecauseiteasilyallowstheimplementationoftheaccelerationtechniquespresentedinthiswork.
6 Nataleetal.
Figure1.Examplesofstarsanddustgeometrieswherethecriteriatoidentifythelimitofasourceinfluencevolumeworks(leftpanel)orcanfail(middle
andrightpanel).Left:TheraycomingfromthedistantsourcedoesnotcontributesignificantlytotheRFEDincell1,whichisilluminatedmainlybyother
sources.Asaconsequence,thesameraydoesnotcontributesignificantlyalsotothecellsbeyondcell1;Middle:Asintheleftpanel,therayfromthedistant
sourcedoesnotcontributesignificantlytotheRFEDincell1.However,inthisparticulargeometry,itscontributiontotheRFEDinthecellsbeyondcell1can
benotnegligiblebecausetheemissionfromtheothersourcesishighlyattenuatedbydust(dashedcells);Right:thedistancesourceispartofalargeemitter
distribution,whoseindividualRFEDcontributionstocell1areverysmallbutthecumulativecontributioncanbenotnegligible.
Figure2.FirstestimateoftheRFED.Theradialray-tracingisperformedonlyuntilalimitopticaldepthordistance.Dashedsquaresdenotecellcontaining
dust.
entiremodelhasalreadybeencreated.AcompletetechnicaldescriptionofthecodewillbegiveninSection§3.
FirstStep:CalculationofaLowerLimitforU
λ
Thefirststepconsistsofray-tracingfromeachemittingcelladoptingarayangulardensitysuchthatallthecellswithinacertainradiusor
acertainoptical depth fromeachemittingcellareintersectedbymultiplerays, asshown inFig.2(hereafter, wewillrefertotheprocess
ofintersectingallcellsinacertainregionwithmultipleraysas“fullysampling”thatregion).Thespecificvaluesforthelimitradiusand
opticaldepthcanbespecifiedintheinput.Theyshouldbelargeenoughtolettherayscoverarelativelylargefractionoftheentiremodelbut
smallenoughtoavoidatoolongcomputationaltimeinthisstep.TheinferredcontributionsδU tothevalueofU ineachcellcrossedby
λ λ
raysareaccumulatedintoanarrayU (LLstaysas“lowerlimit”).Assaidbefore,theinferredRFEDdistributionU representsonly
λ,LL λ,LL
alowerlimitforthefinalvalueU sincetherecouldbeotheremittingcells,atdistanceslargerthanthoseadoptedlimits,whosecontribu-
λ
tiontoU ateachparticularpointhasnotbeenconsideredyet.Inaddition,scatteredlightcontributionshavealsobeenneglectedatthispoint.
λ
SecondStep:ProcessingofSourceDirectLight
Inthesecondsteptheray-tracingprocedureisrepeatedagainfromthebeginningbutthistimetheraysfullysampletheregionsaroundeach
emittingcell,untiltheraycontributionδU foranintersectedcellbecomesmallerthanaverysmallfractionofthelowerlimitU ,that
λ λ,LL
isuntil δU < f ×U with f being avery small number (see §4 for more detailsabout how tochoose an appropriate value for
λ U λ,LL U
f ).Whenthisconditionisrealized,itmeansthatthepositionreachedbytherayisoutsidetheemittingcellinfluencevolume(seeFig.3)
U
DART-RAY 7
Figure3.Ray-tracing calculation fromonesourceuntilδUλ < fU ×Uλ,LL.Thepositionwhenthiscondition isrealized shouldbeoutsidethesource
influencevolume.Thus,onecanstoptheray-tracingcalculationatthatposition.
and the final value of U inthe crossed cell iscontributed mainlyby emitting cellsdifferent fromthe emittingcell originating theray6.
λ
Therefore,tothepurposeofcalculatingtheRFEDdistributionU ,thereisnoreasontoproceedwithafullsamplingofthecellsbeyondthe
λ
limitdeterminedinthisway.ApartfromstoringthevaluesoftheRFEDcontributionsδU ,aftereachraycrossesacell,thescatteredenergy
λ
informationarealsostoredforthatcell.Thatis,theluminosityscatteredwithinadiscretesetofsolidanglesisstoredforeachcellcontaining
dust.Thesevalueswillbeusedinthenextstep.Beforegoingtothethirdstep,thevalueofU isupdatedwiththecurrentestimationforU .
λ,LL λ
ThirdStep:ScatteringIterations
Afterthedirectlighthasbeenprocessedinthesecondstep,athirdstepisstartedwherethereisaseriesofiterationstoprocessthescat-
teredradiation.Eachcellwhichoriginatesscatteredlightistreatedasanemittingcellexactlyinthesamewayasinstep2.Ray-tracingis
performedfromeachdustycellbyfullysamplingallthevolumesurroundingthecelluntilthecontributionδU isnegligiblecomparedtoa
λ
smallfractionofU ,thelowerlimitforU updatedwiththenewvalueofU foundinstep2.Again,scatteringinformationarestoredas
λ,LL λ λ
wellaftereachraycrossingandthishigherorderscatteredradiationintensityisprocessedinsuccessiveiterations.Thesescatteringiterations
continueuntilthevastmajorityoftheluminosityofthesystemhasbeeneitherabsorbedbydustorhasescapedoutsidethemodel.Thatis,
until∆L <f L ,whereL isthetotalstellarluminosityemittedwithinthemodel,∆L istheamountofluminositystilltobeprocessed
λ L λ λ λ
(thatis,notabsorbedorescapedyet)andf isaparametertobesetintheinput.
L
Thedisadvantage ofthisprocedureisthatmanyofthecalculationsperformedinthefirststeparerepeatedonceagaininthesecond
step.Theadvantage isthatthenumber of calculationsavoidedcanbeveryhigh, thusreducingsignificantlythetotalcalculationtimefor
thosegeometrieswheretheinfluencevolumesoftheemittingcellsareonlyasmallfractionofthetotalvolume.Furthercharacteristicsofthe
code,notmentionedinthesimplifiedalgorithmabove,includetheadaptiveanddirectional-dependentangulardensityoftherays(see§3.2)
andthemethodsimplementedtocalculatethespecificintensityoftheradiationescapingoutsidethemodelinafinitesetofdirections(see
Section§3.4).
3 THERAY-TRACING3DCODE:TECHNICALDESCRIPTION
Inthissectionweprovideanextensivedescriptionofthecodewehavedeveloped.Thecodeconsistsoftwomainprogramsperformingthe
adaptive gridcreationand theRTcalculationrespectively. Inthefollowing subsections wedescribethe adaptivegrid creation(§3.1), the
6 Notethat,asshownintheright-handpanelofFig.1,therecouldbeemittingcellswhicharepartofalargeemittingcelldistribution,whoseindividualRFED
contributionsδUλmightbelowerthantheassumedenergydensitythresholdfUUλ,LL,butthecumulativecontributionisnotnegligible.Overlookingthe
presenceofthiskindofcellscanbeavoidedbychoosinganappropriatevaluefortheconstantfU.However,incasessuchasthatofanextendeddistribution
ofuniformvolumeemissivitywithinanopticallythinsystem,thevalueoffU tobechosen,inordertoreachanaccuratesolution,couldbeextremelylow.
Inthosecases,theremightbenosubstantialadvantageintermsofspeedbyusingthepresentedalgorithm,sincethemajorityofemittingcellscontribute
significantlytotheRFEDateachpositionintheentirevolume.
8 Nataleetal.
basicray-tracingroutine(§3.2)andeachofthethreestepsoftheRTalgorithmindetail(§3.3).Finally,wedescribethemethodsimplemented
toderivetheescapingradiationspecificintensity(§3.4).
3.1 AdaptiveGridCreation
Given an input spatial distribution of dust mass and stellar luminosity (either defined by analytical formulae or in a tabulated form), an
adaptivegridiscreatedinawaysuchthatthespatialresolutionishigherinregionswheretheradiationfieldintensityisexpectedtovaryin
amorerapidway,suchasthosewherethedensityofdustishigher.Aparentcellofcubicshape,enclosingtheentiremodel,issubdividedin
3x3x3=27childcubiccellsofequalsize7.Then,theaveragedustdensityandstellarvolumeemissivityarecalculatedwithinthevolumeof
eachnewlycreatedcell,togetherwithanestimationoftheirvariationwithineachcell.Furthercellsubdivisionproceedsforthosechildcells
whichdonotsatisfyuser-definedcriteriaandthosecellsbecomeparentsofevensmallerchildcells.Afterthat,theestimationofthecelldust
density/stellarvolumeemissivityandthecellsubdivisionprocedureareperformedforthenewcellsandsoon.Inthisway,atreeofcellsis
constructediterativelyuntilthechosencriteriaaresatisfiedforallthesmallestcellswithineachoriginalparentcell.Also,inordertoobtain
asmoothvariationofthegridresolution,furthercellsubdivisionisperformedtoavoiddifferencesincellsubdivisionlevelhigherthanone
betweenneighbourcells.Thecellswhichhavenotbeenfurthersubdivided,aftertheentiregridcreationhasbeencompleted,arecalledthe
leafcells.
ThecriteriaforcellsubdivisionshouldbechoseninordertoobtainbothnumericalaccuracyandanadequatecoverageoftheRFED
distributionwithinthemodel.Inordertoachieveagoodnumericalaccuracy,itisimportanttohaveleafcellswithsmalltotalopticaldepths
andwithsmallgradientsofdustdensityandstellarvolumeemissivitywithinthecells.Forexample,typicalconditionsforacelltobealeaf
cellcouldbe:
τ <<1 (4)
λ
∆ρ
<0.5 (5)
ρ
c
whereτ =k ρ l isthetotalcellopticaldepth,∆ρisanestimateofthedustdensityvariationwithinacellandρ istheaveragecelldust
λ λ c c c
density.Anequivalentcondition,astheoneexpressedbyEq.5,hastobefulfilledbythecellstellarvolumeemissivityj .Theseconditions
λ,c
allow us to usethe simple expression given by Eq. 3 tocalculate the variation of thespecific intensity withinacell toa good degree of
accuracy.However,notethatanadditionalconstraintisthemaximumallowednumberofsubdivisionlevelsNLVL MAX,whichhastobe
chosenintheinputaswell.Thevalueofthisparametershouldbesuchtoguaranteethattheinputcellparameterrequirements,suchasthose
above,arevalidforalloratleastthevastmajorityofcellsinthemodeland,atthesametime,avoidtocreatemodelswithtoomanycells.It
isdesirabletokeepthenumberofcellsaslowaspossible,sinceitisoneofthemainparametersaffectingthetotalcomputationaltime.After
creatingthegrid,theprogramcreatesatablewheretheusercanreadthemaximumandaveragecellparametersand,thus,quicklyverifyto
whichdegreetheinputcellrequirementshavebeenfulfilled.
AlltheinformationwhichdefinethegridisprintedonafilethatcanbereadfromtheRTprogram.Thisincludes:
–cellidnumber
–positionofcellcentreinCartesiancoordinates
–cellchildnumber:whichisequalto“-1”ifthecellisaleafcellortotheidnumberofthefirstchildcellcreatedduringthecellsubdivision
–cellindexnumber:binarycodeexpressingthepositionofthecellwithinthecelltree
–cellsizel
c
–celldustdensityρ
c
–cellstellarvolumeemissivityj
λ,c
3.2 ThebasicRay-TracingRoutine
In this subsection we describe the basic ray-tracing calculation for rays departing from an emitting cell and propagating throughout the
model. This is the core routine used by the RT algorithm described in the next subsection. Given an emitting cell, rays are casted into
multipledirectionsdefined by using theHEALPixsphere pixelation scheme (Gorski et al.2005, see alsoe.g. Abel &Wandelt 2002 and
Bisbaset al.2012 for other applications inRTcodes). Inthisscheme asphere isdivided in12 mainsectorsof equal size, whichcan be
furthersubdividedinsmallersphericalpixels(seeFig.4).ByassumingthatanHEALPixsphereiscentredontheemittingcellcentre,the
linesconnectingthecellcentretothecentreoftheHEALPixsphericalpixelsdefinedirectionsalongwhichonecanfollowthevariationof
thespecificintensitywithinthemodel.TheuseoftheHEALPixroutinesPIX2ANG NEST(translatingpixelindexnumbersintospherical
coordinates)anditsinverseANG2PIX NESTisaconvenientwaytohandlesetsofrayswithanadaptiveangulardensity(seebelow).
7 Weuseagridrefinementfactorequalto3,sinceithastheadvantageofpreservingthepositionofcellcentresaftereachcellsubdivisions.Fortechnical
reasons,thisfacilitatestheinclusioninthegridofcentralemittingpointsourcessuchthoseinthesolutionsdescribedinthenextsection.
DART-RAY 9
Figure4.Healpixsphereatdifferentangularresolution.Thesphericalpixelsintheupperleftspherearethosementionedas“HEALPixmainsectors”inthe
text.FigurefromGorskietal.(2005)(reproducedbypermissionoftheAAS).
Foreachray,thefollowingapproximationsareimplemented:
1) the calculation is performed as if all the cell luminosity propagates through solid angles defined by thepixels of an HEALPixsphere
centredonthecellcentre;
2)sinceonecanonlyfollowtheexactvariationofI alongthefinitesetofHEALPixdirections,itisassumedthattheevolutionofI along
λ λ
allthedirectionsincludedwithinthetotalsolidangleassociatedwitharay(determinedbytheadoptedHEALPixschemeangularresolution)
isexactlythesameasforthemaincentralraydirection(seeFig.5).
TheaboveapproximationsimplythatthetotalluminositydensityassociatedwitharayandflowingthroughanHEALPixsolidangle
ΩHP,EMatanydistancerfromthecentreoftheemittingcellisgivenby:
Lλ,ray(r)=Iλ(r)ΩHP,EMAEM (6)
whereAEMistheprojectedareaoftheemittingcell(assumedtobeequaltotheemittingcellsizesquared).
ThevariationofI ,whentheassociatedraycrossesacell,iscalculatedusingthefollowingexpressions.Thevariationduetothecrossed
λ
opticaldepthτ =k ρ ∆randcellvolumeemissivityj withinthecelloriginatingtherayisequalto:
λ λ c λ,c
j ∆r
I = λ,c (1−e−τλ) (7)
λ,0 τ
λ
ifτ >0(cellcontainingdust)8.Ifτ =0,weassumed:
λ λ
I =j ∆r (8)
λ,0 λ,c
which isconsistent withthe previous expression when τ → 0. During the scatteringiterations, the scattered light isconsidered for the
λ
calculationof thecell volumeemissivity, asitwillbeshown later.Instead, thenew valueof I afterthecrossing ofaraythrough acell
λ
differentfromtheoriginaloneissimplygivenby:
I =I e−τλ (9)
λ,i+1 λ,i
ApartfromcalculatinganewvalueforI ,whenaraycrossesacell,thecodedeterminesthecontributionoftheraytotheintersected
λ
cellRFEDandtheamountofenergyscatteredbythedustintheintersectedcell.However,thesecontributionscanbeaccuratelycalculated
onlyifasufficientlyhighrayangulardensityisused.Infact,dependingonthecellsizesandtheadoptedHEALPixangularresolution,at
anydistancefromtheemittingcell,theraybeamasawhole(notjustthemainraydirection)canintersecteitherasinglecellor agroup
of cellsat thesametime.Atdistances largeenough fromtheemittingcell,thebeam, originallyintersecting onlyone cellat atime,will
beginintersectingmorecellssimultaneously(seeFig.5).Assaidbefore,thecodetracestheevolutionofI onlyalongthemaindirectionof
λ
8 Actuallythenumericalimplementationrequirestoassumeasmallthresholdvaluehigherthanzero,suchthate−τλ 6=1whentheexponentialfunctionis
evaluated.
10 Nataleet al.
Figure5.RaybeamassociatedwithanHEALPixpixel,propagatingthroughoutthemodel.DuringtheRTcalculation,thevariationofIλisfollowedonly
alongthemaindirection(boldarrow).Thesamevariation isassumedfortheotherdirections withinthesameraybeam(dashedlines).Notethatatlarge
distancesfromemittingcell,thebeambeginsintersectingmorecellssimultaneously.
theraybeamanditisassumedthatforallthedirectionswithinthebeamtheI variationisexactlythesame.However,farawayfromthe
λ
emittingcell,thebeampropagatesthroughlargerphysicalvolumesandthepreviousapproximationcanbecomeveryinaccurate.Actually,in
ordertoobtainaprecisecalculationoftheRFEDcontributedbyanemittingcellCEMtoacertaincellCINT,itisdesirablethatseveralrays
originatingfromCEMareintersectingCINT.Inthisway,thecalculationofRFEDismoreaccuratebecauseofthebetterestimationofthe
averageI withintheintersectedcell.
λ
IfonedefinesAINT,theprojectedareaoftheintersectedcell(approximatedbythecellsizesquared), andΩINT = ArIN2T,thesolid
anglesubtendedbytheprojectedareaoftheintersectedcellandwithorigininthecentreoftheemittingcell(seeFig.6),inordertohave
moreraysfromCEMcrossingCINTonerequiresthat:
ΩHP,EM< NΩIrNayTs (10)
where Nrays is an input parameter equal to the minimum number of rays which should cross CINT. If this condition is fulfilled by all
theintersectedcellswithinacertainregionaround anemittingcell,wesaythattheraysare“fullysampling”thatregion,usingthesame
terminologyalreadyintroducedin§2.2.
Oncetheaboveconditionisfulfilledforanintersectedcell,thecodedefinesΩ=ΩHP,EMandthecontributionofaraytotheintersected
cellRFEDδU isgivenby
λ,INT
δU = <Iλ >AEMΩ∆cr (11)
λ,INT VINT
wherecisthespeedoflight, ∆cr isequaltothetimeneededbythelighttocrosstheintersectedcell,VINTisthevolumeoftheintersected
cellandtheaveragevalueofI alongthecrossingpathisequalto:
λ
I (1−e−τλ)
<I >= λ,i (12)
λ τ
λ
ifτ >0and
λ
<I >=I (13)
λ λ,i
ifτ =0.
λ
TheraycontributiontothespecificintensityscatteredbytheintersectedcellδI (θ,φ)is(seeFig.6):
λ,sca
δI (θ,φ)= Iλ,extwλ,scaAEMΩΦ(θ,φ) (14)
λ,sca AINTΩHP,INT
where ωλ = kkλλ,,escxat, ΩHP,INT is the solid angle determined by the HEALPixangular resolution adopted to store the scattered radiation
intensityintheintersectedcell,Φ (θ,φ)isatermrepresentingtheintegrationoftheHenyey-Greensteinscatteringphasefunctionoverthe
HG
solidangleΩHP,INT: