# Phase-space integral of n-body decay

Could anyone please hint me what is the best way to calculate numerically the phase-space integral $\Phi$? $$\Phi={\int\limits_{-\infty}^{+\infty}}_{3n\;times}\hspace{-13mm}\dots\hspace{-0.7mm}\dots\hspace{-0.7mm}\dots{\int\limits_{-\infty}^{+\infty}}\delta^{4}\Biggl(\begin{pmatrix}0\\0\\0\\m_0\end{pmatrix}-\sum_{i=1}^{n}\begin{pmatrix}p_{i_{x}}\\p_{i_{y}}\\p_{i_{z}}\\E(m_i,p_{i_{x}},p_{i_{y}},p_{i_{z}})\end{pmatrix}\Biggr)\prod_{i=1}^{n}\frac{dp_{i_{x}}dp_{i_{y}}dp_{i_{z}}}{16\pi^3\times E(m_i,p_{i_{x}},p_{i_{y}},p_{i_{z}})},$$

where $E(m_i,p_{i_{x}},p_{i_{y}},p_{i_{z}}) = \sqrt{m_{i}^2 + p_{i_{x}}^2 + p_{i_{y}}^2 +p_{i_{z}}^2}$. All the parameters $m_i$ are fixed and satisfy: $$m_0>\sum_{i=1}^{n}m_{i}\hspace{15mm}\rm and \hspace{15mm}m_i>0,~~~i=0,1,2,\dots,n$$

The number of parameters $n$ is fixed within single $\Phi$-integral, but I'd prefer to be able to change it easily (as well as $m_i$) without rewriting the script from scratch.

For simplicity let's put $n=3,~m_0=5,~m_1=m_2=m_3=1.$

I tried to calculate it the way I used to calculate integrals over the implicitly declared area, but this method fails in this particular case:

m0 = 5;
m1 = 1;
m2 = 1;
m3 = 1;

energy[px_, py_, pz_, m_] := Sqrt[px^2 + py^2 + pz^2 + m^2];

df1[px_, py_, pz_] := 1./(2.*energy[px, py, pz, m1] (2*Pi)^3);
df2[px_, py_, pz_] := 1./(2.*energy[px, py, pz, m2] (2*Pi)^3);
df3[px_, py_, pz_] := 1./(2.*energy[px, py, pz, m3] (2*Pi)^3);

PhaseSpaceRegion = ImplicitRegion[
(px1 + px2 + px3)^2 +
(py1 + py2 + py3)^2 +
(pz1 + pz2 + pz3)^2 +
(Sqrt[px1^2 + py1^2 + pz1^2 + m1^2] +
Sqrt[px2^2 + py2^2 + pz2^2 + m2^2] +
Sqrt[px3^2 + py3^2 + pz3^2 + m3^2] - m0)^2 == 0,
{px1, py1, pz1, px2, py2, pz2, px3, py3, pz3}];

NIntegrate[
df1[px1, py1, pz1]
df2[px2, py2, pz2]
df3[px3, py3, pz3]
{px1, py1, pz1, px2, py2, pz2, px3, py3, pz3} \[Element] PhaseSpaceRegion]

Unfortunately it fails with the following error:

DiscretizeRegion::cdim: The region given at position 1 in
DiscretizeRegion[ImplicitRegion[(px1+py1+pz1)^2+(px2+py2+pz2)^2+(px3+py3+pz3)^2+
(-5+Sqrt[Plus[<<4>>]]+Sqrt[Plus[<<4>>]]+Sqrt[Plus[<<4>>]])^2==0,
{px1,py1,pz1,px2,py2,pz2,px3,py3,pz3}]] is in dimension 9.
DiscretizeRegion only supports dimensions 1 through 3.

NIntegrate::nsr: Automatic is not a valid specification of an integration
strategy or rule.

P.S. This is quite common construction in High Energy Physics, maybe someone could suggest a package with predefined functions or anything else somehow useful in this case?

