Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

CMB BBN

Authors
Affiliations
Université Paris-Saclay
IJCLab
Université Paris-Saclay
IJCLab
Université Paris-Saclay
IJCLab 
import numpy as np
import matplotlib.pyplot as plt
import astropy.constants as const
import astropy.units as u
from scipy.special import zeta
from scipy.interpolate import interp1d
from astropy.cosmology import Planck18

1Neutron freeze-out

Γn≈H\Gamma_n \approx H
(1)
Γn=σwnνc\Gamma_n = \sigma_w n_\nu c
(2)
nν(Tγ)=34411nγ(Tγ)(onlyνe)n_\nu(T_\gamma) = \frac{3}{4}\frac{4}{11}n_\gamma(T_\gamma)\quad (only \nu_e)
(3)
nγ=2ζ(3)(kBT)3π2ℏ3c3n_\gamma = \frac{2 \zeta(3) (k_B T)^3}{\pi^2 \hbar^3 c^3}
(4)
σw=(GF2E2)/(ℏc)2=GF2(kBT/(ℏc))2=5×10−48(kBT1 MeV)2 m2\sigma_w = (G_F^2 E^2)/(\hbar c)^2 = G_F^2 (k_B T/(\hbar c))^2 = 5\times 10^{-48} \left(\frac{k_B T}{1\,\mathrm{MeV}}\right)^2\,\mathrm{m}^2
(5)
GF=1.166×10−5 GeV−2G_F = 1.166 \times 10^{-5}\,\text{GeV}^{-2}
(6)
gA=1.26(Kolb Turner p.91)g_A = 1.26 \quad\text{(Kolb Turner p.91)}
(7)
H=8π3G90ℏ3c5g⋆1/2(T)(kBT)2H = \sqrt{\frac{8 \pi^3 G}{90 \hbar^3 c^5}} g_\star^{1/2}(T) (k_B T)^2
(8)
g⋆(100 MeV>T>511 keV)=10.75g_\star(100\,\mathrm{MeV} > T > 511\,\mathrm{keV}) = 10.75
(9)
g⋆(T<511 keV)=3.36g_\star(T < 511\,\mathrm{keV}) = 3.36
(10)
T0 = 2.726 * u.K  # K
Qn = 1.29 * u.MeV
Neff = 3.046
gstar = 10.75  # at m_e < T < 100 MeV
gstar_after = 3.36  # at m_e > T
GF = 1.166e-5 / (u.GeV**2)
gA = 1.26
taun = 878.4 * u.s
def ng(T):
    return 2 * zeta(3) * (const.k_B * T)**3 / (np.pi**2 * const.hbar**3 * const.c**3)

ng(T0).to(1/u.cm**3)
Loading...
# 3/4 for fermion nature, only 1 of the 3 flavors
# then 2 polarisations for photons vs g=2 for massless neutrinos
def nnu_e(T):
    return 4*ng(T)/11 * (3/4) #* Neff

nnu_e(T0).to(1/u.cm**3)
Loading...
def H(T, gstar=gstar):
    return (np.sqrt(8*np.pi**3*const.G*gstar/(90*const.hbar**3*const.c**5))*(const.k_B * T)**2).to(1/u.s)

T_1MeV = (1 * u.MeV / const.k_B).to(u.K)
T_10MeV = (10 * u.MeV / const.k_B).to(u.K)
H(T_1MeV)  # OK with Kolb Turner formula p.91
Loading...
def t_to_T(t, gstar=gstar):
    return (1/const.k_B*np.sqrt(1/(2*t)*np.sqrt(90*const.hbar**3*const.c**5/(8*np.pi**3*const.G*gstar) ))).to(u.K)

def T_to_t(T, gstar=gstar):
    return 1/(2*H(T, gstar=gstar))

t_to_T(10*u.s), T_to_t(t_to_T(10*u.s)), t_to_T(1*u.s)
(<Quantity 3.15287605e+09 K>, <Quantity 10. s>, <Quantity 9.97026951e+09 K>)

2Neutrino temperature and time to temperature accurate conversion

See Weinberg p.540 p.529

