Pith. sign in

REVIEW 3 major objections 3 minor 20 references

Bayesian Deep-stacking for High-energy Neutrino Searches

T0 review · 3 major / 3 minor · reviewed 2026-08-09 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Bayesian deep-stacking beats maximum likelihood on faint neutrino sources.

desk verdict Clear, useful framework for Bayesian deep-stacking, but the headline sensitivity gain is inflated by a prior that uses the true signal count; the main result needs a redo with a data-derived prior. read the letter →

arxiv 2502.01452 v2 pith:JMTAYRFE submitted 2025-02-03 astro-ph.HE

I. Bartos , M. Ackermann , M. Kowalski This is my paper · ORCID
classification astro-ph.HE
keywords Bayesianinferencedeep-stackinghigh-energyneutrinossourcecatalogstackinghierarchicalmodelmaximumlikelihoodcomparisonluminosityuncertaintyneutrinoastronomy
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tries to establish that deep-stacking, combining large numbers of individually undetectable neutrino source candidates into one statistical search, is best carried out with Bayesian inference because the dominant uncertainty is the poorly known mapping from a source's observable properties to its expected neutrino flux. In a realistic simulation of roughly 120 cataloged sources, it finds that when per-source flux expectations are scattered by 1 dex relative to the true luminosities, Bayesian deep-stacking keeps the median detection significance above 3σ, while the standard maximum-likelihood stack drops below it; the authors describe the gain as about 2σ of added sensitivity, equivalent to doubling the effective signal strength. The same framework reconstructs population properties more accurately than maximum likelihood in a simplified 1000-source example, because marginalizing over individual source luminosities uses all information without overfitting. If the claim holds, future catalog-driven neutrino searches would gain sensitivity by adopting Bayesian odds-ratio tests rather than frequentist likelihood ratios.

What carries the argument

The machinery is the hierarchical Bayesian odds-ratio test built on the Poisson product likelihood $p(x|\xi)\propto e^{-\hat n_{\rm tot}}\prod_i [\sum_j \hat n_j S_j(x_i) + \hat n_{\rm bg} B(x_i)]$, where $\hat n_j$ is the expected event count of catalog source $j$, $S_j$ is the signal spatial probability, and $B$ is the background probability. Deep-stacking is the strategy of including all such sources, even very distant and faint ones, so that the search targets population detection rather than individual detections. The key move is to put a log-normal prior on each $\hat n_j$ with width $\sigma_S$ characterizing flux-model scatter, marginalize over the $\hat n_j$ inside the odds ratio, and then reinterpret that odds ratio as a frequentist test statistic by calibrating it on background-only simulations. A practical simplification makes the marginalization tractable: when the average source spacing is much larger than the detector angular resolution, each event contributes only to the nearest source's signal term, reducing the computation to a series of one-dimensional integrals. The prior on the total signal scale, taken uniform over 0 to 2 times the true expected signal, is what carries the sensitivity claim.

What would settle it

Rerun the realistic simulation exactly as described but with a prior on the total expected signal that is not centered on the truth, for example uniform over a decade-wide range whose midpoint is off by a factor of three, and compare Bayesian deep-stacking with the source-normalization maximum-likelihood method at 1 dex scatter; if the Bayesian advantage falls below the claimed roughly 2σ, the truth-centered prior is the load-bearing element.

Watch

Extended reading notes

Core claim

The central claim is that a hierarchical Bayesian treatment of source-stacking searches for high-energy neutrinos is more sensitive and more accurate than the standard frequentist likelihood-ratio approach whenever individual source flux expectations are uncertain. The paper writes the unbinned Poisson likelihood for a catalog search and treats each source's expected event count as a nuisance parameter with a log-normal prior whose width encodes the uncertainty of the flux model. It then forms the odds ratio between the signal and background hypotheses, marginalizes over the nuisance parameters, and uses the odds ratio itself as a test statistic calibrated on background-only pseudo-data. In the realistic simulation, the Bayesian method outperforms both a standard sample-normalization maximum-likelihood search and a source-normalization search with per-source nuisance parameters, and it is the only method that reaches average significance above 3σ for 1 dex flux scatter. The paper also shows, in a simplified scenario with 1000 cataloged sources, that the Bayesian posterior on the average source luminosity is less biased and has smaller error bars than maximum-likelihood estimates.

Load-bearing premise

