Speeding up NIntegrate and DensityPlot


I’m making a density plot of a spectral function Sfull[kx, ky, kz, \[CapitalOmega]] along a specified path in {kx,ky,kz}-space. The plot looks as it should, but it takes a long time to run (12+ hours). Sfull contains numerical integrals which I believe are the source of the problem. The evaluation of Sfull along the {kx,ky,kz}-space path could perhaps be made more efficient as well.

==============

Definition of Sfull[kx_, ky_, kz_, \[CapitalOmega]_]:

\[Theta] = 1.119; Qx = 0; Qy = 0; Qz = 0; \[Alpha] = 0.2; S = 0.5; J = 1.2; \[Gamma][kx_, ky_, kz_] := (Cos[kx] + Cos[ky] + \[Alpha] Cos[kz])/(2 + \[Alpha]); \[Omega]k[kx_, ky_, kz_] := Sqrt[(1 + \[Gamma][kx, ky, kz]) (1-Cos[2\[Theta]] \[Gamma][kx, ky, kz])]; \[Epsilon][kx_, ky_, kz_] := 2 J S (2 + \[Alpha]) \[Omega]k[kx, ky,kz]; Ak[kx_, ky_, kz_] := 2 J S (2 + \[Alpha]) (1 + Sin[\[Theta]]^2\[Gamma][kx, ky, kz]); u[kx_, ky_, kz_] := Sqrt[(Ak[kx, ky, kz] + \[Epsilon][kx, ky, kz])/(2 \[Epsilon][kx, ky,kz])]; v[kx_, ky_, kz_] := Sqrt[(Ak[kx, ky, kz] - \[Epsilon][kx, ky, kz])/(2 \[Epsilon][kx, ky,kz])]; \[CapitalPhi]1[kx_, ky_, kz_, qx_,qy_,qz_] := (\[Gamma] [kx, ky, kz] (u[kx, ky, kz] + v[kx, ky,kz]) (u[qx, qy, qz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[qx, qy, qz] u[kx - qx +Qx, ky - qy + Qy,kz - qz + Qz]) + \[Gamma][qx, qy, qz] (u[qx, qy, qz] + v[qx, qy, qz]) (u[kx, ky, kz] u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx, ky, kz] v[kx - qx + Qx, ky - qy + Qy,kz - qz + Qz]) + \[Gamma][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] (u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) (u[kx, ky, kz] u[qx, qy, qz] + v[kx, ky, kz] v[qx, qy, qz])); \[CapitalPhi]2[kx_, ky_, kz_, qx_,qy_,qz_] := \[Gamma][kx, ky, kz] (u[kx, ky, kz] + v[kx, ky, kz])(u[qx, qy, qz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[qx, qy, qz] u[kx - qx + Qx,ky - qy + Qy,kz - qz + Qz]) + \[Gamma][qx, qy, qz] (u[qx, qy, qz] + v[qx, qy, qz]) (u[kx, ky, kz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx, ky, kz] u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) + \[Gamma][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] (u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) (u[kx, ky, kz] v[qx, qy, qz] + v[kx, ky, kz] u[qx, qy, qz]); \[CapitalSigma][kx_, ky_, kz_, \[CapitalOmega]_]:=(1/2)NIntegrate[(-Abs[\[CapitalPhi]1[kx, ky, kz, qx, qy, qz]]^2/(\[CapitalOmega]-\[Epsilon][qx, qy, qz] - \[Epsilon][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + I 0.001) + Abs[\[CapitalPhi]2[kx, ky, kz, qx, qy, qz]]^2/(\[CapitalOmega] + \[Epsilon][qx, qy, qz] + \[Epsilon][kx + qx - Qx, ky + qy - Qy, kz + qz - Qz] - I 0.001)),{qy, 0, 2 \[Pi]}, {qz, 0, 2 \[Pi]}, Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 4, MaxRecursion -> 10^4]; G[kx_, ky_, kz_, \[CapitalOmega]_] := (1/(\[CapitalOmega] - \[Epsilon][kx, ky, kz] + \[CapitalSigma][kx, ky, kz, \[CapitalOmega]])); AA[kx_, ky_, kz_, \[CapitalOmega]_] := ((-1/\[Pi]) Im[G[kx, ky, kz,\[CapitalOmega]]]); Sxx[kx_, ky_, kz_, \[CapitalOmega]_] := (\[Pi] S (u[kx, ky, kz] +v[kx, ky, kz])^2 AA[kx, ky, kz, \[CapitalOmega]]); Syy[kx_, ky_, kz_, \[CapitalOmega]_] := (\[Pi] S (u[kx, ky, kz] -v[kx, ky, kz])^2 AA[kx, ky, kz, \[CapitalOmega]]); Szz[kx_, ky_, kz_, \[CapitalOmega]_] := \[Pi]NIntegrate[((u[qx, qy, qz] v[kx - qx, ky - qy, kz - qz] + v[qx, qy, qz] u[kx - qx, ky -qy,kz - qz])^2 (1/Sqrt[2 \[Pi]] (0.1 J)) Exp[-(\[CapitalOmega] - \[Epsilon][qx,qy, qz] -\[Epsilon][kx - qx, ky - qy, qz - kz])^2/(2 (0.1 J)^2)]), {qx, 0, 2 \[Pi]}, {qy, 0, 2 \[Pi]},{qz,0, 2 \[Pi]}, Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 4, MaxRecursion -> 10^4]; Sxx0[kx_, ky_, kz_, \[CapitalOmega]_] := (Sin[\[Theta]]^2 Sxx[kx,ky,kz, \[CapitalOmega]] + Cos[\[Theta]]^2 Szz[kx - Qx, ky - Qy, kz - Qz, \[CapitalOmega]]); Szz0[kx_, ky_, kz_, \[CapitalOmega]_] := (Cos[\[Theta]]^2 Sxx[kx - Qx, ky - Qy, kz - Qz, \[CapitalOmega]] + Sin[\[Theta]]^2 Szz[kx, ky, kz, \[CapitalOmega]]); Syy0[kx_, ky_, kz_, \[CapitalOmega]_] := (Syy[kx, ky,kz,\[CapitalOmega]]); Sfull[kx_, ky_, kz_, \[CapitalOmega]_] := Sxx0[kx, ky, kz, \[CapitalOmega]] + Syy0[kx, ky, kz, \[CapitalOmega]] + Szz0[kx, ky, kz, \[CapitalOmega]]; 

Definition of {kx,ky,kz}-space path:

kpath = {{0, {0, 0, 0, \[CapitalOmega]}}, {1, {\[Pi], \[Pi], 0, \[CapitalOmega]}}, {2, {\[Pi], 0, 0, \[CapitalOmega]}}, {3, {0, \[Pi], 0, \[CapitalOmega]}}, {4, {0, 0, 0, \[CapitalOmega]}}}; if = Interpolation[kpath, InterpolationOrder -> 1];  

Density plot:

DensityPlot[{Sfull@@if[i]}, {i, 0, 4}, {\[CapitalOmega], 0, 10},ColorFunction -> "SunsetColors", PlotLegends -> Automatic] 

This produces the plot

enter image description here

which takes over 12 hours to run. \[CapitalSigma][kx, ky, kz, \[CapitalOmega]] and Szz[kx, ky, kz, \[CapitalOmega]] contain the numerical integrals (which I would like to perform according to the adaptive quasi Monte Carlo scheme with AccuracyGoal -> 4, MaxRecursion -> 10^4) and I’d like to speed these up. I’m also unsure that using Sfull@@if[i] is the most efficient way to produce the density plot. I would eventually like to make the plot smoother by increasing PlotPoints, but this would of course increase the runtime further.

I have tried SymbolicPreprocessing->0 and defining functions with ?NumericQ variables but this does not speed up the integrals noticably.