Stable computation of entanglement entropy for two-dimensional interacting fermion systems

There is no doubt that the information hidden in entanglement entropy (EE), for example, the nth order R\'enyi EE, i. e. , S₍^A=11-nlnTr (₀^n), where ₀=Tr₀ is the reduced density matrix, can be used to infer the organizing principle of two-dimensional (2D) interacting fermion systems, ranging from spontaneous symmetry-breaking phases and quantum critical points to topologically ordered states. It is far from clear, however, whether EE can be obtained with the precision required to observe these fundamental features---usually in the form of universal finite-sized scaling behavior. Even for the prototypical 2D interacting fermion model---the Hubbard model---to all existing numerical algorithms, the computation of EE has not been successful with reliable data from which the universal scaling regime can be accessed. Here, we explain the reason for these unsuccessful attempts of 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 EE in 2D interacting fermion systems can be achieved and universal scaling information can be extracted. Relevance toward experimental 2D interacting fermion systems is discussed.