  There is no doubt that the information hidden in entanglement entropy (EE),
for example, the $n$-th order R\'enyi EE, i.e., $S^{A}_n=\frac{1}{1-n}\ln \Tr
(\rho_A^n)$ where $\rho_A=\mathrm{Tr}_{\overline{A}}\rho$ is the reduced
density matrix, can be used to infer the organizing principle of 2D interacting
fermion systems, ranging from spontaneous symmetry breaking phases, quantum
critical points to topologically ordered states. It is far from clear, however,
whether the EE can actually be obtained with the precision required to observe
these fundamental features -- usually in the form of universal finite size
scaling behavior. Even for the prototypical 2D interacting fermion model -- the
Hubbard model, to all existing numerical algorithms, the computation of the EE
has not been succeeded with reliable data that the universal scaling regime can
be accessed. Here we explain the reason for these unsuccessful attempts in EE
computations in quantum Monte Carlo simulations in the past decades and more
importantly, show how to overcome the conceptual and computational barrier with
the incremental algorithm, such that the stable computation of the EE in 2D
interacting fermion systems can be achieved and universal scaling information
can be extracted. Relevance towards the experimental 2D interacting fermion
systems is discussed.
