Organizations: Institute for Particle Physics Phenomenology Durham University, Durham, UK · Department of Mathematical Sciences University of Liverpool, Liverpool, UK
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, Neff, 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 ϕϕ↔ϕϕ , calculated semi-analytically and with Monte Carlo using BEST and NBE, with the momentum-averaged relative Monte Carlo uncertainty, σˉ , as precision metric. Right: Number density, nt , and energy density, ρt , over time for different couplings, see Yoon (2026) for the BEST result. The lower panel rescales t→λ2t 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: Neff 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
ΔNeff
Instantaneous decoupling
3.00000
–
+ ν - e± 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 Neff contributions, as discussed in the text.
Quantity
Value
Reference
α
1/137.035999084
Tiesinga et al. (2021)
mephys
0.51099895MeV
Tiesinga et al. (2021)
GFphys
1.1663787⋅10−5GeV−2
Zyla and others (2020)
mWphys
80.379GeV
Zyla and others (2020)
mPlphys
1.22089⋅1019GeV
Zyla and others (2020)
gLe
0.727
Kumar et al. (2013) ; Erler and Su (2013)
Appendix
Table 2: Experimental inputs used in the Neff precision calculation.
Effect
106⋅δNeff
Seed statistics
±0.3
Monte Carlo sampling (batch size)
±0.6
Monte Carlo sampling (collision integral)
±0.6
Moment quadrature
±0.7
Evolution starting point
±0.6
Step size controller
±1.4
Appendix
Table 3: Contributions to the ±2⋅10−6 solver uncertainty in the Neff 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
106⋅δNeff
sin2θ12 ( 3σ band)
±18
sin2θ13 ( 3σ band)
±4
Δm212 ( 3σ band)
±3
Neutrino mass ordering (NO → IO)
+7
Total
±20
Appendix
Table 4: Contributions to the ±2⋅10−5 experimental uncertainty in the Neff 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.
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↔2 processes (left) and also 2↔3 processes (right).
Setting
precision
fast
Batch size
4000→8000
500→1000
Collision integral samples
4000→8000
4000→8000
Neutrino moment quadrature
32→64
same
EM-sector moment quadrature
64→128
same
Photon series truncation
1000→2000
same
QED two-dimensional quadrature
32→64
same
Appendix
Table 5: Solver settings varied for the uncertainty estimate of Table 6 . We write a→b for each entry, with a the default value and b the increased-precision value used for the solver uncertainty estimate.
Effect
106⋅δNeff (precision)
106⋅δNeff (fast)
Seed statistics
±0.3
±1.1
Monte Carlo sampling (batch size)
±0.6
±3.4
Monte Carlo sampling (collision integral)
±0.6
±2.4
Moment quadrature
±0.7
±4.9
Evolution starting point
±0.6
±2.4
Step size controller
±1.4
±6.5
Appendix
Table 6: Contributions to the solver uncertainty in the Neff 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 use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.
In this paper, we extend a weak adversary neural network pushforward method for solving time dependent Boltzmann equation and a weak formulation of the collision operator is proposed where an invertible neural pushforward mapping is used to generating samples given by the distribution governed by the Boltzmann equation. The training of the pushforward mapping is learnt by enforcing the weak form of the Boltzmann equation. Numerical results have demonstrated the effectiveness of the proposed method.
Jenia Fardousi Koly, Andrew Qing He, Wei Cai
Department of Mathematics, Southern Methodist University, Dallas, Texas, USA
Sampling from Boltzmann distributions of molecular systems is an inference problem that has seen significant recent developments fuelled by advances in neural density estimation. We develop a novel sampling algorithm, Neural Transport Nested Sampling (NTNS), which combines the classical strengths of nested sampling with modern neural flow-based methods. NTNS uses a flow matching velocity as the drift in a Metropolis--Hastings corrected Langevin kernel inside a nested sampling outer loop, requiring only evaluations of the target energy function and providing scalable estimation of the full partition function of high-dimensional particle systems. We benchmark NTNS on challenging molecular sampling benchmarks, scaling up to Lennard--Jones clusters of 55 interacting particles, where it reduces both interatomic distance and energy Wasserstein errors to reference MCMC by over an order of magnitude relative to the strongest neural baselines at lower wall-clock cost. To our knowledge, NTNS is also the first neural sampler to return a calibrated, temperature resolved partition function estimate at this scale, recovering the phase structure across temperature from a single run.
David Yallup, Will Handley
Kavli Institute for Cosmology Cambridge Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
Plasma modeling is central to the design of nuclear fusion reactors, yet simulating collisional plasma kinetics from first principles remains a formidable computational challenge: the Vlasov-Maxwell-Landau (VML) system describes six-dimensional phase-space transport under self-consistent electromagnetic fields together with the nonlinear, nonlocal Landau collision operator. A recent deterministic particle method for the full VML system estimates the velocity score function via the blob method, a kernel-based approximation with O(n2) cost. In this work, we replace the blob score estimator with score-based transport modeling (SBTM), in which a neural network is trained on-the-fly via implicit score matching at O(n) cost. We prove that the approximated collision operator preserves momentum and kinetic energy, and dissipates an estimated entropy. We also characterize the unique global steady state of the VML system and its electrostatic reduction, providing the ground truth for numerical validation. On three canonical benchmarks -- Landau damping, two-stream instability, and Weibel instability -- SBTM is more accurate than the blob method, achieves correct long-time relaxation to Maxwellian equilibrium where the blob method fails, and delivers 50% faster runtime with 4× lower peak memory.
Vasily Ilin, Jingwei Hu
Department of Mathematics, University of Washington, Seattle, WA, USA. · Department of Applied Mathematics, University of Washington, Seattle, WA, USA.