T_Weinberg = np.array([1e12, 6e11, 3e11, 2e11, 1e11, 6e10, 3e10, 2e10, 1e10, 6e9, 3e9, 2e9, 1e9, 3e8, 1e8, 1e7, 1e6, 1e5, 1e4, 4e3])
a_Weinberg = np.array([1.9e-12, 3.2e-12, 6.4e-12, 9.6e-12, 1.9e-11, 3.2e-11, 6.4e-11, 9.6e-11, 1.9e-10, 3.1e-10, 5.9e-10, 8.3e-10, 2.6e-9, 9e-9, 2.7e-8, 2.7e-7, 2.7e-6, 2.7e-5, 2.7e-4, 6.3e-4])
Tnu_Weinberg = T_Weinberg / np.array([1, 1, 1, 1, 1, 1, 1.001, 1.002, 1.008, 1.022, 1.081, 1.159, 1.346, 1.401, 1.401,1.401,1.401,1.401,1.401,1.401])
t_Weinberg = np.array([1e-4, 1.94e-4, 1.129e-3, 2.61e-3, 1.078e-2, 3.01e-2, 0.1209, 0.273, 1.103, 3.14, 13.83, 35.2, 1.82e2, 2.08e3, 1.92e4, 1.92e6,1.92e8,1.92e10,1.92e12,1.2e13])
fig = plt.figure()
plt.loglog(t_Weinberg, T_Weinberg, lw=2, label=r"$T_\gamma$")
plt.loglog(t_Weinberg, Tnu_Weinberg, lw=2, color="r", label=r"$T_\nu$")
plt.xlim(1e-1, 1e4)
plt.ylim(1e8, 1e11)
plt.axvline(2, linestyle="--", color="k")
plt.axvline(2e3, linestyle="--", color="k")
plt.legend(fontsize=14)
plt.xlabel("t [s]", fontsize=16)
plt.ylabel("Temperature [K]", fontsize=16)
plt.text(4e0, 5e10, "Era of $e^+ e^-$ annihilation", fontsize=14)
plt.grid()
plt.show()
<Figure size 640x480 with 1 Axes>
def T_to_t_Weinberg(T):
    return np.exp(np.interp(np.log(T), np.log(T_Weinberg[::-1]), np.log(t_Weinberg[::-1])))*u.s

def t_to_T_Weinberg(t):
    return np.exp(np.interp(np.log(t), np.log(t_Weinberg), np.log(T_Weinberg)))*u.K

ttt = np.logspace(-4, 5, 200) * u.s

fig = plt.figure()
plt.loglog(ttt, t_to_T(ttt, gstar=gstar), label="T with $g_*$ before $m_e$")
plt.loglog(ttt, t_to_T(ttt, gstar=gstar_after), label="T with $g_*$ after $m_e$")
plt.loglog(ttt, t_to_T_Weinberg(ttt.value), "k-", lw=2, label="Exact from Weinberg")
plt.grid()
plt.legend(fontsize=14)
plt.xlabel("t [s]", fontsize=16)
plt.ylabel("Temperature [K]", fontsize=16)
plt.show()
<Figure size 640x480 with 1 Axes>
TTT = (np.logspace(3, 8, 50) * u.eV/ const.k_B).to(u.K)

fig = plt.figure()
plt.loglog(TTT, T_to_t(TTT, gstar=gstar), label="t with $g_*$ before $m_e$")
plt.loglog(TTT, T_to_t(TTT, gstar=gstar_after), label="t with $g_*$ after $m_e$")
plt.loglog(TTT, T_to_t_Weinberg(TTT.value), "k-", lw=2, label="Exact from Weinberg")
plt.grid()
plt.legend(fontsize=14)
plt.ylabel("t [s]", fontsize=16)
plt.xlabel("Temperature [K]", fontsize=16)
plt.show()
<Figure size 640x480 with 1 Axes>
T_fine = np.logspace(np.log10(np.min(T_Weinberg)), np.log10(np.max(T_Weinberg)), 200)
t_fine = np.logspace(np.log10(np.min(t_Weinberg)), np.log10(np.max(t_Weinberg)), 200)
a_fine = np.logspace(np.log10(np.min(a_Weinberg)), np.log10(np.max(a_Weinberg)), 200)
H_fine = np.gradient(a_fine, t_fine) / a_fine

def Ht_Weinberg(t):
    return np.interp(t, t_fine, H_fine) / u.s

def HT_Weinberg(T):
    return np.interp(T, T_fine, H_fine) / u.s

fig = plt.figure()
plt.loglog(ttt, Ht_Weinberg(ttt.value), label="Exact from Weinberg")
plt.loglog(ttt, 1/(2*ttt), "r--", label="$1/(2t)$")
plt.loglog(ttt, H(t_to_T_Weinberg(ttt.value), gstar=gstar), "b--", label="$g_*$ before $m_e$")
plt.loglog(ttt, H(t_to_T_Weinberg(ttt.value), gstar=gstar_after), "k--", label="$g_*$ after $m_e$")
plt.grid()
plt.legend(fontsize=14)
plt.xlabel("t [s]", fontsize=16)
plt.ylabel("H(t) [1/s]", fontsize=16)
plt.show()
<Figure size 640x480 with 1 Axes>

3Orders of magnitude

From equilibrium equations and rate with GFG_F and gAg_A :

# at 1 MeV
sigma_np = (((GF**2)  * (1+3*gA**2)  * 1*u.MeV**2).to(1/u.GeV**2)* (const.hbar * const.c)**2).to(u.m**2)
sigma_np
Loading...
def Gamma_n(T):
    return (sigma_np * (const.k_B * T / (1*u.MeV))**2 * nnu_e(T) * const.c).to(1/u.s)

1/Gamma_n(0.8*T_1MeV)
Loading...

Freeze out temperature

from scipy.optimize import fsolve
Tfreeze = fsolve(lambda T: Gamma_n(T*u.K)-H(T*u.K), [1e8, 1e12])[1]
(Tfreeze*u.K * const.k_B).to(u.MeV)
Loading...
TT = (np.logspace(4, 7, 50) * u.eV/ const.k_B).to(u.K)

