Linking High-Energy Cosmic Particles by Black-Hole Jets Embedded in Large-Scale Structures

TL;DR

This paper proposes that black hole jets in galaxy clusters can simultaneously explain the origins of ultrahigh-energy cosmic rays, PeV neutrinos, and gamma rays, unifying these phenomena within a single astrophysical framework.

Contribution

It introduces a model where black hole jets in galaxy clusters produce and accelerate cosmic rays, neutrinos, and gamma rays, linking these high-energy phenomena.

Findings

Model reproduces the observed UHECR spectrum and composition.

Explains IceCube neutrino flux and Fermi gamma-ray background.

Consistent with energy output estimates from black hole jets.

Abstract

The origin of ultrahigh-energy cosmic rays (UHECRs) is a half-century old enigma (Linsley 1963). The mystery has been deepened by an intriguing coincidence: over ten orders of magnitude in energy, the energy generation rates of UHECRs, PeV neutrinos, and isotropic sub-TeV gamma rays are comparable, which hints at a grand-unified picture (Murase and Waxman 2016). Here we report that powerful black hole jets in aggregates of galaxies can supply the common origin of all of these phenomena. Once accelerated by a jet, low-energy cosmic rays confined in the radio lobe are adiabatically cooled; higher-energy cosmic rays leaving the source interact with the magnetized cluster environment and produce neutrinos and gamma rays; the highest-energy particles escape from the host cluster and contribute to the observed cosmic rays above 100 PeV. The model is consistent with the spectrum, composition, and isotropy of the observed UHECRs, and also explains the IceCube neutrinos and the non-blazar component of the Fermi gamma-ray background, assuming a reasonable energy output from black hole jets in clusters.

Equations

$$\Phi (E)=\frac{1}{4\pi }\int \frac{cdz}{H(z)}F(z){\int }_{M_{\min }}^{\infty }dM\frac{dn}{dM}\frac{d\dot{N}}{d{E}^{\prime }}(M,z),$$

$$E_{\gamma}Q_{E_{\gamma}}\approx\frac{2}{3}(E_{\nu}Q_{E_{\nu}})\big{|}_{E_{\nu}=E_{\gamma}/2}\,,$$

$$\frac{dn}{dM}(M,z)=f(\sigma )\frac{{\rho }_m}{M}\frac{d\ln {\sigma }^{-1}}{dM}$$

$${\chi }^2(C_{\mathrm{norm}},{\delta }_E^{\mathrm{CR}},{\delta }_E^{\nu })={\chi }_{\mathrm{CR},\mathrm{spec}}^2(C_{\mathrm{norm}},{\delta }_E^{\mathrm{CR}})+{\chi }_{\mathrm{CR},X_{\max }}^2+{\chi }_{\nu ,\mathrm{spec}}^2(C_{\mathrm{norm}},{\delta }_E^{\nu }).$$

$${\chi }_{CR/\nu ,\mathrm{spec}}^2$$

$${\chi }_{\mathrm{CR},X_{\max }}^2$$

Full text

The origin of ultrahigh-energy cosmic rays (UHECRs) is still unknown Hillas (1984). Measurements by the Pierre Auger Observatory (Auger) Aab et al. (2015) and the Telescope Array (TA) Charles et al. (2015) find a power-law spectrum, $\Phi\propto E^{-2.6}-E^{-2.7}$ (where $E$ is particle energy and $\Phi$ is the diffuse intensity in units of particles per energy, area, time, and solid angle), with a decline above $6\times{10}^{19}\,\rm eV$, probably due to the interaction of UHECRs with cosmic radiation backgrounds such as the cosmic microwave background (CMB) or an upper limit of the particle energy reachable by the accelerator. Small-scale anisotropy in their arrival directions has not been established Aab et al. (2015); Charles et al. (2015).

The IceCube Observatory recently discovered high-energy cosmic neutrinos Halzen (2016); Aartsen et al. (2013), which have been anticipated to provide crucial clues to this age-old mystery. An astrophysical flux in the $0.1-1$ PeV range is found at the level of $10^{-8}\,\rm GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}$ per flavor Halzen (2016); Aartsen et al. (2013, 2016a, 2017a), which is consistent with expectations of cosmic-ray “reservoir” models Murase et al. (2008); Kotera et al. (2009); Loeb and Waxman (2006); Murase et al. (2013). The arrival directions of the observed events present no significant clustering, and indicate that the sources are extragalactic Halzen (2016); Aartsen et al. (2016a, 2017a).

A $\gamma$-ray counterpart is expected from the hadronic processes responsible for neutrino production. If the source environment is transparent to $\gamma$ rays, these side products should show up at $1-100$ GeV energies after cascading in the extragalactic background light (EBL). A significant fraction of the extragalactic $\gamma$-ray background (EGB) Ackermann et al. (2015, 2016) measured by the Fermi Gamma-Ray Space Telescope may be explained by neutrino sources Murase et al. (2013). Despite the unknown origins of these multi-messenger emissions and fine structures in their data (such as a possible excess in the $10-100$ TeV neutrino spectrum Halzen (2016)), it is remarkable that over ten orders of magnitude in energy, the energy generation rates of UHECRs, IceCube neutrinos, and Fermi EGB are all comparable Murase et al. (2013); Murase and Waxman (2016).

