Not yet reviewed by Pith; the record is open.

The nature of correlated states in twisted bilayer graphene (TBG) at the magic angle has received intense attention in recent years. We present a numerical study of an interacting Bistritzer-MacDonald (IBM) model of TBG using a suite of methods in quantum chemistry, including Hartree-Fock, coupled cluster singles, doubles (CCSD), and perturbative triples (CCSD(T)), as well as a quantum chemistry formulation of the density matrix renormalization group method (DMRG). Our treatment of TBG is agnostic to gauge choices, and hence we present a new gauge-invariant formulation to detect the spontaneous symmetry breaking in interacting models. To benchmark our approach, we focus on a simplified spinless, valleyless IBM model. At integer filling ($\nu=0$), all numerical methods agree in terms of energy and $C_{2z} \mathcal{T}$ symmetry breaking. Additionally, as part of our benchmarking, we explore the impact of different schemes for removing ``double-counting'' in the IBM model. Our results at integer filling suggest that cross-validation of different IBM models may be needed for future studies of the TBG system. After benchmarking our approach at integer filling, we perform the first systematic study of the IBM model near integer filling (for $|\nu|< 0.2$). In this regime, we find that the ground state can be in a metallic and $C_{2z} \mathcal{T}$ symmetry breaking phase. The ground state appears to have low entropy, and therefore can be relatively well approximated by a single Slater determinant. Furthermore, we observe many low entropy states with energies very close to the ground state energy in the near integer filling regime.