fig = plt.figure()
plt.plot(TT, Gamma_n(TT), label="$\Gamma_n(T)$")
plt.plot(TT, H(TT), label="$H(T)$")
plt.gca().invert_xaxis()
plt.yscale("log")
plt.xscale("log")
plt.legend()
plt.xlabel("$T$ [K]")
plt.show()
<Figure size 640x480 with 1 Axes>

4Exact computation

The freeze-out temperature obtained with orders of magnitude is too high to explain the helium aboundance (by a factor of 3, but here numerical factors do count). We now use exact (Weinberg integrals) and semi-exact (Dodelson approximation) methods.

Using Weinberg and Dodelson equations

Check equality of constant gVg_V in Weinberg eq 15.7.9 p.547 with with Fermi constant GFG_F

gV = 1.418e-49*u.erg*u.cm**3
(gV/(const.hbar*const.c)**3).to(u.GeV**-2)
Loading...
GF
Loading...
A=(((gV/(const.hbar*const.c)**3)**2+3*(gA*gV/(const.hbar*const.c)**3)**2)/(2*np.pi**3)).to(u.GeV**-4)
A #*7/15*np.pi**4
Loading...

In Weinberg, neutron life time is 1013s but current value if 878s. We rescale its amplitude constant:

# rescale with current neutron life as Weinberg
A = GF**2*(1+3*gA**2)/(2*np.pi**3) * 1013/878
A
Loading...

Computation of interaction rates:

from scipy.integrate import quad
me = 511*u.keV

integrand_lpn = lambda q, T: (A*np.sqrt(1-(me/(Qn+q))**2)*(Qn+q)**2*(q)**2/(1+np.exp(-q/(const.k_B*T)))/(1+np.exp((Qn+q)/(const.k_B * T)))).to(u.dimensionless_unscaled)

integrand_lnp = lambda q, T: (A*np.sqrt(1-(me/(Qn+q))**2)*(Qn+q)**2*(q)**2/(1+np.exp(q/(const.k_B*T)))/(1+np.exp(-(Qn+q)/(const.k_B * T)))).to(u.dimensionless_unscaled)

def lnp(T):
    intp = quad(lambda q: integrand_lnp(q*u.MeV,T), (-Qn+me).to(u.MeV).value, np.inf)[0]
    intm = quad(lambda q: integrand_lnp(q*u.MeV,T), -np.inf, (-Qn-me).to(u.MeV).value)[0]
    return ((intp+intm)*u.MeV/const.hbar).to(1/u.s)

def lpn(T):
    intp = quad(lambda q: integrand_lpn(q*u.MeV,T), (-Qn+me).to(u.MeV).value, np.inf)[0]
    intm = quad(lambda q: integrand_lpn(q*u.MeV,T), -np.inf, (-Qn-me).to(u.MeV).value)[0]
    return ((intp+intm)*u.MeV/const.hbar).to(1/u.s)

def lambda_np_Bernstein(T):  # Dodeslon p.67
    x = Qn / (const.k_B * T)
    return (255/(taun*x**5)*(12+6*x+x**2)).to(1/u.s)

lnp(0.8*T_1MeV), lpn(0.8*T_1MeV) , lambda_np_Bernstein(0.8*T_1MeV), H(0.8*T_1MeV)
(<Quantity 0.65549597 1 / s>, <Quantity 0.13069837 1 / s>, <Quantity 0.64641477 1 / s>, <Quantity 0.4335013 1 / s>)

Check that λpn∼λnpe−Qn/T\lambda_{pn} \sim \lambda_{np} e^{-Qn/T} :

lnp(0.8*T_1MeV) * np.exp(-Qn/(const.k_B * 0.8 * T_1MeV))
Loading...
TT = (np.logspace(5, 7, 50) * u.eV/ const.k_B).to(u.K)

fig = plt.figure()
plt.plot(TT, np.array([lnp(T).value for T in TT])/u.s)
plt.plot(TT, lambda_np_Bernstein(TT))
plt.yscale("log")
plt.xscale("log")
plt.show()
/Users/jneveu/miniforge3/envs/m2-cosmo/lib/python3.11/site-packages/astropy/units/quantity.py:671: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
<Figure size 640x480 with 1 Axes>

Solve

from scipy.integrate import odeint

def dXndt_notaun(Xn, t):
    return (lambda_np_Bernstein(t_to_T(t*u.s))*( (1-Xn)*np.exp(-Qn/(const.k_B*t_to_T(t*u.s))) - Xn)).value

def dXndt(Xn, t):
    # return (lambda_np_Bernstein(t_to_T(t*u.s))*( (1-Xn)*np.exp(-Qn/(const.k_B*t_to_T(t*u.s))) - Xn) - Xn/taun).value
    # return (lnp(t_to_T(t*u.s))*(-Xn + np.exp(-Qn/(const.k_B * t_to_T(t*u.s))) * (1-Xn))).value
    return (-Xn * lnp(t_to_T(t*u.s)) + lpn(t_to_T(t*u.s)) * (1-Xn)).value
    # return (-Xn * lnp(t_to_T_Weinberg(t)) + lpn(t_to_T_Weinberg(t)) * (1-Xn)).value  # slower