Recent UHECR observations have revealed additional characteristic features of extragalactic cosmic rays. First, a hardening in the spectrum of light particles is seen around 100 PeV Apel et al. (2013); Buitink et al. (2016), right in the energy range where a steepening in the spectrum of heavy primary particles is observed (which is often called the “second knee”). Second, a transition from light elements to medium-to-heavy elements around ${10}^{19}$ eV is suggested by the Auger data Aab et al. (2015), and a heavy UHECR composition is also supported by the non-detection of cosmogenic neutrinos Aartsen et al. (2016b). Although the interpretation of the UHECR composition is still debated, direct measurements of indicators of the particle mass seem to be consistent between different experiments. These features were not considered in the simplest convergence theory Murase and Waxman (2016); here we provide a concrete astrophysical model in which black hole jets embedded in large-scale structures reconcile these observations.

Relativistic jets from the accretion onto supermassive black holes provide promising sites for UHECR acceleration. The Hillas condition suggests that an active galactic nucleus (AGN) can accelerate a particle with charge $Z$ to a maximum energy, $E_{\rm max}\sim Z\,10^{19}\,\rm eV$, in jets or at external shocks that are known to be sites of leptonic emissions Hillas (1984); Murase et al. (2012). The energy spectrum of particles accelerated by the Fermi mechanism can be described by a power law, $dN_{\rm acc}/dE\propto E^{-s_{\rm acc}}$, with an index $s_{\rm acc}\sim 2-2.5$. Radio observations often find extended lobes (sometimes referred as bubbles or cocoons, which are plasma cavities inflated by the jet), with $10-100$ kiloparsec scales Kaiser and Best (2007) and $0.1-10$ microgauss-level magnetic fields Kataoka and Stawarz (2005). Particles with energies below $E_{{\rm lobe},c}=Z\,e\,B_{\rm lobe}\,l_{{\rm lobe},c}\sim 1.4\times{10}^{18}\,Z\,\left(B_{\rm lobe}/5\,\mu\rm G\right)\,\left(l_{{\rm lobe},c}/0.3\,\rm kpc\right)\,\rm eV$ have a Larmor radius $r_{L}=E/(Ze\,B)$ that is much smaller than the coherence length of their magnetic structure, where $B_{\rm lobe}$ and $l_{{\rm lobe},c}\sim(0.01-0.1)\,r_{\rm lobe}$ are the magnetic field strength and the coherence length, and $r_{\rm lobe}$ is the lobe size and $c$ is the speed of light. Unlike relativistic electrons that cool inside jets or lobes, high-energy ions diffuse for $t_{\rm diff}^{\rm lobe}\sim 6.1\,\left({r_{\rm lobe}}/{10\,\rm kpc}\right)^{2}\,\left({E}/{Z\,1\,\rm PeV}\right)^{-1/3}\,\left({l_{{\rm lobe},c}}/{0.3\,\rm kpc}\right)^{-2/3}$ $\left({B_{\rm lobe}}/{5\,\mu\rm G}\right)^{1/3}\,{\rm Myr}$ and can enter the intracluster medium (ICM). Meanwhile, these high-energy ions suffer from adiabatic losses due to the expansion of the cocoon. The characteristic cooling time is $t_{\rm ad}\sim 4.9\,\left(r_{\rm lobe}/10\,\rm kpc\right)$ $\left(v_{\rm lobe}/2000\,\rm km\,s^{-1}\right)^{-1}\,\rm Myr$, where $v_{\rm lobe}$ is the typical expansion velocity at source ages of $\sim 0.1-10$ Myr Bordas et al. (2011). Particles with energies above $E_{{\rm lobe},c}$ are less impacted, escaping semi-diffusively with $t_{\rm diff}^{\rm lobe}\propto E^{-2}$. Considering the competition between diffusion and cooling, we approximate the spectrum of cosmic rays leaking into the cluster to be $dN_{\rm inj}/dE\propto E^{-s_{\rm acc}}\,\exp(-t_{\rm diff}^{\rm lobe}/t_{\rm ad})$.

Radio-loud AGN activity that leads to cosmic-ray injections would preferentially reside in the centers of rich clusters Best et al. (2007). A cluster with a halo mass of $M=10^{14}\,M_{14}\,M_{\odot}$ has a virial radius $r_{\rm vir}\sim 1.2\,M_{14}^{1/3}\,\rm Mpc$. The distribution of thermal gas is often described using the $\beta$ model as $n_{\rm ICM}(r)\propto\left[1+\left({r}/{r_{c}}\right)^{2}\right]^{-3\beta/2}$, where $\beta\approx 0.8$ and $r_{c}\sim 0.1\,r_{\rm vir}$ is the core radius Brunetti and Jones (2014). Turbulent magnetic fields in the ICM, which are probably induced by accretion shocks and other cluster dynamics, typically have a strength of a few $\mu$G in the cluster center Brunetti and Jones (2014). Assuming flux conservation and that the field traces the baryon distribution, we adopt a magnetic field profile $B(r)=B_{0}\left[1+\left({r}/{r_{\mathrm{c}}}\right)^{2}\right]^{-\beta}$ with $B_{0}\sim 5\mu G$.

