Quit[];

Formulas from the notes

1.  Setup

params={
T0->0.234meV*a0
,Tovarec->T3*0.3eV
,H0->a0*100*h100*km/sec/Mpc
,h100->h7*0.7
,gastS0->3.91
,gast0->3.36
,gastSeq->gastS0
,gasteq->gast0
,gastSrec->gastS0
,gastrec->gast0
,gastS->g100*106.75
,gast->g100*106.75
(*,arec->a0(1+z0)/(1+zrec)*)
,arec->Trec/Tovarec
,Trec->T0
(*,z0->0*)
(*,zrec->1100*)
,Mpl->2.435*(10^18)GeV
,meV->(10^(-3))eV
,GeV->(10^9)eV
,TeV->(10^12)eV
,Hz->1/sec
,km->(10^5)cm
,GaussHL->0.019535277122580234` (eV^2)
,Mpc->3.0856775814913676`*^24cm
,pc->(10^(-6))Mpc
,Hcal->π*(gastS0^(2/3))/Sqrt[90]*(T0^2)/Mpl*((gast^(1/2))/(gastS^(2/3)))
,Hcaleq->π*(gastS0^(2/3))/Sqrt[90]*(T0^2)/Mpl*((gasteq^(1/2))/(gastSeq^(2/3)))
,ℏ->6.58211956856515`*^-16 eV sec
,ℏc->0.00001973269804310045` cm eV
,cKL->(10^(-2))
,cBS->1.3*(10^(-2))
,ye->Sqrt[2]me/v
,me->511.*(10^(-6))GeV
,v->246.GeV
,αY->α2*0.01
,σY->σ2*100T[η]
,μY5->μ3*(10^(-3))T[η]
};

