Searching for High-energy Neutrino Emission from Galaxy Clusters with IceCube

The IceCube Neutrino Observatory has been observing a steady diffuse flux of TeV–PeV astrophysical neutrinos for over a decade (Aartsen et al. 2020a). Despite recent indications of potential point sources of neutrinos (Aartsen et al. 2018a, 2018b, 2020b), the origin of the majority of the flux remains undiscovered. The mystery is closely tied to the sources of ultrahigh-energy (>10¹⁸ eV) cosmic rays (UHECRs; Abu-Zayyad et al. 2013; Aab et al. 2017), as well as the GeV–TeV extragalactic gamma-ray background (EGB) observed by Fermi-LAT (Ackermann et al. 2015a), both of which have been observed to produce energy fluxes comparable to that of high-energy neutrinos above 100 TeV. This has led several authors to propose a unique astrophysical source class that could explain all three aforementioned fluxes (e.g., Murase et al. 2008; Tamborra et al. 2014; Zandanel et al. 2015; Fang & Murase 2018). Galaxy clusters are one such viable class of sources.

During large-scale structure formation, clusters of galaxies emerge as a result of accretion and mergers of smaller structures, producing megaparsec-scale shock waves (Miniati et al. 2000; Ryu et al. 2003) and a turbulent intracluster medium (ICM). The shocks can potentially accelerate cosmic rays (CRs) to energies as high as 10²⁰ eV (Inoue et al. 2007), depending on the mass of the cluster (Enßlin et al. 2011; Fang & Olinto 2016; Wiener & Zweibel 2019), through the first-order Fermi acceleration mechanism (Blandford & Eichler 1987). Other proposed candidates for CR acceleration within clusters are embedded compact objects such as active galactic nuclei (AGN; Stecker et al. 1991; Winter 2013; Fang & Murase 2018) and magnetars (Murase et al. 2009).

The predicted maximum accelerated CR energy and resultant neutrino flux are highly dependent on the properties of the embedded sources and vary greatly from model to model. For both the megaparsec-scale shock and embedded source scenarios, the accelerated CRs are confined in clusters for times exceeding the age of the universe (Berezinsky et al. 1997). These CRs interact with the ICM to produce pions, which decay, yielding neutrinos and gamma rays.

Galaxy clusters remain unobserved in very high energy gamma rays (>100 GeV). Fermi-LAT (Dutson et al. 2013; Huber et al. 2013) and imaging air Cherenkov telescopes (Abramowski et al. 2012; Arlen et al. 2012; Acciari et al. 2018) have reported upper limits on the gamma-ray emission from several nearby clusters. A detection of astrophysical neutrinos from galaxy clusters would constitute incontrovertible evidence for proton acceleration in these objects. Previous IceCube analyses have reported constraints on the neutrino emission from five nearby clusters (Abbasi et al. 2011a; Aartsen et al. 2013). However, the source class as a whole remains largely uncharted territory in terms of neutrino searches. Some works (Fang & Murase 2018; Hussain et al. 2021) predict that almost the entirety of the diffuse flux observed by IceCube could be explained by the cumulative emission from galaxy clusters. In this work, for the first time, we use cosmological data to test the contribution of galaxy clusters to the diffuse astrophysical neutrino flux. In particular, we focus on testing the scenarios in which the most massive clusters accelerate CRs in accretion shocks and/or cluster mergers (Fang & Olinto 2016; Hussain et al. 2021, 2022).

The mass and redshift range of galaxy clusters with the largest contribution to the neutrino flux varies across models. For this study, we choose a catalog of clusters that extends in redshift coverage between 0.01 and 0.97 and has cluster masses above 7.8 × 10¹³ M_⊙. This selection allows us to test neutrino production in some of the largest and latest-formed objects in the universe. Our catalog of galaxy clusters is taken directly from the Planck 2015 survey (Ade et al. 2016a). The catalog consists of all galaxy clusters detected through the variations observed in the brightness of the cosmic microwave background (CMB). These variations are caused by the inverse Compton scattering of CMB photons off high-energy electrons in the ICM, known as the Sunyaev–Zel'dovich (SZ) effect (Carlstrom et al. 2002). In this work, we solely use the 1094 sources for which a redshift estimate is available. The resulting selection of clusters covers 83% of the sky.

2.1. Catalog Completeness