Cosmic rays leaving the acceleration site and lobe enter the ICM of the host cluster (which functions as a cosmic-ray reservoir Murase et al. (2008); Kotera et al. (2009)). The highest-energy ions travel in a straight line through the ICM. Particles reaching an energy $E_{c}\sim 2\times 10^{19}\,Z\,B_{-6}\,(l_{c}/{20\,\rm kpc})\,\rm eV$ have a gyro-radius comparable to the typical scales of magnetic field fluctuations in massive clusters, with $l_{c}$ about 1-10% of the virial radius Brunetti and Jones (2014). Ions with energies well below $E_{c}$ propagate diffusively in the turbulent magnetic field of the cluster. The confinement, which could last for $\sim 1-10$ Gyr depending on the particle energy, leads to efficient interactions of cosmic-ray nuclei with baryons and infrared background photons in the cluster, producing pions that decay into neutrinos and $\gamma$ rays via $\pi^{\pm}\rightarrow\nu_{e}(\bar{\nu}_{e})+e^{\pm}+\nu_{\mu}+\bar{\nu}_{\mu}$ and $\pi^{0}\rightarrow 2\gamma$, respectively. Finally, particles that leave the cluster propagate to the Earth through the intergalactic medium and extragalactic magnetic fields. UHECRs from sources beyond the energy-loss horizon are depleted via photodisintegration, photomeson production and Bethe-Heitler pair production processes with the CMB and the EBL, producing cosmogenic neutrinos peaked around EeV and $\gamma$ rays that cascade down to GeV-TeV energies.

We numerically simulate the propagation of cosmic rays in the magnetized ICM and from the source to the observer. We assume that a jetted source as a cosmic-ray accelerator can be anywhere in the core of a cluster with equal probability. We inject five representative groups of elements: hydrogen (1H), helium (4He), nitrogen (14N), silicon (28Si) and iron (56Fe) according to the abundances of elements in Galactic cosmic rays (see Supplementary Information for details), and let each group follow the same power-law spectrum with a cutoff above the maximum rigidity, $dN_{\rm inj}/dR\propto R^{-s_{\rm acc}}\,\exp(-R/R_{\rm max})$, where $R=E/Ze$ is the rigidity, $s_{\rm acc}=2.3$, and $R_{\rm max}=2\times 10^{21}/26\,{\rm V}$. We assume that ions are confined up to $t_{\rm inj}=2\,\rm Gyr$, given that the peak period of AGN activity effectively lasts for $\sim 2-3$ Gyr (see Supplementary Information for discussions on model uncertainties and details). The redshift evolution of the source density is taken to be $F(z)=(1+z)^{3}$ up to $z_{c}=1.5$, but its moderate variations barely impact our results. The cumulative flux is obtained by Murase et al. (2008):

[TABLE]

where $n$ is the halo number density, $dn/dM$ is the halo mass function, $H(z)$ is the Hubble parameter at redshift $z$, $d\dot{N}/dE^{\prime}$ is the production rate of neutrinos (or propagated cosmic rays) from a given cluster with a redshifted energy $E^{\prime}=(1+z)\,E$. We consider clusters with a halo mass above $M_{\rm min}=5\times 10^{13}\,M_{\odot}$ (corresponding to $\sim 10^{11}\,M_{\odot}$ for the stellar mass of the main halo), which present higher radio-loud AGN fractions Best et al. (2007). For the intergalactic propagation, we assume that cosmic rays from a galaxy cluster have 50% chance of encountering magnetic structures with an average strength of 2 nG and a coherence length of 1 Mpc.

Figure 1 shows the integrated spectra of UHECRs and neutrinos from overdense regions with black hole jets. The normalization of the spectra is determined by a combined fit to the Auger spectral and $\langle X_{\rm max}\rangle$ data above $10^{18.45}$ eV, and the IceCube data above $2\times{10}^{14}$ eV. The goodness fit results in a $\chi^{2}=44.5$ for 30 degrees of freedom, corresponding to a p-value of 0.043 for this fiducial case. The cosmic-ray confinement in the lobe and the host cluster makes the injection spectrum harder below the second knee Murase et al. (2008, 2013). The spectral shape is agreement with measurements by both Auger and TA above $10^{18}\,\rm eV$. Primary and secondary cosmic-ray particles received by the observer are divided into two composition groups: light (including H and He) and intermediate/heavy (including CNO, Si, Mg, Fe), with the two crossing around $10^{19.5}\,\rm eV$. The mean of the maximum depth of an air shower, $\langle X_{\rm max}\rangle$, which depends on the mass of the UHE nucleon or nucleus, is shown in Figure 2. The trend follows the $\langle X_{\rm max}\rangle$ data measured by Auger. Below $10^{18}\,\rm eV$, accounting for a Galactic contribution with $\Phi\propto E^{-3.4}$, the predicted cosmic-ray spectrum matches the light component of the KASCADE-Grande data Apel et al. (2013).

