Constraining Baryon Loading Efficiency of Active Galactic Nuclei with Diffuse Neutrino Flux from Galaxy Clusters

Active galactic nuclei (AGNs) are widely believed to be one of the promising acceleration sites of ultrahigh-energy cosmic rays (CRs). Essentially, AGNs are powered by the gravitational energy of matter falling into supermassive black holes. However, the conversion efficiency of gravitational to kinetic energy of CRs in AGNs, which is defined as the baryon loading factor η p , is not well known yet. After being accelerated, high-energy CRs could escape the host galaxy and enter the intracluster medium (ICM). These CRs can be confined within the galaxy cluster and produce γ-rays and neutrinos through proton–proton collisions with the ICM. In this paper, we study the diffusion of CRs in galaxy clusters and calculate the diffuse neutrino flux from the galaxy cluster population. Using the latest upper limits on the cumulative unresolved TeV–PeV neutrino flux from galaxy clusters posed by the IceCube Neutrino Observatory, we derive the upper limit of the average baryon loading factor as η p,grav ≲ 2 × 10−3 − 0.1 for the population of galaxy clusters. This constraint is more stringent than the one obtained from γ-ray observation on the Coma cluster.


Introduction
Active galactic nuclei (AGNs) are the most powerful persistent emitters of radiation in the Universe and have been considered as a potential source of extragalactic high-energy cosmic rays (CRs ;Biermann 1988;Takahara 1990;Rachen & Biermann 1993;Berezinsky et al. 2006;Dermer et al. 2009) and neutrinos (Mannheim et al. 1992;Stecker & Salamon 1996;Atoyan & Dermer 2001;Murase et al. 2014).Acceleration of baryonic CRs in jets/outflows of AGNs consumes the kinetic energy or magnetic energy of the jets/outflows, which are essentially fueled by the gravitational energy of matter falling into a supermassive black hole at the center of the nuclei.The efficiency of gravitational energy converting into CRs, which is defined as the baryon loading factor η p , can help us understand the physical mechanism of particle acceleration.
The IceCube Neutrino Observatory has been observing TeV-PeV astrophysical neutrinos for over one decade.All the discovered potential sources are associated with AGNs, such as the blazar TXS 0506+056 (IceCube Collaboration et al. 2018a, 2018b) and the Seyfert II galaxy NGC 1068 (Aartsen et al. 2020) detected by IceCube.Indeed, numerous studies have shown that protons accelerated in AGNs can interact with their intense radiation fields and produce high-energy neutrinos via the photohadronic interactions (Rachen & Mészáros 1998;Atoyan & Dermer 2003;Stecker 2013) or hadronuclear interactions (Fraija et al. 2012;Sahakyan et al. 2013;Li et al. 2022;Xue et al. 2022).On the other hand, whether the all-sky diffuse neutrino background can be accounted for by the AGN population highly depends on the average baryon loading factor of the AGN population, because the expected neutrino flux is proportional to this parameter.For example, Murase et al. (2014) found that blazars may account for the diffuse neutrino background above 100 TeV with η p,rad ≡ L p /L γ ∼ 3-300 assuming an E p 2 -CR proton spectrum.However, the baryon loading mechanism of AGN jets or outflows is not well known yet, which prevents us from drawing a concrete conclusion on the contribution of AGNs to the all-sky diffuse neutrino background (Berezhko 2008;Cuoco & Hannestad 2008;Kadler et al. 2016;Righi et al. 2017;Palladino et al. 2019).
On the other hand, CRs generally lose only a small fraction of energies through the hadronic interactions in AGNs (e.g., Murase et al. 2014;Xue et al. 2019).Thus, they may eventually escape AGNs and host galaxies after being accelerated, and propagate into the intracluster medium (ICM) of the galaxy cluster.Galaxy clusters are the largest gravitational-bound structures in the Universe.They are also known as efficient reservoirs for CRs (Völk et al. 1996;Berezinsky et al. 1997; see Brunetti & Jones 2015 for a recent review).The diffusion timescale of CRs with energies 1 PeV is longer than the Hubble timescale with a reasonable diffusion coefficient in the ICM.The confined CRs would interact with the ICM via pp collision and produce γ-rays and neutrinos.Therefore, measurements on γ-rays and neutrinos from galaxy clusters can serve as constraints on the amount of CRs accelerated in AGNs, which can be translated to the baryon loading factor.
The γ-ray emissions from individual galaxy cluster has been searched in very-high-energy band (>100 GeV) by Fermi Large Area Telescope (LAT; Han et al. 2012;Ackermann et al. 2016) and Imaging Air Cherenkov Telescopes (Perkins et al. 2006;Aleksić et al. 2010Aleksić et al. , 2012;;Arlen et al. 2012) for a long time.Recently, extended GeV γ-ray emission from the direction of the Coma cluster has been reported (Xi et al. 2018;Adam et al. 2021;Baghmanyan et al. 2022).The former estimation suggested that NGC 4869 and NGC 4874, which are the two brightest radio galaxies in the radio band of the cluster, Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.cannot account for the entire observed γ-ray emission (Ackermann et al. 2016;Baghmanyan et al. 2022).However, it still remains uncertain whether the emission is from the diffuse CRs in the ICM or a combination of several unresolved sources in the region.Under the assumption that the γ-ray signal comes entirely from the decay of π 0 produced in pp collisions, the observed γ-ray flux can set an upper limit for the CR content in a galaxy cluster.In our previous work (Shi et al. 2022), we studied the propagation of CRs in the Coma cluster and calculated the radial distribution of generated pionic γ-rays emission.By comparing the γ-ray flux and upper limits obtained by Fermi-LAT and VERITAS for the Coma cluster, we have established an upper limit on the average baryon loading factor for AGNs in the cluster as η p,grav = W p /W grav  0.1 (or η p,rad  1).This limit is found to be lower than the baryon loading factor required for blazars as obtained by Murase et al. (2014).
However, we note that the constraint on the baryon loading factor of AGNs obtained from the Coma cluster may not be generalized to the entire AGN population in the Universe.A more representative constraint would be based on the all-sky diffuse neutrino flux or all-sky γ-ray flux.Recently, Abbasi et al. (2022) performed stacking analysis of 1094 galaxy clusters from Planck using 9.5 yr of muon-neutrino track events and found no evidence for significant neutrino emission.The differential upper limits presented by IceCube in Abbasi et al. (2022) are most constraining in the energy range between 10 TeV and 1 PeV, suggesting that the contribution of the galaxy cluster population cannot exceed 9%-13% of the diffuse neutrino flux.We may calculate the expected flux of diffuse neutrinos produced by CRs escaping from AGNs via pp collisions in ICM.By comparing the expected flux with the measured upper limit, a constraint on the average baryon loading factor of AGN jets or outflows can be obtained.
The rest of this paper is organized as follows.In Section 2, we first review the propagation of CRs escaped from AGN in the galaxy cluster, and then calculate the high-energy neutrino production from galaxy clusters through interactions between these CRs and ICM.In Section 3, we compare the expected diffuse neutrino flux with the upper limit given by IceCube to constrain the amount of injected protons from AGN, which can be translated to the baryon loading factor.We summarize and discuss our results in Section 4.

