Here we follow the method introduced in Refs. [1]}, [2]}, [3]} to treat the phase space integrals. In this method,
all the variables appearing in the phase space integrals are Lorentz scalars.
For a process \(p \rightarrow p_1 + p_2 + p_3 + p_4\) , the corresponding phase space integral is
\(\begin{aligned}\Phi _{4} &=\frac{1}{(2\pi )^{12}}\int \frac{d^{3} \overrightarrow{p_{1}}}{2 E_{1}} \frac{d^{3} \overrightarrow{p_{2}}}{2 E_{2}} \frac{d^{3} \overrightarrow{p_{3}}}{2 E_{3}} \frac{d^{3} \overrightarrow{p_{4}}}{2 E_{4}} \delta ^{(4)}\left(p-p_{1}-p_{2}-p_{3}-p_{4}\right) \\&=\frac{\pi ^{2}}{32(2\pi )^{12} m_{0}^{2}} \int d S_{12} d S_{123} d S_{34} d S_{13} d S_{134} \frac{1}{\sqrt{-\Delta _{4}\left(p_{1}, p_{2}, p_{3}, p_{4}\right)}}\,,\end{aligned}\) 
