Spectra (`hazma.spectra`)¶

Overview¶

The hazma.spectra module contains functions for:

Computing photon,positron and neutrino spectra from the decays of individual SM particles,
Computing spectra from an \(N\)-body final state,
Boosting spectra into new frames,
Computing FSR spectra using Altarelli-Parisi splitting functions.

Below, we demonstrate how each of these actions can be easily performed in hazma.

Computing spectra (General API)¶

To compute the spectra for an arbitrary final state, use hazma.spectra.dnde_photon(), hazma.spectra.dnde_positron() or hazma.spectra.dnde_neutrino(). These functions takes the energies to compute the spectra, the center-of-mass energy and the final states and returns the spectrum.

import numpy as np
from hazma import spectra

cme = 1500.0 # 1.5 GeV
energies = np.geomspace(1e-3, 1.0, 100) * cme

# Compute the photon spectrum from a muon, a long kaon and a charged kaon
dnde = spectra.dnde_photon(energies, cme, ["mu", "kl", "k"])

plt.plot(energies, energies **2 * dnde)
plt.yscale("log")
plt.xscale("log")
plt.xlabel(r"$E_{\gamma} \ [\mathrm{MeV}]$", fontsize=16)
plt.ylabel(r"$E_{\gamma}^2\dv*{N}{E_{\gamma}} \ [\mathrm{MeV}]$", fontsize=16)
plt.tight_layout()
plt.ylim(1e-1, 2e2)
plt.xlim(np.min(energies), np.max(energies))
plt.show()

(Source code)

Under the hood, hazma computes the spectra from an \(N\)-body final state by computing the expectation values of the decay spectra from each particle w.r.t. each particles energy distribution. In the case of photon spectra, we include both the spectra from the decays of final state particles as well as the final state radiation from charged states. The final spectrum is:

\[\frac{dN}{dE} = \sum_{f}\frac{dN_{\mathrm{dec.}}}{dE} + \sum_{f}\frac{dN_{\mathrm{FSR}}}{dE}\]

where the decay and FSR components are:

\[\begin{split}\frac{dN_{\mathrm{dec}}}{dE} &= \sum_{f}\int{d\epsilon_f} P(\epsilon_f) \frac{dN_{f}}{dE}(E,\epsilon_f)\\ \frac{dN}{dE} &= \sum_{i\neq j}\frac{1}{S_iS_j}\int{ds_{ij}} P(s_{ij}) \left( \frac{dN_{i,\mathrm{FSR}}}{dE}(E,s_{ij}) + \frac{dN_{j,\mathrm{FSR}}}{dE}(E,s_{ij})\right)\end{split}\]

where \(P(\epsilon_f)\) is the probability distribution of particle \(f\) having energy \(\epsilon_f\). The sum over f in the decay expression is over all final-state particles. For the FSR expression, \(s_{ij}\) is the squared invariant-mass of particles \(i\) and \(j\). The factors \(S_{i}\) and \(S_{j}\) are symmetry factors included to avoid overcounting.

We can see what the distributions look like using the hazma.phase_space module. For three body final state consisting of a muon, long kaon and charged pion, we can compute the energy distributions and invariant mass distributions in a couple of ways. The first option is using Monte-Carlo integration using Rambo:

from hazma import phase_space
from hazma.parameters import standard_model_masses as sm_masses

cme = 1500.0
states = ["mu", "kl", "k"]
masses = tuple(sm_masses[s] for s in states)
dists = phase_space.Rambo(cme, masses).energy_distributions(n=1<<14, nbins=30)

The second option is to use the specialized ThreeBody class:

from hazma import phase_space
from hazma.parameters import standard_model_masses as sm_masses

cme = 1500.0
states = ["mu", "kl", "k"]
masses = tuple(sm_masses[s] for s in states)
dists = phase_space.ThreeBody(cme, masses).energy_distributions(nbins=30)

(Source code)

For positron and neutrino spectra, the calculation is the same except that we omit the FSR.

Individual Particle Spectra¶

hazma provides access to the individual spectra for all supported particles. All the photon, positron and neutrino spectra follow the naming convention dnde_photon_*, dnde_positron_* and dnde_neutrino_*. Each spectrum function takes in the energies of the photon/positron/neutrino and the energy of the decaying particle are returns the spectrum. As an example, let’s generate all the photon spectra for the particles available in hazma:

import numpy as np
import matplotlib.pyplot as plt
from hazma import spectra

cme = 1500.0
es = np.geomspace(1e-4, 1.0, 100) * cme

# Make a dictionary of all the available decay spectrum functions
dnde_fns = {
   "mu":    spectra.dnde_photon_muon,
   "pi":    spectra.dnde_photon_charged_pion,
   "pi0":   spectra.dnde_photon_neutral_pion,
   "k":     spectra.dnde_photon_charged_kaon,
   "kl":    spectra.dnde_photon_long_kaon,
   "ks":    spectra.dnde_photon_short_kaon,
   "eta":   spectra.dnde_photon_eta,
   "etap":  spectra.dnde_photon_eta_prime,
   "rho":   spectra.dnde_photon_charged_rho,
   "rho0":  spectra.dnde_photon_neutral_rho,
   "omega": spectra.dnde_photon_omega,
   "phi":   spectra.dnde_photon_phi,
}

# Compute the spectra. Each function takes the energy where the spectrum
# is evaluated as the 1st argument and the energy of the decaying
# particle as the 2nd
dndes = {key: fn(es, cme) for key, fn in dnde_fns.items()}

plt.figure(dpi=150)
for key, dnde in dndes.items():
   plt.plot(es, es**2 * dnde, label=key)