Model
In this section, we describe the model we used in our work.In Section 2.1, we review particle diffusion and examine the confinement of CRs in the turbulent magnetic field of galaxy clusters.Then, in Section 2.2, we calculate the neutrino emissivity from an individual galaxy cluster.Finally, in Section 2.3, we integrate our results and obtain the diffuse neutrino flux contributed by AGNs in galaxy cluster populations.

Cosmic-Ray Diffusion
After acceleration, particles propagate through the turbulent magnetic field of the cluster.Their diffusion depends on both the particle's Larmor radius r L and the coherence length l c of the magnetic field.For typical parameters in the galaxy clusters, when r L < l c , corresponding to the particle's energy E p  5 × 10 20 Z(B/5 μG)(l c /0.1 Mpc) eV, the propagation of the CRs enters the diffusive regime with the diffusion coefficient of For Kolmogorov diffusion, the spectral index is w = 5/3; we assumed that B ∼ δB and l c ∼ 10% of the virial radius r vir,200 is the typical magnetic field coherence length in the galaxy clusters.The virial radius r vir,200 of the cluster with mass M vir is defined as , where Δ c = 200.The diffusion timescale of the CRs in a galaxy cluster with mass M vir = 10 15 M e can be estimated as ´-When diffusion timescale is longer than the Hubble time t H ∼ 14 Gyr, the CRs with energy E p  0.9 PeV are confined in the cluster by magnetic fields.
Then, we calculate the radial density distribution of the CRs at present in a typical galaxy cluster.As different CR injection histories have a slight effect on the result (Shi et al. 2022), we assume a time-independent CR injection rate with a single power-law spectrum of injection index α, Neglecting the energy loss of particles, the radial density distribution of the CRs at present can be written as We assume that injection occurs at the center of the galaxy cluster (r = 0) and that the injection duration lasts for ∼10 Gyr.

