Phase-space integral of n-body decay

Question

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?

Answer 1

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.