PaperlayerView exact revision on arXiv

Source-linked AI summary

Neural Boltzmann Equations

Jonas Spinner, Jack Shergold

arXiv:2608.23022v1hep-phcs.LGstat.ML

TL;DR

Conventional early-universe Boltzmann solvers struggle with high-dimensional phase-space integrals, complex processes, and parameter scans. The paper introduces NBEs, combining neural distribution functions, Monte Carlo importance sampling, and natural-gradient evolution. NBEs achieve 1000× more efficient importance sampling, per-mille precision for a complex 2 ↔3 system, agreement with literature for relativistic neutrino degrees of freedom, and a two-dimensional conditioned parameter scan.

  • Problem

    Conventional Boltzmann solvers use fixed momentum grids, quadrature, and process-specific reductions that become costly for complex systems and parameter scans.

  • Method

    NBEs combine physics-inspired neural distribution functions, Monte Carlo collision integrals with importance sampling, and natural-gradient evolution of the Boltzmann equation.

  • Results

    1000× more efficient importance sampling enabled per-mille precision for a complex non-equilibrium 2 ↔3 system; the framework also agreed with literature for relativistic neutrino degrees of freedom and supported a two-dimensional conditioned parameter scan.

  • Takeaways & Limitations

    NBEs provide a precise, flexible, and scalable framework for early-universe calculations involving higher-multiplicity processes and high-dimensional parameter spaces.

  • Takeaways & Limitations

    Natural-gradient optimisation is limited mainly by the overhead of matrix inversion, while collision-operator evaluation dominates computational cost in the reported experiments.

Abstract

from arXiv · show

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 a precision calculation of the effective number of relativistic neutrino degrees of freedom in the early universe.

1 INTRODUCTION

Neural Boltzmann Equations (NBEs) are introduced to make early-universe Boltzmann calculations more flexible and scalable. They combine neural distribution functions, Monte Carlo collision integrals, and natural-gradient training, supporting efficient parameter scans and complex-process calculations.

  • NBEs target early-universe Boltzmann equations, whose conventional solutions rely on momentum discretisation, distribution assumptions, and process-specific interaction-rate reductions.
  • Neural distribution functions provide a discretisation-free, physics-inspired parameterisation that can be sampled efficiently.Their parameters can be predicted from arbitrary conditioning variables, enabling parameter scans with a single distribution-function solve.
  • Monte Carlo collision integrals use importance-sampling techniques developed for collider event generators to evaluate high-dimensional phase-space integrals.
  • Natural-gradient training efficiently trains neural distribution functions on the Boltzmann equation.
  • NBEs are applied to a precision calculation of relativistic neutrino degrees of freedom and to a full parameter scan for a new-physics model after one training.

2 BACKGROUND

The background develops the physical and numerical structure underlying NBEs, from particle distributions and collision operators to phase-space integration and natural-gradient optimisation. It also identifies matrix inversion as the main limitation of the natural-gradient method.

  • Early universe dynamics: The Boltzmann equation couples particle propagation in an expanding universe to microscopic interactions through distribution functions and collision operators.The Liouville operator describes propagation, while the collision operator encodes interactions involving all species.
  • Early universe dynamics: FLRW homogeneity and isotropy reduce each distribution function’s dependence to one momentum variable, while comoving variables concentrate relevant evolution near x = 1.
  • Collision operators: Collision operators combine statistical factors, squared matrix elements, and phase-space integration to describe processes that populate or deplete particle phase space.Statistical factors include Pauli blocking for fermions and Bose enhancement for bosons.
  • Observables: Macroscopic number, energy, and pressure densities are moments of distribution functions, and Boltzmann-equation moments yield fluid equations constraining those quantities.
  • Phase-space integration: Beyond two or three final-state particles, Monte Carlo integration is the only viable option because phase-space integrals are high-dimensional, expensive, and sharply featured.Importance sampling uses proposal distributions over mapped unit-hypercube variables; collider techniques refine these proposals with process-specific mappings and generative networks.
  • Natural-gradient optimisation: Natural-gradient optimisation uses function-space geometry through the Fisher matrix, with direct or Gram-matrix inversion chosen according to network size.The overhead from matrix inversion is the main limitation of the method.

3 RELATED WORK

