Many-body localization is a striking mechanism that prevents interacting quantum systems from thermalizing. The absence of thermalization behavior manifests itself, for example, in a remanence of local particle number configurations, a quantity that is robust over a parameter range. Local particle numbers are directly accessible in programmable quantum simulators, in systems of cold atoms, even in two spatial dimensions. Yet, the classical simulation aimed at building trust in quantum simulations is highly challenging. In this work, we present a comprehensive tensor network simulation of a many-body localized systems in two spatial dimensions using a variant of an infinite projected entangled pair states algorithm. The required translational invariance can be restored by implementing the disorder into an auxiliary spin system, providing an exact disorder average under dynamics. We can quantitatively assess signatures of many-body localization for the infinite system: Our methods are powerful enough to provide crude dynamical estimates for the transition between localized and ergodic phases. Interestingly, in this setting of finitely many disorder values, which we also compare with simulations involving noninteracting fermions and for which we discuss the emergent physics, localization emerges in the interacting regime, for which we provide an intuitive argument, while Anderson localization is absent.
