Dark-Matter-Enhanced Probe of Relic Neutrino Clustering

TL;DR

This paper explores how decaying heavy dark matter can produce ultrahigh-energy neutrinos, which, through their interactions with the cosmic neutrino background, can reveal properties of relic neutrino clustering with future neutrino telescopes.

Contribution

It introduces a novel method using heavy decaying dark matter as a probe for relic neutrino background and clustering effects via ultrahigh-energy neutrino observations.

Findings

Local CνB overdensity of >10^6 can be detected at IceCube-Gen2 within ten years.

Neutrino interactions with CνB can provide insights into neutrino mass and dark matter properties.

Complete neutrino transport modeling enhances the accuracy of flux predictions.

Abstract

We propose heavy decaying dark matter (DM) as a new probe of the cosmic neutrino background (C$\nu$B). Heavy DM, with mass $\gtrsim 10^9$ GeV, decaying into neutrinos can be a new source of ultrahigh-energy (UHE) neutrinos. Including this contribution along with the measured astrophysical and predicted cosmogenic neutrino fluxes, we study the scattering of UHE neutrinos with the C$\nu$B via standard weak interactions mediated by the $Z$ boson. We solve the complete neutrino transport equation, taking into account both absorption and reinjection effects, to calculate the expected spectrum of UHE neutrino flux at future neutrino telescopes, such as the IceCube-Gen2 radio. We argue that such observations can be used to probe the C$\nu$B properties and, in particular, local C$\nu$B clustering. We find that, depending on the absolute neutrino mass and the DM mass and lifetime, a local C$\nu$B overdensity $\gtrsim 10^6$ can be probed at the IceCube-Gen2 radio within ten years of data taking.

Equations

$$n_{\nu ,0}=\frac{3}{4}\frac{\zeta (3)}{{\pi }^2}{T}_{\nu ,0}^3\simeq 56{\mathrm{cm}}^{-3},$$

$$\frac{d{\Phi }_{{\nu }_{\alpha }+{\bar{\nu }}_{\alpha }}^{\text{DM}}}{dE}=\frac{{\Omega }_{\mathrm{DM}}{\rho }_c}{4\pi m_{\mathrm{DM}}{\tau }_{\mathrm{DM}}}{\int }_0^{\infty }\frac{dz}{H(z)}\frac{d{N}_{\alpha }}{dE}(E,z).$$

$$N_{\text{w}(\text{w/o})}=T{\int }_{E_{\nu }^{\min }}^{E_{\nu }^{\max }}dE\int d\Omega A_{\mathrm{eff}}(E,\Omega ){\left(\frac{d{\Phi }_{{\nu }_{\alpha }+{\bar{\nu }}_{\alpha }}}{dE}\right)}_{\text{w}(\text{w/o})},$$

$${\chi }^2=2\left(N_{\mathrm{w}}-N_{w/o}+N_{w/o}\log \frac{N_{w/o}}{{\mathrm{N}}_{\mathrm{w}}}\right).$$

$${\Omega }_{\mathrm{DM}}{h}^2\sim {\left(\frac{m_{\mathrm{DM}}}{10^{11}\mathrm{GeV}}\right)}^2\left(\frac{T_{\mathrm{RH}}}{10^9\mathrm{GeV}}\right).$$

$${\Omega }_{\mathrm{DM}}{h}^2\simeq$$

$$\times \left(\frac{\mathrm{Br}}{10^{-8}}\right).$$

$$\frac{\partial {\tilde{N}}_A(E,t)}{\partial t}$$

$$+{\int }_E^{\infty }d{E}^{\prime }{\tilde{N}}_A(E^{\prime },t){\mathcal{B}}_{A\to A}(E,E^{\prime })$$

$$+\sum _{B\ne A}{\int }_E^{\infty }d{E}^{\prime }{\tilde{N}}_B(E^{\prime },t){\mathcal{C}}_{B\to A}(E,E^{\prime })$$

$$+{\tilde{\mathcal{Q}}}_A(E,t),$$

Full text

††preprint: FERMILAB-PUB-25-0566-T, MITP-2025-045, CETUP-2025-004

I Introduction

Analogous to the cosmic microwave background (CMB), the standard cosmological model predicts the existence of a cosmic neutrino background (C$\nu$B) – a relic population of neutrinos that decoupled from the primordial plasma approximately one second after the big bang, at temperatures on the order of 1 MeV Weinberg (1962); Scott (2024). After decoupling, neutrinos free stream through the Universe, retaining a Fermi-Dirac distribution, with their temperature evolution following an adiabatic expansion: $T_{\nu}\propto a^{-1}$, where $a(t)$ is the scale factor. After electron-positron annihilation, which heats the photon bath but not the neutrinos, the neutrino temperature becomes slightly lower than that of the photons. Consequently, the present-day C$\nu$B temperature can be written in terms of the CMB temperature as $T_{\nu,0}=(4/11)^{1/3}T_{\gamma,0}\simeq 1.95~{\rm K}\simeq 1.68\times 10^{-4}~{\rm eV}$. Assuming that the lightest neutrino mass $m_{1}\gtrsim T_{\nu,0}$, the C$\nu$B today can be thought of as a nonrelativistic gas of fermions with a number density per flavor and per degree of freedom

