phase_space.rst

hazma.phase_space

Phase Space (hazma.phase_space)

• Rambo <rambo section>
  • Integrating over phase-space <rambo integrate>
  • Generating phase-space points <rambo generator>
  • Computing decay widths and cross sections <rambo widths and cs>
  • Energy and Invariant Mass Distributions <rambo dists>
• Three Body Phase Space <tbps section>
  • Integrating over phase-space <tbps integrate>
  • Energy and Invariant Mass Distributions <tbps dists>

Overview

The phase_space module contains the :pyhazma.phase_space.Rambo class for working with N-body phase-space and the specialized :pyhazma.phase_space.ThreeBody class for three-body phase-space. :pyRambo contains methods for generating phase-space points, integrating over N-body phase-space, computing cross-sections and computing decay widths.

Rambo

The :py~hazma.phase_space.Rambo class is instantiated by supplying the center-of-mass energy, the masses of the final-state particles and optionally a function to compute the squared matrix element. The function to compute the squared matrix element must be a vectorized unary function that return the squared matrix elements given four-momenta. The internal momenta generated by :pyRambo have a shape of (4, nfsp, n) where nfsp in the number of final-state particles and n is the number of points generated simultaneously. The leading dimension holds the energy, x, y and z-components of the four-momenta. If no function to compute the squared matrix element is supplied, it will be taken to be 1.

In the following snippet, we create a :pyRambo object for a process with 3 final-state particles, and a squared matrix element equal to p₁ ⋅ p₂. Note the utils module contains a function :pyhazma.utils.ldot to compute the Lorentzian scalar product between two four-vectors.

from hazma import phase_space
from hazma import utils

cme = 10.0
masses = [1.0, 2.0, 3.0]

def msqrd(momenta):
   p1 = momenta[:, 0] # Pick out the 4-momenta of particle 1
   p2 = momenta[:, 1] # Pick out the 4-momenta of particle 2
   return utils.ldot(p1, p2) # Compute p1.p2

phase_space = phase_space.Rambo(cme, masses, msqrd)

Integrating over phase-space

The :pyRambo.integrate method computes the following integral:

$$\Pi_{\mathrm{LIPS}} = \int \qty(\prod_{i=1}^{N}\dfrac{d^{3}\vec{p}_{i}}{(2\pi)^3}\dfrac{1}{2E_{i}}) (2\pi)^4 \delta^{4}\qty(P-\sum_{i=1}^{N}p_{i}) {\left|\mathcal{M}\right|}^2$$

This function takes in an integer specifying the number of points to use for the Monte-Carlo integration and computes the integral using the RAMBO algorithm. As a simple example, let's compute the N-body integral for massless particles. The analytical result is:

$$\Pi^{(n)}_{\mathrm{LIPS}} = \frac{1}{\Gamma(n)\Gamma(n-1)} {(2\pi)}^{4-3n} {\qty(\frac{\pi}{2})}^{n-1} {E_{\mathrm{CM}}}^{2n-4}$$

We can verify this using the :pyRambo.integrate method:

import math
import numpy as np
from hazma import phase_space

def analytic(n):
   fact = math.factorial(n - 2) * math.factorial(n - 1)
   return (0.5 * math.pi) ** (n - 1) * (2 * np.pi) ** (4 - 3 * n) / fact
analytic_integrals = [analytic(n) for n in range(2, 10)]

integrals = []
for n in range(2, 10):
   rambo = phase_space.Rambo(1.0, [0.0] * n)
   integral, error = rambo.integrate(n=10)
   integrals.append(integral)
np.max([abs(i1 - i2) / i2 for i1, i2 in zip(integrals, analytic_integrals)])
# Possible output: 2.2608966769988504e-15

Generating phase-space points

Sometimes it is useful to have access the momenta and weights of N-body phase-space. The are two methods for generating momenta and weights: :pyRambo.generate and func:Rambo.generator. The func:Rambo.generate :py:method will return a specified number of momenta and weights. The :pyRambo.generator method returns a python generator, allowing for iterating over batches of momenta and weights.

As described above, the momenta will have a shape of (4, nfsp, n) where nfsp is the number of final-state particles and n is the number of requested points. The weights will have a shape of (n,). To see this explicitly, consider the following. Here we generate 10 phase-space points:

from hazma import phase_space

rambo = phase_space.Rambo(cme=10.0, masses=[1.0, 2.0, 3.0])
momenta, weights = rambo.generate(n=10)
print(momenta.shape)
print(weights.shape)
# (4, 3, 10)
# (10,)

In some cases, one may not want to generate all points at once (since it is costly memory-wise or maybe one wants to monitor convergence.) We supply the :pyRambo.generator method for generating batches of phase-space points. For example, suppose we want to integrate over phase-space ourselves using 1M points and batches of 50,000 points at a time. To do this, we can use:

