In the article Cosmic abundances of stable particles: Improved analysis, P. Gondolo and G. Gelmini, Nucl. Phys. B 360 (1991), p. 145-179, they convert $\rm{d}^3p_1\rm{d}^3p_2=2\pi^2p_1p_2\rm{d}E_1\rm{d}E_2\rm{d}\cos{\theta}$ (eq.3.2) into $\rm{d}^3p_1\rm{d}^3p_2=2\pi E_1E_2\rm{d}E_+\rm{d}E_-\rm{d}s$ (eq. 3.4) with the following change of variables:
$$E_+=E_1+E_2, \quad E_-=E_1-E_2, \quad s=2m^2+2E_1E_2-2p_1p_2\cos{\theta}, \quad \rm{(eq. 3.3)}$$
When I try to derive the second expresion of $\rm{d}^3p_1\rm{d}^3p_2$ using these new variables $E_+,\ E_-, s$, I always get second order diferentials that in theory should vanish. I don't know how one could get $\rm{d}^3p_1\rm{d}^3p_2=2\pi E_1E_2\rm{d}E_+\rm{d}E_-\rm{d}s$.
They also define new limits of integration due to this change of variables: $\{ E_1>m,\ E_2>m, \ |\cos{\theta}|\leq 1 \}$ changes to $\{ s\geq 4m^2,\ E_+\geq \sqrt{s}, \ |E_-|\leq\sqrt{1-4m^2/s}\sqrt{E_+^2-s} \}$ (eq. 3.5). I get the first two, but I don't know how to compute $E_-$. The limit cases give some idea on how the expresion should be, but this is not enough to know the exact expresion, something like: $$\rm{If} \quad E_+^2=s \implies E_-=0$$ $$\rm{If} \quad s=4m^2 \implies E_-=0$$ $$\rm{Then (?):} \quad |E_-| \leq \sqrt{1-\frac{4m^2}{s}}\sqrt{E_+^2-s}$$