  Many-body localization is a striking mechanism that prevents interacting
quantum systems from thermalizing. The absence of thermalization behaviour
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 iPEPS 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 non-interacting
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.