tt = np.logspace(-3, 3, 200)

Xn_sol = odeint(dXndt, 1/(1+np.exp(Qn/(const.k_B * t_to_T(tt[0]*u.s)))), tt)
Xnt = lambda t: np.interp(t, tt*u.s, Xn_sol[:, 0])
XnT = lambda T: np.interp(T, t_to_T(tt[::-1]*u.s), Xn_sol[:, 0][::-1])

Xn_sol_notaun = odeint(dXndt_notaun, 1/(1+np.exp(Qn/(const.k_B * t_to_T(tt[0]*u.s)))), tt)
Xnt_notaun = lambda t: np.interp(t, tt*u.s, Xn_sol_notaun[:, 0])
XnT_notaun = lambda T: np.interp(T, t_to_T(tt[::-1]*u.s), Xn_sol_notaun[:, 0][::-1])

XnT(T_1MeV), Xnt(0.73*u.s)
(np.float64(0.2546931106785803), np.float64(0.2554279084910855))
dXndt(XnT(1e10*u.K), T_to_t(1e10*u.K).value), XnT(1e10*u.K), H(1e10*u.K),1/(2*T_to_t(1e10*u.K)), T_to_t(1e10*u.K)
(np.float64(-0.059484795266252305), np.float64(0.23618253503009692), <Quantity 0.50298636 1 / s>, <Quantity 0.50298636 1 / s>, <Quantity 0.99406274 s>)
-dXndt(Xnt(1*u.s), 1), Xnt(1*u.s), H(t_to_T(1*u.s)), 1/(2*1*u.s), (const.k_B*t_to_T(1*u.s)).to(u.MeV)
(np.float64(0.0590577858660713), np.float64(0.23588394768970708), <Quantity 0.5 1 / s>, <Quantity 0.5 1 / s>, <Quantity 0.85917135 MeV>)

Differential equation 3.27 from Dodelson p.67

def dXndx(Xn, x):
    T = (Qn / (const.k_B * x)).to(u.K)
    return (x*lnp(T)/H(Qn/const.k_B)*(np.exp(-x) - Xn*(1+np.exp(-x)))).value

