This is the project page for BaBooNs.
Kaitlin Cook: Parameters
Justin Michael Brown: Code
Rob Almus: Analysis
Aim: Investigate the effects of non-thermal baryon injection on BBN.
Introduction
If dark matter particles (e.g. WIMPS) exist, they inevitably did so during the big bang. They also may decay during BBN into non-thermal particles - this may be in the form of photons, leptons (electron/positron pairs and/or neutrinos) or into hadronic channels (n/p, mesons, nuclei). These decays may have two influences into BBN:
1) Change the expansion rate due to the injection of relativistic species.
2) Effect the abundance of light nuclei produced in BBN via reactions with these decay particles.
In case (1), this effect is important only if the decay channels are relativistic, thus adding a lot of entropy into the universe. In this project, we will consider the injection of neutrons, so we will not consider this effect.
It is thus the goal of this project to investigate the effect of the injection of non-thermal neutrons on the light element abundances.
Key parameters:
The mean life of the decaying particle X → 2n (some magic decay. Not very physical, probably. But hey, we don't actually know what WIMPS are).
The energy of these neutrons (although if it is high enough, it may not matter too much?)
The abundance of X. Jedamzik, 2004: 7Be is suppressed if O(10^-5) neutrons per proton are injected. Increases D/H. Lifetime of ~1000 s required.
n interaction rates with every other nucleus.
Mass of the WIMP. Chosen to be 100 GeV, based on the wikipedia WIMP entry, the font of all wisdom.
This figure:
from CDMS in 2004 suggests that this value isn't too far off.
Code:
How to add an Isotope to bigbang
Add pointer identifier, id, in network.dek (e.g. ini56)
In subroutine init_bigbang:
Add 1 to ionmax
Set id number (usually +1 of the last isotope)
id = number
Set ratnam (name e.g. “ni56”)
ionam (id) = isotope name
Set zion (# of protons)
zion (id) = Z
Set aion (mass number)
aion (id) = A
Set bion (binding energy)
bion (id) = binding energy
In net_input
Add initial abundance in xin
xin (id)
In init_isotope_rate_pointers
Set pointer to 0
id = 0
How to add a Reaction Rate to bigbang
Add reaction subroutine
Add reaction and reverse reaction pointers, irid and irrid, in network.dek
In subroutine init_bigbang:
Add 2 to nrat
Set id numbers, usually +1 of last reaction
irid = number
irrid = number + 1
Set names
ratnam (irid) = reaction name
ratnam (irrid) = reverse reaction name
In subroutine init_isotope_rate_pointers
Set pointers to 0
irid = 0
irrid = 0
In subroutine bigbangrat
Call reaction subroutine
call rate_id (btemp,bden,ratraw(irid),dratrawdt(irid),dratrawdd(irid),ratraw(irrid),dratrawdt(irrid),dratrawdd(irrid)
Arguments
temp ⇒ temperature in K
den ⇒ density in g/cm^3
fr ⇒ forward reaction rate, N_a<sigma v>
dfrdt ⇒ derivative of fr with respect to temp
dfrdd ⇒ derivative of fr with respect to den
rr ⇒ reverse reaction rate, N_a<sigma v>
drrdt ⇒ derivative of fr with respect to temp
drrdd ⇒ derivative of fr with respect to den
In subroutine bigbangtab
Add density dependence
dtab(irid) = 1.0d0 for decay
dtab(irid) = bden for 2-body interaction
Likewise for irrid
In subroutine rhs:
Add dy/dt terms
dydt(isotope1) = dydt(isotope1) - y(isotope1) * y(isotope2) * rate(irid)
dydt(isotope2) = dydt(isotope2) - y(isotope1) * y(isotope2) * rate(irid)
dydt(output) = dydt(output) + y(isotope1) * y(isotope2) * rate(irid)
Likewise for reverse reaction
In subroutine bbigbang
Add tree calls for reaction if they don't exist
call tree(isotope1,isotope1,eloc,neloc,nterms,nzo,iloc,jloc,np)
call tree(isotope1,isotope2,eloc,neloc,nterms,nzo,iloc,jloc,np)
call tree(isotope2,isotope1,eloc,neloc,nterms,nzo,iloc,jloc,np)
call tree(isotope2,isotope2,eloc,neloc,nterms,nzo,iloc,jloc,np)
call tree(output,isotope1,eloc,neloc,nterms,nzo,iloc,jloc,np)
call tree(output,isotope2,eloc,neloc,nterms,nzo,iloc,jloc,np)
In subroutine sbigbang
Add Jacobian elements
For example, for d(isotope1)/d(isotope2)
a1 = -y(isotope1)*ratdum(irid)
nt = nt + 1, only if this Jacobian element doesn't already exist
iat = eloc(nt), only if this Jacobian element doesn't already exist
dfdy(iat) = dydy(iat) + a1
xsum(isotope2) = xsum(isotope2) + a1 * mion(isotope1)
Likewise for d(isotope1)/d(isotope1), d(isotope2)/d(isotope1), etc.
For example, for d(output)/d(isotope2)
a1 = -y(isotope1)*ratdum(irid)
nt = nt + 1, only if this Jacobian element doesn't already exist
iat = eloc(nt), only if this Jacobian element doesn't already exist
dfdy(iat) = dydy(iat) + a1
xsum(isotope2) = xsum(isotope2) + a1 * mion(output)
Likewise for d(output)/d(isotope1)
Do likewise for the reverse reaction rate, making sure not to repeat any of the above lines that shouldn't be done if the element already exists
Cross Sections
To add in reaction rates, as above, obviously one needs to find some reaction rates. Unfortunately, these don't exist at the temperatures and energy distributions of interest, as far as the authors are aware. Instead, we used the Ramsauer model, with parameters from Gowder et. al. 2005 to find total scattering cross sections. We then took a nominal value of 10% for the destruction cross sections. See the presentation for details. The cross sections were then put into a reaction rate, using a dirac function distribution of neutron energies.
Presentation