plt.yscale("log")
plt.xscale("log")
plt.ylim(1e-2, 1e4)
plt.xlim(np.min(es), np.max(es))
plt.ylabel(r"$E_{\gamma}^2\dv{N}{dE} \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.xlabel(r"$E_{\gamma} \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.legend()
plt.tight_layout()

(Source code)

You can also produce the same results as above using dnde_photon(). To demonstrate this, let’s repeat the above with the positron and neutrino spectra:

import numpy as np
import matplotlib.pyplot as plt
from hazma import spectra

cme = 1500.0
es = np.geomspace(1e-3, 1.0, 100) * cme

# `dnde_photon`, `dnde_positron` and `dnde_neutrino` all have an
# `availible_final_states` attribute that returns a list of strings for all
# the available particles
e_states = spectra.dnde_positron.availible_final_states
nu_states = spectra.dnde_neutrino.availible_final_states

dndes_e = {key: spectra.dnde_positron(es, cme, key) for key in e_states}
dndes_nu = {key: spectra.dnde_neutrino(es, cme, key) for key in nu_states}

plt.figure(dpi=150, figsize=(12,4))
axs = [plt.subplot(1,3,i+1) for i in range(3)]
for key, dnde in dndes_e.items():
    axs[0].plot(es, es**2 * dnde, label=key)

for key, dnde in dndes_nu.items():
    axs[1].plot(es, es**2 * dnde[0], label=key)
    axs[2].plot(es, es**2 * dnde[1], label=key)

titles = ["positron", "electron-neutrino", "muon-neutrino"]
for i, ax in enumerate(axs):
    ax.set_yscale("log")
    ax.set_xscale("log")
    ax.set_xlim(np.min(es), np.max(es))
    ax.set_ylim(1e-2, 1e3)
    ax.set_xlim(np.min(es), np.max(es))
    ax.set_title(titles[i])
    ax.set_ylabel(r"$E^2\dv{N}{E} \ [\mathrm{MeV}]$", fontdict=dict(size=16))
    ax.set_xlabel(r"$E \ [\mathrm{MeV}]$", fontdict=dict(size=16))
    ax.legend()

plt.tight_layout()

(Source code)

Including the Matrix Element¶

All three of the functions, dnde_photon, dnde_positron and dnde_neutrino allow the user to supply as squared matrix element. The matrix element is then used to correctly determine the energy distributions of the final state particles. If no matrix element is supplied, it is taken to be unity.

To demonstrate the use of a matrix element, we consider the decay of a muon into an electron and two neutrinos. We will compute the photon , the positron and the neutrino spectra. We compute the photon spectrum using only the FSR (no initial state radiation or internal bremstralung from the W.)

import numpy as np
import matplotlib.pyplot as plt
from hazma import spectra
from hazma.parameters import GF
from hazma.parameters import muon_mass as MMU

def msqrd(s, t):
   return 16.0 * GF**2 * t * (MMU**2 - t)

es = np.geomspace(1e-4, 1.0, 100) * MMU
final_states = ["e", "ve", "vm"]

dndes = {
   1: {
      # Photon
      r"$\gamma$ approx": spectra.dnde_photon(es, MMU, final_states, msqrd=msqrd),
      r"$\gamma$ analytic": spectra.dnde_photon_muon(es, MMU),
   },
   2: {
      # Positron
      r"$e^{+}$ approx": spectra.dnde_positron(es, MMU, final_states, msqrd=msqrd),
      r"$e^{+}$ analytic": spectra.dnde_positron_muon(es, MMU),
   },
   3: {
      # Electron and muon neutrino
      r"$\nu_e$ approx": spectra.dnde_neutrino(es, MMU, final_states, msqrd=msqrd, flavor="e"),
      r"$\nu_e$ analytic": spectra.dnde_neutrino_muon(es, MMU, flavor="e"),
      r"$\nu_{\mu}$ approx": spectra.dnde_neutrino(es, MMU, final_states, msqrd=msqrd, flavor="mu"),
      r"$\nu_{\mu}$ analytic": spectra.dnde_neutrino_muon(es, MMU, flavor="mu"),
   },
}


plt.figure(dpi=150, figsize=(12, 4))
for i, d in dndes.items():
   ax = plt.subplot(1, 3, i)
   lss = ["-", "--", "-", "--"]
   for j, (key, dnde) in enumerate(d.items()):
      ax.plot(es, dnde, ls=lss[j], label=key)
   ax.set_yscale("log")
   ax.set_xscale("log")
   ax.set_xlim(1e-1, 1e2)
   ax.set_ylim(1e-5, 1)
   if i == 1:
      ax.set_ylabel(r"$\dv{N}{E} \ [\mathrm{MeV}]$", fontdict=dict(size=16))
   ax.set_xlabel(r"$E \ [\mathrm{MeV}]$", fontdict=dict(size=16))
   ax.legend()

plt.tight_layout()

(Source code)

We can see that the results match the analytic results in large portions of the energy range. For the positron and neutrino spectra, the disagreements for low energies is due to the energy binning (we a linear binning procedure.) If we set nbins to a larger value, we can cover more of the energy range.

Boosting Spectra¶

It’s common that one is interested in a spectrum in a boosted frame. hazma provides facilities to boost a spectrum from one frame to another. See <table boost> for the available functions. As an example, let’s take the photon spectrum from a muon decay and boost it from the muon rest frame to a new frame. There are multiple ways to do this. The first way would be to compute the spectrum in the rest frame and boost the array. This is done using:

import matplotlib.pyplot as plt
import numpy as np
from hazma import spectra
from hazma.parameters import muon_mass as MMU

beta = 0.3
gamma = 1.0 / np.sqrt(1.0 - beta**2)
emu = gamma * MMU

es = np.geomspace(1e-4, 1.0, 100) * MMU
# Rest frame spectrum (emu = muon mass)
dnde_rf = spectra.dnde_photon_muon(es, MMU)
# Analytic boosted spectrum
dnde_boosted_analytic = spectra.dnde_photon_muon(es, emu)
# Approximate boosted spectrum
dnde_boosted = spectra.dnde_boost_array(dnde_rf, es, beta=beta)

plt.figure(dpi=150)
plt.plot(es, dnde_boosted, label="dnde_boost_array")
plt.plot(es, dnde_boosted_analytic, label="analytic", ls="--")
plt.yscale("log")
plt.xscale("log")
plt.ylabel(r"$\dv{N}{E} \ [\mathrm{MeV}^{-1}]$", fontdict=dict(size=16))
plt.xlabel(r"$E \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.legend()
plt.show()

(Source code)

This method will always have issues near the minimum energy, as we have no information about the spectrum below the minimum energy. A second way of performing the calculation is to use make_boost_function(). This method take in the spectrum in the original frame and returns a new function which is able to compute the boost integral.

import matplotlib.pyplot as plt
import numpy as np
from hazma import spectra
from hazma.parameters import muon_mass as MMU

beta = 0.3
gamma = 1.0 / np.sqrt(1.0 - beta**2)
emu = gamma * MMU

es = np.geomspace(1e-4, 1.0, 100) * emu

# Rest frame function:
dnde_rf_fn = lambda e: spectra.dnde_photon_muon(e, MMU)
# New boosted spectrum function:
dnde_boosted_fn = spectra.make_boost_function(dnde_rf_fn)

# Analytic boosted spectrum
dnde_boosted_analytic = spectra.dnde_photon_muon(es, emu)
# Approximate boosted spectrum
dnde_boosted = dnde_boosted_fn(es, beta=beta)

plt.figure(dpi=150)
plt.plot(es, dnde_boosted, label="make_boost_function")
plt.plot(es, dnde_boosted_analytic, label="analytic", ls="--")
plt.yscale("log")
plt.xscale("log")
plt.ylabel(r"$\dv{N}{E} \ [\mathrm{MeV}^{-1}]$", fontdict=dict(size=16))
plt.xlabel(r"$E \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.legend()
plt.show()

(Source code)

The last method is to use dnde_boost(), which similar to make_boost_function() but it returns the boosted spectrum rather than a function. To use it, do the following:

import numpy as np
from hazma import spectra
from hazma.parameters import muon_mass as MMU

beta = 0.3
gamma = 1.0 / np.sqrt(1.0 - beta**2)
emu = gamma * MMU

es = np.geomspace(1e-4, 1.0, 100) * emu

# Analytic boosted spectrum
dnde_boosted_analytic = spectra.dnde_photon_muon(es, emu)
# Approximate boosted spectrum
dnde_boosted = spectra.dnde_boost(
   lambda e: spectra.dnde_photon_muon(e, MMU),
   es,
   beta=beta
)

Approximate FSR Using the Altarelli-Parisi Splitting Functions¶

In the limit that the center-of-mass energy of a process is much larger than the mass of a charged state, the final state radiation is approximately equal to the it’s splitting function (with additional kinematic factors.) The splitting functions for scalar and fermionic states are:

\[\begin{split}P_{S}(x) &= \frac{2(1-x)}{x},\\ P_{F}(x) &= \frac{1 + (1-x)^2}{x}.\end{split}\]

In these expressions, \(x = 2 E / \sqrt{s}\), with \(E\) being the particle’s energy and \(\sqrt{s}\) the center-of-mass energy. Given some process \(I \to A + B + \cdots + s\) or \(I \to A + B + \cdots + f\) (with \(s\) and \(f\) being a scalar and fermionic state), the approximate FSR spectra are:

\[\begin{split}\dv{N_{S}}{E} &= \frac{Q_{S}^2\alpha_{\mathrm{EM}}}{\sqrt{s}\pi} P_{S}(x)\qty(\log(\frac{s(1-x)}{m_{S}^{2}}) - 1)\\ \dv{N_{F}}{E} &= \frac{Q_{F}^2\alpha_{\mathrm{EM}}}{\sqrt{s}\pi} P_{F}(x)\qty(\log(\frac{s(1-x)}{m_{F}^{2}}) - 1)\end{split}\]

where \(Q_{S,F}\) and \(m_{S,F}\) are the charges and masses of the scalar and fermion. To compute these spectra in hazma, we provide the functions dnde_photon_ap_scalar() and dnde_photon_ap_fermion(). Below we demonstrate how to use these:

import matplotlib.pyplot as plt
import numpy as np
from hazma import spectra
from hazma.parameters import muon_mass as MMU
from hazma.parameters import charged_pion_mass as MPI

cme = 500.0 # 500 MeV
es = np.geomspace(1e-4, 1.0, 100) * cme

# These functions take s = cme^2 as second argument and mass of the state as
# the third.
dnde_fsr_mu = spectra.dnde_photon_ap_fermion(es, cme**2, MMU)
dnde_fsr_pi = spectra.dnde_photon_ap_scalar(es, cme**2, MPI)

plt.figure(dpi=150)
plt.plot(es, dnde_fsr_mu, label=r"$\mu$")
plt.plot(es, dnde_fsr_pi, label=r"$\pi$", ls="--")
plt.yscale("log")
plt.xscale("log")
plt.ylabel(r"$\dv{N}{E} \ [\mathrm{MeV}^{-1}]$", fontdict=dict(size=16))
plt.xlabel(r"$E \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.legend()
plt.show()

(Source code)

API Reference¶

Altarelli-Parisi¶

Functions for compute FSR spectra using the Altarelli-Parisi splitting functions:

Function	Description
`dnde_photon_ap_fermion()`	Approximate FSR from fermion.
`dnde_photon_ap_scalar()`	Approximate FSR from scalar.

Boost¶

Function	Description
`boost_delta_function()`	Boost a delta-function.
`double_boost_delta_function()`	Perform two boosts of a delta-function.
`dnde_boost_array()`	Boost spectrum specified as an array.
`dnde_boost()`	Boost spectrum specified as function.
`make_boost_function()`	Construct a boost function.

Photon Spectra¶

hazma includes decay spectra from several unstable particles:

Function	Description
`dnde_photon()`	Photon spectrum from decay/FSR of n-body final states.
`dnde_photon_muon()`	Photon spectrum from decay of \(\mu^{\pm}\).
`dnde_photon_neutral_pion()`	Photon spectrum from decay of \(\pi^{0}\).
`dnde_photon_charged_pion()`	Photon spectrum from decay of \(\pi^{\pm}\).
`dnde_photon_charged_kaon()`	Photon spectrum from decay of \(K^{\pm}\).
`dnde_photon_long_kaon()`	Photon spectrum from decay of \(K_{L}\).
`dnde_photon_short_kaon()`	Photon spectrum from decay of \(K_{S}\).
`dnde_photon_eta()`	Photon spectrum from decay of \(\eta\).
`dnde_photon_eta_prime()`	Photon spectrum from decay of \(\eta'\).
`dnde_photon_charged_rho()`	Photon spectrum from decay of \(\rho^{\pm}\).
`dnde_photon_neutral_rho()`	Photon spectrum from decay of \(\rho^{0}\).
`dnde_photon_omega()`	Photon spectrum from decay of \(\omega\).
`dnde_photon_phi()`	Photon spectrum from decay of \(\phi\).

Positron Spectra¶

Function	Description
`dnde_positron()`	Positron spectrum from decay of n-body final states.
`dnde_positron_muon()`	Positron spectrum from decay of \(\mu^{\pm}\).
`dnde_positron_charged_pion()`	Positron spectrum from decay of \(\pi^{\pm}\).
`dnde_positron_charged_kaon()`	Positron spectrum from decay of \(K^{\pm}\).
`dnde_positron_long_kaon()`	Positron spectrum from decay of \(K_{L}\).
`dnde_positron_short_kaon()`	Positron spectrum from decay of \(K_{S}\).
`dnde_positron_eta()`	Positron spectrum from decay of \(\eta\).
`dnde_positron_eta_prime()`	Positron spectrum from decay of \(\eta'\).
`dnde_positron_charged_rho()`	Positron spectrum from decay of \(\rho^{\pm}\).
`dnde_positron_neutral_rho()`	Positron spectrum from decay of \(\rho^{0}\).
`dnde_positron_omega()`	Positron spectrum from decay of \(\omega\).
`dnde_positron_phi()`	Positron spectrum from decay of \(\phi\).

Neutrino Spectra¶

Function	Description
`dnde_neutrino()`	Neutrino spectrum from decay of n-body final states.
`dnde_neutrino_muon()`	Neutrino spectrum from decay of \(\mu^{\pm}\).
`dnde_neutrino_charged_pion()`	Neutrino spectrum from decay of \(\pi^{\pm}\).
`dnde_neutrino_charged_kaon()`	Neutrino spectrum from decay of \(K^{\pm}\).
`dnde_neutrino_long_kaon()`	Neutrino spectrum from decay of \(K_{L}\).
`dnde_neutrino_short_kaon()`	Neutrino spectrum from decay of \(K_{S}\).
`dnde_neutrino_eta()`	Neutrino spectrum from decay of \(\eta\).
`dnde_neutrino_eta_prime()`	Neutrino spectrum from decay of \(\eta'\).
`dnde_neutrino_charged_rho()`	Neutrino spectrum from decay of \(\rho^{\pm}\).
`dnde_neutrino_omega()`	Neutrino spectrum from decay of \(\omega\).
`dnde_neutrino_phi()`	Neutrino spectrum from decay of \(\phi\).

Spectra (hazma.spectra)¶

Overview¶

Computing spectra (General API)¶

Individual Particle Spectra¶

Including the Matrix Element¶

Boosting Spectra¶

Approximate FSR Using the Altarelli-Parisi Splitting Functions¶

API Reference¶

Altarelli-Parisi¶

Boost¶

Photon Spectra¶

Positron Spectra¶

Neutrino Spectra¶

Spectra (`hazma.spectra`)¶