Related work describes conventional Boltzmann solvers as accurate but increasingly specialised and costly for complex processes, parameter variation, and growing phase-space dimensionality. Neural quantum states provide a related precedent for neural ansätze trained directly from governing equations.

  • Boltzmann equation solvers: Existing Boltzmann solvers are tailored to specific problems through analytic reductions and numerical approximations, such as fixed-grid collision reductions for neutrino decoupling.
  • Boltzmann equation solvers: Momentum-grid and moment-based solvers lose phase-space information or become rapidly intractable as particle, flavour, or momentum-point counts grow.
  • Boltzmann equation solvers: Standard quadrature makes collision integrals prohibitively expensive beyond simple 2 →2 scattering, while each parameter point generally requires a complete system solve.
  • Boltzmann equation solvers: Process-specific analytic reductions and quadrature grids require substantial re-derivation when particles or processes are added.
  • Neural quantum states: Neural quantum states use neural networks as variational ansätze trained directly on governing equations rather than data.

4 METHODS

Neural Boltzmann Equations combine neural distribution functions, Monte Carlo phase-space integration, and natural-gradient evolution to solve early-universe Boltzmann equations flexibly and efficiently.

  • NBEs combine physics-inspired neural distribution functions, Monte Carlo collision integrals, and sequential natural-gradient updates.The framework targets flexible solutions without extensive process-specific assumptions.
  • Neural distribution functions encode equilibrium solutions, avoid momentum discretisation, support efficient sampling, and can be predicted from conditioning variables.A neural network can predict their parameters, allowing cached reuse across evaluations and parameter scans.
  • The distribution parameterisation factorises particle distributions into an invariant normalisation and a normalised momentum density for importance sampling.Mixtures of gamma distributions parameterise the density; fermionic constructions enforce f(y) ≤1.
  • Monte Carlo collision integrals sample initial momenta from neural densities and generate valid final states with MadSpace mappings and importance weights.Normalising flows can refine initial- and final-state sampling, while the generated momenta evaluate matrix elements and gain-loss factors.
  • The solver trains distributions through fixed-point iterations, evaluating collision integrals and applying natural-gradient steps until convergence or a stopping criterion.Previous-step parameters provide a warm start, and conditions are evaluated in batches with cached distribution parameters.
  • The solver’s update rule is invariant under small perturbations of the new solution, supporting application to stiff regimes near thermal equilibrium.Collision-operator evaluation remains the dominant computational cost, typically requiring 10^3 function evaluations per test momentum.

5 RESULTS

NBEs accurately learn equilibrium distributions, improve collision-integral sampling and natural-gradient training efficiency, and support coupling-conditioned evolution and precision neutrino calculations with parameter scans.

  • Detailed balance: 10^-6 relative precision is achieved in the bulk when learning massless and massive Bose-Einstein and Fermi-Dirac equilibrium distributions.The massless Fermi-Dirac distribution reaches even higher precision because one mixture expresses it exactly.
  • Detailed balance: 10^-6 average residual is reached in 10 natural-gradient iterations, compared with 10^-2 for Adam under the same budget.Adam plateaus at 10^-4 after around 10k iterations, while natural-gradient updates are more than an order of magnitude more efficient overall.
  • Thermalisation: 20× lower average Monte Carlo uncertainty than BEST corresponds to a 400× more sample-efficient collision-operator estimate for ϕϕ ↔ϕϕ.Distribution-function sampling, MadSpace mappings, and MadNIS refinement provide a combined efficiency gain of 1000×.
  • Thermalisation: For ϕϕ ↔ϕϕϕ with 10^6 samples, the average uncertainty is 6.8% for BEST and 0.58% for NBE.The comparison uses the same uncertainty hierarchy as the ϕϕ ↔ϕϕ collision operator.
  • Thermalisation: Coupling-conditioned NBE evolution preserves energy density within 0.1% and produces number-density evolution consistent across coupling strengths within 0.1%.After time rescaling, the number density converges to nt=∞/nt=0 ≈0.900 at equilibrium.
  • Neutrino decoupling: The NBE framework reproduces the precision neutrino calculation and identifies mϕ ∈[8 MeV, 13 MeV] where the new-physics model agrees with Neff = 2.81 ± 0.12.The numerical uncertainty varies between 10^-3 and 10^-4 depending on the parameter region.

6 CONCLUSION