In order to extrapolate the results from our source sample to the entire population of galaxy clusters, we need an estimate for the completeness of the Planck catalog. This requires a calculation of the fraction of clusters that would be reliably detected by the Planck mission given a theoretical distribution of galaxy clusters in the universe. We extract Planck's completeness as a function of redshift and cluster mass from the information given in Ade et al. (2016a) using simulations. We describe the steps briefly here and point the reader to Raghunathan et al. (2022) for more details.

We start by simulating maps in individual Planck frequency bands that contain signals from astrophysical foregrounds, the CMB, instrumental noise, and the thermal SZ (tSZ) signal. The astrophysical foregrounds include contributions from diffuse SZ signals (both kinematic and thermal) and also emission from dusty star-forming and radio galaxies. Galactic foregrounds are ignored, since the Planck galaxy cluster catalog is constructed after applying a mask to remove the majority of the galactic dust emission. The tSZ signal from a hypothetical cluster of a given mass and redshift is modeled using a generalized Navarro–Frenk–White profile (Nagai et al. 2007) calibrated with X-ray observations by Arnaud et al. (2010). The dimensionless pressure profile of the ICM is integrated along the line of sight to obtain the Compton y signal y(x), where x is a given position on the sky.

Next, we integrate y(x) over the angular extent of the cluster, R_500c , to obtain the integrated Comptonization parameter, Y_500c . Here the subscript 500c refers to a sphere where the mean mass overdensity of the cluster is at least 500 times the cosmic critical density at the redshift of the cluster. The Y_500c is directly related to the cluster mass (Ade et al. 2016b).

Following the maximum-likelihood approach described in Raghunathan et al. (2022), we combine the simulated maps from all Planck frequency bands and obtain the signal-to-noise ratio (S/N) for each cluster as a function of mass (M_500c ) and redshift (z) in the ranges $\mathrm{log}({M}_{500c}/{M}_{\odot })\in [13,16]$ in logarithmic steps of 0.1 and z ∈ [0.1, 3] in steps of 0.1. For detection, a cluster must have a minimum S/N of 5. Using the detected clusters, we estimate the catalog completeness as (Ade et al. 2016b; Alonso et al. 2016)

$$\begin{eqnarray}&&\chi ({Y}_{500c})=\displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{{Y}_{500c}^{\mathrm{true}}-{q}_{\mathrm{lim}}{\sigma }_{{Y}_{500c}}}{\sqrt{2}{\sigma }_{{Y}_{500c}}}\right)\right],\end{eqnarray} \tag{ 1 }$$

where ${q}_{\mathrm{lim}}=5$ is the detection threshold, ${\sigma }_{{Y}_{500c}}$ is the measurement uncertainty of the integrated Y_500c signal estimated above in simulations, and ${Y}_{500c}^{\mathrm{true}}$ is the true Y_500c flux. Figure 1 (right panel) shows the completeness obtained as a function of redshift and mass.

Zoom In Zoom Out Reset image size

We then create a theoretical distribution of galaxy clusters by drawing ${ \mathcal O }({10}^{5})$ samples from the Tinker 2010 halo mass function (HMF; Tinker et al. 2010) binned in redshift and mass (Figure 1, left panel). We only consider galaxy clusters with masses between 10¹⁴ and 10¹⁵ M_⊙ and redshifts between 0.01 and 2 in our simulation. In our energy range of interest and in models where CR acceleration takes place through shock waves due to the ICM accretion onto clusters, clusters with masses below 10¹⁴ M_⊙ or z > 1 are not expected to produce a significant flux of neutrinos at Earth (Fang & Olinto 2016). In the aforementioned scenarios, the bulk of the flux contribution comes from 0.3 < z < 1, in accordance with the cosmological distribution of sources (Fang & Olinto 2016).

We scale this theoretical distribution with the Planck completeness function obtained above to estimate the fraction of detectable clusters. This fraction is determined by dividing the sum of the scaled distribution by the total number of clusters in the HMF sample. In the case where each redshift and mass bin is equally weighted, we find a catalog completeness of ∼17%. Other weighting schemes are discussed in the following section. We note that our completeness calculations are robust within 2% under different choices of the HMF.

2.2. Source Weighting

