How to reduce an integral in phase space to a one-dimensional form?

[RicardoMP]

I've been trying for a very long time to show that the following integral:
$$ I_D=2{\displaystyle \int} \, {\displaystyle \prod_{i=1}^3} d \Pi_i \, (2\pi
)^4\delta^4(p_H-p_L-p_R) |{\cal M}({e_L}^c e_R \leftrightarrow h^*)|^2
f_{L}^0f_{R}^0(1+f_{H}^0). $$
can be reduced to one dimension:
$$
I_D = \frac{m_H T^3h_e^2\gamma^2}{\pi^3}\int_1^\infty d u \frac{e^{um_H/T}}
{(e^{um_H/T}-1)^2}
\! \ln \left( \frac {\textstyle \cosh (\alpha_{e_L} u + \gamma \sqrt{u^2-1})}
{\textstyle \cosh (\alpha_{e_L} u - \gamma \sqrt{u^2-1})}
% \right. \nonumber \\ & \t \left.
\frac {\textstyle \cosh (\alpha_{e_R} u + \gamma \sqrt{u^2-1})}
{\textstyle \cosh (\alpha_{e_R} u - \gamma \sqrt{u^2-1})} \right).$$
where:

$$\alpha_{e_L}\equiv (m_H^2+m_{e_L}^2-m_{e_R}^2)/4m_HT$$
$$\alpha_{e_R}\equiv (m_H^2+m_{e_R}^2-m_{e_L}^2)/4m_HT$$
$$\gamma\equiv \lambda^{\frac{1}{2}}(m_H^2, m_{e_L}^2, m_{e_R}^2)/4m_HT$$
$$\lambda(x,y,z)\equiv (x-y-z)^2-4yz$$
and $f^0_i=(e^{\beta E_i}\pm 1)^{-1}$, whether the particle species is a fermion(+) or a boson(-).
My attempt started by determining $|{\cal M}({e_L}^c e_R \leftrightarrow h^*)|^2$ which, for the indicated process, gave me the following:
$$I_D=2{\displaystyle \int} \, \frac{d^3p_H}{(2\pi)^32E_H}{\displaystyle \int} \, \frac{d^3p_L}{(2\pi)^32E_L}{\displaystyle \int} \, \frac{d^3p_R}{(2\pi)^32E_R}\, (2\pi
)^4\delta(E_H-E_L-E_R)\delta^3(\vec{p_H}-\vec{p_L}-\vec{p_R}) [2h_e^2(|\vec{p_L}||\vec{p_R}|-\vec{p_L}.\vec{p_R})]\frac{e^{\beta E_H}}{(e^{\beta E_H}-1)(e^{\beta E_L}+1)(e^{\beta E_R}+1)}.$$

I've proceeded by separating the integral and focus on the manifestly Lorentz invariant term $\vec{p_L}.\vec{p_R}$ and on the center of mass frame where the 4-momenta $p_H =(E_H,\vec{0})$:
$$I_D=\frac{h_e^2}{2(2\pi)^5}{\displaystyle \int} \frac{d^3p_H}{E_H}{\displaystyle \int} \frac{d^3p_L}{E_L}{\displaystyle \int} \, \frac{d^3p_R}{E_R}\delta(m_H-E_L-E_R)\delta^3(\vec{p_L}+\vec{p_R})\vec{p_L}.\vec{p_R}\frac{e^{\beta E_H}}{(e^{\beta E_H}-1)(e^{\beta E_L}+1)(e^{\beta E_R}+1)}.$$
By using some of the delta function properties and the conservation of momentum I easily arrive to:
$$\delta(m_H-E_L-E_R) = \delta(m_H-\sqrt{m_L^2 + |\vec{p_L}|}-\sqrt{m_R^2 + |\vec{p_R}|}) \equiv \delta(f(|\vec{p_R}|)) = \left|\frac{df}{d|\vec{p_R}|}\right|^{-1}_{p*}\delta(|\vec{p_r}|-p*) .$$
$$p* = \frac{1}{2m_H}\lambda^{\frac{1}{2}}(m_H^2, m_{e_L}^2, m_{e_R}^2)$$
$$\left|\frac{df}{d|\vec{p_R}|}\right|^{-1}_{p*}=\frac{1}{p*}\frac{E_RE_L}{E_R+E_L}$$
The delta function on the momenta enabled me to remove one of the integrals (in this case, the integral over $\int d^3p_L$. The new delta function $\delta(|\vec{p_r}|-p*)$ will allow me to take out the integral over $\int dp_R$ that comes from $d^3p_R = |\vec{p_R}|^2dp_Rd\Omega$ and so:

$$I_D=\frac{h_e^2}{2(2\pi)^5}{\displaystyle \int} \frac{d^3p_H}{E_H}{\displaystyle \int}\frac{p*^3}{E_R+E_L}\frac{e^{\beta E_H}}{(e^{\beta E_H}-1)(e^{\beta E_L}+1)(e^{\beta E_R}+1)}d\Omega.$$

From here on out, I have no idea on what to do in order to arrive to the one dimensional form, as I've no clue on how to fit the natural logarithm of those hyperbolic cosines in the integrand and neither which substitution I should do.
Any help would be welcome!

 

[RGevo]

Is this a process of the form (1->2)? Cannot one immediately consider the 2body final state phase space of eL and eR evaluated in h CoM frame? Then apply the same procedure as you did and just deal with the dOmega integral.

 

[RicardoMP]

Is this a process of the form (1->2)? Cannot one immediately consider the 2body final state phase space of eL and eR evaluated in h CoM frame? Then apply the same procedure as you did and just deal with the dOmega integral.

It is indeed a decay process of the form $1\rightarrow 2$, in particular, a Higgs decaying to an electron pair. You're suggesting using the same procedure of going into the Higgs CoM frame on the last equation? Nonetheless, I have no idea on how to move on from there.

 

[RGevo]

I am instead suggesting to perform the phase-space sum with i = eL, eR meaning that you will have an integrals over the momenta of those out-going particles only (and not also the Higgs).

This will mean that the manipulations you did to deal with the delta functions before can be applied in the same way. You will then reach the same point with an integral over the solid angle (but not over the Higgs).

Because there is no $Phi$ dependence, d$Omega$ = d$Phi$ d cos $Theta$, the $Phi$ integral can be performed trivially (just 2$Pi$) leaving the d cos $Theta$ integral. I suspect the format of the integral you are trying to find is then obtained by performing a change of variables from cos $Theta$.