BBN-simple: a code for simplified and pedagogical Big Bang Nucleosynthesis#
This code is a pedagogical approach to light-element Big Bang Nucleosynthesis (BBN) stripped to its simplest physics. It is aimed at an advanced undergraduate or early graduate student interested in the most fundamental and necessary physics required to consider the problem of BBN. You can download the code as an .ipynb or this site compiled as a .pdf in the upper right of the page.
This code and its associated publication were written by Aidan Meador-Woodruff and Dragan Huterer of the University of Michigan in the Lienweber Center for Theoretical Physics. Our emails are aidanmw@umich.edu and huterer@umich.edu, respectively.
Section 0: Constants, Functions, and Import Statements: Typical Pythonic Things#
### IMPORT STATEMENTS ###
import numpy as np
import math
from scipy.special import *
from scipy.optimize import *
from scipy import interpolate
from scipy.integrate import *
import matplotlib
from matplotlib import pyplot as plt
# ----------------------------------- #
### CONSTANTS ###
## Numbers ##
pi = np.pi
apery = zeta(3)
## Masses ##
me = 0.510998 # Mass of electron in [MeV/c^2]
mp = 938.272013 # Mass of proton in [MeV/c^2]
mn = 939.565346 # Mass of neutron in [MeV/c^2]
md = 1875.612793 # Mass of Deuterium in [MeV/c^2]
mt = 2808.920906 # Mass of Tritium in [MeV/c^2]
m3 = 2808.391383 # Mass of Helium-3 in [MeV/c^2]
m4 = 3727.379109 # Mass of Helium-4 (Alpha) in [MeV/c^2]
mli6 = 5603.05141 # Mass of Lithium-6 in [MeV/c^2]
mli7 = 6533.833166 # Mass of Lithium-7 in [MeV/c^2]
mbe7 = 6534.184060 # Mass of Beryllium-7 in [MeV/c^2]
Q = 939.5654133 - 938.2720813 # Proton-Neutron Mass Difference in [MeV/c^2]
amu = 931.494003 # AMU in MeV
q = Q/me # Exponential NSE Mass ratio [unitless]
## Energies ##
Bd = 2.224 # Binding Energy of Deuterium in [MeV]
Bt = 8.50 # Binding Energy of Tritium in [MeV]
B3 = 7.718 # Binding Energy of Helium-3 in [MeV]
B4 = 28.295 # Binding Energy of Helium-4 (Alpha) in [MeV]
Bli7 = 39.2 # Binding Energy of Lithum-7 in [MeV]
Bbe7=37.6 # Binding Energy of Beryllium-7 in [MeV]
## General Physical Units ##
G = 6.67e-45 # Gravitational Constant []
c = 2.99e10 # Speed of Light [cm/s]
kB = 8.617e-11 # Boltzmann Constant [MeV/Kelvin]
hbar = 6.582e-22 # Reduced Planck constant h/2pi [MeV s]
NA = 6.022e23 # Avogadro's number [# molecules/mole]
## Particle Physics Units ##
gA = 1.26 # Axial vector corretion for weak charge nucleon rxns between p and n. [scalar]
GF = 1.1663787e-5*1e-6 # Fermi coupling constant in [MeV^-2]
RadConst = (np.pi**2)/15 # Radiation constant for Radiation pressures [unitless]
## Free Parameters ##
Neff = 3.046 # Number of neutrino flavors [unitless]
taun = 880 # Average lifetime of a neutron before decay [s]
η0 = 6.12e-10 # The present day ratio of baryons to photons [unitless, number density]
# ----------------------------------- #
### USEFUL CONVERSIONS ###
def K9toMeV(T):
'Converts T9 to T in MeV.'
K = 8.619e-5 # 1 Kelvin to eV/kB
K9 = K*1e9 # 10^9 Kelvin in eV/kB
MeV = K9/1e6 # 10^9 Kelvin in MeV/kB
return MeV*T
def MeVtoK9(T):
'Converts T in MeV to T9'
eV = 1/(8.619e-5) # 1 eV/kB to Kelvin
MeV = eV*1e6 # Converting MeV to eV/kB
K9 = MeV/1e9 # Converting Kelvin to 10^9 Kelvin
return K9*T
Section 1: Thermodynamics and the Time-Temperature Relationship#
To begin with BBN, we must lay the groundwork of thermodynamics in the early universe. In our cake analogy, this is the equivalent of creating our oven and preheating it to the appropriate temperature. The evolution of photon and neutrino temperature with resepct to time are driven by the expansion of space after the Big Bang. The expansion of space and its thermodynamic consequences are the primary driver of BBN.
To relate the evolution of time after the Big Bang and the temperature of radiation, we must use the Friedmann equations. We will use the fiducial values for a Flat, \(\Lambda\) CDM universe in this notebook. We will assume a basic understanding of cosmology, at least to the degree of Kolb & Turner’s The Early Universe (1994) or Ryden’s Introduction to Cosmology (2016).
Friedmann Equations#
There are two independent Friedmann equations. They are, (for zero curvature):
and $\(\frac{\ddot a}{a} = - \frac{4\pi G}{3} (\rho + 3P)\)$
Where \(a\), \(\rho\), and \(P\) are the scale factor, energy densities, and radiation pressure, respectively. A third crucially important equation is the so-called continuity equation, which lets us calculate the time derivatives of energy densities. It is
By use of the First Law of Thermodynamics and the thermodynamic integrability condition, we arrive at the equivalent equation
Which is explained from fundamentals in K&T and is equivalent to Equation 3.67. With these equations in hand, we have the basis to derive the time-temperature relationships.
Time-Temperature Relationships#
From Table 1 in Wagoner, Fowler, and Hoyle’s 1967 Synthesis of Elements, the Freidmann equations can be used to express the following differential relationship between temperature of photons and time
Which they show can be integrated (from \(t=0\) and \(T = \infty\)) to give the necessary relationship
A common way to express the temperature of photons in the early universe is in units of \(10^9 ~ {\mathrm{K}}\). Temperature in these units is usually denoted as \(T_9\). This equation is in units of \(\mathrm{\frac{MeV}{k_{B}}}\). However, these units can be converted to Kelvin and thus this can be expressed as the approximation
One can also solve the Friedmann equations, noting that
and that $\(\frac{dt}{dT} = -\frac{d \ln(a)}{dT} \frac{1}{H}\)$ Using these two equations, we may evolve the neutrino temperature and time appropriately.
### ELECTRON/POSITRON AND PHOTON TEMPERATURES ####
# Initial and final temperatures
Ti = 1e1
Tf = 1e-2
n = 512
T = np.geomspace(Ti,Tf,n)
T9i = MeVtoK9(Ti)
T9f = MeVtoK9(Tf)
T9 = MeVtoK9(T)
# Equations 18 and 19:
def ElectronPositron(T):
'''
Solves for the electron and positron energy density and pressure, 𝜌e and Pe, as a function of temperature.
INPUT :
Tγ in MeV [1xn ndarray]
OUTPUTS :
𝜌e [1xn ndarray]
Pe [1xn ndarray]
'''
# GL Quadrature shorthands and initialization
N = 64
points, weights = roots_laguerre(N)
y = points
x = np.expand_dims(me/T,-1)
dist = np.sqrt(x**2 + y**2)
rhoe_numerator = y**2 * dist
rhoe_denominator = np.exp(dist) + 1
# Have to reweight for GL quadrature
rhoe_integrand = rhoe_numerator/rhoe_denominator * np.exp(y)
# Approximates the integral by a finite sum in quadrature
rhoe_integral = np.dot(rhoe_integrand, weights)
# Full expression with coefficient
rhoe = T**4 * (2/pi**2) * rhoe_integral
Pe_numerator = y**2
Pe_denominator = dist * (np.exp(dist) + 1)
# Likewise, weight for GL quadrature
Pe_integrand = y**2*Pe_numerator/Pe_denominator * np.exp(y)
Pe_integral = np.dot(Pe_integrand, weights)
# Full expression with coefficient
Pe = T**4 * (2/(3 * pi**2)) * Pe_integral
#1/T Derivative of rhoe with respect to temperature. Specific heat, which is used to compute the derivatives dTnu/dT and dT/dt
drhoedT_integrand = y**2 * dist**2 * np.exp(dist) /(rhoe_denominator) * np.exp(y)/(rhoe_denominator)
drhoedT_integral = np.dot(drhoedT_integrand,weights)
drhoedT = T**4 * (2/pi**2) * drhoedT_integral
return rhoe, Pe, drhoedT
# Equations 11 and 12:
def Photon(T):
'''
Solves for the photon energy density and pressure, 𝜌γ, Pγ, as a function of temperature.
INPUT :
Tγ in MeV [1xn ndarray]
OUTPUTS :
𝜌γ in dynes per cm^3 [1xn ndarray]
Pγ in [pressure...] [1xn ndarray]
'''
rhogamma = pi**2/15 * T**4
Pgamma = rhogamma/3
return rhogamma, Pgamma
def Neutrino(Tnu):
rhoneutrino =7/8*pi**2/15 *Neff* Tnu**4
Pneutrino = rhoneutrino/3
return rhoneutrino, Pneutrino
def Friedmann(T,a):
Tnu = a[0]
rhog, Pg = Photon(T)
rhonu,Pnu = Neutrino(Tnu)
rhoe, Pe,cve = ElectronPositron(T)
drhodT = (cve+4*rhog)/T/3
H = np.sqrt(rhoe + rhonu + rhog) * np.sqrt(8*pi*G/3)/hbar
dtdT =-(drhodT)/(rhoe +rhog + Pe + Pg)/H
dTnudT = drhodT/(rhoe +rhog + Pe + Pg)*Tnu
return [dTnudT,dtdT]
FriedmannIntegration= solve_ivp(Friedmann,T[[0,-1]],[Ti,0], t_eval=T)
T = FriedmannIntegration.t
Tnu = FriedmannIntegration.y[0]
t = FriedmannIntegration.y[1]
T9 = MeVtoK9(T)
Tnu9 = MeVtoK9(Tnu)
#Tnu9 = ((11/4 * (rhogamma+Pgamma))/(rhogamma+rhoe+Pgamma+ Pe))**(-1/3) * T9
#Tnu = K9toMeV(Tnu9)
# Equation 8 in 10^9 K and seconds:
#t= ((Tnu9)/10.4)**(-2)
Show code cell source
### Plot of Neutrino-Photon Decoupling ###
matplotlib.rcParams['mathtext.fontset'] = 'cm'
default_fontfam = matplotlib.rcParams['font.family']
plt.rcParams["font.family"] = "serif"
plt.rcParams["mathtext.fontset"] = "dejavuserif"
## Initialize the figure and axis
fig = plt.figure(figsize=(9,7))#, dpi=200)
ax = plt.subplot(111)
ax.set_xlim(5e-1,1e4)
ax.set_ylim(5e-2,5e1)
## RC Parameters
plt.rcParams['lines.linewidth'] = 4
ax.tick_params(axis='both', direction='in', which='major', length = 9, width=1.3, right=False, labelsize=24, pad=8)
ax.tick_params(axis='x', direction='in', length=7, width=1, top=False, which='minor', labelsize=24, pad=8, reset=True)
ax.tick_params(axis='y', direction='in', length=0, width=0, right=False, which='minor', labelsize=24)
## Labels
#ax.set_title(r'Neutrino Decoupl|ing', fontsize=16)
ax.set_ylabel(r'Temperature $(10^{9} $ K$)$', fontsize=30);
ax.set_xlabel(r'Time (Seconds)', fontsize=30);
## Plot
ax.loglog(t, T9, color='k');
ax.loglog(t, Tnu9, color='red');
## Captions
labels = [r"$T_{\gamma}$ ", r"$T_{\nu}$"]
legend = ax.legend(labels, fontsize=24)
Section 2: The Weak Rates and Nuclear Statistical Equilibrium#
The weak interactions are those between protons and neutrons. Namely, the following reactions:
We denote the rate \(\lambda_{n \rightarrow p}\) as the sum of all of the forward weak reactions and likewise \(\lambda_{p \rightarrow n}\) as the sum of all the reverse weak reactions. In our paper, these are given by
and
Where we have that the variables:
\(q = (m_n - m_p)/m_e\) is the neutron-proton mass difference divided by the electron mass
\(z = m_ec^2/k_B T_\gamma\) is a substitution
\(z_\nu = m_ec^2/k_B T_\nu\) is the same as \(z\) but for the neutrino temperature
\(K\) is the normalization constant, set so that these integrals go to \(1/\tau\) for low temperatures.
\(\xi_e\) is the chemical potential of the electron, taken to be zero.
We elect to solve this by Gauss-Laguerre quadrature below.
### WEAK RATE CALCULATION ###
# Equations 27 and 28:
def WeakRate(T,Tnu):
'''
Solves for the weak rates λnp and λpn by approximating the integral in quadrature.
INPUT :
Tγ in MeV [1xn ndarray]
Tν in MeV [1xn ndarray]
OUTPUTS :
λnp [1xn ndarray]
λpn [1xn ndarray]
'''
# The normalization coefficient
K = 1/(quad(lambda x: np.sqrt(x**2-1)*x*(q-x)**2, 1, q)[0])* (1/taun)
# Quadrature setup
N = 64
points, weights = roots_laguerre(N)
# Shorthands
x = points
z = np.expand_dims(me/(T),-1)
x = points + 1 #Shift integral from 0 to infinity to 1 to infinity.
znu = np.expand_dims(me/(Tnu), -1)
# Define both integrals, sum, and normalize
int1 = (x*(x+q)**2 * (x**2 - 1)**(1/2))
denom1 = ((1+np.exp(-x*z))*(1+np.exp((x+q)*znu)))
int2 = (x*(x-q)**2 * (x**2 - 1)**(1/2))
denom2 = ((1+np.exp(x*z))*(1+np.exp(-(x-q)*znu)))
lpn = K*(int1/denom1+int2/denom2)
# Likewise for reverse
rint1 = (x*(x-q)**2 * (x**2 - 1)**(1/2))
rdenom1 = ((1+np.exp(-x*z))*(1+np.exp((x-q)*znu)))
rint2 = (x*(x+q)**2 * (x**2 - 1)**(1/2))
rdenom2 =((1+np.exp(x*z))*(1+np.exp((-(x+q))*znu)))
lnp = K*(rint1/rdenom1 + rint2/rdenom2)
# Weight for GL Quadrature
lnp = lnp*np.exp(x-1)
lpn = lpn*np.exp(x-1)
# Approximates the integral
lnp = np.dot(lnp, weights)
lpn = np.dot(lpn,weights)
return lnp, lpn
weak_np, weak_pn= WeakRate(T,Tnu);
Show code cell output
/tmp/ipykernel_6007/1673097470.py:33: RuntimeWarning: overflow encountered in exp
denom1 = ((1+np.exp(-x*z))*(1+np.exp((x+q)*znu)))
/tmp/ipykernel_6007/1673097470.py:35: RuntimeWarning: overflow encountered in exp
denom2 = ((1+np.exp(x*z))*(1+np.exp(-(x-q)*znu)))
/tmp/ipykernel_6007/1673097470.py:41: RuntimeWarning: overflow encountered in exp
rdenom1 = ((1+np.exp(-x*z))*(1+np.exp((x-q)*znu)))
/tmp/ipykernel_6007/1673097470.py:44: RuntimeWarning: overflow encountered in exp
rdenom2 =((1+np.exp(x*z))*(1+np.exp((-(x+q))*znu)))
Show code cell source
### Plot of Weak Reactions.
## Initialize the figure and axis
fig = plt.figure(figsize=(9,7))
ax = plt.subplot(111)
ax.tick_params(axis='both', direction='in', which='major', length = 7, width=1.15, right=True, labelsize=24, pad=8)
ax.tick_params(axis='x', direction='in', length=0, width=0.8, right=False, which='minor', labelsize=24)
ax.tick_params(axis='y', direction='in', length=0, width=0.6, right=False, which='minor', labelsize=24)
## Labels
ax.set_ylabel(r'Reaction Rate', fontsize=30); # ($\rm{cm^3 s^{-1}}$)'
ax.set_xlabel(r'Temperature ($10^9 ~ \rm{K}$)', fontsize=30);
## Plot
ax.loglog(T9,weak_np, color='k');
ax.loglog(T9, weak_pn, color='red');
ax.hlines(taun**(-1), 1e2,1e-3, linestyle = 'dashed', color='deepskyblue', linewidth = 3.5)
ax.set_ylim(1.1e-10,5e1)
ax.set_xlim(1e2,1e-1)
## Legend
labels = [r"$n \rightarrow p$", r"$p \rightarrow n$", r"$\tau^{\,-1}$"]
legend = ax.legend(labels,fontsize=24);
#plt.savefig('/home/aidan/WeakRates.pdf', dpi=250, bbox_inches='tight')
In equilibrium (for a photon temperature of >0.8 MeV or so), the number density of nuclear species is governed by the Saha equation. We can express the mass fraction of elements immediately before BBN begins by knowing only the binding energy and spin states of each nucleus along with the equilibrium mass fractions of protons and neutrons. This is because at high temperatures, all species are in equilibrium with protons and neutrons. In NSE, we can express the number density of nuclei of species \(X\) , which would be a light nuclei in our case, such as Deuterium or Tritium, as:
Where
\(n_n\), \(n_p\), \(n_X\) are the number densities of neutrons, protons, and species \(X\) respectively
\(m_{s,X}\) is the number of possible spin states of nucleus \(X\)
\(m_n\), \(m_p\), \(m_X\) are the masses of neutrons, protons, and species \(X\) and
\(T\) is temperature in MeV
This can be reexpressed as a more useful quantity for our case by expressing the number density of a species with respect to the number density of baryons overall. That rexpression, given in Dodelson and Schmidt’s Modern Cosmology (equation 4.16), is
Where \(\eta_b\) is the baryon to photon ratio and \(B_X\) is the binding energy of a nucleus of species \(X\). This relation holds when \(\eta_b\) dominates the exponential term, and thus holds well for immediately pre-BBN temperatures, since they are not too high. From Kolb and Turner equation 4.7, we have that the NSE mass fraction of a species X can be expressed as (the somewhat unwieldy) expression:
These allow us to find our initial conditions.
### NUCLEAR STATISTICAL EQUILIBRIUM (NSE) ###
def NSEAbundances(T,Bind, A, gA):
'''
T - Temperature in MeV
Bind - Binding energy in MeV
A - Mass number of the atom
gA -
'''
coeff = gA * η0**(A-1) * (apery**(A-1) * pow(np.pi, ((1-A)/2)) * pow(2, ((3*A-5)/2))) * A**(5/2)
XA = coeff * (T/mn)**((3*(A-1))/2) * np.exp(Bind/T)
return XA
# Initial conditions for BBN code, NSE abundances of elements.
# Out-front factor is assuming initial conditions at our starting time have Xp ~ Xn, so it is (1/2)^A where A is the mass number.
XDNSE0 = (1/4) * NSEAbundances(T9i, Bd, 2, 3)
XTNSE0 = (1/8) * NSEAbundances(T9i, Bt, 3, 2)
X3NSE0 = (1/8) * NSEAbundances(T9i, B3, 3,2)
X4NSE0 = (1/16) * NSEAbundances(T9i, B4, 4, 1)
XLi7NSE0 = (1/128) * NSEAbundances(T9i, Bli7, 7, 5)
XBe7NSE0 = (1/128) * NSEAbundances(T9i, Bbe7, 7, 5)
X0 = np.zeros(8,)
X0[1] = (np.exp(Q/(T9i))+1)**(-1)
X0[0] = 1 - X0[1]
X0[2] = XDNSE0
X0[3] = XTNSE0
X0[4] = X3NSE0
X0[5] = X4NSE0
X0[6] = XLi7NSE0
X0[7] = XBe7NSE0
Section 3: The Strong Rates and The Reaction Network#
All reactions between nucleides interacting by way of strong reactions, that is atomic reactions, must also be considered in the big bang nucleosynthesis model. These reactions consist of no more than 4 distinct nuclei, which we can express schematically as the nuclear reaction
Where we have two species, \(Y_i\) and \(Y_j\) combining with a rate \([ij]_{kl}\) to produce the two species \(Y_k\) and \(Y_l\). The reverse rate, \([kl]_{ij}\) is the rate at which \(Y_k\) and \(Y_l\) combine and produce \(Y_i\) and \(Y_j\). We are interested in the forward and reverse rates for our reaction network. They can be computed and generally take the form
Where \(\rho_b\) is the baryon density, \(N_A\) is Avogadro’s number, and \( \langle\sigma v\rangle_{jk}\) is an experimentally determined cross-section. The reverse rate is $\([kl]_{{ij}} = \frac{g_{i}g_{j}}{g_{k}g_{l}} \left( \frac{A_{i}A_{j}}{A_{k}A_{l}} \right)^{3/2} \exp\left( \frac{Q}{k_BT} \right) [ij]_{kl}\)$
Or, in the case of one product, $\([k\gamma]_{ij} = \frac{g_{i}g_{j}}{(1+\delta_{{ij}})g_{k}} \left( \frac{A_{i}A_{j}}{A_{k}} \right)^{3/2} \rho_{b}^{-1}T_{9}^{3/2}[ij]_{k\gamma} \exp\left( \frac{Q}{k_BT}\right).\)$
Here we adopt the rates from Smith et al., as in our paper, and ignore any uncertainties in them. This choice provides a straightforward one-stop-shop for the rates. Our emphasis on simplicity rather than precision justifies the ignorance of any updated best-fit values in these rates, or the uncertainties in them.
We will consider the 12 reactions below involving/creating protons, neutrons, Deuterium, Tritium, Helium-3, Helium-4, Lithium-6, Lithium-7, and Beryllium-7
While other elements are involved in BBN, these represent the vast majority of reactants by mass fraction. The reaction network can be constructed as follows
Which allows us to create a coupled system of first-order non-linear ODEs.
def StrongRates(T9, Tnu9):
"""
Approximate the strong reaction rates using Kawano & Smith's 1993 polynomial approximations.
INPUT:
T9 - Temperature array in 10^9 Kelvins. TYPE : [1xn] or [n] ndarray or equivalent.
OUTPUTS:
Forward - Reaction rates in cubic cm per second per mole. TYPE : [11xn] ndarray
Reverse - Reverse reaction rates in the same units. Same type.
"""
# A word of warning to those who want to code this-- Use an Einstein sum. Especially for bigger networks/stellar nucleosynthesis.
# I had already written the reaction network by the time I figured out that trivializes it and makes it much more compact.
#T = K9toMeV(T9)
h= 3.33683e4*η0
rhob = h*(Tnu9)**3 *Neff#* 7/8*(15/np.pi**2)
# Shorthands introduced by Kawano 1993
T9a = T9/(1 + 0.1378*T9)
T9b = T9/(1 + 0.1071*T9)
T9c = T9/(1 + 0.759*T9)
T9d = T9/(1 + 13.08*T9)
# Shorthands
T913 = T9**(1/3)
T953 = T9*(5/3)
lt9 = np.log10(T9)
def fit(a0,a1,a2,a3,a4,a5,a6):
rate = np.exp(a0 + a1/T9 + a2/T913 + a3*T913 + a4*T9 + a5*T953 + a6*np.log(T9))
return rate
T9 = np.asarray(T9)
K = []
if T9.any() < 0.1:
K = 2.60191e+04*(1 -3.1293e-01*lt9**1 +5.5526e-01*lt9**2+5.4925e-01*lt9**3 +1.7919e-01*lt9**4 +2.0267e-02*lt9**5)
elif T9.any() < 10:
K = 2.61952e+04*(1 -2.5918e-01*lt9**1 +6.1972e-01*lt9**2+4.4927e-01*lt9**3 -4.9802e-02*lt9**4 -9.0029e-02*lt9**5)
else:
K = 2.56323e+04*(1 -6.5312e-02*lt9**1 -4.9220e-02*lt9**2+1.5391e+00*lt9**3 -8.4075e-01*lt9**4 +1.2262e-01*lt9**5)
# All updated rate parameters from REACLIB DB.
def fit(a0,a1,a2,a3,a4,a5,a6):
rate = np.exp(a0 + a1/T9 + a2/T913 + a3*T913 + a4*T9 + a5*T953 + a6*np.log(T9))
return rate
# Forward Rates are polynomial fits, pulled from Table 2 of Kawano 1993.
# They are reaction cross-sections, <sigma v> * NA
crosssection = np.asarray([
#1. p + n --> D + γ
4.742e4 * (1 - 0.8504 * T9**(1/2) + 0.49 * T9 - 0.0962 * T9**(3/2) + 8.47e-3 * T9**(2) - 2.8e-4 * T9**(5/2)),
#2. D + p --> He3 + γ
fit(7.528980e+00,0.000000e+00,-3.720800e+00,8.717820e-01,0.000000e+00,0.000000e+00,-6.666670e-01) + fit(8.935250e+00,0.000000e+00,-3.720800e+00,1.986540e-01,0.000000e+00,0.000000e+00,3.333330e-01),
#3. D + D --> n + He3
fit(1.908760e+01,-1.900200e-04,-4.229200e+00,1.693200e+00,-8.555290e-02,-1.357090e-25,-7.345130e-01),
#4. D + D --> p + T
fit(1.880520e+01,4.362090e-05,-4.322960e+00,1.915720e+00,-8.156200e-02,-3.288040e-22,-8.795180e-01),
#5. He3 + n --> p + T
fit(2.037870e+01,0.000000e+00,0.000000e+00,-3.327880e-01,-7.004850e-01,9.765210e-02,0.000000e+00) + fit(1.927620e+01,0.000000e+00,0.000000e+00,4.385570e-02,-2.015270e-01,1.534330e-02,1.000000e+00),
#6. T + D --> n + He4
1.063e11 * T9**(-2/3) * np.exp(-4.559/(T9**(1/3)) - (T9/0.0754)**2) * (1 + 0.092*T9**(1/3) - 0.375*T9**(2/3) - 0.242*T9 + 33.82*T9**(4/3) + 55.42*T9**(5/3)) + 8.047e8*T9**(-2/3) * np.exp(-0.4857/T9),
#fit(3.934570e+01,0.000000e+00,-4.524400e+00,-1.640280e+01,1.731030e+00,-1.229660e-01,2.313040e+00) + fit(2.517940e+01,0.000000e+00,-4.524400e+00,3.503370e-01,5.874700e-01,-8.849090e+00,-6.666670e-01),
#7. He3 + D --> p + He4
5.021e10 * T9**(-2/3) * np.exp(-7.144/(T9**(1/3)) - (T9/0.27)**2) * (1 + 0.058*T9**(1/3) + 0.603*T9**(2/3) + 0.245*T9 + 6.97*T9**(4/3) + 7.19*T9**(5/3)) + 5.212e8*T9**(-1/2) * np.exp(-1.762/T9),
#fit(2.468390e+01,0.000000e+00,-7.182000e+00,4.732880e-01,1.468470e+00,-2.796030e+01,-6.666670e-01) + fit(4.129690e+01,0.000000e+00,-7.182000e+00,-1.713490e+01,1.369080e+00,-8.144230e-02,3.353950e+00),
#8. He3 + He4 --> Be7 + γ
4.817e6 * T9**(-2/3) * np.exp(-14.964/(T9**(1/3))) * (1 + 0.0325*T9**(1/3) - 1.04e-3 * T9**(2/3) - 2.37e-4 * T9 - 8.11e-5 * T9**(4/3) - 4.69e-5 * T9**(5/3)) + 5.938e6*T9b**(5/6)* T9**(-3/2) * np.exp(-12.859/(T9b**(1/3))),
# UHHHHHH
#fit(1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01) + fit(1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00),
#9. T + He4 --> Li7 + γ
fit(1.361620e+01,0.000000e+00,-8.080500e+00,-2.175140e-01,-1.148590e-01,4.700430e-02,-6.666670e-01),
#10. Be7 + n --> p + Li7
fit(2.178990e+01,7.280980e-04,-3.025400e-01,-3.602000e-01,1.747200e-01,-2.230000e-02,-4.581000e-01),
#11. Li7 + p --> He4 + He4
fit(2.044380e+01,0.000000e+00,-8.472700e+00,2.979340e-01,5.823350e-02,-4.133830e-03,-6.666670e-01) + fit(2.189990e+01,-2.615270e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00) + fit(1.195760e+01,0.000000e+00,-8.472700e+00,4.179430e-01,5.345650e+00,-4.868400e+00,-6.666670e-01) + fit(1.425380e+01,-4.478000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00)
])
# Equation 40:
ForwardRate = (crosssection)*rhob#/0.6212711322669096
mass_excess = MeVtoK9(np.array([
2.22457,5.493, 3.2689 ,4.03266,0.7638, 17.589 , 18.353, 1.5861,2.467, 1.644, 17.347
]))
an, ap, ad, at, a3, a4, ali7, abe7 = np.array([mn, mp, md, mt, m3, m4, mli7, mbe7])/amu
conversion = 0.98677e10
prefactors = np.array([
# A = 1 for p & n, = 2 for D, = 3 for T & He3, = 4 for He4, = 7 for Be7 and Li7
# g = 2 for n,p, T, He3, g=3 for D, g=1 for He4, g=4 for Li7 and Be7.
# Below is (g)*(A)^3/2 * exp(-Q/T) and, in the case of photon emission, multiplied by the corrective factor in Eq. 42.
#1. p + n --> D + γ
conversion*((2*2)/(3) * ((an*ap)/(ad))**(3/2)* rhob**(-1) * T9**(3/2) * np.exp(-mass_excess[0]/(T9))),#/((2*np.pi)**1.5*(hbar*c)**3),#((2*2)/(3) * ((mn*mp)/(md))**(3/2)* rhob**(-1) * T**(3/2) * np.exp(-mass_excess[0]/(T)))/((2*np.pi)**1.5*(hbar*c)**3),
#2. D + p --> He3 + γ
conversion*((3*2)/(2) * ((ad*ap)/(a3))**(3/2)* rhob**(-1) * T9**(3/2)* np.exp(-mass_excess[1]/(T9))),
#3. D + D --> n + He3
(3*3)/(2*2) * ((ad*ad)/(an*a3))**(3/2)* np.exp(-mass_excess[2]/(T9)),
#4. D + D --> p + T
(3*3)/(2*2) * ((ad*ad)/(ap*at))**(3/2)* np.exp(-mass_excess[3]/(T9)),
#5. He3 + n --> p + T
(2*2)/(2*2) * ((a3*an)/(ap*at))**(3/2)* np.exp(-mass_excess[4] /(T9)),
#6. T + D --> n + He4
(2*3)/(2*1) * ((at*ad)/(an*a4))**(3/2)* np.exp(-mass_excess[5]/(T9)),
#7. He3 + D --> p + He4
(2*3)/(2*1) * ((a3*ad)/(ap*a4))**(3/2)* np.exp(-mass_excess[6]/(T9)),
#8. He3 + He4 --> Be7 + γ
conversion*((2*1)/(4) * ((a3*a4)/(abe7))**(3/2)* rhob**(-1) * T9**(3/2)* np.exp(-mass_excess[7]/(T9))),
#9. T + He4 --> Li7 + γ
conversion*((2*1)/(4) * ((at*a4)/(ali7))**(3/2)* rhob**(-1) * T9**(3/2)* np.exp(-mass_excess[8]/(T9))),
#10. Be7 + n --> p + Li7
(2*4)/(2*4) * ((an*abe7)/(ap*ali7))**(3/2)* np.exp(-mass_excess[9]/(T9)),
#11. Li7 + p --> He4 + He4
(4*1)/(1*1) * ((ali7*ap)/(a4*a4))**(3/2)* np.exp(-mass_excess[10]/(T9))
])
# Equations 41 and 42:
ReverseRate = prefactors*ForwardRate
return np.asarray([ForwardRate, ReverseRate])
ForwardRate, ReverseRate = StrongRates(T9, Tnu9)
Below is the code for the BBN network. We must input \(t\) for the ODE solver.
### REACTION NETWORK ###
# Allows us to call functions within this function.
interpolation = interpolate.interp1d(t, [T9, Tnu9, weak_np, weak_pn], 'cubic',fill_value='extrapolate')
# Equation 45:
def BBN(t,Y):
Yp, Yn, YD, YT, YHe3, YHe4, YBe7, YLi7 = Y
T9,Tnu9,np,pn = interpolation(t)
F, R = StrongRates(T9,Tnu9)
#1. p + n <--> D + γ
#2. D + p <--> He3 + γ
#3. D + D <--> n + He3
#4. D + D <--> p + T
#5. He3 + n <--> p + T
#6. T + D <--> n + He4
#7. He3 + D <--> p + He4
#8. He3 + He4 <--> Be7 + γ
#9. T + He4 <--> Li7 + γ
#10. Be7 + n <--> p + Li7
#11. Li7 + p <--> He4 + He4
dYpF = -pn*Yp + np*Yn - F[0]*Yp*Yn - F[1]*Yp*YD + (1/2)*F[3]*YD**2 + F[4]*YHe3*Yn + F[6]*YHe3*YD + F[9]*YBe7*Yn - F[10]*YLi7*Yp
dYnF = -np*Yn + pn*Yp - F[0]*Yp*Yn + (1/2) * F[2]*YD**2 - F[4]*YHe3*Yn + F[5]*YT*YD - F[9]*YBe7*Yn
dYDF = F[0]*Yp*Yn - F[1]*Yp*YD - F[2]*YD**2 - F[3]*YD**2 - F[5]*YT*YD - F[6]*YHe3*YD
dYTF = (1/2)* F[3]*YD**2 + F[4]*YHe3*Yn - F[5]*YT*YD - F[8]*YT*YHe4
dYHe3F = F[1]*Yp*YD + (1/2)*F[2]*YD**2 - F[4]*YHe3*Yn - F[6]*YHe3*YD - F[7]*YHe4*YHe3
dYHe4F = F[5]*YT*YD + F[6]*YHe3*YD - F[7]*YHe4*YHe3 - F[8]*YT*YHe4 + F[10]*YLi7*Yp
dYBe7F = F[7]*YHe4*YHe3 - F[9]*YBe7*Yn
dYLi7F = F[8]*YT*YHe4 + F[9]*YBe7*Yn - (1/2)*F[10]*YLi7*Yp
dYpR = R[0]*YD + R[1]*YHe3 - R[3]*Yp*YT - R[4]*Yp*YT - R[6]*YHe4*Yp - R[9]*YLi7*Yp + (1/2) * R[10]*YHe4**2
dYnR = R[0]*YD - R[2]*YHe3*Yn + R[4]*Yp*YT - R[5]*YHe4*Yn + R[9]*YLi7*Yp
dYDR = -R[0]*YD + R[1]*YHe3 + 2*R[2]*YHe3*Yn + 2*R[3]*Yp*YT + R[5]*YHe4*Yp + R[6]*YHe4*Yp
dYTR = - R[3]*Yp*YT - R[4]*Yp*YT + R[5]*YHe4*Yp + R[8]*YLi7
dYHe3R = -R[1]*YHe3 - R[2]*YHe3*Yn + R[4]*Yp*YT + R[6]*YHe4*Yp + R[7]*YBe7
dYHe4R = -R[5]*YHe4*Yn - R[6]*YHe4*Yp + R[7]*YBe7 + R[8]*YLi7 - R[10]*YHe4**2
dYBe7R = -R[7]*YBe7 + R[9]*YLi7*Yp
dYLi7R = -R[8]*YLi7 - R[9]*YLi7*Yp + (1/2) * R[10]*YHe4**2
dYp = dYpF + dYpR
dYn = dYnF + dYnR
dYD = dYDF + dYDR
dYT = dYTF+dYTR
dYHe3 = dYHe3F + dYHe3R
dYHe4 = dYHe4F + dYHe4R
dYBe7 = dYBe7F + dYBe7R
dYLi7 = dYLi7F + dYLi7R
return dYp, dYn, dYD, dYT, dYHe3, dYHe4, dYBe7, dYLi7
### SOLVING THE ODE ###
Y = solve_ivp(BBN, t[[0,-1]], X0, t_eval=t, dense_output=True, method='Radau',rtol=1e-14,atol=1e-15)
print(Y)
/home/aidan/micromamba/envs/physics/lib/python3.12/site-packages/scipy/integrate/_ivp/common.py:47: UserWarning: At least one element of `rtol` is too small. Setting `rtol = np.maximum(rtol, 2.220446049250313e-14)`.
warn("At least one element of `rtol` is too small. "
message: The solver successfully reached the end of the integration interval.
success: True
status: 0
t: [ 0.000e+00 2.022e-04 ... 1.285e+04 1.321e+04]
y: [[ 5.028e-01 5.327e-01 ... 7.565e-01 7.565e-01]
[ 4.972e-01 4.673e-01 ... 1.941e-09 1.838e-09]
...
[ 6.203e-62 3.354e-68 ... 3.583e-10 3.583e-10]
[ 6.118e-62 3.431e-68 ... 3.560e-11 3.562e-11]]
sol: <scipy.integrate._ivp.common.OdeSolution object at 0x78b8541954f0>
t_events: None
y_events: None
nfev: 109159
njev: 218
nlu: 3554
Show code cell source
weights = np.array([1,1,2,3,3,4,7,7])
#weights = np.array([mp/amu, mn/amu, md/amu, mt/amu, m3/amu, m4/amu, mbe7/amu, mli7/amu])
## Initialize the figure and axis
fig = plt.figure(figsize=(15,10), dpi=400)
ax = plt.subplot(111)
ax.grid()
#colors = ['darkblue', 'red', 'palevioletred', 'coral', 'mediumslateblue', 'darkolivegreen','darkturquoise', 'goldenrod', 'firebrick']
colors = ['black', '#e31a1c', 'cornflowerblue', 'sienna', '#33a02c', 'blueviolet', 'mediumvioletred', '#ff7f00', '#cab2d6']
## Plot
for n in range(8):
ax.loglog(T, (Y.y[n])*weights[n], c=colors[n])
## Legend
labels = [r'$\rm p \, (H)$', 'n', r'$\rm D \, (^2 H)$', r'$\rm T \, (^3 H)$', r'$\rm ^3He$', r'$\rm Y_p \, (^4 He)$', r'$\rm ^7Be$', r'$\rm ^7Li$']
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax.legend(labels, loc='center left', bbox_to_anchor=(1,0.5), fontsize=30)
## Tick and Ruling Parameters
from matplotlib.ticker import LogFormatterSciNotation, LogLocator
#ax.set_xlim(1e2,1e-1);
ax.set_xlim(1e0,1e-2)
ax.set_ylim(1e-12,2e0)
ax.tick_params(axis='both', direction='in', which='major', length = 20, width=1.15, right=True, labelsize=36, pad=8)
ax.tick_params(axis='x', direction='in', length=10, width=0.8, right=True, which='minor')
ax.tick_params(axis='y', direction='in', length=20, width=0.6, right=True, which='minor', labelsize=15)
ax.xaxis.set_minor_locator(LogLocator())
ax.yaxis.set_minor_locator(LogLocator())
ax.xaxis.set_minor_formatter(LogFormatterSciNotation())
ax.yaxis.set_minor_formatter(LogFormatterSciNotation())
ax.minorticks_on()
## Labels
#ax.set_title('Mass Fraction of Elements in BBN')
ax.set_ylabel('Mass Fraction', fontsize=36)
ax.set_xlabel(r'Temperature ($\rm{10^9 K}$)', fontsize=36);
textlabels = ['p','n','D','T','3','4','Be','Li']
print(f'The final mass fraction of Hydrogen (protons) is {weights[0]*(Y.y[0][-1]) : 0.4g}')
print(f'The final mass fraction of Helium 4 (Yp) is {(Y.y[5][-1]) *weights[5]: 0.4g}')
print(f'The final percentage of mass fraction of the sum Li+Be is',f' {(Y.y[6][-1] + Y.y[7][-1] )/Y.y[0][-1] : 0.4g}')
print('----------------------------- \n')
for n in [1,2,3,4,6,7]:
print(f'The final percentage of mass fraction of',textlabels[n],f'is {(Y.y[n][-1])/Y.y[0][-1]:0.4g}')
sum = 0
for n in range(8):
sum += Y.y[n][-1]*weights[n]
print(f'The total sum of elements is: {sum*100:0.0f} %.')
The final mass fraction of Hydrogen (protons) is 0.7565
The final mass fraction of Helium 4 (Yp) is 0.2434
The final percentage of mass fraction of the sum Li+Be is 5.207e-10
-----------------------------
The final percentage of mass fraction of n is 2.43e-09
The final percentage of mass fraction of D is 2.152e-05
The final percentage of mass fraction of T is 7.337e-08
The final percentage of mass fraction of 3 is 9.453e-06
The final percentage of mass fraction of Be is 4.736e-10
The final percentage of mass fraction of Li is 4.709e-11
The total sum of elements is: 100 %.
