reference for multidimensional gaussian integral

I was reading on Wikipedia in this article about the n-dimensional and functional generalization of the Gaussian integral. In particular, I would like to understand how the following equations are derived:
$$
\begin{eqnarray}
& {} \quad \int x^{k_1}\cdots x^{k_{2N}} \, \exp\left( – \frac 1 2 \sum_{i,j=1}^{n}A_{ij} x_i x_j \right) \, d^nx \\
& =
\sqrt{\frac{(2\pi)^n}{\det A}} \, \frac{1}{2^N N!} \, \sum_{\sigma \in S_{2N}}(A^{-1})^{k_{\sigma(1)}k_{\sigma(2)}} \cdots (A^{-1})^{k_{\sigma(2N-1)}k_{\sigma(2N)}}
\end{eqnarray}
$$

and

$$
\begin{eqnarray}
\int f(\vec x) \, \exp\left( – \frac 1 2 \sum_{i,j=1}^{n}A_{ij} x_i x_j \right) d^nx
=
\sqrt{(2\pi)^n\over \det A} \, \left. \exp\left({1\over 2}\sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right)f(\vec{x})\right|_{\vec{x}=0}
\end{eqnarray}
$$

The problem is that I can’t find a book or internet source that shows a complete derivation of this .. any reference tip would be highly appreciated!

Another thing that I am not so sure about is what the complete symmetrization means in the above equation, if someone could give a short explanation (or a more detailed equation) of what the following expression means that would be great:

$$
\begin{eqnarray}
\sum_{\sigma \in S_{2N}}(A^{-1})^{k_{\sigma(1)}k_{\sigma(2)}} \cdots (A^{-1})^{k_{\sigma(2N-1)}k_{\sigma(2N)}}
\end{eqnarray}
$$

Many thanks !!

Solutions Collecting From Web of "reference for multidimensional gaussian integral"

The presentation here is typical of those used to model and motivate the infinite dimensional Gaussian integrals encountered in quantum field theory.
I will use subscripts instead of superscripts to indicate components.

I. Wick’s theorem

First consider the integral
$$Z_0 = \int d^n x \exp\left(-\frac{1}{2} x^\mathrm{T} A x\right),$$
where $d^n x = \prod_i d x_i$ and $A$ is symmetric and positive definite.
Diagonalize $A$ with an orthogonal transformation.
Letting $D = M^\mathrm{T} A M = \mathrm{diag}(\lambda_1,\ldots,\lambda_n)$ and $z = M^\mathrm{T} x$, and noting that the Jacobian is the identity, we find
$$\begin{eqnarray}
Z_0 &=& \int d^n z \exp\left(-\frac{1}{2} z^\mathrm{T} D z\right) \\
&=& \prod_i \int d z_i \exp\left(-\frac{1}{2} \lambda_i z_i^2\right) \\
&=& \prod_i \sqrt{\frac{2\pi}{\lambda_i}} \\
&=& \sqrt{\frac{(2\pi)^n}{\det A}}.
\end{eqnarray}$$