The IceCube Neutrino Observatory consists of a cubic kilometer of South Pole ice with 5160 digital optical modules (DOMs) deployed between 1450 and 2450 m below the surface in Antarctica. IceCube is composed of 86 vertical strings, with 60 DOMs per string. The DOMs consist of a photomultiplier tube and onboard readout electronics designed to detect Cherenkov light induced by charged particles that can be produced in charged-current (CC) and neutral-current (NC) neutrino interactions (Abbasi et al. 2009, 2010). These secondary particles result in either track-like (CC interactions of muon neutrinos) or cascade-like (all other CC and NC interactions) optical signatures that are recorded by DOMs. The photon counts and timing information are used to reconstruct both the energy deposition and particle direction. IceCube can detect all flavors of neutrinos and cannot discriminate between neutrinos and antineutrinos. The angular resolution achieved depends on the optical signature, with tracks yielding the best resolution, typically below 1° above 1 TeV (Aartsen et al. 2017a).

Due to the superior angular resolution, the analysis presented in this work uses a 9.5 yr data sample consisting exclusively of muon tracks. The data set, recorded between 2008 April and 2017 November (Aartsen et al. 2017a, 2017b), covers the full sky. The majority of the data set consists of background events induced by atmospheric neutrinos and muons created by CR air showers in the upper atmosphere. A detailed description of the data set and the related reconstruction algorithms used in this work can be found in Aartsen et al. (2014) and Abbasi et al. (2011b).

4.1. Likelihood and Test Statistic

In order to search for the cumulative neutrino emission from our catalog of galaxy clusters in data dominated by atmospheric backgrounds, we perform a stacking analysis using an unbinned maximum-likelihood method (Achterberg et al. 2006; Braun et al. 2008).

