⚛️ phenomenology

Neutrino Direct Simulation Monte Carlo: Accurate modeling of out-of-equilibrium decaying cosmological relics

This paper extends the Neutrino Direct Simulation Monte Carlo ( $\nu$ DSMC) framework to jointly model decaying cosmological relics, neutrino transport, and Big Bang nucleosynthesis, demonstrating that this efficient and comprehensive approach outperforms traditional approximations by accurately capturing nonthermal effects that significantly impact cosmological predictions.

Original authors: Kensuke Akita, Miguel Escudero, Oleksii Ihnatenko, Maksym Ovchynnikov

Published 2026-09-29

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Kensuke Akita, Miguel Escudero, Oleksii Ihnatenko, Maksym Ovchynnikov

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the first few seconds after the Big Bang, the universe was a seething, dense soup of fundamental particles. Among the most abundant of these were neutrinos, ghostly particles that rarely interact with anything else. As the universe expanded and cooled, these neutrinos stopped bumping into the surrounding plasma of electrons and photons, effectively "decoupling" and drifting freely through space. This moment, occurring roughly one second after the beginning of time, set the stage for the formation of the first atomic nuclei in a process known as Big Bang nucleosynthesis. The way these neutrinos behaved during that critical window determined how much energy was available to drive the expansion of the universe and how many neutrons and protons were available to build the first atoms. Today, astronomers can look back at this ancient history by measuring the abundance of light elements like helium and deuterium, as well as by studying the afterglow of the Big Bang. These observations act as a precise fossil record, allowing scientists to test whether the standard story of the early universe holds up or if something unexpected happened.

For decades, the standard model of cosmology has assumed that the universe evolved in a predictable, smooth way during those first few seconds. However, recent theoretical work has suggested that the universe might have been more chaotic. Imagine a heavy, unstable particle that existed in the early universe and decayed later than expected, injecting a fresh burst of energy into the mix. If this happened while neutrinos were still trying to settle into their natural state, it could have pushed them out of equilibrium, creating a jumbled, non-standard distribution of speeds and energies. Such a scenario would leave a distinct fingerprint on the cosmic radiation density and the amounts of light elements created. The problem is that calculating exactly how these neutrinos would behave under such extreme, non-standard conditions is incredibly difficult. Traditional methods often rely on simplifying assumptions, such as averaging out the different speeds of the particles, which can lead to qualitatively wrong answers when the physics gets complicated.

A team of researchers has now developed a new, highly detailed computational tool to solve this problem, allowing them to simulate the chaotic dance of neutrinos with unprecedented accuracy. They call their method "Neutrino Direct Simulation Monte Carlo." Instead of trying to solve complex equations for the average behavior of the entire neutrino population, this approach tracks millions of individual, simulated neutrinos as they zip through the expanding universe. The computer follows each particle, deciding when and how it interacts with the surrounding plasma, when it changes its identity through a process called flavor oscillation, and how it responds to the decay of those hypothetical heavy relics. By following the actual trajectories and collisions of these simulated particles, the researchers can build up a complete picture of the neutrino energy distribution without relying on the simplifying averages that often fail in extreme scenarios.

The researchers tested this new framework against several different cosmological scenarios, including cases where heavy particles decayed into neutrinos or into electromagnetic energy. They compared their results with other sophisticated methods that use different mathematical approaches, such as solving quantum kinetic equations. In cases where the physics was relatively calm and the neutrinos remained close to a standard thermal state, the new simulation agreed perfectly with the established methods, confirming its reliability. However, the true power of the tool emerged in the more extreme scenarios. When the injected energy was high enough to create a highly distorted, non-thermal population of neutrinos, the traditional methods that rely on averaging began to fail. In some cases, these older approaches predicted the wrong sign for the change in the number of effective neutrino species, essentially getting the direction of the effect backwards. They also produced significantly incorrect predictions for the abundance of primordial elements like helium and deuterium.

The study found that the new simulation method is not only more accurate in these difficult regimes but also substantially faster than previous high-precision techniques when dealing with energetic, non-thermal neutrino populations. This speed and accuracy make it practical to run broad scans across a wide range of possible new physics models, something that was previously too computationally expensive to do reliably. The researchers also incorporated the production of charged pions—short-lived particles that can further alter the balance between neutrons and protons—into their simulation, ensuring that the link between the neutrino evolution and the final nuclear abundances was complete. By coupling their transport simulation directly to the calculations for Big Bang nucleosynthesis, they could trace the entire chain of events from the decay of a heavy relic to the final composition of the universe.

The results of this work suggest that when the early universe deviates significantly from thermal equilibrium, the details of the neutrino momentum distribution matter immensely. Approximations that smooth over these details can lead to misleading conclusions about the nature of new physics. The new simulation framework provides a robust way to test these ideas, offering a reliable path to understanding how neutrinos thermalize and what observable consequences that thermalization has for the cosmos. By following the individual histories of these elusive particles, the researchers have opened the door to more rigorous tests of the standard cosmological model and the potential discovery of new physics lurking in the first moments of our universe.

Technical Summary: Neutrino Direct Simulation Monte Carlo (νDSMC)

Problem Statement
Non-standard cosmologies involving decaying relics (unstable particles $X$ ) that inject energy into the universe at MeV temperatures can drive neutrinos out of thermal equilibrium. This departure from equilibrium alters the cosmic radiation density ( $N_{\text{eff}}$ ) and the primordial light element abundances (Big Bang Nucleosynthesis, or BBN). Existing methods for modeling this evolution face significant limitations:

