Big Bang nucleosynthesis

Most Boltzmann codes take the helium fraction from an external table. Here the network is part of the package: the weak $n \leftrightarrow p$ rates with their finite-temperature QED and thermal corrections, a reaction network over the light nuclides with the PRIMAT rate compilation (each of the 63 reaction rates shipped as a data file with provenance, plus detailed-balance closures), and the nuclear evolution integrated through the same background the rest of the code uses — including any exotic components, since the Friedmann equation is assembled from the species list.

This buys three things. First, consistency: cosmology() with the default Yp = :bbn gives every downstream module — recombination, CMB damping tail, diffusion distortions — a helium abundance that actually corresponds to its baryon density and expansion history. Second, reach: deuterium, helium-3, lithium — the observables of BBN — come out of the same call (bbn, primordial_bbn), with Monte-Carlo rate-uncertainty propagation in bbn_mc. Third, a gate for new physics: anything that changes the expansion rate or injects energy (extra radiation from decaying dark matter, a gravitational-wave background contributing to $N_{\rm eff}$, lepton asymmetries) feeds through to abundances that can be checked against observation within the same package.

Two caveats worth knowing. Rate compilations differ: PRIMAT and PArthENoPE disagree at the 2.6% level on D/H and more on lithium-7, so comparisons should name the compilation. And the quoted abundances inherit the network's truncation — the full CNO extension relevant to the lithium problem at the $10^{-4}$ level is not yet included.

Reference

Cosmic.T_nu_over_T_gamma — Method
T_nu_over_T_gamma(T_γ)

Neutrino-to-photon temperature ratio, from entropy conservation in the γ-e± plasma.

a·Tν is constant (neutrinos are decoupled and free), and s{γe}·a³ is constant (photons and e± stay coupled), so Tν/Tγ = [ŝ(T_γ)/ŝ(∞)]^{1/3} with ŝ = s/T³. Taking T → 0 gives back (4/11)^{1/3} without it ever being typed in.

source
Cosmic.bbn — Method
bbn(c::Cosmology; kwargs...)

Run nucleosynthesis for the cosmology c, closing the loop: ωb and Neff are read off the species list rather than passed separately, so the helium that comes out is the helium that that baryon density and radiation content actually make.

Any relativistic species beyond γ and the standard neutrinos – a dark-radiation component, a light relic – shows up in the expansion rate during nucleosynthesis and therefore in Yp, which is the whole reason Yp is a useful probe of N_eff. Passing a friedmann hook is the entry point for modified gravity: BBN feels a modified H(ρ) through the freeze-out of the weak rates and through nothing else.

Matter, curvature, and the cosmology's exact dark-energy ρ(a) are retained too; their standard-ΛCDM contributions are tiny at BBN, but early-dark-energy models are therefore handed to the network without a radiation-era substitution.

source
Cosmic.bbn — Method
bbn(; ω_b = 0.02242, ω_m = 0, N_eff = nothing, τ_n = 878.4, ξ_ν = 0,
    qke_history = nothing, qke_electron_flavor = 1,
    T9_weak = 100, T9_start = 10, T9_end = 0.005)

Integrate the nuclear network and return the primordial abundances.

τ_n is the neutron lifetime in seconds (PDG 2024: 878.4 ± 0.5). It is the single measured number the weak rates need, and Yp is genuinely sensitive to it: ±0.5 s moves Yp by about ±0.0002.

ξ_ν = μ_νe/T_νe is the electron-neutrino degeneracy. It enters the Born, finite-nucleon-mass, and finite-temperature QED weak-rate terms with the neutrino/antineutrino signs fixed by detailed balance. N_eff is the total relativistic energy density, including any ξ-induced enhancement; bbn(c) computes that value from the cosmology automatically.

thermal_weak=false is the PRIMAT-style diagnostic switch that omits only the finite-temperature Brown-Sawyer QED term. It defaults to true; the Born, zero-temperature radiative, Coulomb, and finite-nucleon-mass terms remain active.

qke_history supplies a checked FlavorDensityHistory. Its exact energy density, scale-factor history, and electron-flavour occupation then replace the effective-temperature decoupling backend consistently. Since the public FortEPiaNO history has zero lepton asymmetry, this mode rejects ξ_ν != 0. If N_eff is omitted it is taken from the final QKE history; an explicit different value represents additional or missing free radiation. qke_electron_flavor identifies the weak-interaction electron flavour in a backend with arbitrary flavour ordering.

The integration variable is ln T9 rather than time. Time never appears explicitly: a·Tν is constant, so H = -dln Tν/dt gives dt/dln Tγ = -(dln Tν/dln Tγ)/H, and the whole thermal history follows from the entropy of the γ-e± plasma. `ωmsupplies the present physical matter density to the Friedmann rate;bbn(cosmology)` passes the cosmology's baryon-plus-CDM value automatically.

source
Cosmic.bbn_mc — Method
bbn_mc(; n = 400, στ_n = 0.5, seed = 0x5eed, kwargs...) -> NamedTuple

Monte-Carlo propagation of the nuclear-rate and neutron-lifetime uncertainties through the network. Each draw displaces every forward rate by ξᵢ standard deviations of its own tabulated 1σ envelope (λ = median·σ(T9)^ξᵢ, ξᵢ ~ N(0,1) – the PArthENoPE/PRIMAT convention; the error columns in the rate files finally earn their keep) and draws τn ~ N(τn, στ_n). Reverse rates follow each draw automatically through detailed balance.

Returns means, standard deviations and 16/50/84 percentiles for Y_p, D/H, ³He/H and ⁷Li/H. Draws run threaded; kwargs are passed through to bbn.

source
Cosmic.detailed_balance — Method
detailed_balance(reactants, products) -> (α, β, γ)

Reverse-rate coefficients derived from masses and spins (see the note above). Takes nuclide indices as stored in a Reaction; the photon, if any, is simply absent from the product list.

source
Cosmic.weak_rate — Method
weak_rate(T_γ, T_ν, τ_n; direction)

n→p (direction = +1) or p→n (direction = -1) conversion rate in 1/s.

Γ = (1/(τ_n λ₀)) ∫₀^∞ dp p² [ χ(E) + χ(-E) ],   E = √(p²+1)
χ(E) = E_ν² · g_ν(E_ν, x_ν) · g(-E, x),   E_ν = E - s·q

Temperatures in MeV. See the module docstring for why the χ(E) + χ(-E) form is the right way to write this.

source