Neutrino Production
As CRs propagate through the cluster magnetic field, they interact with the ICMs and produce γ-ray photons and neutrinos.In this section, we calculate the neutrino flux from a typical galaxy cluster.The hot ICMs emit X-rays via bremsstrahlung radiation, with the emissivity proportional to the square of the number density of electrons in the gas (n ICM,e ).Therefore, the density of ICM can be derived from X-ray observations.From stacking the Chandra data of 320 clusters, the mean density profile of electrons in the ICM can be approximated with a form introduced by Patej & Loeb (2015): is the Hubble parameter at redshift z, today's Hubble parameter is referred to as the Hubble constant, H 0 .In a fully ionized gas, the number density ratio of electron and proton is n ICM,e = 1.17nICM,p .
Following the calculation in Kelner et al. (2006), we first calculate the γ-ray emissivity as where σ pp (E p ) is the total inelastic cross section of pp interactions, F γ (E γ /E p , E p ) is the spectrum of the secondary γ-ray in a single collision.
Due to the advantage of the angular resolution of muon track events, source analyses presented with IceCube usually focus on (anti)muon neutrinos.Assuming equal amount of neutrinos of three flavors after the oscillation, we can relate the muonneutrino emissivity to that of γ-rays as . 6 Integrating the total extent of the galaxy cluster, the total muonneutrino flux from an individual galaxy cluster can be calculated as where dV r drd d sin 2 q q j = is the differential volume element of coordinate r = (r, θ, j), ψ is the angular extension of the source, and d L is the luminosity distance of the cluster.
The virial mass of a galaxy cluster can affect the neutrino emissivity in our model.To study this influence, we compare the normalized radial distribution of the neutrino emissivity at E ν = 1 TeV for a cluster of different virial masses.First, the energy budget of CRs in a more massive cluster is higher if the same fraction of gravitational energy is converted to CRs (i.e., the same η p,grav ).Also, a more massive cluster is surrounded by a larger amount of ICM, leading to a higher gas density compared to the same radius r.It results in a higher interaction rate of pp collisions and consequently increases the neutrino emissivity.On the other hand, the magnetic field coherence length l c is related to the virial mass by l r M 0.1 The resulting diffusion coefficient is positively related to the virial mass and hence particles diffuse faster in more massive clusters, leading to a flatter radial distribution of the neutrino emissivity in a more massive cluster than in a less massive one.Note that the virial radius increases with the virial mass as M vir 1 3 , and hence the diffusive escape timescale increases with the virial mass as shown in Equation (2).Therefore, even with the same CR injection rate, the total neutrino luminosity still increases with virial mass.The comparison of the radial distribution of the neutrino emissivity between a 10 14 M e cluster and a 10 15 M e cluster is shown in the left panel of Figure 1.A direct dependence of the total neutrino luminosity on the virial mass of the cluster is shown in the right panel of Figure 1.We observed a monotonic increase of the neutrino luminosity with the virial mass, approximately with a linear relation.
To ensure consistency with the results from Abbasi et al. (2022), we integrate the total neutrino flux over the extent of the cluster from 0 to the virial mass r vir,500 , where Δ c = 500, when calculating the total neutrino flux.As the neutrino flux is concentrated in the central region of the galaxy cluster, the integration range has little effect on the final results.

