.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/astrophysics/astro_dark_matter.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. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_astrophysics_astro_dark_matter.py: Dark matter spatial and spectral models ======================================= Convenience methods for dark matter high level analyses. Introduction ------------ Gammapy has some convenience methods for dark matter analyses in `gammapy.astro.darkmatter`. These include J-Factor computation and calculation the expected gamma flux for a number of annihilation channels. They are presented in this notebook. The basic concepts of indirect dark matter searches, however, are not explained. So this is aimed at people who already know what the want to do. A good introduction to indirect dark matter searches is given for example `here `__ (Chapter 1 and 5). .. GENERATED FROM PYTHON SOURCE LINES 23-28 Setup ----- As always, we start with some setup for the notebook, and with imports. .. GENERATED FROM PYTHON SOURCE LINES 28-47 .. code-block:: Python import numpy as np import astropy.units as u from astropy.coordinates import SkyCoord from regions import CircleSkyRegion, RectangleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt from matplotlib.colors import LogNorm from gammapy.astro.darkmatter import ( DarkMatterAnnihilationSpectralModel, DarkMatterDecaySpectralModel, JFactory, PrimaryFlux, profiles, ) from gammapy.maps import WcsGeom, WcsNDMap .. GENERATED FROM PYTHON SOURCE LINES 48-56 Profiles -------- The following dark matter profiles are currently implemented. Each model can be scaled to a given density at a certain distance. These parameters are controlled by `~gammapy.astro.darkmatter.profiles.DMProfile.LOCAL_DENSITY` and `~gammapy.astro.darkmatter.profiles.DMProfile.DISTANCE_GC` .. GENERATED FROM PYTHON SOURCE LINES 56-74 .. code-block:: Python profiles.DMProfile.__subclasses__() for profile in profiles.DMProfile.__subclasses__(): p = profile() p.scale_to_local_density() radii = np.logspace(-3, 2, 100) * u.kpc plt.plot(radii, p(radii), label=p.__class__.__name__) plt.loglog() plt.axvline(8.5, linestyle="dashed", color="black", label="local density") plt.legend() plt.show() print("LOCAL_DENSITY:", profiles.DMProfile.LOCAL_DENSITY) print("DISTANCE_GC:", profiles.DMProfile.DISTANCE_GC) .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_001.png :alt: astro dark matter :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none LOCAL_DENSITY: 0.3 GeV / cm3 DISTANCE_GC: 8.33 kpc .. GENERATED FROM PYTHON SOURCE LINES 75-82 J Factors --------- There are utilities to compute J-Factor maps that can serve as a basis to compute J-Factors for certain regions. In the following we compute a J-Factor annihilation map for the Galactic Centre region .. GENERATED FROM PYTHON SOURCE LINES 82-85 .. code-block:: Python profile = profiles.NFWProfile(r_s=20 * u.kpc) .. GENERATED FROM PYTHON SOURCE LINES 86-87 Adopt standard values used in H.E.S.S. .. GENERATED FROM PYTHON SOURCE LINES 87-104 .. code-block:: Python profiles.DMProfile.DISTANCE_GC = 8.5 * u.kpc profiles.DMProfile.LOCAL_DENSITY = 0.39 * u.Unit("GeV / cm3") profile.scale_to_local_density() position = SkyCoord(0.0, 0.0, frame="galactic", unit="deg") geom = WcsGeom.create(binsz=0.05, skydir=position, width=3.0, frame="galactic") jfactory = JFactory(geom=geom, profile=profile, distance=profiles.DMProfile.DISTANCE_GC) jfact = jfactory.compute_jfactor() jfact_map = WcsNDMap(geom=geom, data=jfact.value, unit=jfact.unit) plt.figure() ax = jfact_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title(f"J-Factor [{jfact_map.unit}]") .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_002.png :alt: J-Factor [GeV2 / cm5] :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'J-Factor [GeV2 / cm5]') .. GENERATED FROM PYTHON SOURCE LINES 105-106 1 deg circle usually used in H.E.S.S. analyses without the +/- 0.3 deg band around the plane .. GENERATED FROM PYTHON SOURCE LINES 106-118 .. code-block:: Python sky_reg = CircleSkyRegion(center=position, radius=1 * u.deg) pix_reg = sky_reg.to_pixel(wcs=geom.wcs) pix_reg.plot(ax=ax, facecolor="none", edgecolor="red", label="1 deg circle") sky_reg_rec = RectangleSkyRegion(center=position, height=0.6 * u.deg, width=2 * u.deg) pix_reg_rec = sky_reg_rec.to_pixel(wcs=geom.wcs) pix_reg_rec.plot(ax=ax, facecolor="none", edgecolor="orange", label="+/- 0.3 deg band") plt.legend() plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_003.png :alt: astro dark matter :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-121 Note the value quoted by H.E.S.S. in `this paper `__ is :math:`2.67\times10^{21}` .. GENERATED FROM PYTHON SOURCE LINES 121-135 .. code-block:: Python total_jfact = ( pix_reg.to_mask().multiply(jfact).sum() - pix_reg_rec.to_mask().multiply(jfact).sum() ) total_jfact = ( pix_reg.to_mask().multiply(jfact).sum() - pix_reg_rec.to_mask().multiply(jfact).sum() ) print( "J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a " f"{profile.__class__.__name__} is {total_jfact:.3g}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a NFWProfile is 2.36e+21 GeV2 / cm5 .. GENERATED FROM PYTHON SOURCE LINES 136-137 The J-Factor can also be computed for dark matter decay .. GENERATED FROM PYTHON SOURCE LINES 137-151 .. code-block:: Python jfactory = JFactory( geom=geom, profile=profile, distance=profiles.DMProfile.DISTANCE_GC, annihilation=False, ) jfact_decay = jfactory.compute_jfactor() jfact_map = WcsNDMap(geom=geom, data=jfact_decay.value, unit=jfact_decay.unit) plt.figure() ax = jfact_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title(f"J-Factor [{jfact_map.unit}]") .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_004.png :alt: J-Factor [GeV / cm2] :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'J-Factor [GeV / cm2]') .. GENERATED FROM PYTHON SOURCE LINES 152-153 1 deg circle usually used in H.E.S.S. analyses without the +/- 0.3 deg band around the plane .. GENERATED FROM PYTHON SOURCE LINES 153-178 .. code-block:: Python sky_reg = CircleSkyRegion(center=position, radius=1 * u.deg) pix_reg = sky_reg.to_pixel(wcs=geom.wcs) pix_reg.plot(ax=ax, facecolor="none", edgecolor="red", label="1 deg circle") sky_reg_rec = RectangleSkyRegion(center=position, height=0.6 * u.deg, width=2 * u.deg) pix_reg_rec = sky_reg_rec.to_pixel(wcs=geom.wcs) pix_reg_rec.plot(ax=ax, facecolor="none", edgecolor="orange", label="+/- 0.3 deg band") plt.legend() plt.show() total_jfact_decay = ( pix_reg.to_mask().multiply(jfact_decay).sum() - pix_reg_rec.to_mask().multiply(jfact_decay).sum() ) total_jfact_decay = ( pix_reg.to_mask().multiply(jfact_decay).sum() - pix_reg_rec.to_mask().multiply(jfact_decay).sum() ) print( "J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a " f"{profile.__class__.__name__} is {total_jfact_decay:.3g}" ) .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_005.png :alt: astro dark matter :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a NFWProfile is 1.14e+20 GeV / cm2 .. GENERATED FROM PYTHON SOURCE LINES 179-186 Gamma-ray spectra at production ------------------------------- The gamma-ray spectrum per annihilation is a further ingredient for a dark matter analysis. The following annihilation channels are supported. For more info see https://arxiv.org/pdf/1012.4515.pdf .. GENERATED FROM PYTHON SOURCE LINES 186-214 .. code-block:: Python fluxes = PrimaryFlux(mDM="1 TeV", channel="eL") print(fluxes.allowed_channels) fig, axes = plt.subplots(4, 1, figsize=(4, 16)) mDMs = [0.01, 0.1, 1, 10] * u.TeV for mDM, ax in zip(mDMs, axes): fluxes.mDM = mDM ax.set_title(rf"m$_{{\mathrm{{DM}}}}$ = {mDM}") ax.set_yscale("log") ax.set_ylabel("dN/dE") for channel in ["tau", "mu", "b", "Z"]: fluxes = PrimaryFlux(mDM=mDM, channel=channel) fluxes.channel = channel fluxes.plot( energy_bounds=[mDM / 100, mDM], ax=ax, label=channel, yunits=u.Unit("1/GeV"), ) axes[0].legend() plt.subplots_adjust(hspace=0.9) plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_006.png :alt: m$_{\mathrm{DM}}$ = 0.01 TeV, m$_{\mathrm{DM}}$ = 0.1 TeV, m$_{\mathrm{DM}}$ = 1.0 TeV, m$_{\mathrm{DM}}$ = 10.0 TeV :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ['eL', 'eR', 'e', 'muL', 'muR', 'mu', 'tauL', 'tauR', 'tau', 'q', 'c', 'b', 't', 'WL', 'WT', 'W', 'ZL', 'ZT', 'Z', 'g', 'gamma', 'h', 'nu_e', 'nu_mu', 'nu_tau', 'V->e', 'V->mu', 'V->tau'] .. GENERATED FROM PYTHON SOURCE LINES 215-220 Flux maps for annihilation -------------------------- Finally flux maps can be produced like this: .. GENERATED FROM PYTHON SOURCE LINES 220-238 .. code-block:: Python channel = "Z" massDM = 10 * u.TeV diff_flux = DarkMatterAnnihilationSpectralModel(mass=massDM, channel=channel) int_flux = ( jfact * diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV) ).to("cm-2 s-1") flux_map = WcsNDMap(geom=geom, data=int_flux.value, unit="cm-2 s-1") plt.figure() ax = flux_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title( f"Flux [{int_flux.unit}]\n m$_{{DM}}$={fluxes.mDM.to('TeV')}, channel={fluxes.channel}" ) plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_007.png :alt: Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 239-244 Flux maps for decay ------------------- Finally flux maps for decay can be produced like this: .. GENERATED FROM PYTHON SOURCE LINES 244-260 .. code-block:: Python channel = "Z" massDM = 10 * u.TeV diff_flux = DarkMatterDecaySpectralModel(mass=massDM, channel=channel) int_flux = ( jfact_decay * diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV) ).to("cm-2 s-1") flux_map = WcsNDMap(geom=geom, data=int_flux.value, unit="cm-2 s-1") plt.figure() ax = flux_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title( f"Flux [{int_flux.unit}]\n m$_{{DM}}$={fluxes.mDM.to('TeV')}, channel={fluxes.channel}" ) plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_008.png :alt: Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 25.697 seconds) .. _sphx_glr_download_tutorials_astrophysics_astro_dark_matter.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v2.0?urlpath=lab/tree/notebooks/2.0/tutorials/astrophysics/astro_dark_matter.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: astro_dark_matter.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: astro_dark_matter.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: astro_dark_matter.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_