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 :)


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


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


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[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.

| improve this answer | |

Package-X has an internal function X`Utilities`PassarinoVeltmanTensor that generates the polynomial for you. The code below generates the correlation $\langle\eta_{j1}\eta_{j2}\eta_{j3}\eta_{j4}\rangle$:


    /. LTensor -> (corr[#2, #3] &)
corr[j1, j4] corr[j2, j3] + 
     corr[j1, j3] corr[j2, j4] + corr[j1, j2] corr[j3, j4]

You can increase the number of indices to calculate bigger correlations.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.