xx = np.logspace(-3, 3, 200)
XnDodelson_sol = odeint(dXndx, 1/(1+np.exp(xx[0])), xx)
XnDodelson_sol
array([[ 4.99750000e-01], [ 4.99732027e-01], [ 4.99712762e-01], [ 4.99692113e-01], [ 4.99669978e-01], [ 4.99646252e-01], [ 4.99620821e-01], [ 4.99593561e-01], [ 4.99564342e-01], [ 4.99533022e-01], [ 4.99499450e-01], [ 4.99463465e-01], [ 4.99424893e-01], [ 4.99383548e-01], [ 4.99339230e-01], [ 4.99291727e-01], [ 4.99240808e-01], [ 4.99186229e-01], [ 4.99127727e-01], [ 4.99065019e-01], [ 4.98997802e-01], [ 4.98925754e-01], [ 4.98848526e-01], [ 4.98765746e-01], [ 4.98677015e-01], [ 4.98581904e-01], [ 4.98479956e-01], [ 4.98370679e-01], [ 4.98253546e-01], [ 4.98127992e-01], [ 4.97993412e-01], [ 4.97849157e-01], [ 4.97694530e-01], [ 4.97528788e-01], [ 4.97351130e-01], [ 4.97160700e-01], [ 4.96956582e-01], [ 4.96737795e-01], [ 4.96503279e-01], [ 4.96251904e-01], [ 4.95982457e-01], [ 4.95693653e-01], [ 4.95384073e-01], [ 4.95052240e-01], [ 4.94696554e-01], [ 4.94315301e-01], [ 4.93906642e-01], [ 4.93468609e-01], [ 4.92999091e-01], [ 4.92495824e-01], [ 4.91956384e-01], [ 4.91378172e-01], [ 4.90758401e-01], [ 4.90094834e-01], [ 4.89382962e-01], [ 4.88619959e-01], [ 4.87802167e-01], [ 4.86925653e-01], [ 4.85986224e-01], [ 4.84979408e-01], [ 4.83900349e-01], [ 4.82743902e-01], [ 4.81504556e-01], [ 4.80176432e-01], [ 4.78753223e-01], [ 4.77228200e-01], [ 4.75594190e-01], [ 4.73843543e-01], [ 4.71968118e-01], [ 4.69959279e-01], [ 4.67807812e-01], [ 4.65503922e-01], [ 4.63037146e-01], [ 4.60396312e-01], [ 4.57569463e-01], [ 4.54543946e-01], [ 4.51306504e-01], [ 4.47843362e-01], [ 4.44140352e-01], [ 4.40182994e-01], [ 4.35956664e-01], [ 4.31446691e-01], [ 4.26638529e-01], [ 4.21518144e-01], [ 4.16071925e-01], [ 4.10287074e-01], [ 4.04151718e-01], [ 3.97655330e-01], [ 3.90789577e-01], [ 3.83549444e-01], [ 3.75934762e-01], [ 3.67951811e-01], [ 3.59614629e-01], [ 3.50945697e-01], [ 3.41976128e-01], [ 3.32740413e-01], [ 3.23297666e-01], [ 3.13690719e-01], [ 3.03979764e-01], [ 2.94228780e-01], [ 2.84505817e-01], [ 2.74881745e-01], [ 2.65428386e-01], [ 2.56215748e-01], [ 2.47309246e-01], [ 2.38764471e-01], [ 2.30636395e-01], [ 2.22961146e-01], [ 2.15767659e-01], [ 2.09074573e-01], [ 2.02890716e-01], [ 1.97215665e-01], [ 1.92040636e-01], [ 1.87349560e-01], [ 1.83120195e-01], [ 1.79325456e-01], [ 1.75934664e-01], [ 1.72914673e-01], [ 1.70230901e-01], [ 1.67848201e-01], [ 1.65731520e-01], [ 1.63846538e-01], [ 1.62160022e-01], [ 1.60640093e-01], [ 1.59256354e-01], [ 1.57979885e-01], [ 1.56783173e-01], [ 1.55639946e-01], [ 1.54524953e-01], [ 1.53413694e-01], [ 1.52282148e-01], [ 1.51106475e-01], [ 1.49862703e-01], [ 1.48526470e-01], [ 1.47072748e-01], [ 1.45475603e-01], [ 1.43708014e-01], [ 1.41741733e-01], [ 1.39547231e-01], [ 1.37093733e-01], [ 1.34349408e-01], [ 1.31281738e-01], [ 1.27858099e-01], [ 1.24046626e-01], [ 1.19817430e-01], [ 1.15144192e-01], [ 1.10006209e-01], [ 1.04390877e-01], [ 9.82966449e-02], [ 9.17363196e-02], [ 8.47405662e-02], [ 7.73613310e-02], [ 6.96747364e-02], [ 6.17829304e-02], [ 5.38140758e-02], [ 4.59198049e-02], [ 3.82693599e-02], [ 3.10400002e-02], [ 2.44037408e-02], [ 1.85113633e-02], [ 1.34756805e-02], [ 9.35687663e-03], [ 6.15346978e-03], [ 3.80192508e-03], [ 2.18646731e-03], [ 1.15797389e-03], [ 5.57878755e-04], [ 2.41069349e-04], [ 9.19319884e-05], [ 3.03683527e-05], [ 8.50525346e-06], [ 1.97021569e-06], [ 3.66736100e-07], [ 5.34917859e-08], [ 5.72355380e-09], [ 2.79701081e-10], [-4.65247228e-10], [-5.35028485e-10], [-2.08088554e-10], [-1.81120963e-10], [-1.52214636e-10], [-1.21230195e-10], [-8.80182415e-11], [-5.24186372e-11], [-1.42597303e-11], [-1.40219498e-11], [-2.82188292e-11], [-4.34363410e-11], [-5.97478597e-11], [-7.02529061e-11], [-6.79678413e-11], [-6.55185003e-11], [-6.28930728e-11], [-6.00788999e-11], [-5.70624125e-11], [-5.38290657e-11], [-5.03632694e-11], [-4.66483124e-11], [-4.26662823e-11], [-3.83979787e-11]])
XnDodelson_sol[-1]
array([-3.83979787e-11])
fig = plt.figure()
#plt.plot((TT*const.k_B).to(u.MeV), XnT(TT))
plt.plot((const.k_B*t_to_T(tt*u.s)).to(u.MeV), Xn_sol)
plt.plot((Qn / xx).to(u.MeV), XnDodelson_sol, '--')
#plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.grid()
plt.show()
<Figure size 640x480 with 1 Axes>


fig = plt.figure()
#plt.plot((TT*const.k_B).to(u.MeV), XnT(TT))
plt.plot((Qn / xx).to(u.MeV), [-dXndx(XnDodelson_sol[k], x) for k,x in enumerate(xx)])
plt.plot((Qn / xx).to(u.MeV), H((Qn / xx / const.k_B).to(u.K)), label="$H(T)$")
plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.grid()
plt.show()
<Figure size 640x480 with 1 Axes>
fig = plt.figure()
plt.plot((t_to_T(tt*u.s)*const.k_B).to(u.MeV), Xn_sol)
#plt.plot((t_to_T(tt*u.s)*const.k_B).to(u.MeV), XnDodelson_sol)
#plt.plot((Qn / xx).to(u.MeV), XnDodelson_sol)
#plt.plot((t_to_T(tt*u.s)*const.k_B).to(u.MeV), XnT(t_to_T(tt*u.s)))
#plt.plot((TT*const.k_B).to(u.MeV), XnT(TT))
#plt.plot((TT*const.k_B).to(u.MeV), Xnt(T_to_t(TT)))
#plt.plot((Qn / xx).to(u.MeV), Xn_sol)
#plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.grid()
plt.xlabel("$T$ [K]", fontsize=16)
plt.ylabel("$X_n$", fontsize=16)

plt.show()
<Figure size 640x480 with 1 Axes>
TTT = t_to_T(tt*u.s)