T[η_]=(((gastS0/gastS))^(1/3))*T0;
s[η_]=((2(π^2))/45)gastS*((T[η])^3);
ρ[η_]=((π^2)/30)gast*((T[η])^4);
a[η_]=Hcal(η-η1)+a1;
H[η_]=Hcal*((a[η])^(-1));
t[η_]=t1+a1(η-η1)+(1/2)Hcal(((η-η1))^2);
(*Γf[η_]=cKL*((ye^2)/(8π))*T[η]; (* Kamada & Long (2016) *) *)
Γf[η_]=cBS*(ye^2)*T[η]; (* Bodeker & Schroder (2019) *)
(*H[η]==a'[η]/a[η];*)
(*a[η]==t'[η];*)
rule1={a1->Hcal*η1,t1->Hcal*((η1)^2)/2};
rulet=Join[rule1,{η->Sqrt[2tt/Hcal]}];
rulea=Join[rule1,{η->aa/Hcal}];
Print["H(η) = "
,H[η]/.rule1//Simplify," = "
,H[η]/.rulet//Simplify//PowerExpand," = "
,H[η]/.rulea//Simplify//PowerExpand
];
Print["H(η)/a(η) = "
,H[η]/a[η]/.rule1//Simplify," = "
,H[η]/a[η]/.rulet//Simplify//PowerExpand," = "
,H[η]/a[η]/.rulea//Simplify//PowerExpand
];
Print["a(η) = "
,a[η]/.rule1//Simplify," = "
,a[η]/.rulet//Simplify//PowerExpand," = "
,a[η]/.rulea//Simplify//PowerExpand
];
Print["t(η) = "
,t[η]/.rule1//Simplify," = "
,t[η]/.rulet//Simplify//PowerExpand," = "
,t[η]/.rulea//Simplify//PowerExpand
];
Print["Hcal = ",(Hcal/(eV*(a0^2)))//.params," eV (Subscript[a, 0]^2)"];

ηew=η/.First[Solve[(100GeV)==T[η]/a[η]/.rule1//Simplify,η]];
Tew=T[ηew]/.rule1//Simplify;
aew=a[ηew]/.rule1//Simplify;
tew=t[ηew]/.rule1//Simplify;
Hew=H[ηew]/.rule1//Simplify;
Print["(Subscript[η, ew]) = ",((ηew*ℏ)/(sec*(a0^(-1))))//.params//PowerExpand," sec (Subscript[a, 0]^(-1))"];
Print["(Subscript[t, ew]) = ",((tew*ℏ)/sec)//.params//PowerExpand," sec"];
Print["(Subscript[a, ew]) = ",(aew/(a0^1))//.params//PowerExpand," (Subscript[a, 0])"];
Print["(Subscript[H, ew]) = ",(Hew/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(H/a)_ew = ",((Hew/aew)/eV)//.params//PowerExpand," eV"];
Print["(Subscript[T, ew]) = ",(Tew/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(T/a)_ew = ",((Tew/aew)/GeV)//.params//PowerExpand," GeV"];

ηf=η/.First[Solve[Γf[η]η==1,η]];
Tf=T[ηf]/.rule1//Simplify;
af=a[ηf]/.rule1//Simplify;
tf=t[ηf]/.rule1//Simplify;
Hf=H[ηf]/.rule1//Simplify;
Print["(Subscript[η, f]) = ",((ηf*ℏ)/(sec*(a0^(-1))))//.params//PowerExpand," sec (Subscript[a, 0]^(-1))"];
Print["(Subscript[t, f]) = ",((tf*ℏ)/sec)//.params//PowerExpand," sec"];
Print["(Subscript[a, f]) = ",(af/(a0^1))//.params//PowerExpand," (Subscript[a, 0])"];
Print["(Subscript[H, f]) = ",(Hf/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(H/a)_f = ",((Hf/af)/eV)//.params//PowerExpand," eV"];
Print["(Subscript[T, f]) = ",(Tf/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(T/a)_f = ",((Tf/af)/GeV)//.params//PowerExpand," GeV"];

2.  Magnetic field generation

kCPI=2αY*μY5/π;
λCPI=2π/kCPI;
ηCPI=2σY/(kCPI^2);
TCPI=T[ηCPI]/.rule1//Simplify;
aCPI=a[ηCPI]/.rule1//Simplify;
tCPI=t[ηCPI]/.rule1//Simplify;
HCPI=H[ηCPI]/.rule1//Simplify;
Print["(Subscript[k, CPI]) = ",((kCPI*((ℏc)^(-1)))/((cm^(-1))*a0))//.params//PowerExpand," (cm^(-1)) (Subscript[a, 0])"];
Print["(k/a)_CPI = ",((kCPI/aCPI*((ℏc)^(-1)))/(cm^(-1)))//.params//PowerExpand," (cm^(-1))"];
Print["(Subscript[λ, CPI]) = ",((λCPI*((ℏc)^1))/(cm*(a0^(-1))))//.params//PowerExpand," cm (Subscript[a, 0]^(-1))"];
Print["(aλ)_CPI = ",((aCPI*λCPI*((ℏc)^1))/cm)//.params//PowerExpand," cm"];
Print["(Subscript[η, CPI]) = ",((ηCPI*ℏ)/(sec*(a0^(-1))))//.params//PowerExpand," sec (Subscript[a, 0]^(-1))"];
Print["(Subscript[t, CPI]) = ",((tCPI*ℏ)/sec)//.params//PowerExpand," sec"];
Print["(Subscript[a, CPI]) = ",(aCPI/(a0^1))//.params//PowerExpand," (Subscript[a, 0])"];
Print["(Subscript[H, CPI]) = ",(HCPI/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(H/a)_CPI = ",((HCPI/aCPI)/eV)//.params//PowerExpand," eV"];
Print["(Subscript[T, CPI]) = ",(TCPI/(eV*(a0^1)))//.params//PowerExpand," eV (Subscript[a, 0])"];
Print["(T/a)_CPI = ",((TCPI/aCPI)/GeV)//.params//PowerExpand," GeV"];

(TCPI/aCPI)==(Tf/af);
μY5min=μY5/.Last[Solve[%,μY5]];
μY5min//.params;
Print["|(Subscript[μ, Y, 5])|/T > ",μY5min/T[η]//PowerExpand," = ",μY5min/T[η]//.params//PowerExpand," "];

hBCPI=(π/(6αY))μY5*((T[ηCPI])^2);
ξBCPI=λCPI;
BCPI=Sqrt[hBCPI/ξBCPI];
ρBCPI=(BCPI^2)/2;
ρcritCPI=3(aCPI^2)(Mpl^2)(HCPI^2);
ΩBCPI=ρBCPI/ρcritCPI;
Print["(Subscript[h, B, CPI]) = ",hBCPI," = ",((hBCPI*ℏc)/((GaussHL^2)*Mpc*(a0^3)))//.params//PowerExpand," (G^2)Mpc (Subscript[a, 0]^3)"];
Print["(Subscript[ξ, B, CPI]) = ",ξBCPI//PowerExpand," = ",((λCPI*((ℏc)^1))/(cm*(a0^(-1))))//.params//PowerExpand," cm (Subscript[a, 0]^(-1))"];
Print["(Subscript[B, CPI]) = ",BCPI//PowerExpand," = ",((BCPI*((ℏc)^0))/(GaussHL*(a0^2)))//.params//PowerExpand," G (Subscript[a, 0]^2)"];
Print["(Subscript[Ω, B, cpi]) = ((Subscript[ρ, B])/(Subscript[ρ, c]))_CPI = ",ρBCPI/ρcritCPI//PowerExpand," = ",ρBCPI/ρcritCPI*((ℏc)^0)//.params//PowerExpand," "];

Clear[ηrec];
(*ηrec=Trec/(Hcal*Tovarec);*)
ηrec=Trec/(Hcaleq*Tovarec);
(*ηrec=Sqrt[(90/(π^2))]((Mpl*arec)/((gastrec^(1/2))*(Trec^2)))/.{Trec->Tovarec*arec};*)
B0=BCPI(((ηrec/ηCPI))^(-1/3));
ξB0=ξBCPI(((ηrec/ηCPI))^(2/3));
hB0=ξB0*(B0^2);
Print["(Subscript[η, rec]) = ",((ηrec*ℏ)/(sec*(a0^(-1))))//.params//PowerExpand," sec (Subscript[a, 0]^(-1))"];
Print["(B/(a^2))_0 = ",(((B0/(a0^2))*(ℏ^0))/GaussHL)//.params//PowerExpand," G"];
Print["(a (Subscript[ξ, B]))_0 = ",(((a0*ξB0)*((ℏc)^1)*(ℏ^0))/pc)//.params//PowerExpand," pc"];
Print["((Subscript[h, B])/(a^3))_0 = ",(((hB0/(a0^3))*((ℏc)^1))/((GaussHL^2)Mpc))//.params//PowerExpand," (G^2)Mpc"];

{((B0/(a0^2))/((10^(-8))GaussHL)),((ℏc*a0*ξB0)/(1Mpc))}//.params//PowerExpand
%[[1]]/%[[2]]

3.  Gravitational wave radiation

fgw0=1/(a0*λCPI);
Print["(Subscript[f, gw, 0]) = ",fgw0," = ",((fgw0*(ℏ^(-1)))/Hz)//.params//PowerExpand," Hz"];

Ωgw0MAX=ΩBCPI;
Print["(Subscript[Ω, gw, 0]) << (Subscript[Ω, B, cpi]) ",ΩBCPI//PowerExpand," = ",ΩBCPI*((ℏc)^0)//.params//PowerExpand," "];

hc=Sqrt[((3(((H0/a0))^2))/(2(π^2)))((Ωgw0MAX)/(fgw0^2))];
Print["(Subscript[h, c]) = ",hc//PowerExpand," = ",hc*ℏ//.params//PowerExpand," "];

4.  Baryogenesis and baryon over-production

5.  Dark radiation constraints

6.  Summary of viable parameter space

paramLO=Join[params,{σ2->1,α2->1,g100->1,T3->1,μ3->1}];
paramHI=Join[params,{σ2->1,α2->1,g100->1,T3->1,μ3->10}];
Print["|(Subscript[μ, Y, 5])|/T:  ",{μY5/TCPI//.paramLO,μY5/TCPI//.paramHI}//PowerExpand//MatrixForm];
Print["(T/a)_CPI:  ",{(((TCPI/aCPI))/GeV)//.paramLO,(((TCPI/aCPI))/GeV)//.paramHI}//PowerExpand//MatrixForm," GeV"];
Print["(B/(a^2))_0:  ",{(((B0/(a0^2))*(ℏ^0))/GaussHL)//.paramLO,(((B0/(a0^2))*(ℏ^0))/GaussHL)//.paramHI}//PowerExpand//MatrixForm," G"];
Print["((Subscript[aξ, B]))_0:  ",{(((a0*ξB0)*((ℏc)^1)*(ℏ^0))/pc)//.paramLO,(((a0*ξB0)*((ℏc)^1)*(ℏ^0))/pc)//.paramHI}//PowerExpand//MatrixForm," pc"];
Print["(Subscript[Ω, gw, 0]^((max))):  ",{Ωgw0MAX//.paramLO,Ωgw0MAX//.paramHI}//PowerExpand//MatrixForm];
Print["(Subscript[f, gw, 0]):  ",{((fgw0*(ℏ^(-1)))/Hz)//.paramLO,((fgw0*(ℏ^(-1)))/Hz)//.paramHI}//PowerExpand//MatrixForm," Hz"];

CPI equation of motion

Brule={
BR'[η]->-((k^2)/(σY))BR[η]+(2/π)αY((μY5)/(σY))k BR[η]
,BRc'[η]->-((k^2)/(σY))BRc[η]+(2/π)αY((μY5)/(σY))k BRc[η]
,BL'[η]->-((k^2)/(σY))BL[η]-(2/π)αY((μY5)/(σY))k BL[η]
,BLc'[η]->-((k^2)/(σY))BLc[η]-(2/π)αY((μY5)/(σY))k BLc[η]
};
ρBk[η_]=((k^2)/(2(π^2)))×(1/2)(BRc[η]BR[η]+BLc[η]BL[η]);
hBk[η_]=((k^2)/(2(π^2)))×(k^(-1))(BRc[η]BR[η]-BLc[η]BL[η]);
Simplify[ρBk'[η]/.Brule]/.{BR[η]->((((2(π^2))/(k^2))×2×ρB[η]+((2(π^2))/(k^2))×k×hB[η])/(2BRc[η])),BL[η]->((((2(π^2))/(k^2))×2×ρB[η]-((2(π^2))/(k^2))×k×hB[η])/(2BLc[η]))}//Simplify//Expand
Simplify[hBk'[η]/.Brule]/.{BR[η]->((((2(π^2))/(k^2))×2×ρB[η]+((2(π^2))/(k^2))×k×hB[η])/(2BRc[η])),BL[η]->((((2(π^2))/(k^2))×2×ρB[η]-((2(π^2))/(k^2))×k×hB[η])/(2BLc[η]))}//Simplify//Expand