Momentum-averaged approaches (e.g., thermal Boltzmann equations) assume neutrino spectra remain close to Fermi-Dirac distributions. These approximations can yield qualitatively incorrect predictions, including the wrong sign for changes in $N_{\text{eff}}$ , when spectral distortions are large.
Momentum-resolved approaches (e.g., Quantum Kinetic Equations or Quasi-classical Boltzmann equations) are computationally expensive, particularly when energetic neutrinos are injected, and often struggle to incorporate complex processes like inverse decays, quantum statistics, and hadronic production simultaneously.
Previous Monte Carlo implementations (early versions of $\nu$ DSMC) lacked explicit evolution of unstable relics, did not include inverse decays or charged-pion production, and suffered from algorithmic errors related to the ordering of expansion and interaction operators and timestep-dependent flavor updates.

Methodology
The authors present an extended Neutrino Direct Simulation Monte Carlo (νDSMC) framework that jointly evolves the relic population, neutrino spectra, electromagnetic plasma, and cosmic expansion. The method replaces the direct evaluation of collision kernels with Monte Carlo sampling of particle interactions. Key technical extensions include:

Strang Splitting: The algorithm replaces the asymmetric ordering of expansion and interaction operators ( $O_{\text{exp}} \to O_{\text{int}}$ ) with a symmetric Strang splitting scheme ( $O_{\text{exp}}(\Delta t/2) \to O_{\text{int}} \to O_{\text{exp}}(\Delta t/2)$ ). This reduces the splitting error from $O(\Delta t)$ to $O((\Delta t)^2)$ .
Majorant Collision Frequency (MCF): The previous No-Time-Counter (NTC) scheme is replaced by an MCF scheme. This allows for continuous-time event sampling, enabling the simulation of decays and scatterings on a single event clock without artificial splitting between decay and scattering operators. It supports large timesteps where $\Gamma_{\text{max}}\Delta t \gg 1$ .
Relic and Inverse Decay Evolution: The framework explicitly evolves the phase-space distribution of the unstable relic $X$ . It includes both direct decays ( $X \to \text{SM}$ ) and inverse decays ( $\text{SM} \to X$ ), accounting for thermal masses and quantum-statistical factors (Pauli blocking and Bose stimulation) in the plasma.
Flavor Conversion: Neutrino oscillations are treated using adiabatic propagation in the propagation basis with averaged oscillation phases (ATAO). Flavor labels are sampled based on in-medium projectors at the moment of interaction, rather than globally reassigning flavors at fixed timesteps.
Quantum Statistics and Pion Production: The method calculates local occupation numbers in momentum bins to determine Pauli blocking and Bose enhancement factors dynamically. It also includes a deterministic source counter for charged-pion production ( $\nu e \to \pi$ ) induced by non-thermal neutrinos, providing inputs for BBN calculations.
Zero-Mode Handling: For light relics where a significant fraction of the population may occupy the zero-momentum state, a specialized treatment handles the inverse decay into this state and the associated Bose stimulation.

Key Contributions and Results
The paper validates the extended framework through equilibrium tests and comparisons with independent state-of-the-art calculations across four cosmological scenarios:

Validation: The framework reproduces the standard cosmological prediction $N_{\text{eff}} \approx 3.044$ and passes tests for relativistic time dilation, cosmological expansion, and chemical equilibrium in fixed volumes.
Comparison with Other Methods:
- In electromagnetic reheating scenarios (Setup I), $\nu$ DSMC shows excellent agreement with Quantum Kinetic Equations (QKE) and momentum-resolved Boltzmann solvers (qcbe) for $N_{\text{eff}}$ and electron-neutrino spectra.
- In neutrinophilic decay scenarios (Setups II and III), the authors demonstrate that momentum-averaged "Thermal" approximations can fail qualitatively. When departures from equilibrium are large, averaging over momenta can predict the wrong sign for $\Delta N_{\text{eff}}$ and yield substantially incorrect primordial abundances.
- The extended $\nu$ DSMC is found to be substantially faster than momentum-resolved Boltzmann solvers for strongly non-thermal populations, making broad parameter scans practical.
Physical Insights:
- Inverse Decays: For light relics ( $m_X \lesssim T$ ), inverse decays can replenish the relic population, reshaping the history of neutrino injection and affecting thermalization.
- Pion Production: While the yield of charged pions from neutrino thermalization is small for the masses considered ( $m_X \lesssim 10$ GeV), the framework quantifies their impact on BBN, showing they can modify the neutron-to-proton ratio.
- Spectral Distortions: The study confirms that accurate predictions for $N_{\text{eff}}$ and BBN abundances in non-standard cosmologies require resolving the full momentum distribution of neutrinos, as thermal approximations are insufficient when spectral distortions are significant.

Significance and Claims
The paper claims that the extended $\nu$ DSMC framework enables reliable cosmological tests of new physics by following neutrino thermalization and its observable consequences with high fidelity. Its primary significance lies in:

Robustness: It includes a broader range of physical processes (direct/inverse decays, flavor conversion, quantum statistics, pion production) than previous Monte Carlo or Boltzmann approaches.
Efficiency: It offers a computational speed advantage over momentum-resolved kinetic solvers for non-thermal populations, facilitating comprehensive parameter scans.
Accuracy: It avoids the qualitative errors associated with momentum-averaged approximations in regimes of large non-equilibrium, ensuring correct predictions for $N_{\text{eff}}$ and light element abundances.

The authors conclude that this framework is a necessary tool for interpreting future precision cosmological data (such as CMB and BBN measurements) in the context of decaying relics and non-standard neutrino interactions. The implementation is released in Mathematica, with a C++ version planned for public release.

More like this