The sensitivity demonstration assumes the analysis uses a prior on the total expected signal that is centered on the true simulated value (uniform from zero to twice that value), information a real search would not have; if that prior is misspecified, the reported 2σ advantage may shrink.

Editorial extensions

If this is right

  • Catalog-based neutrino stacking searches should prefer Bayesian odds-ratio tests over sample-normalization maximum likelihood when per-source flux predictions are unreliable, because the standard method's significance degrades as flux scatter grows.
  • A real search can report posterior credible intervals on population parameters such as average luminosity, giving physical interpretation of the stacked population rather than only a detection p-value.
  • The gain is expected to grow with the number of sources in the catalog, since marginalization avoids overfitting the many per-source nuisance parameters that maximum-likelihood searches must fit.
  • If source density becomes high enough that point-spread functions overlap, the nearest-source approximation fails and full multi-dimensional marginalization would be required, which the paper identifies as the current computational bottleneck.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • A natural extension, not explored in the paper, is to apply the same odds-ratio construction to other astrophysical transient searches or to gravitational-wave source catalogs, where per-source distance and mass uncertainties play a role similar to neutrino flux-model uncertainty.
  • A practical prescription suggested but not tested is to set the total-signal prior from an independent measurement such as the diffuse neutrino flux, which would remove the truth-centered-prior dependence; testing this would show how much of the reported gain survives in a real search.
  • Treating the spectral index as a hyperparameter with its own prior, rather than fixing it to the injected value, is a natural extension that could change the energy-bin weights and hence the sensitivity comparison.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 3 minor

Summary. This paper proposes a Bayesian framework for 'deep-stacking' searches for high-energy neutrinos from large catalogs of faint candidate sources. The method combines a Poisson likelihood over reconstructed events with priors on per-source expected event counts, marginalizes over nuisance parameters, and uses a Bayes factor / odds ratio as a test statistic calibrated with background-only simulations. The authors demonstrate the approach in two settings: a simplified 1000-source example reconstructing the average source luminosity, and a more realistic simulation comparing Bayesian deep-stacking with two maximum-likelihood stacking variants under increasing scatter between the injected source fluxes and the fluxes assumed in the analysis. They report that the Bayesian method outperforms the frequentist approaches, with a factor-50 improvement in p-value and roughly 2σ added sensitivity at 1 dex scatter.

Significance. If the reported sensitivity advantage survives a fully fair comparison, the paper would make a useful contribution to neutrino multi-messenger searches and population inference. The hierarchical Bayesian formulation is natural for this problem, the use of calibrated test-statistic distributions to compare Bayesian and frequentist methods is commendable, and the authors are transparent about computational limitations. The main methodological ideas — incorporating source-flux uncertainties as priors and marginalizing rather than profiling — are clearly presented and potentially valuable. However, the central quantitative claim in Section 4 is currently supported by a prior that uses the true simulated signal count as its normalization anchor, and the simplified example in Section 3 contains a prior expression that is not a uniform distribution as claimed. These issues are load-bearing for the paper's main demonstrations, so the significance of the reported sensitivity gain is not yet established.

major comments (3)
  1. [Section 4.4] The prior on the total signal normalization is informative in a way that leaks the simulation truth into the analysis. The text states: 'A uniform prior was used for 0 ≤ Σ_j n̂*_j / n̂_sig ≤ 2', where n̂_sig is defined in Section 4.2 as the expected number of signal events from the source population in the simulated data. An analyst performing a real search would not know n̂_sig; at best the total normalization could be bounded by the observed diffuse flux, which would give a far wider and less informative prior. Because the claimed advantage at 1.0 dex scatter (p ≈ 8.9×10^-5 for Bayesian versus p ≈ 1.4×10^-2 for sample-normalization maximum likelihood in Fig. 2) occurs precisely where the total-scale prior matters most, the reported sensitivity gain may be driven by this truth-centered prior rather than by Bayesian marginalization. The authors should either re-run the comparison with a prior anchored to observables (e.g., a broad prior derived from the diffuse astrophysical neutrino flux) or demonstrate that the reported gains are robust to the choice of prior scale.
  2. [Section 3.2] The prior for the average luminosity is stated as 'We adopted a non-informative uniform distribution for the luminosity: π(ξvar) = 2 exp(L0 + 2σL).' This expression is not uniform in L0: it grows exponentially with L0 and depends on the parameter being estimated and on σL. As written it is not a valid probability density over L0 (nor a uniform density over log L0). This is not a minor typographical issue, because the simplified example in Section 3 is one of the paper's two demonstrations of Bayesian parameter reconstruction and of superiority over maximum likelihood. The authors should correct the prior specification, state the actual support of the prior, and re-run the Monte Carlo results in Fig. 1 to verify that the claimed unbiasedness and improved performance still hold.
  3. [Section 4.4] The scatter width σS used in the log-normal per-source priors is set to the same value as the scatter injected into the simulated data. The text says, 'we computed various scenarios of scattering between the injected source flux and the n̂_j used in the analysis model using the same set of widths, i.e., σS = (0, 0.3, 0.5, 1.0)'. In a real analysis σS is not known a priori; it would need to be estimated from data, marginalized over, or treated with a hyperprior. Since this parameter controls how much the likelihood can adapt to source-by-source deviations, the comparison as presented is optimistic for both the Bayesian and the source-normalization methods. The authors should show that the reported ranking of methods is robust to a mismatch between the assumed and true σS, for example by evaluating the Bayesian analysis with a moderately mis-specified σS.