[TABLE]

and a total number density of $6\,n_{\nu,0}$ (valid irrespective of whether neutrinos are Dirac or Majorana Long et al. (2014); Bauer and Shergold (2023)), considering three neutrino flavors and both neutrinos and antineutrinos.

The presence of C$\nu$B has only been indirectly confirmed through its imprints on (i) the big bang nucleosynthesis via the expansion rate and neutron-to-proton ratio Steigman (2012), (ii) the anisotropies in the CMB power spectrum Aghanim et al. (2020), (iii) the formation and distribution of large-scale cosmological structures Elbers et al. (2025), and (iv) acoustic oscillation phase shifts Follin et al. (2015). Despite several proposed strategies Stodolsky (1975); Weinberg (1962); Irvine and Humphreys (1983); Gelmini (2005); Yanagisawa (2014); Baracchini et al. (2018); Apponi et al. (2022); Long et al. (2014); Bauer and Shergold (2023); del Castillo et al. (2025), the direct detection of C$\nu$B remains elusive. The extremely low kinetic energy of C$\nu$B and the weakly interacting nature of neutrinos are the biggest obstacles to observing these relic neutrinos directly in the terrestrial laboratories via their scattering.

Therefore, an interesting alternative is to flip the script and consider the impact of the C$\nu$B as scatterers instead of scattering particles. The first such proposal is the so-called $Z$ burst Weiler (1982) (see also Refs. Roulet (1993); Yoshida et al. (1997); Fargion et al. (1999); Weiler (1999); Yoshida et al. (1998); Eberle et al. (2004); Barenboim et al. (2005)), where ultrahigh-energy (UHE) neutrinos scatter off C$\nu$B through a resonant process in which an on-shell $Z$ boson is exchanged. Such a process could, in principle, leave a distinct absorption feature in the cosmogenic neutrino flux measured at Earth. However, the neutrino energy required to produce an on-shell $Z$ boson, $E_{\nu}=m_{Z}^{2}/2m_{\nu}(1+z)\gtrsim 10^{13}$ GeV ($z$ being the redshift at which the $\nu\bar{\nu}$ annihilation occurs in the sky), far exceeds the Greisen–Zatsepin–Kuzmin (GZK) cutoff Greisen (1966); Zatsepin and Kuzmin (1966), which sets a theoretical upper bound on the energy of cosmic rays that are the progenitors of the cosmogenic neutrinos Berezinsky and Zatsepin (1969). The resonance energy gets smaller for the on-shell production of less massive particles compared to the $Z$ boson. This was employed in Ref. Brdar et al. (2022) for the case of $\rho$ meson exchange within the Standard Model (SM), where the resonance occurs at $E_{\nu}=m_{\rho}^{2}/2m_{\nu}(1+z)\gtrsim 10^{9}$ GeV. By numerical coincidence, the position of the $\rho$ resonance overlaps with the peak of the predicted cosmogenic neutrino flux, which offers an enhanced absorption effect at higher redshifts. In presence of new (light) scalars or gauge bosons, similar resonances might occur, even at lower energies Ng and Beacom (2014); Ibe and Kaneta (2014); Araki et al. (2015); Kamada and Yu (2015); DiFranzo and Hooper (2015); Altmannshofer et al. (2016). However, there are two major challenges in searching for resonance-induced absorption features in the cosmogenic neutrino spectrum: (i) The width of the resonance must be comparable to or greater than the resolution of the experiment in order to resolve the dip feature. This happens to be the case for the $\rho$ meson which is the lightest vector meson in the SM and has a width-to-mass ratio of 19% Navas et al. (2024). In contrast, almost all other SM resonances have narrow widths Navas et al. (2024); Dev and Soni (2022), thus making the $\rho$ meson possibly the only SM resonance in the $\nu\bar{\nu}$ channel relevant at future neutrino telescopes. (ii) There is a large uncertainty on the cosmogenic neutrino flux predictions Berezinsky and Zatsepin (1969); Engel et al. (2001); Hooper et al. (2005); Allard et al. (2006); Ahlers et al. (2010); Kotera et al. (2010); Kampert and Unger (2012); Ahlers and Halzen (2012); Aloisio et al. (2015); Alves Batista et al. (2019); Møller et al. (2019); van Vliet et al. (2019); Heinze et al. (2019); Muzio et al. (2019, 2022); Valera et al. (2023); Ehlert et al. (2024); Muzio et al. (2023), which is mainly due to the current uncertainty in the cosmic ray composition at the highest energies and the unknown cosmic ray source evolution history. The recent measurements from Auger Aab et al. (2019) and IceCube Abbasi et al. (2025a) have already excluded some of the cosmogenic flux models mentioned above; however, there still remains almost an order of magnitude uncertainty in the flux prediction among the allowed models, thus making it further difficult to search for a dip feature.