The neutrino spectrum is composed of two parts. Between $10^{14}\,\rm eV$ and $10^{17}\,\rm eV$, it is mostly contributed by particle interactions in the ICM. It agrees with the IceCube measurements above $10^{14}\,\rm eV$. The low-energy neutrino spectrum is harder than that of accelerated cosmic rays, and the spectral steepening above $10^{15}\,\rm eV$ results from the faster escape of higher-energy cosmic rays. Above $10^{18}\,\rm eV$, the neutrino flux is dominated by the cosmogenic neutrinos produced when UHECRs interact with the CMB and the EBL, and is consistent with the IceCube constraints at extremely high energies Aartsen et al. (2016b). Likewise, the observed sub-TeV $\gamma$ rays are produced both in the ICM and during intergalactic propagation Murase and Waxman (2016). Thanks to the hard injection spectrum, the total $\gamma$-ray flux largely originates from electromagnetic cascades, and is consistent with the non-blazar component of the EGB Ackermann et al. (2016). In addition to the hard $\gamma$-ray spectrum, our model also predicts a dominance of low-mass clusters, and the $\gamma$-ray and radio limits from individual clusters Zandanel et al. (2015) can be satisfied.

The chance of previously or currently having active jets in a cluster, $f_{\rm jet}$, and the average cosmic-ray luminosity of contained active galaxies per cluster, $L_{\rm CR}$, are left as free parameters. Assuming $L_{\rm CR}\sim 10^{44}-10^{45}\,\rm erg\,s^{-1}$, we obtain $f_{\rm jet}\sim 10-100\%$. This is consistent with duty cycles of the accretion-driven evolution of black holes Ma et al. (2011). The number density of clusters and groups with a mass above $5\times 10^{13}\,M_{\odot}$ is $\sim{\rm a~{}few}\times 10^{-5}\,\rm Mpc^{-3}$. This satisfies the muon neutrino limits on the neutrino source density derived from the absence of multiplets Murase and Waxman (2016), as well as the lower bounds on the UHECR source density derived from the lack of strong anisotropy in the UHECR data Abreu et al. (2013). Our model predicts an association between directions of neutrino events and low-mass clusters that have past and ongoing jet activities of supermassive black holes. Alongside multiplet signals from nearby candidate sources, this correlation can be tested Murase and Waxman (2016) with future neutrino observations by experiments such as IceCube-Gen2.

**Data availability

**The authors declare that the data supporting the plots within this paper and other findings of this study are available from the authors upon reasonable request. The data of Figures 1 and 2 of the main text can be found at https://figshare.com/s/8216c3831633e29dace3.

**Acknowledgements

**We thank Rafael Alves Batista, Mauricio Bustamante, M. Coleman Miller, Chris Reynolds and Michael Unger for helpful comments. This work made use of supercomputing resources at the University of Maryland. We gratefully acknowledge support from the Eberly College of Science of Penn State University and the Institute for Gravitation and the Cosmos. The work of K.M. is supported by Alfred P. Sloan Foundation and NSF grant No. PHY-1620777.

**Author contributions

**K.F. performed simulations and produced the figures. K.M. designed the research and contributed to the calculations. Both authors equally edited the manuscript.

**Competing financial interests

**The authors declare no competing financial interests.

**Additional information

**Reprints and permissions information is available at www.nature.com/reprints. Correspondence and requests for materials should be addressed to K.M.

**Note added

**After this Letter was submitted to Nature Physics, the preprint arXiv:1704.06893 appeared with an incorrect quote about our work. Our reservoir model indeed explains diffuse spectra of neutrinos, $\gamma$ rays and UHECRs “dominantly”, as well as the UHECR composition and sub-ankle cosmic rays.

Supplementary Information

Methods

Cosmic rays can be accelerated by black hole jets by mechanisms such as shock acceleration, shear acceleration, and magnetic reconnection. After escaping from the jet and lobe/bubble, cosmic rays rectilinearly or diffusively propagate in the turbulent magnetic field of the host cluster (or group). In our grand-unified model, while we consider AGN jets as the principal cosmic-ray accelerators, we treat galaxy clusters and groups as cosmic-ray reservoirs, in which PeV neutrinos and a fraction of sub-TeV $\gamma$ rays are produced Murase et al. (2008, 2009); Kotera et al. (2009). The semi-diffusive propagation is implemented based on the public UHECR propagation code CRPropa 3 Alves Batista et al. (2016). The diffusive propagation is computed semi-analytically by letting particles random walk with a step size equal to the coherence length of the magnetic field. The numerical step size is then converted to an actual trajectory length through the on-site diffusion coefficient Kotera and Lemoine (2008); Fang and Olinto (2016). Interactions between nuclei and target nucleons of the ICM gas are computed using pre-tabulated cross sections and products that were calculated with the hadronic interaction model EPOS Werner et al. (2006). Photomeson productions (photodisintegration) due to interactions between cosmic-ray protons (nuclei) and the infrared background of clusters are computed based on SOPHIA Mücke et al. (2000) (TALYS Koning and Duijvestijn (2013)) through CRPropa 3.