minor comments (3)
  1. [Throughout] There are several typographical errors that should be corrected: 'a called deep-stacking' in the abstract, 'frequentiest' in Section 4, 'nuissance' in Section 4.3, 'Incoporating' and 'strenght' in the Conclusion, and 'Jeffries prior' in Section 2.3 (should be Jeffreys).
  2. [Section 4.3] The sentence 'To match this challenge, the likelihood calculations for this manuscript have been optimized using the PyTorch package' appears without a preceding full stop in the submitted text; please fix the punctuation.
  3. [Section 2.6] In Eq. (2.6), the notation p(x|Ltot > 0, ξ) is used, but the right-hand side should be a marginal likelihood with the source luminosities integrated out; this is presumably implied by the context but should be stated explicitly for clarity.

Circularity Check

1 steps flagged · score 7.0 of 10

Realistic-simulation sensitivity claim is anchored to a prior normalized by the true simulated signal count n_sig.

  1. self definitional [Section 4.2 (Eq. 4.1) and Section 4.4 (Eqs. 4.4-4.5)]
    "where the star in ˆn∗ sig indicates that this parameter is a fit parameter (to distinguish it from ˆnsig, which is the expected number of signal events from the source in the simulated population) ... A uniform prior was used for 0≤ P j ˆn∗ j /ˆnsig ≤2, while a log-normal shape was used for the individual priors for ˆnj."

    The Bayesian prior on the total fitted signal is defined in units of n_sig, the true expected number of signal events in the simulated population. This leaks the simulation truth into the analysis: the prior is centered on the actual total signal and restricts it to at most twice that value. In a real search the analyst would not know n_sig; the paper's own Section 2.3 suggests using the observed diffuse flux as an upper bound, not the injected true count. The reported factor-of-50 p-value advantage at 1.0 dex scatter (p=8.9e-5 vs 1.4e-2) is therefore partly an artifact of giving the Bayesian method the correct total normalization, while the frequentist sample-normalization approach must infer that scale from the data.

full rationale

The main circularity is confined to the realistic-simulation comparison in Section 4. Equation 4.1 explicitly defines n_sig as the expected signal count in the simulated population, and Section 4.4 places a uniform prior on the total fitted signal normalized by that same true n_sig. This is a truth-informed prior: it supplies the Bayesian analysis with the overall signal scale that a real search would have to estimate or bound from observables. Consequently, the headline sensitivity advantage in Fig. 2—especially the factor of 50 in p-value for 1.0 dex scatter—is partly built into the prior rather than earned by Bayesian marginalization. The simplified reconstruction example in Section 3 is self-contained and does not share this flaw; it uses a standard hierarchical prior. No self-citation chain or imported uniqueness theorem is load-bearing, and the computational-limitation note in Section 4.4 is not a circularity issue. The score reflects that the central realistic-case claim is compromised by the truth-anchored prior, while the general framework retains independent content.

Assumptions & free parameters 3 free parameters · 4 assumptions · 0 invented entities

The central claim rests on standard Bayesian inference axioms, plus domain-specific assumptions about the source luminosity distribution and telescope response. The main concern is the free parameter setting the prior upper bound to the true signal count, which leaks information into the analysis. No new physical entities are postulated.

