A minimal model for extragalactic cosmic rays and neutrinos

TL;DR

This paper presents a unified minimal model explaining ultrahigh energy cosmic rays and neutrinos by combining hadronic and photo-hadronic interactions near sources, successfully matching multiple observational data sets.

Contribution

It introduces a minimal model incorporating both gas and photon backgrounds to explain cosmic ray and neutrino data simultaneously, consistent with gamma-ray background limits.

Findings

Good fit to cosmic ray flux, $X_{max}$, and $RMS(X_{max})$ data above 10^{17} eV

Predicted neutrino flux aligns with IceCube measurements

Contribution to diffuse gamma-ray background is about 30%

Abstract

We aim to explain in a unified way the experimental data on ultrahigh energy cosmic rays (UHECR) and neutrinos, using a single source class and obeying limits on the extragalactic diffuse gamma-ray background (EGRB). If UHECRs only interact hadronically with gas around their sources, the resulting diffuse CR flux can be matched well to the observed one, providing at the same time large neutrino fluxes. However, air showers in the Earth's atmosphere induced by UHECRs with energies $E>3\times 10^{18}$ eV would reach in such a case their maxima too high. Therefore additional photo-hadronic interactions of UEHCRs close to the accelerator have to be present, in order to modify the nuclear composition of CRs in a relatively narrow energy interval. We include thus both photon and gas backgrounds, and combine the resulting CR spectra with the high-energy part of the Galactic CR fluxes predicted by the escape model. As result, we find a good description of experimental data on the total CR flux, the mean shower depth $X_\max$ and its width $RMS(X_\max)$ in the whole energy range above $E\simeq 10^{17}$ eV. The predicted high-energy neutrino flux matches IceCube measurements, while the contribution to the EGRB is of order 30%.

Tables

Table 1. Table 1: Parameters used for the two cases shown in Fig. 1.

source type	only gas	photons+gas
power-law injection spectrum $\alpha$	$1.8$	$1.5$
maximal energy $E_{\max }/$eV	$3\times {10}^{18}$	$6\times {10}^{18}$
evolution	BL Lac	AGN
interaction depth ${\tau }_0^{pp}$	0.035	0.035
interaction depth ${\tau }_0^{p\gamma }$	0	0.29
photon temperature $T/$K	–	850
power-law index diffusion ${\delta }_{Ap}$	$0.5$	$0.5$
power-law index diffusion ${\delta }_{A\gamma }$	–	$0.77$

Equations

$$j_{\mathrm{inj}}(E)\propto E^{-\alpha }\exp [-E/(Z{E}_{\max })].$$

$$\rho (z)=\{\begin{array}{ll}(1+z{)}^m& \text{for }z<z_c,\\ (1+z_c{)}^m& \text{for }{z}_c<z<z_d,\\ (1+z_c{)}^{m}10^{k(z-z_d)}& \text{for }z>z_d,\end{array}$$

$$N_c(z)\propto {\int }_{L_{\gamma }^{min}}^{L_{\gamma }^{max}}\rho (z,L_{\gamma })L_{\gamma }d{L}_{\gamma },$$

$$\rho (z,L_{\gamma })=\rho (L_{\gamma })e(z,L_{\gamma }),$$

$$\rho (L_{\gamma })=\frac{A}{\log (10)L_{\gamma }}{\left[{\left(\frac{L_{\gamma }}{L_c}\right)}^{{\gamma }_1}+{\left(\frac{L_{\gamma }}{L_c}\right)}^{{\gamma }_2}\right]}^{-1},$$

$$e(z,L_{\gamma })={\left[{\left(\frac{1+z}{1+z_c(L_{\gamma })}\right)}^{p_1}+{\left(\frac{1+z}{1+z_c(L_{\gamma })}\right)}^{p_2}\right]}^{-1},$$

$$z_c(L_{\gamma })=z_c^{\star }{\left(\frac{L_{\gamma }}{10^{48}{\mathrm{erg}\mathrm{s}}^{-1}}\right)}^{\alpha }.$$

Full text

I Introduction

A major motivation for the construction of km3 neutrino telescopes has been the goal to identify the sources of ultrahigh energy cosmic rays (UHECR): While deflections of charged CRs in magnetic fields have prevented so far a successful correlation of their arrival directions with potential sources even at the highest energies, photons and neutrinos point back to their sources. These neutral secondaries are produced by UHECRs interacting with gas or photons in their sources, and with cosmic microwave and other background photons during propagation. Any process involving hadronization leads mainly to the production of pions, and isospin symmetry fixes then the ratio of charged to neutral pions produced. The production of neutrinos is thus intimately tied to the one of photons, and both depend in turn on the flux of primary cosmic rays. This is the basis of the multi-messenger approach to UHECR physics Gaisser:1994yf .

In this work, we will address the question if a single source class can explain $i)$ the extragalactic CR flux, $ii)$ its nuclear composition and $iii)$ the observed neutrino flux in IceCube. Moreover, we will require that $iv)$ the accompanying photon flux is only a subdominant contribution to the extragalactic gamma-ray background (EGRB) measured by Fermi-LAT EGRB . Finally, $v)$ the model should be consistent with an early galactic to extragalactic transition. No model has been developed yet which satisfies all five requirements. Since the neutrino flux measured by IceCube is high nuflux ; nuflux2 , i.e. close to the cascade limit casc , combining $iii)$ and $iv)$ is challenging for many source classes. Moreover, existing models aiming to reproduce the observed nuclear composition of UHECRs fail to produce sizable neutrino fluxes others1 ; Unger:2015laa . In contrast, models leading to large neutrino fluxes in the 0.1–1 PeV energy range use typically proton primaries with 10–100 PeV energies, without a direct connection to measurements of UHECR composition others2 ; more .

