Proof of this version of Wick’s theorem? I can’t prove/find the following version of Wick’s theorem:
Say we have a system of free fermions with Hamiltonian
$$
H = \sum_{ij} t_{ij}c^{\dagger}_ic_j\quad \longrightarrow \quad H = \sum_k E_k d^{\dagger}_kd_k.
$$
An eigenstate  $|\psi \rangle$ of $H$ is given by acting on $|0\rangle$ with $d^{\dagger}$s as usual:
$$
|\psi^{N_p}\rangle = \prod_{a \in N_p} d_a | 0 \rangle \quad \text{where} \quad d_a = \sum_i\phi^a_ic_i,
$$
for some suitable $\phi$s.
How can one prove that $2n$-point correlators factorise into products of 2-point correlators? I.e.
$$
\langle c_i^{\dagger}c_j^{\dagger}c_kc_l\rangle_{\psi}\equiv \langle \psi | c_i^{\dagger}c_j^{\dagger}c_kc_l|\psi\rangle =\langle c_i^{\dagger}c_l\rangle_{\psi} \langle c_j^{\dagger}c_k\rangle_{\psi} − \langle c^{\dagger}_ic_k \rangle_{\psi}\langle c^{\dagger}_jc_l\rangle _{\psi},
$$
for any $\psi$.
Edit: Found above statement at https://arxiv.org/abs/cond-mat/0212631 . "Because $|\psi \rangle$ is a determinant, all higher order correlations can be written in terms of $C$"[where $C_{ij}=\langle c^{\dagger}_i c_j\rangle$].
 A: On a product states such as $\lvert \psi \rangle = \prod \limits_{a \in N_p} d^\dagger_a \lvert 0 \rangle$, for $\hat n_i := d^\dagger_i d_i$, and $\langle \hat n_i \rangle = n_i$, $\langle \hat n_i \hat n_j \rangle = n_i n_j = \langle \hat n_i \rangle \langle \hat n_j \rangle$ since $n_i \in \{0, 1\}$. Also, $\delta_{ij} \langle n_i \rangle = \langle d^\dagger_i d_j \rangle$ (I will use this later).
Now consider the matrix element $\langle d^\dagger_i d^\dagger_j d_k d_l \rangle$. Off-diagonal elements vanish, so the only non-zero terms are when $i=k, j=l$ and $i=l, j=k$. Therefore
$$\begin{aligned}
\langle d^\dagger_i d^\dagger_j d_k d_l \rangle 
& = \delta_{ik} \delta_{jl} \langle d^\dagger_i d^\dagger_j d_k d_l \rangle
+
\delta_{il} \delta_{jk} \langle d^\dagger_i d^\dagger_j d_k d_l \rangle \\ 
&= 
\delta_{ik} \delta_{jl} \langle d^\dagger_i (-d_k d^\dagger_j + \delta_{jk}) d_l \rangle 
-
\delta_{il} \delta_{jk} \langle d^\dagger_j d^\dagger_i d_k d_l \rangle \\
&= 
\delta_{ik} \delta_{jl} \delta_{jk} \langle d^\dagger_i d_l \rangle 
-
\delta_{ik} \delta_{jl} \langle d^\dagger_i d_k d^\dagger_j d_l \rangle 
-
\delta_{il} \delta_{jk} \langle d^\dagger_j (- d_k d^\dagger_i + \delta_{ik}) d_l \rangle \\
&= 
\delta_{ik} \delta_{jl} \delta_{jk} \delta_{il} \langle d^\dagger_i d_l \rangle 
-
\delta_{ik} \delta_{jl} \langle d^\dagger_i d_k d^\dagger_j d_l \rangle 
+
\delta_{il} \delta_{jk} \langle d^\dagger_j d_k d^\dagger_i d_l \rangle
-
\delta_{il} \delta_{jk} \delta_{ik} \delta_{jl} \langle d^\dagger_j d_l \rangle,
\end{aligned}$$
where I added redundant kronecker deltas in the last line to show that the first and last term cancel each other. Now we can note that $\delta_{ik} \delta_{jl} \langle d^\dagger_i d_k d^\dagger_j d_l \rangle = \delta_{ik} \delta_{jl} \langle \hat n_i \hat n_j \rangle = \langle d^\dagger_i d_k \rangle \langle d^\dagger_j d_l \rangle$, and similarly for the third term as well. So we have
$$ \langle d^\dagger_i d^\dagger_j d_k d_l \rangle = \langle d^\dagger_j d_k \rangle \langle d^\dagger_i d_l \rangle - \langle d^\dagger_i d_k \rangle \langle d^\dagger_j d_l \rangle.$$
We can substitute $c^\dagger_i = \sum_{a} \phi^a_i d^\dagger_a$, which transfers the results to the original basis:
$$\begin{aligned}
\langle c^\dagger_i c^\dagger_j c_k c_l\rangle 
&=
\sum_{abcd} \phi^a_i \phi^b_j \phi^c_k \phi^d_l \langle d^\dagger_a d^\dagger_b d_c d_d \rangle \\ & =  \sum_{abcd} \phi^a_i \phi^b_j \phi^c_k \phi^d_l \left(\langle d^\dagger_a d_d \rangle \langle d^\dagger_b d_c \rangle - \langle d^\dagger_a d_c \rangle \langle d^\dagger_b d_d \rangle \right) \\
&=
\left\langle \left(\sum_a \phi^a_i d^\dagger_a \right) \left(\sum_d \phi^d_l d_d \right) \right\rangle \left\langle \left(\sum_b \phi^b_j d^\dagger_b \right) \left(\sum_c \phi^c_k d_c \right) \right\rangle
-
\left\langle \left(\sum_a \phi^a_i d^\dagger_a \right) \left(\sum_c \phi^c_k d_c \right) \right\rangle \left\langle \left(\sum_b \phi^b_j d^\dagger_b \right) \left(\sum_d \phi^d_l d_d \right) \right\rangle \\
&=
\langle c^\dagger_i c_l \rangle \langle c^\dagger_j c_k \rangle - \langle c^\dagger_i c_k \rangle \langle c^\dagger_j c_l \rangle,
\end{aligned}$$
as required. Note that this holds for any other basis related via a unitary transformation, since the order of the operators do not change.
Edit: I forgot to consider the case where $i=j=k=l$, but in that case, the matrix element vanishes, and in cases where at least one of the indices are different, the above argument holds.
A: I think the easiest way to understand it is to reduce to "conventional" Wick theorem (I assume that the proof of this one is known). Any state from fermion Fock space basis $\prod \limits_{a \in N_p} d^\dagger_a |0\rangle$ can be thought of as a vacuum for system of "particles" (with creation operators $d^\dagger_i, i \notin N_p$) and "holes" (with creation operators $a^\dagger_i \equiv d_i, i \in N_p$). There is a symmetry between $d$ and $d^\dagger$ in fermionic algebra, of which there is, of course, no analogue in the bosonic case.
