Phase Space¶

Overview¶

The rambo module contains the hazma.rambo.PhaseSpace class for working with N-body phase-space. In addition, it contains older functions that use Cython and multiproccessing which tend to be slower. hazma.rambo.PhaseSpace contains methods for generating phase-space points, integrating over N-body phase-space, computing cross-sections and computing decay widths.

The hazma.rambo.PhaseSpace 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 uniary function that return the squared matrix elements given four-momenta. The internal momenta generated by hazma.rambo.PhaseSpace have the 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 hazma.rambo.PhaseSpace object for a process with 3 final-state particles, and a squared matrix element equal to \(p_{1}\cdot p_{2}\). Note the Utilities (hazma.utils) module contains a function hazma.utils.ldot() to compute the Lorentzian scalar product between two four-vectors.

from hazma import rambo
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 = PhaseSpace(cme, masses, msqrd)

Integrating over phase-space¶

The hazma.rambo.PhaseSpace.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}) \qty|\mathcal{M}|^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 hazma.rambo.PhaseSpace.integrate() method:

Integrating N-body phase space with massless particles¶

import math
import numpy as np
from hazma import rambo

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):
   phase_space = rambo.PhaseSpace(1.0, [0.0] * n)
   integral, error = phase_space.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 the momenta and weights of N-body phase-space. The are two methods for generating momenta and weights: hazma.rambo.PhaseSpace.generate() and hazma.rambo.PhaseSpace.generator(). The hazma.rambo.PhaseSpace.generate() method will return a specified number of momenta and weights. The hazma.rambo.PhaseSpace.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 explicity, consider the following. Here we generate 10 phase-space points:

from hazma import rambo

phase_space = rambo.PhaseSpace(cme=10.0, masses=[1.0, 2.0, 3.0])
momenta, weights = phase_space.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 hazma.rambo.PhaseSpace.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 rambo

phase_space = rambo.PhaseSpace(cme=10.0, masses=[1.0, 2.0, 3.0])
n = int(1e6)
batch_size = 50_000
integrals = []
for momenta, weights in phase_space.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 hazma.rambo.PhaseSpace.generate() with n=int(1e6) and batch_size=50_000. All methods accept a batch_size argument that can be used to split up the computation if memory is limited.

Computing decay widths and cross sections¶

The most common use of the hazma.rambo.PhaseSpace() is computing cross sections or decay widths. The functions hazma.rambo.PhaseSpace.cross_section() and hazma.rambo.PhaseSpace.decay_width() can be used for these purposes. These methods just wrap hazma.rambo.PhaseSpace.integrate() and append the appropriate prefactors. As an example, let’s compute the muon decay width for the process \(\mu\to e\nu_{e}\nu_{\mu}\). The squared matrix element (ignoring the electron mass) is given by:

\[\qty|\mathcal{M}|^2 = 16 G_{F}^{2} t (m_{\mu}^2 - t)\]

where \(G_{F}\) is the Fermi constant and \(t=(p_{e}+p_{\nu_{\mu}})^2\). 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 rambo
from hazma import utils

GF = 1.166379e-05
MMU = 1.056584e-01


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


phase_space = rambo.PhaseSpace(cme=mmu, masses=[0.0, 0.0, 0.0], msqrd=msqrd)
width, error = phase_space.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 Utilities (hazma.utils) module contains a couple of functions useful for dealing with four-vectors, namely, hazma.utils.ldot() for computing scalar products and hazma.utils.lnorm_sqr() for computing the squared norm.