.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_gallery/plot_relic_density.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_gallery_plot_relic_density.py: Relic density calculations ========================== Calculate the relic density of the model. .. GENERATED FROM PYTHON SOURCE LINES 7-12 .. code-block:: default import matplotlib.pyplot as plt import numpy as np from singletscalar_dm import * .. GENERATED FROM PYTHON SOURCE LINES 13-15 In order to calculate the relic density it is possible to use the function `interpolate_Omega`, which uses the code DRAKE near the Higgs resonance and MicrOMEGAs elsewhere. For example, let's compute the relic density, expressed as :math:`\Omega h^2` for the QCDA model. .. GENERATED FROM PYTHON SOURCE LINES 15-21 .. code-block:: default mS = 50 # GeV lambda_hs = 0.1 omegah2 = interpolate_Omega(mS,lambda_hs,'QCDA',True) print(omegah2) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.12801745167589393 .. GENERATED FROM PYTHON SOURCE LINES 22-23 In the following, we compute the relic density using only MicrOMEGAs. .. GENERATED FROM PYTHON SOURCE LINES 23-29 .. code-block:: default mS = 50 # GeV lambda_hs = 0.1 omegah2_micromegas = interpolate_Omega_MicrOMEGAs(mS,lambda_hs) print(omegah2_micromegas) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.12796691577839414 .. GENERATED FROM PYTHON SOURCE LINES 30-32 As expected the results obtained with the functions `interpolate_Omega` and `interpolate_Omega_MicrOMEGAs` are similar for the mass point tested above. However, if we test masses close to the Higgs resonance the differences emerge: .. GENERATED FROM PYTHON SOURCE LINES 32-40 .. code-block:: default mS = 60 # GeV lambda_hs = 0.001 omegah2 = interpolate_Omega(mS,lambda_hs,'QCDA',True) omegah2_micromegas = interpolate_Omega_MicrOMEGAs(mS,lambda_hs) print(omegah2) print(omegah2_micromegas) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.06396689881423608 0.047495662241396217 .. GENERATED FROM PYTHON SOURCE LINES 41-43 The calculation of the relic density close to the resonance with DRAKE has been done for `lambda_hs` values between 1 and two orders of magnitude below and above the value that gives the correct relic abundance. If one tries to calculate the relic density with DRAKE beyond this region a warning will be printed. .. GENERATED FROM PYTHON SOURCE LINES 43-49 .. code-block:: default mS = 60 # GeV lambda_hs = 0.1 omegah2 = interpolate_Omega(mS,lambda_hs,'QCDA',True) print(omegah2) .. rst-class:: sphx-glr-script-out .. code-block:: none Warning, extrapolating..... The relic density with DRAKE has been calculated for lambda_HS between 7.32e-05 0.01464 Following Value of lambda_HS is given 0.1 46.92619215842866 .. GENERATED FROM PYTHON SOURCE LINES 50-53 Note that the precomputed results with DRAKE have been done to cover the region which gives the right relic abundance, so investigating region beyond the tested one is not recommended. It is also possible to find the :math:`\lambda_{HS}` value which provides a give relic density. In order to do that you can use the function `interpolate_lambda` .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: default mS = 60 # GeV Omegah2 = 0.1 lambda_hs = interpolate_lambda(mS,Omegah2,'QCDB',True) print(lambda_hs) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0013339697597695931 .. GENERATED FROM PYTHON SOURCE LINES 60-62 Note that for dark matter masses below 70 GeV and :math:`\lambda_{HS}` values larger than 1 the relic density have a turning point. See the following plot. .. GENERATED FROM PYTHON SOURCE LINES 62-69 .. code-block:: default mass = 45 lambda_vec = np.logspace(-4.,2.,50) omegah2_vec = np.zeros(len(lambda_vec)) for t in range(len(lambda_vec)): omegah2_vec[t] = interpolate_Omega_MicrOMEGAs(mass,lambda_vec[t]) .. GENERATED FROM PYTHON SOURCE LINES 70-85 .. code-block:: default fig = plt.figure(figsize=(8,6)) plt.plot(lambda_vec,omegah2_vec,lw=1.5,ls='-',color='black',label='$m_S=45$ GeV') plt.ylabel(r'$\Omega h^2$', fontsize=18) plt.xlabel(r'$\lambda_{HS}$', fontsize=18) plt.axis([1e-4,100,1e-4,1e2]) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid(True) plt.yscale('log') plt.xscale('log') plt.legend(loc=4,prop={'size':14},numpoints=1, scatterpoints=1, ncol=1) fig.tight_layout(pad=0.5) plt.show() .. image-sg:: /examples_gallery/images/sphx_glr_plot_relic_density_001.png :alt: plot relic density :srcset: /examples_gallery/images/sphx_glr_plot_relic_density_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 86-89 This implies that for :math:`m_S < 70` GeV there could be two possible for the :math:`\lambda_{HS}` which gives a value of the relic density. The code prints a warning if this is the case. See example below. .. GENERATED FROM PYTHON SOURCE LINES 89-95 .. code-block:: default mS = 30 # GeV Omegah2 = 0.001 lambda_hs = interpolate_lambda(mS,Omegah2,'QCDB',True) print(lambda_hs) .. rst-class:: sphx-glr-script-out .. code-block:: none Warning, extrapolating. For this mass pick a range of Omegah^2 between 0.001075 984200.0 Following Value of Omegah^2 is given 0.001 0 .. GENERATED FROM PYTHON SOURCE LINES 96-102 .. code-block:: default mS = 70 # GeV lambda_hs = 0.80952387 lambda_hs = interpolate_Omega(mS,lambda_hs,'QCDA',True) print(lambda_hs) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0009975472734606318 .. GENERATED FROM PYTHON SOURCE LINES 103-104 In the following, we will generate a plot with the values of :math:`m_S` and :math:`\lambda_{HS}` which provide 100% or 30% of the relic density, showing the calculations performed with the QCDA and QCDB models and with MicrOMEGAs. .. GENERATED FROM PYTHON SOURCE LINES 104-117 .. code-block:: default mass_vec = massz_vec lambda_QCDA_100_vec = np.zeros(len(mass_vec)) lambda_QCDB_100_vec = np.zeros(len(mass_vec)) lambda_QCDB_30_vec = np.zeros(len(mass_vec)) lambda_Micro_100_vec = np.zeros(len(mass_vec)) Omegah2 = 0.120 for t in range(len(mass_vec)): lambda_QCDA_100_vec[t] = interpolate_lambda(mass_vec[t],Omegah2,'QCDA',False) lambda_QCDB_100_vec[t] = interpolate_lambda(mass_vec[t],Omegah2,'QCDB',False) lambda_QCDB_30_vec[t] = interpolate_lambda(mass_vec[t],0.1*Omegah2,'QCDB',False) lambda_Micro_100_vec[t] = interpolate_lambda_MicrOMEGAs(mass_vec[t],Omegah2,False) .. rst-class:: sphx-glr-script-out .. code-block:: none Warning, mass should be between 2.0 20000.0 Following Value of mass is given 2.0 Warning, mass should be between 2.0 20000.0 Following Value of mass is given 2.0 Warning, mass should be between 2.0 20000.0 Following Value of mass is given 2.0 Warning, mass should be between 2.0 20000.0 Following Value of mass is given 2.0 .. GENERATED FROM PYTHON SOURCE LINES 118-153 .. code-block:: default fig = plt.figure(figsize=(8,6)) plt.plot(mass_vec,lambda_QCDA_100_vec,lw=1.5,ls='--',color='black',label='DRAKE QCDA, $\Omega h^2=0.12$') plt.plot(mass_vec,lambda_QCDB_100_vec,lw=1.5,ls='-',color='blue',label='DRAKE QCDB, $\Omega h^2=0.12$') plt.plot(mass_vec,lambda_QCDB_30_vec,lw=1.5,ls=':',color='red',label='DRAKE QCDB, $\Omega h^2=0.012$') plt.plot(mass_vec,lambda_Micro_100_vec,lw=1.5,ls='-.',color='green', label=r'MicroOMEGA, $\Omega h^2=0.12$') plt.ylabel(r'$\lambda_{HS}$', fontsize=18) plt.xlabel(r'$m_{S}$ [GeV]', fontsize=18) plt.axis([10.,1000,1e-4,2e0]) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid(True) plt.yscale('log') plt.xscale('log') plt.legend(loc=4,prop={'size':14},numpoints=1, scatterpoints=1, ncol=1) fig.tight_layout(pad=0.5) plt.show() fig = plt.figure(figsize=(8,6)) plt.plot(mass_vec,lambda_QCDA_100_vec,lw=1.5,ls='--',color='black',label='DRAKE QCDA, $\Omega h^2=0.12$') plt.plot(mass_vec,lambda_QCDB_100_vec,lw=1.5,ls='-',color='blue',label='DRAKE QCDB, $\Omega h^2=0.12$') plt.plot(mass_vec,lambda_QCDB_30_vec,lw=1.5,ls=':',color='red',label='DRAKE QCDB, $\Omega h^2=0.012$') plt.plot(mass_vec,lambda_Micro_100_vec,lw=1.5,ls='-.',color='green', label=r'MicroOMEGA, $\Omega h^2=0.12$') plt.ylabel(r'$\lambda_{HS}$', fontsize=18) plt.xlabel(r'$m_{S}$ [GeV]', fontsize=18) plt.axis([40.,80,1e-4,1e0]) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid(True) plt.yscale('log') plt.xscale('log') plt.legend(loc=3,prop={'size':14},numpoints=1, scatterpoints=1, ncol=1) fig.tight_layout(pad=0.5) plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_gallery/images/sphx_glr_plot_relic_density_002.png :alt: plot relic density :srcset: /examples_gallery/images/sphx_glr_plot_relic_density_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples_gallery/images/sphx_glr_plot_relic_density_003.png :alt: plot relic density :srcset: /examples_gallery/images/sphx_glr_plot_relic_density_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 154-155 The parameters :math:`m_S` and :math:`\lambda_{HS}` which provide the right relic abundance are reported in the file `Omega_MicroOMEGAs_DRAKE_QCDB_QCDA.dat`, which can be imported through the function `import_data_file`. .. GENERATED FROM PYTHON SOURCE LINES 155-176 .. code-block:: default omega_drake_micromegas = import_data_file('Omega_MicroOMEGAs_DRAKE_QCDB_QCDA.dat') table = np.loadtxt(omega_drake_micromegas) mass_RD = table[:,0] lambda_RD_FBQCDA = table[:,1] lambda_RD_FBQCDB = table[:,2] fig = plt.figure(figsize=(8,6)) plt.fill_between(mass_RD,lambda_RD_FBQCDA,lambda_RD_FBQCDB,color='cyan',alpha=0.5) plt.plot(mass_RD,lambda_RD_FBQCDA,lw=2.5,ls='-.',color='black', label=r'DRAKE fBE, $QCD_A$') plt.plot(mass_RD,lambda_RD_FBQCDB,lw=2.5,ls=':',color='grey', label=r'DRAKE fBE, $QCD_B$') plt.ylabel(r'$\lambda_{HS}$', fontsize=18) plt.xlabel(r'$m_{S}$ [GeV]', fontsize=18) plt.axis([50.,70,1e-4,1e-1]) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid(True) plt.yscale('log') plt.xscale('linear') plt.legend(loc=4,prop={'size':12},numpoints=1, scatterpoints=1, ncol=1) fig.tight_layout(pad=0.5) plt.show() .. image-sg:: /examples_gallery/images/sphx_glr_plot_relic_density_004.png :alt: plot relic density :srcset: /examples_gallery/images/sphx_glr_plot_relic_density_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 11 minutes 37.779 seconds) .. _sphx_glr_download_examples_gallery_plot_relic_density.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_relic_density.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_relic_density.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_