Many-body-localization protection of eigenstate topological order in two dimensions

Many-body localization (MBL) has been proposed to enable and protect topological order in all eigenstates, vastly expanding the traditional ground-state setting. However, for the most intriguing case of two-dimensional (2D) systems with anyons and topology-dependent degeneracies, the dense many-body spectrum challenges studying this MBL protection numerically. Here we use large-scale full-spectrum variational ans\"atze to demonstrate MBL-protected topological order in the disordered 2D toric code perturbed by magnetic fields. We show that the system has topological local integrals of motion (tLIOMs) for magnetic field strengths below $h_c\approx0.1$ times the toric code coupling scale. Combining tLIOMs with exact diagonalization, we also identify high-energy topological multiplets in the dense many-body spectrum. The phase diagram we find is consistent with toric-code and trivial MBL phases being separated by an intervening thermal phase.