Galaxy clusters and groups, in which AGN and other cosmic-ray accelerators are embedded, provide promising sites for neutrino and gamma-ray production Berezinsky et al. (1997); Murase and Beacom (2013). The density profile of the intracluster medium gas in a cluster is normalized by $f_{\rm b}M/(\mu\,m_{p})=\int n_{\rm ICM}(r)dV$, where $f_{b}$ is the gas fraction in galaxy clusters, which is compatible with the cosmic mean baryon fraction $\Omega_{b}/\Omega_{m}=0.167$ Komatsu et al. (2011); Planck Collaboration et al. (2013), $\mu\approx 0.61$ is the mean molecular weight, and $m_{p}$ is the mass of a proton. The infrared background photons in the cluster is modeled following reference Takami and Murase (2012) (scaled up by a factor of 3 to be consistent with reference Kotera et al. (2009)), with a spectral energy distribution resulting from the superposition of the emission of 100 giant elliptical galaxies, and a density profile following the ICM gas distribution. The infrared background of the cluster add to the CMB and the EBL as target radiation fields. While not only radio-loud AGN but also radio-quiet AGN have been suggested as UHECR accelerators Pe’er et al. (2009), the former objects are known to be more powerful. Radio-loud AGN are more common as a central galaxy in cool core clusters than in non-cool core clusters. To be conservative, we have not taken into account the high densities at the center of the subset of clusters that have cool cores Kotera et al. (2009).

Our simulations allow us to approximately calculate spectra of neutrinos and cosmic rays for a steady injection over a duration $t_{\rm inj}\sim 1-10$ Gyr. The history of injections from AGN is time-dependent and is dominated by their past activities at $z\sim 1-2$. Considering this effect, we have adopted $t_{\rm inj}=2$ Gyr. Two effects induce a spectral steepening in our model. First, a spectral break due to the escape of cosmic rays occurs at sufficiently high energies when $t_{\rm diff}(E)<t_{\rm inj}$ Murase et al. (2008, 2013). In addition, a break in the injection spectrum is caused by the confinement of cosmic rays in the cocoon. Figure S1 demonstrates that $t_{\rm inj}=t_{H}$, where $t_{H}=13.75~{}{\rm Gyr}$ is the age of the Universe, leads to similar results. Thus for our fiducial parameters, our results are insensitive to changes in $t_{\rm inj}$ which affects the confinement ability of the cluster. In essence, since we achieve $s<s_{\rm acc}$ (where $s$ is the effective injection spectral index), our model is consistent with model-independent bounds on cosmic-ray reservoir scenarios, in which the effective injection spectral index, $s$, is constrained by the $\gamma$-ray background to be smaller than $s\sim 2.1-2.2$ Murase et al. (2013).

All cosmic-ray ions and neutrinos, including the primary cosmic-ray particles injected into the simulation and their secondary and higher-order products, are tracked down to 1 TeV. The neutrino spectrum below this energy is extrapolated based on a spline fit to the simulation results between 2.2 TeV and 10 TeV. From the neutrino spectrum we obtain the spectrum of their $\gamma$-ray counterparts by Murase et al. (2016)

[TABLE]

where $EQ_{E}$ is the energy generation rate density per logarithmic energy, and $\gamma$-ray and neutrino energies are related as $E_{\gamma}\approx 2E_{\nu}$. Then we take into account electromagnetic cascades during the intergalactic propagation by solving transport equations Murase et al. (2012). Note that clusters and groups are expected to be transparent up to $\sim 100$ TeV energies Inoue et al. (2005); Murase and Beacom (2013). Either leptonic or leptohadronic, $\gamma$-ray emission can be observed directly from jets and lobes of radio-loud AGNs. In particular, the emission from blazars jets significantly contributes to the EGB. Primary relativistic electrons cool inside the AGN, and are not expected to diffusively enter the ICM like cosmic-ray ions.

The diffuse flux computation convolves the mass function of reservoirs (that are clusters and groups) and luminosity density of all possible cosmic-ray accelerators (that are mainly AGN jets) as shown by equation 1. Since it is assumed that all radio-loud AGNs, including BL Lac objects and their misaligned counterparts, can produce UHECRs (see reference Murase and Beacom (2013) for the discussion on the inner jet model for UHECRs), we use a luminosity density obtained by integrating the luminosity function. Note that the photon luminosity densities of radio galaxies, blazars, and AGN (rather than the number densities) show a positive redshift evolution Willott et al. (2001); Ajello et al. (2014); Ueda et al. (2014), and uncertainty in the redshift evolution can affect the diffuse fluxes only by a factor of $\sim 3$ (see reference Murase and Waxman (2016) for more details). Then, we simulate host clusters in 10 mass bins from $M_{\rm min}=5\times 10^{13}\,M_{\odot}$ to $M_{\rm max}=10^{16}\,M_{\odot}$, and 14 redshift bins from $z_{\rm min}=0.005$ to $z_{\rm max}=5$. The result, however, barely depends on $M_{\rm max}$ and $z_{\rm max}$ since the population of massive clusters at high redshifts is negligible. The minimum redshift $z_{\rm min}=0.005$ corresponds to the distance of the Virgo cluster at $\sim 20$ Mpc. The halo mass function is calculated by

[TABLE]