In light of these difficulties, we complement Ref. Brdar et al. (2022) in the present work in two aspects. First, instead of resonant production, we focus on off-shell $Z$ boson-induced scattering of UHE neutrinos with C$\nu$B Roulet (1993) and hence consider a wide energy range rather than resonance energies only. Second, in addition to astrophysical and cosmogenic neutrinos, we add the neutrino flux contribution from the decay of superheavy dark matter (SHDM) particles Gelmini and Kusenko (2000). Specifically, we consider decaying neutrinophilic dark matter (DM) particles with mass $m_{\text{DM}}\gtrsim 10^{9}$ GeV. Such an alternative DM scenario is well motivated Chung et al. (1998a), especially in the era of the waning of the weakly interacting massive particle (WIMP) paradigm Arcadi et al. (2025). SHDM easily evades the direct detection and collider constraints, unless they are strongly interacting Albuquerque and Baudis (2003). Moreover, the standard indirect detection signals from DM annihilation are highly suppressed, because of the small annihilation cross section $\sigma v$ and small DM number density $n_{\rm DM}$. The indirect detection prospects for SHDM improve if the DM is not absolutely stable as the decay signal scales like $n_{\rm DM}$ (instead of $n_{\rm DM}^{2}$ as in annihilation). This is conceivable because discrete symmetries typically invoked to stabilize DM are expected to be broken at high-energy scales by quantum gravity effects. Note that the SHDM production in the early Universe usually requires a nonthermal mechanism to obtain the correct relic abundance, since the thermal freeze-out due to DM annihilation in a standard cosmological history results in DM overproduction for masses larger than $\sim 100$ TeV–the so-called unitarity bound Griest and Kamionkowski (1990).111For possible ways to relax or evade the unitarity bound, see, e.g., Ref. Nussinov (2025). Explicit constructions for SHDM production and decay have been considered in inflationary models Chung et al. (1998b, 2001); Kannike et al. (2017); Kolb and Long (2017); Li et al. (2019); Haro and Aresté Saló (2019); Babichev et al. (2020); Ling and Long (2021); Verner (2025); Murase et al. (2025), in string theory Benakli et al. (1999); Allahverdi et al. (2021, 2024); Ling et al. (2025), and in seesaw framework Uehara (2001); Huong and Dong (2017); Dev et al. (2016); Sui and Dev (2018); Deligny (2024). Here we take a model-independent, phenomenological approach and simply treat the SHDM mass $m_{\rm DM}$ and lifetime $\tau_{\rm DM}$ as free parameters. We will choose representative benchmark points such that these quantities are consistent with the existing limits, as summarized in Appendix A. We assume that the correct DM relic abundance can always be obtained for these benchmarks, depending on the details of the explicit model construction and cosmological history; see Appendix B for a short discussion and parametric estimate indicating the relevant cosmological scales.

Including the two above-mentioned effects, namely, (i) the off-shell $Z$ contributions and (ii) the DM-induced neutrino flux, we find that future neutrino telescopes like IceCube-Gen 2 radio Aartsen et al. (2021) can probe the C$\nu$B, but only in the presence of a local overdensity. We will denote by $\xi$ the dimensionless ratio of the testable C$\nu$B number density to the standard cosmological value given in Eq. (1).

As neutrinos are massive, the C$\nu$B can become clustered because of the gravitational potential felt from both baryonic matter and DM. Hence, the local neutrino density is expected to be enhanced; however, such an enhancement can at most be at the level of few tens of percent in the Milky Way Ringwald and Wong (2004); de Salas et al. (2017); Mertsch et al. (2020); Hotinli et al. (2023); Zimmer et al. (2023); Holm et al. (2024). Nevertheless, very large overdensities of C$\nu$B can be achieved in the presence of new physics. For example, attractive Yukawa interactions of neutrinos with a new light scalar boson can lead to the formation of stable bound states and neutrino clusters Wise and Zhang (2015); Smirnov and Xu (2022) with maximal central density up to $10^{9}~{\rm cm}^{-3}$, which translates to an overdensity parameter $\xi\sim 10^{7}$.