Neural Boltzmann Equations provide a precise, flexible framework for early-universe Boltzmann equations without momentum discretisation. Their Monte Carlo integration, neural distribution functions, and natural-gradient training support complex processes and parameter scans.

  • NBEs solve early-universe Boltzmann equations without momentum discretisation, using neural distribution functions conditioned on theory parameters.Natural-gradient training supports both conditional and unconditional distributions to high precision.
  • 1000× more efficient importance sampling than a recently proposed Monte Carlo Boltzmann solver enabled per-mille evolution of a system with number-changing 2 ↔3 processes.
  • The framework produced a precision calculation of relativistic neutrino degrees of freedom with excellent agreement with the literature.
  • A single conditioned NBE performed a two-dimensional parameter scan of an extended Neff scenario.
  • Together, the results establish NBEs as a precise, flexible, and scalable framework for higher-multiplicity processes and high-dimensional parameter spaces.

A METHOD DETAILS

The method details define neural distribution functions and their sampling strategy for equilibrium and learned distributions. The construction combines gamma-mixture representations, rejection sampling, quadrature, phase-space mappings, and adaptive Monte Carlo allocation.

  • Equilibrium solutions: The neural distribution ansatz represents equilibrium distributions of arbitrary mass using kinetic energy ε and gamma-mixture components.The exact equilibrium density is represented as a mixture with a bounded correction factor, while learned distributions absorb that factor into mixture weights.
  • Equilibrium solutions: Bose-Einstein distributions are expanded in a geometric series, with truncation at M terms leaving relative error O(e−Mε).
  • Equilibrium solutions: In the non-relativistic limit, the Maxwell-Boltzmann distribution becomes a single gamma distribution, peq(y) dy = γ(ε|3/2, 1).
  • Distribution function sampling: For fermions, rejection sampling maps proposals from the Maxwell-Boltzmann distribution to Fermi-Dirac distributions; the massless zero-chemical-potential efficiency is N = π2/12 ≈0.82.
  • Moments of neural distribution functions: Moments of neural distribution functions reduce to expectations under the normalised density and can be evaluated with Gauss-Laguerre quadrature when needed.
  • Equilibrium solutions: The implementation switches from Laguerre to Hermite nodes for a = m/T > 0.1 after mapping a square-root branch point to a Gaussian half-line integral.
  • Distribution function sampling: Factorising f(y) into invariant density η and normalised density p(y) enables efficient importance sampling of collision-integral initial states.Sampling angles uniformly and momenta from p(y) causes phase-space factors to cancel, leaving η/f(y) to combine with statistical factors.
  • Collision integral evaluation: Neyman allocation dynamically assigns Monte Carlo samples across processes using their previous-step standard deviations rather than equal sample counts.

A.3 BOLTZMANN EQUATION TRAINING

Boltzmann-equation training combines multistep integrators, approximate collision-operator derivatives, and neutrino-decoupling equations under explicit cosmological assumptions. The formalism reduces oscillation-averaged quantum kinetic evolution to tractable distribution-function equations.

  • Integrator schemes: The training implementation supports Euler, BDF2, and Crank-Nicolson schemes, with Crank-Nicolson as the default choice using trapezoidal collision-integral refinement.
  • Collision operator derivative estimate: The collision-operator derivative approximation avoids the impractical cost of exact differentiation while improving on the unstable naive choice M = 1.The approximation is exact on the diagonal y = y′ after neglecting specified f1 occurrences.
  • Standard Model neutrino decoupling: Neutrino decoupling is modeled in FLRW spacetime with no lepton asymmetry, no CP-violating lepton-sector phase, massless neutrino collision kinematics, and truncated weak collision rates.
  • Quantum kinetic equation: Neutrino distributions are represented by a Hermitian density matrix whose diagonal occupations and off-diagonal components encode flavour populations and coherence.
  • Adiabatic transfer of averaged oscillations: Adiabaticity and rapid oscillations justify neglecting the second commutator and averaging off-diagonal phases, yielding a diagonal mass-basis QKE with Boltzmann-equation form.
  • Flavour-diagonal limit: In the flavour-diagonal limit, mixing changes Neff only indirectly through interaction-rate redistribution, with a reported correction ΔNeff = 0.0005.
  • Flavour-diagonal limit: Identical initial conditions reduce the three neutrino flavours to two independent distributions, fνe and fνβ with β ∈ {µ, τ}.
  • Equilibrium electromagnetic bath: Electrons and photons remain in equilibrium, so the electromagnetic sector is described by its temperature TEM without chemical potentials.

C EXPERIMENT DETAILS

