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.