Note

Go to the end to download the full example code or to run this example in your browser via Binder.

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).

Setup#

As always, we start with some setup for the notebook, and with imports.

import astropy.units as u

# %matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import SkyCoord
from matplotlib.colors import LogNorm
from regions import CircleSkyRegion, RectangleSkyRegion

from gammapy.astro.darkmatter import (
    ContinuumPrimaryFlux,
    DarkMatterAnnihilationSpectralModel,
    DarkMatterDecaySpectralModel,
    JFactory,
    profiles,
)
from gammapy.maps import WcsGeom, WcsNDMap

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 LOCAL_DENSITY

profile_classes = [
    profiles.NFWProfile,
    profiles.EinastoProfile,
    profiles.IsothermalProfile,
    profiles.BurkertProfile,
    profiles.MooreProfile,
]

profile_labels = {
    profiles.NFWProfile: "NFW",
    profiles.EinastoProfile: "Einasto",
    profiles.IsothermalProfile: "Isothermal",
    profiles.BurkertProfile: "Burkert",
    profiles.MooreProfile: "Moore",
}

for profile_class in profile_classes:
    p = profile_class()
    p.scale_to_local_density(10 * u.kpc)
    radii = np.logspace(-3, 2, 100) * u.kpc
    plt.plot(radii, p(radii), label=profile_labels[profile_class])


plt.loglog()
plt.xlabel("Radius [kpc]")
plt.ylabel(r"Density [GeV cm$^{-3}$]")
plt.axvline(8.33, linestyle="dashed", color="black", label="local density")
plt.legend()
plt.show()

print("LOCAL_DENSITY:", profiles.DMProfile.LOCAL_DENSITY)
astro dark matter
LOCAL_DENSITY: 0.3 GeV / cm3

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

profile = profiles.NFWProfile(r_s=20 * u.kpc)

Adopt standard values used in H.E.S.S.

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=10 * u.kpc, rmax=4 * u.kpc)
jfact = jfactory.compute_jfactor()

jfact_map = WcsNDMap(geom=geom, data=jfact.value, unit=jfact.unit)

Plot the J-factor map

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()
J-Factor [GeV2 / cm5]
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/astrophysics/astro_dark_matter.py:130: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  plt.legend()

Note the value quoted by H.E.S.S. in this paper is \(2.67\times10^{21}\)

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}"
)
total_jfact = u.Quantity(float(total_jfact.value), unit=total_jfact.unit)

J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a NFWProfile is 1.76e+22 GeV2 / cm5

The J-Factor can also be computed for dark matter decay

jfactory = JFactory(
    geom=geom,
    profile=profile,
    distance=10 * u.kpc,
    annihilation=False,
    rmax=4 * u.kpc,
)
jfact_decay = jfactory.compute_jfactor()

jfact_map = WcsNDMap(geom=geom, data=jfact_decay.value, unit=jfact_decay.unit)

Plot the J-factor map

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()
J-Factor [GeV / cm2]
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/astrophysics/astro_dark_matter.py:182: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  plt.legend()

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}"
)

total_jfact_decay = u.Quantity(
    float(total_jfact_decay.value), unit=total_jfact_decay.unit
)

J-factor in 1 deg circle without the +/- 0.3 deg band around GC assuming a NFWProfile is 2.8e+20 GeV / cm2

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

fluxes = ContinuumPrimaryFlux(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 = ContinuumPrimaryFlux(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()
m$_{\mathrm{DM}}$ = 0.01 TeV, m$_{\mathrm{DM}}$ = 0.1 TeV, m$_{\mathrm{DM}}$ = 1.0 TeV, m$_{\mathrm{DM}}$ = 10.0 TeV
['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', 'aZ', 'HZ', 'd', 'u', 's']
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/modeling/models/spectral.py:649: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  ax.plot(energy.center, flux.quantity[:, 0, 0], **kwargs)

Flux maps for annihilation#

Finally flux maps can be produced like this:

channel = "Z"
massDM = 10 * u.TeV

diff_flux = DarkMatterAnnihilationSpectralModel(
    mass=massDM, channel=channel, jfactor=total_jfact
)
int_flux = 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 * jfact / total_jfact).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()
Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/astrophysics/astro_dark_matter.py:249: GammapyDeprecationWarning: The DarkMatterAnnihilationSpectralModel class is deprecated and may be removed in a future version.
        Use DarkMatterSpectralModel instead.
  diff_flux = DarkMatterAnnihilationSpectralModel(
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/astro/darkmatter/spectra.py:895: GammapyDeprecationWarning: "mass" was deprecated in version 2.2 and will be removed in a future version. Use argument "mDM" instead.
  super().__init__(*args, **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/astropy/utils/decorators.py:620: GammapyDeprecationWarning: "jfactor" was deprecated in version 2.2 and will be removed in a future version. Use argument "factor" instead.
  return function(*args, **kwargs)

Flux maps for decay#

Same approach for decay, using the D-factor instead of the J-factor.

channel = "Z"
massDM = 10 * u.TeV

diff_flux_decay = DarkMatterDecaySpectralModel(
    mass=massDM, channel=channel, jfactor=total_jfact_decay
)
int_flux_decay = diff_flux_decay.integral(
    energy_min=0.1 * u.TeV, energy_max=10 * u.TeV
).to("cm-2 s-1")

flux_map_decay = WcsNDMap(
    geom=geom,
    data=(int_flux_decay * jfact_decay / total_jfact_decay).value,
    unit="cm-2 s-1",
)
plt.figure()
ax = flux_map_decay.plot(cmap="viridis", norm=LogNorm(), add_cbar=True)
plt.title(
    f"Flux [{flux_map_decay.unit}]\n m$_{{DM}}$={fluxes.mDM.to('TeV')},\
          channel={channel}"
)
plt.show()
Flux [1 / (s cm2)] m$_{DM}$=10.0 TeV, channel=Z
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/astrophysics/astro_dark_matter.py:280: GammapyDeprecationWarning: The DarkMatterDecaySpectralModel class is deprecated and may be removed in a future version.
        Use DarkMatterSpectralModel instead.
  diff_flux_decay = DarkMatterDecaySpectralModel(
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/gammapy/astro/darkmatter/spectra.py:905: GammapyDeprecationWarning: "mass" was deprecated in version 2.2 and will be removed in a future version. Use argument "mDM" instead.
  super().__init__(*args, **kwargs)
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.12/site-packages/astropy/utils/decorators.py:620: GammapyDeprecationWarning: "jfactor" was deprecated in version 2.2 and will be removed in a future version. Use argument "factor" instead.
  return function(*args, **kwargs)

Gallery generated by Sphinx-Gallery