BBN 0.0.1

Big Bang Nucleosynthesis

BBN: Big Bang Nucleosynthesis

BBN(eta, index, n_step=256, **kw):
input:
  eta = baryon to photon ratio (Kolb and Turner eq.3.104)
  index = list of element's name (string)
          (available indices are given in nuclear.py)
  n_step = number of steps at which data are evaluated
           with logarithmic interval in T
  kw = keyword arguments passed to solve_ivp, e.g.
    atol = tolerance for absolute error (default 1e-6)
    rtol = tolerance for relative error (default 1e-3)
return T,X where
  T = temperature / MeV (shape(n_step,)) in [T_init, T_final]
  X = mass fraction of elements (shape(len(index),n_step))
reference:
  E. W. Kolb and M. S. Turner
    "The Early Universe" chapter 4
-------------------------------------------------------------

initialize(T_init=1e1, T_final=1e-2, N_nu=3, tau_n=885.7):
input:
  T_init,T_final = temperature / MeV
  N_nu = number of neutrino generation
  tau_n = neutron mean lifetime / sec

example code:

import matplotlib.pyplot as plt
from BBN import BBN

# all available indices
index = ['neutron', 'proton', 'deutron', 'tritium',
         'helium3', 'helium4', 'lithium7', 'beryllium7']
label = ['n', 'p', 'd', 't', r'He$^3$',
         r'He$^4$', r'Li$^7$', r'Be$^7$']

# initialize() has been executed with default arguments
T,X = BBN(5e-10, index, atol=1e-13)
plt.axis([2, 1e-2, 1e-13, 2])
plt.loglog(T, X.T)
plt.xlabel('T = temperature / MeV')
plt.ylabel('X = mass fraction')
plt.legend(label)
plt.show()

example code:

import numpy as np
import matplotlib.pyplot as plt
from BBN import BBN,initialize

index = ['helium4']
label = [r'$N_{\nu}$ = 2', '3', '4']

eta = np.geomspace(3e-11, 1e-8, 50)
X = []
for N_nu in [2,3,4]:# number of neutrino generagion
    initialize(N_nu=N_nu)
    X1 = []
    for e in eta:
        T,X2 = BBN(e, index, 2, rtol=1e-6, atol=1e-9)
        X1.append(X2[0,-1])
        print(e, X1[-1])

    X.append(X1)

plt.axis([eta[0], eta[-1], 0.17, 0.27])
plt.semilogx(eta, np.asarray(X).T)
plt.xlabel(r'$\eta$ = baryon to photon ratio')
plt.ylabel(r'$X_4$ = mass fraction of He$^4$')
plt.legend(label, markerfirst=False)
plt.show()