Assume we are given are fermionic many-particle system with creation operators $c_1^\dagger,c_2^\dagger,…$ acting on some Hilbert space $\mathcal{H}$. That is, the $c_i^\dagger$ and their corresponding annihilation operators $c_i=(c_i^\dagger)^\dagger$ act as bounded operators on $\mathcal{H}$ and satisfy the CAR relations $$\{c_i,c_j\}=\{c^\dagger_i,c_j^\dagger\}=0\quad\text{and}\quad\{c_i^\dagger,c_j\}=\delta_{ij}\quad\forall i,j$$
In the bosonic case, the Wick theorem can be elegantly written in the form
$$\prod_{j=1}^nc_{i_j}^{\sigma_j}=\sum_{Q\subset\mathbb{N}_n}\langle\prod_{j\in Q}^nc_{i_j}^{\sigma_j}\rangle:\prod_{j\not\in Q}^nc_{i_j}^{\sigma_j}:$$
for any $i_1,…,i_n\in\mathbb{N}$ and $\sigma_1,…,\sigma_n\in\{-1,1\}$, where $c^-_i:=c_i$, $c^+_i:=c^\dagger_i$ for any $i\in\mathbb{N}$, $\mathbb{N}_n:=\{1,\dots,n\}$, $\langle A\rangle:=\langle\Omega\mid A\Omega\rangle$ tenotes the vacuum expectation value and $:\cdots:$ denotes the Normal ordering. However, for fermions this formula has to be modified by attaching an appropriate sign to every of the summands. What is this sign?
Best Answer
Here's a trick to recover the fermionic signs, knowing the bosonic formulas. First, introduce a sufficient number of algebraically independent "c-numbers" $\epsilon_k$, that are only required to anti-commute among themselves and with the $c_i$, $c_i^\dagger$ generators. Then, according to the usual sign rules, the products $\epsilon c_i$ and $\epsilon c_j^\dagger$, where $\epsilon$ is any one of the new generators, will obey the bosonic canonical commutation relations with extra anti-commuting factors: $$ [\epsilon c_i, \eta c_j^\dagger] = -\epsilon \eta \{c_i, c_j^\dagger\} = -\epsilon \eta \delta_{ij} = -\epsilon \eta \langle c_i c_j^\dagger \rangle = \langle (\epsilon c_i) (\eta c_j^\dagger) \rangle, $$ for example. The extra intermediate $-$ sign appeared because the $c_i$ and $\eta$ exchanged positions while moving $\eta$ out of the commutator or moving it inside the expectation value. Note that $\epsilon$ and $\eta$ need to be chosen independently, to avoid writing down tautologies like $0=0$ due to the nilpotency $\epsilon^2 = 0$ of the new generators.
So, the bosonic Wick formula still works in the following way $$ \prod_{j=1}^n \epsilon_{k_j} c_{i_j}^{\sigma_j}=\sum_{Q\subset\mathbb{N}_n}\langle\prod_{j\in Q}^n \epsilon_{k_j} c_{i_j}^{\sigma_j}\rangle:\prod_{j\not\in Q}^n \epsilon_{k_j} c_{i_j}^{\sigma_j}: , $$ where the $\epsilon_{k_j}$ multipliers are all distinct for distinct $j$. I'll leave it as an exercise to obtain the exact sign factor that you were looking for, which you can get by factoring out $\prod_{j=1}^n \epsilon_{k_j}$ from both sides of the above equation.