where $\rho_{m}(z)=\rho_{m}(0)\,(1+z)^{3}$ is the mean density of the universe at given redshift, and $\sigma(M,z)$ corresponds to the variance of the linear density field smoothed on a top-hat window function Tinker et al. (2008) $r_{t}=(3M/4\pi\rho_{m})^{1/3}$. For $f(\sigma)$, we adopt the Sheth-Tormen mass function multiplicity Sheth and Tormen (1999), which is consistent with results from N-body cosmological simulations Jenkins et al. (2001); Warren et al. (2006); Tinker et al. (2008).

A galaxy cluster and group naturally hosts supermassive black holes with powerful jets as cosmic-ray accelerators, which may include not only radio-loud AGN but also radio-quiet AGN Murase (2015). The emission from the accelerators themselves is a very different subject beyond the scope of this work. As noted above, synchrotron emission from primary electrons in the jets (or lobes) of blazars and radio galaxies have been observed. Since the injected cosmic-ray ion flux may be $\sim 100$ times higher than the electron flux, as indicated by the Galactic cosmic-ray data, it is also possible to have a hadronic component as often discussed in the leptohadronic model for blazars. Our model is consistent with either leptonic or leptohadronic model for the high-energy emission from the cosmic-ray accelerators, and the blazars may make a subdominant contribution to the PeV neutrino flux. This is a generic feature of cosmic-ray reservoir scenarios. In principle, cosmic-ray accelerators may even be powerful transients in galaxies including $\gamma$-ray bursts, supernova explosions, fast-rotating pulsars, tidal disruption events Waxman (2013); Senno et al. (2015), and blazar flares Dermer et al. (2014). In addition, particle acceleration can also happen in the accretion shocks that are commonly found at the outskirts of large-scale structures Kang et al. (1996); Inoue et al. (2007). The accretion shock contribution is expected to be more important for massive clusters Murase et al. (2008); Fang and Olinto (2016); Murase and Waxman (2016), and the corresponding neutrino production has been subject to various constraints Kushnir et al. (2009); Fang and Olinto (2016); Zandanel et al. (2015); Murase and Waxman (2016). In contrast, the contribution from internal accelerators such as AGN can be more important for less massive clusters Murase et al. (2008); Kotera et al. (2009); Murase et al. (2013), since the population of such clusters is higher and AGN are found to have a strong redshift evolution. The distribution of higher-energy cosmic rays is expected to be more uniform Keshet (2010).

For the extragalactic propagation, we use CRPropa 3 and take into account the photomeson and photodisintegration interactions between cosmic rays and background photons including the CMB and the EBL (where the best-fit model of reference Kneiske et al. (2004) is used), as well as the decay of unstable intermediate products and the Bethe-Heitler production of electron-positron pairs. We use the transmission factor Kotera and Lemoine (2008) to take into account the magnetic horizon effect due to structured EGMFs Lemoine (2005); Das et al. (2008); Alves Batista and Sigl (2014). The confinement in the structured EGMFs may also lead to additional secondary production especially below the ankle, which is not included to save computation time. Cosmogenic neutrinos Berezinsky and Zatsepin (1969) are generated by photomeson interactions of protons above the Greisen-Zatsepin-Kuzmin energy Greisen (1966); Zatsepin and Kuz’min (1966) ($E_{\rm GZK}\sim 6\times 10^{19}$ eV) during their intergalactic propagation. The flux of cosmogenic neutrinos is lower if UHECRs are dominated by intermediate or heavy nuclei, as opposed to pure protons Murase and Beacom (2010); Kotera and Olinto (2011); Ahlers and Halzen (2012). Cosmogenic $\gamma$ rays from intermediate and heavy nuclei are dominated by the Bethe-Heitler pair production Ahlers and Salvado (2011); Decerprit and Allard (2011); Murase et al. (2012). It is known that cascaded $\gamma$ rays follow a universal spectral shape, and that their flux is basically determined by the energy injection rate of cosmic rays at ultrahigh energies Berezinsky and Smirnov (1975); Wang et al. (2011); Murase et al. (2012). For demonstration purposes, the cluster contribution to the EGB is estimated by adding the cosmogenic $\gamma$-ray flux of a mixed composition scenario Decerprit and Allard (2011) to the source $\gamma$-ray flux obtained in this work. As in the simplest convergence theory Murase and Waxman (2016), the total $\gamma$-ray flux from our model is compatible to the non-blazar component of the EGB Ackermann et al. (2016); Lisanti et al. (2016).

The combined fit to the Auger spectral data, the $\langle X_{\rm max}\rangle$ data, and the IceCube spectral data is performed as follows. The goodness-of-fit is assessed with a total $\chi^{2}$ contributed by three separate fits, given that the UHECR energy spectrum, the $X_{\rm max}$ distributions, and the neutrino energy spectrum are independent measurements:

[TABLE]

