forked from MiguelEA/nudec_BSM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
BasicModules.m
157 lines (125 loc) · 10 KB
/
BasicModules.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
(* ::Package:: *)
(************************************************************************)
(* This file was generated automatically by the Mathematica front end. *)
(* It contains Initialization cells from a Notebook file, which *)
(* typically will have the same name as this file except ending in *)
(* ".nb" instead of ".m". *)
(* *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs. Doing so is equivalent *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end. *)
(* *)
(* DO NOT EDIT THIS FILE. This entire file is regenerated *)
(* automatically each time the parent Notebook file is saved in the *)
(* Mathematica front end. Any changes you make to this file will be *)
(* overwritten. *)
(************************************************************************)
(* ::Input::Initialization:: *)
(* Load the tabulated leading order finite temperature QED corrections to the electromagnetic energy density and pressure *)
Clear[Pint,dPintdT,d2PintdT2];
SetDirectory[NotebookDirectory[]];
data=Import["QED/QED_P_int.cvs","Table",HeaderLines-> 3];
auxxvec1=Take[data,{1,Length[data]},{1,2}];
auxxvec2=Take[data,{1,Length[data]},{2,3}];
auxxvec2[[All,1]]=0;
Pint:=Interpolation[auxxvec1+auxxvec2];(* For LO QED use only auxxvec1 *)
data=Import["QED/QED_dP_intdT.cvs","Table",HeaderLines-> 3];
auxxvec3=Take[data,{1,Length[data]},{1,2}];
auxxvec4=Take[data,{1,Length[data]},{2,3}];
auxxvec4[[All,1]]=0;
dPintdT:=Interpolation[auxxvec3+auxxvec4];(* For LO QED use only auxxvec3 *)
data=Import["QED/QED_d2P_intdT2.cvs","Table",HeaderLines-> 3];
auxxvec5=Take[data,{1,Length[data]},{1,2}];
auxxvec6=Take[data,{1,Length[data]},{2,3}];
auxxvec6[[All,1]]=0;
d2PintdT2:=Interpolation[auxxvec5+auxxvec6];(* For LO QED use only auxxvec5 *)
Clear[data];
(* Uncomment in order to remove QED thermal corrections *)
(*
Clear[Pint,dPintdT,d2PintdT2];
Pint[x_]:=0;dPintdT[x_]:=0;d2PintdT2[x_]:=0;*)
(* ::Input::Initialization:: *)
(* Load the corrections to the rates as a result of a non-negligible electron mass *)
SetDirectory[NotebookDirectory[]];
Clear[geL,geR,g\[Mu]L,g\[Mu]R,fann\[Nu]e,fsca\[Nu]e,fann\[Nu]\[Mu],fsca\[Nu]\[Mu],fscan\[Nu]e,fscan\[Nu]\[Mu]];
data=Import["SM_Rates/nue_ann.dat","Table",HeaderLines-> 1];
fann\[Nu]e=Interpolation[data];
data=Import["SM_Rates/nue_scatt.dat","Table",HeaderLines-> 1];
fsca\[Nu]e=Interpolation[data];
data=Import["SM_Rates/numu_ann.dat","Table",HeaderLines-> 1];
fann\[Nu]\[Mu]=Interpolation[data];
data=Import["SM_Rates/numu_scatt.dat","Table",HeaderLines-> 1];
fsca\[Nu]\[Mu]=Interpolation[data];
data=Import["SM_Rates/sigmav_nue.dat","Table",HeaderLines-> 1];
fscan\[Nu]e=Interpolation[data];
data=Import["SM_Rates/sigmav_numu.dat","Table",HeaderLines-> 1];
fscan\[Nu]\[Mu]=Interpolation[data];
(* Uncomment if one wants to consider me = 0 in the rates *)
(*Clear[fann\[Nu]e,fsca\[Nu]e,fann\[Nu]\[Mu],fsca\[Nu]\[Mu],fscan\[Nu]e,fscan\[Nu]\[Mu]];
fann\[Nu]e[x_]:=1;fsca\[Nu]e[x_]:=1;fann\[Nu]\[Mu][x_]:=1;fsca\[Nu]\[Mu][x_]:=1;fscan\[Nu]e[x_]:=1;fscan\[Nu]\[Mu][x_]:=1;*)
(* Spin-Statistics Corrections to the MB interaction rates. Set to 1 for the MB case *)
faFD = 0.884;fsFD = 0.829;fnFD = 0.852;
(* Low energy couplings desciribin NC and CC interactions of neutrinos and electrons *)
(* See Table 10.3 of the EW review of the PDG and Table 6 of 1303.5522 *)
geL=0.727;geR=0.233;g\[Mu]L=-0.273;g\[Mu]R=0.233;
Func[T1_,T2_,\[Mu]1_,\[Mu]2_]:=32 *faFD*(Exp[2 \[Mu]1/T1]T1^9-Exp[2 \[Mu]2/T2]T2^9)+ 56 *fsFD*T1^4 T2^4 Exp[\[Mu]1/T1]Exp[\[Mu]2/T2]*(T1-T2);
(* Energy Density Transfer Rates: \[Delta]\[Rho]/\[Delta]t *)
\[CapitalDelta]\[Rho]SM\[Nu]e[T\[Gamma]_,T\[Nu]e_,T\[Nu]\[Mu]_,\[Mu]\[Nu]e_,\[Mu]\[Nu]\[Mu]_]:=FacMeVtos1 (GF^2)/\[Pi]^5 (4 *(geL^2+geR^2)(32*faFD*fann\[Nu]e[T\[Gamma]]*(T\[Gamma]^9-Exp[2 \[Mu]\[Nu]e/T\[Nu]e]T\[Nu]e^9)+56*fsFD* T\[Gamma]^4 *T\[Nu]e^4*fsca\[Nu]e[T\[Gamma]]*Exp[\[Mu]\[Nu]e/T\[Nu]e]*(T\[Gamma]-T\[Nu]e))+2Func[T\[Nu]\[Mu],T\[Nu]e,\[Mu]\[Nu]\[Mu],\[Mu]\[Nu]e])/.GF-> 1.1663787 10^-11;
\[CapitalDelta]\[Rho]SM\[Nu]\[Mu][T\[Gamma]_,T\[Nu]e_,T\[Nu]\[Mu]_,\[Mu]\[Nu]e_,\[Mu]\[Nu]\[Mu]_]:=FacMeVtos1 (GF^2)/\[Pi]^5 (4 *(g\[Mu]L^2+g\[Mu]R^2)(32*faFD*fann\[Nu]\[Mu][T\[Gamma]]*(T\[Gamma]^9-Exp[2 \[Mu]\[Nu]\[Mu]/T\[Nu]\[Mu]]T\[Nu]\[Mu]^9)+56 *fsFD*T\[Gamma]^4 *T\[Nu]\[Mu]^4*fsca\[Nu]\[Mu][T\[Gamma]]*Exp[\[Mu]\[Nu]\[Mu]/T\[Nu]\[Mu]]*(T\[Gamma]-T\[Nu]\[Mu]))-Func[T\[Nu]\[Mu],T\[Nu]e,\[Mu]\[Nu]\[Mu],\[Mu]\[Nu]e])/.GF-> 1.1663787 10^-11;
\[CapitalDelta]\[Rho]SM\[Nu][T\[Gamma]_,T\[Nu]_,\[Mu]\[Nu]_]:=1/3 (\[CapitalDelta]\[Rho]SM\[Nu]e[T\[Gamma],T\[Nu],T\[Nu],\[Mu]\[Nu],\[Mu]\[Nu]]+2\[CapitalDelta]\[Rho]SM\[Nu]\[Mu][T\[Gamma],T\[Nu],T\[Nu],\[Mu]\[Nu],\[Mu]\[Nu]]);
(* Number Density Transfer Rates: \[Delta]n/\[Delta]t *)
\[CapitalDelta]nSM\[Nu]e[T\[Gamma]_,T\[Nu]e_,T\[Nu]\[Mu]_,\[Mu]\[Nu]e_,\[Mu]\[Nu]\[Mu]_]:=FacMeVtos1*8*fnFD * (GF^2)/\[Pi]^5 (4 *(geL^2+geR^2)*fscan\[Nu]e[T\[Gamma]]*(T\[Gamma]^8-Exp[2 \[Mu]\[Nu]e/T\[Nu]e]T\[Nu]e^8)+2(T\[Nu]\[Mu]^8 Exp[2 \[Mu]\[Nu]\[Mu]/T\[Nu]\[Mu]]-T\[Nu]e^8 Exp[2 \[Mu]\[Nu]e/T\[Nu]e]) )/.GF-> 1.1663787 10^-11;
\[CapitalDelta]nSM\[Nu]\[Mu][T\[Gamma]_,T\[Nu]e_,T\[Nu]\[Mu]_,\[Mu]\[Nu]e_,\[Mu]\[Nu]\[Mu]_]:=FacMeVtos1*8*fnFD * (GF^2)/\[Pi]^5 (4*(g\[Mu]L^2+g\[Mu]R^2) *fscan\[Nu]\[Mu][T\[Gamma]]*(T\[Gamma]^8-Exp[2 \[Mu]\[Nu]\[Mu]/T\[Nu]\[Mu]]T\[Nu]\[Mu]^8)-(T\[Nu]\[Mu]^8 Exp[2 \[Mu]\[Nu]\[Mu]/T\[Nu]\[Mu]]-T\[Nu]e^8 Exp[2 \[Mu]\[Nu]e/T\[Nu]e]))/.GF-> 1.1663787 10^-11;
\[CapitalDelta]nSM\[Nu][T\[Gamma]_,T\[Nu]_,\[Mu]\[Nu]_]:=1/3 (\[CapitalDelta]nSM\[Nu]e[T\[Gamma],T\[Nu],T\[Nu],\[Mu]\[Nu],\[Mu]\[Nu]]+2\[CapitalDelta]nSM\[Nu]\[Mu][T\[Gamma],T\[Nu],T\[Nu],\[Mu]\[Nu],\[Mu]\[Nu]]);
(* ::Input::Initialization:: *)
(* All the formulae is written for particles with g = 1. *)
(* Numerical integrals are performed with respect to Ered = En * T. *)
(* Massless Species : n, \[Rho], p *)
nBE[T_,\[Mu]_]:=(T^3 PolyLog[3,E^(\[Mu]/T)])/\[Pi]^2;
nFD[T_,\[Mu]_]:=-((T^3 PolyLog[3,-E^((\[Mu]/T))])/\[Pi]^2);
\[Rho]BE[T_,\[Mu]_]:=(3 T^4 PolyLog[4,E^(\[Mu]/T)])/\[Pi]^2;
\[Rho]FD[T_,\[Mu]_]:=-((3 T^4 PolyLog[4,-E^((\[Mu]/T))])/\[Pi]^2);
pBE[T_,\[Mu]_]:=1/3 \[Rho]BE[T,\[Mu]];
pFD[T_,\[Mu]_]:=1/3 \[Rho]FD[T,\[Mu]];
(* Massless Species : \!\(
\*SubscriptBox[\(\[PartialD]\), \(T\)]\ \
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]\) *)
dnBEdT[T_,\[Mu]_]:=(T (-\[Mu] PolyLog[2,E^(\[Mu]/T)]+3 T PolyLog[3,E^(\[Mu]/T)]))/\[Pi]^2;
dnBEd\[Mu][T_,\[Mu]_]:=(T^2 PolyLog[2,E^(\[Mu]/T)])/\[Pi]^2;
d\[Rho]BEdT[T_,\[Mu]_]:=(3 T^2 (-\[Mu] PolyLog[3,E^(\[Mu]/T)]+4 T PolyLog[4,E^(\[Mu]/T)]))/\[Pi]^2;
d\[Rho]BEd\[Mu][T_,\[Mu]_]:=(3 T^3 PolyLog[3,E^(\[Mu]/T)])/\[Pi]^2;
dnFDdT[T_,\[Mu]_]:=(T (\[Mu] PolyLog[2,-E^((\[Mu]/T))]-3 T PolyLog[3,-E^((\[Mu]/T))]))/\[Pi]^2;
dnFDd\[Mu][T_,\[Mu]_]:=-((T^2 PolyLog[2,-E^((\[Mu]/T))])/\[Pi]^2);
d\[Rho]FDdT[T_,\[Mu]_]:=(3 T^2 (\[Mu] PolyLog[3,-E^((\[Mu]/T))]-4 T PolyLog[4,-E^((\[Mu]/T))]))/\[Pi]^2;
d\[Rho]FDd\[Mu][T_,\[Mu]_]:=-((3 T^3 PolyLog[3,-E^((\[Mu]/T))])/\[Pi]^2);
(* Massive Species : n, \[Rho], p *)
nBEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] 1/(Exp[Enred-\[Mu]/T]-1),{Enred,m/T,\[Infinity]}];
nFDM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] 1/(Exp[Enred-\[Mu]/T]+1),{Enred,m/T,\[Infinity]}];
nMBM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] 1/Exp[Enred-\[Mu]/T],{Enred,m/T,\[Infinity]}];
\[Rho]BEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^4 NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] 1/(Exp[Enred-\[Mu]/T]-1),{Enred,m/T,\[Infinity]}];
\[Rho]FDM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^4 NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] 1/(Exp[Enred-\[Mu]/T]+1),{Enred,m/T,\[Infinity]}];
pBEM[T_,\[Mu]_,m_]:=1/(6 \[Pi]^2) T^4 NIntegrate[(Enred^2-m^2/T^2)^(3/2) 1/(Exp[Enred-\[Mu]/T]-1),{Enred,m/T,\[Infinity]}];
pFDM[T_,\[Mu]_,m_]:=1/(6 \[Pi]^2) T^4 NIntegrate[(Enred^2-m^2/T^2)^(3/2) 1/(Exp[Enred-\[Mu]/T]+1),{Enred,m/T,\[Infinity]}];
(* Massive Species : \!\(
\*SubscriptBox[\(\[PartialD]\), \(T\)]\ \
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]\) *)
dndTBEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] ((Enred T-\[Mu]) Csch[1/2 (Enred-\[Mu]/T)]^2)/(4 T^2),{Enred,m/T,\[Infinity]}];
dnd\[Mu]BEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] Csch[1/2 (Enred-\[Mu]/T)]^2/(4 T),{Enred,m/T,\[Infinity]}];
d\[Rho]dTBEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^4 NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] ((Enred T-\[Mu]) Csch[1/2 (Enred-\[Mu]/T)]^2)/(4 T^2),{Enred,m/T,\[Infinity]}];
d\[Rho]d\[Mu]BEM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^4 NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] Csch[1/2 (Enred-\[Mu]/T)]^2/(4 T),{Enred,m/T,\[Infinity]}];
dndTFDM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^2 Quiet[NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] ((Enred -\[Mu]/T) Sech[(Enred -\[Mu]/T)/2]^2)/4,{Enred,m/T,\[Infinity]}]];
d\[Rho]dTFDM[T_,\[Mu]_,m_]:=
1/(2 \[Pi]^2) T^3 Quiet[NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] ((Enred -\[Mu]/T) Sech[(Enred -\[Mu]/T)/2]^2)/4,{Enred,m/T,\[Infinity]}]];
dnd\[Mu]FDM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^2 Quiet[NIntegrate[Enred Sqrt[Enred^2-m^2/T^2] 1/(2(Cosh[Enred -\[Mu]/T]+1)),{Enred,m/T,\[Infinity]}]];
d\[Rho]d\[Mu]FDM[T_,\[Mu]_,m_]:=1/(2 \[Pi]^2) T^3 Quiet[NIntegrate[Enred^2 Sqrt[Enred^2-m^2/T^2] 1/(2(Cosh[Enred -\[Mu]/T]+1)),{Enred,m/T,\[Infinity]}]];
(* Entropy Densities *)
sFD[T_,\[Mu]_]:= (\[Rho]FD[T,\[Mu]]+ pFD[T,\[Mu]]-\[Mu] nFD[T,\[Mu]])/T;
sFDM[T_,\[Mu]_,m_]:= (\[Rho]FDM[T,\[Mu],m]+ pFDM[T,\[Mu],m]-\[Mu] nFDM[T,\[Mu],m])/T;
sBE[T_,\[Mu]_]:= (\[Rho]BE[T,\[Mu]]+ pBE[T,\[Mu]]-\[Mu] nBE[T,\[Mu]])/T;
sBEM[T_,\[Mu]_,m_]:= (\[Rho]BEM[T,\[Mu],m]+ pBEM[T,\[Mu],m]-\[Mu] nBEM[T,\[Mu],m])/T;
(* ::Input::Initialization:: *)
m0=0.511; (* Used to generate a dimensionless scale factor *)
me = 0.511; (* Electron Mass in MeV *)
FacMeVtos1=1/(6.58212 10^-16 10^-6); (* This is used to convert MeV to s^-1 *)
Mpl=1.2209 10^19 10^3 (* Planck Mass in MeV *);