The experiments use fixed neural-distribution hyperparameters, natural-gradient stabilization, Crank-Nicolson integration, adaptive step sizes, and specified iteration and optimizer settings.

  • Hyperparameters: All experiments use M = 50 gamma-mixture components, Nquad = 16 quadrature nodes, and a three-layer, 64-channel tanh multilayer perceptron.
  • Natural-gradient training: Natural-gradient matrix inversion uses diagonal Tikhonov damping of 10−8, increased to 10−5 for conditional-training Gram matrices.
  • Time integration: Crank-Nicolson refines the solver trajectory, while an uncertainty estimate drives a PI step-size controller with experiment-specific minimum and maximum step sizes.
  • Training convergence: Fixed-point iteration typically converges within 10 iterations for unconditional training and up to 50 iterations for conditional training.
  • Normalising-flow training: Normalising flows use 10 bins, a three-layer 32-channel MLP, and Adam with learning rate γ = 10−3.
  • Hyperparameters: Adam trainings use learning rate γ = 0.02 and batch size B = 1000; 100 iterations take around five seconds on a GPU.

C.2 THERMALISATION

The thermalisation study tests NBE evolution on a challenging massive-particle model with both 2 ↔2 and 2 ↔3 processes. The known asymptotic number density provides an analytic benchmark for assessing the evolution and coupling-conditional training.

  • Thermalisation: The model contains a massive particle with m = 1 and includes both 2 ↔2 and 2 ↔3 processes sharing coupling λ.The collision operator scales as λ^2, allowing the coupling to be absorbed into a rescaled time variable.
  • Thermalisation: The evolution in rescaled time t′ is independent of λ, providing a test of coupling-conditional training.The model's coupling dependence is used to assess the precision of conditional training in Figure 2.
  • Thermalisation: The non-thermal initial condition evolves toward the Bose-Einstein equilibrium configuration through momentum-mode coupling in the collision operator.Energy density is conserved, while the 2 ↔3 process can change the number density.
  • Thermalisation: n_t→∞ ≈ 0.900272 is the analytically determined asymptotic number density used as a benchmark in Figure 2.The limit follows from the known initial and final distributions together with energy conservation.
  • Thermalisation: The semi-analytical baseline is biased at y → 0, and refined quadrature reduces the resulting shift.The tilt toward larger NBE deviations at y = 10^-1 is attributed to the suboptimal quadrature result.
  • Thermalisation: A full training takes 4 hours on an A30 GPU under the stated thermalisation settings.The calculation uses batch size B = 1000, N_s = 10^4 Monte Carlo samples, and fixed step size Δt = 20.

C.3 NEUTRINO DECOUPLING

The neutrino-decoupling study applies NBE approximations to calculate Neff while modeling the electromagnetic plasma in equilibrium and neutrinos with general distributions. It also evaluates training costs and uncertainty, including a conservative method uncertainty.

  • Neutrino decoupling: The calculation follows Froustey et al. (2020) and targets 10^-4 precision in Neff.The electromagnetic plasma is assumed to remain in thermal equilibrium, while the three neutrino species use general distribution functions.
  • Neutrino decoupling: The neutrino-decoupling trainings range from 20 to 90 minutes on an A30 GPU, depending on the distribution ansatz and flavour-oscillation treatment.The general ansatz without flavour oscillations takes 50 minutes, and including flavour oscillations takes 90 minutes.
  • Neutrino decoupling: Conditional training takes 3 hours on the same GPU and does not include flavour oscillations.
  • Neutrino decoupling: Random seed variations are 10^-5 for unconstrained distributions and 5 · 10^-5 for equilibrium distributions under default settings.The authors suspect the larger equilibrium uncertainty reflects a mismatch between the true target and the chosen ansatz.
  • Neutrino decoupling: The method reports a conservative uncertainty of 10^-4 because no systematic uncertainty treatment was performed.

D ADDITIONAL RESULTS

The parameter-scan study measures method uncertainty across the full new-physics parameter space using ten independent seeds. The uncertainty is larger than for the unconditional solve and varies with scan complexity.

  • Additional results: The predicted Neff uncertainty across the full parameter space is measured using ten independent seeds.The resulting uncertainty on the mean prediction is reported in Figure 4.
  • Additional results: 10^-4 to 10^-3 is the observed uncertainty range for the new-physics parameter scan.
  • Additional results: The parameter-scan uncertainty is significantly larger than the O(10^-5) uncertainty of the unconditional solve.
  • Additional results: The uncertainty correlates with the size of the probed parameter space and the additional iterations required to cover it fully.
Loading 2608.23022v1…