I'm trying to solve numerically the following double integral in Mathematica, which I know, a priori, has a solution for values of $mh \in[0.4,2.0]$ and tmer = 0.017, tmel = 0.044.
ID[mh_] :=
1/(2*Pi^3)*
NIntegrate[(mh*tmer^(1/
2)*(tmer^(1/2)*Sqrt[w^2 - 1]*
Sqrt[(u*mh - w*tmer^(1/2))^2 - tmel] +
1/2*(mh^2 - tmer - tmel) + (w*tmer^(1/2))^2 -
u*w*mh*tmer^(1/2)))*
Exp[u*mh]/((Exp[u*mh] - 1)*(Exp[w*tmer^(1/2)] +
1)*(Exp[u*mh - w*tmer^(1/2)] + 1)), {u, 1, +\[Infinity]}, {w,
1, +\[Infinity]}]
If you run, for instance, ID[0.4] you'll see that the numerical integration is converging too slowly, I suppose due to the Exp[u*mh] term in the numerator. Furthemore, it returns a complex number, thanks to Sqrt[(u*mh - w*tmer^(1/2))^2 - tmel]. However, I want to impose that (u*mh - w*tmer^(1/2))^2 >= tmel. If done correctly, I would expect ID[0.4] to output 0.000025385 or close to it. Is there a way to deal with the mentioned slow convergence and impose the above condition?