• A couple suggestions: (1) Mathematica supports DiracDelta, so you may be able to implement that directly. (2) You'll probably have more luck if you describe your implicit region as a list of separate equations $f_1 = 0$, $f_2 = 0$, etc. rather than as a single equation $f_1^2 + f_2^2 + ... = 0$. While these are mathematically equivalent, Mathematica usually finds it easier to deal with the former than the latter (it has to do with the gradient vanishing on the constraint surface.) Commented Apr 12, 2018 at 20:17
• Also, I'm pretty sure those functions should be px1 + px2 + px3 = 0 (and similarly for y and z); right now you have px1 + py1 + pz1 = 0 (and similarly for 2 and 3.) Commented Apr 12, 2018 at 20:24
• Oh, sure, you're right. Thanks! I've edited this typo in the question. I've also tried to separate equations but it still doesn't work. Commented Apr 12, 2018 at 21:14
• I'm doing this exact same thing now. Did you ever find anything? Commented Feb 25, 2021 at 19:45

There's always Monte Carlo integration. It helps a lot to do the trivial integrals delta-function integrals beforehand, so we do the momentum integrals explicitly with respect to $\vec{p}_1$, setting $\vec{p}_1 = - \vec{p}_2 - \vec{p}_3$ throughout and leaving us with a 6-D integral over $\vec{p}_2$ and $\vec{p}_3$. This reduced integral still contains a single delta function that enforces energy conservation, which we "thicken" into a "pulse" functions of width $\epsilon$ and height $1/\epsilon$ in order to get a non-zero chance of selecting a point. We also note that no component of $\vec{p}_i$ can have a magnitude greater than $\sqrt{m_0^2 - m_i^2}$; so our region of integration is bounded within $\pm (\sqrt{m_0^2 - m_i^2} + \epsilon)$ in all six dimensions.

eps = 0.1;

momentumconstraints = {px1 -> -(px2 + px3), py1 -> -(py2 + py3), pz1 -> -(pz2 + pz3)}
integrand[{px2_, py2_, pz2_, px3_, py3_, pz3_}] = df1[px1, py1, pz1] df2[px2, py2, pz2] df3[px3, py3, pz3] /. momentumconstraints
energyconstraint[{px2_, py2_, pz2_, px3_, py3_, pz3_}] := Abs[Sqrt[px1^2 + py1^2 + pz1^2 + m1^2] +
Sqrt[px2^2 + py2^2 + pz2^2 + m2^2] +
Sqrt[px3^2 + py3^2 + pz3^2 + m3^2] - m0 /. momentumconstraints] <= eps/2;

ptssel = 0;
sum = 0;
sqsum = 0;
vol = (2 m0)^6;
npts = 10^7;
bound = Sqrt[m0^2 - Min[m2,m3]^2] + eps;
For[i = 1, i <= npts, i++,
{px2, py2, pz2, px3, py3, pz3} = RandomReal[{-bound, bound}, 6];
If[energyconstraint[{px2, py2, pz2, px3, py3, pz3}],
sum = sum + integrand[{px2, py2, pz2, px3, py3, pz3}];
sqsum = sqsum + integrand[{px2, py2, pz2, px3, py3, pz3}]^2;
ptssel++;
]
]

integralestimate = (vol sum)/(eps npts);
integralvariance = 1/eps vol/Sqrt[npts] Sqrt[sqsum/npts - (sum/npts)^2]

This takes about 10 minutes to run on my (not terribly fast) desktop machine, and results in an estimate of $1.23 \pm 0.05 \times 10^{-6}$.

This is all probably horrendously unclear, and it may be partly or entirely wrong (I'm by no means an expert in Monte Carlo integration.) Please ask questions about what I've done here, and I'll either edit or delete this answer as appropriate.

• Sorry for the late response. In your code there are few errors in terms of formulas, but I get the main idea and rewrote it. Unfortunately, in adequate time limits (~10 minutes) of calculation the results I got were unstable, i.e. dispersion is of the order of the mean value. Therefore I’ve written similar code in C++ to achieve more efficient performance. And it worked well for 3-body decays. But unfortunately for 4 and 5-body decays it turned out to be again too cpu consuming calculation as the complexity grows as n^m, where m is dimension of phase-space (number of final-state particles). Commented Nov 19, 2018 at 14:39