free parameters (3)
  • Prior upper bound on total signal normalization = 2 × n̂_sig (true value)
    The Bayesian analysis uses a uniform prior 0 ≤ Σ_j n̂*_j / n̂_sig ≤ 2 in Section 4.4. Since n̂_sig is the simulated truth, this prior is centered on the true signal count and biases the sensitivity comparison.
  • Prior scatter width σ_S = 0, 0.3, 0.5, 1.0 dex
    This width of the log-normal prior on individual source normalizations is chosen by the authors for each scenario; it controls the amount of source-flux uncertainty and is not inferred from data, but it is a modeling choice that affects the comparison.
  • Prior on average luminosity in simplified example = π(ξvar)=2 exp(L0+2σL)
    In Section 3.2 the authors claim to use a non-informative uniform prior, but the expression given is neither uniform nor normalizable; this is a mis-specified prior that may affect the toy demonstration.
assumptions (4)
  • standard math Poisson likelihood for detected neutrino counts from sources and background.
    Used throughout Section 2.4 and in Eq. 2.3; this is a standard statistical model for rare-event counting.
  • domain assumption Spatial and energy signal PDFs are known and separable.
    Eq. 2.4 assumes a Gaussian spatial term and a power-law energy term with common spectral index γ. This is a modeling assumption about neutrino telescope response.
  • domain assumption Source luminosities are drawn from a log-normal distribution with known σ_L (or σ_S).
    Used in Sections 3.1 and 4.4; the log-normal form and its width are assumed, not derived from data.
  • ad hoc to paper Only the closest source to each event contributes significantly to the signal PDF.
    Stated in Section 4.4 to reduce the integration dimension; this approximation is valid only when sources are well separated relative to the angular resolution, limiting scalability.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Bayesian Deep-stacking for High-energy Neutrino Searches." pith.science (2026). https://pith.science/paper/JMTAYRFE

@misc{pith2026250201452,
  author       = {Pith},
  title        = {Pith review of: Bayesian Deep-stacking for High-energy Neutrino Searches},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/JMTAYRFE}},
  note         = {Machine review of arXiv:2502.01452}
}
read the original abstract

Following the discovery of the brightest high-energy neutrino sources in the sky, the further detection of fainter sources is more challenging. A natural solution is to combine fainter source candidates, and instead of individual detections, aim to identify and learn about the properties of a larger population. Due to the discreteness of high-energy neutrinos, they can be detected from distant very faint sources as well, making a statistical search benefit from the combination of a large number of distant sources, a called deep-stacking. Here we show that a Bayesian framework is well-suited to carry out such statistical probes, both in terms of detection and property reconstruction. After presenting an introductory explanation to the relevant Bayesian methodology, we demonstrate its utility in parameter reconstruction in a simplified case, and in delivering superior sensitivity compared to a maximum likelihood search in a realistic simulation.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

20 extracted references · 2 canonical work pages

  1. [1]

    IceCube Collaboration,Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,Science342(2013) 1242856 [1311.5238]

    arXiv 2013

  2. [2]

    Aartsen, M

    IceCube Collaboration, M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al.,Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A,Science361(2018) eaat1378 [1807.08816]

    arXiv 2018

  3. [3]

    Aartsen, M

    IceCube Collaboration, M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al.,Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert,Science361(2018) 147 [1807.08794]. – 13 –

    arXiv 2018

  4. [4]

    Abbasi, M

    IceCube Collaboration, R. Abbasi, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al., Evidence for neutrino emission from the nearby active galaxy NGC 1068,Science378(2022) 538 [2211.09972]

    arXiv 2022

  5. [5]

    Kadler, F

    M. Kadler, F. Krauß, K. Mannheim, R. Ojha, C. Müller, R. Schulz et al.,Coincidence of a high-fluence blazar outburst with a PeV-energy neutrino event,Nature Physics12(2016) 807 [1602.02012]

    arXiv 2016

  6. [6]

    Neronov, D

    A. Neronov, D. Savchenko and D.V. Semikoz,Neutrino Signal from a Population of Seyfert Galaxies,Phys. Rev. Lett.132(2024) 101002 [2306.09018]

    arXiv 2024

  7. [7]

    E. Kun, I. Bartos, J.B. Tjus, P.L. Biermann, A. Franckowiak, F. Halzen et al.,Possible correlation between unabsorbed hard x rays and neutrinos in radio-loud and radio-quiet active galactic nuclei,Phys. Rev. D110(2024) 123014 [2404.06867]

    arXiv 2024

  8. [8]

    Abbasi, M

    R. Abbasi, M. Ackermann, J. Adams, S.K. Agarwalla, J.A. Aguilar, M. Ahlers et al.,Search for neutrino emission from hard X-ray AGN with IceCube,arXiv e-prints(2024) arXiv:2406.06684 [2406.06684]

    arXiv 2024

