Wick's theorem for Gaussian stochastic variables

Question

I wonder whether there exists a clever way to implement Wick's theorem for Gaussian stochastic variables $\eta_{j_{i}}$ (with $\langle \eta_{j_{i}}\rangle=0$ for $\forall i$) which in general states:

$$ \begin{align} \langle \eta_{j_1}\eta_{j_2}\dots\eta_{j_{2n+1}}\rangle&=0\\ \langle \eta_{j_1}\eta_{j_2}\dots\eta_{j_{2n}}\rangle&=\sum_\mathrm{P_d}\sigma_{k_1, k_2}\sigma_{k_3, k_4}\dots\sigma_{2n-1, 2n}\ \end{align} $$ where one has to sum over only those $(2n)!/(2^n n!)$ permutations $(j_1\dots j_{2n})\Rightarrow (k_1,\dots,k_{2n})$ which lead to different expressions for $\sigma_{k_1, k_2}\dots\sigma_{k_{2n-1}, k_{2n}}$. And $\sigma_{jk}=\langle \eta_j\eta_k\rangle$ is defined as the covariance matrix.

I think computationally I would prefer a functional form where I define a $\langle \dots\rangle$ operator.

Problems that occur to my mind:

• Need to distinguish between random variables and deterministic variables (the latter ones can be pulled out of the operator).

• How do I cover all the cases at the same time, i.e. four-point correlation functions, six-point correlation functions etc.

• How do I impose the rules to make sure that only permutations which lead to different expressions remain?

Maybe anyone of you has already tried to implement this theorem/problem? Any help is highly appreciated as always :)

Answer 1

First define a function that tells you whether a variable is random or not

constantQ[x_?NumericQ]:=False

And, according to your program you might want to give some up values as follows

x1/:constantQ[x1]:=True;
y1/:constantQ[y1]:=False;
(*etc...*)

Then define a function ExpValue that encloses all your variables. Make it so that variables that evaluate to True when constantQ is applied are taken out

ExpValue[a___,x_?constantQ,b___]:= x ExpValue[a,b];

And apply the actual Wick contraction rule when all variables have been taken out

ExpValue[a__] /; ! (constantQ /@ And @@ {a}) := Wcontract[{a}];

Now the function Wcontract can be defined recursively

Wcontract[{}]:=1;
Wcontract[a_List] := Sum[sigma[a[[1]], a[[j]]] Wcontract[Delete[a, {{1}, {j}}]], {j, 2, Length[a]}]

Where I called sigma the two point correlator.