Add a source term
$$Z_J = \int d^n x \exp\left(-\frac{1}{2} x^\mathrm{T} A x + J^\mathrm{T} x\right).$$
Complete the square to eliminate the cross term.
Letting $x = y + A^{-1}J$, we find
$$\begin{eqnarray}
Z_J &=& \int d^n y \exp\left(-\frac{1}{2} {y}^\mathrm{T} A y
+ \frac{1}{2} J^\mathrm{T}A^{-1}J\right) \\
&=& \sqrt{\frac{(2\pi)^n}{\det A}}
\exp\left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right).
\end{eqnarray}$$
There is always a factor of $\sqrt{\frac{(2\pi)^n}{\det A}}$.
For convenience, let’s define
$$\langle x_{k_1} \cdots x_{k_{2N}}\rangle = \frac{1}{Z_0}
\int d^n x \ x_{k_1} \cdots x_{k_{2N}} \exp\left(-\frac{1}{2} x^\mathrm{T} A x\right).$$
(Notice that $\langle x_{k_1} \cdots x_{k_{2N+1}}\rangle = 0$ since the integral is odd.)
Roughly speaking, these are free field scattering amplitudes.
This integral can be found by taking derivatives of $Z_J$,
$$\langle x_{k_1} \cdots x_{k_{2N}}\rangle =
\left.\frac{1}{Z_0} \frac{\partial}{\partial J_{k_1}} \cdots \frac{\partial}{\partial J_{k_{2N}}} Z_J\right|_{J=0}.$$
(Notice that every factor of $\partial/\partial_{J_{k_i}}$ brings down one factor of $x_{k_i}$ in the integral $Z_J$.)
Using this formula it is a straightforward exercise to work out, for example, that
$$\begin{eqnarray}
\langle x_{k_1} x_{k_{2}}\rangle
&=& \left.\frac{\partial}{\partial J_{k_1}} \frac{\partial}{\partial J_{k_{2}}}
\exp\left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right)\right|_{J=0} \\
&=& \frac{\partial}{\partial J_{k_1}} \frac{\partial}{\partial J_{k_{2}}}
\frac{1}{2} J^\mathrm{T}A^{-1}J \\
&=& \frac{1}{2} ( A^{-1}_{k_1 k_2} + A^{-1}_{k_2 k_1}) \\
&=& A^{-1}_{k_1 k_2},
\end{eqnarray}$$
which agrees with the formula for $N=1$.
This is the “free field propagator.”
(In the last step we have used the fact that $A^{-1}$ is symmetric.)
It is possible to see by inspection that the formula for general $N$ is
$$\langle x_{k_1}\cdots x_{k_{2N}}\rangle = \frac{1}{2^N N!}
\sum_{\sigma \in S_{2N}}(A^{-1})_{k_{\sigma(1)}k_{\sigma(2)}} \cdots (A^{-1})_{k_{\sigma(2N-1)}k_{\sigma(2N)}}.$$
In fact
$$\begin{eqnarray}
\langle x_{k_1}\cdots x_{k_{2N}}\rangle &=&
\left.\frac{\partial}{\partial J_{k_1}} \cdots \frac{\partial}{\partial J_{k_{2N}}}
\exp \left(\frac{1}{2} J^\mathrm{T}A^{-1}J \right) \right|_{J=0} \\
&=& \frac{\partial}{\partial J_{k_1}} \cdots \frac{\partial}{\partial J_{k_{2N}}}
\frac{1}{N!} \left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right)^{N} \\
&=& \frac{1}{2^N N!} \frac{\partial}{\partial J_{k_1}} \cdots \frac{\partial}{\partial J_{k_{2N}}}
\left(J^\mathrm{T}A^{-1}J\right)^{N}.
\end{eqnarray}$$
Note that the derivative $\partial/\partial J_{k_1}$ will operate on all possible $J$s.
Likewise for the other derivatives.
Thus, we will get a sum over all $(2N)!$ possible permutations of the $k_i$.
These permutations are denoted by $\sigma$ and they live in the symmetric group $S_{2N}$.
Thus, we arrive at the result, known as Wick’s theorem.
(If this seems too vague, it is straightforward to find the result by induction on $N$.
We have the formula for $N=1$ above.)

Let’s unwind the scattering amplitude for $N=2$.
We have
$$\langle x_{k_1}x_{k_2}x_{k_3}x_{k_4}\rangle = \frac{1}{2^2 2!}
\sum_{\sigma \in S_{4}}(A^{-1})_{k_{\sigma(1)}k_{\sigma(2)}} (A^{-1})_{k_{\sigma(3)}k_{\sigma(4)}}.$$
There are $4! = 24$ permutations in $S_4$ and thus $24$ terms in the sum.
For example $(12)$ takes $1$ to $2$ and $2$ to $1$.
Not all permutations give independent terms.
In fact, the degeneracy is $2^N N!$, since $A^{-1}$ is symmetric and the order of the $A^{-1}$s doesn’t matter.
Thus, there will only be
$$\frac{(2N)!}{2^N N!} = (2N-1)!!$$
independent terms.
For $N=2$ there are three independent terms,
$$\langle x_{k_1}x_{k_2}x_{k_3}x_{k_4}\rangle
= A^{-1}_{k_1 k_2}A^{-1}_{k_3 k_4}
+ A^{-1}_{k_1 k_3}A^{-1}_{k_2 k_4}
+ A^{-1}_{k_1 k_4}A^{-1}_{k_2 k_3}.$$