import numpy as np
from hazma import phase_space

rambo = phase_space.Rambo(cme=10.0, masses=[1.0, 2.0, 3.0])
n = int(1e6)
batch_size = 50_000
integrals = []
for momenta, weights in rambo.generator(n, batch_size, seed=1234):
   integrals.append(np.nanmean(weights))
np.average(integrals)
# Output: 0.0036118278252665406

Note

The above code is equivalent to using :pyhazma.phase_space.Rambo.generate with n=int(1e6) and batch_size=50_000. All methods accept a batch_size argument used can to split the computation into chunks in cases where the user has limited memory.

Computing decay widths and cross sections

The most common use of the :pyRambo is computing cross sections or decay widths. The functions :pyhazma.phase_space.Rambo.cross_section and :pyhazma.phase_space.Rambo.decay_width can be used for these purposes. These methods just wrap :pyhazma.phase_space.Rambo.integrate and append the appropriate prefactors. As an example, let's compute the muon decay width for the process μ → eν_eν_μ. The squared matrix element (ignoring the electron mass) is:

|ℳ|² = 16G_F²t(m_μ² − t)

where G_F is the Fermi constant and t = (p_e + p_νμ)². The analytic result is

$$\Gamma = \frac{G_{F}^{2}m_{\mu}^{5}}{192\pi^3} \sim 3 \times 10^{-19}$$

To compute this, we use the following:

from hazma import phase_space
from hazma import utils
from hazma.parameters import GF
from hazma.parameters import muon_mass as MMU

def msqrd(momenta):
   p1 = momenta[:, 0]
   p3 = momenta[:, 2]
   t = utils.lnorm_sqr(p1 + p3)
   return 16.0 * GF**2 * t * (MMU**2 - t)


rambo = phase_space.Rambo(MMU, [0.0, 0.0, 0.0], msqrd=msqrd)
width, error = rambo.decay_width(n=50_000, seed=1234)

analytic = GF**2 * MMU**5 / (192 * np.pi**3)
actual_error = abs(width - analytic)
print(f"width = {width:.2e} +- {error:.2e}")
print(f"actual error = {actual_error:.2e} = {actual_error / analytic * 100:.2f} %")
# Output:
#  width = 3.02e-19 +- 5.99e-22
#  actual error = 9.50e-22 = 0.32 %

Note

The utils module contains a couple of functions useful for dealing with four-vectors, namely, :pyhazma.utils.ldot for computing scalar products and :pyhazma.utils.lnorm_sqr for computing the squared norm.

Energy and Invariant Mass Distributions

The :pyRambo class can also compute energy distributions as well as the invariant mass distributions of pairs of final state particles. To compute energies distributions, use :pyRambo.energy_distributions:

from hazma import phase_space from hazma.parameters import standard_model_masses as sm_masses import matplotlib.pyplot as plt

states = ["pi", "e", "mu", "k"] masses = [sm_masses[s] for s in states] cme = 3 * sum(masses) rambo = phase_space.Rambo(cme, masses)

energy_dists = rambo.energy_distributions(n=1<<16, nbins=25)

plt.figure(dpi=150) labels=[r"$pi$", r"$e$", r"$mu$", r"$K$"] for i, dist in enumerate(energy_dists): plt.plot(dist.bin_centers, dist.probabilities, label=labels[i])

plt.ylabel(r"$P(epsilon) [mathrm{MeV}^{-1}]$", fontdict=dict(size=16)) plt.xlabel(r"$epsilon [mathrm{MeV}]$", fontdict=dict(size=16)) plt.tight_layout() plt.legend()

We note the 'choppy' behavior of the curves. Since we are performing naive Monte-Carlo integration, we need a large number of points to properly sample phase-space. The invariant mass distributions are generated in a similar fashion. Note that the return value of the invariant mass distributions is a dictionary with keys given by a pair of integers that specify the pair of particles the distribution corresponds to. See :pyhazma.phase_space.PhaseSpaceDistribution1D for more information on the distribution objects.

Three Body Phase Space (:pyhazma.phase_space.ThreeBody)

Attention

This class is meant for cases where the squared matrix element only depends on the total momentum and the momenta of the final state particles. If this isn't the case, then this class is not suitable.

As an example of a case where this class does apply, recall that the squared matrix element (summed over spins) for the decay of an unstable particle into a three-body final state only depends on the momenta of the final-state particles.

For three-body phase space where the squared matrix element only depends on the final-state momenta, the phase space integral simplifies. The result is:

$$\Phi_3 = \frac{1}{16(2\pi)^2Q^2} \int_{s_{-}}^{s_{+}}\dd{s} \int_{t_{-}(s)}^{t_{+}(s)}\dd{t} {\left|\mathcal{M}\right|}^2$$

where Q is the center-of-mass energy, s = (p₂ + p₃)² and t = (p₁ + p₃)². The limits of integration are:

