# multiple Numeric integral calculation lasts forever

I'm plotting some multiple numeric integral and after 16-17 hours(i7-6700k 4.2 GHz) I didn't get the result. I tried to avoid all secondary calculations, input functions directly, disabled Symbolic Processing, made scaling of the problem so integrand value is magnitude of 1. The precision in 8 digits forced because of value of parameter $e$, function $\Sigma$ inside the integral diverges if it goes zero. If you uncomment plots before last one, you see that all functions have nice behaviour in integration region. Further I will need to integrate this again with some distribution over $p$. Any suggestions and recommendations? Code:

int = -((E^(-2 r (Sqrt[u (-4 k^2 (-1 + u) + x - y) + y] +
I k (1 - 2 u) Cos[θ])) (-1 + u) u (3 +
2 r (3 Sqrt[u (-4 k^2 (-1 + u) + x - y) + y] +
2 r (u (-4 k^2 (-1 + u) + x - y) + y))))/(
2 π^3 Sqrt[
u (-4 k^2 (-1 + u) + x - y) +
y] (4 k^2 (-1 + u) u - u x + (-1 + u) y)^2))
x = 1
y = 1
W[r_?NumericQ, k_?NumericQ, θ_?NumericQ] :=
NIntegrate[Re[int]/scaleW, {u, 0, 1},
Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
PrecisionGoal -> 8, WorkingPrecision -> 8]
Σ[p_?NumericQ, k_?NumericQ, θ_?NumericQ] :=
1/scaleΣ*(2*a^2*m^2*
Pi*α^4*(2*
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])] + ((k^2 - (e - m)^2)^2 +
m^4 -
2*(-k^2 + (e - m)^2)*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ]) +
2*(-k^2 + (e - m)^2 + m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])^2 +

6*m^2*(-m^2 - 2*(-e + m)*Sqrt[m^2 + p^2] +
2*k*p*Cos[θ]))/(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[θ] -

Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]) - ((k^2 - (e -
m)^2)^2 + m^4 -
2*(-k^2 + (e - m)^2)*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ]) +
2*(-k^2 + (e - m)^2 + m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])^2 +

6*m^2*(-m^2 - 2*(-e + m)*Sqrt[m^2 + p^2] +
2*k*p*Cos[θ]))/(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[θ] +

Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]) +
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] - k*p*Cos[θ])*

Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[θ] -

Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] -
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] - k*p*Cos[θ])*

Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[θ] +

Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] - (1/(m^2 + (e -
m)*Sqrt[m^2 + p^2] + k*p*Cos[θ]))*
((12*m^4 -
8*m^2*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[θ]) +
(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[θ])^2)*
(-Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[θ] -
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) -
2*k*p*Cos[θ]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/
(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] +

Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[θ] +

Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] +
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] -

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] -
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] +

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]])) + (2*
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]*(e^4 - 2*e^2*k^2 +
k^4 - 4*e^3*m + 4*e*k^2*m + 8*e^2*m^2 - 4*k^2*m^2 -
8*e*m^3 + m^4 -
4*e^3*Sqrt[m^2 + p^2] +
4*e*k^2*Sqrt[m^2 + p^2] + 12*e^2*m*Sqrt[m^2 + p^2] -
4*k^2*m*Sqrt[m^2 + p^2] -
6*e*m^2*Sqrt[m^2 + p^2] -
2*m^3*Sqrt[m^2 + p^2] + 9*e^2*(m^2 + p^2) -
18*e*m*(m^2 + p^2) + 9*m^2*(m^2 + p^2) -
4*e^2*k*p*Cos[θ] +
4*k^3*p*Cos[θ] + 8*e*k*m*p*Cos[θ] +
2*k*m^2*p*Cos[θ] +
18*e*k*p*Sqrt[m^2 + p^2]*Cos[θ] -
18*k*m*p*Sqrt[m^2 + p^2]*Cos[θ] +
9*k^2*p^2*
Cos[θ]^2 - (e^2 - k^2 - 2*m^2 +
2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) -
2*k*p*Cos[θ])*(-m^2 + (m^2 + (-e + m)*
Sqrt[m^2 + p^2] - k*p*Cos[θ])^2/
(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ]))) +
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] - k*p*Cos[θ])*
(m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] -

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])])*(m^2 +
e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] +

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])])*(-Log[
m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] -

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]] +
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] +

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*
Cos[θ])]]))/((m^2 +
e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] -

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])])*(m^2 +
e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[θ] +

Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[θ]]*

Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[θ])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])]))))/((-k^2 + (e -
m)^2 - 3*m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ])*
(-k^2 + (e - m)^2 + m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[θ]))