The current most stringent local constraint on $\xi$ comes from the Karlsruhe tritium neutrino (KATRIN) experiment: $\xi<1.1\times 10^{11}$ at 95% C.L. Aker et al. (2022). The high-precision tracking of astronomical objects in Solar System places comparable constraints of the same order Tsai et al. (2024). These two probes are sensitive to overdensities located closer than 10 A.U. from Earth. For other astrophysical probes including overdensities across larger distances, see Refs. Bauer et al. (2021); Brdar et al. (2022); Císcar-Monsalvatje et al. (2024); Franklin et al. (2025); De Marchi et al. (2025); Herrera et al. (2025); Dev et al. (2025); Chauhan (2025); Das et al. (2025a); Zhang et al. (2025), where the derived constraints vary between $\xi\sim 10^{6}$ and $\xi\sim 10^{14}$. The strongest astrophysical bound comes from relic neutrino upscattering with ultrahigh-energy cosmic rays (UHECRs) Herrera et al. (2025); Zhang et al. (2025), but this crucially depends on the UHECR composition and source evolution, both of which are largely uncertain at the moment. In this work, by comparing expected neutrino event rates induced by heavy DM decay with and without overdensity at the future IceCube-Gen2 radio experiment Aartsen et al. (2021); Abbasi et al. (2021); Glaser (2025), we derive a projected sensitivity of $\xi\gtrsim 10^{6}$, assuming a single overdense C$\nu$B cloud centered at redshift $z=0$. The radius of the cloud $R$ is calculated assuming that the total number of relic neutrinos matches the one in standard cosmology and is given by $R=\frac{c}{H_{0}}\xi^{-1/3}$, where $H_{0}$ is the Hubble constant today and $c$ is the speed of light in vacuum. Our projected sensitivity is comparable to those from other astrophysical probes and is significantly better than the current laboratory constraint.

The rest of this paper is organized as follows. In Section II, we discuss the production of UHE neutrinos from the decay of heavy DM particles and from astrophysical environments, as well as their evolution to Earth. In Section III, we use these UHE neutrino fluxes to obtain IceCube-Gen2 radio event rates and we derive the corresponding sensitivity to a C$\nu$B overdensity. We conclude in Section IV. In Appendix A, we summarize the existing constraints on SHDM mass and lifetime. In Appendix C, we present the neutrino transport equation.

II Ultrahigh-Energy Neutrino Flux

We consider a heavy scalar DM, with masses above $10^{9}$ GeV, that can decay into neutrino-antineutrino pairs of all three flavors ($\text{DM}\rightarrow\nu_{\alpha}\bar{\nu}_{\alpha}$). The differential flux of neutrinos and antineutrinos produced in this way has two components: galactic and extragalactic. Extragalactic neutrinos encounter a significantly larger fraction of the C$\nu$B cloud on their way to Earth, and hence we focus exclusively on this contribution. Specifically, including the galactic component would worsen the sensitivity to C$\nu$B, simply because this component of the flux is hardly impacted through scattering with C$\nu$B given the much shorter propagation distances involved. Since radio experiments such as IceCube-Gen2 radio are expected to have exceptional angular resolution below a few degrees Hallgren (2021); Bouma and Nelles (2023), neutrinos coming from the direction of the Galactic Center can be isolated, as recently demonstrated by IceCube with the detection of neutrinos from the Galactic Plane Abbasi et al. (2023).

The differential flux of extragalactic neutrinos and antineutrinos of flavor $\alpha$ reads Esmaili et al. (2012); Brdar et al. (2017)

[TABLE]

Here, $\rho_{c}$ is the critical energy density of the Universe, $\Omega_{\mathrm{DM}}$ is the relic abundance of DM, $z$ is the cosmological redshift and $H(z)=H_{0}\sqrt{\Omega_{\Lambda}+\Omega_{\rm m}(1+z)^{3}}$ is the Hubble expansion rate with the dark energy and matter density components $\Omega_{\Lambda}=0.685$, $\Omega_{\rm m}=0.315$, and the Hubble constant $H_{0}=100\,h$ km.s*-1*.Mpc*-1* with $h=0.673$ Aghanim et al. (2020). The quantity $dN_{\alpha}/dE$ in the absence of neutrino interactions with C$\nu$B would simply be equal to the prompt energy spectrum $dN_{\alpha}^{\rm P}/dE$ (appropriately redshifted to account for cosmic expansion) that we calculate using HDMSpectra Bauer et al. (2021) and shown in Fig. 1.

However, UHE neutrinos produced from DM decays can scatter while propagating in the C$\nu$B, and this is what is captured by the neutrino transport equation Ng and Beacom (2014); Bustamante et al. (2020); Creque-Sarbinowski et al. (2021); Esteban et al. (2021); Das et al. (2025b) (see Appendix C). Specifically, this equation accounts for energy redshift due to the expansion of the Universe, as well as neutrino absorption and reinjection effects from scattering on the C$\nu$B. We adopt the cross sections for $\nu$-$\nu$, $\bar{\nu}$-$\bar{\nu}$, and $\nu$-$\bar{\nu}$ interactions through $Z$ boson exchange from Ref. Roulet (1993) (see also Ref. Dev and Soni (2022)). For brevity, in what follows we will just refer to all these processes simply by $\nu$-$\nu$. While we have investigated all the aforementioned terms in the transport equation, to keep the numerical treatment efficient and the runtime manageable, we eventually solved the transport equation by ignoring the reinjection terms, which we found not to significantly impact our final results. This is consistent with what was shown in Ref. Das et al. (2025b). With this procedure, we calculate the neutrino flux at Earth from extragalactic DM decays.

