ELMAG: A Monte Carlo simulation of electromagnetic cascades on the extragalactic background light and in magnetic fields

0.16 MB·English

by  M. Kachelriess

ELMAG: A Monte Carlo simulation of electromagnetic cascades on the
extragalactic background light and in magnetic fields
M.Kachelrießa,S.Ostapchenkoa,b,andR.Toma`sc
aInstituttforfysikk,NTNU,Trondheim,Norway
bD.V.SkobeltsynInstituteofNuclearPhysics,MoscowStateUniversity,Russia
cII.Institutfu¨rTheoretischePhysik,Universita¨tHamburg,Germany
2
1
0
Abstract
2

n AMonteCarloprogramforthesimulationofelectromagneticcascadesinitiatedbyhigh-energyphotonsandelectrons
a interactingwithextragalacticbackgroundlight(EBL)ispresented. PairproductionandinverseComptonscattering
J
on EBL photonsas well as synchrotronlosses and deflections of the chargedcomponentin extragalactic magnetic

1
fields (EGMF) are includedin the simulation. Weighted sampling of the cascade developmentis applied to reduce
1
thenumberofsecondaryparticlesandtospeedupcomputations. Asfinalresult, thesimulationprocedureprovides

] the energy, the observationangle, and the time delay of secondarycascade particles at the presentepoch. Possible
E applicationsare the study of TeV blazars and the influence of the EGMF on their spectra or the calculation of the
H contributionfromultrahighenergycosmicraysordarkmattertothediffuseextragalacticgamma-raybackground.As
. anillustration,wepresentresultsfordeflectionsandtime-delaysrelevantforthederivationoflimitsontheEGMF.
h
p
Keywords: Electromagneticcascades,extragalacticbackgroundlight,extragalacticmagneticfields.
-
o
r
t PROGRAMSUMMARY
s
a Manuscript Title: ELMAG: A Monte Carlo simulation of electromagnetic cascades on the extragalactic background light and in
[ magneticfields

ProgramTitle:ELMAG 1.01
2
v JournalReference:
8 Catalogueidentifier:
0 Licensingprovisions:
5 Programminglanguage: Fortran95
5 Computer:AnycomputerwithFortran95compiler
6. Operatingsystem:AnysystemwithFortran95compiler
0 RAM:4Mbytes
1 Numberofprocessorsused:arbitraryusingtheMPIversion
1 Supplementarymaterial:seehttp://elmag.sourceforge.net/
: Keywords:Electromagneticcascades,extragalacticbackgroundlight,extragalacticmagneticfields
v
Classification:11.3CascadeandShowerSimulation,11.4QuantumElectrodynamics
i
X Natureofproblem:Calculationofsecondariesproducedbyelectromagneticcascadesontheextragalacticbackgroundlight(EBL)
r Solutionmethod: MonteCarlosimulationofpairproductionandinverseComptonscatteringonEBLphotons; twoparametrisa-
a tionsfromRef.[1]canbechosenasEBL;weightedsamplingofthecascadingsecondaries;recordingofenergy,observationangle
andtimedelayofsecondaryparticlesatthepresentepoch.
Restrictions:Deflectionsandtime-delaysarecalculatedinthesmall-angleapproximation.
Unusualfeatures:
Additionalcomments:
Runningtime: 400secondsfor103 photonsinjectedatredshiftz=0.2withenergyE =100TeVusingoneIntel(R)Core(TM)i7
CPUwith2.8GHz.
PreprintsubmittedtoComputerPhysicsCommunications January12,2012
References
[1] T.M.KneiskeandH.Dole,Astron.Astrophys.515(2010)A19[arXiv:1001.2132[astro-ph.CO]].
1. Introduction
TheUniverseisopaquetothepropagationofγ-rayswithenergiesintheTeVregionandabove[1]. Suchphotons
are absorbed by pair production on the extragalactic background light (EBL) [2–5], consisting mainly of infrared
lightandthecosmicmicrowavebackground(CMB).AsaresultthephotonfluxatenergiesE >10TeVfromdistant
∼
sourcesase.g.blazarsissignificantlyattenuatedonthewayfromthesourcetotheEarth. High-energyphotonsare
howevernotreallyabsorbedbutinitiateelectromagneticcascadesintheintergalacticspace,viathetwoprocesses
γ+γ →e++e− (1)
b
e±+γ →e±+γ. (2)
b
The cascade develops very fast until it reaches the pair creation threshold at1 s = 4E ε = 4m2 with ε as the
min γ γ e γ
characteristicenergyofthebackgroundphotonsγ . Electrons2 continuetoscatteronEBLphotonsintheThomson
b
regimewithaninteractionlengthofafewkpc,producingphotonswithaverageenergy
E = 4εγEe2 ≈3GeV Ee 2 (3)
γ 3 m2 1TeV
e
(cid:18) (cid:19)
usingε =2.7T ≈6.3×10−4eVasthetypicalenergyofCMBphotons.
γ CMB
The resulting shapeof the energyspectrumof the diffuse photonflux J can be estimated analytically[6] for a
γ
monochromaticbackgroundas
K(E /E )−3/2 for E ≤E ,
γ x γ x
J (E )= K(E /E )−2 for E ≤ E ≤ E ,
γ γ γ x x γ min

 0 for Eγ >Emin.
Here,E =m2/ε isthethresholdenergyforpair-production,whileE ≤ E istheenergyregionwherethenumber
min e γ γ x
of electrons remains constant. Since the last generation of e+e− pairs produced share the initial energy equally,
E = E /2,thistransitionenergyisgivenbyE = 4ε E2/(3m2) = E /3. Thusforamonochromaticbackground
e min x γ e e min
theplateauregioncharacterisedbyan1/E2spectrumextendsonlyoveronethirdofanenergydecade.
γ
AbetteranalyticaldescriptionofthecascadedevelopmentintheEBLusesadichromaticphotongas,withε =
CMB
6.3×10−4eVandε =1eVastypicalenergiesfortheCMBandthe(second)peakoftheIRbackground,respectively.
IR
Below one half of the thresholdenergyof pair productionon the IR, E ≈ E /2 = m2/(2ε ) ≈ 1.3×1011eV,
e min,IR e IR
thenumberofelectronsremainsconstant. Intheintermediateregime,E <E <E ,electronsareCompton
min,IR∼ ∼ min,CMB
scatteringonCMBphotonsintheThomsonregime,whilephotonsarestillproducinge+e−pairsonIRphotons.Thus
theenergyE belowwhichnoadditionalelectronsareinjectedinthecascadeisgivenby
x
4 ε E2 1 ε
E = CMB e = CMBE ≈50MeV. (4)
x 3 m2 3 ε min,IR
e IR
LetusnowcomparehowwellthisqualitativepictureagreeswiththecascadespectrumcalculatedwithourMonte
Carlo simulation. Figure 1 shows the results obtained with ELMAG for the diffuse spectrum of secondary photons
producedbyphotonsinjectedwithenergyE = 1014eVatredshiftz = 0.02andz = 0.15. Wenotefirstthattheslope
ofthephotonspectrumbelowE andtheobtainedvalueofE agreeverywellforbothdistanceswiththeprediction
x x
inthe simple “dichromaticmodel,”althoughthe latter assumesan infinitenumberof interactions. Incontrast, both
theextensionandtheshapeofthe plateauregionareless universal: Thesmallerthe distancetothe source,theless
pronouncedoccursthesteepeningofthephotonspectrumfromtheE−1.5 Thomsonslopetowardsthepredicted1/E2
2
1
E-1.5
10-1 E-1.9
)
E z=0.02
2J( 10-2
E
10-3 z=0.15
10-4
106107 108 109 10101011101210131014
E (eV)
Figure1: The(normalised)diffusephotonfluxE2J(E)fortwosourcesinjectingphotonswithenergyE =1014eVatredshiftz=0.02(—)and
z=0.15(—),respectively.
plateau. However,thedeviationoftheslopefromtheprediction,E−1.9 versusE−2,isminoralreadyfordistancesof
∼500Mpc.
Animportantapplicationofelectromagneticcascadesisthecalculationofvariouscontributionstotheextragalac-
tic diffuse gamma-raybackground(EGRB). Since the Universe acts as a calorimeter for electromagneticradiation,
accumulatingitin the MeV–TeVrange, themeasuredEGRB limitsallprocessesduringthe historyofthe Universe
thatinjectelectromagneticenergyabovethepaircreationthreshold.Examplesforsuchprocessesarephoto-pionand
p+γ → p+e++e−pair-productionofUHECRprotonsinteractingwiththecosmicmicrowavebackground(CMB)
CMB
[7],thedecayorannihilationof(superheavy)darkmatteroroftopologicaldefects[8]. Anotherimportantapplication
of electromagnetic cascades is the calculation of spectra from point sources as TeV blazars. If the spectra of such
sourcesextendtosufficientlyhighenergies,emittedphotonsinteractwiththeEBL.Thechargedcomponentofthese
cascadesisdeflectedbyextragalacticmagneticfields(EGMF),leadingpotentiallytohalosaroundpointsources[9–
12],todelayedechosofflaringemission[13]andinfluencestheobservedenergyspectrum[14]. Adetailedmodelling
of the electromagnetic cascade process is thus not only necessary to connect the observed energy spectra of TeV
sourceswiththeirintrinsicspectra,butprovidesalsoinformationaboutEGMFs.
TheextremelysmallinteractionlengthcomparedtotypicalsourcedistancesfromhundredsofMpctoGpcmeans
thatalargenumbernofinteractionstepshastobesimulatedusingaMonteCarloapproach.Theexponentialgrowth
of the number N = 2n of secondaries aggravates the computational load in a brute-force Monte Carlo approach.
The MonteCarlo programpresentedhere uses weightedsampling of the cascade developmentto reduce efficiently
the number of secondary particles which are traced explicitly. For maximally weighted sampling, the number of
secondaries stays on average constant as function of interaction steps. Synchrotron losses and deflections of the
chargedcomponentinextragalacticmagneticfields(EGMF)areincludedinthesimulationtoo.Theversionpresented
hereisrestrictedtothelimitofsmalldeflections.
2. Modellingofthecascadeprocess
2.1. Interactionrateofphotonsandelectrons
TheinteractionrateR (E,z)ofphotonswithenergyE atredshiftzcanbeconnectedtothepair-productioncross
γ
sectionσ (s)andthespectraldensityofbackgroundphotonsn (E,z)as
pair γ
1 ∞ 1
R (E,z) = dE′ n (E′,z) dµ(1−µ)σ (s)Θ(s−s )
γ γ pair min
2
Z0 Z−1
1Weusenaturalunits,~=c=kB=1,throughoutthetext.
2Wecallfromnowonelectronsandpositronscollectivelyelectrons.
3
1 smax(E) s
= ds sσ (s) I ,z , (5)
8E2 pair γ 4E
Zsmin (cid:18) (cid:19)
whereweintroducedtheauxiliaryfunction
EmaxdE′
I (E ,z)= n (E′,z). (6)
γ min E′2 γ
ZEmin
HerewehavealsoassumedthattheEBL,asanytrulydiffusebackground,isisotropic.
Thec.m.energysquaredin a γγ interactionis s = 2EE′(1−µ) withµ = cosϑ, whilethe integrationlimits are
givenbythepairproductionthresholds =4m2ands (E)=4EE withE ∼14eVasthehighenergycutoff
min e max max max
oftheEBLbackground.Thewell-knownpair-productioncrosssectionσ (s)isgivenby
pair
3 m2 1+β
σ (s)= σ e (3−β4)ln −2β(2−β2) , (7)
pair 4 Th s 1−β
" #
withσ =8πα2/(3m2)asThomsoncrosssectionandβ= 1−4m2/s.
Th e e
Electronsemitin theThomsonregimemainlysoftphotons,cf.Eq.(3). Tospeedupthe simulation,we include
p
therefore as discrete interactions only those which produce secondary photons above an arbitrary energy threshold
E . Theremainingsoftinteractionsareintegratedoutandincludedascontinuousenergyloss. Thuswedefinethe
thr
interactionrateR (E,z)ofanelectronwithenergyEatredshiftzas
e
1 ∞ 1
R (E,z) = dE′n (E′,z) dµ(1−βµ)σ (s,ε)Θ(s−s (ε))
e γ C min
2
Z0 Z−1
1 smax(E) s−m2
= ds(s−m2)σ (s,ε)I e ,z (8)
8βE2 e C γ 2E(1+β)
Zsmin(ε) !
withε= E /E,β= 1−m2/E2,ands=m2+2EE′(1−βµ).Theintegrationlimitsaregivenbys (ε)=m2/(1−ε)
thr e e min e
ands (E)=m2+2EE (1+β),whiletheComptonscatteringcrosssectionσ (s,ε)integratedabovethethreshold
max e p max C
εisgivenby
3 y −y ln(y /y ) 4y (1+y )
σ (s,ε) = σ y max min max min 1− min min
C 4 Th min 1−y y −y (1−y )2
min " max min min !
4(y /y +y ) y +y
+ min max min + max min . (9)
(1−y )2 2
min #
Here,y = m2/sandy = 1−εarerespectivelytheminimalandthemaximalenergyfractionsofthesecondary
min e max
electron.
In turn, for the electron energy loss per unit distance due to the emission of photons of energies E < E one
thr
obtains
dE 3 1−y ln(1/y ) 4y (1+2y )
ICS/thr(s,ε) = σ y max max −1 1− min min
dx 4 Th min 1−y 1−y (1−y )2
min " max ! min !
1 2y (1+2y /y )(1−y )
+ (1−y )(1+2y )+ min min max max . (10)
6 max max (1−y )2
min #
IntheleftpanelofFig.2weshowtheinteractionratesR ofelectronsandphotonsatthepresentepochasfunction
i
of energy. The difference between the “best-fit” and the “lower-limit” EBL from Ref. [5] becomes visible in the
interactionrate R ofphotonsonlyin theenergyrangebelow1014eV.TheinteractionrateR ofelectronsisshown
γ e
only for the “best-fit” EBL but for two different values of the threshold E used in the Compton scattering cross
thr
section, E = 3×104eV and E = 3×106eV. Note that while E = 3×106eV leads already below 1011eV
thr thr thr
to strong deviationsfrom the Thomson scattering cross section, the resulting photonspectrum is influenced by the
thresholdmainlyatenergiesintheMeVrangeandbelow,cf.Eq.(3).
4
103 e, Eth=3×104 eV 1) -11 Fermi
-s
102 2
) -m
-1pc 10 e, Eth=3×106 eV g c -12 HESS
M r
e
R ( 1 F/
2 f=0.10
10-1 (E0-13 ff==00..5700
10-2 g, EBL lower limit g1 f=0.80
g, EBL best fit o f=0.90
10-3 l
-14
109 1011 1013 1015 1017 1019 8 9 10 11 12 13 14
log (E/eV) log (E/eV)
10 10
Figure2: Left:InteractionrateRatz=0asfunctionoftheenergyEforelectronswithEth =3×104eVandEth =3×106eVandforphotons
withthe“best-fit”andthe“lower-limit”EBLfromRef.[5].Right:Fluencecontainedinsidethe95%confidencecontourofthePSFofFermi-LAT
asfunctionofenergyforEGMFwithtop-hatprofileandfillingfactor f varyingfrom f =0.1to f =0.9withEmax=20TeV.
2.2. Interactions
Themodellingofγγandeγinteractionsstartsfromsamplingthec.m. energysquaredsofthecollisionaccording
to the integrandsof Eqs. (5) and (8), respectively. Technically, the rejection method is used as in most other cases
to choose s according to its probability distribution: The value of s is first sampled logarithmically in the interval
[s ,s ],thenthechoiceisacceptedwiththeprobabilityproportionaltostimestheintegrandofEq.(5)and(8),or
min max
otherwiserejected.
Forgivens,theenergyfractionyofthelowestenergysecondarylepton(electronorpositron)inthepairproduction
processissampledaccordingtothecorrespondingdifferentialcrosssection
dσ (s,y) 1 y2 1−β2 (1−β2)2
pair ∝ +1−y+ − 1+2β2(1−β2) , (11)
dy y 1−y 1−y 4y(1−y)2
" #(cid:30)h i
withβ= 1−4m2/s. Theothersecondaryleptonhasthentheenergyfraction1−y. Similarly,theenergyfractiony
e
ofthesecondaryelectronintheinverseComptonprocessissampledaccordingtothedifferentialcrosssection
p
dσ (s,y) 1 1+y2 2y (y−y )(1−y)
ICS ∝ − min min , (12)
dy y 2 y(1−y )2
" min #
withy =m2/s.
min e
2.3. Stackingandweightedsampling
Theproducedsecondaryparticlesarethen subjectto a weightedsamplingprocedure: A secondaryparticle car-
ryingthefractionyoftheparentenergyisdiscardedwiththeprobability(1−yαsample),oraddedwiththeprobability
yαsample tothe stack. Dependingonthechoiceofthesamplingparameter(0 ≤ αsample ≤ 1)eitherallthe secondaries
arekeptinthecascade(α = 0)oronlysomerepresentativeonesareretained. Inparticular,onesecondaryper
sample
interaction is retained on average for the default value α = 1. As compensation, each particle in the cascade
sample
acquiresaweightwwhichisaugmentedaftereachinteractionasw→w/yαsample. Theparticlesinthestackareordered
according to their energies. After the interaction, the lowest energy particle is extracted from the stack and traced
furtherinthecascadeprocess.
Theoptimalvalueofα dependsonthetypicalenergyoftheinjectedphotons. Ifthelatterissolowthatthe
sample
cascades consists on average of only few steps, reducingα may be advantageousbecause the fluctuations are
sample
therebyreduced.
5
2.4. Synchrotronlosses
Synchrotronenergylossesofelectronsareaccountedforinthe continuousenergylossapproximationusingthe
interpolationformula[23],
dE m2χ2
≈ e , (13)
dx [1+4.8(1+χ) ln(1+1.7χ)+3.44χ2]2/3
withχ=(p /m )(B/B ),wherep denotesthemomentumperpendiculartothemagneticfieldandB =4.14×1013G
⊥ e cr ⊥ cr
thecriticalmagneticfield.
2.5. Angulardeflectionandtimedelay
For the energies considered, E > MeV and E >10GeV, secondaryparticles are emitted in the forwarddirec-
γ∼ e∼
tion. Thustheangulardeflectionofthecascadeparticlesresultsfromthedeflectionsofelectronsintheextragalactic
magneticfield(EGMF).Ifthecoherencescale oftheEGMFismuchlargerthanelectronmeanfreepathλ = R−1,
e e
anelementarydeflectionangleofi-thelectroninthecascadechaincanbecalculatedassumingthefieldtoberegular
overthedistanced travelledbytheelectron,
i
p −1 d B
β ≃0.52◦ ⊥ i , (14)
i TeV 10kpc 10−15G
(cid:18) (cid:19) !(cid:18) (cid:19)
with p beingthemomentumcomponentperpendiculartothelocaldirectionofthemagneticfield. Insidethepatch j
⊥
ofachosencoherencelength,thedeflectionanglesβ perelectronpatharesummedupcoherently,β = β. Thede-
i j i i
flectionsanglesβ percoherentmagneticfieldpatcharethensummedquadraticallyintherandom-walkapproximation
i
P
toobtainastotaldeflection,i.e.theangleβbetweentheinitialandfinalphotonsinthecascade,
β= β2, (15)
i
s i
X
Inthesmall-angleapproximationandassumingsphericalsymmetry,theangleβbetweentheinitialandfinalphotons
inthecascadeisrelatedtotheemissionangleαandtheobservationangleϑasα=β−ϑ,cf.Ref.[11].
As the energy of the cascade particles quickly degrades along the cascade chain, the largest contribution to β
comesfromthelastelectroninthechain. Thisallowsustoapproximatethecorrespondinggeometrybyatriangular
configurationandtoobtainasrelationbetweenβandϑ
x
sinϑ= sinβ. (16)
L
Here xreferstothedistancefromthesourceStothepointPwherethefinalphotoninthecascadebranchhasbeen
createdandListhetotaldistancebetweenthesourceandtheobserverO.Forsmallϑ,wethushave
x
ϑ= sinβ. (17)
L
Thetimedelay∆t ofphotonswithrespecttothestraightlinepropagationfromthesourceisthen
geo
∆t ≃ x(1+sinα/sinϑ)−L≃2x(1−x/L) sin2β. (18)
geo
We add to this geometricaltime delay ∆t the kinematicaltime delay ∆t due to velocity v < c of the electron,
geo kin
althoughthelatterisusuallynegligible.
2.6. Cosmology
The connectionbetween redshiftz, comovingdistance r and light-traveltime t calculated for a flat Friedmann-
Robertson-WalkeruniversewithΩ =0.7andΩ =0.3iscontainedinthefileredshift.
Λ m
6
3. Programmestructure
Theprogrammeisdistributedamongthefilesmodules101.f90,user101.f90,init101.f90,elmag101.f90
andaux101.f90.Thefilemodules101.f90containsthedefinitionofinternalvariables,mathematicalandphysical
constants; for standard applications of the programme no changes by the user are needed. The file user101.f90
contains the input/output subroutines developed by the user for the desired task. An example file is discussed in
Sec. 5. Data files of the used EBL backgrounds and the cosmological evolution of the universe are provided in
the directory Tables. They are read by the subroutines init EBL(myid), init arrays(myid) and the func-
tion aintIR(E,z)inside the file init101.f90. Then the function w EBL density tab(emin,zz)tabulates the
weightedbackgroundphotondensityI definedinEq.(6),followedbythetabulationoftheinteractionrateinthethe
γ
subroutinerate EBL tab(e0,zz,icq)andoftheelectronenergylossesduetotheemissionofphotonswithenergy
belowthethresholdineloss thr tab(e0,zz,icq).
Wediscussnowinmoredetailthesubroutinesandfunctionsofthefileelmag101.f90whichconstitutethecore
oftheprogramme:
• subroutine cascade(icq,e00,weight0,z in)
Follows the evolution of the cascade initiated by a photon (icq = 0) or an electron/positron (icq = ±1)
injectedatredshiftz inwithenergye00andweightweight0untilallsecondaryparticleshaveenergiesbelow
theenergythresholdethrorreachedtheobserveratz=0.
• subroutine angle delay(the2,xx,rcmb,theta,dt)
Determinesthephotontime-delaydtandtheobservationanglethetafromthermscascadedeflectionangle
the2andthephotonemissionpointxxbytheparentelectron/positron.
• subroutine interaction(e0,x0,zz,t,weight,the1,the2,xxc,xx,dt,icq)
Handlesoneinteractionwithbackgroundphotons: determinesthec.m.energysgamofthereactionviaacall
to sample photon or sample electron(e0,zz,sgam,ierr), the energy fraction z of secondaries via a
calltothefunctionszpair(sgam)orzics(e0,sgam),andstoresthenthesecondariescallingthesubroutine
store particle.
• subroutine store particle(e0,x0,zz,t,ze,weight,the1,the2,xxc,xx,dt,icq)
Decidesifa producedsecondaryis storedusingweightedsampling;if yes, it addsthesecondarytothe array
eventandre-ordersthearrayaccordingtotheparticleenergies.
• subroutine get particle(e0,x0,zz,t,weight,the1,the2,xxc,xx,dt,icq)
Readsthe secondarywiththe lowestenergyoutofthe arrayeventandreducesthe particlecounterjcmbby
one.
• subroutine sample photon(e0,zz,sgam,ierr)andsample electron(e0,zz,sgam,ierr)
Determinesthecmsenergysgamofaninteractionatredshiftzz.
• double precision function w EBL density(emin,zz)
DeterminestheweightedbackgroundphotondensityI definedinEq.(6).
γ
• function int point(e0,x0,zz,icq)
FindsthenextinteractionpointforinteractionwithEBLphotons.
• function sigpair(sgam)andsigics(e0,sgam)
CalculatethepairproductionandinverseComptoncrosssection,respectively.
• function zpair(sgam)andzics(e0,sgam)
DeterminetheenergydistributioninpairproductionandinverseComptonscattering,respectively.
• function zsigics(e0,sgam)
Calculates the electronenergylosses per unitdistance due to photonemission belowthe chosenthresholdin
Comptonscattering.
7
1) -11 Fermi 1) -11 Fermi
-s t < 10 yr -s
2 2
-g cm -12 103 < t(yr) < 107 HESS -g cm -12 t < 107 yr HESS
er er t < 10 yr
F/ F/ t < 102 yr
2 2
(E0-13 (E0-13 tt << 110034 yyrr
g1 g1 t < 105 yr
lo 10 < t(yr) < 103 lo t < 106 yr
-14 -14
7 8 9 10 11 12 13 14 7 8 9 10 11 12 13 14
log (E/eV) log (E/eV)
10 10
Figure3: Fluencecontainedinsidethe95%confidencecontourofthePSFofFermi-LATasfunctionofthetime-delayforB=10−17G;leftfor
individualtimelayers,rightcumulativetimes.
• function zloss(e0,zz)
Interpolatestheintegratedenergylossduetoemissionofphotonsbelowthethreshold.
• subroutine rate EBL(e0,zz,icq)
InterpolatestheinteractionratesR onEBLphotons.
i
• function eloss syn(E,begmf)
CalculatesthesynchrotronlossesaccordingtoEq.(13).
• function themf(e0,dx,begmf)
DeterminesthedeflectionangleintheEGMF.
Thefileaux101.f90containsauxiliaryfunctions,e.g.therandomnumbergeneratorpsranfromRef.[24].
4. Exampleinputandoutput
Thefileuser101.f90isanexamplefilefortheinput/outputsubroutineswhichshouldbedevelopedbytheuser
forthedesiredtask. Wediscussnowtheexamplecontainedinthedistribution.
4.1. Exampleinput
The input variables specified in the module user variablesare: the choice of the EBL model (model), the
numberofinjectedparticles(nmax),thejetopeningangleofthesourceindegrees(th jet),thesamplingparameter
a smp, the energy threshold ethr for Compton scattering, and the maximal photon energy egmax. The last two
parameters serve also as minimal and maximal energy in the energy spectra produced as output. In subroutine
user main(myid,nmax)theinitialredshiftzandtheparticletypeicqoftheinjectedparticlesisfixed.
z = 0.14d0 ! initial redshift
do nl=1,nmax
call initial_particle(e0,weight) ! generate initial energy
icq = 0 ! (0 - gamma, +-1 - e+-)
call cascade(icq,e0,weight,z) ! starts e/m cascade
enddo
Thesubroutineinitial particle(e0,weight)choosestheenergyandtheweightofoneinitialparticleinthe
energyrange[emin,egmax]accordingtoabrokenpower-lawwithexponentsgam1belowebreak,andgam2above.
The magnetic field B is modeled as patches of uniform field-strength |B| of size l . The value of the coher-
coh
encelengthisfixedbytheparametercohlnthinmodule user variables,thefield-strengthperpendiculartothe
propagationdirectioninthefunction bemf(r).
8
4.2. Exampleoutput
The energy e0, the observationangle theta and the time delay dt of secondarycascade particles with weight
weightoftypeicqreachingtheobserveratz=0arerecordedbythesubroutineregister(e0,theta,dt,weight,icq)
andbinnedin variousdata arraysdefinedin the module user result. All dataarraysexistsin two versions, e.g.
spec(n bin,0:1) and spec tot(n bin,0:1). Using MPI [15], the former arrays contain the result of a single
process,whicharesummedbycall MPI REDUCEintospec tot(n bin,0:1),
n_array = 2*n_bin
call MPI_REDUCE(spec,spec_tot,n_array,MPI_DOUBLE_PRECISION,MPI_SUM,0, & \\
MPI_COMM_WORLD,ierr) ! sum individal arrays spec
Finally,thesubroutine user output(n max,n proc)writesthedataarrayswiththeresultsinthefilescontained
inthesubdirectoryData.
Thefilespec diffcontainsthenormalised(diffuse)energyspectraofphotonsandelectrons,intheformatE/eV,
E2dN /dE andE2dN /dE. Thefilespec 95includestheenergyspectraofphotonsinsideandoutsidethe95%area
γ e
ofthepoint-spreadfunctionofFermi-LAT,in theformatE/eV, E2dN /dE(ϑ < ϑ )and E2dN /dE(ϑ < ϑ ). An
γ 95 γ 95
approximationtothepoint-spreadfunctionofFermi-LATisdefinedinthefunction thereg en(en).
The rightpanelof Fig. 2 shows the energyspectra of photonsarrivingwithin the 95% area of the point-spread
functionof Fermi-LATfordifferentfillingfactorsoftheEGMFwhichcanbechosenbythe parameterfracin the
function bemf. Otherwise the default values contained in the distributed file user101.f90are used. The photon
fluenceiscomparedtoH.E.S.S.data[16]andupperlimitsfromFermi-LAT[17]fortheTeVblazar1ES0229+200.
The files spec 95 t and spec 95 c contain the energyspectra of photons arriving within the 95% area of the
point-spread function, with the time-delay binned in seven time intervals, t < 10yr, 10yr < 102yr,...t > 106yr.
The file spec 95 c is the cumulative version of the distribution in spec 95 t. Figure 3 shows the fluence inside
the 95% PSF of Fermi-LAT for an injection spectrum dN /dE ∝ E−2/3 at redshift z = 0.14 with maximal energy
γ
E =20TeV.
max
5. Possibleextensions
We discuss four possible applications of the simulation ELMAG and the required extensions to perform them,
orderedbythecomplexityofthenecessarychangesandadditions.
EGRBfromdarkmatterdecaysorannihilations. High-energyelectronsandphotonscanbegeneratedbydecaysor
annihilationsofsufficientlyheavydarkmatterparticles. Forinstance,theannihilationmodeXX → e+e− ofthedark
matterparticle X withmassm wouldcorrespondtotheinjectionoftwoelectronswithenergyE = m . Theonly
X e X
necessaryadditionfor the calculationof the resulting EGRB is a subroutinechoosingthe injectionpointaccording
to the so-called boost factor B(z) which accounts for the redshift dependent clustering of dark matter in galaxies.
Additionally,the desiredfragmentationfunctionsdN /dE of the X particlesshouldbe includedintothe subroutine
γ
initial particleinthecaseofphotonsfromhadronicdecayorannihilationmodes.
EGRBfromUHECRsandcosmogenicneutrinos. TheGreisen-Zatsepin-Kuzmincutoffisasteepeningoftheproton
spectrum at the energy E ≈ (4 − 5) × 1019 eV, caused by photo-pion production on the CMB. An additional
GZK
signatureforthepresenceofextragalacticprotonsinthecosmicrayfluxandtheirinteractionwithCMBphotonsisthe
existenceofultrahighenergycosmogenicneutrinosproducedbychargedpiondecays[18],whilethecorresponding
fluxofcosmogenicneutrinosfromultrahighenergynucleiissuppressed. Photonsandelectronsfrompiondecayand
p+ γ → p +e+ + e− pair-production lead to a contribution to the EGRB which can be used to limit cosmic
CMB
rays(CR)modelsandfluxesofcosmogenicneutrinos. To performthistask, ELMAGhastobecoupledtoa program
performingthe propagationof ultrahigh energy cosmic rays which providessecondary electrons and photonsfrom
CRinteractionsasinput.Asthecommunicationbetweenthetwoprogrampartsisrestrictedtocallsofthesubroutine
cascade(icq,e0,weight,z),such a combination should be straightforward. For an example where ELMAG was
usedinthiscontextseeRef.[19].
9
Extensionto3-dimensionalcascades. Goingbeyondthesmall-angleapproximationrequiresto calculatetheactual
trajectoryofelectronssolvingtheLorentzequation.Additionally,scalarquantitieslikex,xxc, e0, begmf,...have
to be changed into three-dimensional vectors and the type one event in the module stack has to be adjusted.
Finally,theimageofathree-dimensionalcascadecanbecalculatedusingthemethoddescribedinRef.[20]. Foran
illustrationofpossibleapplicationsofELMAGtothisproblemseeRef.[21].
Treatmentofinteractionsinsources. Photonsandelectronsaregeneratedofteninsourcescontainingdensephoton
fields, as e.g. near the cores of active galactic nuclei. In this case, electromagnetic cascades take place on non-
thermal,anisotropicphotonbackgroundsinsidethesourcebeforetheescapingparticlescascadeontheEBL.Inorder
todescribebothcascadesontheEBLandinsidethesource,subroutinesase.g.init EBLhavetobedoubled,addinga
correspondingsubroutineinit sourceforthephotonfieldinsidethesource.Allexistingsubroutineswhichdepend
onthechosenEBL(i.e.use EBL fit)havetoadapted. Inparticular,therejectionmechanisminsubroutinesase.g.
sample photonhastoadjusted. Asaresult,suchanextensionrequiresconsiderableworkandthoroughtestsofthe
changedcode.AshortdiscussionofanisotropicphotonfieldsisgiveninRef.[22].
6. Summary
We presented a Monte Carlo program for the simulation of electromagnetic cascades initiated by high-energy
photons and electrons interacting with the extragalactic background light. The program uses weighted sampling
of the cascade developmentand treats Thomson scattering below a chosen threshold in the continuousenergy loss
approximationinordertospeedupcomputations.
Possible applicationsare the studyof TeV blazarsand the influenceof the EGMFon theirspectra orthe calcu-
lation of the contribution from ultrahigh energy cosmic rays or dark matter to the diffuse extragalactic gamma-ray
background.Asanillustrationforpossibleapplicationswepresentedresultsfordeflectionsandtime-delaysrelevant
for the derivation of limits on the EGMF studying the spectra of TeV blazars. Other possible applications include
e.g.thecalculationofthecontributionfromultrahighenergycosmicraysordarkmatterannihilationstothediffuse
extragalacticgamma-raybackground.
Acknowledgements
We aregratefultoVenyaBerezinskyforvaluablediscussionsandtoAndrewTaylorforcross-checkingsomeof
ourresults. ThisworkwaspartiallysupportedbytheprogramRomforskningoftheNorwegianResearchCouncil.
References
[1] R.J.GouldandG.P.Schreder,Phys.Rev.155,1404(1967);ibid.,1408(1967);F.W.Stecker,Astrophys.SpaceSci.6,377-389(1970).
[2] F.W.Stecker,M.A.MalkanandS.T.Scully,Astrophys.J.648,774(2006)[arXiv:astro-ph/0510449].
[3] J.R.Primack,R.C.GilmoreandR.S.Somerville,AIPConf.Proc.1085,71(2009)[arXiv:0811.3230[astro-ph]].
[4] A.Franceschini,G.RodighieroandM.Vaccari,Astron.Astrophys.487,837(2008)[arXiv:0805.1841[astro-ph]].
[5] T.M.KneiskeandH.Dole,Astron.Astrophys.515,A19(2010)[arXiv:1001.2132[astro-ph.CO]].
[6] V.S.BerezinskyandA.Yu.Smirnov,Astrophys.Sp.Sci.32461(1975); seealsoV.S.Berezinskyetal.,AstrophysicsofCosmicRays
(Elsevier,Amsterdam(1990))andV.S.Berezinsky,Nucl.Phys.B380,478(1992).
[7] K.Greisen,Phys.Rev.Lett.16,748(1966);G.T.ZatsepinandV.A.Kuzmin,JETPLett.4,78(1966)[PismaZh.Eksp.Teor.Fiz.4,114
(1966)];V.S.BerezinskyandS.I.Grigor’eva,Astron.Astrophys.199,1(1988).
[8] ForareviewseeM.Kachelrieß,ComptesRendusPhysique5,441(2004)[arXiv:hep-ph/0406174].
[9] F.A.Aharonian,P.S.CoppiandH.J.Vo¨lk,Astrophys.J.423,L5(1994)[arXiv:astro-ph/9312045].
[10] A.NeronovandD.V.Semikoz,JETPLett.85,473(2007)[arXiv:astro-ph/0604607].
[11] K.Dolag,M.Kachelrieß,S.OstapchenkoandR.Toma`s,Astrophys.J.703,1078(2009)[arXiv:0903.2842[astro-ph.HE]].
[12] A.Elyiv,A.NeronovandD.V.Semikoz,Phys.Rev.D80,023010(2009)[arXiv:0903.3649[astro-ph.CO]].
[13] R.Plaga,1995,Nature,374,430;K.Muraseetal.,Astrophys.J.686,L67(2008).
[14] P.d’Avezac,G.DubusandB.Giebels,Astron.Astrophys.469,857(2007)[arXiv:0704.3910[astro-ph]];A.NeronovandI.Vovk,Science
328,73(2010)[arXiv:1006.3504[astro-ph.HE]];F.Tavecchioetal.,Mon.Not.Roy.Astron.Soc.406,L70(2010).
[15] For information on MPI see e.g. “Message Passing Interface Forum. MPI: A Message Passing Interface Standard, June 1995” on
http://www.mpi-forum.org.
[16] F.Aharonianetal.[H.E.S.S.Collaboration],Astron.Astrophys.475,L9(2007)[arXiv:0709.4584[astro-ph]].
[17] A.A.Abdoetal.[FermiLATCollaboration],Astrophys.J.707,707(2009)[arXiv:0910.4881[astro-ph.HE]].
10