Neural Boltzmann Equations

Jonas Spinner; Jack D. Shergold

hep-phAug 24, 2026

Neural Boltzmann Equations

Authors: Jonas Spinner, Jack D. Shergold

Organizations: Institute for Particle Physics Phenomenology Durham University, Durham, UK · Department of Mathematical Sciences University of Liverpool, Liverpool, UK

Abstract

The dynamics of particles in the early universe are described by Boltzmann equations, which involve high-dimensional phase-space integrals. Classical approaches use quadrature integration and evolve the system on a fixed momentum grid, which scales poorly to complicated systems and parameter scans, severely limiting the complexity of processes that can be studied. We introduce Neural Boltzmann Equations (NBEs), which combine three coupled concepts to overcome these limitations. First, particle properties are encoded in physics-inspired neural distribution functions, with parameters that can be predicted using neural networks, enabling efficient parameter scans. Second, phase-space integrals are evaluated with Monte Carlo, using importance sampling tools from collider physics. Third, we use the natural gradient method to evolve the system. After demonstrating the individual benefits of NBEs, we use the framework to perform the most numerically precise calculation to date of the effective number of relativistic neutrino degrees of freedom, $N_\mathrm{eff}$ , in the Standard Model.

Figures & tables

Figure 1: Natural gradient training for neural distribution functions. Left: Learned distribution functions for massless and massive bosons and fermions, compared to equilibrium targets. Right: Average residual over training time for massless bosons and fermions, comparing the natural gradient method with the Adam optimiser using the mean and standard deviation of ten independent trainings.

Figure 2: Efficient collision operator evaluation and coupling-conditioned training. Left: Collision operator for

\phi\phi\leftrightarrow\phi\phi

, calculated semi-analytically and with Monte Carlo using BEST and NBE, with the momentum-averaged relative Monte Carlo uncertainty,

\bar{\sigma}

, as precision metric. Right: Number density,

n_{t}

, and energy density,

\rho_{t}

, over time for different couplings, see Yoon (2026) for the BEST result. The lower panel rescales

t\to\lambda^{2}t

to collapse the curves for different couplings.

Figure 3: Neutrino decoupling and new physics impact. Left: Evolution of electromagnetic and neutrino energy densities during neutrino decoupling. Middle: Neutrino spectral distortions after decoupling, comparing the equilibrium solution to the full solution including spectral distortions. Right:

N_{\mathrm{eff}}

for the new physics model of Escudero et al. (2026b) , from a single NBE training.

Appendix figures & tables10 assets

Supplementary material from the paper’s appendix.

Appendix

Description	NBE result	$\Delta N_{\mathrm{eff}}$
Instantaneous decoupling	$3.00000$	–
+ $\nu$ - $e^{\pm}$ interactions	$3.03414$	$+0.03414$
+ spectral distortions	$3.03407$	$-0.00007$
+ finite-temperature QED	$3.04336$	$+0.00929$
+ flavour oscillations	$3.04392$	$+0.00056$

Appendix

Table 1: Decomposition of

N_{\mathrm{eff}}

contributions, as discussed in the text.

Quantity	Value	Reference
$\alpha$	$1/137.035999084$	Tiesinga et al. (2021)
$m_{e}^{\mathrm{phys}}$	$0.51099895\,\mathrm{MeV}$	Tiesinga et al. (2021)
$G_{F}^{\mathrm{phys}}$	$1.1663787\cdot 10^{-5}\,\mathrm{GeV}^{-2}$	Zyla and others (2020)
$m_{W}^{\mathrm{phys}}$	$80.379\,\mathrm{GeV}$	Zyla and others (2020)
$m_{\mathrm{Pl}}^{\mathrm{phys}}$	$1.22089\cdot 10^{19}\,\mathrm{GeV}$	Zyla and others (2020)
$g_{L}^{e}$	$0.727$	Kumar et al. (2013) ; Erler and Su (2013)

Appendix

Table 2: Experimental inputs used in the

N_{\mathrm{eff}}

precision calculation.

Effect	$10^{6}\cdot\delta N_{\mathrm{eff}}$
Seed statistics	$\pm 0.3$
Monte Carlo sampling (batch size)	$\pm 0.6$
Monte Carlo sampling (collision integral)	$\pm 0.6$
Moment quadrature	$\pm 0.7$
Evolution starting point	$\pm 0.6$
Step size controller	$\pm 1.4$

Appendix

