= Solution
The probabilistic form of <Wick's theorem> is <Isserlis's theorem>. If $(Y_1,\ldots,Y_k)$ is a centered <random vector> with a <multivariate Gaussian distribution> and <covariance matrix> $C$, its odd product <moments> vanish, and, for even $k$,
$$
\boxed{\mathbb E\prod_{a=1}^kY_a
=\sum_{\pi\in\mathcal P_2(k)}\prod_{\{a,b\}\in\pi}C_{ab}.}
$$
Here $\mathcal P_2(k)$ is the set of <perfect matchings> of the $k$ labelled positions. The formula allows repeated coordinates and singular <covariance matrices>. To see the pairing rule, differentiate the <moment-generating function> $\exp(t^TCt/2)$ once in each coordinate and set $t=0$: only terms pairing all differentiated coordinates survive, each with the displayed product of <covariances>.
Take $Y_a=M_{i_ai_{a+1}}$, with $i_{k+1}=i_1$. These variables are linear combinations of the original <independent> <normal random variables>, hence are jointly <Gaussian>. They have zero <means>, so
$$
\boxed{\mathbb E[M_{i_1i_2}\cdots M_{i_ki_1}]
=\sum_{\pi\in\mathcal P_2(k)}\prod_{\{a,b\}\in\pi}
\left(\delta_{i_ai_b}\delta_{i_{a+1}i_{b+1}}
+\delta_{i_ai_{b+1}}\delta_{i_{a+1}i_b}\right).}
$$
Substituting the entry <covariance> from part (a)(i) identifies each parenthesis with $\mathbb E[M_{i_ai_{a+1}}M_{i_bi_{b+1}}]$, giving the requested identity.
Back to article page