We also consider UHE neutrinos produced in the astrophysical environments through baryonic interactions. This additional flux consists of two types–cosmogenic and astrophysical neutrinos. The cosmogenic neutrinos are created when the cosmic ray protons or nuclei, during their propagation, interact with the CMB Berezinsky and Zatsepin (1969); Engel et al. (2001); Hooper et al. (2005); Allard et al. (2006); Ahlers et al. (2010); Kotera et al. (2010); Kampert and Unger (2012); Ahlers and Halzen (2012); Aloisio et al. (2015); Alves Batista et al. (2019); Møller et al. (2019); van Vliet et al. (2019); Heinze et al. (2019); Muzio et al. (2019, 2022); Valera et al. (2023); Ehlert et al. (2024); Muzio et al. (2023). The astrophysical neutrinos, on the other hand, are produced when the cosmic ray protons or heavier nuclei that are generated at the cosmic reservoirs or accelerators undergo interactions with each other or with the ambient photons, and produce neutrinos as by-products Meszaros (2017). Similarly as for DM flux, we employ the transport equation in order to propagate cosmogenic and astrophysical neutrinos through the C$\nu$B cloud.

In Fig. 2, we show both DM-induced and astrophysical plus cosmogenic neutrino fluxes at Earth, comparing scenarios with and without scattering off an overdense C$\nu$B cloud. We take $\xi=10^{11}$ (which saturates the current KATRIN bound Aker et al. (2022)) for illustration and assume normal neutrino mass ordering with the lightest neutrino mass $m_{1}=0.01$ eV. For the DM-induced flux component, we show two benchmark points with $m_{\text{DM}}=10^{11}$ GeV, one that saturates the lifetime constraint from $\gamma$ rays (green) and another one that saturates the lifetime constraint from neutrinos (orange) –see Fig. 5. The cosmogenic neutrino flux with the associated uncertainties (as shown by the gray hatched band) is taken from Ref. Muzio et al. (2023), derived by fitting the observed cosmic ray spectrum and composition data from Pierre Auger Aab et al. (2020). Here, the upper and lower lines correspond to the active galactic nucleus (AGN) and star formation rate (SFR) source evolution models, respectively. For the astrophysical flux component, we use a recent IceCube best-fit spectrum using a broken power law  Abbasi et al. (2025b). The light yellow band shows the astrophysical plus cosmogenic flux after scattering with the C$\nu$B cloud. Clearly, the scattering in regions with C$\nu$B overdensity alters the predicted spectrum, and this is precisely what our analysis is based on. In Fig. 2, we also show the flux upper limits from Auger Aab et al. (2019) and IceCube Abbasi et al. (2025a) (shaded regions in the upper half-plane), as well as the sensitivities of a number of future neutrino telescopes (taken from Ref. Ackermann et al. (2022)), including IceCube-Gen2 radio Aartsen et al. (2021); Abbasi et al. (2021) which we will use as a representative example in our analysis.

III Event Rates and Sensitivity to C$\nu$B Overdensity

The number of predicted neutrino events in IceCube-Gen2 radio experiment can be written as

[TABLE]

where w (w/o) in the subscript denotes with (without) including scattering off the overdense neutrino cloud. Further, in Eq. 3, $T=10\;\mathrm{yr}$ is the assumed data taking time and $A_{\mathrm{eff}}(E,\Omega)$ is the effective area of the IceCube-Gen2 radio detector which is a function of the neutrino energy and solid angle. Given the available solid-angle average effective area of IceCube-Gen2, we have $\int d\Omega A_{\textrm{eff}}(E,\Omega)=4\pi A_{\textrm{eff}}(E)$. The energy integration limits in Eq. 3 are taken to be $E_{\nu}^{\mathrm{min}}=10^{7.5}\;$ and $E_{\nu}^{\mathrm{max}}=10^{15}\;$GeV. The IceCube-Gen2 radio effective area is available for the energy range between $10^{7.5}\;$GeV and $10^{11}\;$GeV Abbasi et al. (2021); above $10^{11}\;$GeV, we conservatively assume the effective area to be constant, although a $\log(E)$ growth is expected.

We calculate the projected sensitivity on the local C$\nu$B overdensity, $\xi$, assuming that the likelihood ratio test statistic follows the $\chi$-squared distribution with 1 degree of freedom,

[TABLE]

For the IceCube-Gen2 radio, to detect the signature of $\nu$-$\nu$ scattering at 90% C.L., the value of $\chi^{2}$ needs to exceed $2.7$. By utilizing this condition, we calculate the projected sensitivity on $\xi$.