II Constraints

Let us now explain these conditions in more detail. Additionally to the all-particle CR spectrum, data on the primary composition have become available in the last years: The Auger collaboration derived the fraction of four different elemental groups above $6\times 10^{17}$ eV Aab:2014aea , while the KASCADE-Grande experiment measured the composition up to $2\times 10^{17}$ eV dataKG . These measurements can be summarized as follows: First, the proton fraction amounts to $\sim 40$–60% in the energy range between $7\times 10^{17}$ eV and $7\times 10^{18}$ eV and decreases afterwards, while the fraction of intermediate nuclei increases. Second, the iron fraction in the energy range between $7\times 10^{17}$ eV and $2\times 10^{19}$ eV is limited by \;<$$\sim\;15–20%. Despite both theoretical and experimental uncertainties, the following conclusions can be drawn: First, the Galactic contribution to the observed CR spectrum has to die out around $7\times 10^{17}$ eV. This inference is supported by limits on the CR dipole anisotropy which require a transition below $\sim 10^{18}$ eV in case of a light composition dipole . Consequently, we demand an early Galactic-extragalactic transition and have to explain therefore the ankle as a feature in the extragalactic CR spectrum. Second, the composition measurements are inconsistent with a strong dominance of either protons or iron nuclei. We assume therefore that a mixture of nuclei is injected in the source, with a rigidity dependent maximal energy.

The main part of the EGRB is attributed to unresolved blazars Neronov:2011kg ; dimauro ; Ajello:2015mfa ; Fermiun . This makes blazars and in particular BL Lacs attractive neutrino sources, since their contribution to the EGRB is much larger than the one from other sources. An attempt to connect the observed UHECR proton flux with neutrino and gamma-ray data was performed in Ref. escape3 . However, correlation studies of arrival directions of muon neutrinos with Fermi blazars showed that blazars cannot be the main source of IceCube neutrinos bllac_fermi . Thus leptonic models are favored to explain the main part of the photon flux from blazars. As a result, neutrino sources should give a subdominant contribution to the EGRB.

III Source model

To be specific, we will assume in the following that UHECRs are accelerated by (a subclass of) active galactic nuclei (AGN). Acceleration could proceed either via shock acceleration in accretion shocks KE86 or via acceleration in regular fields Bl76 ; Lo76 ; BZ77 ; MT82 ; Neronov:2007mh close to the supermassive black hole (SMBH). Alternatively, UHECRs could be accelerated in large-scale radio jets RB or via a two-step acceleration process in the jet and a radio-lobe lobes . We neglect the details of the acceleration process, and assume that the energy spectra of nuclei follow a power-law with a rigidity dependent cutoff

[TABLE]

Subsequently, the CR nuclei diffuse first through a zone dominated by photo-hadronic interactions, before they escape into a second zone dominated by hadronic interactions with gas. The propagation in both zones is modeled as a one-dimensional process, determined by the ratio $\tau$ of the interaction rate $R_{\rm int}$ to the escape rate $R_{\rm esc}$.

