  Given a specific interacting quantum Hamiltonian in a general spatial
dimension, can one access its entanglement properties, such as, the
entanglement entropy corresponding to the ground state wavefunction? Even
though progress has been made in addressing this question for interacting
bosons and quantum spins, as yet there exist no corresponding methods for
interacting fermions. Here we show that the entanglement structure of
interacting fermionic Hamiltonians has a particularly simple form -- the
interacting reduced density matrix can be written as a sum of operators that
describe free fermions. This decomposition allows one to calculate the Renyi
entropies for Hamiltonians which can be simulated via Determinantal Quantum
Monte Carlo, while employing the efficient techniques hitherto available only
for free fermion systems. The method presented works for the ground state, as
well as for the thermally averaged reduced density matrix.
