forked from MiguelEA/nudec_BSM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
WIMP_generic.py
198 lines (166 loc) · 11.3 KB
/
WIMP_generic.py
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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
# Last modified 27/05/2020. Miguel Escudero Abenza, [email protected]
# Check ArXiv:1812.05605 for the relevant equations
# Temperature evolution for a generic WIMP taking into account the energy exchange between electrons and neutrinos
# from the WIMP annihilation.
# The command to run should look like
# python WIMP_generic.py 3.4 2 F e 1e-2 test_generic
# or
# python WIMP_generic.py 10.4 2 F nu 1e-4 test_generic
# where 10.4 corresponds to the mass in MeV, 2 are the NDOF, F corresponds to the statitcs
# nu or e will corresdond to the sector to which the WIMP annihilates the most
# 1e-4 will be the BR of the annihilation of the annihilation that proceeds to the other sector
# for instance, nu 1e-4 means that the WIMP annihilates 1 every 10^4 times to electrons
# This will output the temperature evolution to nameoffile.dat and nameoffile_Neff.dat the Neff value
#####################################################################
import os
import sys
import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import quad, odeint
from scipy.special import kv
# Flag to tell how to run the program
if len(sys.argv) < 5:
print("ERROR"+"\n"+"NEED MORE INFORMATION, you should run for example:"+"\n"+" python WIMP_generic.py 10.4 2 F nu 1e-4 test_generic")
sys.exit()
else:
MDM, gDM, SPIN, MODE, BR, file = float(sys.argv[1]), float(sys.argv[2]), str(sys.argv[3]), str(sys.argv[4]), float(sys.argv[5]),str(sys.argv[6])
if BR > 0.5:
print("ERROR, if you are in e mode BR cannot be > 0.5")
sys.exit()
# WIMP Thermodynamics
if SPIN == 'B' or SPIN == 'BOSE' or SPIN == 'Bose':
def rho_DM(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(2*np.pi**2)*T**4*quad(lambda E: E**2*(E**2-(MDM/T)**2)**0.5/(np.exp(E)-1.) ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
def p_DM(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(6*np.pi**2)*T**4*quad(lambda E: (E**2-(MDM/T)**2)**1.5/(np.exp(E)-1.) ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
def drho_DMdT(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(2*np.pi**2)*T**3*quad(lambda E: 0.25*E**3*(E**2-(MDM/T)**2)**0.5*np.sinh(E/2.0)**-2 ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
if SPIN == 'F' or SPIN == 'Fermi' or SPIN == 'Fermion':
def rho_DM(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(2*np.pi**2)*T**4*quad(lambda E: E**2*(E**2-(MDM/T)**2)**0.5/(np.exp(E)+1.) ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
def p_DM(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(6*np.pi**2)*T**4*quad(lambda E: (E**2-(MDM/T)**2)**1.5/(np.exp(E)+1.) ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
def drho_DMdT(T,MDM):
if T < MDM/30.0: return 0.0
else: return gDM/(2*np.pi**2)*T**3*quad(lambda E: 0.25*E**3*(E**2-(MDM/T)**2)**0.5*np.cosh(E/2.0)**-2 ,MDM/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
#####################################################################
#################### RELEVANT CODE STARTS HERE ####################
#####################################################################
# QED O(e^3)
P_int =interp1d(np.loadtxt("QED/QED_P_int.cvs")[:,0],np.loadtxt("QED/QED_P_int.cvs")[:,1]+np.loadtxt("QED/QED_P_int.cvs")[:,2],bounds_error=False,fill_value=0.0,kind='linear')
dP_intdT =interp1d(np.loadtxt("QED/QED_dP_intdT.cvs")[:,0],np.loadtxt("QED/QED_dP_intdT.cvs")[:,1]+np.loadtxt("QED/QED_dP_intdT.cvs")[:,2],bounds_error=False,fill_value=0.0,kind='linear')
d2P_intdT2 =interp1d(np.loadtxt("QED/QED_d2P_intdT2.cvs")[:,0],np.loadtxt("QED/QED_d2P_intdT2.cvs")[:,1]+np.loadtxt("QED/QED_d2P_intdT2.cvs")[:,2],bounds_error=False,fill_value=0.0,kind='linear')
## Uncomment in order to remove the QED corrections
#P_int, dP_intdT, d2P_intdT2 = lambda x: 0, lambda x: 0, lambda x: 0
# All in MeV Units!
GF = 1.1663787e-5*1e-6 #in MeV^{-2}
me = 0.511
Mpl = 1.22091e19*1e3
# Conversion factor to convert MeV^-1 into seconds
FAC = 1./(6.58212e-22)
# Left and Right nu-e couplings as relevant for E < 10 MeV. From the EW review of the PDG
geL, geR, gmuL, gmuR = 0.727, 0.233, -0.273, 0.233
# Thermodynamics
def rho_nu(T): return 2 * 7./8. * np.pi**2/30. * T**4
def rho_gam(T): return 2 * np.pi**2/30. * T**4
def rho_e(T):
if T < me/30.0: return 0.0
else: return 4./(2*np.pi**2)*T**4*quad(lambda E: E**2*(E**2-(me/T)**2)**0.5/(np.exp(E)+1.) ,me/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
def p_e(T):
if T < me/30.0: return 0.0
else: return 4./(6*np.pi**2)*T**4*quad(lambda E: (E**2-(me/T)**2)**1.5/(np.exp(E)+1.) ,me/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
# Derivatives
def drho_nudT(T): return 4*rho_nu(T)/T
def drho_gamdT(T): return 4*rho_gam(T)/T
def drho_edT(T):
if T < me/30.0: return 0.0
else: return 4./(2*np.pi**2)*T**3*quad(lambda E: 0.25*E**3*(E**2-(me/T)**2)**0.5*np.cosh(E/2.0)**-2 ,me/T,100,epsabs=1e-12,epsrel = 1e-12)[0]
if MODE == 'e':
def rho_tot(T_gam,T_nue,T_numu,MDM):
return rho_gam(T_gam) + rho_e(T_gam) + rho_DM(T_gam,MDM) + rho_nu(T_nue) + 2*rho_nu(T_numu) - P_int(T_gam) + T_gam*dP_intdT(T_gam)
def p_tot(T_gam,T_nue,T_numu,MDM):
return 1./3. * rho_gam(T_gam) + p_e(T_gam) + p_DM(T_gam,MDM) + 1./3. * rho_nu(T_nue) + 1./3. * 2*rho_nu(T_numu) + P_int(T_gam)
if MODE == 'nu':
def rho_tot(T_gam,T_nue,T_numu,MDM):
return rho_gam(T_gam) + rho_e(T_gam) + rho_DM(T_nue,MDM) + rho_nu(T_nue) + 2*rho_nu(T_numu) - P_int(T_gam) + T_gam*dP_intdT(T_gam)
def p_tot(T_gam,T_nue,T_numu,MDM):
return 1./3. * rho_gam(T_gam) + p_e(T_gam) + p_DM(T_nue,MDM) + 1./3. * rho_nu(T_nue) + 1./3. * 2*rho_nu(T_numu) + P_int(T_gam)
def Hubble(T_gam,T_nue,T_numu,MDM):
return FAC * (rho_tot(T_gam,T_nue,T_numu,MDM)*8*np.pi/(3*Mpl**2))**0.5
# Neff Definition
def Neff_func(T_gam,T_nue,T_numu):
return 8./7.*(11./4.)**(4./3.)* (rho_nu(T_nue) + 2*rho_nu(T_numu)) / rho_gam(T_gam)
# Suppression of the rates as a result of a non-negligible electron mass
f_nue_s = interp1d(np.loadtxt("SM_Rates/nue_scatt.dat")[:,0],np.loadtxt("SM_Rates/nue_scatt.dat")[:,1],kind='linear')
f_numu_s = interp1d(np.loadtxt("SM_Rates/numu_scatt.dat")[:,0],np.loadtxt("SM_Rates/numu_scatt.dat")[:,1],kind='linear')
f_nue_a = interp1d(np.loadtxt("SM_Rates/nue_ann.dat")[:,0],np.loadtxt("SM_Rates/nue_ann.dat")[:,1],kind='linear')
f_numu_a = interp1d(np.loadtxt("SM_Rates/numu_ann.dat")[:,0],np.loadtxt("SM_Rates/numu_ann.dat")[:,1],kind='linear')
##Uncomment in order to remove the effect of m_e in the rates.
#f_nue_s, f_numu_s, f_nue_a, f_numu_a = lambda T : 1, lambda T : 1, lambda T : 1, lambda T : 1
def Ffunc_nue_e(T1,T2):
return 32* 0.884 *(T1**9 - T2**9) * f_nue_a(T1) + 56 * 0.829 * f_nue_s(T1) *T1**4*T2**4*(T1-T2)
def Ffunc_numu_e(T1,T2):
return 32* 0.884 *(T1**9 - T2**9) * f_numu_a(T1) + 56 * 0.829 * f_numu_s(T1) *T1**4*T2**4*(T1-T2)
def Ffunc(T1,T2):
return 32* 0.884 *(T1**9 - T2**9) + 56* 0.829 *T1**4*T2**4*(T1-T2)
# Energy Transfer Rates
def DeltaRho_nue(T_gam,T_nue,T_numu):
return FAC * GF**2/np.pi**5 * ( 4* (geL**2 + geR**2) * Ffunc_nue_e(T_gam,T_nue) + 2*Ffunc(T_numu,T_nue) )
def DeltaRho_numu(T_gam,T_nue,T_numu):
return FAC * GF**2/np.pi**5 * ( 4* (gmuL**2 + gmuR**2) * Ffunc_numu_e(T_gam,T_numu) - Ffunc(T_numu,T_nue) )
# Energy Transfer Rate as induced by WIMP annihilations
if MODE == 'e':
def DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR):
sigmav = 2.57*1e-9*1e-6 # 3*10^{-26} cm^3/s in MeV^{-2}, matches equations 4.7 and 4.8 of 1812.05605
return - gDM**2 * FAC * sigmav * BR * MDM**5/(4*np.pi**4) * (T_gam**2 * kv(2,MDM/T_gam)**2 - T_numu**2 * kv(2,MDM/T_numu)**2 )
if MODE == 'nu':
def DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR):
sigmav = 2.57*1e-9*1e-6 # 3*10^{-26} cm^3/s in MeV^{-2}, matches equations 4.7 and 4.8 of 1812.05605
return gDM**2 * FAC * sigmav * BR * MDM**5/(4*np.pi**4) * (T_gam**2 * kv(2,MDM/T_gam)**2 - T_numu**2 * kv(2,MDM/T_numu)**2 )
# Temperature Evolution Equations
if MODE == 'e':
def dTnu_dt(T_gam,T_nue,T_numu,MDM):
return -(Hubble(T_gam,T_nue,T_nue,MDM) * (3 * 4 * rho_nu(T_nue) ) - (2*DeltaRho_numu(T_gam,T_nue,T_numu) + DeltaRho_nue(T_gam,T_nue,T_numu)) + DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR) )/ (3*drho_nudT(T_nue) )
def dTgam_dt(T_gam,T_nue,T_numu,MDM):
return -(Hubble(T_gam,T_nue,T_numu,MDM)*( 4*rho_gam(T_gam) + 3*(rho_e(T_gam)+p_e(T_gam)) + 3*(rho_DM(T_gam,MDM)+p_DM(T_gam,MDM)) + 3 * T_gam * dP_intdT(T_gam)) + DeltaRho_nue(T_gam,T_nue,T_numu) + 2*DeltaRho_numu(T_gam,T_nue,T_numu) - DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR))/( drho_gamdT(T_gam) + drho_edT(T_gam) + drho_DMdT(T_gam,MDM) + T_gam * d2P_intdT2(T_gam) )
if MODE == 'nu':
def dTnu_dt(T_gam,T_nue,T_numu,MDM):
return -(Hubble(T_gam,T_nue,T_nue,MDM) * (3 * 4 * rho_nu(T_nue) + 3*(rho_DM(T_nue,MDM)+p_DM(T_nue,MDM))) - (2*DeltaRho_numu(T_gam,T_nue,T_numu) + DeltaRho_nue(T_gam,T_nue,T_numu) + DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR)))/ (3*drho_nudT(T_nue) + drho_DMdT(T_nue,MDM) )
def dTgam_dt(T_gam,T_nue,T_numu,MDM):
return -(Hubble(T_gam,T_nue,T_nue,MDM)*( 4*rho_gam(T_gam) + 3*(rho_e(T_gam)+p_e(T_gam)) + 3 * T_gam * dP_intdT(T_gam)) + DeltaRho_nue(T_gam,T_nue,T_numu) + 2*DeltaRho_numu(T_gam,T_nue,T_numu) + DeltaRho_DM(T_gam,T_nue,T_numu,MDM,BR) )/( drho_gamdT(T_gam) + drho_edT(T_gam) + T_gam * d2P_intdT2(T_gam) )
def dT_totdt(vec,t,MDM):
T_gam, T_nu = vec
return [dTgam_dt(T_gam,T_nu,T_nu,MDM),dTnu_dt(T_gam,T_nu,T_nu,MDM)]
#Start the integration at a common temperature of 20 MeV, which corresponds to t ~ 2*10^{-3} s
T0 = 10.0
t0 = 1./(2*Hubble(T0,T0,T0,MDM))
# Finish the calculation at t = 5*10^4 seconds, which will correspond to T ~ 5*10^{-3} MeV
t_max = 5e4
# Calculate
tvec = np.logspace(np.log10(t0),np.log10(t_max),num=300)
sol = odeint(dT_totdt, [T0,T0], tvec, args=(MDM,), rtol = 1e-8, atol= 1e-8)
# Display Results
print("Neff = ", round(Neff_func(sol[-1,0],sol[-1,1],sol[-1,1]),5))
print("Tgam/Tnu = ", round(sol[-1,0]/sol[-1,1],5))
# Output Results
with open("Results/WIMPS/"+file+"_Neff.dat", 'w') as f:
f.write(str(round(Neff_func(sol[-1,0],sol[-1,1],sol[-1,1]),5)))
# Calculate some Thermodynamic quantities
rho_vec, p_vec = np.zeros(len(tvec)),np.zeros(len(tvec))
for i in range(len(tvec)):
rho_vec[i], p_vec[i] = rho_tot(sol[i,0], sol[i,1], sol[i,1],MDM), p_tot(sol[i,0], sol[i,1], sol[i,1],MDM)
MDM, gDM, SPIN, MODE, BR
# Store results
print("Thermodynamics output to : Results/WIMPS/"+file+".dat")
with open("Results/WIMPS/"+file+".dat", 'w') as f:
f.write("# WIMP coupled to "+MODE+". With "+"m = "+str(round(MDM,2))+", g = "+str(int(gDM))+", SPIN = "+SPIN+", and BR to other = "+"{:.2E}".format(BR)+"\n")
f.write("# Neff = "),f.write("{:<12}".format(round(Neff_func(sol[-1,0],sol[-1,1],sol[-1,1]),5))),f.write("\n")
f.write("# Tgam/Tnu = "),f.write("{:<12}".format(round(sol[-1,0]/sol[-1,1],5))),f.write("\n")
f.write("{:<13}".format("#T_gamma (MeV)")), f.write("{:<13}".format("T_gam/T_nu")),f.write("{:<14}".format("R_tot/T_gam^4 ")),f.write("{:<14}".format("P_tot/T_gam^4"+"\n"))
for i in range(len(tvec)):
f.write("{:.7E}".format(1e3*sol[i,0])),f.write(" "), f.write("{:<12}".format(round(sol[i,0]/sol[i,1],6))), f.write("{:.7E}".format(rho_vec[i]/sol[i,0]**4)),f.write(" "),f.write("{:.7E}".format(p_vec[i]/sol[i,0]**4)),f.write("\n")