Tomographic constraints on the high-energy cosmic neutrino emission rate

Despite growing efforts to find the sources of high energy neutrinos measured by IceCube, the bulk of the neutrinos remain with unknown origins. We aim to constrain the emissivity of cosmic high-energy neutrinos from extragalactic sources through their correlation with the large-scale structure. We use cross-correlations between the IceCube 10-year dataset and tomographic maps of the galaxy overdensity to place constraints on the bias-weighted high-energy neutrino emissivity out to redshift $z\sim3$. We test two different models to describe the evolution of neutrino emissivity with redshift, a power law model $\propto (1+z)^a$, and a model tracking the star formation history, assuming a simple power law model for the energy injection spectrum. We also consider a non-parametric reconstruction of the astrophysical neutrino emissivity as a function of redshift. We do not find any significant correlation, with our strongest results corresponding to a $1.9 σ$ deviation with respect to a model with zero signal. We use our measurements to place upper bounds on the bias-weighted astrophysical high-energy neutrino emission rate as a function of redshift for different source models. This analysis provides a new probe to test extragalactic neutrino source models. With future neutrino and galaxy datasets we expect the constraining and detection power of this type of analysis analysis to increase.

💡 Research Summary

This paper presents a novel tomographic cross‑correlation analysis between the IceCube ten‑year muon‑track data set and four wide‑area galaxy and quasar catalogs (2MPZ, WIxSC, DECaLS, and Quaia) to constrain the bias‑weighted emissivity of high‑energy astrophysical neutrinos up to redshift z ≈ 3. The authors first formulate the neutrino intensity in terms of a comoving emissivity and derive the angular cross‑power spectrum Cℓ νg using the Limber approximation, which links the neutrino field to the three‑dimensional matter overdensity through a radial kernel that incorporates the neutrino energy spectrum (parameterised by a spectral index β) and the source evolution.

Three parametrisations of the redshift dependence of the bias‑weighted neutrino emission rate are considered: (i) a simple power‑law model ˙ n̄ν(z) = N bν (1+z)^a, where a describes the evolution; (ii) a model that tracks the cosmic star‑formation rate density fSFR(z) as given by Madau & Dickinson (2014); and (iii) a non‑parametric tomographic model in which the emission rate is independently measured at the effective redshifts of the four galaxy samples. In each case the combination N bν (or Nρ bν for the energy‑weighted emissivity) is treated as a free parameter because the source bias bν is uncertain but expected to be of order unity to a few.

The IceCube events are turned into full‑sky maps of neutrino number density and energy density. Each event is weighted by the inverse of the detector’s effective area Aeff(δ, E) and by its reconstructed energy, and a Gaussian angular kernel approximates the per‑event point‑spread function. Atmospheric muon and neutrino backgrounds are mitigated by restricting to declinations > –5°. The galaxy maps are constructed from the catalog redshift distributions dp/dz and linear galaxy biases bg taken from the literature; the corresponding radial kernels qg(χ) are then inserted into the cross‑power calculation.

Cross‑power spectra are estimated with NaMaster on HEALPix maps (Nside = 128) over multipoles ℓ ≈ 10–200, where the Limber approximation is valid and large‑scale anisotropies dominate. Covariance matrices are generated from 10 000 Gaussian simulations, and Bayesian parameter inference is performed with an MCMC sampler. The analysis finds no statistically significant detection; the strongest excess corresponds to only 1.9 σ above the null hypothesis.

Upper limits are derived for the bias‑weighted emission rates. In the power‑law case the data constrain the evolution index to a ≲ –1.5 (95 % confidence) and the normalisation to N bν ≲ 10⁻⁹ Mpc⁻³ yr⁻¹ (for number density) or Nρ bν ≲ 5 × 10⁻¹⁰ erg Mpc⁻³ yr⁻¹ (for energy density). The star‑formation‑rate‑tracking model yields N bν ≲ 5 × 10⁻¹⁰ Mpc⁻³ yr⁻¹ at z = 0, implying that the contribution of sources following the SFR history is limited to a few percent of the total IceCube flux. The tomographic reconstruction provides redshift‑bin‑by‑bin limits of order 10⁻⁹–10⁻¹⁰ Mpc⁻³ yr⁻¹, again indicating that the cumulative neutrino emissivity is well below the observed diffuse flux.

The authors discuss the implications of these limits. Because the cross‑correlation signal scales with the product b × N, the constraints are robust against uncertainties in the source bias as long as bν remains within the plausible range (1–5). The non‑detection is attributed primarily to the limited statistics of the IceCube muon‑track sample and the modest depth of the current galaxy catalogs. However, the methodology demonstrates that future facilities—IceCube‑Gen2, KM3NeT, and next‑generation optical/infrared surveys such as LSST, Euclid, and the Roman Space Telescope—will dramatically improve sensitivity. With larger neutrino event samples, better angular resolution, and deeper, more complete tomographic tracers, the signal‑to‑noise of the cross‑power spectrum could increase by an order of magnitude, potentially leading to a definitive detection.

In summary, this work establishes tomographic cross‑correlation as a powerful, model‑independent probe of the high‑energy neutrino source population. It provides the first quantitative upper bounds on the bias‑weighted neutrino emissivity as a function of redshift and outlines a clear path toward future detections that will illuminate the long‑standing mystery of the origins of the diffuse astrophysical neutrino background.