.. 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-102 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 103-104 Plot the J-factor map .. GENERATED FROM PYTHON SOURCE LINES 104-121 .. code-block:: Python plt.figure() ax = jfact_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title(f"J-Factor [{jfact_map.unit}]") # 1 deg circle usually used in H.E.S.S. analyses without the +/- 0.3 deg band around the plane 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_002.png :alt: J-Factor [GeV2 / cm5] :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 122-124 Note the value quoted by H.E.S.S. in `this paper `__ is :math:`2.67\times10^{21}` .. GENERATED FROM PYTHON SOURCE LINES 124-138 .. 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 139-140 The J-Factor can also be computed for dark matter decay .. GENERATED FROM PYTHON SOURCE LINES 140-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) .. GENERATED FROM PYTHON SOURCE LINES 152-153 Plot the J-factor map .. GENERATED FROM PYTHON SOURCE LINES 153-169 .. code-block:: Python plt.figure() ax = jfact_map.plot(cmap="viridis", norm=LogNorm(), add_cbar=True) plt.title(f"J-Factor [{jfact_map.unit}]") # 1 deg circle usually used in H.E.S.S. analyses without the +/- 0.3 deg band around the plane 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: J-Factor [GeV / cm2] :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 170-183 .. code-block:: Python 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}" ) .. 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 184-191 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 191-220 .. code-block:: Python fluxes = PrimaryFlux(mDM="1 TeV", channel="eL") print(fluxes.allowed_channels) fig, axes = plt.subplots(2, 2, figsize=(10, 9)) axes = axes.flatten() 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() fig.tight_layout() plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_004.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_004.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 221-226 Flux maps for annihilation -------------------------- Finally flux maps can be produced like this: .. GENERATED FROM PYTHON SOURCE LINES 226-244 .. 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_005.png :alt: Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 245-250 Flux maps for decay ------------------- Finally flux maps for decay can be produced like this: .. GENERATED FROM PYTHON SOURCE LINES 250-266 .. 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_006.png :alt: Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z :srcset: /tutorials/astrophysics/images/sphx_glr_astro_dark_matter_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 26.949 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/main?urlpath=lab/tree/notebooks/dev/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 `_