Figure 3 (left panel) shows our projected IceCube-Gen2 radio sensitivity to $\xi$ in the case where only neutrino fluxes from DM decay are considered. The sensitivity is shown as a function of the lightest neutrino mass $m_{1}$ for two values of DM mass: $m_{\mathrm{DM}}=10^{9}\;\mathrm{GeV}$ (red) and $m_{\mathrm{DM}}=10^{15}\;\mathrm{GeV}$ (orange). For the lighter of these, we show two lines corresponding to different DM-induced neutrino fluxes. The solid (dashed) line represents the case where neutrino fluxes are computed such that the $\gamma$-ray (neutrino) DM lifetime constraints in Fig. 5 are saturated. For $m_{\text{DM}}=10^{9}$ GeV these two DM lifetime constraints are comparable; therefore, the expected sensitivity to $\xi$ is of the same order. The derived sensitivity limits on $\xi$ for $m_{\text{DM}}=10^{9}$ GeV are weaker than the KATRIN overdensity bound Aker et al. (2022) (purple horizontal line), except for $\mathcal{O}$(eV) values of the lightest neutrino mass; that region, however, is in tension with KATRIN’s direct mass limit Aker et al. (2025) (brown vertical line). For $m_{\text{DM}}=10^{15}$ GeV, sensitivity down to $\xi\gtrsim 10^{6}$ is reachable for the case where the DM-induced neutrino flux is produced such that bounds on the DM lifetime from neutrino experiments are saturated. For this value of DM mass, we do not show the scenario with DM-induced neutrino flux that is in accord with DM lifetime limits from $\gamma$-ray experiments, since such neutrino flux would be too small to yield any events at the IceCube Gen-2 radio experiment. For comparison, we also show the cosmological upper limit on the lightest neutrino mass derived from the Planck measurements Aghanim et al. (2020). It should, however, be noted that the exact value of this limit depends on the cosmological datasets used Elbers et al. (2025); Chebat et al. (2025). Moreover, it is strictly applicable at high redshifts, and there are ways to relax it at lower redshifts, e.g. by assuming a nonstandard cosmology Bellomo et al. (2017); Esteban and Salvado (2021); Esteban et al. (2022); Craig et al. (2024), in presence of nonstandard neutrino interactions Beacom et al. (2004); Hannestad (2005); Farzan and Hannestad (2016); Escudero et al. (2023); Benso et al. (2025); Das et al. (2025c), decaying neutrinos Chacko et al. (2020); Escudero et al. (2020); Franco Abellán et al. (2022), time-varying neutrino mass Fardon et al. (2004); Dvali and Funcke (2016); Lorenz et al. (2019, 2021); de Gouvêa et al. (2022); Sen and Smirnov (2025), or neutrinos with a modified distribution function Oldengott et al. (2019); Alvey et al. (2022). Therefore, we do not shade the region disfavored by standard cosmology.

The right panel of Fig. 3 exhibits the computed projected upper limit on $\xi$ based on IceCube-Gen2 radio sensitivity as a function of the DM mass $m_{\text{DM}}$ for four representative values of the lightest neutrino mass: $0.01$ (red), $0.05$ (blue), $0.1$ (orange), and $0.45$ eV (brown). For each value of the DM mass, the DM lifetime, which controls the magnitude of the induced neutrino fluxes, is taken such that the DM lifetime constraints from neutrino searches are saturated. In this case, we are able to probe overdensities down to $\xi\gtrsim 10^{6}$.

Note that generating DM-induced neutrino fluxes saturating the neutrino constraint assumes significant relaxation of the $\gamma$-ray limits. The weakening of these limits by an order of magnitude is possible due to uncertainties in the galactic magnetic field modeling Munbodh and Profumo (2024). Another promising possibility is discussed in Ref. Hiroshima et al. (2018), where DM decays into invisible hidden states in addition to SM neutrinos; this results in relaxed $\gamma$-ray limits as well as in the broadening of the neutrino spectrum.

Now, we consider the scenario in which all available neutrino fluxes–from both astrophysical sources and DM decay–are included in the calculation of the sensitivity limits. This is shown in Fig. 4, where in the left panel we present the expected sensitivity for $\xi$ as a function of the lightest neutrino mass $m_{1}$. The DM-induced part of the total flux is computed for two DM masses ($10^{9}$ and $10^{15}$ GeV), where for each case the flux is calculated by either saturating neutrino or $\gamma$-ray constraints on the DM lifetime. We observe that both benchmark points computed with $\gamma$-ray constraints on the DM lifetime, as well as the benchmark point $m_{\text{DM}}=10^{9}$ GeV ($\nu$ constraint), yield similar results, represented by slightly shifted red, green, and blue bands. This signifies that the sensitivity is only mildly impacted by DM, being dominated by cosmogenic and astrophysical neutrinos. The projected sensitivity for C$\nu$B overdensity in this regime varies between $\xi\sim 10^{8}$ and $\xi\sim 10^{10}$, depending on the value of the lightest neutrino mass. The band illustrates the uncertainty in the theoretical modeling of the cosmogenic fluxes (see, e.g., Ref. Muzio et al. (2023)). In the left panel of  Fig. 4, we also observe that a significant relaxation of the $\gamma$-ray constraints could further improve these predictions, as the DM-induced neutrino flux would not only become competitive with but even surpass the astrophysical and cosmogenic flux. This occurs for large DM masses, as represented by the orange band calculated using a DM mass of $10^{15}$ GeV and a lifetime taken according to neutrino constraints. For this case, we find a projected sensitivity on $\xi$ reaching down to approximately $10^{6}$, thereby converging to the respective DM-only case shown in the left panel of Fig. 3.

