I’m trying to numerically integrate for $ \omega > 0$ the following nine dimensional integral $ $ \iiint_{-\infty}^\infty d^3\vec{k} d^3\vec{k_1} d^3\vec{k_2}\, e^\left(\frac{1}{8} \left(-k_{y2}^2-(-k_{y}+k_{y1}+k_{y2})^2-2 \left(k_{y}^2+k_{y1}^2+k_{z}^2+k_{z1}^2\right)-k_{z2}^2-(-k_{z}+k_{z1}+k_{z2})^2\right)\right) \cdot (\vec{k} \cdot \vec{k_1}) \sqrt{\left(k_{y}^2+k_{z}^2\right) \left(k_{y1}^2+k_{z1}^2\right)} \left(k_{x2}^2+3 \left(k_{y2}^2+k_{z2}^2\right)\right) \cdot \left(\frac{2 \left((-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2\right)}{\left((k_{x}+k_{x1}-k_{x2})^2+(-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2\right)^2}+\frac{1}{(k_{x}+k_{x1}-k_{x2})^2+(-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2}+\frac{2 \left((-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2\right)}{\left((-k_{x}+k_{x1}+k_{x2})^2+(-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2\right)^2}+\frac{1}{(-k_{x}+k_{x1}+k_{x2})^2+(-k_{y}+k_{y1}+k_{y2})^2+(-k_{z}+k_{z1}+k_{z2})^2}\right) \cdot \left(\left(k_{x}^2+k_{y}^2+k_{z}^2\right) \cdot \left(k_{x1}^2+k_{y1}^2+k_{z1}^2\right)+\omega ^2\right) \cdot \frac{1}{\left(k_{x}^2+k_{y}^2+k_{z}^2\right) \left(k_{x1}^2+k_{y1}^2+k_{z1}^2\right) \left(k_{x2}^2+k_{y2}^2+k_{z2}^2\right)^2 \pi ^2 \left(\left(k_{x}^2+k_{y}^2+k_{z}^2\right)^2+\omega ^2\right) \left(\left(k_{x1}^2+k_{y1}^2+k_{z1}^2\right)^2+\omega ^2\right)}. $ $

The code for this integral is as follows:

Int[\[Omega]_] :=  NIntegrate[(E^(     1/8 (-ky2^2 - (-ky + ky1 + ky2)^2 -         2 (ky^2 + ky1^2 + kz^2 + kz1^2) -         kz2^2 - (-kz + kz1 + kz2)^2)) (kx kx1 + ky ky1 +        kz kz1) Sqrt[(ky^2 + kz^2) (ky1^2 + kz1^2)] (kx2^2 +        3 (ky2^2 + kz2^2)) ((       2 ((-ky + ky1 + ky2)^2 + (-kz + kz1 + kz2)^2))/((kx + kx1 -            kx2)^2 + (-ky + ky1 + ky2)^2 + (-kz + kz1 + kz2)^2)^2 +        1/((kx + kx1 - kx2)^2 + (-ky + ky1 + ky2)^2 + (-kz + kz1 +           kz2)^2) + (       2 ((-ky + ky1 + ky2)^2 + (-kz + kz1 + kz2)^2))/((-kx + kx1 +            kx2)^2 + (-ky + ky1 + ky2)^2 + (-kz + kz1 + kz2)^2)^2 +        1/((-kx + kx1 + kx2)^2 + (-ky + ky1 + ky2)^2 + (-kz + kz1 +           kz2)^2)) ((kx^2 + ky^2 + kz^2) (kx1^2 + ky1^2 +           kz1^2) + \[Omega]^2))/((kx^2 + ky^2 + kz^2) (kx1^2 + ky1^2 +        kz1^2) (kx2^2 + ky2^2 +        kz2^2)^2 \[Pi]^2 ((kx^2 + ky^2 + kz^2)^2 + \[Omega]^2) ((kx1^2 +          ky1^2 +          kz1^2)^2 + \[Omega]^2)),      {kx, -\[Infinity], \[Infinity]}, {ky, -\[Infinity], \[Infinity]},      {kz, -\[Infinity], \[Infinity]}, {kx1, -\[Infinity], \[Infinity]},      {ky1, -\[Infinity], \[Infinity]}, {kz1, -\[Infinity], \[Infinity]},      {kx2, -\[Infinity], \[Infinity]}, {ky2, -\[Infinity], \[Infinity]},      {kz2, -\[Infinity], \[Infinity]},      PrecisionGoal -> 4, WorkingPrecision -> MachinePrecision] 

Due to the high number of dimensions, and the singular nature of the integral, I’m finding it difficult to make it converge. Furthermore, I don’t see any symmetry that I could exploit to reduce the number of dimensions. Is there any way to make NIntegrate more tractable in this case?

Thanks.