Stable computation of entanglement entropy for 2D interacting fermion systems

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.