fig = plt.figure()
#plt.plot(tt, np.array([-dXndt(Xnt(t*u.s), t) for t in tt])/u.s, label="$\Gamma_n(T)$")
plt.plot(tt, np.array([lnp(T).value for T in TTT])/u.s, label="$\Gamma_n(T)$ exact")
plt.plot(tt, lambda_np_Bernstein(TTT), label="$\Gamma_n(T)$ no free neutron decay")
#plt.plot(tt, np.array([lnp(T).value+lpn(T).value for T in TTT])/u.s, label="$\Gamma_n(T)$")
plt.plot(tt, 1/(2*tt), label="$H(T)$")
plt.plot(tt, Ht_Weinberg(tt), "--", label="$H(t)$ Weinberg")
#plt.gca().invert_xaxis()
plt.axvline(1.7)
plt.yscale("log")
plt.xscale("log")
plt.legend()
plt.grid()
plt.xlabel("$t$ [s]")
plt.show()
<Figure size 640x480 with 1 Axes>
lnp(T_1MeV), T_1MeV
(<Quantity 1.77604491 1 / s>, <Quantity 1.16045181e+10 K>)
from scipy.optimize import fsolve
Tfreeze = fsolve(lambda T: lambda_np_Bernstein(T*u.K).value-H(T*u.K).value, [1e8, 1e12])[1] * u.K
#Tfreeze = fsolve(lambda T: lnp(T*u.K).value-H(T*u.K).value, [1e8, 1e12])[1] * u.K
(Tfreeze * const.k_B).to(u.MeV), Tfreeze
(<Quantity 0.67481545 MeV>, <Quantity 7.83090811e+09 K>)
Tfreeze = t_to_T_Weinberg(fsolve(lambda t: lambda_np_Bernstein(t_to_T_Weinberg(t)).value-Ht_Weinberg(t).value, [1e-2, 1e2])[0])
(Tfreeze * const.k_B).to(u.MeV), Tfreeze
(<Quantity 0.64293845 MeV>, <Quantity 7.46099084e+09 K>)
T_to_t(Tfreeze)
Loading...
def ratio_np(T):
    return XnT(T) / (1-XnT(T))

def ratio_np_notaun(T):
    return XnT_notaun(T) / (1-XnT_notaun(T))

npfreeze = ratio_np(Tfreeze)
npfreeze
np.float64(0.25901803917730154)

5Deutérium ratio

BD = 2.22 * u.MeV
mn = 939.6 * u.MeV
Omegab0 = Planck18.Ob0
H0 = Planck18.H0
age = Planck18.lookback_time(np.inf)
rhoc0 = (3 * Planck18.H0**2 / (8 * np.pi * const.G) / const.m_p).to(1/u.m**3)
def nb(z):
    return Omegab0 * rhoc0 * (1+z)**3

eta = nb(0) / ng(T0)
eta
Loading...
def ratio_Dn(T, eta=eta):
    # return (6*0.2*np.pi**3/2*eta*(const.k_B * T/(mn*const.c**2))**(3/2)*np.exp(BD/(const.k_B*T))).to(u.dimensionless_unscaled)
    return (6*(1-XnT(T))*2*zeta(3)*(np.pi/2)*eta*(const.k_B * T/mn)**(3/2)*np.exp(BD/(const.k_B*T))).to(u.dimensionless_unscaled)

ratio_Dn(Tfreeze)
Loading...
from scipy.optimize import brentq
Tnuc = brentq(lambda T: ratio_Dn(T*u.K, eta=eta)-1, 1e8, 1e10) * u.K
Tnuc
Loading...
(const.k_B * Tnuc).to(u.MeV)
Loading...
tnuc = T_to_t(Tnuc, gstar=gstar_after)
tnuc = T_to_t_Weinberg(Tnuc.value)
tnuc
Loading...
T_to_t(0.1*u.MeV/const.k_B, gstar=gstar_after), T_to_t_Weinberg((0.1*u.MeV/const.k_B).to(u.K).value)
(<Quantity 132.03654981 s>, <Quantity 127.90464621 s>)
npnuc = ratio_np(Tnuc)
Xnnuc = XnT(Tnuc)
npfreeze, npnuc, Xnnuc
(np.float64(0.25901803917730154), np.float64(0.14734740090432036), np.float64(0.12842439943489092))

Maximum helium 4 abondance (if all neutrons end in helium nuclei, no D nor Li nor...):

def Yp(np):  # Kolb Turner p.95
    return 2*np/(1+np)

def Yp_eta(eta):
    Tnuc = brentq(lambda T: ratio_Dn(T*u.K, eta)-1, 1e8, 1e10) * u.K
    npnuc = ratio_np(Tnuc)
    return 2*npnuc/(1+npnuc)

Yp(npnuc), Yp_eta(eta)
(np.float64(0.25684879886978185), np.float64(0.25684879886978185))

6Plots

TT = (np.logspace(np.log10(7e4), np.log10(5e7), 100) * u.eV/ const.k_B).to(u.K)
TTafter = (np.logspace(np.log10(7e4), np.log10(5e6), 100) * u.eV/ const.k_B).to(u.K)
npeq = np.exp(-Qn/(const.k_B * TTafter))

