Proof that the reduced density matrix of free fermions is thermal? I found this question here but it was partly unanswered. The question remains, namely:
Given a free theory of fermions in a bi-partite system $S=A\cup B$ with Hamiltonian
$$
H = \sum_{ij} t_{ij}a^{\dagger}_ia_j\quad \longrightarrow \quad H = \sum_k E_k c^{\dagger}_kc_k
$$
can anybody help me with proving that given any eigenstate $|\psi\rangle$ of $H$ with its associated density matrix $\rho=|\psi\rangle\langle\psi|$ of $H$, the reduced density matrix $\rho_B=\text{tr}_A(\rho)$ has a "thermal" form:
$$
\rho_B \sim \exp(-H_B) \quad \text{with} \quad H_B=\sum_i h_{ij}c^{\dagger}_ic_j.
$$

The proof I’m reading sketches that since $|\psi\rangle$ is a determinant (?), then its correlators factorise
$$
\langle c_i^{\dagger}c_j^{\dagger}c_kc_l\rangle_{\psi}= \langle c_n^{\dagger}c_l\rangle_{\psi} \langle c_m^{\dagger}c_k\rangle_{\psi} − \langle c^{\dagger}_nc_k \rangle_{\psi}\langle c^{\dagger}_mc_l\rangle _{\psi}
$$
($\langle \cdot \rangle_{\psi} := \langle \psi | \cdot | \psi \rangle$), therefore
$$
C_{ij} := \text{Tr}(\rho_B c_i^{\dagger} c_j)
$$
must factorise as well. According to Wick’s theorem, this property holds if (why? is this if an if and only if?) $\rho$ is the exponential of a free-fermion operator
$$
\rho_B = K\exp(-H_B) \quad \text{with} \quad H_B=\sum_i h_{ij}c^{\dagger}_ic_j.
$$
My problem is that this might not be the unique form of $\rho_B$. The above to me is a bit meaningless without the if and only if.
 A: The answer assumes that you know that for Gaussian states, i.e. states of the form $\exp(-H_B)$, with $H_B$ quadratric in the fermionic operators, Wick's theorem holds, that is, $2n$-point correlators $\langle f_1\cdots f_{2n}\rangle$ can be expressed in terms of two-point correlators, where the $f_i$ are some fermionic operators (creation or annihiliation).
The question remains why this implies that a state which satisfies Wick's theorem must be a Gaussian state. The reason is that any quantum state is entirely determined by all its expectation values $\langle O\rangle$. For fermionic states, all admissible operators must have even parity, and thus, knowing all $\langle f_1\cdots f_{2n}\rangle$ means that you know all $\langle O\rangle$, and thus you know everything about the state: That is, the state is uniquely determined by all $\langle f_1\cdots f_{2n}\rangle$.
So if you have a state which satisfies Wick's theorem, you know it must be Gaussian, since there is a Gaussian state with the same $\langle f_1\cdots f_{2n}\rangle$, and those uniquely determine the state.