The likelihood function for M stacked clusters is defined as

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }\left({n}_{s},\gamma \right)\\ \ \ =\prod _{i}^{N}\left(\displaystyle \frac{{n}_{s}}{N}\sum _{k}^{M}\displaystyle \frac{{r}_{k}{w}_{k}}{\sum _{j}^{M}{r}_{j}{w}_{j}}{S}_{i}^{k}\left(\gamma \right)+\left(1-\displaystyle \frac{{n}_{s}}{N}\right){B}_{i}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where n_s is the number of signal events, γ is the spectral index of an assumed unbroken power-law spectrum that is common to all sources, N is the total number of data events, S_i ^k is the signal probability distribution function (PDF) for the ith neutrino event for the kth source, B_i is the background PDF for the ith neutrino event, and r_k and w_k are stacking weights described below.

The background PDF is a function of both the reconstructed event energy, E_i , and the reconstructed event decl., δ_i , and is calculated by randomizing the R.A. of events in a given decl. band. The signal PDF is a weighted superposition of the signal PDFs of all sources in our selection of galaxy clusters. A potential source of astrophysical neutrinos should be distinguishable from the background due to the spatial clustering of high-energy events in its vicinity. Therefore, each S_i ^k consists of a Gaussian spatial term and an energy term that favors high-energy events,

$$\begin{eqnarray}&&{S}_{i}^{k}=\displaystyle \frac{1}{2\pi {\sigma }_{i}^{2}}\exp \left(-\displaystyle \frac{| {{\rm{\Psi }}}_{{ik}}{| }^{2}}{2{\sigma }_{i}^{2}}\right)\times { \mathcal E }({E}_{i},{\delta }_{i},\gamma ).\end{eqnarray} \tag{ 3 }$$

The term Ψ_ik is the angular separation between event i and source k, σ_i ⁶⁴ is the angular uncertainty of event i, and ${ \mathcal E }({E}_{i},{\delta }_{i},\gamma )$ is a simulated energy PDF, where γ is assumed to have a value between 1 and 4. Nearby clusters with extended profiles are not part of our sample; therefore, the likelihood function is independent of source extension.

The source weight w_k in Equation (2) parameterizes the relative contribution of each cluster to the total signal. It is determined by the weighting schemes described in Section 2.2. The model-independent term, r_k , depends on the detector acceptance and is determined by the effective area of the detector at the source location and the assumed neutrino spectrum. For more details, see Abbasi et al. (2022a) and Aartsen et al. (2017c).

In order to determine the number of signal events and the spectral index that best describes the data, the likelihood in Equation (2) is maximized with respect to n_s and γ and tested against the null hypothesis, which corresponds to the background-only scenario (n_s = 0). For that purpose, we define a test statistic (TS) as

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{log}\displaystyle \frac{{ \mathcal L }\left({\hat{n}}_{s},\hat{\gamma }\right)}{{ \mathcal L }\left({n}_{s}=0\right)},\end{eqnarray} \tag{ 4 }$$

where ${\hat{n}}_{s}$ and $\hat{\gamma }$ are the best-fit values. We note that each choice of weights will result in a different fitted ${\hat{n}}_{s}$, $\hat{\gamma }$, and TS.

4.2. Sensitivity

We perform a stacking analysis on our catalog as described in Section 4 for the different weighting schemes described in Section 2.2. None of the scenarios we investigate show a significant excess of signal events over background. The best-fit number of signal events, as well as the corresponding TS value, is zero for each weighting scheme.

In order not to report overly constraining limits due to potential underfluctuations in the data, we set the upper limits at the 90% sensitivity level described above. Table 1 shows the upper limits on the unresolved flux of galaxy clusters at 100 TeV for each weighting scenario after scaling for the effective completeness of the catalog and correcting for the 83% sky coverage of Planck.

Table 1. Upper Limits at 90% CL on the Diffuse Muon Neutrino Flux of Planck SZ Galaxy Clusters after Scaling for the Effective Catalog Completeness Corresponding to Each Weighting Scheme Considered

	γ = 2.0		γ = 2.5		γ = 3.0
Weighting Scheme	${E}_{0}^{2}\tfrac{{dN}}{{dE}}$	Diffuse Contribution	${E}_{0}^{2}\tfrac{{dN}}{{dE}}$	Diffuse Contribution	${E}_{0}^{2}\tfrac{{dN}}{{dE}}$	Diffuse Contribution
Equal	5.31 × 10−12	37%	1.11 × 10−11	77%	4.70 × 10−12	32.7%
Distance	6.53 × 10−13	4.5%	8.96 × 10−13	6.2%	2.41 × 10−13	1.7%
Mass	6.12 × 10−13	4.3%	6.58 × 10−13	4.6%	2.21 × 10−13	1.5%

Note. For each weighting scheme, upper limits are quoted for three E^−γ power-law spectra. The upper limits are presented in terms of the diffuse galaxy cluster flux at E₀ = 100 TeV, ${E}_{0}^{2}\tfrac{{dN}}{{dE}}$, as well as its relative contribution to the diffuse muon neutrino observations of Abbasi et al. (2022b) at 100 TeV. All fluxes are given in units of TeV cm⁻² s⁻¹ sr^−1.

Download table as:  ASCIITypeset image

The upper limits set in this analysis can then be compared with the diffuse muon neutrino flux measured by IceCube (Abbasi et al. 2022b) to estimate the maximum contribution of the galaxy cluster population to the diffuse observations. This fraction is also reported in Table 1 for each scenario and assumed spectral index. In addition, Figure 2 shows the upper limits on the neutrino emission from clusters for an assumed E^−2.5 power-law spectrum. This limit is directly compared to the diffuse muon neutrino observations over the full sensitive energy range of IceCube, which have a best-fit E^−2.37 power-law spectrum (Abbasi et al. 2022b).

Zoom In Zoom Out Reset image size

We find the equal-weighting case to provide the weakest constraints across all spectral indices. This is expected because this weighting scheme is meant to act as a check on our methods without corresponding to a realistic physical scenario. The distance and mass weighting schemes provide similar constraints to one another, which is also expected, as they are similar in terms of both the physics cases they are testing and their effective completeness. Both cases assume that the neutrino flux falls off with the cluster distance as d⁻². Distance weighting assumes that galaxy clusters are neutrino standard candles. On the other hand, the mass weighting assumes linear proportionality between neutrino luminosity and cluster mass. For a spectral index of 2.5, the results exclude a maximum contribution of 77% at 100 TeV from galaxy clusters assuming equal weights. Under the more realistic distance and mass weighting schemes, this contribution is constrained down to 6.2% and 4.6% at 100 TeV, respectively.

Another useful way of comparing the contribution of galaxy clusters to diffuse neutrino observations is by evaluating the upper limits in quasi-differential energy bins. We calculate the upper limits in one-decade energy bins, assuming an E⁻² spectrum across the bin. Figure 3 shows the quasi-differential limits for the distance weighting scheme along with comparable theoretical predictions from the mass and redshift range used in this work. The IceCube upper limits obtained in this work are most constraining between 10 TeV and 1 PeV compared to the diffuse neutrino observations, where we find that the population of galaxy clusters in our mass and redshift range of interest cannot contribute to more than 9%–13% of the diffuse neutrino flux.

Zoom In Zoom Out Reset image size

Our results above 100 TeV constrain the model based on magnetohydrodynamic (MHD) simulations by Hussain et al. (2021), which predict that most of the contribution to the diffuse neutrino flux can be attributed to clusters with masses above 10¹⁴ M_⊙. For the aforementioned model, our limits exclude a CR injection spectrum of 1.5 and a conservative energy conversion rate of 0.5%. Our limits are robust between 100 TeV and 1 PeV to the choice of the injection spectrum and maximum CR energy. These conclusions can change for a different choice of kinetic energy conversion rates and assumptions about source evolution. We note that certain models of structure formation followed by hadronic interactions in galaxy clusters suggest that the total flux from clusters is indeed very low and not able to explain more than a few percent of the total diffuse flux (Zandanel et al. 2015; Fang & Olinto 2016). Our analysis is not sensitive enough to probe such models. As explained in Section 2.2, models that assume AGN as the embedded sources of CRs in clusters (Fang & Murase 2018) are not constrained using the mass range of clusters tested in this work.

In this work, we presented a search for neutrinos produced by a large catalog of galaxy clusters, which are expected to produce neutrinos through embedded sources, accretion shocks, or both. We performed a stacking analysis with 1094 SZ-selected galaxy clusters from Planck using 9.5 yr of muon-track IceCube data. Three different weighting schemes were investigated, i.e., equal weights, distance weights (1/d²), and mass weights (M/d²). We find no evidence for neutrino emission in any of the scenarios considered. Assuming an unbroken E^−γ power-law flux common among all galaxy clusters, we compute upper limits at 90% CL on the diffuse neutrino flux from galaxy clusters with masses between 10¹⁴ and 10¹⁵ M_⊙ and a redshift between 0.01 and 2 using a (simulated and extrapolated in z) complete catalog of galaxy clusters. We constrain the galaxy cluster contribution to the diffuse muon neutrino observations of Abbasi et al. (2022b) at 100 TeV. In particular, for γ = 2.5, galaxy clusters can contribute no more than 67%, 6.2%, and 4.6% for the equal, distance, and mass weighting scenarios, respectively.

We also compute quasi-differential limits at 90% CL on the diffuse neutrino flux from galaxy clusters. These limits exclude the MHD model of Hussain et al. (2021) as neutrino production mechanisms in galaxy clusters. Accretion shock model predictions by Zandanel et al. (2015) and Fang & Olinto (2016) are not constrained by our limits, making our results consistent with scenarios in which accretion shocks make a subdominant contribution to the diffuse neutrino flux (Murase & Waxman 2016).

An absence of GeV–TeV gamma rays from nearby clusters, along with our new IceCube limits, constrains the shock acceleration efficiencies in galaxy clusters. Galaxy clusters have long been proposed as the single unified solution to explain the diffuse fluxes of UHECRs, neutrinos, and gamma rays (Fang & Murase 2018). Our results encourage a revisiting of the multimessenger connections between different classes of these high-energy particles. The next generation of observatories, such as IceCube Gen-2 (Aartsen et al. 2021), CTA (Pérez-Romero 2022), and AMEGO-X (Fleischhack 2021), offers a promising prospect for the detection of smoking-gun MeV–PeV multimessenger signatures of CR acceleration in galaxy clusters.

The IceCube Collaboration acknowledges the significant contributions made to this manuscript by Mehr Un Nisa, Andrew Ludwig, and Srinivasan Raghunathan. We also acknowledge support from: USA - US National Science Foundation Office of Polar Programs, US National Science Foundation Physics Division, US National Science Foundation EPSCoR, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin–Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), the Frontera computing project at the Texas Advanced Computing Center, US Department of Energy National Energy Research Scientific Computing Center, the particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and the astroparticle physics computational facility at Marquette University; Belgium - Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany - Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden - Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia - Australian Research Council; Canada - Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark - Villum Fonden and Carlsberg Foundation; New Zealand - Marsden Fund; Japan - Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea - National Research Foundation of Korea (NRF); Switzerland - Swiss National Science Foundation (SNSF); United Kingdom - Department of Physics, University of Oxford.