fig = plt.figure()
plt.plot((TT*const.k_B).to(u.MeV), XnT(TT), label=rf"exact with $\tau_n={taun:.1f}$", lw=2)
plt.plot((TT*const.k_B).to(u.MeV), XnT_notaun(TT), linestyle="-.", label="approximation without\nfree neutron decay", zorder=-2, color="k")
plt.plot((TTafter*const.k_B).to(u.MeV), npeq/(1+npeq), linestyle="dotted", label=r"$X_n^{eq}(T)$", zorder=-3, color="k")
#plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.xlabel("$T$ [MeV]", fontsize=16)
plt.ylabel("$X_n(T)$", fontsize=16)
plt.axvline((Tfreeze * const.k_B).to(u.MeV).value, color="b", linestyle="--", label="$T_\mathrm{freeze}$")
plt.axvline((Tnuc * const.k_B).to(u.MeV).value, color="r", linestyle="--", label="$T_\mathrm{nuc}$")
plt.legend(fontsize=14, ncol=1)
plt.axvspan((Tnuc * const.k_B).to(u.MeV).value, 5e-2, alpha=0.2, color="gray")  #horizontal shading
plt.xlim((50, 5e-2))
plt.ylim(0.0, 0.52)
plt.grid()

secax = plt.gca().twiny()
ttt = T_to_t_Weinberg(TT.value)
secax.plot(ttt, Xnt(ttt), linestyle="none")
secax.set_xscale("log")
secax.set_xlabel('$t$ [s]', fontsize=16)
fig.tight_layout()
plt.show()
<Figure size 640x480 with 2 Axes>
def ratio_mass_nb(T):
    return ratio_np(T)/(1+ratio_np(T))

def ratio_mass_Db(T):
    return 2*ratio_Dn(T)/(1+ratio_np(T))

def ratio_mass_Hb(T):
    if T > Tnuc:
        return 1/(1+ratio_np(T))
    else:
        return 1/(1+ratio_np(Tnuc)) - Yp(npnuc)/2  # half the mass of helium must be subtracted (2 protons)

def ratio_mass_Heb(T):
    if T > Tnuc:
        return 0
    else:
        return Yp(npnuc)
fig = plt.figure()
plt.plot((TT*const.k_B).to(u.MeV), ratio_mass_nb(TT), label="n")
p = plt.plot((TT*const.k_B).to(u.MeV), [ratio_mass_Hb(T) for T in TT], linestyle="-", label="H")
TTHe = (np.logspace(4, np.log10(3*(const.k_B * Tnuc).to(u.eV).value), 20) * u.eV/ const.k_B).to(u.K)
plt.plot((TTHe*const.k_B).to(u.MeV), [ratio_mass_Hb(T) for T in TTHe], linestyle="-", color=p[0].get_color())
plt.plot((TT*const.k_B).to(u.MeV), ratio_mass_Db(TT), linestyle="-", label="D")
plt.plot((TTHe*const.k_B).to(u.MeV), [ratio_mass_Heb(T) for T in TTHe], linestyle="-", label="$^4$He")
plt.axvline((Tfreeze * const.k_B).to(u.MeV).value, color="b", linestyle="--", label="$T_\mathrm{freeze}$")
plt.axvline((Tnuc * const.k_B).to(u.MeV).value, color="r", linestyle="--", label="$T_\mathrm{nuc}$")
plt.xscale("log")
#plt.yscale("log")
plt.xlabel(r"$T$ [MeV]", fontsize=16)
plt.ylabel(r"Element mass fractions $\rho_X / \rho_b$", fontsize=16)
plt.gca().invert_xaxis()
#plt.ylim(1e-2, 1)
plt.grid()
plt.legend(fontsize=14, loc="upper left", ncol=2)
plt.show()
<Figure size 640x480 with 1 Axes>
TT = (np.logspace(np.log10(7e4), np.log10(5e7), 100) * u.eV/ const.k_B).to(u.K)

E = (TT * const.k_B).to(u.MeV)
T_before_me = np.array([T.value for T in TT if (T * const.k_B >= me)])*u.K
T_after_me = np.array([T.value for T in TT if (T * const.k_B <= me)])*u.K

