Origin of the ankle in the ultra-high energy cosmic ray spectrum and of the extragalactic protons below it
OOrigin of the ankle in the ultrahigh energy cosmic ray spectrum,and of the extragalactic protons below it
Michael Unger,
1, 2, ˚ Glennys R. Farrar, : and Luis A. Anchordoqui
3, 4, 5, ; Center for Cosmology and Particle Physics, Department of Physics, New York University, NY 10003, USA Karlsruher Institut f¨ur Technologie, Institut f¨ur Kernphysik, Postfach 3640, 76021 Karlsruhe, Germany Department of Physics and Astronomy, Lehman College, City University of New York, NY 10468, USA Department of Physics, Graduate Center, City University of New York, 365 Fifth Avenue, NY 10016, USA Department of Astrophysics, American Museum of Natural History, Central Park West 79 St., NY 10024, USA (Dated: Aug 28th 2015)The sharp change in slope of the ultrahigh energy cosmic ray (UHECR) spectrum around 10 . eV(the ankle), combined with evidence of a light but extragalactic component near and below the ankleand intermediate composition above, has proved exceedingly challenging to understand theoretically,without fine-tuning. We propose a mechanism whereby photo-disintegration of ultrahigh energynuclei in the region surrounding a UHECR accelerator accounts for the observed spectrum andinferred composition at Earth. For suitable source conditions, the model reproduces the spectrumand the composition over the entire extragalactic cosmic ray energy range, i.e. above 10 . eV.Predictions for the spectrum and flavors of neutrinos resulting from this process are also presented. I. INTRODUCTION
The cosmic ray spectrum spans roughly eleven decadesof energy, 10 eV À E À eV and has three major fea-tures: the steepening of the spectrum dubbed the “knee”at « . eV [1], a pronounced hardening of the spec-trum at E « . eV, the so-called “ankle” feature [2–4], and finally a cutoff around 10 . eV [3, 5]. Three ad-ditional more subtle features have been reported betweenthe knee and the ankle: A hardening of the spectrum ataround 2 ˆ eV [6–9] followed by two softenings at „ . eV [6, 7] and and 10 . eV [2, 8–11]. The latteris traditionally referred to as the “second knee”.The variations of the spectral index reflect various as-pects of cosmic ray production, source distribution andpropagation. The first and second knee have straightfor-ward explanations, as reflecting the maximum energy ofGalactic magnetic confinement or acceleration capabilityof the sources, both of which grow linearly in the charge Z of the nucleus; the first knee being where protons dropout and the second knee where the highest- Z Galacticcosmic rays drop out. As the energy increases above thesecond knee to the ankle, the composition evolves fromheavy to light [12] while the cosmic ray arrival directionsare isotropic to high accuracy throughout the range [13–15]. Finally, as the energy increases above the ankle,not only does the spectrum harden significantly, butthe composition gradually becomes heavier (interpretingthe data using conventional extrapolations of accelerator-constrained particle physics models) [16, 17].This observed evolution in the extragalactic cosmic raycomposition and spectral index presents a major conun-drum. A pure proton composition might be compati- ˚ [email protected] : [email protected] ; [email protected] ble with the observed spectrum of extragalactic cosmicrays [18] when allowance is made for experimental un-certainties in the energy scale and the fact that the reallocal source distribution is not homogeneous and contin-uous [19] (although the sharpness of the ankle is difficultto accommodate), but a pure proton composition is in-compatible with the depth-of-shower-maximum ( X max )distributions observed by Auger [16, 17] unless currentextrapolations of particle physics are incorrect. More-over, a fit of the spectrum with a pure proton compositionseems to require a very strong source evolution [20] whichleads to a predicted neutrino flux in excess of experimen-tal limits [21]. On the other hand, models which fit thespectrum and composition at highest energies, predict adeep gap between the end of the Galactic cosmic raysand the onset of the extragalactic cosmic rays [22–27].Models can be devised to fill this gap, but fine-tuning isrequired to position this new population so as to just fitand fill the gap [28–30].Here we offer a resolution to this conundrum, byshowing that “post-processing” of UHECRs via photo-disintegration in the environment surrounding the source,can naturally explain the entire spectrum and composi-tion. In our model, extragalactic cosmic rays below theankle are predominantly protons from nucleons knockedoff higher energy nuclei in the region surrounding the ac-celerator, and the spectrum and composition above theankle are predominantly dictated by the accelerator andpropagation to Earth. The model makes distinctive pre-dictions about the spectrum and flavor ratios of neutri-nos, which should enable it to be tested. If the ankle andthe protons below it arise on account of our mechanism,we obtain a new constraint on UHECR sources beyondthe Hillas criterion and total-energy-injection require-ments, namely that the environment around the sourcehas the conditions giving rise to the required amount ofphoto-disintegration.Up until now, photo-disintegration (PD) has beenmainly considered as a danger inside the accelerator, a r X i v : . [ a s t r o - ph . H E ] N ov c o s m i c r a y source environment EBL/CMB detection FIG. 1. Illustration of our model calculation: Sources (yellow stars) inject cosmic rays with a power law in energy, into asurrounding region of radiation and turbulent magnetic fields. After propagation through this local environment and thenintergalactic space, these cosmic rays and their spallation products are detected at Earth. The photon energies in the sourceenvironment are characteristically of much higher energy than in the extragalactic background light. as it would cut off the cosmic ray spectrum at energiessuch that the PD interaction length and the accelerationlength are comparable. Since the acceleration length in-creases with energy, whereas the PD interaction lengthgenerally decreases with energy, photo-dissociation actsas a low-pass filter . The insight underlying the mech-anism we propose, is that if the primary locus of PDis outside the accelerator, PD generally acts as a high-pass filter , permitting the highest energy cosmic rays toescape unscathed while the lower energy ones are disinte-grated inside the source region, generating nucleons withenergy 1 { A of the original nucleus of mass A . As we shallsee, these spallated nucleons naturally produce the anklefeature, explain why extragalactic cosmic rays below theankle are protonic, and account for the spectral indexbelow the ankle. Examples of systems in which the ac-celerator is embedded in a photon field and the cosmicrays are trapped by magnetic fields in that environmentcould be the dusty torus surrounding an active galacticnucleus or the interstellar medium of the star-formingregion surrounding most young pulsars; see also [31–38].The basic setup of our phenomenological model is illus-trated in Fig. 1.The layout of the paper is as follows. In Sec. II weintroduce our model and in Sec. III we compare its pre-dictions with experimental data. Details about particlepropagation and the calculation of multi-messenger sig-natures are given in the appendices. Section IV containsour conclusions. II. FORMATION OF THE ANKLE
To illustrate the mechanism we have identified to cre-ate the ankle and generate protons below, consider asystem in which the accelerator (also referred to as thesource) is embedded in an environment in which the cos-mic rays are confined for some time by magnetic fieldswhile interacting with the ambient radiation field. Ouressential simplifications are: (i) a fast acceleration mech- anism and/or a low photon density inside the accelerator, (ii) no energy is lost except through an interaction, andwhenever a nucleus interacts it loses one or more nucleonsby photo-disintegration or photo-pion production (in thiscase the nucleus loses a fraction of its energy correspond-ing to the reduction in its nuclear mass); (iii) a cosmicray either escapes without changing energy, with a rate τ esc , or the cosmic ray interacts one or more times beforeescaping; (iv) τ esc and τ int are independent of positionin the source environment and depend only on t E, A, Z u of the nucleus. In this approximation the number of nu-clei in a given energy range and with a specified t A, Z u decreases exponentially with time, with τ “ p τ ´ ` τ ´ q ´ . (1)A fraction η esc “ p ` τ esc { τ int q ´ (2)of the particles escape without interaction and the restinteract before escaping, so η int “ ´ η esc . Note that η esc and η int depend only on the ratio of the escape andinteraction times, but not on the absolute value of eitherof them.A simple analytic treatment is instructive. To illus-trate the low/high-pass filter mechanism, consider thecase that the escape and interaction times are both powerlaws in energy, τ esc “ a p E { E q δ and τ int “ b p E { E q ζ . (3)Then η esc p E q “ ` ` R p E { E q δ ´ ζ ˘ ´ , (4)where R “ a { b is the ratio of the escape and inter-action time at reference energy E . When δ ą ζ , thesource environment acts as a low-pass filter on the parti-cles injected from the accelerator, leading to a cutoff inthe escaping spectrum at high energies. This situation -9 -8 -7 photo-disintegration [ a . u .] t c -9 -8 -7 photo-pion production lg(E/eV) -9 -8 -7 total narrow width approximationnumerical integration FIG. 2. Interaction times of Si in a broken power-law pho-ton field with parameters α “ , β “ ´ ε “ .
11 eV.Top panel: photo-disintegration, middle panel: photo-pionproduction, bottom panel: sum of the two processes. Theresults of numerical integration using detailed cross sectionsare shown as thick solid lines, while those of the narrow-resonance-approximation (detailed in Appendix B) are dis-played with thin dashed lines. is typical of leaky box models of diffuse acceleration attime-independent shocks [39–41] where δ ą δ ă ζ leading to a high-pass filter on the energy spectrum of injected nuclei:the lower the energy, the more time the nuclei have tointeract before escaping, leading to a hardening of thespectrum and lightening of the composition of nuclei es-caping the region surrounding the source. The spallatednucleons have energies of E “ E A { A ; these nucleons aremost abundant at low energies and have a steeper spec-trum ∝ p ´ η esc p E ˚ A qq . Thus the high-pass scenarioleads naturally to an ankle-like feature separating thenucleonic fragments from the remaining nuclei. The nor-malization and slope of the spectrum of spallated nucle-ons relative to that of the primary nuclei is determinedby how thoroughly the primary nuclei are disintegrated,which is governed by the ratio of escape and interactionlengths of the most abundant primaries. To obtain a more realistic treatment of the interac-tion time, we must specify the shape of the spectrum ofthe target photons. In our work to date we have con-sidered: (i) a broken power-law (BPL), characterized byits peak energy (cid:15) and lower and upper spectral indices α, β (this is a simplified representative of non-thermalemission that allows for analytic calculation as discussedbelow and in Appendix B); (ii) a black-body spectrum; (iii) two types of modified black-body spectrum, whichresult from a reprocessed black-body in a dusty environ-ment [42]. Details are given in Appendix A. For suchpeaky photon spectra the interaction time does not havethe simple representation of (3) but it does have a ratheruniversal structure. In our actual calculations we adopt anumerical integration of Talys [43, 44] and
Sophia [46]cross sections using [47], but the analytic expression for τ int derived in Appendix B for the BPL in the narrow-width approximation for the interaction cross sections,is qualitatively similar and useful for understanding. Ascan be seen in Fig. 2 the folding of a single resonance witha broken power-law spectrum leads to a “V” shape curvefor τ int in a log-log plot for both photo-disintegration (toppanel) and photopion production (middle panel). Com-bining both processes in narrow-resonance approxima-tion yields an interaction time with a “W” shape, whilenumerical integration including the plateau for multi-pion production softens the “W” to what we shall refer toas an “L” shape for brevity, a shown in the bottom panelof Fig. 2. As evident from Fig. 2, below the inflectionpoint for photodisintegration E b , the narrow-resonanceapproximation is good, while from the full numerical in-tegration in the high-energy region τ int is roughly con-stant, so using the BPL spectrum, we have the approxi-mate representation: τ int p E q « τ b " p E { E b q β ` E ď E b E ą E b , (5)where formulae for τ b and E b are given in Appendix B,and the parameter values for photodisintegration are tobe used.Returning to the discussion of τ int in (3) with (5) yieldsthe fraction of nuclei which escape without interaction ina peaky photon spectrum. It is straightforward to seethat if δ ă η esc has the properties of a high-pass filter. These conclusions do not depend on the exactshape of the photon spectrum. As can be seen in Fig. 13of Appendix A, the interaction times flatten to an L-curve as well if the photon density is assumed to followa (modified ) black body spectrum. III. COMPARISON WITH EXPERIMENTA. Fiducial Model
As our fiducial example, we adopt a broken power-lawphoton spectrum as a simplified representative of non-thermal emission given by n p ε q “ n BPL0 p ε { ε q α ε ă ε p ε { ε q β otherwise . (6)where ε is the photon energy, the maximum photon num-ber density is at an energy of ε and following [39] we takethe slope parameters α “ ` and β “ ´
2. As we shallsee later, any peaky spectrum gives similar results, withthe position of the peak, ε , being the most importantparameter besides the peak photon density.Inspired by the energy dependence of the diffusion co-efficient for propagation in a turbulent magnetic field, wemodel τ esc as a power law in rigidity E { Z , τ esc “ τ p EZ ´ { E q δ . (7)Since only the ratio of escape and interaction times mat-ters, and the t E, A, Z u dependence of this ratio is entirelydetermined once the spectral index of the escape time δ is specified, the remaining freedom in characterizing thesource environment can be encoded by specifying the ra-tio of escape to interaction time for a particular choiceof t E, A, Z u , which we take at 10 eV for iron nuclei,denoted R Fe19 . In application to a particular source can-didate, R Fe19 depends on the density of photons and theproperties of the turbulent magnetic field that delays theescape of the UHECRs from the environment of theirsource.Figure 3, upper panel, shows the escape and interac-tion times in the fiducial source environment, as a func-tion of the cosmic ray energy, for proton, He, N, Si andFe; the interaction times are calculated including bothphoto-disintegration and photo-pion production. Thegross features of the energy dependence of the interactiontimes can be understood in the approximation of reso-nant interactions in the nucleus rest frame ε res . At lowcosmic-ray energies, reaching ε res requires high photonenergy ( ε ą ε ), so that the interaction time decreaseswith increasing cosmic-ray energy as τ ∝ E β ` . Howeverfor high enough cosmic ray energy, the resonance can bereached in collisions with photons of ε ă ε . From here,as the cosmic ray energy increases, the photon density atthe resonant energy decreases as ε α , and correspondinglythe interaction times increase. The laboratory energy ofthe inflection point of the interaction times for a cosmicray nucleus of mass Am p is at E “ Am p ε res {p ε q . Theinflection point of the photo-dissociation times can beseen as a dip in the plot in the upper panel of Fig. 3, e.g.,at around 10 . eV for iron nuclei. At slightly higherenergy, photo-pion production becomes important, withthe result that the energy dependence of the interactiontime is roughly speaking an L-shaped curve in a log-logpresentation.Using these energy-dependent interaction and escapetimes, we propagate nuclei through the source environ-ment with the procedure described in Appendix C. Cos-mic rays of some given composition are injected from theaccelerator into the source environment with a power law lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 interactionescape = -1.00 inj g – /eV) = 18.5 pmax lg(E 0.04 – ) = 2.64 Fe19esc lg(R 0.04 – = -0.764 esc d – = 0.569 gal f 0.04 – = -4.18 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = 0 max (X sys = 0, n sys lgE D /ndf = 386.74/61 c spec: 206.149/30, lnA: 136.076/18, VLnA: 44.5156/18 evolution: SFR2, IRB: Gilmore12 f(28)= 1.0e+00 = 9.2e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 interactionescape = -1.00 inj g – /eV) = 18.5 pmax lg(E 0.04 – ) = 2.64 Fe19esc lg(R 0.04 – = -0.764 esc d – = 0.569 gal f 0.04 – = -4.18 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = 0 max (X sys = 0, n sys lgE D /ndf = 386.74/61 c spec: 206.149/30, lnA: 136.076/18, VLnA: 44.5156/18 evolution: SFR2, IRB: Gilmore12 f(28)= 1.0e+00 = 9.2e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) FIG. 3.
Top:
Interaction and escape times for A “ , , ,
28 and 56 (bottom to top for escape; vice versa forinteraction) for the fiducial model photon field with ε “ . Bottom:
Injected Si flux (bold dashed) and escapingfluxes: thin black solid line denotes the sum of all escap-ing nuclei and solid curves give contribution of different massgroups with low energy intercept increasing with mass. Nu-cleons from photo-dissociation and photo-pion-production areshown with thin-dashed and dotted curves, respectively. spectrum and an exponential cutoff at some maximumrigidity. To keep the complexity of the fiducial model toa minimum, we inject only a single nuclear species and fixthe injection spectral index γ “ ´
1, as expected for accel-eration in young neutron stars [48]. The particles escap-ing the source environment are then propagated throughthe intergalactic medium using the procedure explainedin Appendix D.In total the fiducial model has 14 parameters, with 8parameters allowed to float freely in the fit, as indicatedin Table I. The spectral index and normalization of theGalactic spectrum are free “nuisance” parameters withthe best fit giving a spectral index of ´ .
2. This shouldbe understood as an effective spectral index describingthe cutoff of the Galactic cosmic ray population, andhence cannot be directly compared with the parameterreported by the KASCADE-Grande Collaboration [49],because their single-power law fit is driven by the “low-energy” data. The fraction of Galactic cosmic rays at10 . eV is 55%.The best description of the data is obtained with Siof maximum energy Z . eV “ . ˆ eV; the im-pact of allowing other parameters to vary is discussedin following sections. Normalizing this model to the ob-served flux at Earth, we infer a comoving volumetric en-ergy injection rate in CRs at z “
0, above 10 . eV, of . (cid:15) . “ . ˆ erg Mpc ´ yr ´ .The unmodified injection spectrum and the spectrumof escaping nuclei for this fiducial model are shown inthe lower panel of Fig. 3. At low energies, the nuclei aredepleted relative to the injected flux because τ esc " τ int ,but the escaping nuclei follow the original spectral indexbecause in this example the interaction and escape timesare parallel, as to be expected for δ “ β `
1. Once thecorner of the L-shape is reached, the fraction of escapingnuclei grows, leading to an apparent hardening of thespectral index.Even for the simple case in which a single nuclearspecies is injected into the source environment, we ob-tain a complex evolution of the mass composition withenergy. At low energies the composition is dominated byknock-off nucleons whereas at high energies the composi-tion becomes heavier as the ratio of escape to interactiontime drops and more heavy nuclei can escape before in-teracting.This fiducial model of interactions in the source envi-ronment is a very simple one, yet even so it offers a re-markably good accounting for the flux and compositionat Earth as determined by the Pierre Auger Observatory.(Data from the Telescope Array (TA) are consistent withthe Auger results within systematic and statistical un-certainties [51, 52] and also can be well-fit; we come toTA separately below.) In Fig. 4 we compare the fidu-cial model prediction to the Auger measured flux, from10 . eV to above 10 eV [50] and to the mean and vari-ance of the distribution of the logarithm of mass on topof the atmosphere, x ln A y and V p ln A q [16, 53, 54]. Thereis a good overall agreement between the model and thedata. The shape of the spectrum is described well, in-cluding the ankle and the flux suppression. The modelalso qualitatively reproduces the increase of the averagelogarithmic mass with energy and the decrease of its vari-ance.The neutrino signals of the fiducial model are shownin Fig. 5; details of the calculation are given in AppendixE. An exciting aspect of our model for the ankle is thepresence of a detectable anti-electron-neutrino flux fromneutron β -decay, with a rate consistent with the na¨ıve es-timate of [57]. The right panel of Fig. 5 shows the numberof events as a function of energy predicted in ten years ofIC86, using the IceCube acceptance for different neutrinoflavors given in [56]. In total, the fiducial model predicts3.5 events in the range 10 ´ eV after 10 years ofoperation of IceCube (corresponding to about one year ofoperation for an upgraded IceCube-Gen2 detector [58]).We emphasize the distinctive ¯ ν e enrichment due to betadecay of spallated neutrons.The associated photon flux from nuclear de-excitationin our model is well below the Fermi-LAT data (seeAppendix E for more detail). Photo-pion interactionsat the source and during propagation produce an addi- tional flux of photons via π -decay; this is consistent withFermi-LAT data, as follows: If the origin of the photonsmeasured by Fermi-LAT is exclusively from these inter-actions, then from [59] the associated diffuse neutrinoflux saturates the IceCube upper limit [56]. Since theneutrino flux in the fiducial model is below the IceCubelimit, it follows that also the associated photon flux isconsistent with Fermi-LAT data. A more sophisticatedrealization of our mechanism than in the fiducial modelmust also respect the IceCube limits, and therefore theFermi-LAT data as well. B. Model Variations
In this section we discuss the impact of theoretical andexperimental uncertainties on our model, as well as dif-ferent choices for the fiducial parameters.
1. Experimental Uncertainties
To study the influence of the experimental systematicuncertainties on our fit, we have repeated the fit for allcombinations of altering the measurements by ` ` ´ σ sys . of the quoted uncertainties on the energyand composition scale. We find that the best fit is ob-tained within the experimental systematics when shiftingthe energy scale up by ` σ sys . “ `
15% and by shifting x ln A y and V p ln A q corresponding to a shift of the showermaximum by ´ σ sys . « ´
10 g { cm . The best-fit valuesafter the application of these shifts are shown in bracketsin Table I. Most notably, the peak energy of the photonspectrum decreases from 110 to 70 meV and the best-fitvalue of the spectral index of the escape time decreasesfrom „ ´ { ´
1. The neutrino flux at Earthobtained for this fit is about 30% smaller than in case ofthe fiducial model. This is mainly due to the difference inthe best-fit peak energy of the photon field in the sourceenvironment. The sensitivity of the neutrino flux to ε will be further discussed in Sec. III B 5.The overall description of the spectrum and composi-tion is considerably improved, as can be seen in Fig. 6.The model variations discussed below will therefore beperformed based on shifted data.
2. Hadronic Interactions in Air Showers
The interpretation of experimental air shower data interms of mass composition relies on the validity of extrap-olations of the properties of hadronic interactions to ul-trahigh energies. Using alternative models for this inter-pretation (
Sibyll2.1 [62] or
QGSJetII-04 [63] insteadof
Epos-LHC [64]), decreases the value of the x ln A y datapoints by about x ln A y “ ´ . lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 interactionescape = -1.00 inj g – /eV) = 18.5 pmax lg(E 0.04 – ) = 2.64 Fe19esc lg(R 0.04 – = -0.764 esc d – = 0.569 gal f 0.04 – = -4.18 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = 0 max (X sys = 0, n sys lgE D /ndf = 386.74/61 c spec: 206.149/30, lnA: 136.076/18, VLnA: 44.5156/18 evolution: SFR2, IRB: Gilmore12 f(28)= 1.0e+00 = 9.2e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (a) Flux at Earth lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 interactionescape = -1.00 inj g – /eV) = 18.5 pmax lg(E 0.04 – ) = 2.64 Fe19esc lg(R 0.04 – = -0.764 esc d – = 0.569 gal f 0.04 – = -4.18 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = 0 max (X sys = 0, n sys lgE D /ndf = 386.74/61 c spec: 206.149/30, lnA: 136.076/18, VLnA: 44.5156/18 evolution: SFR2, IRB: Gilmore12 f(28)= 1.0e+00 = 9.2e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (b) Composition at Earth FIG. 4. Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50], error bars denotethe statistical uncertainties and the shaded boxes illustrate the experimental systematic uncertainties of the composition. Thecomposition estimates are based on an interpretation of air shower data with
Epos-LHC ; the lines denote the predictions ofour fiducial model. lg(E/eV)
13 14 15 16 17 18 19 ] - s - s r - ( E ) [ G e V c m F E -12 -11 -10 -9 -8 -7 e n m n t n sum IC2013 e n m n t n IC2014 lg(E/eV) e v en t s / . l g ( E ) / I C - y ea r s e n + e n m n + m n t n + t n sum events = 3.5 S lg(E/eV)
13 14 15 16 17 18 19 ] - s - s r - ( E ) [ G e V c m F E -12 -11 -10 -9 -8 -7 e n m n t n sum IC2013 e n m n t n IC2014 lg(E/eV) e v en t s / . l g ( E ) / I C - y ea r s e n + e n m n + m n t n + t n sum events = 3.5 S FIG. 5. Neutrino spectrum (left) and expected number of events in 10 IC86-years (right) for the fiducial model. The measuredflux of low-energy extragalactic neutrinos from IceCube [55] is shown in the left panel (purple lines) as well as the 90% CLupper limit on the flux of high-energy neutrinos (dashed area) [56]. The peak in the electron neutrino flux at about 10 . eVseen in the right panel is due to the increased interaction probability of anti-electron neutrinos at the Glashow resonance. in both directions, σ theo px ln A yq “ ˘ .
6, then a hadronicinteraction model that leads to a heavier interpretationof Auger data than
Epos-LHC would make the fit withthe fiducial model even better, similar to the systematicshift in the composition scale discussed in the previoussection.
3. Mass Composition at the Source
It is remarkable that a good description of both thespectrum and mass composition at Earth is possible by assuming only a single injected species at the source asassumed for simplicity in the fiducial model. However,depending on the astrophysical scenario, this might bean unrealistic assumption.In Fig. 7 we explore the capability of our model to in-corporate additional flux components of mass A belowand above the mass A „
29 that gives the best fit for thefiducial single-mass model. As can be seen, our calcula-tion allows for an additional proton or helium componentas large as 80% and up to 70% for nitrogen.For an additional flux component with a heavy mass,the model is more restrictive as illustrated in the lower source parameters power law index of injected nuclei γ fix ´ A free 28 (29)maximum energy E p max free 10 . p . q eVcosmic ray power density, E ą . eV . (cid:15) . free 9.2 (13) ˆ erg Mpc ´ yr ´ evolution ξ p z p t qq fix star formation rate [60] source environment energy of maximum of photon field density ε free 0 . p . q eVpower law index of photon spectrum ( ε ă ε ) α fix ` power law index of photon spectrum ( ε ě ε ) β fix ´ δ free ´ .
77 ( ´ . R Fe19 free 4.4 (3.7) ˆ propagation to Earth infra-red photon background – fix Gilmore12 [61] spectrum of Galactic cosmic rays power law index at Earth γ gal free ´ . ´ . A gal fix 56flux fraction at 10 . eV f gal free 57 (72) %TABLE I. Parameters of the fiducial model. Values in brackets denote the parameters of the best-fit obtained when shiftingthe data by its systematic uncertainties (see text). left panel of Fig. 7 using A “
56. In this case, the de-scription of the data considerably deteriorates for frac-tions above 10%. The reason for this behavior is twofold.Firstly, the injection of too much iron at the source leadsto a too heavy composition at Earth as compared to theestimates from the Pierre Auger Observatory. Secondly,if the end of the cosmic ray spectrum is to be describedby the maximum rigidity of iron nuclei, then the energyof secondary nucleons needed to populate the flux at andbelow the ankle is too small to describe the data (themaximum energy of secondary nucleons is 1 { A of themaximum energy of nuclei).If the cut-off of the flux is at higher energies, as sug-gested by the measurement of TA [65], then a larger frac-tion of iron primaries at the source can be incorporated,as shown in the lower right panel of Fig. 7. When usingthe TA data in the fit, as shown in Fig. 8, the spectrumcan be described reasonably well even for an injectedflux consisting of 100% iron nuclei. But in this scenariothe composition at Earth at ultrahigh energies is heavierthan suggested by the interpretation of the X max data ofAuger.As an illustration of a more complex compositionmodel, we use the abundances of Galactic nuclei at anucleus energy of 1 TeV, which we read from Fig. 28.1in [66]. The flux fractions are 0.365, 0.309, 0.044, 0.077,0.019, 0.039, 0.039, 0.0096, 0.014, 0.084 for H, He, C,O,Ne, Mg, Si, S, Ar+Ca, Fe, respectively. The result-ing fit is shown in Fig. 9 ( γ “ ´ .
25 and δ “ ´
4. Source Evolution and Spectral Index
To have a concrete fiducial model, we needed to specifyhow the production of UHECRs varied over cosmologicaltime scales. This is known as the source evolution, whichwe took to be in direct proportion to the star-formation-rate – as would be expected in a source scenario such asyoung magnetars. In this section, we consider alternativeevolutions of the source luminosity density described bythe simple one-parameter functional form ξ p z q “ p ` z q m z ă z p ` z q m exp p´p z ´ z qq otherwise (8)with z “ m ranging from ´ ` m “ m “ ` γ of the injected flux: ´
1, as in thefiducial model, ´ γ float freely in the fit. As can beseen in Fig. 10(a), γ “ ´ m Á
0, but is a viable choice for closeby sources,in accordance to the findings of [67]. For positive valuesof m , a fixed value of γ “ ´ γ , but the latter converges to val-ues larger than ´ m ą m ě γ “ ´ m , with the exception lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape = -1.00 inj g – /eV) = 18.6 pmax lg(E 0.05 – ) = 2.57 Fe19esc lg(R 0.05 – = -0.942 esc d – = 0.718 gal f 0.02 – = -3.72 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.07 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 174.85/61 c spec: 89.615/30, lnA: 35.6525/18, VLnA: 49.5829/18 evolution: SFR2, IRB: Gilmore12 f(29)= 1.0e+00 = 1.3e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (a) τ int and τ esc . lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape = -1.00 inj g – /eV) = 18.6 pmax lg(E 0.05 – ) = 2.57 Fe19esc lg(R 0.05 – = -0.942 esc d – = 0.718 gal f 0.02 – = -3.72 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.07 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 174.85/61 c spec: 89.615/30, lnA: 35.6525/18, VLnA: 49.5829/18 evolution: SFR2, IRB: Gilmore12 f(29)= 1.0e+00 = 1.3e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (b) Injected (dashed line) and escaping (solid lines) fluxes. lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape = -1.00 inj g – /eV) = 18.6 pmax lg(E 0.05 – ) = 2.57 Fe19esc lg(R 0.05 – = -0.942 esc d – = 0.718 gal f 0.02 – = -3.72 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.07 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 174.85/61 c spec: 89.615/30, lnA: 35.6525/18, VLnA: 49.5829/18 evolution: SFR2, IRB: Gilmore12 f(29)= 1.0e+00 = 1.3e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (c) Flux at Earth lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -2 -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape = -1.00 inj g – /eV) = 18.6 pmax lg(E 0.05 – ) = 2.57 Fe19esc lg(R 0.05 – = -0.942 esc d – = 0.718 gal f 0.02 – = -3.72 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.07 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 174.85/61 c spec: 89.615/30, lnA: 35.6525/18, VLnA: 49.5829/18 evolution: SFR2, IRB: Gilmore12 f(29)= 1.0e+00 = 1.3e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (d) Composition at Earth FIG. 6. Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50] shifted by plusone sigma of systematic uncertainty for the energy scale and minus one sigma for the X max scale. The lines denote the best-fitwithin our fiducial model. of the power-law index of the escape time δ (Fig. 10(e))and the power density . (cid:15) . (Fig. 10(e)).We conclude that our model for the ankle does notcritically depend on the choice of the source evolution,but that for a given choice of m we can constrain theallowed values of γ , δ and . (cid:15) . .
5. Photon Spectrum
We repeated the model fits using alternative energydistributions of the photon density instead of the brokenpower law used in the fiducial model: a black body spec-trum and two modified black body spectra. All four spec-tra are normalized to the same integral photon densityand depend only on one parameter, the peak energy ε (see Appendix A). The resulting fit results are shown inFig. 11 for a freely floating spectral index γ and for sourceevolutions with m ě
0. As can be seen, all four photonspectra describe the data equally well (Fig. 11(a)). The best-fit values of the free model parameters are very sim-ilar and in particular the obtained peak values are within ˘
20 meV. We conclude that as long as the photon spec-trum is “peaky”, the particular details of its shape donot influence the parameters of our model.The sensitivity of the fit to the peak energy is shown inthe left panel of Fig. 12. As can be seen, the χ deterio-rates very quickly at low values of ε , but it is almost flatabove the minimum. This feature can be easily under-stood recalling ε in the “L-curve” approximation intro-duced in Sec. II: The smaller ε , the larger is the energyof inflection of the interaction length, E b . For too-smallvalues of ε , the interaction and escape times are parallelover the full energy range and thus no high-pass filter iscreated. On the other hand, once E b is small enough,a further decrease changes only the flux at low energy,where the escaping spectrum is dominated by low-massnuclei from spallation (see e.g. Fig. 3) which can be com-pensated by adjusting other parameters such as R Fe19 . = 1 flux fraction of A A / N d f, ) c m i n ( (a) Proton. = 4 flux fraction of A A / N d f, ) c m i n ( (b) Helium. = 14 flux fraction of A A / N d f, ) c m i n ( (c) Nitrogen. = 56 flux fraction of A A / N d f, ) c m i n ( (d) Iron, Auger flux. = 56 flux fraction of A A / N d f, ) c m i n ( (e) Iron, TA flux. FIG. 7. Injection of two mass components. The first mass value, A , is fixed and contributes the fraction indicated on thex-axis to the total flux. The second mass value, A , is varied as shown on the y-axis. The fit quality is indicated by the colors. To first order, our model can therefore only give a lowerlimit on the peak energy of the photon flux in the sourceenvironment. However, future limits or observations ofneutrinos in the 10-100 PeV range will help to constrainthis important source property, because the number ofpredicted neutrinos strongly depends on ε , as shown inthe left panel of Fig. 12 by the superimposed open sym-bols. A larger peak energy of the ambient photon en-vironment increases neutrino production at the sourcein two ways. Firstly, shifting E b to lower energies (andcompensating as necessary by adjustment of R ) movesthe interaction times of protons closer to the escape timeand correspondingly additional neutrinos are producedvia photo-pion production of protons (compare e.g. thered curves at around 10 eV in the upper panel of Fig. 3( ε “
110 meV) to the ones in Fig. 6(a) ( ε “
70 meV)).Secondly, increasing ε moves the minimum of the inter-action time for photo-pion production of nuclei to lowerenergies. Since the neutrinos from photo-pion productioncarry a larger fraction of the nucleon energy than the neu-trinos from neutron decay after photo-dissociation, thisincreases the neutrino flux as well. It is tempting to give a quantitative interpretationof the χ -curve of Fig. 12 in terms of a lower limit on ε and the number of neutrinos. However, the mini-mum of χ is far away from χ { N df “ S “ p χ { N df q to bring the rescaled χ { N df to 1. This rescales the χ value of any given model sothat the number of standard-deviations it is from theminimum is given by N σ “ S ´ a χ ´ χ . Thisyields an approximate lower limit on ε at N σ “ ε ą
34 meV and N ν p ˆ IC86 q ą . σ “ lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV)
18 18.5 19 19.5 20 20.5 ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
TA 2013 lg(E/eV) [ a . u .] t c -9 -8 -7 -6 interactionescape = -1.00 inj g – /eV) = 18.7 pmax lg(E 0.08 – ) = 2.38 Fe19esc lg(R 0.09 – = -0.489 esc d – = 0.768 gal f 0.07 – = -3.79 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = -1 max (X sys = 0, n sys lgE D /ndf = 220.705/55 c spec: 42.2301/24, lnA: 119.446/18, VLnA: 59.0293/18 evolution: SFR2, IRB: Gilmore12 f(56)= 1.0e+00 = 7.4e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (a) τ int and τ esc . lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV)
18 18.5 19 19.5 20 20.5 ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
TA 2013 lg(E/eV) [ a . u .] t c -9 -8 -7 -6 interactionescape = -1.00 inj g – /eV) = 18.7 pmax lg(E 0.08 – ) = 2.38 Fe19esc lg(R 0.09 – = -0.489 esc d – = 0.768 gal f 0.07 – = -3.79 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = -1 max (X sys = 0, n sys lgE D /ndf = 220.705/55 c spec: 42.2301/24, lnA: 119.446/18, VLnA: 59.0293/18 evolution: SFR2, IRB: Gilmore12 f(56)= 1.0e+00 = 7.4e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (b) Injected (dashed line) and escaping (solid lines) fluxes. lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV)
18 18.5 19 19.5 20 20.5 ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
TA 2013 lg(E/eV) [ a . u .] t c -9 -8 -7 -6 interactionescape = -1.00 inj g – /eV) = 18.7 pmax lg(E 0.08 – ) = 2.38 Fe19esc lg(R 0.09 – = -0.489 esc d – = 0.768 gal f 0.07 – = -3.79 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = -1 max (X sys = 0, n sys lgE D /ndf = 220.705/55 c spec: 42.2301/24, lnA: 119.446/18, VLnA: 59.0293/18 evolution: SFR2, IRB: Gilmore12 f(56)= 1.0e+00 = 7.4e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (c) Flux at Earth lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV)
18 18.5 19 19.5 20 20.5 ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
TA 2013 lg(E/eV) [ a . u .] t c -9 -8 -7 -6 interactionescape = -1.00 inj g – /eV) = 18.7 pmax lg(E 0.08 – ) = 2.38 Fe19esc lg(R 0.09 – = -0.489 esc d – = 0.768 gal f 0.07 – = -3.79 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.11 eV, e s ) = -1 max (X sys = 0, n sys lgE D /ndf = 220.705/55 c spec: 42.2301/24, lnA: 119.446/18, VLnA: 59.0293/18 evolution: SFR2, IRB: Gilmore12 f(56)= 1.0e+00 = 7.4e+44 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (d) Composition at Earth FIG. 8. Spectrum and composition at Earth. The data points are from the TA [65] (flux) and the Pierre Auger Observatory [16](composition). The latter have been shifted in energy to match the energy scale of TA and the X max scale is shifted down by1 sigma. The lines denote the fit with our model assuming a pure iron composition at the source.
6. Hadronic Interactions in the Source Environment
In addition to interactions with the background pho-ton field, nucleons and nuclei can also scatter off hadronsin the source environment. In this paper we assumethat the density of hadronic matter in the source en-vironment is low enough that such hadronic interactionscan be neglected. For any concrete astrophysical real-ization of our scenario, one must check and if necessaryinclude hadronic interactions in the source environment.Production of π ˘ ’s and π ’s in hadronic collisions couldsignificantly increase the fluxes of neutrinos and photonsemitted in the EeV energy range. Fast-spinning newbornneutron stars provide a particular example [69]. Pre-cise estimates of the impact of hadronic collisions on thepredictions of our model will be presented in a separatepublication. The results presented here are valid for allastrophysical systems in which the interactions are dom-inated by photo-nuclear processes. IV. CONCLUSIONS
In this paper we have proposed a new explanation forthe ankle in the cosmic ray spectrum, and for the evo-lution with energy of the composition of extragalacticcosmic rays: from light below the ankle to increasinglyheavy above. When nuclei are trapped in the turbulentmagnetic field of the source environment, their escapetime can decrease faster with increasing energy than doestheir interaction time. Under these conditions, only thehighest energy particles can escape the source environ-ment unscathed, and the source environment acts as ahigh-pass filter on UHECRs. Nuclei below the crossoverenergy such that τ esc ą τ int interact with photons in theenvironment around the source, with ejection of nucle-ons or alpha particles and consequent production of asteep spectrum of secondary nucleons. The superpositionof this steeply falling nucleon spectrum with the harderspectrum of the surviving nuclear fragments creates anankle-like feature in the total source emission spectrum.1 lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape – = -1.25 inj g – /eV) = 18.6 pmax lg(E 0.02 – ) = 2.57 Fe19esc lg(R 0.03 – = -1.01 esc d – = 0.686 gal f 0.02 – = -3.71 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.09 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 244.504/60 c spec: 120.825/30, lnA: 84.6523/18, VLnA: 39.0268/18 evolution: SFR2, IRB: Gilmore12 f(1) = 3.6e-01f(4) = 3.1e-01f(12)= 4.5e-02f(16)= 7.7e-02f(20)= 1.9e-02f(24)= 3.9e-02f(28)= 3.9e-02f(32)= 9.6e-03f(40)= 1.4e-02f(56)= 8.4e-02 = 1.5e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (a) Injected fluxes in 5 mass groups. lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape – = -1.25 inj g – /eV) = 18.6 pmax lg(E 0.02 – ) = 2.57 Fe19esc lg(R 0.03 – = -1.01 esc d – = 0.686 gal f 0.02 – = -3.71 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.09 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 244.504/60 c spec: 120.825/30, lnA: 84.6523/18, VLnA: 39.0268/18 evolution: SFR2, IRB: Gilmore12 f(1) = 3.6e-01f(4) = 3.1e-01f(12)= 4.5e-02f(16)= 7.7e-02f(20)= 1.9e-02f(24)= 3.9e-02f(28)= 3.9e-02f(32)= 9.6e-03f(40)= 1.4e-02f(56)= 8.4e-02 = 1.5e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (b) Escaping fluxes (sum of injection shown as dashed line). lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape – = -1.25 inj g – /eV) = 18.6 pmax lg(E 0.02 – ) = 2.57 Fe19esc lg(R 0.03 – = -1.01 esc d – = 0.686 gal f 0.02 – = -3.71 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.09 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 244.504/60 c spec: 120.825/30, lnA: 84.6523/18, VLnA: 39.0268/18 evolution: SFR2, IRB: Gilmore12 f(1) = 3.6e-01f(4) = 3.1e-01f(12)= 4.5e-02f(16)= 7.7e-02f(20)= 1.9e-02f(24)= 3.9e-02f(28)= 3.9e-02f(32)= 9.6e-03f(40)= 1.4e-02f(56)= 8.4e-02 = 1.5e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (c) Flux at Earth lg(E/eV) d N / d l g E / d t [ a . u .] n -1 injected lg(E/eV) d N / d l g E / d t [ a . u .] n -1
10 110 injected £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £ lg(E/eV) ] - y r - s r - k m J ( E ) [ e V E · £ A £
1 6 £ A £
3 19 £ A £
7 39 £ A £
20 56 £ A £
40 galactic (A=56)
Auger 2013 prel. lg(E/eV) [ a . u .] t c -9 -8 -7 -6 -5 -4 interactionescape – = -1.25 inj g – /eV) = 18.6 pmax lg(E 0.02 – ) = 2.57 Fe19esc lg(R 0.03 – = -1.01 esc d – = 0.686 gal f 0.02 – = -3.71 gal g lg(fphot) = 0.00 =-2 b =1.5, a = 0.09 eV, e s ) = -1 max (X sys = +0.1, n sys lgE D /ndf = 244.504/60 c spec: 120.825/30, lnA: 84.6523/18, VLnA: 39.0268/18 evolution: SFR2, IRB: Gilmore12 f(1) = 3.6e-01f(4) = 3.1e-01f(12)= 4.5e-02f(16)= 7.7e-02f(20)= 1.9e-02f(24)= 3.9e-02f(28)= 3.9e-02f(32)= 9.6e-03f(40)= 1.4e-02f(56)= 8.4e-02 = 1.5e+45 e yr Mpcerg lg(E/eV)
18 18.5 19 19.5 æ l n A Æ lg(E/eV)
18 18.5 19 19.5 V ( l n A ) (d) Composition at Earth FIG. 9. Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50] shifted by theirsystematic uncertainty as in Fig. 6. The injected composition follows a Galactic mixture with 10 elements (see text).
Above the ankle, the spectrum emerging from the sourceenvironment exhibits a progressive transition to heaviernuclei, as the escape of non-interacting nuclei becomesefficient. Abundant production of ¯ ν e ’s is a signature ofthis mechanism.We illustrated the high quality of the fit which canbe obtained to the Auger data, with a fiducial modelin which nuclei are accelerated up to a maximum rigid-ity found to be « . V, with spectrum E ´ , andare then subject to photo-disintegration in the vicinityof the accelerator before escaping for their journey toEarth. We showed that the details of the photon spec-trum around the accelerator are unimportant, except forits peak energy. The other important characteristic ofthe environment is the photon density relative to themagnetic diffusivity, which we characterized in a verysimplistic way (through a single parameter) in this ini-tial study. We studied the sensitivity of the mechanismto the energy-scale uncertainty and hadronic-interaction-modeling uncertainty, which affects the composition in- ferred from the atmospheric shower observations, andalso used the TA spectrum instead of the Auger spec-trum. The conclusion of these studies is that a goodquality fit can be obtained in most cases, but details ofthe fit parameters such as the composition and maximumenergy characterizing the accelerator change. A corollaryis that until these systematic uncertainties in the obser-vations and their interpretation are reduced, such detailsof the accelerator cannot be reliably inferred from thedata. The fiducial model parameters needed in the fitsare such that the scenario can be reasonably achieved inat least one type of proposed astrophysical source, as willbe discussed in a future publication.Our mechanism has two predictions beyond fitting theshape of the spectrum and composition evolution, whichare independent of many environmental variables and canbe used to test the validity of this scenario for produc-tion of the ankle. i) The spectral cutoff of spallated nu-cleons emerging from the source environment is R max ,where R max is the rigidity cutoff of the accelerator, be-2 evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR / nd f c =-1 g =-2 g =free g (a) fit quality. evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR [ e V ] e (b) peak energy. evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR g s pe c t r a l i nde x -2-1.5-1-0.50 (c) spectral index of injected spectrum. evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR ) F e 19 l g ( R (d) ratio of escape and interaction time. evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR e sc ape d -1-0.9-0.8-0.7-0.6-0.5-0.4-0.3 (e) power-law index of escape time. evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR p ) m a x l g ( E (f) maximum energy for Z “ evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR m a ss (g) injected mass evolution -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 SFR )) - y r - / ( e r g M p c e l g ( (h) UHECR energy injection rate FIG. 10. Fit results as a function of source evolution for different spectral indices of the injected flux: γ fixed to ´ ´ m of the source evolution is shown; thelast bin reports the values for the fiducial model (SFR) evolution from [60], Eq. (D6). evolution / nd f c BPLBB =1 s MBB, =2 s MBB, (a) fit quality. evolution [ e V ] e (b) peak energy. evolution g s pe c t r a l i nde x -2-1.5-1-0.500.5 (c) spectral index of injected spectrum. evolution ) F e 19 l g ( R (d) ratio of escape and interaction time. evolution e sc ape d -1-0.9-0.8-0.7-0.6-0.5-0.4-0.3 (e) power law index of escape time. evolution p ) m a x l g ( E (f) maximum energy for Z “ evolution m a ss (g) injected mass evolution )) - y r - / ( e r g M p c e l g ( (h) UHECR energy injection rate FIG. 11. Fit results as a function of source evolution for different photon spectra: Broken power law (BPL, open squares),black body spectrum (BB, open circles), modified black body spectrum (MBB) with σ “ σ “ m of the source evolution is shown and in the last bin the fit values for the fiducial evolutionfrom [60], Eq. (D6), is shown. [meV] e / nd f c nu m be r o f neu t r i no e v en t s i n I C y ea r s fit quality number of neutrinos BPL BB =1 s MBB =2 s MBB lg(E/eV)
13 14 15 16 17 18 19 ] - s - s r - ( E ) [ G e V c m F E -12 -11 -10 -9 -8 -7 e n m n t n sum IC2013 e n m n t n IC2014 lg(E/eV) e v en t s / . l g ( E ) / I C - y ea r s e n + e n m n + m n t n + t n sum events = 0.4 S FIG. 12. Left: Fit quality of the fiducial model (closed symbols) and number of neutrinos (open symbols) as a function ofpeak energy ε of the photon spectrum in the source environment. Four types of photon spectra are shown: Broken power law(BPL), black body spectrum (BB) and two modified black body spectra (MBB). The minimum χ of BPL corresponds to theresult shown in Fig. 6. Right: Lower limit on the neutrino flux obtained for a modified black body spectrum with σ “
2, and ε “
34 meV ( T “
100 K). The lines and hatched area at the top of the figure are the measured neutrino flux and upper limitfrom IceCube [55, 56] cause E max , spal . nuc . “ E max ,A { A while E max ,A “ Z R max ,and finally Z { A “ , largely independent of composition.This relation holds prior to the extragalactic propaga-tion from the source, thus giving complementary infor-mation on the accelerator to that obtained from the spec-trum and composition above the ankle alone. ii) Thereis a one-to-one relation between the spectrum of spal-lated nucleons and the anti-electron-neutrinos producedby beta decay of neutrons, unless the spallated nucle-ons lose energy by interacting with hadronic material inthe source environment. Independent of other proper-ties of the environment or the source evolution, ¯ ν e ’s willhave an identical spectral shape, shifted down by a factor „ { n Ñ p e ´ ¯ ν e and reducedby a factor-2 in normalization because only half the nu-cleons are neutrons. This follows because propagationenergy losses are small for nucleons of such low energy,and redshift impacts both nucleons and neutrinos identi-cally. Thus, detailed comparison of the ¯ ν e and spallatednucleon spectra will reveal if hadronic interactions in thesource environment are important, which would imply acorrelated production of photo-pion produced neutrinos. NOTE ADDED
After this work was presented at the IceCube ParticleAstrophysics Symposium a paper appeared on the arXiv exploring another mechanism for producing the ankle,arising in the context of gamma-ray bursts [70].
ACKNOWLEDGMENTS
We would like to acknowledge many useful discussionswith our colleagues of the Pierre Auger Collaboration.Furthermore we thank David Walz for his support regard-ing questions about
CRPropa and Benoit Marchand forhis help running the calculations on the BuTinah highperformance computing cluster of NYU Abu Dhabi. MUacknowledges the financial support from the EU-fundedMarie Curie Outgoing Fellowship, Grant PIOF-GA-2013-624803. The research of GRF is supported in part by theU.S. National Science Foundation (NSF), Grant PHY-1212538 and the James Simons Foundation; she thanksKIPAC/SLAC for their hospitality. The research of LAAis supported by NSF (Grant CAREER PHY-1053663)and NASA (Grant NNX13AH52G); he thanks the Centerfor Cosmology and Particle Physics at New York Univer-sity for its hospitality.5 † [eV] d n / d † [ e V − c m − ] BB, T =364 KMBB, T =206 K, σ =1 MBB, T =148 K, σ =2 BPL, † =0 . eV, α =3 / , β = − . lg(E/eV) [ y r ] t c -1 FIG. 13. Left: Comparison of photon spectra. BPL: Broken power law (solid), BB: black body spectrum (dashed), MBB:modified black body spectrum (dotted and dash-dotted). The curves are normalized to match the integral of the black bodyspectrum and the temperatures are chosen to match the peak energy of the broken power law. Right: Interaction timescorresponding to the four photon spectra.
Appendix A: Photon spectra
In this paper we explore the propagation effects of the four types of photon spectra shown in Fig. 13. The firstconsists of broken power-law (see e.g. [39]) as a simplified representative of non-thermal emission given by n p ε q “ n BPL0 p ε { ε q α ε ă ε p ε { ε q β otherwise . (A1)where ε is the photon energy and the maximum of the number density is at an energy of ε .We also consider modified black-body spectra using the functional form n p ε q “ n MBB0 π p hc q ε e εkT ´ ˆ εε ˙ σ (A2)where T denotes the temperature in the case of pure black-body, and the absorption factor is given by p εε q σ (seee.g. [42]). h , k and c are the Planck constant, Boltzmann constant and speed of light respectively. For σ “ n MBB0 “ ε “ “ W ` ´ e ´ b b ˘ ` b ‰ k T, (A3)where W p x q is the Lambert function (see e.g. [71]) and b “ σ ` I BPL “ n BPL0 ε ˆ α ` ´ β ` ˙ (A4)and for Eq. (A2) it is I σ MBB “ n MBB0 π p hc q p kT q ζ p σ ` , q Γ p σ ` q , (A5)where ζ p x q denotes the Riemann zeta function and Γ p x q is the Gamma function. Choosing the photon density of theunmodified black body spectrum as reference we use the following normalization constants, n BPL0 “ I { I BPL (A6)6and n MBB0 p σ q “ I { I σ MBB . (A7)An example of the four photon spectra after normalization and for the same peak energy of ε “
50 meV is shownin the left panel of Fig. 13. The corresponding interaction time for the sum of photo-dissociation and photo-pionproduction is shown in the right panel.
Appendix B: Photo-Nuclear Interactions
The interaction between photons and high energy nuclei has been extensively discussed in the literature [72–82]. Inthis appendix, we describe how we implement the photon-nucleus collisions in our analysis. The interaction time fora highly relativistic nucleus with energy E “ γAm p (where γ is the Lorentz factor) propagating through an isotropicphoton background with energy ε and spectrum n p ε q , normalized so that the total number density of photons is ş n p ε q dε , is given by [72] 1 τ int “ c ż dε n p ε q γ ε ż γε dε ε σ p ε q , (B1)where σ p ε q is the photo-nuclear interaction cross section of a nucleus of mass Am p by a photon of energy ε in therest frame of the nucleus.Detailed tables of σ p ε q are available in CRPropa [83, 84]. We use the numerical tools provided at [47] to calculatethe interaction times for the photon field given by Eqs. (A1) and (A2).For illustrative purposes, the cross section can be approximated by a single pole in the narrow-width approximation, σ p ε q “ π σ res Γ res δ p ε ´ ε res q , (B2)where σ res is the resonance peak, Γ res its width, and ε res the pole in the rest frame of the nucleus. The factor of 1 { τ int p E q « c π σ res ε res Γ res γ ż dεε n p ε q Θ p γε ´ ε res q“ c π σ res ε res Γ res γ ż (cid:15) res { γ dεε n p ε q . (B3)Substituting (A1) into (B3) yields:1 τ int p E q “ τ b p E b { E q β ` E ď E b p ´ β q{p ´ α q ” p E b { E q α ` ´ p E b { E q ı ` p E b { E q E ą E b , (B4)where τ b “ E b p ´ β q c π σ res A m p Γ res n and E b “ ε res A m p ε . (B5)The parameters characterizing the photo-disintegration cross section are: σ res « . ˆ ´ cm A , Γ res “ (cid:15) res “ . A ´ . p . A . q MeV , for A ą A ď
4) [74]. The parameters for the photo-pion productioncross section are: σ res » . ˆ ´ cm A , Γ res “
150 MeV, and ε res “ p m ´ m p q{p m p q »
340 MeV [66].
Appendix C: Propagation in the Source Environment
In our simple model we consider interactions and escape of particles treating the source environment as a leakybox. If at a given time, t , N p t “ t q “ N particles are injected at random into the source environment, then thenumber of particles N remaining in the source at any later time t changes as dNdt “ ´ τ esc N ´ τ int N, (C1)7where τ esc and τ int are the escape and interaction times respectively. Integration yields the time evolution of N as N p t q “ N e ´ t ´ t τ , (C2)where τ “ τ esc τ int τ esc ` τ int . (C3)The total number of escaping particles is given by N esc “ ż t τ esc N p t q dt “ N τ int τ esc ` τ int “ N ` τ esc { τ int ” N f esc (C4)and likewise the number of particles suffering interactions is N int “ N τ esc τ esc ` τ int “ N ` τ int { τ esc ” N f int , (C5)with N esc ` N int “ N . As can be seen, N esc and N int depend only on the ratio of the escape and interaction times,but not on the absolute value of either of them.In the following we consider sources at steady state, i.e. sources which are active long enough to justify integratingto infinity in Eq. (C4) and for which the injected flux equals the escaping flux. Interacting particles constitute thesource for secondary particles of lower mass number.Since the particle trajectory in the source is treated as a random walk starting from a random position, the escapetime of a secondary does not depend on the time it was produced. Therefore we can apply Eqs. (C1) to (C5) also tothe secondary particle production, which greatly simplifies the equations with respect to previous analytic approachesthat had been developed for the extra-galactic propagation of cosmic-ray nuclei [79–82].
1. Single-Nucleon Emission
The basic principle of the analytic calculation can be best illustrated by firstly describing the case where interactionswith the photon field lead to the knock-out of a single nucleon, A ` γ Ñ p A ´ q ` n { p, (C6)and the nucleon carries away a fraction of 1 { A of the initial energy of the nucleus. This approach has been successfullyapplied to the photo-disintegration (PD) during the extra-galactic propagation of nuclei (see e.g. [79]). It can also serveas a good approximation for the losses due to photo-pion production (PP) if nuclei are treated as the superposition of A individual nucleons (see e.g. [83, 84]). The interaction time is therefore the combination of the two processes, i.e. τ int “ τ PDint ` τ PPint τ PDint τ PPint . (C7)In this simplified propagation scheme, secondaries with mass A and energy E ˚ originate from nuclei with energy E “ A ` A E ˚ and mass A `
1. They are produced at a rate Q p E ˚ , A q “ Q int ˆ A ` A E ˚ , A ` ˙ ˇˇˇˇ dE dE ˚ ˇˇˇˇ “ Q ˆ A ` A E ˚ , A ` ˙ η int ˆ A ` A E ˚ , A ` ˙ A ` A , (C8)where the factor ˇˇˇ dE dE ˚ ˇˇˇ is the Jacobian determinant needed to transform the differential injection rate from the primaryto secondary energy. In analogy to Eq. (C4), a fraction of the secondaries escapes the source environment, Q esc p E ˚ , A q “ Q p E ˚ , A q η esc p E ˚ , A q , (C9)and the remaining particles interact again at a rate of Q int p E ˚ , A q “ Q p E ˚ , A q η int p E ˚ , A q . (C10)8This assumes that the escape probability of a secondary is independent of the time or position it got produced inthe source environment. This calculation can be iterated to obtain the escape rate of any remnant with mass A ˚ produced during the propagation of a nucleus of mass A : Q remesc p E ˚ , A ˚ , A q “ Q ˆ A A ˚ E ˚ , A ˙ A A ˚ η esc p E ˚ , A ˚ q A ź A ˛ “ A ˚ ` η int ˆ A ˛ A ˚ E ˚ , A ˛ ˙ . (C11)The rate of nucleons being knocked out of nuclei during propagation via either of the considered processes i “ PD { PPis Q i esc p E ˚ , n ` p, A q “ Q ˆ A E ˚ κ i , A ˙ A κ i A ÿ A ¯ “ f i ˆ A ¯ E ˚ κ i , A ¯ ˙ A ź A ˛ “ A ¯ η int ˆ A ˛ E ˚ κ i , A ˛ ˙ (C12)with elasticities of the knock out nucleon given by κ PD “ κ PP “ . f PD “ ` τ PD { τ PP , and f PP “ ´ f PD “ ` τ PP { τ PD . (C13)The total escape rate of particles of mass A ˚ from injected nuclei of mass A is Q totesc p E ˚ , A ˚ , A q “ Q remesc p E ˚ , A ˚ , A q ` δ A ˚ “ Q PDesc p E ˚ , A ˚ , A q ` Q PPesc p E ˚ , A ˚ , A q ‰ , (C14)where δ A is the Kronecker delta.
2. Branching Ratios from Photo-disintegration
The propagation scheme described in the last section can be easily extended to take into account the emission ofseveral nucleons or light nuclei in photo-nuclear reactions. We use the total interaction time and branching ratios forphoto-dissociation from
Talys [43, 44] as available in
CRPropa and neglect multi-nucleon emission for photo-pionproduction since it can be safely neglected at the energies relevant here.Instead of the closed formulae derived above for the single-nucleon case, we now have the following recursive relationfor the rate of produced remnant nuclei of mass A ˚ . Q p E ˚ , A ˚ q “ A ´ A ˚ ÿ i “ b p E i , A ˚ , A i q η int p E i , A i q Q p E i , A i q ˇˇˇˇ dE i dE ˚ ˇˇˇˇ (C15)with A i “ A ˚ ` i , E i “ A i { A ˚ E ˚ and ˇˇ dE i dE ˚ ˇˇ “ A i { A ˚ . b p E i , A ˚ , A i q is the branching ratio that gives the probabilitythat a nucleus of mass A i with energy E i will have a remnant mass of A ˚ after the interaction. The n knocked outnucleons and nuclei are calculated the same way but replacing the branching fraction by n b p E i , n A ˚ , A i q , i.e. theprobability to produce n fragments of mass A ˚ , and summing over n . The rest of the calculation proceeds as in thecase of single-nucleon emission.
3. Proton Interactions
Once nucleons are generated in photo-disintegration they are assumed to either escape immediately in the case ofneutrons, or to interact further via photo-pion production. The average elasticity of this process is κ PP “ . E “ lg κ PP « .
1. Since we perform the calculation in logarithmic bins ofthis width, proton interactions can be treated similarly to photo-disintegration as a “trickle-down” of particles fluxessubsequently shifted by one energy bin. For the reaction p ` γ Ñ π ` ` n , the neutron escapes and the interactionchain is finished. In case of p ` γ Ñ π ` p , the secondary proton has a reduced energy; it may also interact again.The neutron, proton and positive pion fluxes in an energy bin k are calculated from the recursive relations Q p E ˚ k , n q “ Q p E ˚ k , n q ` p ´ b pp q η int p E ˚ k ` , p q Q p E ˚ k ` , p q { κ (C16) Q p E ˚ k , p q “ Q p E ˚ k , p q ` b pp η int p E ˚ k ` , p q Q p E ˚ k ` , p q { κ (C17)9and Q p E ˚ k , π ` q “ p ´ b pp q η int p E ˚ k ` , p q Q p E ˚ k ` , p q {p ´ κ q , (C18)where b pp « . p ` γ Ñ π ` p and the un-primed fluxes are the sum of theknocked-out nucleons from Eq. (C15) and primary protons. The offset of 7 in the equation for the pion flux is due tothe energy shift of the pions, ∆ lg E “ lg p ´ κ PP q « . Appendix D: Cosmic Ray Production and Propagation in an Expanding Universe
To compare the spectra obtained in the last section, the particles need to be propagated to Earth. The numberof cosmic rays per unit volume and energy in the present universe is equal to the number of particles accumulatedduring the entire history of the universe and is comprised of both primary particles emitted by the sources andsecondaries produced in the photo-disintegration process. Herein, the variable t characterizes a particular age of theuniverse and t H indicates its present age. We adopt the usual concordance cosmology of a flat universe dominatedby a cosmological constant, with Ω Λ « .
69 and a cold dark matter plus baryon component Ω m « .
31 [86]. TheHubble parameter as a function of redshift z is given by H p z q “ H r Ω m p ` z q ` Ω Λ s , normalized to its valuetoday, H “ h km s ´ Mpc ´ , with h » .
68 [86]. The dependence of the cosmological time with redshift can beexpressed via dz “ ´ dt p ` z q H p z q . The co-moving space density of cosmic rays n CR of mass A from a population ofuniformly distributed sources with (possibly age-dependent) emission rate per volume Q p E , A , t q is given by n CR p E, A, A q ” dN CR dE dV “ ż E ż t H d P AA p E , E, t q dE Q p E , A , t q ξ p t q dE dt , (D1)where d P AA { dE is the expectation value for the number of nuclei of mass A in the energy interval ( E, E ` dE q whichderive from a parent of mass A and energy E emitted at time t , and ξ p t q is the ratio of the product of co-movingsource density and Q p E , A , t q , relative to the value of that product today. Note that d P AA { dE includes propagationeffects both at the source environment and en route to Earth.We assume that the emission rate of cosmic rays is the same for all sources and the spectrum and composition isindependent of the age of the universe, so that evolution of the volumetric emission rate with cosmological time canbe described by an overall source evolution factor, ξ p t q discussed below. (It need not be specified whether this is dueto an evolution of the number of sources or their intrinsic power.) We further assume, as per usual practice, thatemission rate is fairly well described by a power-law spectrum. Under these general assumptions the source emissionrate per volume takes the form Q p E , A q “ Q ˆ E E ˙ γ exp ˆ ´ E Z E p max ˙ , (D2)where E p max is the maximal energy of emitted protons, i.e., maximum rigidity of the accelerator, Z is the nucleus’atomic number, E is some reference energy, and Q “ $&% . n dN A dE ˇˇˇ E “ E , for bursting sources n dN A dE dt ˇˇˇ E “ E , for steady sources . (D3)Here, . n is the number of bursts per unit volume per unit time and dN A { dE is the spectrum of particles producedby each burst, or for a steady source n is the number density of sources at z “
0, and dN A { dE dt is the UHECRproduction rate per unit energy per source. The cosmic ray power density above a certain energy E min is given by . (cid:15) E p A q “ ż E min E Q p E , A q dE “ Q ż E min E ˆ E E ˙ γ exp ˆ ´ E Z E p max ˙ dE “ Z E p max ˆ Z E p max E ˙ γ ` ż E min {p Z E p max q t γ ` e ´ t dt “ Q E ˆ Z E p max E ˙ γ ` Γ ˆ γ ` , E min Z E p max ˙ , (D4)0where Γ denotes the upper incomplete gamma function.The cosmological evolution of the source density per co-moving volume is parametrized as n s p z q “ n ξ p z q (D5)with ξ p z “ q “
1. We adopt for the fiducial model that the evolution of sources follows the star formation rate with ξ p z q “ p ` z q a ` rp ` z q{ b s c (D6)where a “ . ˘ . b “ . ˘ .
14 and c “ . ˘ .
19 [60]. Additionally we consider the family of evolutionmodels parameterized as ξ p z q “ p ` z q m .To propagate the particles escaping the source environment to Earth we use the CRPropa framework [83, 84].For this purpose, we generate a library of propagated nuclei with A ˚ “ . . . A ˚ max injected uniformly in light-traveldistance. The latter corresponds to a non-evolving source distribution in comoving distance after accounting for thecosmological time dilation. We simulated particles up to A ˚ max “
56. Given this library of simulated particles, we canconstruct the propagation matrix M ijµν for arbitrary source evolutions for each nuclear mass A ˚ µ escaping the sourceand secondary mass A ν at Earth. The elements of the propagation matrix give the expected number of secondaries inan energy interval r lg E j , lg E j ` ∆ s at Earth originating from nuclei at the source at an energy r lg E ˚ i , lg E ˚ i ` ∆ s fora given source evolution ξ p z p t qq and a uniform logarithmic spacing in energy with ∆ “ .
1. Numerically, the elementsare constructed via discretization of Eq. (D1) n CR p E j , A ν , A q “ A ÿ A ˚ µ “ A ν n i ÿ i “ j n a ÿ a “ ∆ P ijµνa ∆ E j Q totesc p E ˚ i , A ˚ µ , A q ξ p t a q ∆ t a ∆ E ˚ i , (D7)where ∆ t a “ t H { n a and∆ P ijµνa ∆ E j “ E j N Earth ijµνa p E ˚ i , E ˚ i ` ∆ E ˚ i ; E j , E j ` ∆ E j ; A ˚ µ ; A ν ; t a , t a ` ∆ t a q N gen iµa p E ˚ i , E ˚ i ` ∆ E ˚ i ; A ˚ µ ; t a , t a ` ∆ t a q (D8)For a non-evolving injection rate per unit volume, the number of generated events per bin is constant, N gen iµa “ K gen iµ .Then, for any source evolution ξ r z p t qqs , (D7) can be rewritten as n CR p E j , A ν , A q “ A ÿ A ˚ µ “ A ν n i ÿ i “ j Q totesc p E ˚ i , A ˚ µ , A q ∆ E ˚ i ∆ E j n a ÿ a “ N Earth ijµνa N gen iµa ξ r z p t a qs ∆ t a “ A ÿ A ˚ µ “ A ν n i ÿ i “ j Q totesc p E ˚ i , A ˚ µ , A ν q ∆ E ˚ i ∆ E j t H ř n a a “ N Earth ijµνa ξ r z p t a qs n a K gen iµ “ A ÿ A ˚ µ “ A ν n i ÿ i “ j Q totesc p E ˚ i , A ˚ µ , A q ∆ E ˚ i ∆ E j t H ř n a a “ N Earth ijµνa ξ r z p t a qs ř n a a “ N gen iµa “ A ÿ A ˚ µ “ A ν n i ÿ i “ j Q totesc p E ˚ i , A ˚ µ , A q ∆ E ˚ i ∆ E j t H ř N Earth ijµν p “ ξ r z p t p qs N gen iµ “ A ÿ A ˚ µ “ A ν n i ÿ i “ j ∆ E ˚ i ∆ E j t H M ijµν Q iµ , (D9)where ř N Earth ijµν p “ ξ r z p t p qs denotes the ξ -weighted sum over all events generated with p A ˚ µ , E ˚ i q arriving at Earth with p A ν , E j q and N gen iµ is the total number of generated events with p A ˚ µ , E i q . Note that if binned in ∆ z b “ z max { b , then N gen iµb | ∆ z b { ∆ t b | “ constant, and hence (D9) can be rewritten as n CR p E i , A ν , A q “ A ÿ A ˚ µ “ A ν n i ÿ i “ j Q totesc p E ˚ i , A ˚ µ , A q ∆ E ˚ i ∆ E j z max ř n a a “ N Earth ijµνa ξ p z b q ř n a a “ N gen iµa ˇˇˇ ∆ z a ∆ t a ˇˇˇ , (D10)1where | ∆ z b { ∆ t b | “ p ` z b q H p z b q and z max “ ∆ z b n b .For a given spectrum of injected nuclei of mass A we obtain the space density of cosmic rays at Earth with energy E and mass A , n CR p E, A, A q “ dN CR dE dV . (D11)For an isotropic arrival direction distribution (which is an excellent approximation based on current observations) therelation between the spectrum and the cosmic ray density is J p E, A, A q ” dN CR dE dA dt d Ω “ c π n CR p E, A, A q . (D12)The total flux at earth of particles of mass A ν is J p E, A ν q “ A max ÿ A µ “ A ν f p A µ q J p E, A ν , A µ q , (D13)where f p A µ q denotes the fraction of particles of mass A µ injected at the source. Appendix E: Neutrino and Photon Production
The results of the last two sections can be readily applied to obtain the flux of neutrinos at Earth from the decayof neutrons and charged pions. We approximate the emission rate of pions from photo-pion production by using κ π “ ´ κ PP in Eq. (C12). The energies of neutrinos escaping from the source are given by the kinematics of thetwo-body decay of pions and the subsequent muon decay which we treat approximately by assigning a third of themuon energy to each of the decay products. In this way we construct a propagation matrix for pions, M ijµν with ν “ p ν µ , ¯ ν µ , ν e p ¯ ν e qq and µ “ π ˘ . Similarly, a propagation matrix for neutrons is obtained with ν “ p ¯ ν e q and µ “ n .Neutrino oscillation over astronomical distances modifies the initial flavor distribution of fluxes, Φ e : Φ µ : Φ τ , incalculable ways. The relevant parameters for such a calculation are the three Euler rotations ( θ , θ , θ ) and the CP -violating Dirac phase δ . The current best-fit values as well as the allowed ranges of the mixing parameters atthe 1 σ level are: θ { ˝ “ . ` . ´ . , θ { ˝ “ . ` . ´ . ‘ . ` . ´ . , θ { ˝ “ . ` . ´ . , δ { ˝ “ ` ´ [87]. The mixingprobabilities are given by P ν µ Ñ ν µ “ c s ` p c c ` s s s ´ c c s s s c δ q ` p c s ` c s s ` c c s s s c δ q , (E1) P ν e Ø ν µ “ c (cid:32) c s c ` ` c ` s ˘ s s ` c s c s c δ p c ´ s q s ( , (E2) P ν e Ø ν τ “ P ν e Ø ν µ p θ Ñ θ ` π { q , P ν τ Ñ ν τ “ P ν µ Ñ ν µ p θ Ñ θ ` π { q , (E3)and the unitarity relations P ν e Ñ ν e “ ´ P ν e Ø ν µ ´ P ν e Ø ν τ P ν τ Ø ν µ “ ´ P ν e Ø ν µ ´ P ν µ Ñ ν µ P ν τ Ñ ν τ “ ´ P ν e Ø ν τ ´ P ν µ Ø ν τ , (E4)where c ij “ cos θ ij , s ij “ sin θ ij , and c δ “ cos δ [88]. The measurable neutrino flux at Earth is given by ¨˝ Φ e Φ µ Φ τ ˛‚ “ ¨˝ .
55 0 .
24 0 . .
24 0 .
37 0 . .
21 0 .
38 0 . ˛‚ ¨˝ Φ e Φ µ Φ τ ˛‚ . (E5)In addition to neutrinos, photons are produced from π production and decay [89], and by photo-disintegrationof high-energy nuclei followed by immediate photo-emission from the excited daughter nuclei [90]. The γ -rays, elec-trons, and positrons produced in the decay of π and π ˘ trigger an electromagnetic (EM) cascade on the cosmicmicrowave background, which develops via repeated e ` e ´ pair production and inverse Compton scattering. Othercontributions to the cascade are provided by Bethe-Heitler production of e ` e ´ pairs and γ -rays emitted during the2photo-disintegration process, after the photo-dissociated nuclear fragments de-excite. The net result is a pile up of γ -rays at GeV À E γ À TeV, just below the threshold for further pair production on the diffuse optical backgrounds.The EM energy then gets recycled into the so-called Fermi-LAT region, which is bounded by observation [91, 92]to not exceed ω cas „ . ˆ ´ eV { cm [93]. The latest Fermi-LAT limts [92] do not significantly influence thedetermination of the ω cas upper bound of [93], because that bound is not very sensitive to the high energy bins addedby [92].The photons coming from photo-pion production in the source environment were shown to be below the Fermi-LATbound in Section III A. To place a bound on the contribution of photons from nuclear de-excitation to the Fermi-LAT diffuse gammas, without performing an explicit calculation, we can turn to the estimate of [57] which found ω cas « . ˆ ´ eV { cm assuming the high-energy IceCube spectrum to be entirely due to neutron beta-decay. Sincein our model the neutrino flux from neutron decay is significantly below the IceCube spectrum, the correspondingde-excitation photon contribution to the Fermi-LAT data must be far below the limit. [1] T. Antoni et al. (KASCADE), Astropart. Phys. , 1(2005), arXiv:astro-ph/0505413 [astro-ph].[2] D. J. Bird et al. (HiRes), Phys. Rev. Lett. , 3401(1993).[3] R. U. Abbasi et al. (HiRes), Phys. Rev. Lett. , 101101(2008), arXiv:astro-ph/0703099 [astro-ph].[4] J. Abraham et al. (Pierre Auger), Phys. Lett. B685 , 239(2010), arXiv:1002.1975 [astro-ph.HE].[5] J. Abraham et al. (Pierre Auger), Phys. Rev. Lett. ,061101 (2008), arXiv:0806.4302 [astro-ph].[6] W. D. Apel et al. , Astropart. Phys. , 183 (2012).[7] M. G. Aartsen et al. (IceCube), Phys. Rev. D88 , 042004(2013), arXiv:1307.3795 [astro-ph.HE].[8] S. P. Knurenko, Z. E. Petrov, R. Sidorov, I. Ye. Sleptsov,S. K. Starostin, and G. G. Struchkov (Yakutsk), (2013),arXiv:1310.1978 [astro-ph.HE].[9] V. V. Prosin et al. (Tunka), Nucl. Instrum. Meth.
A756 ,94 (2014).[10] T. Abu-Zayyad et al. (HiRes-MIA), Astrophys. J. ,686 (2001), arXiv:astro-ph/0010652 [astro-ph].[11] D. R. Bergman and J. W. Belz, J. Phys.
G34 , R359(2007), arXiv:0704.3721 [astro-ph].[12] K.-H. Kampert and M. Unger, Astropart. Phys. , 660(2012), arXiv:1201.0018 [astro-ph.HE].[13] P. Abreu et al. (Pierre Auger), Astropart. Phys. , 627(2011), arXiv:1103.2721 [astro-ph.HE].[14] P. Abreu et al. (Pierre Auger), Astrophys. J. Suppl. ,34 (2012), arXiv:1210.3736 [astro-ph.HE].[15] A. Aab et al. (Pierre Auger), Astrophys. J. , 111(2015), arXiv:1411.6953 [astro-ph.HE].[16] A. Aab et al. (Pierre Auger), Phys. Rev. D90 , 122005(2014), arXiv:1409.4809 [astro-ph.HE].[17] A. Aab et al. (Pierre Auger), Phys. Rev.
D90 , 122006(2014), arXiv:1409.5083 [astro-ph.HE].[18] V. Berezinsky, A. Z. Gazizov, and S. I. Grigorieva, Phys.Rev.
D74 , 043005 (2006), arXiv:hep-ph/0204357 [hep-ph].[19] M. Ahlers, L. A. Anchordoqui, and A. M. Taylor,Phys. Rev.
D87 , 023004 (2013), arXiv:1209.5427 [astro-ph.HE].[20] M. Fukushima (Telescope Array),
Proceedings, 18th In-ternational Symposium on Very High Energy CosmicRay Interactions (ISVHECRI 2014) , EPJ Web Conf. ,04004 (2015), arXiv:1503.06961 [astro-ph.HE].[21] R. Aloisio, D. Boncioli, A. di Matteo, A. F. Grillo, S. Pe-trera, and F. Salamida, (2015), arXiv:1505.04020 [astro- ph.HE].[22] D. Allard, E. Parizot, and A. V. Olinto, Astropart. Phys. , 61 (2007), arXiv:astro-ph/0512345 [astro-ph].[23] D. Allard, A. V. Olinto, and E. Parizot, Astron. Astro-phys. , 59 (2007), arXiv:astro-ph/0703633 [ASTRO-PH].[24] D. Allard, N. G. Busca, G. Decerprit, A. V. Olinto,and E. Parizot, JCAP , 033 (2008), arXiv:0805.4779[astro-ph].[25] C. De Donato and G. A. Medina-Tanco, Astropart. Phys. , 253 (2009), arXiv:0807.4510 [astro-ph].[26] A. M. Taylor, Astropart. Phys. , 48 (2014),arXiv:1401.0199 [astro-ph.HE].[27] O. Deligny, Comptes Rendus Physique , 367 (2014),arXiv:1403.5569 [astro-ph.HE].[28] T. K. Gaisser, T. Stanev, and S. Tilav, Front.Phys.China , 748 (2013), arXiv:1303.3565 [astro-ph.HE].[29] R. Aloisio, V. Berezinsky, and P. Blasi, JCAP , 020(2014), arXiv:1312.7459 [astro-ph.HE].[30] G. Giacinti, M. Kachelriess, and D. Semikoz, Phys. Rev. D91 , 083009 (2015), arXiv:1502.01608 [astro-ph.HE].[31] D. Allard and R. J. Protheroe, Astron. Astrophys. ,803 (2009), arXiv:0902.4538 [astro-ph.HE].[32] K. Kotera, D. Allard, K. Murase, J. Aoi, Y. Dubois,T. Pierog, and S. Nagataki, Astrophys. J. , 370(2009), arXiv:0907.2433 [astro-ph.HE].[33] A. Pe’er, K. Murase, and P. Meszaros, Phys. Rev.
D80 ,123018 (2009), arXiv:0911.1776 [astro-ph.HE].[34] K. Fang, K. Kotera, and A. V. Olinto, Astrophys. J. , 118 (2012), arXiv:1201.5197 [astro-ph.HE].[35] K. Fang, K. Kotera, and A. V. Olinto, JCAP , 010(2013), arXiv:1302.4482 [astro-ph.HE].[36] N. Globus, D. Allard, R. Mochkovitch, and E. Pari-zot, Mon. Not. Roy. Astron. Soc. , 5270 (2015),arXiv:1409.1271 [astro-ph.HE].[37] E. Parizot,
Proceedings, Cosmic Ray Origin - beyond thestandard models , Nucl. Phys. Proc. Suppl. , 197(2014), arXiv:1410.2655 [astro-ph.HE].[38] K. Kotera, E. Amato, and P. Blasi, JCAP , 026(2015), arXiv:1503.07907 [astro-ph.HE].[39] A. P. Szabo and R. J. Protheroe, Astropart. Phys. , 375(1994), arXiv:astro-ph/9405020 [astro-ph].[40] R. J. Protheroe and T. Stanev, Astropart. Phys. , 185(1999), arXiv:astro-ph/9808129 [astro-ph].[41] L. O. Drury, P. Duffy, D. Eichler, and A. Mastichiadis, Proceedings, 26th International Cosmic Ray Conference, August 17-25, 1999, Salt Lake City , (1999), [Astron. As-trophys.347,370(1999)], arXiv:astro-ph/9905178 [astro-ph].[42] E. Kr¨ugel,
The Physics of Interstellar Dust (Institute ofPhysics Publishing, 2003).[43] .[44] We use the “TALYS-1.6 (restored)” cross sections as de-scribed in [45].[45] R. A. Batista, D. Boncioli, A. di Matteo, A. van Vliet,and D. Walz, (2015), arXiv:1508.01824 [astro-ph.HE].[46] A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe, andT. Stanev, Comput. Phys. Commun. , 290 (2000),arXiv:astro-ph/9903478 [astro-ph].[47] https://github.com/CRPropa/CRPropa3-data .[48] P. Blasi, R. I. Epstein, and A. V. Olinto, Astrophys. J. , L123 (2000), arXiv:astro-ph/9912240 [astro-ph].[49] W. Apel et al. , Phys.Rev.
D87 , 081101 (2013),arXiv:1304.7114 [astro-ph.HE].[50] A. Aab et al. (Pierre Auger), (2013), arXiv:1307.5059[astro-ph.HE].[51] B. R. Dawson et al. (Pierre Auger, Yakutsk, Tele-scope Array), EPJ Web Conf. , 01005 (2013),arXiv:1306.6138 [astro-ph.HE].[52] R. Abbasi et al. (Pierre Auger, Telescope Array), in (2015) arXiv:1503.07540 [astro-ph.HE].[53] P. Abreu et al. (Pierre Auger), JCAP , 026 (2013),arXiv:1301.6637 [astro-ph.HE].[54] E. J. Ahn for the Pierre Auger Collaboration, (2013),proc. 33rd ICRC (2013) arXiv:1307.5059 [astro-ph.HE].[55] M. G. Aartsen et al. (IceCube), Phys. Rev. D91 , 022001(2015), arXiv:1410.1749 [astro-ph.HE].[56] M. G. Aartsen et al. (IceCube), Phys. Rev.
D88 , 112008(2013), arXiv:1310.5477 [astro-ph.HE].[57] L. A. Anchordoqui, Phys. Rev.
D91 , 027301 (2015),arXiv:1411.6457 [astro-ph.HE].[58] M. G. Aartsen et al. (IceCube), (2014), arXiv:1412.5106[astro-ph.HE].[59] M. Ahlers, L. A. Anchordoqui, M. C. Gonzalez-Garcia,F. Halzen, and S. Sarkar, Astropart. Phys. , 106(2010), arXiv:1005.2620 [astro-ph.HE].[60] B. E. Robertson, R. S. Ellis, S. R. Furlanetto, and J. S.Dunlop, Astrophys. J. , L19 (2015), arXiv:1502.02024[astro-ph.CO].[61] R. C. Gilmore, R. S. Somerville, J. R. Primack, andA. Dominguez, Mon. Not. Roy. Astron. Soc. , 3189(2012), arXiv:1104.0671 [astro-ph.CO].[62] E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari,and T. Stanev, Phys. Rev. D80 , 094003 (2009),arXiv:0906.4113 [hep-ph].[63] S. Ostapchenko, Phys. Rev.
D83 , 014018 (2011),arXiv:1010.1869 [hep-ph].[64] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, andK. Werner, (2013), arXiv:1306.0121 [hep-ph].[65] T. Abu-Zayyad et al. (Telescope Array), Astrophys. J. , L1 (2013), arXiv:1205.5067 [astro-ph.HE].[66] K. A. Olive et al. (Particle Data Group), Chin. Phys.
C38 , 090001 (2014).[67] A. M. Taylor, M. Ahlers, and D. Hooper, (2015),arXiv:1505.06090 [astro-ph.HE].[68] A. H. Rosenfeld, Ann. Rev. Nucl. Part. Sci. , 555(1975). [69] K. Fang, K. Kotera, K. Murase, and A. V. Olinto,Phys. Rev. D90 , 103005 (2014), arXiv:1311.2044 [astro-ph.HE].[70] N. Globus, D. Allard, and E. Parizot, Phys. Rev.
D92 ,021302 (2015), arXiv:1505.01377 [astro-ph.HE].[71] D. Veberic, Comput. Phys. Commun. , 2622 (2012),arXiv:1209.0735 [cs.MS].[72] F. W. Stecker, Phys. Rev. , 1264 (1969).[73] J. L. Puget, F. W. Stecker, and J. H. Bredekamp, As-trophys. J. , 638 (1976).[74] S. Karakula and W. Tkaczyk, Astropart. Phys. , 229(1993).[75] L. A. Anchordoqui, M. T. Dova, L. N. Epele, andJ. D. Swain, Phys. Rev. D57 , 7103 (1998), arXiv:astro-ph/9708082 [astro-ph].[76] L. N. Epele and E. Roulet, JHEP , 009 (1998),arXiv:astro-ph/9808104 [astro-ph].[77] F. W. Stecker and M. H. Salamon, Astrophys. J. ,521 (1999), arXiv:astro-ph/9808110 [astro-ph].[78] E. Khan, S. Goriely, D. Allard, E. Parizot, T. Suomi-jarvi, A. J. Koning, S. Hilaire, and M. C. Duijvestijn,Astropart. Phys. , 191 (2005), arXiv:astro-ph/0412109[astro-ph].[79] D. Hooper, S. Sarkar, and A. M. Taylor, Phys. Rev. D77 , 103007 (2008), arXiv:0802.1538 [astro-ph].[80] R. Aloisio, V. Berezinsky, and S. Grigorieva, Astropart.Phys. , 73 (2013), arXiv:0802.4452 [astro-ph].[81] R. Aloisio, V. Berezinsky, and S. Grigorieva, Astropart.Phys. , 94 (2013), arXiv:1006.2484 [astro-ph.CO].[82] M. Ahlers and A. M. Taylor, Phys. Rev. D82 , 123005(2010), arXiv:1010.3019 [astro-ph.HE].[83] K.-H. Kampert, J. Kulbartz, L. Maccione, N. Niersten-hoefer, P. Schiffer, G. Sigl, and A. R. van Vliet, As-tropart. Phys. , 41 (2013), arXiv:1206.3132 [astro-ph.IM].[84] R. A. Batista, M. Erdmann, C. Evoli, K.-H. Kampert,D. Kuempel, et al. , (2013), arXiv:1307.2643 [astro-ph.IM].[85] L. A. Anchordoqui, J. F. Beacom, H. Goldberg,S. Palomares-Ruiz, and T. J. Weiler, Phys. Rev. D75 ,063001 (2007), arXiv:astro-ph/0611581 [astro-ph].[86] P. A. R. Ade et al. (Planck), (2015), arXiv:1502.01589[astro-ph.CO].[87] M. C. Gonzalez-Garcia,