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