fig = plt.figure()
plt.axvline((Tfreeze * const.k_B).to(u.MeV).value, color="b", linestyle="--", label="$T_\mathrm{freeze}$")
plt.axvline((Tnuc * const.k_B).to(u.MeV).value, color="r", linestyle="--", label="$T_\mathrm{nuc}$")
plt.axvline(me.to(u.MeV).value, linestyle=":", label="$m_e$", color="gray")
plt.axhline(1/taun.value, linestyle="--", label=r"$\tau_n^{-1}$", color="green")
#plt.plot(tt, np.array([-dXndt(Xnt(t*u.s), t) for t in tt])/u.s, label="$\Gamma_n(T)$")
p = plt.plot(E, np.array([lnp(T).value for T in TT])/u.s, lw=2, label=r"$\Gamma_{n\to p}(T)$ exact")
plt.plot(E, lambda_np_Bernstein(TT), label="approximation without\nfree neutron decay", color=p[0].get_color(), linestyle="--")
p = plt.plot(E, np.array([lpn(T).value for T in TT])/u.s, lw=2, label=r"$\Gamma_{p\to n}(T)$ exact")
#plt.plot(E, np.array([lnp(T).value+lpn(T).value for T in TTT])/u.s, label="$\Gamma_n(T)$")
#plt.plot((T_before_me * const.k_B).to(u.MeV), H(T_before_me, gstar=gstar), lw=2, label="$H(T)$", color="k")
#plt.plot((T_after_me * const.k_B).to(u.MeV), H(T_after_me, gstar=gstar_after), lw=2,color="k")
plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.grid()
plt.xlabel("$T$ [MeV]", fontsize=16)
plt.ylabel("Rates [1/s]", fontsize=16)
plt.axvspan((Tnuc * const.k_B).to(u.MeV).value, 5e-2, alpha=0.2, color="gray")  #horizontal shading
plt.xlim((50, 5e-2))
plt.ylim((1e-1/taun.value, 5e4))
plt.legend(fontsize=12, ncol=1, loc="lower left")

secax = plt.gca().twiny()
ttt = T_to_t_Weinberg(TT.value)
secax.plot(ttt, 1/(2*ttt), linestyle="-", color="k", lw=2, label="H(T)")
#secax.plot(ttt, H(t_to_T_Weinberg(ttt.value), gstar=gstar_after), linestyle="-", color="r", lw=2, label="H(T)")
#secax.plot(tt_before_me, H(t_to_T(tt_before_me, gstar=gstar), gstar=gstar), linestyle="-", color="k", lw=2, label="H(T)")
#secax.plot(tt_after_me, H(t_to_T(tt_after_me, gstar=gstar_after), gstar=gstar_after), linestyle="-", color="k", lw=2, label="H(T)")
#secax.plot(tt*u.s, H(t_to_T(tt*u.s, gstar=gstar), gstar=gstar), linestyle="none")
#secax.plot(tt*u.s, H(t_to_T(tt*u.s, gstar=gstar_after), gstar=gstar_after), linestyle="none")
secax.set_xscale("log")
secax.set_xlabel('$t$ [s]', fontsize=16)
secax.legend(fontsize=12, loc="upper left")

fig.tight_layout()
plt.show()
<Figure size 640x480 with 2 Axes>
(13.6*u.eV / const.k_B).to(u.K)
Loading...

Taux de désintégration du neutron (plateau de la courbe)`

1/(lnp(TTT[-1]).value/u.s)
Loading...
TT = (np.logspace(np.log10((Tnuc * const.k_B).to(u.eV).value), np.log10(1e5), 100) * u.eV/ const.k_B).to(u.K)

fig = plt.figure()
plt.plot((TT*const.k_B).to(u.MeV), ratio_np(TT), lw=2, label=rf"$n_n / n_p$")
plt.plot((TT*const.k_B).to(u.MeV), ratio_Dn(TT), lw=2, label=rf"$n_D / n_n$")
#plt.yscale("log")
plt.xscale("log")
plt.gca().invert_xaxis()
plt.xlabel("$T$ [MeV]", fontsize=16)
plt.ylabel("Density ratios", fontsize=16)
#plt.axvline((Tfreeze * const.k_B).to(u.MeV).value, color="b", linestyle="--", label="$T_\mathrm{freeze}$")
plt.axvline((Tnuc * const.k_B).to(u.MeV).value, color="r", linestyle="--", label="$T_\mathrm{nuc}$")
plt.legend(fontsize=14, ncol=1)
plt.axvspan((Tnuc * const.k_B).to(u.MeV).value, 6.5e-2, alpha=0.2, color="gray")  #horizontal shading
plt.xlim((np.max(TT)*const.k_B).to(u.MeV).value, 6.5e-2)
#plt.ylim(0.0, 0.52)
plt.grid()

secax = plt.gca().twiny()
ttt = T_to_t_Weinberg(TT.value)
secax.plot(ttt, Xnt(ttt), linestyle="none")
secax.set_xscale("log")
secax.set_xlabel('$t$ [s]', fontsize=16)
fig.tight_layout()
plt.show()
<Figure size 640x480 with 2 Axes>
def Ob0(eta):
    return eta / (2.73e-8 * (H0.value / 100)**2)

etas = np.linspace(1e-11, 1e-8, 50)
Ob0s = Ob0(etas)

fig = plt.figure()
plt.plot(etas, [Yp_eta(eta) for eta in etas])
plt.xscale("log")
plt.grid()
plt.xlabel("$\eta$", fontsize=16)
plt.ylabel("$Y_p$", fontsize=16)

secax = plt.gca().twiny()
secax.plot(Ob0s, [Yp_eta(eta) for eta in etas], linestyle="none")
secax.set_xscale("log")
secax.axvline(Planck18.Ob0, color="k", linestyle="--", label="$\Omega_b^0$ from Planck 2018")
secax.set_xlabel('$\Omega_b^0$', fontsize=16)
secax.legend(fontsize=14)

fig.tight_layout()
plt.show()
<Figure size 640x480 with 2 Axes>
Exercices
Thermal history
Exercices
Large-scale structures