$C_{\rm norm}$ is a universal factor that normalizes the predicted fluxes of UHECRs and neutrinos, $\delta_{E}^{\rm CR}$ and $\delta_{E}^{\nu}$ denote the energy-scale displacements in the UHECR and neutrino spectra, defined so that $E\rightarrow E^{\prime}=(1+\delta_{E})\,E$. The allowed displacements are capped by the systematic energy uncertainties in the Auger spectrum, $\Delta_{E}^{\rm CR}=\pm 14\%$ Aab et al. (2015), and the deposited energy resolution of cascade events of IceCube The IceCube Collaboration et al. (2015), $\Delta_{E}^{\nu}\sim\pm 15\%$. We do not allow an energy-scale displacement in the fits to the $\langle X_{\rm max}\rangle$ data since the energy resolution of the Fluorescence Detector (FD) of Auger is narrower than the width of the energy bins used here Aab et al. (2017). Specifically, the $\chi^{2}$ functions are defined as (see e.g., reference Baerwald et al. (2015)):

[TABLE]

Here $\Phi_{\rm CR}$ and $\Phi_{\nu}$ are the flux of UHECRs measured by the Auger Surface Detector (SD) Aab et al. (2015) and the flux of astrophysical neutrinos in the six-year IceCube data Aartsen et al. (2017a), $\langle X_{\rm max}\rangle$ is the mean of $X_{\rm max}$ measured by the Auger FD Aab et al. (2014, 2015), and $\Delta\left(\Phi_{\rm CR}\right)$, $\Delta\left(\Phi_{\nu}\right)$, and $\Delta\langle X_{\rm max}\rangle$ are the corresponding uncertainties defined via the quadrature sum of the statistical and systematic uncertainties of the measurements. We use $E^{3}\Phi_{\rm CR}$ and $E^{2}\Phi_{\nu}$ to evaluate the fits to cosmic-ray and neutrino data respectively as these are the forms of data provided by the corresponding experiments. The systematic uncertainty of the UHECR spectrum is derived from that of vertical Auger SD data sets recorded by the 750 m and 1500 m arrays, as well as inclined events recorded by the 1500 m array Aab et al. (2015), weighted by the event number from each data set in each energy bin. In addition, as suggested by a comparison of the TA and Auger spectra Verzi et al. (2017), we increase the energies of Auger spectrum by 5% (so that $-19\%<\delta_{E}^{\rm CR}<9\%$) to determine the absolute flux of UHECRs. The data used in the fit consists of the Auger spectrum in 18 bins with 0.1 increment of $\log_{10}(E/\rm eV)$ from 18.45 to 20.15, $\langle X_{\rm max}\rangle$ in 12 bins with 0.1 increment from 18.45 to 19.50, and one last bin between 19.50 and 20.00, and the IceCube spectrum of high-energy starting events above $2\times{10}^{14}$ eV. A total of 33 non-zero data points are used for the fit. When the upper and lower errors of a data point are different, the upper error is used when the model prediction is higher than the data, and vice versa.

The default parameters are set to be $E_{\rm max}=2\times 10^{21}/Z$ eV and $s_{\rm acc}=2.3$ as motivated by black hole jets of radio-loud and radio-quite AGN Murase et al. (2012); Pe’er et al. (2009). Our results are however not very sensitive to small variations of these parameters. Generally a softer cosmic-ray spectrum also leads to a softer neutrino spectrum, and a larger $E_{\rm max}$ leads to higher flux of cosmogenic neutrinos. We inject a chemical composition based on the abundance of elements in Galactic cosmic rays. Specifically, we adopt the elemental abundances at PeV$/Z$ based on cosmic-ray measurements from 10 GeV to 100 PeV Hörandel (2003). The iron abundance is enhanced by a factor of 2, since elliptical galaxies may have higher metallicity than star-forming galaxies such as the Milky Way. The energy flux ($E^{2}dN_{\rm inj}/dE$) at the same rigidity scales to (0.625 H, 0.252 He, 0.053 CNO, 0.009 Si, 0.124 Fe), which is also compatible with the cosmic-ray composition at 100 PeV. By minimizing $\chi^{2}$ through the fits with the Auger and the IceCube data, we find a best-fit $\chi^{2}_{\rm dof}=44.5/30=1.48$ (corresponding to a non-chance occurrence probability of 4.3% for the fiducial case to explain the data) at $\hat{\delta}_{E}^{\rm CR}=-0.12$ and $\hat{\delta}_{E}^{\nu}=-0.01$. The best-fit $C_{\rm norm}$ corresponds to $f_{\rm jet}=23.7\%$, assuming $L_{\rm CR}=10^{45}\,\rm erg\,{\rm s}^{-1}$ above 10 GeV. We obtain similar results, $\chi^{2}_{\rm dof}\sim 1.5-2$, if changing the starting energies to $10^{18.45\pm 0.2}$ eV in the UHECR fits and $6\times{10}^{13}$ eV in the neutrino fit. The fit could be improved to $\chi^{2}_{\rm dof}=1.3$ with the non-blazar component of the Fermi data above 50 GeV, although more careful analyses are required to properly estimate the uncertainties of the $\gamma$-ray data. A fit with $\chi^{2}_{\rm dof}\sim 3-5$ is obtained without applying the EGMF. A more dedicated fitting allowing arbitrary abundance, $E_{\rm max}$, $s_{\rm acc}$, and varying EGMF parameters can potentially provide a better fit to data Aab et al. (2017); Wittkowski (2017), but such a study is beyond the scope of this work. Given that there are uncertainties in the acceleration and propagation models as well as photonuclear interactions, our work is sufficient for demonstrating that the physical connection among multi-messengers is achieved by a concrete model with the physically-motivated source parameters.