m = 0.5109989461
α = 1/137
e = (α^2 m)/2
a = 1/(α*m)
scaleΣ = 10^6
scaleW = 10^-3
(*Plot[Σ[p=0.5,k/a,θ=0],{k,0,5}]
Plot[Σ[p,1/a,θ=0],{p,0.05,1.1}]
Plot[Σ[p=0.5,1/a,θ],{θ,0,Pi}]
Plot3D[r^2*k^2*W[r,k,θ=0],{r,0,5},{k,0,5},PlotRange\[Rule]All,\
AxesLabel\[Rule]Automatic]
Plot3D[r^2*k^2*W[r,k,θ=Pi/2],{r,0,5},{k,0,5},PlotRange\[Rule]\
All,AxesLabel\[Rule]Automatic]*)
Plot[(*scaleΣ*scaleW*2Pi*a^-2**)
NIntegrate[Σ[p, k/a, θ]*W[r, k, θ]*
Sin[θ]*k^2, {θ, 0, Pi}, {k, 0, 5}, {r, 0, 5},
Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
PrecisionGoal -> 8, WorkingPrecision -> 8], {p, 0.055, 1.1}]
Clear[x, y, int, W, k, r, θ, c, u, m, e, a, p, σ, \
Σ, t1, t0, M, t, s]


Edit: I noticed that integral can be done analytically by $r$ so integration dimension reduces to 3. Now it looks like $\int_{0}^{5}\int_{0}^{\pi}\Sigma(p,k,\theta)\int_0^1duf(k,\theta,u)dkd\theta$. It can be done numerically one time by $u$ and I plotted integrand by $k$ and $\theta$ with fixed $p$. It is still very smooth function and should be simple to plot over $p$ but it does not. I have removed precision requirements and decreased digits in $m$. Computation in a fixed $p$ required about 100 sec. For $p=0.5$ I got {100.716, 22.1118}. Any estimates how much time I need to wait to get plot over $p$?

f[k_?NumericQ, \[Theta]_?NumericQ, u_?NumericQ] :=
Re[-(1/(scaleW*2 \[Pi]^3 (1 - 4 k^2 (-1 + u) u)^(5/2))) (-1 +
u) u* (-((-16 + 3 k^2 - 76 k^2 u + 76 k^2 u^2 +
18 I k (-1 + 2 u) Sqrt[1 - 4 k^2 (-1 + u) u] Cos[\[Theta]] +
3 k^2 (1 - 2 u)^2 Cos[2 \[Theta]])/(
4 (Sqrt[1 - 4 k^2 (-1 + u) u] -
I k (-1 + 2 u) Cos[\[Theta]])^3)) +
E^(-10 Sqrt[1 - 4 k^2 (-1 + u) u] +
10 I k (-1 + 2 u) Cos[\[Theta]]) ((
50 I (-1 + 4 k^2 (-1 + u) u))/(
I Sqrt[1 - 4 k^2 (-1 + u) u] + k (-1 + 2 u) Cos[\[Theta]]) + (
5 (5 (-1 + 4 k^2 (-1 + u) u) +
3 I k (-1 + 2 u) Sqrt[1 - 4 k^2 (-1 + u) u]
Cos[\[Theta]]))/(Sqrt[1 - 4 k^2 (-1 + u) u] -
I k (-1 + 2 u) Cos[\[Theta]])^2 + (-16 + 3 k^2 - 76 k^2 u +
76 k^2 u^2 +
18 I k (-1 + 2 u) Sqrt[1 - 4 k^2 (-1 + u) u] Cos[\[Theta]] +
3 k^2 (1 - 2 u)^2 Cos[2 \[Theta]])/(
4 (Sqrt[1 - 4 k^2 (-1 + u) u] -
I k (-1 + 2 u) Cos[\[Theta]])^3)))]
\[CapitalSigma][p_?NumericQ, k_?NumericQ, \[Theta]_?NumericQ] :=
1/scale\[CapitalSigma]*(2*a^2*m^2*
Pi*\[Alpha]^4*(2*
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])] + ((k^2 - (e - m)^2)^2 + m^4 -
2*(-k^2 + (e - m)^2)*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[\[Theta]]) +
2*(-k^2 + (e - m)^2 + m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])^2 +
6*m^2*(-m^2 - 2*(-e + m)*Sqrt[m^2 + p^2] +
2*k*p*Cos[\[Theta]]))/(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]) - ((k^2 - (e - m)^2)^2 + m^4 -

2*(-k^2 + (e - m)^2)*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[\[Theta]]) +
2*(-k^2 + (e - m)^2 + m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])^2 +
6*m^2*(-m^2 - 2*(-e + m)*Sqrt[m^2 + p^2] +
2*k*p*Cos[\[Theta]]))/(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[\[Theta]])]) +
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] - k*p*Cos[\[Theta]])*
Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[\[Theta]])]] -
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] - k*p*Cos[\[Theta]])*
Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]] - (1/(m^2 + (e - m)*
Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]]))*((12*m^4 -
8*m^2*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]]) + (-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])^2)*(-Log[-m^2 -
e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]] +

