Table Of Content

MNRAS000,1–18(4January2017) Preprint4January2017 CompiledusingMNRASLATEXstylefilev3.0 Physics of cosmological cascades and observable properties T. Fitoussi1,2(cid:63), R. Belmont1,2†, J. Malzac1,2, A. Marcowith3, J. Cohen-Tanugi3, and P. Jean1,2 1 Université de Toulouse; UPS-OMP; IRAP; Toulouse, France 2 CNRS; IRAP; 9 Av. colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France 3 Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, Montpellier, France 7 1 AcceptedXXX.ReceivedXXX.;inoriginalform4January2017 0 2 ABSTRACT n TeV photons from extragalactic sources are absorbed in the intergalactic medium a J and initiate electromagnetic cascades. These cascades offer a unique tool to probe the properties of the universe at cosmological scales. We present a new Monte Carlo 3 codededicatedtothephysicsofsuchcascades.Thiscodehasbeentestedagainstboth ] publishedresultsandanalyticalapproximations,andismadepubliclyavailable.Using E thisnumericaltool,weinvestigatethemaincascadeproperties(spectrum,haloexten- H sion, time delays), and study in detail their dependence on the physical parameters . (extra-galactic magnetic field, extra-galactic background light, source redshift, source h spectrum and beaming emission). The limitations of analytical solutions are empha- p sised. In particular, analytical approximations account only for the first generation of - o photons and higher branches of the cascade tree are neglected. r t Key words: gamma-rays: general, astroparticles physics, extra-galactic magnetic s a field, scattering, relativistic processes [ 1 v 4 1 INTRODUCTION Observation (or non-detection) of electromagnetic cas- 5 cades is crucial to several astrophysical issues. It offers a 6 The Universe is opaque to gamma-rays. Very high-energy unique tool to probe the intergalactic medium, especially 0 photonsfromextragalacticsourcesareabsorbedbytheam- theExtragalacticBackgroundLight(EBL)andtheEGMF. 0 bient soft radiation and converted into electron-positron The background photons involved in the cascades have two . pairs(Gould&Schréder1967;Wdowczyketal.1972).These 1 distinct origins. Inverse Compton scattering mainly occurs leptons are deflected by the ExtraGalactic Magnetic Field 0 onphotonsfromtheCosmicMicrowaveBackground(CMB) 7 (EGMF)andcoolthroughinverseComptonscattering,pro- while high energy gamma-rays are mostly absorbed by the 1 ducing new gamma-rays that may, in turn, be absorbed. EBL of stars and dust, which extends from infrared to ul- : The observable properties of the resulting electromagnetic v traviolet. Our knowledge of the EBL is limited. Direct ob- cascade depend on the characteristics of the intergalactic i servations at these wavelengths are very inaccurate due to X medium. The development of cascades has three main ob- strong contamination from the zodiacal light. The predic- servableeffects.First,thesourcespectrumisalteredbecause r tions of the different models proposed in the literature can a each high energy TeV photon is reprocessed into thousands differ by up to an order of magnitude, depending on wave- of GeV photons (Protheroe 1986; Roscherr & Coppi 1998; length and redshift (Franceschini et al. 2008; Dom´ınguez Aharonian et al. 2002; Neronov & Vovk 2010). Second, due etal.2011;Finkeetal.2010;Kneiske&Dole2010;Gilmore tothedeflectionofleptonsbytheEGMF,newgamma-rays et al. 2012, see Fig. 1). Absorption of the gamma ray spec- are emitted along different lines of sight, so that a point trumofhigh-energysourcesprovidesunequalconstrainson sourcemayappearasextended(Aharonianetal.1994;Eu- the EBL (Stecker et al. 1992). ngwanichayapant & Aharonian 2009). Third, as leptons are deflected,cascadephotonstravelalongerdistanceandarrive Recently,cosmologicalcascadeswerealsousedtoprobe with a significant time delay, as compared to unabsorbed, thepropertiesoftheextragalacticmagneticfield,theorigin primary photons (Kronberg 1995; Plaga 1995; Ichiki et al. of which is still debated (Durrer & Neronov 2013). A pri- 2008; Murase et al. 2008; Takahashi et al. 2008). mordialmagneticfieldcouldhavebeengeneratedduringin- flationorduringthephasetransitionwhenelectroweakand (cid:63) E-mail:thomas.fi[email protected] QCD forces decoupled. This field would have remained un- † E-mail:[email protected] affected during the evolution of the extragalactic medium. (cid:13)c 4January2017TheAuthors 2 T. Fitoussi et al. Alternatively magnetic fields generated by galaxies during thatispubliclyavailable1.Thiscodeisdedicatedtocascades largescalestructureformationcouldhavepropagatedinthe induced by high-energy photons (or leptons) and does not intergalacticmediumthroughplasmajets.Dependingonthe take into account hadronic processes (see e.g. Oikonomou properties of the field generation and evolution, its value B et al. 2014; Essey et al. 2011, for results on cosmic-ray in- is expected to lie in the range 10−17 to 10−9 Gauss (Essey duced cascades). It computes the physics of leptonic cas- et al. 2011; Finke et al. 2015), with coherence length λ cades at the highest level of precision and with the fewest B (scale of de-correlation of two nearby field lines) between approximations. Using this code, we present a systematic 10−6 and 104 Mpc. Electromagnetic cascades represent a exploration of the parameter space. unique tool to probe the intergalactic magnetic field (Aha- In Section 2 we present the basic analytical theory of ronian & Atoyan 1985) when conventional methods such as cosmological cascades and simple analytical estimates of Faraday rotation cannot be applied. Neronov & Semikoz their observables. The results of our code are presented in (2007) and Elyiv et al. (2009) have suggested to measure Section 3 and are tested against analytical approximations the extension of pair halos to probe the EGMF. Indeed and other published numerical results. The last sections of Neronovetal.(2010)demonstratedthatforastrongenough the paper is devoted to an exploration of the parameter magnetic field, halos in the GeV energy band can remain space. We study the impact of the source properties (red- longaftertheTeVblazar’sendofactivity.Alternatively,the shift,spectrum,anisotropy)inSection4.Then,inSection5, spectralanalysiscanalsoprovideconstraintsontheEGMF we explore the effects of the intergalactic medium (EBL, (e.g.D’Avezacetal.2007;Neronov&Vovk2010;Kachelrieß EGMF).TechnicalaspectsofthecodearepresentedinAp- 2010).Althoughmoststudiesfocusontheaverageintensity pendix A. andcoherencelengthofthefield,ithasbeenshownrecently thatanisotropiesintheimagesofpairhaloscouldalsopro- vide crucial information on the magnetic helicity (Long & Vachaspati 2015; Batista et al. 2016). 2 PHYSICS OF COSMOLOGICAL CASCADES Cosmological electromagnetic cascades involves three main Inthepastyears,allthreeeffectshavebeensearchedin- processes: pair production through photon-photon anni- tenselyinthedataofthespacegamma-raytelescopeFermi hilation, inverse Compton scattering, and propagation of and of Cherenkov air telescopes such as MAGIC, H.E.S.S, charged particles in a magnetised, expanding universe. All or VERITAS. However none of the methods has provided other processes are negligible. In particular, as long as pri- undisputed evidence yet. Cascade contribution to the GeV mary photons do not exceed 100 TeV and the extragalactic spectrumhasmostlyprovidedupperlimits,andmostBlazar magnetic field remains below B = 10−10 G, synchrotron observations remain compatible with no cascade emission cooling is orders of magnitude weaker than Compton cool- (Arlen et al. 2014). Such constraints however provide lower ing,andsynchrotronphotonsonlycontributeatlowenergy, limits on the EGMF intensity (Neronov & Vovk 2010). No below the infrared range (<0.02 eV). This section presents time delay has been clearly detected either, which also pro- a simple analytical view of the cascade physics. vides lower limits on the amplitude of the random compo- nentofthemagneticfield(Neronovetal.2011).Detectionof pair halos requires a very accurate modelling of the instru- 2.1 Propagation of particles in a magnetised, ment point spreadfunction(PSF)andhasnotgivenundis- expanding Universe puted results either (Krawczynski et al. 2000; Aharonian etal.2001;Abramowskietal.2014;Prokhorov&Moraghan CosmologicalcascadesdeveloponkpctoGpcscales.Onthe 2016). A detection in Fermi-LAT data sets was claimed re- largest scales, the geometry and the expansion of the uni- cently (Chen et al. 2015), but has not been confirmed by versemustbetakenintoaccount.Throughoutthispaperwe the Fermi collaboration yet. Much better constraints are assumeaΛ-CDMmodel.Acompletedescriptionofparticle expected from CTA (Meyer et al. 2016). trajectories can be found in appendix A1. However a few important points must be noted here. First, in an expanding Universe and in absence of any interaction, the particle Regardlessofthedetectionmethod,adeepunderstand- (photons and leptons) momentum p evolves with redshift z ing of the cascade physics is crucial to interpret observa- as p∝(1+z), also meaning that their energy continuously tionaldata.Inthepastdecade,thecascadephysicshasbeen decreases with time. In the specific case of photons, the en- investigated through fast, analytical (or semi-analytical) ergy scales with momentum: E ∝ p ∝ (1+z), providing methodsthatallowtoquicklycoveralargeparameterspace, γ thewell-knowncosmologicalredshift.Thecosmologicalevo- and Monte Carlo simulations. Although much slower, the lution of lepton energy is slightly more complex. However, latter have proven to be mandatory to derive quantitative in the limit of highly relativistic particles, it also scales as results and to interpret precise observations. Several codes (1+z). have been developed over the years but only the most re- The propagation of leptons is also affected by the ex- centincludethecosmologicalexpansionontheparticletra- tragalacticmagneticfield(EGMF).Inthiswork,weassume jectory (Taylor et al. 2011; Kachelriess et al. 2012; Arlen that no field is created or dissipated in the cosmological et al. 2014; Settimo & De Domenico 2015). To our knowl- voids,andthatitissimplydilutedastheuniverseexpands: edge, only one is publicly available (ELMAG Kachelriess et al. 2012) but the lepton trajectories in the magnetised, intergalactic medium is treated in a simple, 1D, diffusion approach.InthispaperwepresentanewMonteCarlocode 1 https://gitlab.com/tfitoussi/cascade-simulation MNRAS000,1–18(4January2017) Impact of cosmological cascades 3 Figure 2. Gamma-ray annihilation mean free path λγγ = ctγγ Figure 1. Comoving spectral energy distribution of target pho- (where tγγ is the mean cosmic time between two interactions) tons at a redshift z = 2, including the Cosmological Microwave as a function of their initial energy, and for different emission Background(blackline)anddifferentmodelsofEBL(colorlines). redshifts (solid lines). For comparison, the mean free path for a static universe (with properties frozen at their values at the initialredshift)isshownindashedlines.Theblacklineshowsthe B(z) ∝ (1+z)2 (see Durrer & Neronov 2013, eq. 22). In thin/thicktransitionwherethemeanfreepathequalsthesource that case, the Larmor radius evolves as: distance Ds and where the energy equals the absorption energy E . abs E (cid:18) E (cid:19)(cid:18) B (cid:19)−1 R = e ≈1.1(1+z)−1 e Mpc, L e B(1+z) 1TeV 10−15G c (1) The photon mean annihilation distance is plotted in Fig. 2 as a function of the initial photon energy, assuming where e is the lepton charge, E and B are respectively c e the EBL model from Dom´ınguez et al. (2011). The solid the lepton energy and the EGMF intensity at z =0. In co- lines show the results for an expanding universe while the moving coordinates, this means that the comoving Larmor dashed lines show the results for a static universe. Below 1 radius R (1+z) is constant and that for a uniform mag- L TeV, the absorption mean free path quickly becomes larger netic field, the perpendicular comoving motion of leptons than the typical distance of targeted sources (100 Mpc to is a pure circle. When Compton losses are included, lepton Gpc), so that only TeV photons are significantly absorbed. trajectories become converging spirals. Below 200-300 GeV, photons tend to travel over such large Inpractice,thecosmologicalmagneticfieldisexpected distances that two cosmological effects work in concert to to be highly turbulent (Caprini & Gabici 2015; Durrer & produce a diverging mean free path. First, the target pho- Neronov 2013). Although it should be described by a full ton density vanishes as (1+z)3 as the universe expands turbulentspectrum,itspropertiesareoftencharacterisedby (horizonevent).Second,gamma-raysphotonsaremoreand its intensity B and its coherence length λ . In this paper, B moreredshiftedbeforereachingthenextannihilationpoint, we will consider that the field structure can be modelled requiring higher and higher energy target photons. As the byuniformmagneticcellsofsameintensityandsizeλ but B EBLphotondensitydropsat10eV,themeanfreepathdi- withrandomorientations.Insideaparticularcell,thelepton verges at low energy. Photons emitted at low redshift, with comovingtrajectoriesbecomesimplehelicoidaltrajectories. anenergyof1TeVtravelafewhundredMpcbeforeproduc- ingpairs.Thisdistancedecreasesquicklyathigherenergyto reachfewMpctofewkpc.ConsideringablazarlikeMrk421 2.2 Photon absorption by the EBL (z = 0.0308, 135 Mpc) emitting very high-energy photons High-energy photons (of energy E ) annihilate with soft, upto100TeV,primarygamma-raysaretypicallyabsorbed γ ambientphotons.Theannihilationcrosssectionbeingmax- over a distance of a few Mpc. imal close to the threshold, the interaction is most efficient Since the density of the EBL photons decreases with with soft photons of energy ∼(m c2)2/E where m is the theirenergy,gamma-rayabsorptionismoreefficientathigh e γ e lepton mass and c is the speed of light. This explains why energy.Thetransitionfromopticallythintoopticallythick TeV photons are absorbed preferentially by eV photons of photon-photonabsorptionwheretheradiationbecomesfully the EBL (Gould & Schréder 1967). Fig. 1 shows six mod- absorbedoccursataninitialenergyE =(1+z)E where abs cut elsofEBLthatcanbefoundintheliterature(Franceschini E isthecorrespondingenergycut-offobservedinthespec- cut etal.2008;Dom´ınguezetal.2011;Finkeetal.2010;Kneiske tra at z = 0. Fig. 3 shows the cut-off energy as a function &Dole2010;Gilmoreetal.2012)andillustratestheuncer- of source redshift. Distant sources are more absorbed and tainty of EBL intensity and spectrum. It can be seen that their absorption occurs at lower energy. Significant differ- at z =2, the EBL photon densities can differ by one order ences are observed between the different EBL models. At of magnitude from one model to the other. large redshift (z > 1), the different EBL models are very MNRAS000,1–18(4January2017) 4 T. Fitoussi et al. After a flight of length x, the lepton energy is E (x) = e E0/(1+x/D0) where e ic 3m2c4 (cid:18) E0 (cid:19)−1 D0 = ≈367 (1+z)−4 e kpc, (7) ic 4σ ρ E0 1TeV T cmb e istheinitialComptoncoolingdistance,andE0isthelepton e energy at the production site. 2.4 Observables Basedonthesesimpleproperties,someanalyticalestimates fortypicalcascadeobservablescanbederived.Herewemake the following standard assumptions: - The source of primary gamma-rays is isotropic and mono-energetic with an energy E0. γ - The primary gamma-rays all annihilate at exactly the same distance λ as given in Fig. 2. γγ Figure3.Cut-offenergyEcut observedatz=0asafunctionof - The leptons are produced by photon-photon annihila- thesourceredshiftfordifferentEBLmodels. tion in the direction of the parent gamma-ray photon and have exactly half of its energy. different, and differences up to a factor 6 are observed in - ComptoninteractionsoccurintheThomsonregimeand the cutoff energy. At lower redshifts the EBL models are the leptons travel exactly over one mean free path λic be- consistent with each other and only differ by a factor 2 at tween two Compton scatterings. thelowestenergiesas,atsuchlowdistance,gamma-gamma - The leptons are deflected by a uniform magnetic field absorption occurs mainly above 10 TeV. The effects of the perpendicular to the motion. EBLmodelonthecascadeswillbediscussedinsection5.1. - Magnetic deflections occur locally on scales much smaller than the photon annihilation length and the source distance. 2.3 Compton scattering by the CMB - The scattered photons all get exactly the average en- High-energy leptons Compton up-scatter soft ambient pho- ergygivenbyequation5andareemittedinthepropagation tons to gamma-ray energies. In the Thomson regime, the direction of the scattering lepton. Compton cross-section does not depend on the energy, and - The absorption of these scattered photons (hereafter the scattering rate scales linearly with the target photon refereed to as first generation photons) and the associated numberdensity,henceComptonscatteringmostlyoccurson pairproductionisneglected.Thecontributionofhighergen- CMB photons. The CMB is modelled by a blackbody with eration particles is not considered. atemperatureT =(1+z)T (BlackcurveinFig.1), - Cosmological effects are neglected. cmb cmb,0 where T = 2.725 K is the temperature at z = 0. The cmb,0 Withtheaboveassumptionswecanderivethefollowing associatedCMBmeandensity,averageenergyandmeanen- cascade geometrical and distribution properties: ergy density are respectively: (cid:18)k T (cid:19)3 ncmb =16πζ(3) Bhccmb ≈411(1+z)3 cm−3, (2) Geometry: π4 Geometrical effects (halo effects and time delay) are due to (cid:15) = k T ≈6.34×10−4(1+z) eV, (3) cmb 30ζ(3) B cmb theleptondeflectionintheextragalacticmagneticfield.Two cases can be studied: ρ =n (cid:15) ≈0.26(1+z)4 eV.cm−3, (4) cmb cmb cmb - If the coherence length is large (λ (cid:29) D0) the mag- where ζ(3) ≈ 1.202, k , and h are the Apéry’s, the Boltz- B ic B netic field can be considered as uniform and the lepton de- mann’s, and the Planck’s constants respectively. In the flection after a travel distance x is: δ = (cid:82)xds/R (s). It Thomson regime, leptons of local energy E up-scatter soft 0 L e means that a lepton of initial energy E0 = E0/2 having photons to typical energy: e γ cooled down to energy E has been deflected from its origi- e 4(cid:15) E2 (cid:18) E (cid:19)2 nal trajectory by an angle: E = cmb e ≈3.23(1+z) e GeV. (5) γ 3m2ec4 1TeV D0 (cid:34)(cid:18)E0(cid:19)2 (cid:35) Between two Compton scatterings, the leptons travel a δ= 2Ri0c Ee −1 , (8) Comptonmeanfreepathλ =1/(n σ )≈1.19kpc(cor- L e ic cmb T responding to a scattering time of tic = λic/c = 3870 yr), whereRL0 (Eq.1)istheinitialLarmorradiusofthelepton. where σT is the Thomson cross section. And on average, As soon as the lepton has lost a significant fraction of its they loose energy at rate: energy, the previous equation reduces to: dE 4 (cid:18) E (cid:19)2 D e Bλ e = cσ ρ e . (6) δ≈ ic = c ic, (9) dt 3 T cmb m c2 2R 2E e L γ MNRAS000,1–18(4January2017) Impact of cosmological cascades 5 equations combined to Eq. 5 and 9 yield respectively to: (cid:16) τ (cid:17)−1(cid:18) B (cid:19)(cid:18) E (cid:19)−1 θ≈0.79◦ γγ γ , (13) 397.4 10−14G 1GeV (cid:18) λ (cid:19)(cid:18) B (cid:19)2(cid:18) E (cid:19)−2 ∆t≈65 γγ γ yrs, (14) 1.32Mpc 10−17G 1GeV where τ = D /λ is the annihilation optical depth. All γγ s γγ values are calculated for a source at z = 0.13 emitting primary photons at E0 = 100 TeV. These equations show γ the complementarity in the search for pair halos and pair Figure 4.Geometryoftheone-generationmodel. echoes. Pair halos can only be observed if they are larger than the instrument point spread function (PSF). For a typical PSF of 0.1◦, the above values correspond to mag- neticfieldslargerthan10−14 G.Hence,largemagneticfield where D and R are no longer the initial values but are strengthscanbeconstrainedthroughdetectionofpairhalos. ic L nowevaluatedlocallyatenergyE (cid:28)E0,correspondingto In this case, time delays are as long as 108 yrs and echoes e e photons scattered to energy E (Eq. 5). cannot be observed. On the other hand, still with the val- γ - If the coherence length is short (λ (cid:28) D0), the lep- uesabove,observableechoesshorterthan5yrrequiremag- B ic tons travel across many zones of a highly turbulent field. netic fields lower than 10−18-10−17 G. Hence low magnetic We assume that the field is composed of many cells of size field strengths can be constrained though detection of pair λ , with uniform field and random directions. In each cell, echoes. In that case, pair halos are typically smaller than B theleptonsaredeflectedbyanangle∼λ /R inarandom 0.0001◦ and cannot be resolved. B L direction. Considering this process as a random walk leads Astheleptonenergydecreases,boththedetectionangle to: andthetimedelayincreaseuntilthemaximaldeflectionδ= π/2 is reached for which θ ≈ 1/τ = 5.7◦(τ /10)−1. max γγ γγ δ= (cid:112)2DRi0c0λB (cid:34)(cid:18)EEe0(cid:19)2−1(cid:35) Ttohnisecnoerrrgeys,pothnedsLtaormthoermraadxiiumsablehcaolmoesiszes.mAatllearlotwhaenr ltehpe- √ L e magnetic coherence length. The leptons are trapped by the ≈ DicλB, (10) magnetic field and cannot travel farther. This leads to the 2RL formationofacloudofe+−e−pairsaroundthesource.The size of the observed halo then corresponds to the physical wherethelastapproximationisobtainedsimilarlytoEq.9. extension of the pair cloud, i.e. λ . γγ In both cases, the magnetic deflection is a function of the lepton energy δ(E ) and of the secondary gamma-ray en- Distributions: e ergy δ(E ). In the following, we will concentrate on large- γ Inthefollowing,unlessotherwisespecified,alldistributions coherencecase(λ (cid:29)D0).However,similarconstraintsare B ic are normalised to one single primary photons emitted. easily obtained for short coherence lengths. The cascade spectrum produced by one sin- The geometrical properties of the cascade (extension, gle high-energy photon is computed as dN/dE = time delay) can be derived from the magnetic deflection γ 2(dN/dt) /(dE /dE )/(dE /dt), where the factor 2 ac- angle (e.g. Neronov & Semikoz 2007; Dermer et al. 2011). ic γ e e counts for the two leptons produced by one single primary TheyareillustratedinFig.4.Intheone-generationapprox- photon. Noting that the number of photons up-scattered imation, halo photons observed with a finite angle θ were byoneleptonperunittimeis(dN/dt) =c/λ ,andusing emitted out of the line of sight and then deflected back to ic ic Eq. 5 and 6 to derive the last terms leads to: theobserver.Assumingtheleptondeflectionoccursonvery short distances compared to the photon absorption length, dN m c2(cid:114)3E (cid:18) E (cid:19)1/2 the detection angle and time delay are: E2 = e γ ≈556 γ (1+z)−1GeV. γdE 2 (cid:15) 1GeV γ cmb (cid:18) (cid:19) (15) λ λ θ=arcsin γγ sinδ ≈ γγδ, (11) Ds Ds It is a simple power-law dN/dEγ ∝ Eγ−Γ with index Γ = δ2 3/2. If unabsorbed, this spectrum extends up to the en- c∆t=λ (1−cosδ)−D (1−cosθ)≈λ , (12) γγ s γγ 2 ergy of photons scattered by the highest-energy leptons, i.e. E = 0.8(E /1 TeV)2 GeV. However, for distant γ,max γ,0 where D is the distance to the source, λ is the annihi- sources,suchapower-lawspectrumistypicallycutatlower s γγ lation distance of the primary photons, and where the last energies by photon absorption (Fig. 3). In principle, a new approximations were obtained in the small-angle approxi- generationofparticlesisthenproduced.Higherphotonsgen- mation (λ (cid:28) D ,δ (cid:28) 1). These relations are illustrated eration have often been neglected although they may con- γγ s in Figs. 5 in the case of large coherence length. tribute significantly to the overall spectrum (see section 3). Ashigh-energyleptonstravelandgetdeflected,thede- TheangulardistributioncanbecomputedasdN/dθ= tection angle and time delay decrease as energy increases. (dN/dE )(dE /dδ)(dδ/dθ).UsingEq.15,9,and11toesti- γ γ In this regime, the small angle approximation is well satis- mate the three factors respectively allows to derive the an- fied and for large coherence length (λ (cid:29) D0) the latter gular distribution of first-generation photons produced by B ic MNRAS000,1–18(4January2017) 6 T. Fitoussi et al. one single primary photon, for large coherence length and cross-sections. Particles with updated parameters and new small angle: particles are then stored for later treatment. The physical parameters of the particles reaching the Earth without any dN (cid:18) 3 τ θ(cid:19)1/2 θ =m c2 γγ interaction are stored for post-processing. dθ e 2λic(cid:15)cmb ecB - Theextragalacticmagneticfieldismodeledbydividing =1972.3(cid:16) τγγ (cid:17)1/2(cid:18) B (cid:19)−1/2(cid:18) θ (cid:19)1/2. (16) the comoving space into a number of cells with size λB, 397.4 10−15G 1◦ defining regions of uniform magnetic field. For each cube a random magnetic field direction is computed and is kept Identically,writingdN/dt=(dN/dEγ)(dEγ/dδ)(dδ/dt)and for the entire simulation. The EGMF is set by its strength using Eq. 12 gives the following distribution of time delays: at redshift z = 0. For very short coherence lengths, the dN (cid:18) 3 (cid:19)1/2(cid:18) c∆t (cid:19)1/4 leptoninteractiondistance(λIC ≈1kpc)canbecomelarger ∆t dt =mec2 4λ (cid:15) e B 2λ than λB. In such a case, the motion to the next interaction ic cmb c γγ location is divided in shorter steps of length a fraction of =437.2(cid:18) λγγ (cid:19)−1/4(cid:18) ∆t (cid:19)1/4(cid:18) B (cid:19)−1/2. λB, to ensure that the leptons are deflected by all cells. 1.32Mpc 106yr 10−15G Thefollowingsubsectionspresenttheresultsoftestsim- (17) ulations.Forthesakeofcomparison,wefirstconsiderasim- These distributions are pure power-laws and do not show ple canonical model consisting in a mono-energetic source any specific angular (or time) scale because they are com- (100 TeV) emitting isotropically at z=0.13 (around 557 posedbythecontributionsofphotonsdetectedwithallpos- Mpc). We use the EBL model of Dom´ınguez et al. (2011). sibleenergies.However,asdiscussedbelow,realobservations With this set-up the characteristic annihilation distance of areobtainedinalimitedenergyrange,withagivenaperture primary photons is λ = 1.32 Mpc. The EGMF is set to γγ angle,andwithafiniteobservationduration.Thisproduces B =3×10−16 G with a coherence length of λ =1 Mpc. B characteristic scales in energy, arrival angle and time delay, We focus on the three main observables characterising thatappearascutsorbreaksinthecorrespondingdistribu- the photons detected on Earth: their energy, their arrival tions. angle,andtheirtimedelay.Wederivethe3Dphotondistri- bution in this space of parameter. We also concentrate on thecontributionofthedifferentgenerationsofparticles,i.e. 3 CODE DESCRIPTION AND TEST CASES their rank in the cascade tree. In the following, the zeroth generation corresponds to the primary photons emitted at Theanalyticalapproximationspresentedintheprevioussec- thesource.Thesephotonsannihilate,producepairs,thatin tions are useful to understand the physics of the cascades turn produce the photons of first generation. Again, these and to obtain orders of magnitude estimates, but accurate can produce pairs which up-scatter target photons produc- calculationsrequirenumericalsimulations.Forthispurpose, ingthesecondphotongeneration,andsoon.Asweareonly wehavedevelopedanewMonteCarlocode.Inthissection, interested in the detection of photons, their energy will be we outline the main features of this code and compare its writtenE,wherethesubscriptγ hasbeendroppedforsim- resultstotheanalyticalpredictionspresentedinsection2.4 plicity. as well as to other published results. 3.2 Correlations between observables 3.1 Code and simulation set-up In the framework of the approximations used in the one- Our Monte Carlo code is designed to track the propagation generation analytical model, the photon energy, detection ofallparticlesandtoreproducethepropertiesofthecascade angleandtimedelayarelinkedexactlythroughsimplerela- as precisely as possible. A complete description of the code tions. When these approximations are relaxed a significant can be found in appendix A. However, the main algorithm scatter is expected. Fig. 5 show the three possible correla- is summarised here (see also Fig. A1): tions between the photon energy E, the detection angle θ, - Primary particles (photons or leptons) are launched at and the time delay ∆t. In addition, the average behaviour agivenredshiftandwithagivenenergy.Particleenergycan is plotted in solid line for each generation by binning the alsobesampledfromapower-lawdistribution.Inthispaper x-axis (with 4 bins per decade) and averaging the y-values. will focus only on cascades initiated by photons. Theanalyticalestimatesoftheone-generationmodel(from - Interaction distances are generated randomly accord- Eq. 11 and 12) are shown for comparison as blue, dashed ing to the probability distribution given by the exact cross- lines. The expected trends are recovered: the smaller the sections (e+ − e− pair production for a photon, Klein- energy,thelargertheobservationangleandthearrivaltime Nishima cross-section for a lepton) and taking into account delay. As the leptons cool down, they produce more and the cosmological evolution of the target particles. moregamma-rayphotonsperunittime.Asaresult,theto- - The particles are propagated taking into account all tal cascade emission is dominated by a very large number cosmological effects (redshift, expansion). In particular, the of photons with low-energy, large angle, long time delay, as transportofleptonsintheEGMFiscomputedasdescribed canbeseeninFig.5.However,severalimportantpointscan inappendixA1.Iftheparticlesinteractbeforereachingthe be made. Earth,theinteractionoutcomeiscomputed.Energyanddi- The averaged values for the first-generation photons rectionoftheoutcomingparticlesaregeneratedaccordingto (blue line) are consistent with the analytical results of sec. the probability distributions given by the exact differential 2.4: both the saturation at low energy and the slopes in MNRAS000,1–18(4January2017) Impact of cosmological cascades 7 the small-angle regime are well recovered. The saturations observed in time correspond to the maximal delay (photons that are emitted away from, and scattered back to the observer) while the saturation observed in detection angle corresponds to photons deflected by an angle of δ = π/2. The power-law regimes correspond to Eq. 13 and 14. Most of the observed deviations at high energy (e.g. fluctuations and peaks) come from the limited statistics of the simulations: averaged values can be contaminated by a few photons with very large values (typically photons scattered by EBL targets instead of CMB photons). Physical processes that are not taken into account in the analytical estimates (Klein-Nishina regime, dispersion in the annihilation distance,non-uniformmagneticfield,energydispersionaround the averaged value ...) are also responsible for deviations to the approximations, but the effects are weaker. These results also show that the one-generation model underestimates the detection angle and the time delay by at least two orders of magnitude when the contribution of second-generation photons is significant. Indeed the highest energy, second-generation photons are almost all produced at the location where the parent leptons which emitted them where produced. As these photons have lower energy than the primary photons, their annihilation distance is larger λgen=1 (cid:29) λgen=0. As a result, the highest energy, γγ γγ second-generation photons are typically produced at a distance λgen=1 from the source. The geometry remains how- γγ ever similar so that an estimate for the second-generation observable quantities can be found by substituting λ by γγ λgen=1 intheresultsofsec.2.4.Forprimaryphotonsat100 γγ TeV, the highest-energy, first-generation photons (8 TeV) have mean free path λgen=1 = 117 Mpc. This is shown in γγ Fig. 5 as green dashed lines. As can be seen, this new estimate matches well the average results of second-generation dominated cascades. Weemphasisethatasinglegenerationneverdominates atallenergies,anglesandtimedelays.Inourcanonicalsim- ulation for instance, only the low-energy spectrum is domi- natedbysecond-generationphotonswhilethehighestenergy photons are mostly first-generation photons (see Fig. 6). Thisiswhytheaveragetimedelaydropsbelowthesecond- generation estimates at high energy. Inprinciple,theratiooffirsttosecondgenerationpho- tons depends on the energy of primary photons and the sourcedistance.Ifprimaryphotonshaveenergysmallerthan or close to the absorption energy E , then the absorption abs is so weak that the production of second-generation photons is quenched, and the cascade is dominated by first- generation photons (this will be illustrated in Fig. 12 for instance). However, the results presented here remain gen- Figure 5. Top panel: Correlation between the time de- eralassoonastheenergyofprimaryphotonsissignificantly lay ∆t and energy E of detected photons. The density map larger than the absorption energy (see also Berezinsky & shows the number of photons per unit energy and time delay Kalashev 2016). (E∆t)d2N/(dEd∆t) with a log color scale. The blue, green and Inanycase,thereisalargedispersionaroundtheaver- blacksolidlinesshowtheaveragetimedelayforfirst-generation age quantities, showing that the latter might be of limited photonsonly,second-generationphotons,andallphotons,respec- practical use. tively. Theblue and greendashed lines show the analyticalesti- mates for the first and second generations only (see sec. 2.4). Withthesamenotations,themiddleandbottompanelsshow 3.3 Photon distributions theenergy-detectionangleanddetectionangle-timedelaycorre- lationsrespectively. The energy, angle and time distributions are obtained by integrating over the other two parameters. In this section we focus on photonswithenergyabove1GeV whichcorre- MNRAS000,1–18(4January2017) 8 T. Fitoussi et al. Figure 6. Full spectrum for an isotropic and mono-energetic Figure 7.DetectionangledistributionforE >1GeVphotons. source (100 TeV) at z=0.13 (black line). The blue, green and ThecolorsarethesameasinFig.6.Theblue,dashedlineshows redlinesshowthecontributionsofgenerations1,2and3respec- theanalyticalestimateofEq.16. tively. The analytic expression (Eq. 15) is shown in blue dashed line. The dotted-dashed, and dotted lines show the results from Elmag (Kachelriess et al. 2012) and Taylor et al. (2011) respectively. sponds to the typical energy range of gamma-ray observa- tories. Thespectrumofourfiducialmodel,integratedoverthe entiresky,overallpossiblearrivaltimes,andnormalisedto one primary photon, is shown in Fig. 6. It is compared to Fig.1ofTayloretal.(2011)2andtotheresultsoftheElmag code (Kachelriess et al. 2012) using the same set-up. Spec- tral shapes provided by the different codes are compatible each other. All primary photons at 100 TeV are absorbed and produce the first-generation spectrum (blue line). Be- low E analytical estimates (Eq. 15) are well reproduced cut by the first-generation population, in spite of the approximations made. Above 1 TeV the photon of the first generation are absorbed (Fig. 3) and produce a second generation whichdominatesthespectrumbelowabout100GeV(green Figure 8. Distribution of time delays for E > 1 GeV photons. line).Thespectrumofthesecondgenerationissimilartothe spectrumofthefirstgeneration(dN/dE ∝E−3/2)exceptit ThecolorsarethesameasinFig.6.Theblue,dashedlineshows theanalyticalestimateofEq.17. is softer in the energy range shown in this figure (it can be shownthatdN/dE ∝E−7/4).Afewsecond-generationpho- tonsarealsoabsorbedandproduceaweakthird-generation Assecond-generationphotonshavealargermeanfree-path, population(redline)whichdoesnotcontributetothetotal they typically arrive with larger angles. Interestingly, the spectrum. distribution is dominated by second-generation photons at The angle distribution integrated over energies E > 1 observable scales (θ > 0.1◦). This result remains general GeV and over all arrival times, normalised to one primary as long as the energy of primary photons is significantly photon is shown in Fig. 7. The emission is peaked at the larger than the absorption energy E (i.e. photons have center and decreases as a power-law with increasing angle. abs absorptiondepthτ (E0)>>1).Primaryphotonswithen- Atsmallangles,thedistributionoffirst-generationphotons γγ ergycomparabletotheabsorptionenergy(τ (E0)≈1)are iswellapproximatedbytheanalyticalestimategiveninEq. γγ weakly absorbed and do not produce high-generation pho- 16: dN/dθ ∝ θ−1/2. However, only photons above 1 GeV tons. For sources with an extended intrinsic spectrum, the are considered here. Hence, low-energy photons with large contribution of first- and second-generation photons to the angles are not observed. And, as compared to the analyti- angular distribution is more complex (see sec. 4.2). calestimate,theangulardistributiondropsatatypicalsize The time delay distribution integrated over energies thatdependsontheminimalenergyandthemagneticfield. E > 1 GeV and all detection angles, and normalised to one primary photon, is shown in Fig. 8. As time evolves af- 2 Thenormalisationsofthepublishedresultswherechosenarbi- ter a source flare, less and less photons are observed. The trary one-generation model (Eq. 17) provides a good estimate of MNRAS000,1–18(4January2017) Impact of cosmological cascades 9 also shows the instantaneous spectra observed at present time, if the source (such as an AGN) has been active for a timet =τ inthepast,withconstantluminosity(Dermer act etal.2011).Astheactivityperiodincreases,weareableto detectsecondariesproducedbyprimariesemittedearlier.As leptonshadmoretimetocooldown,thesesecondarieshave lowerenergy.Asaresult,longactivitysourceshavespectra that extend to lower energy. 4 SOURCE PROPERTIES The simple case presented in section 3 allows us to understand general behavior of electromagnetic cascades. How- ever, several effects must be included to produce realistic cascades. Gamma-ray sources (AGNs, GRBs) do not emit photons at a single energy or isotropically but instead produce non-thermal, beamed radiation. Here we investigate Figure 9. Total spectrum of an flaring event integrated over the following intrinsic properties of the source: a finite exposure time (t =τ), or equivalently: instantaneous obs spectrumofasourcethathasbeenactiveforagivetime(tact=τ) - Redshift z. inthepast. - Intrinsicspectrum:hereweconsiderpower-lawspectra in the form dN/dE ∝E−Γ for 100 MeV(cid:62)E (cid:62)E . max - Emissionprofile:hereweassumeadiskemission,i.e.an thefirst-generationdistributionfortimedelaysrangingfrom axisymmetric angular distribution dN /dΩ uniform up to e e about1monthto100years.Shortertimedelayscorrespond a given half-opening angle θ observed with an angle θ jet obs tohighenergyphotonsthatareabsorbed,producingadrop away from its axis (see Fig. 4). below the analytical estimate. Longer time delays correspond to low-energy photons below the selection criterion Inthefollowingweusesimulationparameterscorresponding E >1 GeV, producing a cut of the distribution above 1000 to the Blazar 1ES0229+200 (Tavecchio et al. 2009; Taylor years.Interestingly,accessibletimedelays(∆t<1yr)after et al. 2011; Vovk et al. 2012) with a redshift z = 0.14 cor- aflaringevent(suchasaGRB)areshortenoughtobedom- responding to a distance of 599 Mpc, and a hard spectrum inatedbyfirst-generationphotonsonly,allowingfortheuse with Γ=1.2. The unconstrained maximal energy of the in- of simple formulae to derive constraints from the potential trinsicspectrumissettoEmax =100TeVandtheemission detection of pair echoes. is assumed to be isotropic. The EGMF has an averaged in- The three distributions presented on Figs. 6, 7 and 8 tensity of B = 10−15 G, and a coherence length λB = 1 (dN/dE,dN/dθ,anddN/dt)areglobaldistributionsinthe Mpc. We use the EBL model of Dom´ınguez et al. (2011). sensethattheyareintegratedoverlargerangesofthe2oth- In this section, the spectra are normalised to L0, the ersquantities.Forinstance,thespectrumshowninFig.6is intrinsic luminosity of the source as observed at z =0 (i.e. theintegratedoveralldetectionanglesandallarrivaltimes. the intrinsic luminosity decreased by a factor 1+z). However,inamorerealisticsituation,onlylimitedrangesof thesequantitiesareaccessible.Aslow-energyphotonsarrive 4.1 Source redshift with large angle, large time delay, and are products of pri- maryphotonsemittedfarawayfromthelineofsight,anyof The spectral evolution with redshift is shown in Fig. 10 for the following effects will damp the spectrum at low energy, z = 0.04 (D ∼ 175 Mpc), z = 0.14 (D ∼ 599 Mpc), s s while cutting the large angle and large time delay part of z = 1 (D ∼ 3.4 Gpc), and z = 2 (D ∼ 5.3 Gpc). As ex- s s the associated distributions: pected,thecut-offenergydecreaseswithincreasingdistance as the column density of target photons increases. At high (i) Iftheinstrumentisonlysensitiveaboveagivenenergy redshift, we find that the absorption depth goes simply as (see discussion before). (ii) If the instrument aperture is limited. τγγ ∝E2producingasuper-exponentialcutoff∝e−E2/Ec2ut. However the spectra of nearby sources show a more com- (iii) Iftheexposuretimeisfiniteafteranimpulsiveflaring plex and harder absorption cutoff. In our setup, the maxi- event. malenergyofprimaryphotons(100TeV)islargeenoughto (iv) If the source emission is beamed within a limited generate an efficient cascade. As a result, almost all shown opening angle. spectra are dominated by second-generation photons. Only TheeffectoffiniteexposuretimeisillustratedinFig.9(see whenthesourceiscloseenough(z=0.04),theabsorptionis alsoIchikietal.2008).Thiscanbeinterpretedintwoways. weakenoughtoquenchtheproductionofsecond-generation Ifthesourceproducesastrong,impulsiveflaringevent(such photons. As discussed by Berezinskii & Smirnov (1975) the asaGRBforinstance),thisfigureshowstheintegratedspec- spectrumsoftensasthegenerationorderincreases,whichis tra as data is accumulated from the detection of the unab- consistent with the harder spectrum observed at z = 0.04. sorbed,primaryphotonsuptotimet =τ.Astimeevolves The intrinsic spectrum used here is hard (Γ=1.2), so that obs lower energy photons are detected and the low-energy part mostoftheintrinsicluminosityisconcentratedatthehigh- of the spectrum builds up slowly. Alternatively, this figure estenergies(E ∼E ).Mostofthespectrumisthenfully 0 max MNRAS000,1–18(4January2017) 10 T. Fitoussi et al. Figure 10. Full spectrum of sources at different redshifts (z = 0.04, 0.14, 0.4, 1 and 2), normalised to the intrinsic luminosity attenuatedbyafactor1+ztoaccountfortheuniverseexpansion. Figure 12. Top panel: Observed spectra for Emax =100 TeV and different spectral indices Γ (solid lines) . The contributions of the primary source and the cascade are shown in dotted and dashed lines respectively. Spectra are normalised to L0, the in- trinsicluminosityattenuatedby1+z.Bottom panel:Spectra forΓ=1.2anddifferentmaximalenergiesEmax. asmallerangleasseenfromamoredistantobserver(seefor instance Eq. 11 for first-generation photons: θ ∝ λ /D ). γγ s In contrast, the time delay does not suffer from any geometrical dependence on distance (see for instance Eq. 12 Figure 11.Averagearrivalangle(top)andtimedelay(bottom) forfirst-generationphotons)andshowsonlylittleevolution of photons with energy E > 1 GeV, as a function of the source withredshift.Tofirstorderhowever,thecosmologicalevolu- redshift, for different EGMF strengths (B = 10−13,10−15, and tionoftheuniversealsoinfluencestheangularsizeandtime 10−17 G). delayofsecondaryphotons(namelythroughλ (z),λ (z), γγ ic B(z), and E(z)), explaining the remaining evolution. absorbedandredistributedascascadecontribution.Asare- sult,theamplitudeoftheobservedspectraisalmostinsensi- 4.2 Source spectrum tivetotheabsorptionenergy,i.e.alsotothesourceredshift. The evolution of the angular distribution and time de- Fig. 12 shows the observed spectra when the source in- lays of photons with energy E > 1 GeV, are illustrated by trinsic spectrum is changed. The top and bottom panel theiraveragevaluesinFig.11fordifferentEGMFstrengths. show the results for different spectral indices (Γ=1.2, 2 Both the halo extension and the time delay increase with and 2.2, E = 100 TeV) and different maximal energies max magneticfield,asleptonsofgivenenergyaremoredeflected (E =10, 50, 100 TeV and 1 PeV, Γ=1.2) respectively. max by stronger fields. Their evolution with redshift is the re- At source distance z = 0.14, photons with energy higher sultofseveraleffects.Thehaloextensiondecreaseswithdis- thanafewTeVareabsorbedandredistributedtowardslow tance.Tozerothorder,itissimplyduetogeometricaleffects: energies. For hard spectra (Γ < 2), many primary photons thesameannihilationdistancetothesourcecorrespondsto are absorbed. This induces a strong cascade which domi- MNRAS000,1–18(4January2017)