This effect can also be seen in the right panel of Fig. 4, where for $m_{\mathrm{DM}}\gtrsim 10^{11}$ GeV the sensitivity limits on the local C$\nu$B overdensity calculated using the saturated neutrino DM lifetime constraints (densely hatched region) start to dominate over the limits calculated using the saturated $\gamma$-ray DM lifetime constraints (sparsely hatched region). Note that this effect is present irrespective of the value of the lightest neutrino mass. In this regime, where the DM-induced neutrino flux dominates over the cosmogenic flux, there is less impact of the astrophysical uncertainty due to UHECR modeling on the overdensity sensitivity.

Once the UHE neutrino flux is observed, distinguishing whether it corresponds to an attenuated spectrum due to a local C$\nu$B overdensity or an unattenuated flux with a different normalization is crucial. We outline several ways to approach this problem. First, the UHE neutrino fluxes in these two cases would exhibit different neutrino energy dependencies, which future precise spectral measurements could disentangle Aartsen et al. (2021). Second, independent $\gamma$-ray observations could reveal the presence of a DM component, thereby favoring an overdensity interpretation. Third, future results from direct probes of local neutrino overdensities in experiments such as KATRIN Aker et al. (2022), Project-8 Asner et al. (2015), HOLMES Alpert et al. (2025), or PTOLEMY Baracchini et al. (2018) may provide complementary evidence. Last but not least, improved modeling of the cosmogenic neutrino flux, driven by a more precise determination of the UHE cosmic ray composition, would reduce theoretical uncertainties and help discriminate between the two cases.

IV Summary and Conclusions

The decay of heavy neutrinophilic DM can produce UHE neutrinos that may scatter off the C$\nu$B. Astrophysical and cosmogenic neutrinos can also undergo such scattering. In the presence of a large local C$\nu$B overdensity, such interactions can significantly modify the UHE neutrino flux, making it possible to test for these overdensities.

In this work, we have calculated projected constraints on the local C$\nu$B overdensity using the future observations of the UHE neutrinos in the next-generation experiment IceCube-Gen2 radio. We considered two scenarios: (i) the UHE neutrino flux originates solely from DM decay, and (ii) the flux includes contributions from DM decay, as well as astrophysical and cosmogenic neutrinos. In case (i), if the neutrino flux is computed in a way that respects $\gamma$-ray limits on DM decay, the projected sensitivity to the C$\nu$B overdensity is generally weaker than current limits, except for large neutrino masses. However, the situation improves significantly, especially for DM masses above approximately $10^{11}$ GeV, if those $\gamma$-ray constraints are relaxed, since a larger DM-induced neutrino flux becomes allowed. In such case, we find that the projected sensitivity to the local C$\nu$B overdensity reaches as low as $\sim 10^{6}$. In case (ii), where we considered UHE neutrinos produced from both DM decays and in astrophysical settings, we generally find the projected sensitivity between $\xi\sim 10^{8}$ and $\xi\sim 10^{10}$, depending on the value of the lightest neutrino mass, if the $\gamma$-ray constraints are respected. The sensitivity further increases and converges to scenario (i) when neutrino fluxes from heavy DM are considered in a scenario with relaxed $\gamma$-ray constraints. Overall, our results point to an interesting novel direction for probing the C$\nu$B clustering in the near future. This is complementary to other direct and indirect searches for the C$\nu$B. At the same time, it could shed some light on the nature of dark matter.

Acknowledgements.

We would like to thank Carlos Argüelles, Jose Carpio, Alex Chen, Bryce Cyr, Saikat Das, Kohta Murase, and Jack Shergold for useful discussions. The work of VB is supported by the U.S. Department of Energy Grant No. DE-SC0025477. The work of WM and BD was partly supported by the U.S. Department of Energy under Grant No. DE-SC0017987. BD was also supported in part by a Humboldt Fellowship from the Alexander von Humboldt Foundation. The work of AMS was supported in part by Department of Energy Grant No. DE-AC02-07CH11359: Neutrino Theory Network Program. We would like to thank the Center for Theoretical Underground Physics and Related Areas (CETUP*) and the Institute for Underground Science at Sanford Underground Research Facility (SURF) for hospitality and for providing a stimulating environment during the 2023, 2024 and 2025 summer workshops, where parts of this work were conducted. We would also like to thank the organizers of PHENO 2022 at University of Pittsburgh, where the initial idea of this project was generated.

