LoganAMorrison
commented
5 years ago hazma.relic_density.relic_density
Okay, I think I've done the best I can without cythonizing. We can now compute the relic densities of any model. As long as the model implements model.mx and
self.annihilation_cross_sections (which must return a dictionary with key total), then the function relic_density will compute the relic density. Alternatively, the model can just implement a function called model.thermal_cross_section and it will still work (still needs model.mx though.)
Semi-Analytic + Boltzmann Methods
I've implemented the calculation in two different ways:
- full numeric integration of the Boltzmann equation and
- a semi-analytic approximation using similar results as arXiv:1204.3622v3.
The semi-analytic calculation takes about 0.5s while the full integration of the Boltzmann equation takes ~10 seconds (yikes...) They both give very similar answers so I made the default use the semi-analytical calculation.
I added specific functions in the scalar and vector mediator model for the thermal cross-sections which only differer by the built-in methods by adding specific breakpoints for the integration.
We can make this calculation even faster if we cythonize the annihilation cross-section functions. I'm sure we could bring the calculation down to ms levels.
Usage
For the scalar and vector mediator models, you can use:
from hazma.scalar_mediator import HiggsPortal
from hazma.vector_mediator import KineticMixing
from hazma.relic_density import relic_density
hp = HiggsPortal(100.0, 123.5, 1.0, 1e-3)
km = KineticMixing(100.0, 125.2, 1.0, 1e-4)
# semi-analytical
relic_density(hp) # 0.11847057586537633
relic_density(km) # 0.11922204452307583
# integration of Boltzmann
relic_density(hp, semi_analytic=False) # 0.12104612783073646
relic_density(km, semi_analytic=False) # 0.13167106960636926One can also make their own model to compute the relic density:
from hazma.relic_density import relic_density
class MyModel(object):
def __init__(self, mx, sigmav):
self.mx = mx
self.sigmav = sigmav
def thermal_cross_section(self, x):
"""x = mx / T"""
return self.sigmav
# These give the correct relic density for DM
mx = 104.74522360006331e3 # ~104 GeV in MeV
sigmav = 1.7597967261428258e-15
model = MyModel(mx, sigmav)
relic_density(model, semi_analytic=True) # 0.11444106905299978
relic_density(model, semi_analytic=False) # 0.11963483643866808These computations are significally faster since no thermal integrals need to be computed (4 ms for semi-analytic and 40 ms for solving Boltzmann equation.)
Tests
I added tests verifying that both the semi-analytical method and Botlzman integration work using the following points:
- mx1, sigmav1 = 10.313897683787216e3, 1.966877938634266e-15
- mx2, sigmav2 = 104.74522360006331e3, 1.7597967261428258e-15
- mx3, sigmav3 = 1063.764854316313e3, 1.837766552668581e-15
- mx4, sigmav4 = 10000.0e3, 1.8795945459427076e-15
Each of these should give a relic density of Omega*h^2 ~ 0.119.
simps and similar functions don't work?
Lastly, I had a lot of issues trying to use different quadrature rules to compute the thermal cross-section integrals. I tried SciPy's simps, romb, and quadrature as well as NumPy's laggauss fo Gauss-Laguerre quadrature. All of these gave an answer that was two orders of magnitude larger than what SciPy's quad gave. I have no idea what the issue is. Add the end, I just used SciPy's quad. @adam-coogan any idea why there would be such a large discrepancy? Here is a MWE of the discrepancy:
from hazma.scalar_mediator import HiggsPortal
from scipy.integrate import simps, quad
from scipy.special import kn, k1
def thermal_cross_section(x, model, use_quad=True):
# need this to avoid divide by zero warnings
if x > 300.0:
return 0.0
pf = x / (2.0 * kn(2, x))**2
def integrand(z):
sig = model.annihilation_cross_sections(model.mx * z)['total']
kernal = z**2 * (z**2 - 4.0) * k1(x * z)
return sig * kernal
if use_quad:
return pf * quad(integrand, 2.0, np.inf)[0]
else:
zs = np.linspace(2.0, 200.0, 1000)
integrands = integrand(zs)
return pf * simps(integrands, zs)
# This doesn't work
hp = HiggsPortal(100.0, 1000.0, 1.0, 1e-3)
print(thermal_cross_section(1.0, hp, use_quad=False)) # 4.972292358963628e-15
print(thermal_cross_section(1.0, hp, use_quad=True)) # 5.252226660601945e-17
# but this works?
hp = HiggsPortal(100.0, 150.0, 1.0, 1e-3)
print(thermal_cross_section(1.0, hp, use_quad=False)) # 2.469852823226098e-07
print(thermal_cross_section(1.0, hp, use_quad=True)) # 2.519329618702442e-07On second thought, I wonder if it's quad that's wrong? Maybe we should come up with an example where we can get an analytic result and compare. But I have a feeling that it's the resonances/thresholds that give the issue. Here is a plot of what the integrand looks like for the scalar:
from hazma.scalar_mediator import HiggsPortal
from scipy.integrate import simps, quad
from scipy.special import k1
import matplotlib.pyplot as plt
def thermal_cross_section_integrand(z, x, model):
sig = model.annihilation_cross_sections(model.mx * z)['total']
kernal = z**2 * (z**2 - 4.0) * k1(x * z)
return sig * kernal
hp = HiggsPortal(100.0, 1000.0, 1.0, 1e-3)
x = 1.0
zs = np.linspace(2.0, 80.0, 1000)
integrands = thermal_cross_section_integrand(zs, x, hp)
plt.figure(dpi=100)
plt.plot(zs, integrands, lw=2)
plt.vlines(hp.ms / hp.mx, 0.0, 1e-10, ls='--', label='Resonance')
plt.vlines(2.0 * hp.ms / hp.mx, 0.0, 1e-10, ls='-.', label='Threshold')
plt.yscale('log')
plt.ylabel('Thermal Cross Section Integrand', fontsize=14)
plt.xlabel(r'$z = E_{\mathrm{CM}} / m_{\mathrm{DM}}$', fontsize=14)
plt.ylim(1e-35, 1e-10)
plt.legend(fontsize=14)thermal_cross_section_integrand.pdf
Review
@adam-coogan wanna review some of the code to make sure things look alright? Also, do you have any other suggestions before a merge can take place?
adam-coogan
We should add submodules inside the scalar and vector mediator classes to compute the relic density of the dark matter in these frameworks.