# AST3220, Cosmo I, project 2 -- Big Bang Nucleosynthesis # Jakob Borg, UiO, 2025 # Module driving the BBN calculations. # Imports: import numpy as np from scipy.integrate import solve_ivp # Own imports: from reaction_rates import ReactionRates # Class with equations for each reaction from cosmology import BackgroundBBN # Class with background equations class BBN: def __init__(self, NR_interacting_species=2, Omega_b0=0.05, N_eff=3) -> None: """ Initialize the BBN calculations including particles up to a given number, NR_interacting_species. This is a number between 2 and 8, determining which particles to be considered in the model. The particle species are n, p, D, T, He3, He4, Li7, Be7. The class also initializes the reaction rates and background cosmology, with the option to change the baryon density and effective number of neutrinos. Parameters ---------- NR_interacting_species : int, optional The number of particle species to be considered, by default 2 Keyword Arguments: Omega_b0 : float, optional The relative baryon density, by default 0.05 N_eff : float, optional The effective number of neutrinos, by default 3 """ self.NR_species = NR_interacting_species ast_msg = "NR_interacting_species must be between 2 and 8" assert self.NR_species <= 8 and self.NR_species >= 2, ast_msg # Define the particle species order and their mass numbers: self.species_labels = ["n", "p", "D", "T", "He3", "He4", "Li7", "Be7"] self.mass_number = [1, 1, 2, 3, 3, 4, 7, 7] # Initialize the reaction rates and background cosmology: self.RR = ReactionRates() self.cosmo = BackgroundBBN(Omega_b0=Omega_b0, N_eff=N_eff) def get_ODE(self, lnT, Y): """ Compute the right hand sides of the ODEs for Boltzmanns equation for particle species in Y, up to the number of interacting species. Parameters ---------- lnT : float The natural log of temperature. Y : np.ndarray Array of relative number densities at current lnT for each species. With all 8 included species Y takes the form: n, p, D, T, He3, He4, Li7, Be7 = Y Returns ------- np.ndarray The differential for each species following the shape of Y. E.g. dY[0] corresponds to dY_n, dY[1] to dY_p, dY[2] to dY_D, etc. Initialized to zero for all species and add to the change for each species. """ # Set current temperature, Hubble parameter and baryon density: T = ... T9 = T / 1e9 Hubble = self.cosmo.get_Hubble(T) rho_b = self.cosmo.get_rho_b(T) # Initialize the differential for each species following the shape of Y, # i.e. dY[0] corresponds to dY_n. Initialized to zero for all species. dY = np.zeros_like(Y) # Example of how to include each interactions: # ~ (n <-> p) (a 1-3) # Relative number densities of neutrons and protons: Y_n, Y_p = Y[0], Y[1] rate_n, rate_p = self.RR.get_weak_force(T9) # Get the weak force rates # The change to the particles on the left hand side of the reaction equation: change_LHS = Y_p * rate_p - Y_n * rate_n # Update the changes to each side: dY[0] += change_LHS # neutron dY[1] -= change_LHS # proton (opposite sign) # Include deuterium if self.NR_species > 2: Y_D = Y[2] # Relative number density of deuterium # (n+p <-> D+gamma) (b.1) Y_np = Y_n * Y_p rate_np, rate_D = self.RR.get_np_to_D(T9, rho_b) change_LHS = ... # Update changes to each side: dY[0] += change_LHS # neutron dY[1] += change_LHS # proton dY[2] -= change_LHS # deuterium (opposite sign) # Include tritium if self.NR_species > 3: Y_T = Y[3] # Relative number density of tritium # (n+D <-> T+gamma) (b.3) Y_nD = Y_n * Y_D rate_nD, rate_T = self.RR.get_nD_to_T(T9, rho_b) change_LHS = ... dY[0] += ... # neutron dY[2] += ... # deuterium dY[3] -= ... # tritium # TWO EQUAL PARTICLES -- MIND THE FACTORS OF 2 # (D+D <-> p + T) (b.8) Y_DD = Y_D * Y_D Y_pT = Y_p * Y_T rate_DD, rate_pT = self.RR.get_DD_to_pT(T9, rho_b) # Different relative changes to each side: change_LHS = 2 * Y_pT * rate_pT - Y_DD * rate_DD change_RHS = 0.5 * Y_DD * rate_DD - Y_pT * rate_pT dY[2] += change_LHS # deuterium dY[1] += change_RHS # proton dY[3] += change_RHS # tritium if self.NR_species > 4: # Include He3 ... if self.NR_species > 5: # Include He4 ... if self.NR_species > 6: # Include Li7 ... if self.NR_species > 7: # Include Be7 ... # Each reaction equation are missing a factor(-1/Hubble) dY = -dY / Hubble return dY def get_IC(self, T_init): """ Defines the initial condition array used in solve_BBN. The only nonzero values are for neutrons and protons, but the shape of Y_init is determined by the number of particles included. Parameters ---------- T_init : float The initial temperature. Returns ------- np.ndarray The initial condition array for the relative number densities of each species. """ # Initialize the array with zeros for all included species: Y_init = np.zeros(self.NR_species) # Calculate the initial number densities for neutrons and protons: Yn_init = ... Yp_init = ... # Assign the initial number densities to the correct indices in Y_init: Y_init[0] = Yn_init Y_init[1] = Yp_init return Y_init def solve_BBN(self, T_init=100e9, T_end=0.01e9, NR_sol_points=1001): """ Solves the BBN-system for a given range of temperature values Keyword Arguments: T_init {float} -- the initial temperature (default: {100e9}) T_end {float} -- the final temperature (default: {0.01e9}) """ # Get initial conditions for the solver: Y_init = self.get_IC(T_init) # Define the logarithmic temperature range for the solver # as our equations are defined over ln(T): lnT_range = [np.log(T_init), np.log(T_end)] lnT_eval = np.linspace(lnT_range[0], lnT_range[-1], NR_sol_points) # Solve the ODE-system using scipy.solve_ivp: sol = solve_ivp( self.get_ODE, t_span=lnT_range, y0=Y_init, t_eval=lnT_eval, method="Radau", rtol=1e-12, atol=1e-12, ) # Store the results: self.T = np.exp(sol.t) self.Y_i = sol.y if __name__ == "__main__": import matplotlib.pyplot as plt # Example use: # Initiate the system including 3 species (neutrons, protons and deuterium): bbn = BBN(NR_interacting_species=3) bbn.solve_BBN(T_end=0.1e9) # solve the system until end temperature 0.1*10^9 K # Plot the mass fraction for each species: fig, ax = plt.subplots() for i, y in enumerate(bbn.Y_i): ax.plot(bbn.T, bbn.mass_number[i] * y, label=bbn.species_labels[i]) ax.legend()