Table 3: Contributions to the

\pm 2\cdot 10^{-6}

solver uncertainty in the

N_{\mathrm{eff}}

calculation in Eq. ( 17 ). We quantify the effects by increasing the precision for every parameter individually, and combine the resulting uncertainties as independent Gaussian effects.

Effect	$10^{6}\cdot\delta N_{\mathrm{eff}}$
$\sin^{2}\theta_{12}$ ( $3\sigma$ band)	$\pm 18$
$\sin^{2}\theta_{13}$ ( $3\sigma$ band)	$\pm 4$
$\Delta m^{2}_{21}$ ( $3\sigma$ band)	$\pm 3$
Neutrino mass ordering (NO $\to$ IO)	$+7$
Total	$\pm 20$

Appendix

Table 4: Contributions to the

\pm 2\cdot 10^{-5}

experimental uncertainty in the

N_{\mathrm{eff}}

calculation in Eq. ( 17 ). We quantify the effects by varying the input parameters within their reported experimental uncertainty band, and for the neutrino mass ordering we report the effect of flipping the choice.

Figure 4: Convergence speed of the gamma mixture series for massless bosons. We compare the exact Bose-Einstein distribution with the shapes obtained after truncating the series after

M

terms.

Panel of Figure 5 — Figure 5: Trained neural distribution functions with different numbers of mixture components $M$ , for massless and massive fermions and bosons.

Figure 6: Left: Longer Adam training for massless equilibrium distribution functions. Right: Benefit of neural importance sampling for the thermalisation experiment in Section 5.2 , comparing NBE with and without the neural importance sampling refinement using MadNIS.

Figure 7: Effect of the moment conservation trick of Appendix A.3 for scalar thermalisation, including only

2\leftrightarrow 2

processes (left) and also

2\leftrightarrow 3

processes (right).

Setting	precision	fast
Batch size	$4000\to 8000$	$500\to 1000$
Collision integral samples	$4000\to 8000$	$4000\to 8000$
Neutrino moment quadrature	$32\to 64$	same
EM-sector moment quadrature	$64\to 128$	same
Photon series truncation	$1000\to 2000$	same
QED two-dimensional quadrature	$32\to 64$	same

Appendix

Table 5: Solver settings varied for the uncertainty estimate of Table 6 . We write

a\to b

for each entry, with

a

the default value and

b

the increased-precision value used for the solver uncertainty estimate.

Effect	$10^{6}\cdot\delta N_{\mathrm{eff}}$ (precision)	$10^{6}\cdot\delta N_{\mathrm{eff}}$ (fast)
Seed statistics	$\pm 0.3$	$\pm 1.1$
Monte Carlo sampling (batch size)	$\pm 0.6$	$\pm 3.4$
Monte Carlo sampling (collision integral)	$\pm 0.6$	$\pm 2.4$
Moment quadrature	$\pm 0.7$	$\pm 4.9$
Evolution starting point	$\pm 0.6$	$\pm 2.4$
Step size controller	$\pm 1.4$	$\pm 6.5$

Appendix

Table 6: Contributions to the solver uncertainty in the

N_{\mathrm{eff}}

calculation, for the two solver settings of Table 5 . We quantify the effects by increasing the precision for every parameter individually and combine the resulting uncertainties as independent Gaussian effects.

Neural PDE Solvers

Neural PDE solvers use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.

See all 269 papers

Weak Adversarial Neural Pushforward Method for Boltzmann Equation
Aug 7, 2026Jenia Fardousi Koly, Andrew Qing He, Wei CaiNeural PDE Solvers
Department of Mathematics, Southern Methodist University, Dallas, Texas, USA
Neural Transport Nested Sampling
Sep 24, 2026David Yallup, Will HandleyMarkov Chain Monte Carlo Boltzmann Sampling
Kavli Institute for Cosmology Cambridge Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
A Neural Score-Based Particle Method for the Vlasov-Maxwell-Landau System
Mar 26, 2026Vasily Ilin, Jingwei HuNeural PDE Solvers
Department of Mathematics, University of Washington, Seattle, WA, USA. · Department of Applied Mathematics, University of Washington, Seattle, WA, USA.

Abstract

Neural PDE Solvers

Explore similar work

Weak Adversarial Neural Pushforward Method for Boltzmann Equation

Neural Transport Nested Sampling

A Neural Score-Based Particle Method for the Vlasov-Maxwell-Landau System