REVIEW 6 minor 25 references
On the difficulty of capturing the distribution function of neutrinos in neutron star merger simulations
T0 review · 0 major / 6 minor · reviewed 2026-08-16 · deepseek-v4-flash
Pith's one-line read Current Monte Carlo merger codes cannot resolve the neutrino distribution function fν; a single packet can make it ~10^5.
desk verdict A solid, honest methods paper that quantifies why current Monte Carlo neutrino transport in mergers cannot resolve f_nu and gives concrete packet-count targets, with the main caveats flagged by the author himself. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing identity is the packet-count relation $N_{\mathrm{target}} = (f_{\mathrm{true}}/\sigma_f)^2$. It follows from writing the maximum energy in a phase-space bin as $E_{\max}(D) = \frac{1}{(hc)^3} \Delta\Omega \Delta V \frac{\epsilon_1^4 - \epsilon_0^4}{4}$ and using the Monte Carlo shot-noise estimate $\sigma_f = \sqrt{\langle N\rangle} e_p/E_{\max} = \sqrt{f_{\mathrm{true}} e_p/E_{\max}}$. This identity does the work of the whole paper: it converts a target absolute error in $f_\nu$ into a required number of packets per bin, exposes the inefficiency of constant-energy or constant-number packet weights, and is inverted to give the optimal packet weight $n_p \propto (\epsilon_1^3 - \epsilon_0^3)/f_{\mathrm{true}}$ so that the absolute error is roughly constant across bins. The paper also introduces the helper estimate $E_{\mathrm{true}} = \eta/(\kappa_a + \kappa_{\mathrm{floor}}) \Delta V/c$ for the equilibrium neutrino energy in a cell, which sets the normalization for the concrete packet-count figures.
What would settle it
Run a Monte Carlo merger simulation, freeze the fluid and metric evolution, and repeatedly resample the neutrino emission with a different random seed while keeping all other inputs identical. In each run, count the packets in a fixed low-energy phase-space bin of a hot cell. If the empirical variance of that count over realizations is close to the mean, the paper's $N_{\mathrm{target}} = (f_{\mathrm{true}}/\sigma_f)^2$ requirement stands; if it is substantially larger, as it should be immediately after packet splitting when sibling packets are correlated, then reaching $\sigma_f = 0.1$ requires more than 100 packets per bin and the paper's feasibility claim is too optimistic.
Extended reading notes
Core claim
The paper's central claim is that the shot noise inherent in Monte Carlo transport, not any missing reaction physics, is what blocks access to $f_\nu$. In a Monte Carlo code the distribution function is represented as a sum of infinitely narrow spikes, one per packet, so $f_\nu$ can only be defined by averaging over a phase-space domain $D$. If the expected number of packets in $D$ is $\langle N\rangle$, the sampling noise is roughly $\sigma_N \approx \sqrt{\langle N\rangle}$, which translates into an error $\sigma_f = \sqrt{f_{\mathrm{true}} e_p/E_{\max}}$ in the inferred occupation number, where $e_p$ is the energy carried by one packet and $E_{\max}$ is the total energy capacity of $D$. Requiring $\sigma_f = 0.1$ in a bin where $f_{\mathrm{true}} \approx 1$ gives $N_{\mathrm{target}} = (f_{\mathrm{true}}/\sigma_f)^2 \approx 100$ packets per bin. Because real merger codes put tens to hundreds of packets in a whole cell and concentrate those packets in high-energy bins, the low-energy bins where $f_\nu$ is largest are left with zero or one packet, producing $f_\nu$ estimates of 0 or about $10^5$. The paper then works out the cost of fixing this: with the optimal per-bin weighting $n_p \propto (\epsilon_1^3 - \epsilon_0^3)/f_{\mathrm{true}}$, hot remnant regions still need about $10^3$ packets per cell for $\sigma_f = 0.1$, while the fixed-number-of-neutrinos weighting used today would need $10^5$ to $10^7$ packets per cell; even the optimal route requires smoothing $f_\nu$ over neighboring cells or coarser energy bins.
Load-bearing premise
The packet-count estimates assume that the number of neutrino packets in a phase-space region scatters like the square root of its average count, with packets acting independently; the paper concedes in its appendix that this is not guaranteed, especially after packet splitting, and if the scatter is larger, every required count in the paper goes up.
Editorial extensions
If this is right
- Reaching $\sigma_f = 0.1$ where $f_\nu \approx 1$ requires about 100 packets in every energy or angular bin; with 16 energy bins and 3 neutrino species, that is about $10^3$ packets per cell, above what current three-dimensional merger codes supply.
- Any fixed packet weight, whether constant energy or constant number of neutrinos, leaves low-energy neutrinos grossly underresolved, and existing $f_\nu$ estimates in low-energy bins are therefore the worst offenders, jumping between 0 and about $10^5$.
- Adopting $n_p \propto (\epsilon_1^3 - \epsilon_0^3)/f_{\mathrm{true}}$ lowers the packet count for a fixed $\sigma_f$ by several orders of magnitude in hot regions, but it increases shot noise in the fluid coupling because fewer, heavier packets carry the high-energy neutrino signal.
- Smoothing $f_\nu$ over neighboring cells, about 27 cells in three dimensions, and merging low-energy bins is the cheapest path toward $\sigma_f \leq 0.1$ without raising total packet numbers.
- Until these changes are made, simulations should keep avoiding explicit $f_\nu$ in interaction rates, which means pair production and inelastic electron scattering remain approximate in merger codes.
Reading between the lines
- My inference: the paper's packet-count arithmetic implies that the quoted 10 to 20 percent agreement between Monte Carlo and moment codes on luminosities and ejecta does not test whether $f_\nu$ is correct; observables that weight the tail of the distribution, such as high-energy neutrino spectra or pair-annihilation heating, could be far more sensitive to the missing resolution.
- My inference: the optimal weighting $n_p \propto (\epsilon_1^3 - \epsilon_0^3)/f_{\mathrm{true}}$ resembles an importance-sampling prescription and could be combined with a control-variate approach that uses the equilibrium distribution $f_{\mathrm{eq}}$ as a baseline, potentially reducing the variance below the paper's $\sqrt{N}$ estimate while keeping the same packet counts.
- My inference: a direct empirical test of the paper's variance assumption is easy to design by running the same fluid snapshot with several independent Monte Carlo emission sequences and measuring the scatter in $f_\nu$ per bin; if the scatter is dominated by correlations inherited from packet splitting, the required counts rise and the paper's feasibility conclusion weakens.
- My inference: the paper's normalization $E_{\mathrm{true}} = \eta/(\kappa_a+\kappa_{\mathrm{floor}}) \Delta V/c$ is most reliable in near-equilibrium regions, so the estimated packet counts are least trustworthy in free-streaming and semi-transparent regions, exactly where $f_\nu$ matters most for the kilonova ejecta; the practical bottleneck may therefore be even tighter than the paper's central-
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. This manuscript argues that Monte Carlo neutrino transport in neutron star merger simulations cannot currently provide useful estimates of the neutrino distribution function f_nu from individual grid cells. It begins with phase-space counting of neutrino quantum states and derives the maximum neutrino number and energy in a region (Eqs. 3-7). It then adopts a Poisson shot-noise model, sigma_N ~ sqrt(N) for N >> 1, to obtain the central scaling sigma_f ~ sqrt(f_true e_p / E_max) (Eq. 15) and the target packet number N_target = (f_true / sigma_f)^2 (Eq. 17), with a correction for low occupancy (Eq. 20). These requirements are evaluated for merger-relevant thermodynamic conditions with the SFHo equation of state and 16 energy bins, leading to the conclusion that O(10^3) packets per cell are needed under f-adapted weighting, while fixed-number-of-neutrinos weighting requires 10^5-10^7 packets per cell (Figs. 1-2). The paper proposes a weighting scheme n_p proportional to (epsilon_1^3 - epsilon_0^3) / f_true (Eq. 23) and discusses smoothing over space and energy. Appendix A candidly discusses the limitations of the sqrt(N) variance assumption, including correlated packets after splitting and fluid back-reaction.
Significance. The paper addresses an important and timely problem: the fermionic nature of neutrinos and the practical inability of Monte Carlo merger simulations to resolve f_nu. The central diagnostic---that current simulations are orders of magnitude away from resolving f_nu from single-cell information---is robust because it depends only on order-of-magnitude packet counts and basic phase-space counting. The specific scaling formulas and target packet counts are useful design targets for future code development. A notable strength is the paper's transparency: it derives the requirements from simple, self-contained arguments, uses no fitted parameters, and clearly lists the limitations of its variance ansatz in Appendix A. The paper does not overclaim; it explicitly frames the estimates as guidance rather than as substitutes for convergence studies. The main residual concern, that the packet-count requirements are ideal-sampling lower bounds, slightly tempers the positive feasibility claim but does not affect the negative diagnostic.
minor comments (6)
- [Sec. III.A and Appendix A] Equation (17) and the packet counts in Figs. 1-2 are presented as requirements, but Appendix A correctly notes that sigma_N ~ sqrt(N) requires independent packet draws and can fail for correlated packets after splitting. Please state explicitly in Section III.A that all quoted packet counts are lower bounds under ideal independent sampling, especially since Section III.B advocates adaptive packet splitting; this caveat should appear near the main quantitative results, not only in the appendix.
- [Sec. III.B, Eq. (24)] The normalization E_true = eta / (kappa_a + kappa_floor) Delta V / c is an approximation that combines local equilibrium and free-streaming limits. Because Fig. 2 focuses on low densities where this approximation is least controlled, please add a sentence quantifying the expected bias or noting that the low-density packet counts inherit this approximation; as written the reader cannot tell how robust the '~100 packets' estimate for disks is to the choice of kappa_floor.
- [Sec. I and Abstract] The motivating claim that a single Monte Carlo packet makes f_nu jump from 0 to ~10^5 in the author's simulations is not demonstrated in this manuscript. A short quantitative example using the stated packet energy, grid resolution, and Eq. (7) would make the motivation self-contained and would let readers verify the order of magnitude without consulting reference [22].
- [Sec. II.A, footnote 1] The notation distinguishing lower-case n/e, upper-case N/E, and script N for packet number is easy to lose, particularly because Eq. (14) uses N for the number of packets while Eq. (3) uses N for the number of neutrinos. Using a distinct symbol such as N_p for packet number throughout would remove ambiguity.
- [Sec. III.B, Fig. 1 caption] The color scale labels should specify whether the plotted quantity is the total number of packets per grid cell summed over energy bins and species or the number per energy bin; the text at the end of Section III.B ('we need at least 100 packets per cell') is also unclear in this regard and appears inconsistent with the earlier O(1000) estimate.
- [Throughout] There are several typos and formatting issues: 'probabily' in Sec. II.B, 'manusript' in Sec. III.A, and 'weighing' where 'weighting' is meant in Sec. IV; these should be corrected in a final pass.
Circularity Check
No significant circularity: the packet-count requirements follow from phase-space counting and a stated accuracy target, with no fitted parameter or self-citation forcing the result.
full rationale
The paper's central derivation is self-contained. The required packet count N_target = (f_true/σ_f)^2 (Eq. 17) follows directly from the Poisson variance estimate σ_N ~ sqrt(⟨N⟩), the definition of the packet estimate σ_f = sqrt(N) e_p / E_max (Eq. 15), and a chosen target absolute error σ_f. No parameter is fitted to data and then renamed as a prediction; the target accuracy is an external specification. The proposed weighting np ∝ (ε1^3−ε0^3)/f_true (Eq. 23) does require knowledge of f_true, but the paper explicitly acknowledges this: 'The scheme additionally requires an estimate of f_true; feq may be a reasonable choice here.' This is an importance-sampling design requirement, not a circular derivation of fν from fν. The variance ansatz is admittedly approximate; Appendix A states that σ_N ∼ sqrt(N) 'is only rigorously defined under very narrow constraints' and warns about correlated packets after splitting. That caveat weakens the precision of the quantitative feasibility estimate but does not make the argument circular. Self-citations (e.g., [9], [11], [22], [25]) are used as concrete examples of current simulation packet counts, grid spacings, energy binning, and reaction rates; they are motivational context, not load-bearing justification of the counting formulas. The headline observation that a single Monte Carlo packet can give fν ~ 10^5 is an empirical statement about the author's existing simulations, and the theoretical section reproduces the difficulty from independent phase-space and sampling considerations. No load-bearing step reduces to its own input by construction, and no external result is replaced by a self-citation chain.
Assumptions & free parameters
free parameters (2)
- sigma_f = 0.1 =
0.1
- kappa_floor = c^2/(G M_sun) =
c^2/(G M_sun)
assumptions (4)
- domain assumption The number of Monte Carlo packets in a region of phase space has variance sigma_N approximately sqrt(⟨N⟩) for ⟨N⟩ much larger than 1, and sigma_N less than or about 1 for ⟨N⟩ of order 1.
- standard math Neutrinos are massless, have fixed chirality, and satisfy 0 < fν < 1; spacetime is locally orthonormal within the phase-space region D.
- domain assumption E_true can be approximated by eta/(kappa_a + kappa_floor) Delta V / c, with eta the emissivity and kappa_a the absorption opacity.
- domain assumption The true distribution function f_true can be approximated by the equilibrium Fermi-Dirac distribution f_eq when setting packet weights.
Cite this review
Pith. "Pith review of On the difficulty of capturing the distribution function of neutrinos in neutron star merger simulations." pith.science (2026). https://pith.science/paper/SL66WZ4K
@misc{pith2026250421822,
author = {Pith},
title = {Pith review of: On the difficulty of capturing the distribution function of neutrinos in neutron star merger simulations},
year = {2026},
howpublished = {\url{https://pith.science/paper/SL66WZ4K}},
note = {Machine review of arXiv:2504.21822}
}
abstract
The collision of two neutron stars is a rich source of information about nuclear physics. In particular, the kilonova signal following a merger can help us elucidate the role of neutron stars in nucleosynthesis, and informs us about the properties of matter above nuclear saturation. Approximate modeling of neutrinos remains an important limitation to our ability to make predictions for these observables. Part of the problem is the fermionic nature of neutrinos. By the exclusion principle, the expected value $f_\nu$ for the number of neutrinos in a quantum state is at most 1. Any process producing neutrinos is suppressed by a blocking factor $(1-f_\nu)$. Recent simulations focused on neutrino physics mostly use a gray two-moment scheme to evolve neutrinos. This evolves integrals of $f_\nu$ over momentum space, preventing direct calculations of blocking factors. Monte Carlo methods may be an attractive alternative, providing access to the full distribution of neutrinos. Their current implementation is however inadequate to estimate $f_\nu$: in our most recent simulations, a single Monte Carlo packet causes, in the worst cases, estimates of $f_\nu$ to jump from $f_\nu=0$ to $f_\nu\sim 10^5$. While this is concerning, this brazen violation of the fermionic nature of neutrinos has been largely inconsequential as the interactions used in simulations avoid direct calculations of $f_\nu$. We are however reaching a level of modeling at which this problem can no longer be ignored. Here, we discuss the relatively simple origin of this issue. We then show that very rough estimates of $f_\nu$ can in theory be obtained in merger simulations, but that they will require a combination of unintuitive weighting schemes for Monte Carlo packets and smoothing of the neutrino distribution at coarser resolution than what the merger simulation uses.
Figures
Reference graph
Works this paper leans on
-
[1]
Y. Sekiguchi, An implementation of the microphysics in full general relativity: a general relativistic neutrino leak- age scheme, Class. Quantum Grav. 27, 114107 (2010), arXiv:1009.3358 [astro-ph.HE]
work page Pith review arXiv 2010
-
[2]
M. B. Deaton, M. D. Duez, F. Foucart, E. O’Connor, C. D. Ott, L. E. Kidder, C. D. Muhlberger, M. A. Scheel, and B. Szil´ agyi, Black Hole-Neutron Star Mergers with a Hot Nuclear Equation of State: Outflow and Neutrino- Cooled Disk for a Low-Mass, High-Spin Case, Astrophys. J. 776, 47 (2013), arXiv:1304.3384 [astro-ph.HE]. 11
arXiv 2013
- [3]
-
[4]
J. Barnes and D. Kasen, Effect of a High Opacity on the Light Curves of Radioactively Powered Transients from Compact Object Mergers, Astrophys. J. 775, 18 (2013), arXiv:1303.5787 [astro-ph.HE]
arXiv 2013
-
[5]
Foucart, Neutrino transport in general relativistic neu- tron star merger simulations, Liv
F. Foucart, Neutrino transport in general relativistic neu- tron star merger simulations, Liv. Rev. Comput. Astro- phys. 9, 1 (2023), arXiv:2209.02538 [astro-ph.HE]
arXiv 2023
-
[6]
M. Shibata, K. Kiuchi, Y. Sekiguchi, and Y. Suwa, Trun- cated Moment Formalism for Radiation Hydrodynamics in Numerical Relativity, Progress of Theoretical Physics 125, 1255 (2011), arXiv:1104.3937 [astro-ph.HE]
arXiv 2011
-
[7]
F. Foucart, E. O’Connor, L. Roberts, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Impact of an im- proved neutrino energy estimate on outflows in neutron star merger simulations, Phys. Rev. D94, 123016 (2016), arXiv:1607.07450 [astro-ph.HE]
arXiv 2016
- [8]
Show all 25 references
-
[9]
Foucart, M
F. Foucart, M. D. Duez, F. Hebert, L. E. Kidder, P. Ko- varik, H. P. Pfeiffer, and M. A. Scheel, Implementation of Monte Carlo Transport in the General Relativistic SpEC Code, Astrophys. J. 920, 82 (2021), arXiv:2103.16588 [astro-ph.HE]
2021 arXiv
-
[10]
J. M. Miller, B. R. Ryan, J. C. Dolence, A. Burrows, C. J. Fontes, C. L. Fryer, O. Korobkin, J. Lippuner, M. R. Mumpower, and R. T. Wollaeger, Full Transport Model of GW170817-Like Disk Produces a Blue Kilo- nova, Phys. Rev. D100, 023008 (2019), arXiv:1905.07477 [astro-ph.HE]
2019 arXiv
-
[11]
Kawaguchi, S
K. Kawaguchi, S. Fujibayashi, and M. Shibata, Long- term Monte Carlo-based neutrino-radiation hydrody- namics simulations for a black hole-torus system, Phys. Rev. D 111, 023015 (2025), arXiv:2410.02380 [astro- ph.HE]
2025 arXiv
-
[12]
Foucart, P
F. Foucart, P. C.-K. Cheong, M. D. Duez, L. E. Kid- der, H. P. Pfeiffer, and M. A. Scheel, Robustness of neu- tron star merger simulations to changes in neutrino trans- port and neutrino-matter interactions, Phys. Rev. D110, 083028 (2024), arXiv:2407.15989 [astro-ph.HE]
2024 arXiv
-
[13]
Y. Qiu, D. Radice, S. Richers, and M. Bhattacharyya, Neutrino Flavor Transformation in Neutron Star Merg- ers, arXiv (2025), arXiv:2503.11758 [astro-ph.HE]
2025 arXiv
-
[14]
H. H.-Y. Ng, C. Musolino, S. D. Tootle, and L. Rez- zolla, Accurate muonic interactions in neutron-star merg- ers and impact on heavy-element nucleosynthesis, arXiv (2024), arXiv:2411.19178 [astro-ph.HE]
2024 arXiv
-
[15]
Fujibayashi, K
S. Fujibayashi, K. Kiuchi, N. Nishimura, Y. Sekiguchi, and M. Shibata, Mass Ejection from the Remnant of a Binary Neutron Star Merger: Viscous-Radiation Hydrodynamics Study, Astrophys. J. 860, 64 (2018), arXiv:1711.02093 [astro-ph.HE]
2018 arXiv
-
[16]
P. C.-K. Cheong, F. Foucart, M. D. Duez, A. Of- fermans, N. Muhammed, and P. Chawhan, Energy- dependent and Energy-integrated Two-moment General- relativistic Neutrino Transport Simulations of a Hyper- massive Neutron Star, Astrophys. J. 975, 116 (2024), arXiv:2407.16017 [astro-ph.HE]
2024 arXiv
-
[17]
Foucart, Monte Carlo closure for moment-based trans- port schemes in general relativistic radiation hydrody- namic simulations, Mon
F. Foucart, Monte Carlo closure for moment-based trans- port schemes in general relativistic radiation hydrody- namic simulations, Mon. Not. Roy. Astr. Soc. 475, 4186 (2018), arXiv:1708.08452 [astro-ph.HE]
2018 arXiv
-
[18]
M. R. Izquierdo, J. F. Abalos, and C. Palenzuela, Guided moments formalism: A new efficient full-neutrino treat- ment for astrophysical simulations, Phys. Rev. D 109, 043044 (2024), arXiv:2312.09275 [astro-ph.HE]
2024 arXiv
-
[19]
Burrows, S
A. Burrows, S. Reddy, and T. A. Thompson, Neutrino opacities in nuclear matter, Nuc. Phys. A 777, 356 (2006), astro-ph/0404432
2006 arXiv
-
[20]
Banerjee, A
A. Banerjee, A. Dighe, and G. Raffelt, Linearized flavor- stability analysis of dense neutrino streams, Phys. Rev. D 84, 053013 (2011)
2011
-
[21]
M.-R. Wu, I. Tamborra, O. Just, and H.-T. Janka, Imprints of neutrino-pair flavor conversions on nucle- osynthesis in ejecta from neutron-star merger remnants, Phys. Rev. D96, 123015 (2017), arXiv:1711.00477 [astro- ph.HE]
2017 arXiv
-
[22]
Foucart, M
F. Foucart, M. D. Duez, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and E. Spira-Savett, General rel- ativistic simulations of collapsing binary neutron star mergers with monte carlo neutrino transport, Phys. Rev. D 107, 103055 (2023)
2023
-
[23]
A. W. Steiner, M. Hempel, and T. Fischer, Core-collapse Supernova Equations of State Based on Neutron Star Ob- servations, Astrophys. J.774, 17 (2013), arXiv:1207.2184 [astro-ph.SR]
2013 arXiv
-
[24]
O’Connor, An Open-source Neutrino Radiation Hy- drodynamics Code for Core-collapse Supernovae, Astro- phys
E. O’Connor, An Open-source Neutrino Radiation Hy- drodynamics Code for Core-collapse Supernovae, Astro- phys. J. Suppl. Ser. 219, 24 (2015), arXiv:1411.7058 [astro-ph.HE]
2015 arXiv
-
[25]
Foucart, M
F. Foucart, M. D. Duez, F. Hebert, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Monte-Carlo neutrino trans- port in neutron star merger simulations, Astrophys. J. Lett. 902, L27 (2020), arXiv:2008.08089 [astro-ph.HE]
2020 arXiv
Reviewed August 16, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.