The spectra of AGN show a characteristic blue bump, produced by UV photons emitted from the accretion disk surrounding the SMBH bluebump . At lower energies, the spectra are dominated by thermal emission from dust surrounding the SMBH. We assume that these IR photons provide the main interaction target, and approximate their energy spectrum by a thermal distribution. In the next section, we will see that interactions on IR photons will result in a suppression of the heavier nuclei fluxes at the desired energy range. We parametrise the interaction depth $\tau$ as $\tau_{0}^{p\gamma}=R_{\rm int}/R_{\rm esc}$, using protons with energy $E_{0}=10$ EeV as reference. We assume that the accelerated CR nuclei escape from the region filled with thermal photons in a diffusive way, so that their escape rate $R_{\rm esc}$ is proportional to $(E/Z)^{\delta_{p\gamma}}$. For the numerical simulations, we use the open source code Kalashev:2014xna which is based on kinetic equations in one dimension. Since the propagation time is the same for all particles in the code, diffusion of charged particles is taken into account by multiplying the interaction rates by the escape times and setting the propagation time equal to one. Neutrons escape freely.

In the 2nd zone, we model both interactions with the gas and the escape of CRs as a Monte Carlo process, using QGSJET-II-04 qgsjetII to describe nucleus-proton collisions. We assume now that charged CRs diffuse in an extended halo and model escape and interactions in the leaky-box picture. The produced neutrons escape again freely. The source is then fully described specifying the interaction depth $\tau_{0}^{pp}=R_{\rm int}/R_{\rm esc}$ of protons at the reference energy $E_{0}=10$ EeV and the energy dependence of the escape rate, $R_{\rm esc}=R_{0}[E/(ZE_{0})]^{\delta_{pp}}$.

The spectrum of particles exiting the source is then used in the third step as an “effective injection spectrum”, from which we calculate using again the code Kalashev:2014xna the resulting diffuse flux, taking into account the distribution of sources as well as the interaction of protons, electrons and photons with the EBL and the CMB. We employ as our baseline EBL model the one of Ref. EBL , but have verified that other models, e.g. the more recent one Stecker:2016fsg , lead only to small variations in the CR flux and composition.

For the cosmological evolution of AGNs we use the parameterization obtained in Ref. Hasinger:2005sb ,

[TABLE]

choosing as parameters $m=3.4$, $z_{c}=1.2$, $z_{d}=1.2$, and $k=-0.32$. They correspond to the ones derived in Ref. Hasinger:2005sb for $\log(L_{X}/{\rm erg})=43.5$.

As second option, we use the parameterization for BL Lac/FR1 evolution presented in dimauro . We describe the cosmological evolution of BL Lac/FR I sources as

[TABLE]

where $\rho(z,L_{\gamma})$ is the $\gamma$-ray luminosity function, i.e. the number of sources per comoving volume and luminosity. For $\rho(z,L_{\gamma})$ we adopt the Luminosity-Dependent Density Evolution (LDDE) model of Ref. dimauro given by

[TABLE]

with

[TABLE]

and

[TABLE]

We use the numerical values given in Table 3 of Ref. dimauro for the free parameters of this model. The evolution of the effective source density as function of redshift is shown in Fig. 3 of Ref. unified . This figure illustrates that, in contrast to average AGNs, the number density of BL Lac and FR I galaxies peaks at low redshift, $z\lesssim 1$. Thus their evolution is similar to that of galaxy clusters. In fact, most of the FR I sources, which are the parent population of BL Lacs reside in the centres of the dominant central elliptical galaxies of galaxy clusters (cD galaxies). Thus these two models serve as templates for an evolution peaking early or late, respectively.

We perform steps 1–3 injecting pure p, He, N, Si or Fe primaries in the initial step and obtain then results for a mixed injection spectrum performing a linear combination of the final spectra. The model parameters used for the calculations are summarized in Table 1. The interaction depths $\tau_{0}^{pp}$ and $\tau_{0}^{p\gamma}$ given there are those for proton primaries at the reference energy $E_{0}=10^{19}$ eV.

IV Results

