  We consider an efficient scheme to simulate fermionic Hubbard models with
nonlocal density-density interactions in two dimensions, based on bond-centered
auxiliary-field quantum Monte Carlo. The simulations are shown to be
sign-problem free within a finite, restricted parameter range. Using this
approach, we first study the Hubbard model on the half-filled square lattice
bilayer, including an interlayer repulsion term in addition to the local
repulsion, and present the ground state phase diagram within the accessible
parameter region. Starting from the antiferromagnetically ordered state in the
absence of interlayer repulsion, the interlayer interactions are found to
destabilize the antiferromagnetic order, leading to a band insulator state.
Moreover, for sufficiently strong interlayer tunneling, we also observe the
emergence of a direct dimer product state of mixed D-Mott and S-Mott character
along the equal coupling line. We discuss the stability range of this state
within strong-coupling perturbation theory. Furthermore, we consider the
Hubbard model on the honeycomb lattice with next-nearest-neighbor interactions.
Such an interaction is found to enhance both charge density and spin-current
correlations within the semimetallic region. However, inside the accessible
parameter region, they do not stabilize long-ranged charge density wave order
nor a quantum spin Hall state, and the only insulating state that we observe
exhibits long-range antiferromagnetism.