Log[-m^2 - e*Sqrt[m^2 + p^2] + m*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 - 2*m^2 + 2*m*Sqrt[m^2 + p^2] -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]] +
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]] -
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]])) + (2*
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]*(e^4 - 2*e^2*k^2 + k^4 -
4*e^3*m + 4*e*k^2*m + 8*e^2*m^2 - 4*k^2*m^2 - 8*e*m^3 +
m^4 - 4*e^3*Sqrt[m^2 + p^2] + 4*e*k^2*Sqrt[m^2 + p^2] +
12*e^2*m*Sqrt[m^2 + p^2] - 4*k^2*m*Sqrt[m^2 + p^2] -
6*e*m^2*Sqrt[m^2 + p^2] - 2*m^3*Sqrt[m^2 + p^2] +
9*e^2*(m^2 + p^2) - 18*e*m*(m^2 + p^2) +
9*m^2*(m^2 + p^2) - 4*e^2*k*p*Cos[\[Theta]] +
4*k^3*p*Cos[\[Theta]] + 8*e*k*m*p*Cos[\[Theta]] +
2*k*m^2*p*Cos[\[Theta]] +
18*e*k*p*Sqrt[m^2 + p^2]*Cos[\[Theta]] -
18*k*m*p*Sqrt[m^2 + p^2]*Cos[\[Theta]] +
9*k^2*p^2*
Cos[\[Theta]]^2 - (e^2 - k^2 - 2*m^2 +
2*m*Sqrt[m^2 + p^2] - 2*e*(m + Sqrt[m^2 + p^2]) -
2*k*p*Cos[\[Theta]])*(-m^2 + (m^2 + (-e + m)*
Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]]))) +
4*(m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])*(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])])*(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])])*(-Log[
m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]] +
Log[m^2 + e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]]))/((m^2 +
e*Sqrt[m^2 + p^2] - m*Sqrt[m^2 + p^2] +
k*p*Cos[\[Theta]] -
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])])*(m^2 + e*Sqrt[m^2 + p^2] -
m*Sqrt[m^2 + p^2] + k*p*Cos[\[Theta]] +
Sqrt[e^2 - k^2 + 2*m*(-m + Sqrt[m^2 + p^2]) -
2*e*(m + Sqrt[m^2 + p^2]) - 2*k*p*Cos[\[Theta]]]*
Sqrt[-m^2 + (m^2 + (-e + m)*Sqrt[m^2 + p^2] -
k*p*Cos[\[Theta]])^2/(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])]))))/((-k^2 + (e - m)^2 -
3*m^2 + 2*(-e + m)*Sqrt[m^2 + p^2] -
2*k*p*Cos[\[Theta]])*(-k^2 + (e - m)^2 + m^2 +
2*(-e + m)*Sqrt[m^2 + p^2] - 2*k*p*Cos[\[Theta]]))
m = 0.5109
\[Alpha] = 1/137
e = (\[Alpha]^2 m)/2
a = 1/(\[Alpha]*m)
scale\[CapitalSigma] = 10^6
scaleW = 10^-3
Plot3D[Sin[\[Theta]]*k^2*\[CapitalSigma][p = 0.5, k/a, \[Theta]]*
NIntegrate[f[k, \[Theta], u], {u, 0, 1}], {k, 0, 5}, {\[Theta], 0,
Pi}, PlotRange -> All]
NIntegrate[(*scale\[CapitalSigma]*scaleW*2Pi*
a^-2*)\[CapitalSigma][p = 0.5, k/a, \[Theta]]*f[k, \[Theta], u]*
Sin[\[Theta]]*k^2, {\[Theta], 0, Pi}, {k, 0, 5}, {u, 0,
1}] // AbsoluteTiming
Plot[NIntegrate[(*scale\[CapitalSigma]*scaleW*2Pi*
a^-2*)\[CapitalSigma][p, k/a, \[Theta]]*f[k, \[Theta], u]*
Sin[\[Theta]]*k^2, {\[Theta], 0, Pi}, {k, 0, 5}, {u, 0, 1}], {p,
0.05, 1.1}]
Clear[x, y, int, W, k, r, \[Theta], c, u, m, e, a, p, \[Sigma], \
\[CapitalSigma], t1, t0, M, t, s, f, \[Alpha]]

• have you got that NIntegrate to evaluate at all? If it works at all you should generate a table over p and use ListPlot. You can often get a reasonable plot with far fewer evaluations than Plot will need. – george2079 Apr 11 '17 at 20:46
• I tried to calculate it in a single point instead of plotting and it still running for about 40 mins... for single point... Tomorrow I will update my question, looks like things are going really bad. – satoru Apr 11 '17 at 21:56
• You have this options PrecisionGoal -> 8, WorkingPrecision -> 8 -- no wonder your computations are slow. Use machine precision for faster results. To start investigating, how quickly you get results with PrecisionGoal -> 2, WorkingPrecision -> MachinePrecision ? – Anton Antonov Apr 12 '17 at 15:47
• @AntonAntonov The value of parameter $e$ is 10^-6. If precision is too small this will be zero and function $\Sigma$ inside the integral diverges – satoru Apr 12 '17 at 15:58
• do you realize WorkingPrecision->8 is much smaller precision than default? What you are saying about singularity makes no sense. That said the integral doesn't evaluate in any reasonable time at MachinePrecision either. – george2079 Apr 12 '17 at 18:40