Show all 20 references
  1. [9]

    Bartos, D

    I. Bartos, D. Veske, M. Kowalski, Z. Márka and S. Márka,The IceCube Pie Chart: Relative Source Contributions to the Cosmic Neutrino Flux,Astrophys. J921(2021) 45 [2105.03792]

    2021 arXiv

  2. [10]

    Aartsen, M

    M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, M. Ahrens et al., Time-Integrated Neutrino Source Searches with 10 Years of IceCube Data,Phys. Rev. Lett.124 (2020) 051103 [1910.08488]

    2020

  3. [11]

    Kowalski, M

    M. Kowalski, M. Ackermann and I. Bartos,The promise of deep-stacking for neutrino astronomy,arXiv e-prints(2025) arXiv:2501.10213 [2501.10213]

    2025 arXiv

  4. [12]

    Mandel, W.M

    I. Mandel, W.M. Farr and J.R. Gair,Extracting distribution parameters from multiple uncertain observations with selection biases,MNRAS486(2019) 1086 [1809.02063]

    2019 arXiv

  5. [13]

    Abbott, T.D

    R. Abbott, T.D. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams et al.,Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog,Astrophys. J Lett.913(2021) L7 [2010.14533]

    2021 arXiv

  6. [14]

    Kovetz, I

    E.D. Kovetz, I. Cholis, P.C. Breysse and M. Kamionkowski,Black hole mass function from gravitational wave measurements,Phys. Rev. D95(2017) 103010 [1611.01157]

    2017 arXiv

  7. [15]

    W.M. Farr, S. Stevenson, M.C. Miller, I. Mandel, B. Farr and A. Vecchio,Distinguishing spin-aligned and isotropic black hole populations with gravitational waves,Nature548(2017) 426 [1706.01385]

    2017 arXiv

  8. [16]

    Capel, D.J

    F. Capel, D.J. Mortlock and C. Finley,Bayesian constraints on the astrophysical neutrino source population from IceCube data,Phys. Rev. D101(2020) 123017 [2005.02395]

    2020 arXiv

  9. [17]

    Capel, J

    F. Capel, J. Kuhlmann, C. Haack, M. Ha Minh, H. Niederhausen and L. Schumacher,A Hierarchical Bayesian Approach to Point-source Analysis in High-energy Neutrino Telescopes, Astrophys. J976(2024) 127 [2406.14268]. [18]IceCubecollaboration,The contribution of Fermi-2LAC blazars ...

    2024 arXiv

  10. [20]

    Piessens, E

    R. Piessens, E. de Doncker-Kapenga, C.W. Überhuber and D. Kahaner,QUADPACK: A Subroutine Package for Automatic Integration, Springer-Verlag (1983). [21]IceCubecollaboration,Characteristics of the diffuse astrophysical electron and tau neutrino flux with six years of IceCube hi...

    1983

  11. [22]

    Aartsen, M

    M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, M. Ahrens et al.,The IceCube Neutrino Observatory: instrumentation and online systems,Journal of Instrumentation12(2017) P03012 [1612.05093]. [23]KM3Netcollaboration,Letter of intent for KM3NeT 2.0,J. Phys. G43(201...

    2017 arXiv

  12. [24]

    Abbasi et al.,Improved Characterization of the Astrophysical Muon–neutrino Flux with 9.5 Years of IceCube Data,Astrophys

    R. Abbasi et al.,Improved Characterization of the Astrophysical Muon–neutrino Flux with 9.5 Years of IceCube Data,Astrophys. J.928(2022) 50 [2111.10299]. – 15 – 0 50 100 150 200 100 101 102 103 104 Number of realizationsp-value = 1.5 × 10 8 TSmed = 30.7 no scatter max. LH: sam...

    2022

Pith tools

Reviewed August 9, 2026 · model on record in the stance chip above.