Constants and “constants”

Constants

Mpl=1.22*10^19(*GeV*);GeVInvTos=1/(1.52*10^24)(*GeV*s*);etaB=6.2*10^-10;mpMeV=938;mpGeV=0.938;mpiGeV=0.135;
GeVInvTocm=2*10^-14(*cm*GeV*);

BD=2.23(*MeV*);BHe=28.3(*MeV*);BHe3=7.72;BH3=8.48;BLi6=31.99;BLi7=39.24;BC12=92.2;
GF=1.17*10^(-5)(*GeV^-2*);
ccms=29979245800(*cm/s*);

Thermodynamic quantities in the early Universe (SM EoS)

(*it's important to upload the file with number of relativistic degrees of freedom in the SM! Change the path below!*)

dataEoS=Import["/Users/kolesova/Work/Bern_projects/data/standardmodel.dat","Table"]

dataEoS[[7]]

datagstar=Table[{dataEoS[[j,1]],dataEoS[[j,3]]*30/Pi^2},{j,8,100}]

ListLogLogPlot[datagstar,Joined->True,AxesLabel->{"T[MeV]","(SubStar[g])"},GridLines->{{0.08,1.4},{}}]

gstarInt=Interpolation[datagstar]

datagstarfull=Table[{dataEoS[[j,1]],dataEoS[[j,3]]*30/Pi^2},{j,8,Length[dataEoS]}];

ListLogLogPlot[datagstarfull,PlotRange->{{0.01,10^4},{3,110}},Joined->True]

datarho3p=Table[{dataEoS[[j,1]],(dataEoS[[j,3]]-3*dataEoS[[j,2]])/dataEoS[[j,3]]},{j,8,Length[dataEoS]}];

ListLogLogPlot[datarho3p,PlotRange->{{0.01,10^4},{10^-4,0.4}},Joined->True]

Neutron freeze out

(*effective number of relativistic degrees of freedom around 1MeV*)
gstar=10.75;
(*Hubble rate at x=1*)
Hx1=1.13(*s^-1*);
(*effective number of relativistic degrees of freedom around 80keV*)
Neff=3.;gstarNS=2+7/8*2*Neff*(4/11)^(4/3);
Hx1NS=Hx1*Sqrt[gstarNS/gstar];
(*neutron lifetime*)
taun=886.7(*s*);
(*interaction rate for neutron-proton conversion*)
lamnp[x_]:=(255/(taun*x^5))(12+6*x+x^2);XnEQ[x_]:=1/(1+Exp[x]);
(*difference of mass between proton and neutron, assigned by Q in the lectures*)
deltam=1.293(*MeV*);

(*temperature in terms of x*)Tem[x_]:=deltam/x;

(*it's reasonable to solve the equation until x=18*)

Tem[18]

gstarNS

xmax=18;

(*solution of the Boltzmann equation with gstar=10.75= const., without taking into account neutron decay*)

solXn=NDSolve[Xn'[x]==x/Hx1*lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])&&Xn[0.1]==XnEQ[0.1],Xn,{x,0.1,xmax}]

(*solution of the Boltzmann equation with gstarNS=3.36=const., without taking into account neutron decay*)

solXnNS=NDSolve[Xn'[x]==x/Hx1NS*lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])&&Xn[0.1]==XnEQ[0.1],Xn,{x,0.1,xmax}]

(*solution of the Boltzmann equation with realistic evolution of gstar, without taking into account neutron decay*)

solXngstNum=NDSolve[Xn'[x]==x/Hx1*Sqrt[gstar/gstarInt[Tem[x]]]*lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])&&Xn[0.1]==XnEQ[0.1],Xn,{x,0.1,xmax}]

(*how does the realistic gstar differ from 10.75*)

LogLogPlot[Sqrt[gstar/gstarInt[Tem[x]]],{x,0.1,16}]

(*taking into account neutron decay, but keeping gstar=10.75= const.*)

solXnDec=NDSolve[Xn'[x]==x/Hx1*(lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])-1/taun*Xn[x])&&Xn[0.1]==XnEQ[0.1],Xn,{x,0.1,xmax}]

(*most realistic equation: taking into account neutron decay, realistic gstar*)

solXnFull=NDSolve[Xn'[x]==x/Hx1*Sqrt[gstar/gstarInt[Tem[x]]]*(lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])-1/taun*Xn[x])&&Xn[0.1]==XnEQ[0.1],Xn,{x,0.1,xmax}]

LogLogPlot[{Xn[x]/.solXn[[1]],Xn[x]/.solXnNS[[1]],Xn[x]/.solXngstNum[[1]],Xn[x]/.solXnDec[[1]],Xn[x]/.solXnFull[[1]],XnEQ[x],1/(1+1/0.18),1/(1+1/0.14)},{x,0.1,xmax},PlotRange->{0.1,0.5},PlotStyle->{Thickness[.003],Thickness[.003], Thickness[.003],Thickness[.003],Thickness[.006],Dotted,Dashed,Dashed},GridLines->{{deltam/0.08},{}},Frame->True,FrameLabel->{"x=Q/T","(Subscript[X, n])=(Subscript[n, n])/(Subscript[n, B])"},PlotLegends->{"(Subscript[X, n]): (SubStar[g])=10.75=const., no neutron decay","(Subscript[X, n]): (SubStar[g])=3.36=const., no neutron decay","(Subscript[X, n]): (SubStar[g]) realistic, no neutron decay","(Subscript[X, n]): (SubStar[g])=10.75, including neutron decay","(Subscript[X, n]): (SubStar[g]) realistic, including neutron decay","(Subscript[X, n, EQ])","(Subscript[n, n])/(Subscript[n, p])=0.18","(Subscript[n, n])/(Subscript[n, p])=0.14"},LabelStyle->20,ImageSize->600]

(*what if the weak interaction is 10x stronger?*)

solXnMod=NDSolve[Xn'[x]==x/Hx1*10*lamnp[x]*((1-Xn[x])*Exp[-x]-Xn[x])&&Xn[0.05]==XnEQ[0.05],Xn,{x,0.05,xmax}]

(*if the weak interaction is stronger, the neutrons stay longer in equilibirum and their final abundance is lower*)

LogLogPlot[{Xn[x]/.solXn[[1]],Xn[x]/.solXnMod[[1]],XnEQ[x]},{x,0.1,xmax},PlotRange->{0.0001,0.5},PlotStyle->{Thickness[.003],Thickness[.003],Dotted},Frame->True,FrameLabel->{"x=Δm/T","(Subscript[X, n])=(Subscript[n, n])/(Subscript[n, B])"},LabelStyle->20,ImageSize->600,PlotLegends->{"(Subscript[X, n]) as before, no neutron decay","(Subscript[X, n]) with 10x larger cross section","(Subscript[X, n, EQ])"},GridLines->{{xFOnum1,xFOnum2},{}}]

(*when H~Γ? It acutally tells us when the abundance stops changing a lot..*)

xFOnum1=FindRoot[Hx1/x^2==lamnp[x],{x,2}][[1,2]]

xFOnum2=FindRoot[Hx1/x^2==10*lamnp[x],{x,3}][[1,2]]