Diffuse Neutrino Flux
The diffuse muon-neutrino flux integrated from the entire galaxy clusters can be estimated as where the differential number density dn/dM of clusters with mass M at redshift z can be obtained from the halo mass function and ρ m is the mean density of the Universe at the epoch of analysis, ( ) is the rms variance of the linear density field smoothed on scale , and f (σ) describes the σ-weighted distribution. ( = +is the angular diameter distance.In this study, we adopted the same mass halo function from Tinker et al. (2010) as used in Abbasi et al. (2022) to ensure consistency.As galaxy clusters with masses below 10 14 M e or z > 1 are not expected to produce a significant flux of neutrinos at Earth (Fang & Olinto 2016), our calculation only considers clusters with masses between 10 14 M e and 10 15 M e and a redshift between 0.01 and 2.

Constraining the Baryon Loading Factor
In this section, we calculate the average baryon loading factor η p,grav for the population of galaxy clusters.Although the value of η p,grav may vary for each cluster, our primary concern lies in the constraints on the total population of galaxy clusters.Therefore, we focus on calculating the average value of η p,grav for the entire cluster population.We use the same definition of baryon loading factor η p,grav in Section 4.1 from Shi et al. (2022) to constrain the total efficiency of releasing gravitational potential energy W g loading into total injected energy of the baryons W p , written as where the total injected energy of the baryons W p can be calculated from the total CR injection rate ò .The releasing gravitational potential energy W g is estimated as W g = 0.2M BH c 2 with an intermediate mass-to-energy conversion efficiency between standard accretion disk model ≈0.1 (Shakura & Sunyaev 1976) and an extreme Kerr black hole ≈0.3 (Thorne 1974).The total black hole mass M BH in the cluster can be estimated as a fraction of the virial mass M vir , M BH = η BH M vir .Having noted that, this fraction is linearly correlated to the total gravitational energy W g and hence would linearly affect the obtained value of η p,grav .For the Coma cluster, the galaxy mass can be integrated with the help of the mass-to-light (V-band) ratio, M gal = 2.03 × 10 13 M e while the virial mass is measured to be M vir,500 = 6.0 × 10 14 M e .Therefore, the ratio of galaxy mass to the virial mass is fixed at η gal = M gal /M vir = 0.034 in our study, which is also consistent with the standard ΛCDM model.
Abbasi et al. (2022) obtained the neutrino flux upper limit from the galaxy cluster population for two different weighting methods for the expected neutrino flux from a galaxy cluster.One is the so-called distance weighting, assuming neutrino luminosity is the same among all clusters and hence the neutrino flux is proportional to d 1 L 2 .The other is the mass weighting, assuming the neutrino luminosity scales linearly with the virial mass of the cluster and hence the neutrino flux from an individual cluster scales with M d L vir 2 .According to our discussion in the previous section and the result shown in Figure 1, the mass weighting is more consistent with our model.However, in the analysis, they fixed the neutrino spectrum to be an unbroken power-law function with a slope of −2.5.As we want to test different values of the spectral index α for the injected protons, it is not appropriate to directly compare our results with the flux upper limit derived with the mass weighting method.On the other hand, Abbasi et al. (2022) also provided the quasi-differential flux upper limits of 90% confidence level in one-decade energy bins with the distance weighting, we choose to constrain the value of η p,grav by comparing our results with this differential upper limit.Note that the neutrino flux upper limit obtained with the mass weighting is stricter than that obtained with the distance weighting.So, such a comparison leads to conservative constraints on the baryon loading factor η p,grav .
Using the same cluster mass ranging from 10 14 M e to 10 15 M e and the redshift between 0.01 and 2, we calculate the results and compare them with the differential upper limits from Abbasi et al. (2022) in Figure 2. Since the modeled neutrino flux linearly depends on η p,grav , we can find out the maximally allowed value of η p,grav that makes the modeled neutrino flux saturating only one of the five bins for the flux upper limit.For reference, we also display the isotropic diffuse γ-ray background (IGRB) observed by Fermi-LAT (Ackermann et al. 2015) in the figure.Furthermore, we explore the maximum η p,grav with different CR spectral injection index α, and obtain the relation between α and the upper limit of η p,grav as shown in Figure 3.For a flat injection CR spectrum α = 2, which is expected under the canonical shock acceleration theory, a quite strict constraint η p,grav = 2 × 10 −3 can be obtained.The constraint becomes less stringent for a softer spectral index, and the upper limit of η p,grav increases up to 0.12 for α = 2.5.
We note that the diffuse neutrino flux from the galaxy cluster highly depends on model parameters such as the baryon loading factor η p,grav and others.In the previous literature (Fang & Olinto 2016;Fang & Murase 2018;Hussain et al. 2023), the authors aim to explore the potential of galaxy clusters as the major sources of the all-sky diffuse high-energy neutrino flux, so they tune the model parameters to make the predicted neutrino flux match the measured one.The latest observational upper limits of the diffuse neutrino flux from galaxy clusters presented by IceCube in Abbasi et al. (2022) suggest that the contribution of the galaxy cluster population with masses between 10 14 M e and 10 15 M e cannot exceed 9%-13% of the diffuse flux, which can thus be also used to constrain the model parameters in previous literature, as we do in the present study.On the other hand, it should be noted that the IceCube's result does not necessarily rule out the entire galaxy cluster population as the main contributor of the diffuse neutrino background, because a galaxy cluster of smaller mass could make an important contribution.For example, in our model, the diffuse neutrino flux emitted by a galaxy cluster with a mass below 10 14 M e is comparable to that emitted by a galaxy cluster with a mass between 10 14 M e and 10 15 M e at 100 TeV.Adding this part would make the contribution of the entire galaxy cluster population reach 30%-40% of the all-sky flux at most.While this still suggests galaxy clusters to be subdominant contributors in our model, the ratio of the neutrino flux from low-mass galaxy clusters to that from high-mass galaxy clusters could vary from model to model.

Cumulative γ-Ray Flux
The CR-ICM interaction can also produce γ-ray photons via the decay of π 0 s.There are suggestions that the observed diffuse neutrino/γ-ray flux could be completely explained by the cumulative emission from galaxy clusters Figure 1.Left: the normalized radial density distribution of neutrino emissivity with energy E ν = 1 TeV for a cluster with virial mass M vir = 10 14 M e (red), 10 15 M e (green), and redshift z = 0.The flux is normalized by fixing it to unity at the center of the cluster with M vir = 10 14 M e .The total neutrino emissivity derived from l c = 0.1r vir,200 under different CR injection rates, regarding either the same total CR injection W const p = (dashed line) or a constant baryon loading factor const p,grav h = (solid line).The latter is the model employed in our study.Right: the relation between the virial mass of a galaxy cluster and the total neutrino emissivity at redshift z = 0.The emissivity is normalized by fixing the value to be unity for the cluster of M vir = 10 14 M e .The line types correspond to the same assumptions as in the left panel.
Figure 2. Integrated muon-neutrino flux (red lines) from galaxy clusters in this work.Different injection spectral index α = 2 (dashed), 2.2 (solid), and 2.5 (dotted) are illustrated separately.It is compared with the differential upper limits in one-decade energy bins for the distance weighting (1/d 2 ) scheme from IceCube analysis (Abbasi et al. 2022).The IGRB observed by Fermi-LAT (Ackermann et al. 2015) is also shown for comparison.
Figure 3.The upper limits of baryon loading factor η p,grav constrained by diffuse muon-neutrino and γ-ray flux with different indices α of the injected protons from 2.0 to 2.5.(Fang & Murase 2018;Hussain et al. 2023).However, the origin of the IGRB is still under debate; besides the galaxy clusters, the observed IGRB is possibly superimposed by different populations of γ-ray emitters, such as star-forming galaxies, starburst galaxies, and AGNs.
We calculate the diffuse γ-ray flux from the clusters of galaxies in the same way as for neutrinos, taking into account the effect of extragalactic background light (EBL) attenuation using the model from Saldana-Lopez et al. (2021).The dominant contribution to the total flux of γ-rays comes from sources at low redshifts (z  0.3), where the effect of EBL attenuation is less pronounced.
Considering the IGRB observed by Fermi-LAT as an upper limit, we can also derive the upper limit of the baryon loading factor from the integrated γ-ray flux.Figure 4 compares the diffuse flux of γ-ray from the clusters of galaxies and the IGRB measured by Fermi-LAT.The upper limits of neutrino flux obtained by IceCube is also shown for reference.For considered injection indices, i.e., 2 α 2.5, the constraints derived from the diffuse γ-ray flux are 1-2 orders of magnitude more stringent than those from neutrinos.
However, there may be some potential uncertainties in the constraints obtained from the diffuse γ-ray flux.This is mainly because the corresponding energies of CRs responsible for these γ-ray emissions are relatively low.Shi et al. (2022) suggested that the propagation of these CRs may be influenced by the streaming instability (Kulsrud & Pearce 1969;Skilling 1971).As a result, the CR spatial distribution may be different from what is predicted in our current model and the CR energies may be dissipated through self-excited Alfvén waves in the ICM.In addition, these relatively low-energy CRs may not be so easy to escape their acceleration sites and could lose energy adiabatically (Fang & Murase 2018).These processes are not considered in our model.On the other hand, we should also note that the most constraining energy bin of the Fermi-LAT data is the highest-energy one in [580, 820] GeV (as shown in Figure 4), corresponding to the CR proton energy of 10 TeV.As a result, the aforementioned uncertainties may not be so severe.Another uncertainty arises from the accuracy of the EBL model, which may affect the attenuated γ-ray spectrum (Hussain et al. 2023).

Distinguish from the Acceleration Caused by Cluster Mergers
In galaxy clusters, CRs can also be accelerated by shock waves arising from cluster merger processes.These shocks would also accelerate particles to relativistic energies and produce γ-ray and neutrino emission (Colafrancesco & Blasi 1998;Ryu et al. 2003).The gravitational energy released by two galaxy clusters with a total mass of 10 15 M e merging from an infinity distance to 1 Mpc is ∼10 64 erg.Most of this gravitational energy is converted into the kinetic energy of dark matter; only about 10% is dissipated into the ICM.The gravitational energy released from central black hole accretions in our study is estimated as ∼10 64 erg for a galaxy cluster with virial mass 10 15 M e .Due to the lack of knowledge of the baryon loading efficiency of these two acceleration mechanisms, it is hard to tell which process dominates CR accelerations in the galaxy cluster.
The future observed γ-ray morphology may help us distinguish which process would dominate the acceleration.In the central AGN injection model, the γ-ray profile is a halolike structure that clusters in the center region and decreases with radius.Numerical simulations modeling the formation of large-scale structures show the presence of strong accretion shocks in the outer regions of galaxy clusters (Miniati et al. 2000;Ryu et al. 2003;Vazza et al. 2012); the radial density profile of CRs and γ-ray predicted by the cluster mergers would be a shell-like structure and significantly different from the central injection model.Studies have shown that the CR-to-thermal pressure ratio from the injection of CRs at cosmological shocks should slightly increase with radius (Vazza et al. 2012).However, in our model considering the CR injection from central AGNs, this ratio decreases with the radius to the cluster center unless the energy of the injected proton exceeds ∼1 PeV.The CR-tothermal energy ratio has been limited to 4%-10% for photon indices α = 2.0-3.2 using stacked Fermi-LAT count maps in the 1-300 GeV band (Huber et al. 2013).In our work, the upper limit on the average CR-to-thermal pressure ratio for proton indices of 2.0-2.5 (corresponding to a photon index of ∼2.2-2.7)from AGN contributions is estimated to be 0.3%-1.2%,which is more stringent than (and consistent with) the limitation obtained from stacking Fermi-LAT count maps.
Moreover, the radial density profile of CRs accelerated by AGNs is clustered in the central region of the galaxy cluster, where the gas density of ICMs is higher.This enhances the overall interaction rates of the pp collision and the neutrino/γray production efficiency.It should be noted that our main purpose is to constrain the upper limit of the baryon loading efficiency.Therefore, considering the γ-ray emission additionally contributed by other processes would not invalidate the results, but would only make the constraints tighter.

Conclusions
In this work, we calculated the diffuse neutrino and γ-ray flux from the galaxy cluster population and compared the results with the upper limit using 9.5 yr of muon track IceCube data and the IGRB observed by Fermi-LAT.In order to compare with the upper limits of the diffuse neutrino flux presented by Abbasi et al. (2022), the same mass range of galaxy clusters from 10 14 M e to 10 15 M e and the redshift between 0.01 and 2 are considered in our calculations.Our best constraint for the upper limits of the average baryon loading factor is η p,grav  2 × 10 −3 , assuming a flat injection spectrum index of CRs (α = 2).When varying the injected power-law spectrum index from 2 to 2.5, we derive the upper limit of the average baryon loading factor as η p,grav  0.1.The constraints using the IGRB observed by Fermi-LAT are about 1-2 orders of magnitude stricter than the constraints derived from the diffuse neutrino flux from the galaxy cluster population, assuming that the cumulative γ-ray flux from clusters is the dominant component of the IGRB.We note that the constraint from IGRB is based on the propagation model of lower-energy (a few TeV) CRs in the ICM.There may be additional physical effects on these lower-energy CRs and hence the obtained upper limits may not be valid.The constraints derived from both the upper limit on cumulative neutrino flux by IceCube analysis and the IGRB observed by Fermi-LAT are more robust than the one inferred from the γ-ray observations of the Coma cluster in Shi et al. (2022).Finally, we should bear in mind that the constraints on the baryon loading factor obtained here are valid for the entire source population over a long period of time comparable to the Hubble timescale.It is possible for an individual galaxy cluster or AGNs to exceed this limit, in particular during a short period of time such as during AGN flares.

Figure 4 .
Figure 4. Cumulative γ-ray flux from the galaxy clusters for different proton injection spectral index α = 2.0 (dashed), 2.2 (solid), 2.5 (dotted).The flux is compared with the IGRB observations from Fermi-LAT.The differential upper limits for the distance weighting (1/d 2 ) scheme from IceCube analysis are also shown in the figure for comparison.