$$\begin{aligned} s_{\mathrm{min}} &= (m_2+m_3)^2\\\ s_{\mathrm{max}} &= (Q-m_1)^2\\\ t_{\pm} &= \frac{-s^2 +(Q^2+m_1^2+m_2^2+m_3^2) s -\qty(M^2-m_1^2)(m_2-m_3)^2 \pm p_1p_2}{2s}\\\ p_{1} &= \lambda^{1/2}(s,m_1^2,Q^2)\\\ p_{2} &= \lambda^{1/2}(s,m_2^2,m_3^2) \end{aligned}$$

Here, λ(a, b, c) = a² + b² + c² − 2ab − 2ac − 2bc.

Integrating over phase-space

The interface for :pyThreeBody is similar to :pyRambo. To integrate over phase space, use :pyThreeBody.integrate. The is an important difference between :pyThreeBody and :pyRambo: the matrix element must be a binary function taking the variables s = (p₂ + p₃)² and t = (p₁ + p₃)².

As an example, let's consider the muon decay again. To integrate over phase space, we use:

import numpy as np
from hazma import phase_space
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)

integral, error = phase_space.ThreeBody(MMU, [0.0, 0.0, 0.0], msqrd=msqrd).integrate()
width = integral / (2.0 * MMU)
error = error / (2.0 * MMU)

analytic = GF**2 * MMU**5 / (192 * np.pi**3)
actual_error = abs(width - analytic)
print(f"analytic = {analytic:.8e}")
print(f"width = {width:.8e} +- {error:.2e}")
print(f"actual error = {actual_error:.2e} = {actual_error / analytic * 100:.2e} %")
# Output:
# analytic = 3.00917842e-16
# width = 3.00917842e-16 +- 6.68e-30
# actual error = 4.93e-32 = 1.64e-14 %

Energy and Invariant Mass Distributions

The interface for computing energy and invariant-mass distributions using :pyThreeBody is nearly identical to :pyRambo. You can simply construct your :pyThreeBody instance and call either :pyThreeBody.energy_distributions or :pyThreeBody.invariant_mass_distributions. The main difference is the signature of the msqrd parameter. As before, it must be a binary function accepting s = (p₂ + p₃)² and t = (p₁ + p₃)².

The hazma.form_factors module uses ThreeBody for the 3-body form factors. As an example of how to generate distributions, let's use the hazma.form_factors.vector.VectorFormFactorPiKK0 to look at the energy distributions of the final state mesons.

import numpy as np import matplotlib.pyplot as plt from hazma import phase_space import hazma.form_factors.vector as ffv

def msqrd(s, t, q, form_factor: ffv.VectorFormFactorPiKK0, gvuu, gvdd, gvss):

ff = form_factor.form_factor(q, s, t, gvuu=gvuu, gvdd=gvdd, gvss=gvss) lor = form_factor.squared_lorentz_structure(q, s, t) return np.abs(ff) ** 2 * lor

# Specify parameters gvuu, gvdd, gvss = 2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0 pk0_ff = ffv.VectorFormFactorPiKK0() q = 1.5e3 # 1.5 GeV # Construct ThreeBody object tb = phase_space.ThreeBody( q, pk0_ff.fsp_masses, msqrd=lambda s, t: msqrd(s, t, q, pk0_ff, gvuu, gvdd, gvss), )

# Construct distributions energy_dists = tb.energy_distributions(nbins=100) inv_mass_dists = tb.invariant_mass_distributions(nbins=100)

# Make plot plt.figure(dpi=150, figsize=(9, 3)) ax1 = plt.subplot(1, 2, 1) ax2 = plt.subplot(1, 2, 2)

labels = [r"$K^0$", r"$K$", r"$pi$"] lss = ["-", "--", "-."] for i, dist in enumerate(energy_dists): ax1.plot(dist.bin_centers, dist.probabilities, label=labels[i], ls=lss[i]) ax1.set_ylabel(r"$P(epsilon)$", fontsize=16) ax1.set_xlabel(r"$epsilon$", fontsize=16)

labels = {

(0, 1): r"$(K^0,K)$", (0, 2): r"$(K^0,pi)$", (1, 2): r"$(K,pi)$",

} lss = {(0, 1): "-", (0, 2): "--", (1, 2): "-."} for key, dist in inv_mass_dists.items(): ax2.plot(dist.bin_centers, dist.probabilities, label=labels[key], ls=lss[key]) ax2.set_ylabel(r"$P(s{ij})$", fontsize=16) ax2.set_xlabel(r"$s{ij}$", fontsize=16)

ax1.legend() ax2.legend() plt.tight_layout() plt.show()

API Reference

hazma.phase_space::Rambo

hazma.phase_space::ThreeBody

hazma.phase_space::PhaseSpaceDistribution1D