Although the original question is a little too narrow, I'm going to interpret it as asking about general multidimensional integrals in which Gaussians appear together with arbitrary factors in the integrand.

Then you may be interested in this observation which speeds up such integrals tremendously:

Take for example the Fourier transform of a Gaussian,

$$\iiint \exp(- i\, \vec{k}\cdot \vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

Here $\Sigma^{-1}$ could be a symmetric positive definite matrix, or in the simplest case just $\sigma^2$ times the identity matrix.

Instead of feeding this integral to Integrate or NIntegrate with infinite integration boundaries, I found that it's much more efficient to evaluate such Gaussian integrals as follows:

For Gaussian integrals over all space (or momentum space, as in the question), this approach is complementary to whuber's solution, and in fact both faster and simpler.

Although the original question is a little too narrow, I'm going to interpret it as asking about general multidimensional integrals in which Gaussians appear together with arbitrary factors in the integrand:

$$\iiint f(\,\vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

Here $\Sigma^{-1}$ could be a symmetric positive definite matrix, or in the simplest case just $\sigma^2$ times the identity matrix. The dimension can be anything, but here I'll just specialize to three.

Then you may be interested in this observation which speeds up such integrals tremendously in comparison with straightforward multiple integration:

Edit: An integrand that doesn't factor

Start by defining the Gaussian in three dimensions:

f = PDF[MultinormalDistribution[{0, 0, 
    0}, \[Sigma] IdentityMatrix[3]], {x, y, z}]

$\frac{e^{\frac{1}{2} \left(-\frac{x^2}{\sigma }-\frac{y^2}{\sigma }-\frac{z^2}{\sigma }\right)}}{2 \sqrt{2} \pi ^{3/2} \sqrt{\sigma ^3}}$

Try the following as a triple integral:

Timing[Integrate[
  Integrate[Integrate[
    (x^2 + y^3)/(1 + z^2)
      f, {x, -Infinity, Infinity}], {y, -Infinity, \
Infinity}], {z, -Infinity, Infinity}]
 ]

$\left\{1.48519,\text{ConditionalExpression}\left[\frac{\sqrt{\frac{ \pi }{2}} e^{\frac{1}{2 \sigma }} \sqrt{\sigma ^3} \text{erfc}\left(\frac{1}{\sqrt{ 2} \sqrt{\sigma }}\right)}{\sigma },\Re(\sigma )>0\right]\right\}$

Now do the same using Expectation, keeping in mind the different normalization of the result:

Timing[Expectation[(x^2 + y^3)/(
  1 + z^2), {x, y, z} \[Distributed] 
   MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]]]
 ]

$\left\{0.775868,\sqrt{\frac{\pi }{2}} e^{\frac{1}{2 \sigma }} \sqrt{\sigma } \text{erfc}\left(\frac{1}{\sqrt{ 2} \sqrt{\sigma }}\right)\right\}$

So the result is obtained about twice as fast using Expectation.

Fourier transforms as a special case

My initial example was the Fourier transform of a Gaussian,

$$\iiint \exp(- i\, \vec{k}\cdot \vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

where the speed advantage of Expectation compared to multiple Integrate with infinite integration boundaries is more significant than above:

Although my initial claim is correct that this is faster than the straightforward triple integral, one can do a lot better in this special case by just using FourierTransform directly:

Timing[FourierTransform[f, {x, y, z}, {kx, ky, kz}]]

$\left\{0.007332,\frac{e^{-\frac{1}{2} \sigma \left(\text{kx}^2+\text{ky}^2+\text{kz}^2\right)} }{2 \sqrt{2} \pi ^{3/2} \left(\frac{1}{\sigma }\right)^{3/2} \sqrt{\sigma ^3}}\right\}$

Interestingly, though, we can even beat this really fast result using the multinormal distribution again:

Timing[CharacteristicFunction[
  MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]], {kx,
    ky, kz}]
 ]

$\left\{0.000185,e^{\frac{1}{2} \left(\text{kx}^2 (-\sigma )-\text{ky}^2 \sigma -\text{kz}^2 \sigma \right)}\right\}$

Conclusion

For Gaussian integrals over all space (or momentum space, as in the question), the approach using MultinormalDistribution is complementary to whuber's solution: general Gaussian integrals can be evaluated by using Expectation and similar tools for probability distributions, such as CharacteristicFunction.