  In the absence of a fermion sign problem, auxiliary field (or determinantal)
quantum Monte Carlo (DQMC) approaches have long been the numerical method of
choice for unbiased, large-scale simulations of interacting many-fermion
systems. More recently, the conceptual scope of this approach has been expanded
by introducing ingenious schemes to compute entanglement entropies within its
framework. On a practical level, these approaches however suffer from a variety
of numerical instabilities that have largely impeded their applicability. Here
we report on a number of algorithmic advances to overcome many of these
numerical instabilities and significantly improve the calculation of
entanglement measures in the zero-temperature projective DQMC approach,
ultimately allowing to reach similar system sizes as for the computation of
conventional observables. We demonstrate the applicability of this improved
DQMC approach by providing an entanglement perspective on the quantum phase
transition from a magnetically ordered Mott insulator to a band insulator in
the bilayer square lattice Hubbard model at half filling.