Having calculated for a specific set of source parameters the diffuse flux of UHECR nuclei and their secondaries, we combine them with the high-energy part of the Galactic CR fluxes predicted in the escape model escape2 . We aim to reproduce the total CR flux Aab:2014aea ; dataKG ; Aab:2013ika , the average maximum depth $X_{\max}$ of observed CR-induced air showers and the width of the $X_{\max}$ distributions XRMS . Since the experimental and theoretical systematic uncertainties in both the $X_{\max}$ and the RMS($X_{\max})$ measurements are difficult to quantify, we do not perform a standard $\chi^{2}$ analysis. Instead, we determine first the nuclear composition of extragalactic CRs which fits best the observed total CR flux, since this is the most reliable quantity. After that we check if the composition is allowed by the EGRB constraint and results in a sufficiently large neutrino flux.

We consider first the case that photo-hadronic interactions are negligible. Requiring both a small contribution to the EGRB and a large contribution to the neutrino signal observed by IceCube restricts the allowed range of slopes $\alpha$ strongly, \alpha\raise 1.29167pt\hbox{;<\kern-7.5pt\raise-4.73611pt\hbox{\sim;}}2.1. We choose $\delta_{Ap}=0.5$ as the energy-dependence of the escape rate $R_{\rm esc}=R_{0}[E/(ZE_{0})]^{1/2}$, corresponding to Kraichnan turbulence. The interaction and escape rates have been normalized such that $\tau^{pp}_{0}=0.035$ for proton primaries at $E_{0}=10$ EeV. Since the interaction depth decreases with increasing energies due to the faster CR escape, the UHECR flux is dominated by primary nuclei. In order to reproduce the light component in the KASKADE-Grande data, the extragalactic CR flux has to remain light up to $10^{18}$ eV. Insisting to reproduce the ankle requires a relatively low cutoff energy, $E_{\max}=3\times 10^{18}$ eV, and a small contribution of intermediate nuclei. Thus the wish to describe the ankle including the sub-ankle region by extragalactic CRs drives the composition towards a two-component model, consisting mainly of protons and iron. Note that in this scenario the spectra of intermediate CNO nuclei are cut off around the ankle energy, and hence their contribution is insignificant, unless the proton flux is strongly reduced. A reduction of the proton flux would in turn reduce the neutrino flux and the model will fall short of explaining the IceCube data.

The upper left panel of Fig. 1 shows the resulting diffuse fluxes of CR nuclei and their secondaries. Both the ankle, which corresponds in this scenario to the transition between the proton and iron dominated components in the extragalactic CR flux, and the proton (light) component observed by KASCADE-Grande are reproduced well. The contribution to the diffuse EGRB is low except in the TeV range, while the neutrino flux is somewhat below the level indicated by IceCube observations. The two lower left panels of Fig. 1 show the predicted shower maximum $X_{\max}$ and the corresponding distribution width RMS($X_{\max})$, respectively. Since the composition above the ankle is heavy, the predicted $X_{\max}$ coincides with the one of Fe, in contradiction to observations. Also the predicted width RMS($X_{\max})$ is smaller than the observed one.

In our second scenario, we add photo-hadronic interactions close to the source. In contrast to hadronic interactions, nuclear scattering on IR photons can give rise to relatively large interaction depths, leading to the dominance of secondary nuclei in the ankle region. The parameter space of this case, however without including hadronic interactions on gas, was extensively studied in Ref. Unger:2015laa . We employ therefore simply their base case, using $T=850$ K, $\tau^{p\gamma}=0.29$ and $\delta_{p\gamma}=0.77$ for diffusion close to the source. In contrast to Ref. Unger:2015laa , we assume however a compact acceleration region and use therefore an exponential attenuation. Moreover, we choose as injection slope $\alpha=1.5$ and $E_{\max}=6\times 10^{18}$ eV. The interaction depth for hadronic interactions is kept as before, i.e. normalised as $\tau_{0}^{pp}=0.035$. The resulting diffuse fluxes of CR nuclei and their secondaries are shown in the upper right panel of Fig. 1. The spectra of intermediate nuclei show a narrow dip: The low-energy end of this dip is determined by the threshold energy for photo-dissociation, which is at higher energies partly filled by secondary nuclei generated by heavier primary nuclei. The resulting neutrino flux matches now IceCube data, while the contribution to the diffuse EGRB is of order 30%. While neutrinos from interactions on gas dominate at \raise 1.29167pt\hbox{;<\kern-7.5pt\raise-4.73611pt\hbox{\sim;}}10^{16} eV, $A\gamma$ interactions become important at higher energies. The two lower right panels of Fig. 1 show the predicted shower maximum $X_{\max}$ and the corresponding distribution width RMS($X_{\max})$, respectively. Accounting for the systematic uncertainties of the experimental data and of the hadronic interaction models, the two distributions are well reproduced. Note that the peak in RMS($X_{\max})$ could be shifted to lower energies, reducing the extragalactic proton flux somewhat. This may indicate that not all neutrons escape freely from their sources.