Appendix A DM decay constraints

In Fig. 5, we summarize the existing limits in the heavy DM mass–lifetime parameter space, coming from $\gamma$-ray Das et al. (2023); Kachelriess et al. (2018); Chianese et al. (2021); Sarmah et al. (2025), $p$-$\bar{p}$ Das et al. (2023); Sarmah et al. (2025), and neutrino Esmaili et al. (2012); Kachelriess et al. (2018); Argüelles et al. (2023) observations. In all these cases, we have assumed that the DM primarily decays into neutrinos. However, other final states appear unavoidably as secondaries through the electroweak bremsstrahlung. These limits are used to set the benchmarks studied in the main text.

Appendix B SHDM production

There are various nonthermal production mechanisms for a SHDM that can yield the correct relic density. Here we briefly mention two commonly discussed mechanisms:

B.1 Gravitational production

This is the most model-independent mechanism, which only requires the DM coupling to gravity. SHDM is produced by vacuum fluctuations in the expanding spacetime background during or at the end of inflation Chung et al. (1998b); Kuzmin and Tkachev (1999), with a relic density today of order Chung et al. (2001)

[TABLE]

Thus, for a given SHDM mass, an appropriate reheating temperature $T_{\rm RH}$ can be chosen to yield the observed relic density of $\Omega_{\rm DM}h^{2}=0.12$ Aghanim et al. (2020). The exact numerical coefficients in Eq. (5) depend on the model of inflation (e.g., chaotic versus hybrid), but this scaling relation is robust, so long as $m_{\rm DM}<H_{I}$, the Hubble expansion rate during inflation. For $m_{\rm DM}>H_{I}$, the particle production is exponentially suppressed and the result will depend on whether the particle was relativistic or nonrelativistic at the end of inflation Chung (2003).

B.2 Inflaton decay

If the inflaton $\phi$ couples to the SHDM $X$, its decay $\phi\to XX$ can also produce a nonthermal population of DM Allahverdi and Drees (2002a, b); Dev et al. (2014); Bernal and Xu (2021); Bernal et al. (2025), with the current relic density given by

[TABLE]

Thus, a tiny decay branching ratio (Br) of $\phi\to XX$ is enough to generate the observed DM relic density. Note that Eq. (6) assumes perturbative inflaton decay and negligible annihilation or thermalization of DM after production. For the modified expression when backreaction, scatterings, or nonperturbative preheating are important, see Ref. Bernal et al. (2025).

This mechanism can be generalized to the decay of any matterlike component that drives an early-matter-domination phase. A well-known example is moduli decay in supersymmetric theories Moroi and Randall (2000); Gelmini and Gondolo (2006). In string theory, generic features of string compactifications, such as high-scale supersymmetry breaking that gives rise to successful inflation and epochs of early matter domination driven by string moduli, also provide a suitable framework to obtain the correct relic density for SHDM Allahverdi et al. (2013, 2021, 2024).

Appendix C Neutrino transport equation

The neutrino transport equation captures the cosmological evolution of the comoving number of neutrinos per unit energy at time $t$, $\tilde{N}_{A}(E,t)=dN_{A}/dE$ for a given neutrino species $A$. The neutrino source located at a particular redshift $z_{s}$ injects the neutrinos at time $t=0$. The comoving number of neutrinos per unit energy at the redshift $z_{s}$ is obtained using HDMSpectra Bauer et al. (2021). These neutrinos are then propagated until the flux reaches Earth. The complete neutrino transport equation reads Das et al. (2025b)

[TABLE]

where $A\in\{\nu_{1},\nu_{2},\nu_{3},\bar{\nu}_{1},\bar{\nu}_{2},\bar{\nu}_{3}\}$ are the neutrino mass eigenstates, which are connected to the flavor eigenstates of neutrinos $\nu_{\alpha}$ by the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) mixing matrix, assumed to be unitary. The first term on the right-hand side of Eq. (7) is the adiabatic loss term due to the expansion of the Universe. The last term $\tilde{\mathcal{Q}}_{A}(E,t)$ is the source term, accounting for the number of neutrinos of species $A$ with energy $E$ injected by the source per unit time per unit energy at a particular time $t$. We assume $\tilde{\mathcal{Q}}_{A}(E,t)=\delta(t){dN_{A}}/{dE}$. The $\mathcal{A}_{A}(E,t)$-term in Eq. (7) is the total absorption rate of neutrino species $A$ with energy $E$, $\mathcal{B}_{A\to A}$ is the self-production rate of species $A$ with energy $E$ being created from the same particle type $A$ but with different energy $E^{\prime}$, and $\mathcal{C}_{B\to A}$ is the production rate of species $A$ with energy $E$ from other neutrino species $B\neq A$ with energy $E^{\prime}$. Both $\mathcal{B}$ and $\mathcal{C}$ terms provide the reinjection of neutrinos into the flux. For details, see Ref. Das et al. (2025b).