II. Central identity

Consider
$$I_J = \frac{1}{Z_0} \int d^n x \exp\left(-\frac{1}{2}x^\mathrm{T}A x + J^\mathrm{T}x \right) f(x).$$
We are interested in $I_0$.
The presence of the source allows us to take $f$ out of the integral if we replace its argument with $\partial/\partial J$,
$$\begin{eqnarray}
I_0 &=& \left.f(\partial_J) \frac{1}{Z_0} \int d^n x \exp\left(-\frac{1}{2}x^\mathrm{T}A x + J^\mathrm{T}x \right)\right|_{J=0} \\
&=& \left.f(\partial_J) \exp\left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right)\right|_{J=0}.
\end{eqnarray}$$
This is a typical trick.
In fact, it is equivalent to what Anthony Zee calls the “central identity of quantum field theory.”
Usually $f(x) = \exp[-V(x)]$, where $V(x)$ is the potential.
There is a nice graphical interpretation of the formula in the form of Feynman diagrams.
The process of calculating
$$\left.e^{-V(\partial_J)} \exp\left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right)\right|_{J=0}$$
is equivalent to “tying together” the propagators (the $A^{-1}$) with vertices represented by the operator $-V(\partial_J)$.
Of course, there is a lot more to it than that!

Now we wish to show
$$I_0 = \left.\exp\left(\frac{1}{2}\partial_x^\mathrm{T} A^{-1}\partial_x\right) f(x)\right|_{x=0}.$$
If we consider a Taylor expansion for $f(\partial_J)$ and $f(x)$ we can see this is equivalent to
$$\left.\partial_{J_{k_1}}\cdots \partial_{J_{k_{2N}}} \exp\left(\frac{1}{2} J^\mathrm{T}A^{-1}J\right)\right|_{J=0}
= \left.\exp\left(\frac{1}{2}\partial_x^\mathrm{T} A^{-1}\partial_x\right) x_{k_1}\cdots x_{k_{2N}}\right|_{x=0}.$$
The left hand side is $\langle x_{k_1}\cdots x_{k_{2N}}\rangle$, by previous arguments.
But
$$\begin{eqnarray}
\left.\exp\left(\frac{1}{2}\partial_x^\mathrm{T} A^{-1}\partial_x\right) x_{k_1}\cdots x_{k_{2N}}\right|_{x=0} &=&
\frac{1}{2^N N!} \left(\partial_x^\mathrm{T} A^{-1}\partial_x\right)^N x_{k_1}\cdots x_{k_{2N}} \\
&=& \frac{1}{2^N N!} \sum_{\sigma \in S_{2N}}(A^{-1})_{k_{\sigma(1)}k_{\sigma(2)}} \cdots (A^{-1})_{k_{\sigma(2N-1)}k_{\sigma(2N)}} \\
&=& \langle x_{k_1}\cdots x_{k_{2N}}\rangle.
\end{eqnarray}$$
For example, for $N=1$,
$$\begin{eqnarray}
\left.\exp\left(\frac{1}{2}\partial_x^\mathrm{T} A^{-1}\partial_x\right) x_{k_1} x_{k_{2}}\right|_{x=0}
&=& \left(\frac{1}{2} \partial_x^\mathrm{T} A^{-1}\partial_x\right) x_{k_1} x_{k_2} \\
&=& \frac{1}{2} A^{-1}_{ij}(\delta_{i k_1}\delta_{j k_2} + \delta_{i k_2}\delta_{j k_1}) \\
&=& \frac{1}{2} (A^{-1}_{k_1 k_2} + A^{-1}_{k_2 k_1}) \\
&=& A^{-1}_{k_1 k_2}.
\end{eqnarray}$$

References

Many texts on quantum field theory deal with such finite dimensional integrals on the way to treating the infinite dimensional case. See, for example,

A. Zee. Quantum field theory in a nutshell

J. Zinn-Justin. Quantum Field Theory and Critical Phenomena