V Discussion

Let us now comment in which astrophysical environments the considered UHECR interactions could be realized. Interaction depths of order one for proton-gamma interaction arise naturally from scattering on IR photons PASA . These photons may be either emitted by the dust torus of few parsecs extension or, as considered here, from a more compact source region M . Similarly, the dust and gas in the accretion disc surrounding the SMBH provides a target for hadronic interactions of UHECRs. In both cases, UHECRs have to be accelerated close to the SMBH, excluding e.g. acceleration sites as radio lobes of AGN jets. The composition of injected CRs has to be strongly enhanced towards intermediate nuclei compared to the solar composition. This may indicate a connection to the tidal ignition of white dwarf stars close the the SMBH, as suggested recently in Ref. TDE . Next we comment on the implication of our results for the EGRB. Since the AGN evolution peaks early, the spectral shape of the photon flux is close to the “universal shape” obtained after many cascade generations Berezinsky:2016feh , which in turn reproduces well the observed shape of the EGRB. Moreover, reproducing the large neutrino flux observed by IceCube requires that in our model unresolved AGN contribute around 30% to the EGRB.

Finally, we comment on the dependence of our results on the used experimental data and the hadronic interaction model. Using the data of the Telescope Array (TA) experiment Ai is complicated, since the effects of the experimental cuts and of the reconstruction bias have to be included. A preliminary analysis, using the method of Ref. Aii to mimic these effects, showed that the TA data favor a lighter composition than obtained here, represented by a mixture of protons and helium nuclei. On the other hand, such a mixture is disfavored by the PAO analysis using the correlations between the maximum depth $X_{\max}$ and the ground detector signal Aiii . Employing in the analysis the EPOS-LHC model would result in a somewhat heavier UHECR composition: because of the deeper $X_{\max}$ and smaller RMS($X_{\max})$ predicted by that model. However, as discussed in Refs. Av ; Avi , the deeper $X_{\max}$ values of EPOS are disfavored by the PAO analysis of the maximal muon production depth $X^{\mu}_{\max}$ Aiv . In addition, Ref. Avi suggested that the smaller RMS($X_{\max})$ for primary iron in EPOS results from deficiencies in its treatment of the fragmentation of nuclear spectators.

VI Conclusions

We have studied, if a single source class can explain both the flux and the composition of extragalactic CRs including the sub-ankle region and the high-energy neutrinos observed by IceCube. Using only hadronic interactions of UHECRs with gas around their sources, we have obtained a good fit to the CR energy spectrum, however, only if intermediate nuclei are sub-dominant. Therefore the predicted maximum $X_{\max}$ of CR-induced air showers lies higher in the atmosphere than observed. Adding photo-nuclear interactions, with a relatively large interaction depth, we were able to reduce significantly the fraction of heavy nuclei in the primary fluxes and, consequently, to fit satisfactorily both the spectrum and the composition data on UHECRs. At the same time, the high-energy neutrino flux obtained matches IceCube measurements, while the contribution of unresolved UHECR sources to the EGRB is of order 30%. The large interaction depth for photo-nuclear interactions suggests that UHECRs are accelerated close to SMBHs.

Acknowledgements.

S.O. acknowledges support from project OS 481/1 of the Deutsche Forschungsgemeinschaft. This research was supported in part with computational resources at NTNU provided by NOTUR, http://www.sigma2.no, and at the Theoretical Physics Division of the INR RAS with support by the RFBR, grant 16-29-13065 ofi-m. Note added: While we were preparing this work for submission, the preprint arXiv:1704.00015 appeared. The authors explain the neutrino flux measured by IceCube by CR interactions in the galaxy clusters surrounding UHECR sources. While the basic scheme may be considered as a specific realisation of our first secenario, our attitude differs. Since we require that the neutrino sources explain the observed UHECR flux, we discard the case $\alpha\simeq 1.9-2.1$ and AGN evolution considered in arXiv:1704.00015, because it does not fit the UHECR flux. In contrast to our assumptions, the scenario of arXiv:1704.00015 implies that neutrino sources are only a subdominant class of UHECR sources.