Gamma ray spectra¶
Overview¶
This page discusses how to compute the gamma-ray spectrum for a particular
particle physics process. Since computing gamma-ray spectra is a
model-dependent process, we include in hazma
tools for computing gamma-ray
spectra from both FSR and the decay of final state-particles.
The gamma_ray
module contains two
functions called gamma_ray_decay
and gamma_ray_fsr
. The
gamma_ray_decay
function accepts a list of the final-state particles,
the center-of-mass energy, the gamma-ray energies to compute the spectrum
at and optionally the matrix element. Currently, the final-state particles
can be \(\pi^{0}\), \(\pi^{\pm}\), \(\mu^{\pm}\),
\(e^{\pm}\), \(K^{\pm}\), \(K_{L}\) and \(K_{S}\) where
\(K\) stands for kaon.
We caution that when including many final-state mesons, one needs to take
care to supply the properly unitarized matrix element. The gamma_ray_decay
functions works by first computing the energies distributions of all the
final-state particles and convolving the energy distributions with the decay
spectra of the final-state particles. The gamma_ray_decay
function can be
used as follows:
>>> from hazma.gamma_ray_decay import gamma_ray_decay
>>> import numpy as np
# specify the final-state particles
>>> particles = np.array(['muon', 'charged_kaon', 'long_kaon'])
# choose the center of mass energy
>>> cme = 5000.
# choose list of the gamma-ray energies to compute spectra at
>>> eng_gams = np.logspace(0., np.log10(cme), num=200, dtype=np.float64)
# compute the gamma-ray spectra assuming a constant matrix element
>>> spec = gamma_ray_decay(particles, cme, eng_gams)
The gamma_ray_fsr
function computes the gamma-ray spectrum from
\(X\to Y\gamma\), i.e.:
where \(X\) and \(Y\) are any particles excluding the photon. This function takes in as input a list of the initial state particle masses (either 1 or 2 particles), the final state particle masses, the center-of-mass energy, a function for the tree-level matrix element (for \(X\to Y\)) and a function for the radiative matrix element (\(X\to Y\gamma\)). The functions for the matrix elements must take is a single argument which is a list of the four-momenta for the final state particles. As an example, we consider the process of two dark-matter particles annihilating into charged pions, \(\bar{\chi}\chi\to \pi^{+}\pi^{-}(\gamma)\) using the model from Advanced Usage. In Advanced Usage, we gave the analytic expressions for the gamma-ray spectra. The tree-level and radiative matrix elements for this process are:
where \(Q\) is the center-of-mass energy, \(e\) is the electromagnetic coupling, \(\mu_{\pi,\chi} = m_{\pi,\chi}/Q\) and
with \(p_{\pi,1,2}\) are the four-momenta of the two final-state pions and \(k\) is the four-momenta of the final-state photon. Below, we create a class to implement functions for the tree and radiative matrix elements. Note that these functions take in an array of four-momenta.
from hazma.field_theory_helper_functions.common_functions import \
minkowski_dot as MDot
class Msqrd(object):
def __init__(self, mx, c1, lam):
self.mx = mx # DM mass
self.c1 = c1 # effective coupling of DM to charged pions
self.lam = lam # cut off scale for effective theory
def tree(self, momenta):
ppi1 = momenta[0] # first charged pion four-momentum
ppi2 = momenta[1] # second charged pion four-momentum
Q = ppi1[0] + ppi2[0] # center-of-mass energy
return -((self.c1**2 * (4 * self.mx**2-Q**2)) / (2 * self.lam**2))
def radiative(self, momenta):
ppi1 = momenta[0] # first charged pion four-momentum
ppi2 = momenta[1] # second charged pion four-momentum
k = momenta[2] # photon four-momentum
Q = ppi1[0] + ppi2[0] + k[0] # center-of-mass energy
mux = self.mx / Q
mupi = mpi / Q
s = MDot(ppi1 + ppi2, ppi1 + ppi2)
t = MDot(ppi1 + k, ppi1 + k)
u = MDot(ppi2 + k, ppi2 + k)
return ((2*self.c1**2*(-1 + 4*mux**2)*Q**2*qe**2 *
(s*(-(mupi**2*Q**2) + t)*(mupi**2*Q**2 - u) +
(-2*mupi**3*Q**3 + mupi*Q*(t + u))**2)) /
(self.lam**2*(-(mupi**2*Q**2) + t)**2*
(-(mupi**2*Q**2) + u)**2))
Next, we can compute the gamma-ray spectrum for \(\bar{\chi}\chi\to\pi^{+}\pi^{-}\gamma\) using:
>>> from hazma.gamma_ray import gamma_ray_fsr
# specify the parameters of the model
>>> params = {'mx': 200.0, 'c1':1.0, 'lam':1e4}
# create instance of our Msqrd class
>>> msqrds = Msqrd(**params)
# specify the initial and final state masses
>>> isp_masses = np.array([msqrds.mx, msqrds.mx])
>>> fsp_masses = np.array([mpi, mpi, 0.0])
# choose the center-of-mass energy
>>> cme = 4.0 * msqrds.mx
# compute the gamma-ray spectrum
>>> spec = gamma_ray_fsr(isp_masses, fsp_masses, cme, msqrds.tree,
msqrds.radiative, num_ps_pts=500000, num_bins=50)
# plot the spectrum
>>> import matplotlib.pyplot as plt
>>> plt.figure(dpi=100)
>>> plt.plot(spec[0], spec[1])
>>> plt.yscale('log')
>>> plt.xscale('log')
>>> plt.ylabel(r'$dN/dE_{\gamma} \ (\mathrm{MeV}^{-1})$', fontsize=16)
>>> plt.xlabel(r'$E_{\gamma} \ (\mathrm{MeV})$', fontsize=16)