Animated BBN abundance evolution¶
Creates two animated GIFs of the small-network abundance evolution \(A_i Y_i(t)\):
anim_deltaNeff.gif— sweeps \(\Delta N_{\rm eff}\) from \(-1\) to \(+1\), keeping \(\Omega_b h^2\) at its standard value.anim_omegabh2.gif— sweeps \(\Omega_b h^2\) from \(0.005\) to \(0.05\), keeping \(\Delta N_{\rm eff} = 0\).
All runs use network='small' (8 nuclides, 12 reactions). The data are precomputed
before the animation loop so that GIF rendering is fast.
import importlib.util, subprocess, sys
# In Colab (or any environment without a local primat checkout), pip-install
# the released package; local dev runs already have it importable (editable
# install / repo on sys.path via the cell below), so this is a no-op there.
if importlib.util.find_spec("primat") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "primat"], check=True)
import sys, os
sys.path.insert(0, os.path.abspath('..'))
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
from matplotlib.animation import FuncAnimation, PillowWriter
from primat.backend import run_bbn
from primat.evolution import Y_interpolator
from primat.config import PRIMATConfig
# Uses primat.backend.run_bbn() only (works on either backend); per-frame
# Y(t) interpolators come from primat.evolution.Y_interpolator on each run's
# EvolutionResult (output_time_evolution=True), the same backend-agnostic
# building block primat.plotting.abundance_evolution_curves() uses.
_A_SMALL = {name: NZ[0] + NZ[1] for name, NZ in PRIMATConfig({'network': 'small'}).Nuclides.items()}
# Time grid (1 s → 10^6 s) and species style shared by both animations.
# Each entry: (species key, line color, linestyle, legend label).
t = np.logspace(0, 6, 500)
SPECIES = [
('n', 'black', 'solid', 'n'),
('p', 'red', 'solid', 'p'),
('H2', 'red', 'dashed', 'D'),
('H3', 'red', 'dotted', 'T'),
('He3', 'blue', 'dashed', r'$^3$He'),
('He4', 'blue', 'solid', r'$^4$He'),
('Li7', 'green', 'dashed', r'$^7$Li'),
('Be7', 'purple', 'solid', r'$^7$Be'),
]
# Fixed axis limits that encompass all frames in both animations.
XLIM = (1, 1e6)
YLIM = (1e-12, 2)
FPS = 5 # frames per second in the output GIF
os.makedirs('plots', exist_ok=True)
Helper: precompute abundance tracks¶
Running PRIMAT inside a FuncAnimation update callback would be slow.
We instead precompute the masked \((t, A_i Y_i)\) arrays for every frame
and cache them in a list of dicts.
def precompute_runs(param_values, param_name, base_params=None):
"""
Run the small-network BBN solve (via primat.backend.run_bbn) for each
value in *param_values*, evaluating all species on *t*.
Parameters
----------
param_values : array-like
Sequence of parameter values (one BBN solve per entry).
param_name : str
PRIMATConfig key to vary (e.g. ``'DeltaNeff'`` or ``'Omegabh2'``).
base_params : dict, optional
Fixed parameters shared by all runs (defaults to small network, silent).
Returns
-------
frames : list of dict
Each entry maps species name -> (t_masked, AY_masked) tuple.
"""
base = {'network': 'small', 'verbose': False,
'output_time_evolution': True, 'output_file': None}
if base_params:
base.update(base_params)
frames = []
n = len(param_values)
for i, val in enumerate(param_values):
print(f'\r {param_name} = {val:.4g} ({i+1}/{n})', end='', flush=True)
result = run_bbn({**base, param_name: val})
evolution = result['evolution']
frame = {}
for name, *_ in SPECIES:
AY = _A_SMALL[name] * Y_interpolator(evolution, name)(t)
mask = AY > 0
frame[name] = (t[mask], AY[mask])
frames.append(frame)
print() # newline after progress
return frames
Animation 1 — varying \(\Delta N_{\rm eff}\)¶
\(\Delta N_{\rm eff}\) parametrises extra relativistic energy density beyond the Standard Model value. Increasing it raises the Hubble rate, shifting freeze-out earlier and increasing \(Y_P\); decreasing it has the opposite effect.
delta_neff_values = np.linspace(-2.0, 2.0, 21) # 41 frames, step 0.2
print(f'Precomputing {len(delta_neff_values)} runs for ΔNeff sweep …')
frames_dneff = precompute_runs(delta_neff_values, 'DeltaNeff')
Precomputing 21 runs for ΔNeff sweep …
DeltaNeff = -2 (1/21)
DeltaNeff = -1.8 (2/21)
DeltaNeff = -1.6 (3/21)
DeltaNeff = -1.4 (4/21)
DeltaNeff = -1.2 (5/21)
DeltaNeff = -1 (6/21)
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT.
DeltaNeff = -0.8 (7/21)
DeltaNeff = -0.6 (8/21)
DeltaNeff = -0.4 (9/21)
DeltaNeff = -0.2 (10/21)
DeltaNeff = 0 (11/21)
MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT.
DeltaNeff = 0.2 (12/21)
DeltaNeff = 0.4 (13/21)
DeltaNeff = 0.6 (14/21)
DeltaNeff = 0.8 (15/21)
MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
DeltaNeff = 1 (16/21)
DeltaNeff = 1.2 (17/21)
DeltaNeff = 1.4 (18/21)
DeltaNeff = 1.6 (19/21)
DeltaNeff = 1.8 (20/21)
DeltaNeff = 2 (21/21)
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT.
MT. LT. done.
frames_dneff[20]['n'][1]
array([2.38345489e-01, 2.36982013e-01, 2.35629761e-01, 2.34277292e-01,
2.32982715e-01, 2.31685664e-01, 2.30399491e-01, 2.29145525e-01,
2.27892972e-01, 2.26672028e-01, 2.25474149e-01, 2.24277528e-01,
2.23115454e-01, 2.21962680e-01, 2.20815380e-01, 2.19717771e-01,
2.18620555e-01, 2.17538284e-01, 2.16486126e-01, 2.15436927e-01,
2.14419071e-01, 2.13417634e-01, 2.12418331e-01, 2.11461635e-01,
2.10511215e-01, 2.09569446e-01, 2.08661196e-01, 2.07755800e-01,
2.06870433e-01, 2.06012858e-01, 2.05158900e-01, 2.04329675e-01,
2.03514836e-01, 2.02703969e-01, 2.01932081e-01, 2.01163621e-01,
2.00405234e-01, 1.99675341e-01, 1.98949616e-01, 1.98245537e-01,
1.97557981e-01, 1.96872973e-01, 1.96214283e-01, 1.95566474e-01,
1.94923992e-01, 1.94311267e-01, 1.93700951e-01, 1.93101457e-01,
1.92526415e-01, 1.91954564e-01, 1.91398898e-01, 1.90856426e-01,
1.90317082e-01, 1.89802428e-01, 1.89293962e-01, 1.88790758e-01,
1.88310127e-01, 1.87831588e-01, 1.87363164e-01, 1.86912022e-01,
1.86463673e-01, 1.86029856e-01, 1.85605241e-01, 1.85183120e-01,
1.84779952e-01, 1.84380563e-01, 1.83986852e-01, 1.83610891e-01,
1.83236196e-01, 1.82869252e-01, 1.82513509e-01, 1.82159771e-01,
1.81818844e-01, 1.81483425e-01, 1.81149556e-01, 1.80831346e-01,
1.80514555e-01, 1.80202572e-01, 1.79901536e-01, 1.79601523e-01,
1.79309133e-01, 1.79023372e-01, 1.78738489e-01, 1.78463487e-01,
1.78191394e-01, 1.77920593e-01, 1.77660337e-01, 1.77400354e-01,
1.77144075e-01, 1.76894584e-01, 1.76645188e-01, 1.76401415e-01,
1.76161168e-01, 1.75920829e-01, 1.75687656e-01, 1.75455215e-01,
1.75223475e-01, 1.74998114e-01, 1.74772045e-01, 1.74548282e-01,
1.74328056e-01, 1.74106927e-01, 1.73889416e-01, 1.73672963e-01,
1.73455370e-01, 1.73242377e-01, 1.73028316e-01, 1.72814010e-01,
1.72602735e-01, 1.72389692e-01, 1.72177447e-01, 1.71966066e-01,
1.71752656e-01, 1.71540849e-01, 1.71328004e-01, 1.71112838e-01,
1.70899838e-01, 1.70684147e-01, 1.70466964e-01, 1.70250025e-01,
1.70030084e-01, 1.69809265e-01, 1.69586984e-01, 1.69361359e-01,
1.69135281e-01, 1.68906224e-01, 1.68673706e-01, 1.68440483e-01,
1.68203173e-01, 1.67962836e-01, 1.67720240e-01, 1.67473152e-01,
1.67223317e-01, 1.66969809e-01, 1.66711373e-01, 1.66450325e-01,
1.66184305e-01, 1.65913273e-01, 1.65638882e-01, 1.65358691e-01,
1.65073620e-01, 1.64783846e-01, 1.64487753e-01, 1.64186804e-01,
1.63879882e-01, 1.63566087e-01, 1.63247370e-01, 1.62921457e-01,
1.62588553e-01, 1.62249645e-01, 1.61902817e-01, 1.61548911e-01,
1.61187729e-01, 1.60817974e-01, 1.60440993e-01, 1.60055474e-01,
1.59660731e-01, 1.59258485e-01, 1.58846466e-01, 1.58425032e-01,
1.57994782e-01, 1.57553990e-01, 1.57103550e-01, 1.56642941e-01,
1.56170962e-01, 1.55689060e-01, 1.55195556e-01, 1.54689981e-01,
1.54173709e-01, 1.53644478e-01, 1.53102701e-01, 1.52548462e-01,
1.51980063e-01, 1.51398181e-01, 1.50801500e-01, 1.50189019e-01,
1.49560306e-01, 1.48913078e-01, 1.48246733e-01, 1.47554803e-01,
1.46836243e-01, 1.46083994e-01, 1.45275246e-01, 1.44413017e-01,
1.43434839e-01, 1.42275093e-01, 1.40965749e-01, 1.39009088e-01,
1.36366091e-01, 1.32832757e-01, 1.26586267e-01, 1.18376740e-01,
1.06463817e-01, 9.12241647e-02, 7.43143189e-02, 5.70128790e-02,
4.16549723e-02, 2.78594189e-02, 1.90876780e-02, 1.19203226e-02,
7.17780469e-03, 4.46114035e-03, 2.47467112e-03, 1.46663323e-03,
8.46586952e-04, 4.44389188e-04, 2.75190177e-04, 1.56334334e-04,
8.75964972e-05, 5.70531855e-05, 3.45854498e-05, 2.25553165e-05,
1.56118598e-05, 1.06331153e-05, 7.91796385e-06, 5.97619717e-06,
4.52978608e-06, 3.68422736e-06, 2.98070937e-06, 2.44750137e-06,
2.07883664e-06, 1.75997873e-06, 1.51636337e-06, 1.32320220e-06,
1.15361780e-06, 1.02346270e-06, 9.10729344e-07, 8.10283288e-07,
7.33186004e-07, 6.62133678e-07, 6.00226983e-07, 5.48938287e-07,
5.01451701e-07, 4.60833884e-07, 4.25160055e-07, 3.91997212e-07,
3.64047473e-07, 3.38338507e-07, 3.14513203e-07, 2.94322073e-07,
2.75159893e-07, 2.57766142e-07, 2.42376406e-07, 2.27745452e-07,
2.14701814e-07, 2.02731960e-07, 1.91336684e-07, 1.81306963e-07,
1.71790534e-07, 1.62815935e-07, 1.54836986e-07, 1.47145630e-07,
1.40018685e-07, 1.33477400e-07, 1.27165184e-07, 1.21415335e-07,
1.15986216e-07, 1.10738299e-07, 1.06015071e-07, 1.01446760e-07,
9.70960569e-08, 9.31050996e-08, 8.92193527e-08, 8.55764771e-08,
8.21577299e-08, 7.88267259e-08, 7.57431256e-08, 7.27815111e-08,
6.98959741e-08, 6.72526700e-08, 6.46628245e-08, 6.21789657e-08,
5.98584880e-08, 5.75847445e-08, 5.54375559e-08, 5.33880887e-08,
5.13771753e-08, 4.94969681e-08, 4.76703389e-08, 4.58885331e-08,
4.42306306e-08, 4.25993904e-08, 4.10312986e-08, 3.95447731e-08,
3.80806461e-08, 3.66895991e-08, 3.53464132e-08, 3.40234044e-08,
3.27914738e-08, 3.15795756e-08, 3.03919654e-08, 2.92770517e-08,
2.81778423e-08, 2.71217875e-08, 2.61144974e-08, 2.51210315e-08,
2.41762085e-08, 2.32551598e-08, 2.23443727e-08, 2.15012441e-08,
2.06666769e-08, 1.98493687e-08, 1.90818581e-08, 1.83234567e-08,
1.75918780e-08, 1.68902527e-08, 1.61981395e-08, 1.55422379e-08,
1.49035179e-08, 1.42753955e-08, 1.36886068e-08, 1.31089927e-08,
1.25467950e-08, 1.20134633e-08, 1.14870349e-08, 1.09836412e-08,
1.05002418e-08, 1.00245969e-08, 9.57854582e-09, 9.14245358e-09,
8.71542808e-09, 8.31606449e-09, 7.92299363e-09, 7.54553003e-09,
7.18645718e-09, 6.83295158e-09, 6.49892231e-09, 6.17856343e-09,
5.86460328e-09, 5.57224049e-09, 5.28657593e-09, 5.00958587e-09,
4.75107561e-09, 4.49764337e-09, 4.25710732e-09, 4.02978971e-09,
3.80723361e-09, 3.59957332e-09, 3.39965129e-09, 3.20398875e-09,
3.02521482e-09, 2.85064738e-09, 2.68340409e-09, 2.52858583e-09,
2.37744563e-09, 2.23546202e-09, 2.10145532e-09, 1.97088535e-09,
1.85122150e-09, 1.73632197e-09, 1.62503027e-09, 1.52445659e-09,
1.42652840e-09, 1.33375656e-09, 1.24886746e-09, 1.16630254e-09,
1.08887512e-09, 1.01656914e-09, 9.46822403e-10, 8.83858453e-10,
8.23497811e-10, 7.65788725e-10, 7.13827535e-10, 6.63523701e-10,
6.16573324e-10, 5.73485093e-10, 5.31920323e-10, 4.94151830e-10,
4.58580511e-10, 4.24257746e-10, 3.93893109e-10, 3.64751613e-10,
3.37230866e-10, 3.12608566e-10, 2.88887720e-10, 2.67031073e-10,
2.47003321e-10, 2.27775946e-10, 2.10530399e-10, 1.94345668e-10,
1.78858305e-10, 1.65318762e-10, 1.52297284e-10, 1.40156342e-10,
1.29287006e-10, 1.18865444e-10, 1.09408436e-10, 1.00735011e-10,
9.24444564e-11, 8.51036775e-11, 7.81977383e-11, 7.16636285e-11,
6.59579510e-11, 6.04905446e-11, 5.54533593e-11, 5.09418411e-11,
4.66329681e-11, 4.27714611e-11, 3.92177114e-11, 3.58331753e-11,
3.28826903e-11, 3.00944867e-11, 2.74836610e-11, 2.52015019e-11,
2.30222397e-11, 2.10385262e-11, 1.92556382e-11, 1.75582743e-11,
1.60561319e-11, 1.46682471e-11, 1.33509329e-11, 1.22174800e-11,
1.11410021e-11, 1.01437870e-11, 9.26955513e-12, 8.43724227e-12,
7.68830744e-12, 7.01263001e-12, 6.37134599e-12, 5.81043663e-12,
5.28981342e-12, 4.79843313e-12, 4.37877431e-12, 3.97904297e-12,
3.61276223e-12, 3.29048264e-12, 2.98452429e-12, 2.71229242e-12,
2.46559852e-12, 2.23219893e-12, 2.03055772e-12, 1.84232173e-12,
1.66645650e-12, 1.51579038e-12, 1.37267299e-12, 1.24297578e-12,
1.12831442e-12, 1.01978030e-12, 9.24400646e-13, 8.37480345e-13,
7.55466574e-13, 6.85494992e-13, 6.19810416e-13, 5.59121510e-13,
5.06858677e-13, 4.57361628e-13, 4.13045215e-13, 3.73680433e-13,
3.36523554e-13, 3.04232822e-13, 2.74663392e-13, 2.46863271e-13,
2.23420751e-13, 2.01281253e-13, 1.81058152e-13, 1.63575978e-13,
1.47062832e-13, 1.32437789e-13, 1.19387887e-13, 1.07112246e-13,
9.65783714e-14, 8.68686653e-14, 7.78064753e-14, 7.01996757e-14,
6.30069928e-14, 5.65082585e-14, 5.08634924e-14, 4.55497106e-14,
4.09059848e-14, 3.67349622e-14, 3.28221305e-14, 2.95122685e-14,
2.64434957e-14, 2.36130267e-14, 2.12224745e-14, 1.89706945e-14,
1.69644084e-14, 1.52112901e-14, 1.35656683e-14, 1.21454579e-14,
1.08638556e-14, 9.66644336e-15, 8.66524778e-15, 7.73145961e-15,
6.88109008e-15, 6.16034404e-15, 5.48328584e-15, 4.88680192e-15,
4.36400373e-15, 3.87535862e-15, 3.45896868e-15, 3.08075465e-15,
2.72855782e-15, 2.43888908e-15, 2.16669379e-15, 1.92164595e-15,
1.71309375e-15, 1.51796774e-15, 1.34852969e-15, 1.19893481e-15])
fig, ax = plt.subplots(figsize=(9, 6))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Time [s]')
ax.set_ylabel(r'$A_i\,Y_i$')
ax.set_xlim(*XLIM)
ax.set_ylim(*YLIM)
ax.legend(
handles=[
plt.Line2D([0], [0], color=color, ls=ls, label=label)
for _, color, ls, label in SPECIES
],
ncol=2, loc='center left'
)
# One Line2D artist per species; updated each frame via set_data.
lines = {
name: ax.plot([], [], color=color, ls=ls)[0]
for name, color, ls, _ in SPECIES
}
# Parameter label in the upper-right corner.
param_text = ax.text(
0.97, 0.77, '', transform=ax.transAxes,
ha='right', va='bottom', fontsize=12, fontweight='bold',
path_effects=[pe.withStroke(linewidth=3, foreground='white')]
)
def update_dneff(i):
"""Update all line artists and the annotation for frame *i*."""
frame = frames_dneff[i]
val = delta_neff_values[i]
ax.set_title(
fr'BBN abundance evolution',
fontsize=11
)
param_text.set_text(fr'$\Delta N_{{\rm eff}} = {val:+.2f}$')
for name, *_ in SPECIES:
t_m, AY_m = frame[name]
lines[name].set_data(t_m, AY_m)
return list(lines.values()) + [param_text]
anim_dneff = FuncAnimation(
fig, update_dneff,
frames=len(frames_dneff),
interval=1000 // FPS,
blit=False
)
out_dneff = 'plots/anim_deltaNeff.gif'
anim_dneff.save(out_dneff, writer=PillowWriter(fps=FPS))
plt.close(fig)
print(f'Saved → {out_dneff}')
Saved → plots/anim_deltaNeff.gif
Preview the \(\Delta N_{\rm eff}\) animation:
from IPython.display import Image
Image(out_dneff)
Animation 2 — varying \(\Omega_b h^2\)¶
The baryon density \(\Omega_b h^2\) controls the baryon-to-photon ratio \(\eta_b\). Higher \(\Omega_b h^2\) accelerates nuclear burning, increasing \(Y_P\) slightly, reducing D/H (more efficient deuterium burning), and raising \(^7\)Li/\(^7\)Be.
omegabh2_values = np.linspace(0.005, 0.05, 19) # 19 frames, step ~0.0025
print(f'Precomputing {len(omegabh2_values)} runs for Ωbh² sweep …')
frames_obh2 = precompute_runs(omegabh2_values, 'Omegabh2')
Precomputing 19 runs for Ωbh² sweep …
Omegabh2 = 0.005 (1/19)
Omegabh2 = 0.0075 (2/19)
Omegabh2 = 0.01 (3/19)
Omegabh2 = 0.0125 (4/19)
Omegabh2 = 0.015 (5/19)
Omegabh2 = 0.0175 (6/19)
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
Omegabh2 = 0.02 (7/19)
Omegabh2 = 0.0225 (8/19)
Omegabh2 = 0.025 (9/19)
Omegabh2 = 0.0275 (10/19)
Omegabh2 = 0.03 (11/19)
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT.
Omegabh2 = 0.0325 (12/19)
Omegabh2 = 0.035 (13/19)
Omegabh2 = 0.0375 (14/19)
Omegabh2 = 0.04 (15/19)
Omegabh2 = 0.0425 (16/19)
MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
Omegabh2 = 0.045 (17/19)
Omegabh2 = 0.0475 (18/19)
Omegabh2 = 0.05 (19/19)
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
[primat] HT. MT. LT. done.
fig, ax = plt.subplots(figsize=(9, 6))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Time [s]')
ax.set_ylabel(r'$A_i\,Y_i$')
ax.set_xlim(*XLIM)
ax.set_ylim(*YLIM)
ax.legend(
handles=[
plt.Line2D([0], [0], color=color, ls=ls, label=label)
for _, color, ls, label in SPECIES
],
ncol=2, loc='center left'
)
lines = {
name: ax.plot([], [], color=color, ls=ls)[0]
for name, color, ls, _ in SPECIES
}
param_text = ax.text(
0.97, 0.7, '', transform=ax.transAxes,
ha='right', va='bottom', fontsize=12, fontweight='bold',
path_effects=[pe.withStroke(linewidth=3, foreground='white')]
)
def update_obh2(i):
"""Update all line artists and the annotation for frame *i*."""
frame = frames_obh2[i]
val = omegabh2_values[i]
ax.set_title(
fr'BBN abundance evolution',
fontsize=11
)
param_text.set_text(fr'$\Omega_b h^2 = {val:.4f}$')
for name, *_ in SPECIES:
t_m, AY_m = frame[name]
lines[name].set_data(t_m, AY_m)
return list(lines.values()) + [param_text]
anim_obh2 = FuncAnimation(
fig, update_obh2,
frames=len(frames_obh2),
interval=1000 // FPS,
blit=False
)
out_obh2 = 'plots/anim_omegabh2.gif'
anim_obh2.save(out_obh2, writer=PillowWriter(fps=FPS))
plt.close(fig)
print(f'Saved → {out_obh2}')
Saved → plots/anim_omegabh2.gif
Preview the \(\Omega_b h^2\) animation:
from IPython.display import Image
Image(out_obh2)
