Lots of integrals in life look like Gaussians.

\[ \int_{\mathbb{R}^n} \frac{e^{-\frac{|y|^2}{2}}}{(2\pi)^{n/2}} \mathrm{d}y = 1. \] In the heat kernel, we often have some rescaling $y \to \sqrt{\lambda} y$, scaling the integral to \[ \int_{\mathbb{R}^n} \frac{e^{-\frac{|y|^2}{2\lambda}}}{(2\pi \lambda)^{n/2}} \mathrm{d}y = 1, \] where perhaps we might take $\lambda = \frac{1}{2\pi}$ to get the familiar \[ \int_{\mathbb{R}^n} e^{-\pi |y|^2} \mathrm{d}y = 1. \] More generally, we might change coordinates by an invertible linear transformation $Sx = y$, whence $|y|^2 = |Sx|^2$ and \[ \int_{\mathbb{R}^n} \frac{e^{-\frac{1}{2} |S^{-1} x|^2}}{(2\pi)^{n/2}}\det(S) \mathrm{d}y = 1; \] that is, \[ \int_{\mathbb{R}^n} \frac{e^{-\frac{1}{2} |Sx|^2}}{\det(\sqrt{2\pi} S^{-1})} \mathrm{d}y = 1; \] Note that $|Sx|^2 = x^T A x$ with $A = S^T S$, and $\det(A) = \det(S)^2$; thus, for any symmetric $A$ \[ \int_{\mathbb{R}^n} \frac{e^{-\frac{1}{2} x^T A x}}{\sqrt{\det(2\pi A^{-1})}} \mathrm{d}x = 1; \] taking $A = \lambda I$ recovers our earlier formulae.

Each of these integrals is a measure on $\mathbb{R}^n$. Let $\mu_A$ denote the measure in this last formulation; a Gaussian integral is an integral $\int f(x) \mu_A(x)$ with respect to such a Gaussian measure. Let $\langle f \rangle_A$ denote this integral, or just $\langle f \rangle$ if $A$ is understood. Let $\mu_{\lambda}$ denote $\mu_{\lambda I}$, and let $\mu$ denote $\mu_1$.

Gaussian integrals are particularly nice to evaluate. There's the usual Calculus 1 problem \[ \langle x \rangle = 0 \] since $x e^{-\frac{x^2}{2}} = \frac{\mathrm{d}}{\mathrm{d}x} e^{-\frac{x^2}{2}}$ is a total derivative. We likewise compute \begin{align*} \langle x^2 \rangle &= \int_{-\infty}^\infty x \frac{\mathrm{d}}{\mathrm{d}x} e^{-\frac{x^2}{2}} \mathrm{d} x\\ &= \langle 1 \rangle = 1. \end{align*} Continuing this computation, clearly for all $n \ge 1$, $\langle x^{2n-1} \rangle = 0$ and \[ \langle x^{2n} \rangle = (2n-1) \langle x^{2n-2} \rangle = (2n-1)!! \] is a double factorial.

The goal is to extend this formula into the multivariate setting, and therefore find a rule to evaluate for any monomial the quantity $\langle x_{i_1} \cdots x_{i_k}\rangle$.