The connection between PeV neutrinos and UHECRs in cosmic-ray reservoir scenarios was suggested before the discovery of IceCube neutrinos Murase et al. (2008); Kotera et al. (2009), and their link to the diffuse isotropic $\gamma$-ray background has also been discussed in references Murase et al. (2013); Murase and Waxman (2016); Wang and Loeb (2017). Our model is the first concrete model that simultaneously explains the measured diffuse fluxes of all three messengers without violating observational constraints. The model avoids ad hoc tuning of parameters, in the sense that the source population (including the number density and redshift evolution), the profiles of target gas and radiation, and the injection composition are all determined by observations, and the choice of free parameters ($E_{\rm max}$ and $s_{\rm acc}$) is physically motivated. Our numerical approach takes into account impacts of magnetic structures within sources and contributions from secondary particles. Note that non-flaring blazar models Murase et al. (2014); Giacinti et al. (2015) as the dominant sources of IceCube neutrinos are disfavored by the stacking analyses Aartsen et al. (2017b) and multiplet analyses Murase and Waxman (2016) of the IceCube data. High-redshift sources or rapidly evolving sources such as quasars Wang and Loeb (2017); Xiao et al. (2016) do not explain the Fermi data above 10 GeV due to the strong EBL absorption. Unlike the blazar and high-redshift source models, our results respect the updated constraints on the non-blazar component of the Fermi EGB Ackermann et al. (2016) with the important contribution from cosmogenic $\gamma$ rays Murase and Waxman (2016), the neutrino flux above 10 PeV Aartsen et al. (2016b) by including cosmogenic neutrinos, and the heavy-rich composition of UHECRs indicated by Auger Aab et al. (2015), as well as the sub-ankle component above ${10}^{17}$ eV.

Light composition scenario

So far we have mainly focused on a mixed composition scenario, but the interpretation of the chemical composition of UHECRs is under debate. While three interaction models, EPOS Pierre Auger Collaboration (2013), QGSJET Ostapchenko (2008) and SIBYLL Ahn et al. (2009) are commonly used, hadronic interaction models at the highest energies are based on the extrapolation from experimental data at lower energies, and thus still suffer from non-negligible uncertainties. A simultaneous fit to the Auger spectrum, $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$ suggests a dominance of intermediate or heavier nuclei at the injection site Taylor et al. (2011); Aab et al. (2017), which is quite challenging for known astrophysical sources. References Aab et al. (2017); Wittkowski (2017) find a main minimum and a second minimum in their fitting. The main minimum requires a hard injection spectral index $s=0.9$ (or 1.61 with the EGMF), which is different from typical expectations from the diffusive shock acceleration mechanism (although such a hard spectrum could be reached by the shock acceleration at ultrarelativistic shocks Meli et al. (2008); Aoi et al. (2008) or other mechanisms such as the shear acceleration Kimura et al. (2018)). The second minimum suggests $s=2.0$ (or 2.3 with the EGMF) and is similar to the mixed composition scenario shown here, but with a $\sigma(X_{\rm max})$ lighter than the Auger measurements between $10^{18.8}$ eV and $10^{19.4}$ eV. The UHECR composition can also be modified by interactions inside cosmic-ray accelerators Murase et al. (2012); Fang et al. (2012); Murase et al. (2014). More quantitative discussion on the composition is beyond the scope of this work, and we do not attempt to fit to the $\sigma(X_{\rm max})$ data.

Below we consider an alternative scenario, where cosmic rays from radio-loud AGN are dominated by pure protons (aka the “proton scenario”). Such a possibility is shown in Figure S2 for the completeness of our demonstration, since it has been argued that $X_{\rm max}$ measurements by TA may be consistent with the light composition Charles et al. (2015). Except for a different composition, we adopt the same injection spectrum, $dN_{\rm inj}/dE\propto E^{-2.3}$ and the same maximum energy, $E_{\rm max}=10^{19.6}$ eV, as in our mixed composition case. Like in the mixed-composition scenario, the neutrino products from the ICM is comparable to the $0.1-1$ PeV measurements by IceCube. The difference between a nucleus and a proton in neutrino production can be very roughly estimated as follows. Considering that the inelastic cross section of a nucleus with mass number $A$ scales as $\sigma_{Np}\sim A^{2/3}\,\sigma_{pp}$, and that a nucleus with $Z\times E$ experiences the same level of diffusion as a proton with $E$, the interaction between the nucleus and the ICM gas has an effective optical depth $f_{Np}(ZE)\sim(\kappa_{Np}/\kappa_{pp})A^{2/3}\,f_{pp}(E)$, where $\kappa_{Np}$ and $\kappa_{pp}$ are inelasticities of the projectile particle in each case. Therefore, a neutrino product from a nucleus and that from a proton with $Z$ times lower energy have an energy ratio $Z/A\sim 0.5$ and an energy flux ratio $\sim(\kappa_{Np}/\kappa_{pp})A^{2/3}$.