Transition to a many-body localized regime in a two-dimensional disordered quantum dimer model

TL;DR

This paper provides evidence for a transition to a many-body localized phase in a two-dimensional disordered quantum dimer model, advancing understanding of localization phenomena beyond one dimension through large-scale simulations.

Contribution

It demonstrates the existence of a many-body localization transition in a 2D quantum dimer model, a significant step beyond prior 1D studies, using extensive numerical analysis.

Findings

Evidence of a transition to a localized regime in 2D

Large-scale simulations of static and dynamical properties

Highlighting unusual quantum dynamics in constrained systems

Abstract

Many-body localization is a unique physical phenomenon driven by interactions and disorder for which a quantum system can evade thermalization. While the existence of a many-body localized phase is now well-established in one-dimensional systems, its fate in higher dimension is an open question. We present evidence for the occurrence of a transition to a many-body localized regime in a two-dimensional quantum dimer model with interactions and disorder. Our analysis is based on the results of large-scale simulations for static and dynamical properties of a consequent number of observables. Our results pave the way for a generic understanding of occurrence of a many-body localization transition in dimension larger than one, and highlight the unusual quantum dynamics that can be present in constrained systems.

Figures

Tables

Table 1. Table 1: Two-dimensional lattices used in this work, with their shape (regular or tilted square lattices, or rectangular) and the size of their Hilbert space in the ( W x = 0 , W y = 0 ) formulae-sequence subscript 𝑊 𝑥 0 subscript 𝑊 𝑦 0 (W_{x}=0,W_{y}=0) winding sector.

Lattice	$4\times 6$	$32$	$4\times 8$	$6\times 6$	$40$	$4\times 10$	$6\times 8$	$52$	$8\times 8$
Lattice type	Rectangular	Tilted sq.	Rectangular	Regular sq.	Tilted sq.	Rectangular	Rectangular	Tilted sq.	Regular sq.
Hilbert space size $\mathscr{H}$	$\mathrm{1\hspace{0.17em}456}$	$\mathrm{13\hspace{0.17em}348}$	$\mathrm{17\hspace{0.17em}412}$	$\mathrm{44\hspace{0.17em}176}$	$\mathrm{139\hspace{0.17em}212}$	$\mathrm{216\hspace{0.17em}016}$	$\mathrm{1\hspace{0.17em}504\hspace{0.17em}896}$	$\mathrm{4\hspace{0.17em}572\hspace{0.17em}468}$	$\mathrm{153\hspace{0.17em}722\hspace{0.17em}916}$

Equations

$$H$$

$$\mathcal{I}=\frac{1}{N_p}\sum _{{\square }_p}{e}^{\mathrm{i}{\phi }_p}\left(2\left(|\text{leavevmode}to11.55pt{\text{vbox}to12.75pt\text{pgfpicture}\text{makeatletter}\text{lower}-1.39497ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}0.0pt9.95863pt\text{pgfsys@lineto}9.95863pt9.95863pt\text{pgfsys@lineto}9.95863pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt9.95863pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}4.97931pt9.95863pt\text{pgfsys@moveto}10.5561pt9.95863pt\text{pgfsys@curveto}10.5561pt10.6186pt8.05933pt11.15361pt4.97931pt11.15361pt\text{pgfsys@curveto}1.89929pt11.15361pt-0.59749pt10.6186pt-0.59749pt9.95863pt\text{pgfsys@curveto}-0.59749pt9.29866pt1.89929pt8.76366pt4.97931pt8.76366pt\text{pgfsys@curveto}8.05933pt8.76366pt10.5561pt9.29866pt10.5561pt9.95863pt\text{pgfsys@closepath}\text{pgfsys@moveto}4.97931pt9.95863pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}4.97931pt0.0pt\text{pgfsys@moveto}10.5561pt0.0pt\text{pgfsys@curveto}10.5561pt0.65997pt8.05933pt1.19498pt4.97931pt1.19498pt\text{pgfsys@curveto}1.89929pt1.19498pt-0.59749pt0.65997pt-0.59749pt0.0pt\text{pgfsys@curveto}-0.59749pt-0.65997pt1.89929pt-1.19498pt4.97931pt-1.19498pt\text{pgfsys@curveto}8.05933pt-1.19498pt10.5561pt-0.65997pt10.5561pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}4.97931pt0.0pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}_p⟩⟨\text{leavevmode}to11.55pt{\text{vbox}to12.75pt\text{pgfpicture}\text{makeatletter}\text{lower}-1.39497ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}0.0pt9.95863pt\text{pgfsys@lineto}9.95863pt9.95863pt\text{pgfsys@lineto}9.95863pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt9.95863pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}4.97931pt9.95863pt\text{pgfsys@moveto}10.5561pt9.95863pt\text{pgfsys@curveto}10.5561pt10.6186pt8.05933pt11.15361pt4.97931pt11.15361pt\text{pgfsys@curveto}1.89929pt11.15361pt-0.59749pt10.6186pt-0.59749pt9.95863pt\text{pgfsys@curveto}-0.59749pt9.29866pt1.89929pt8.76366pt4.97931pt8.76366pt\text{pgfsys@curveto}8.05933pt8.76366pt10.5561pt9.29866pt10.5561pt9.95863pt\text{pgfsys@closepath}\text{pgfsys@moveto}4.97931pt9.95863pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}4.97931pt0.0pt\text{pgfsys@moveto}10.5561pt0.0pt\text{pgfsys@curveto}10.5561pt0.65997pt8.05933pt1.19498pt4.97931pt1.19498pt\text{pgfsys@curveto}1.89929pt1.19498pt-0.59749pt0.65997pt-0.59749pt0.0pt\text{pgfsys@curveto}-0.59749pt-0.65997pt1.89929pt-1.19498pt4.97931pt-1.19498pt\text{pgfsys@curveto}8.05933pt-1.19498pt10.5561pt-0.65997pt10.5561pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}4.97931pt0.0pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}_p|+|\text{leavevmode}to12.95pt{\text{vbox}to11.55pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.79749ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}0.0pt9.95863pt\text{pgfsys@lineto}9.95863pt9.95863pt\text{pgfsys@lineto}9.95863pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt9.95863pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt4.97931pt\text{pgfsys@moveto}1.29466pt4.97931pt\text{pgfsys@curveto}1.29466pt8.05933pt0.71503pt10.5561pt0.0pt10.5561pt\text{pgfsys@curveto}-0.71503pt10.5561pt-1.29466pt8.05933pt-1.29466pt4.97931pt\text{pgfsys@curveto}-1.29466pt1.89929pt-0.71503pt-0.59749pt0.0pt-0.59749pt\text{pgfsys@curveto}0.71503pt-0.59749pt1.29466pt1.89929pt1.29466pt4.97931pt\text{pgfsys@closepath}\text{pgfsys@moveto}0.0pt4.97931pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}9.95863pt4.97931pt\text{pgfsys@moveto}11.2533pt4.97931pt\text{pgfsys@curveto}11.2533pt8.05933pt10.67366pt10.5561pt9.95863pt10.5561pt\text{pgfsys@curveto}9.2436pt10.5561pt8.66397pt8.05933pt8.66397pt4.97931pt\text{pgfsys@curveto}8.66397pt1.89929pt9.2436pt-0.59749pt9.95863pt-0.59749pt\text{pgfsys@curveto}10.67366pt-0.59749pt11.2533pt1.89929pt11.2533pt4.97931pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt4.97931pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}_p⟩⟨\text{leavevmode}to12.95pt{\text{vbox}to11.55pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.79749ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}0.0pt9.95863pt\text{pgfsys@lineto}9.95863pt9.95863pt\text{pgfsys@lineto}9.95863pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt9.95863pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt4.97931pt\text{pgfsys@moveto}1.29466pt4.97931pt\text{pgfsys@curveto}1.29466pt8.05933pt0.71503pt10.5561pt0.0pt10.5561pt\text{pgfsys@curveto}-0.71503pt10.5561pt-1.29466pt8.05933pt-1.29466pt4.97931pt\text{pgfsys@curveto}-1.29466pt1.89929pt-0.71503pt-0.59749pt0.0pt-0.59749pt\text{pgfsys@curveto}0.71503pt-0.59749pt1.29466pt1.89929pt1.29466pt4.97931pt\text{pgfsys@closepath}\text{pgfsys@moveto}0.0pt4.97931pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}[named]pgffillcolorrgb\mathrm{1,.5,0}\text{pgfsys@color@rgb@fill}1.50\text{pgfsys@invoke}\text{pgfsys@moveto}9.95863pt4.97931pt\text{pgfsys@moveto}11.2533pt4.97931pt\text{pgfsys@curveto}11.2533pt8.05933pt10.67366pt10.5561pt9.95863pt10.5561pt\text{pgfsys@curveto}9.2436pt10.5561pt8.66397pt8.05933pt8.66397pt4.97931pt\text{pgfsys@curveto}8.66397pt1.89929pt9.2436pt-0.59749pt9.95863pt-0.59749pt\text{pgfsys@curveto}10.67366pt-0.59749pt11.2533pt1.89929pt11.2533pt4.97931pt\text{pgfsys@closepath}\text{pgfsys@moveto}9.95863pt4.97931pt\text{pgfsys@fillstroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}_p|\right)-\hat{I}\right).$$

$$S_p(|n⟩)=-\sum _i{p}_i\log (p_i)$$

$$S(|n⟩)=-{\mathrm{Tr}}_A{\rho }_A\log ({\rho }_A)$$

Full text

I Introduction

The emergence of thermalization in isolated quantum systems, which are not in contact with a bath and solely follow unitary Hamiltonian dynamics, has long been a central issue in statistical mechanics Neumann (1929); Goldstein et al. (2010); Gogolin and Eisert (2016). This question has been revived thanks to an extraordinary increase in the degree of control and isolation reached in cold-atom experiments Bloch et al. (2008); Gross and Bloch (2017). The central concept of Eigenstate Thermalization Hypothesis (ETH) Deutsch (1991); Srednicki (1994), a conjecture on expectation values of local operators in the eigenstates of the isolated system, allows to provide a physical understanding and justification of this thermalization. The ETH is found to be valid for the generic majority of many-body systems Rigol et al. (2008); Deutsch (2018). Recent years have witnessed a vast body of research concentrating on a situation where the ETH fails, namely the phenomenon of many-body localization (MBL) Basko et al. (2006); Gornyi et al. (2005), which occurs most noticeably in the presence of strong disorder.

In a many-body localized system, eigenstates at the same energy density do not form a statistical ensemble, and local expectation values strongly differ from one eigenstate to the other. MBL states sustain a low level of quantum entanglement in the form of an area law Bauer and Nayak (2013) (in contrast with ETH states which obey a volume law for entanglement entropy), and have been found to be able to sustain long-range order in cases where equilibrium states are forbidden to Huse et al. (2013); Bahri et al. (2015); Pekker et al. (2014); Kjäll et al. (2014); Chandran et al. (2014). Dynamics in MBL systems markedly differ from those of thermalizing systems, with e.g. long-time memory of the initial state, absence of transport for conserved quantities (localization), and a slow propagation of quantum information with logarithmic growth of entanglement Žnidarič et al. (2008); Bardarson et al. (2012). The latter example allows to distinguish MBL from Anderson localization Anderson (1958) (where entanglement is bounded) which occurs in non-interacting systems.

The existence of many-body localization in one dimension is now established and reasonably well-understood in one-dimensional (1d) quantum lattice systems with short-range interactions, thanks to a concerted effort of approaches, including numerical simulations Pal and Huse (2010); Luitz et al. (2015); Wahl et al. (2017), phenomenological renormalization Vosk et al. (2015); Potter et al. (2015); Dumitrescu et al. (2017); Thiery et al. (2018); Zhang and Yao (2018) approaches, rigorous results Imbrie (2016a, b) and cold-atom experiments Schreiber et al. (2015); Smith et al. (2016); Lukin et al. (2018). Crucial to this understanding is the provable Imbrie (2016a) existence of local integral of motions Huse et al. (2014); Serbyn et al. (2013a), denoted as *l-*bits (for localized bits): these emergent localized operators, which diagonalize the Hamiltonian of the system at strong disorder, explain most of the salient features of the MBL phase Imbrie et al. (2017). In the typical setting of lattice models, the MBL phase indeed emerges in the presence of strong enough disorder, and is separated from a ETH phase at low disorder by a many-body localization transition involving an exponential number of eigenstates Pal and Huse (2010). The reviews Nandkishore and Huse (2015); Altman and Vosk (2015); Abanin and Papić (2017); Parameswaran et al. (2017); Alet and Laflorencie (2018); Abanin et al. (2019) provide an introduction and recent insights on this alive field of research.

I.1 Many-body localization in higher dimensions

An outstanding question which remains open is the fate of many-body localization in dimension higher than one. A basic observation in this case is that localization is more difficult to achieve as there are more transport channels than in one dimension. The Anderson localization survives in two dimensions (2d) for all disorder strengths, but with a much higher localization length, and needs an increasingly stronger critical disorder in higher dimensions Evers and Mirlin (2008).

In the many-body case, analytical arguments were presented De Roeck and Huveneers (2017); De Roeck and Imbrie (2017); Potirniche et al. (2019) that suggest the intrinsic instability of a MBL phase in dimension higher than one. They are based on the inevitable existence of thermal ‘bubbles’ (created by e.g. local regions of space with low disorder), which grow in time and which eventually thermalize any MBL sample. The thermalization conveyed in this scenario is asymptotic, meaning that it may require an asymptotically large time before signs of localization (for instance memory of the initial state) disappear. In these works, the MBL region is represented in terms of the *l-*bits, the exponentially localized integral of motions which form a complete set of operators spanning the Hilbert space.

In spite of these arguments, recent remarkable numerical simulations Wahl et al. (2018); Thomson and Schiró (2018); De Tomasi et al. (2019); Doggen et al. (2020) analyze possible signatures of a MBL phase in two-dimensional lattice models of spins or bosons, in either static Wahl et al. (2018) or dynamical quantities Thomson and Schiró (2018); De Tomasi et al. (2019); Doggen et al. (2020). While results of Refs. Wahl et al., 2018; Thomson and Schiró, 2018; De Tomasi et al., 2019 are in favor of the existence of a true thermodynamic MBL phase at finite disorder strenght in 2d, Ref. Doggen et al., 2020 concluded opppositely (albeit evidencing slow dynamics even at moderate disorder strength). These results are very intriguing, but the methods used in these works are not unbiased: in Refs. Wahl et al., 2018; Thomson and Schiró, 2018; De Tomasi et al., 2019 an approximate unitary diagonalization transformation based on the lowly-entangled nature or localized nature of MBL eigenstates is used, which misses by construction the ETH phase (we note nevertheless that Ref. Wahl et al., 2018 presents evidence for the transition in 2d, but without a finite-size or finite-entanglement scaling analysis). The finite-time dynamical results of Ref. Doggen et al., 2020 also rely on a method which is entanglement limited. Other exact diagonalization studies in 2d Geraedts et al. (2016); van Nieuwenburg et al. (2019); Wiater and Zakrzewski (2018) are limited to small systems (lattices with at most 16 lattice sites), such that finite-size scaling is not possible. We also note a recent interesting investigation using exact numerics on a 2d model in the continuum compatible with MBL, albeit with a transition point shifting with the number of particles Krishna et al. (2019). A quantum Monte Carlo treatment of a Hamiltonian constructed from *l-*bits also argues for the existence of MBL in a 2d spin system Inglis and Pollet (2016). On the experimental side, recent advances in cold-atom experiments present clear dynamical signs of localization (in terms of memory of an initial imbalance of atoms) on the longest accessible time scale in two dimensional systems, either for bosonic atoms in random potential Choi et al. (2016) or fermionic gases in quasi-periodic potentials Bordia et al. (2017).

In this work, we present a comprehensive study of a two-dimensional strongly interacting disordered model, based on extensive numerical simulations on fairly large systems (square lattices with up to 64 sites). These simulations are exact, unbiased, and allow to study both ETH and MBL regimes and the transition between them. We study an important range of clusters of different sizes and shapes, allowing to follow the evolution of many ETH or MBL indicators with size. This is made possible thanks to two main specificities of our work: we push large-scale parallel numerical methods to their limits (see Sec. II.2), and we work on a constrained model (described in Sec. II.1) that presents specific features appropriate for localization and numerical studies (Sec. I.2).

I.2 Searching for MBL in 2d constrained models

Constraints exist or emerge in the description of a wide variety of physical systems, ranging from gauge theories, surface physics (adsorption of molecules Fowler and Rushbrooke (1937)), glass physics Fredrickson and Andersen (1984); Ritort and Sollich (2003); Chandler and Garrahan (2010), quantum or classical magnets Lacroix et al. (2011); Castelnovo et al. (2012) or atomic physics (such as with Rydberg blockade Bernien et al. (2017)). Theoretical descriptions are written in terms of models with adapted degrees of freedom, such as dimer, loop or ice models, which exhibit physical features not present in unconstrained counterparts. Consider for instance the dimer model (Fig. 1) that will be used in this work, where the degrees of freedom are dimers located on bonds of a given lattice (or graph). The constraint is given by the following rule: each site of a lattice must be touched by one and only one dimer. Already at the classical static level, this local rule enforces long-range, algebraic correlations, between dimers on bipartite lattices Fisher and Stephenson (1963); Huse et al. (2003). Taking again the dimer example, it is not possible to obtain another dimer configuration by simply moving one dimer: only constrained moves (the smallest of which is represented in Fig. 1) are possible. This can result in very complex dynamics (see Refs. Oakes et al., 2016; Das et al., 2005; Cépas, 2014 for classical examples).

Recent work highlighted the also unusual quantum dynamics exhibited by quantum constrained models, either in theory Chamon (2005); Lan et al. (2018); Turner et al. (2018); Feldmeier et al. (2019); Pancotti et al. (2020) or experiments Bernien et al. (2017), as well as for related lattice gauge theories Smith et al. (2017a); Brenes et al. (2018); Smith et al. (2017b, 2018); Karpov et al. (2020); Papaefstathiou et al. (2020); Halimeh and Hauke (2020a, b). On one hand, we would thus expect that constrained models are good natural candidates to exhibit localization due to their enhanced slow dynamics. On the other hand, the local constraint is such that these models are already in a strongly interacting limit (even without interactions encoded in a Hamiltonian), which would rather favor delocalization and thermalization. This competition provides strong motivation to investigate the possibility of MBL in quantum constrained models. The interplay between constraints, interactions and disorder has been explored in recent works on 1d models Chen et al. (2018); Ostmann et al. (2019) (for Ref. Chen et al., 2018, equivalent to a dimer model on a ladder), which constraints are relevant for the physics of cold-atomic chains in the regime of Rydberg blockade Bernien et al. (2017), resulting in rich phase diagrams, with e.g. different mechanisms to produce MBL Chen et al. (2018).

Another important aspect (crucial for the numerical approach taken in this work) worth mentioning is the fact that the constraints overall reduce the size of the phase space for allowed configurations. For instance, the size of the Hilbert space scales as $\sim 1.34^{N}$ for dimer coverings of square lattices with $N$ sites, instead of $2^{N}$ for spin 1/2 models (or $4^{N}$ if no constraints were present). This greatly helps numerics based on exact approaches, as this enlarges the set of available finite samples in 2d. Given that exact diagonalization simulations have been instrumental in the study and understanding of MBL, this is also an important positive feature in the search for MBL in quantum constrained models.

The plan of the paper is as follows. We first describe the model and the 2d lattices used in this work in Sec. II.1, along with a brief description of the exact numerical methods employed (Sec. II.2). We then describe in detail physical properties of the eigenstates of this model in Sec. III, including spectral and eigenstate statistics, expectation values of local observables as well as half-system entanglement entropy. Sec. IV is devoted to a machine-learning analysis of entanglement spectra. Sec. V then presents an analysis of dynamical properties after a quench. Finally, we summarize our results and present perspectives in Sec. VI.

II Model and numerical methods

II.1 Model

This work focuses on a quantum dimer model on a square lattice with random potential terms. Quantum dimer models have a long history in condensed matter as low-energy effective models. They have been first introduced by Rokhsar and Kivelson Rokhsar and Kivelson (1988) in the context of the Resonating Valence Bond theory of high-temperature superconductivity, an interest which has been revived with the use of QDMs to describe experiments on underdoped cuprates Punk et al. (2015). QDMs are also prominent representations of highly frustrated magnets where dimers represent nearest-neighbor spin singlets 111See R. Moessner, K. S. Raman, Chapter 17 in Ref. Lacroix et al., 2011 for a review. Dimer models are indeed often more tractable than frustrated spin models from which they can be systematically derived Schwandt et al. (2010), while capturing the essential physics. For instance, the existence of 2d quantum spin liquids exhibiting topological order was first demonstrated in a QDM on a non-bipartite lattice Moessner and Sondhi (2001).

Quantum dimer models are defined in a configuration space where only configurations which fulfill the dimer constraint are allowed (dimer coverings, see Fig. 1 (a)). Their Hamiltonian contains a potential term that attributes an energy to every dimer covering and a kinetic term that flips dimers along a loop, allowing to pass from one dimer configuration to another. The specific QDM Hamiltonian that we consider in this work is based on the square lattice and contains a potential term with a random component, taking the following form

[TABLE]

The sums run over all plaquettes ${\square_{p}}$ of the square lattice, $V_{p}$ is a random potential different for each plaquette which is drawn uniformly from a box distribution $V_{p}\in[-V,V]$. The kinetic energy scale is set to $t=1$. The various terms are illustrated in Fig. 1 (b). At $V=0$ (no disorder), Ref. Lan and Powell, 2017 showed that eigenstates in the middle of the spectrum obey the eigenstate thermalization hypothesis. Quite importantly, there is no limit where Eq. 1 can be mapped to a model of free particles moving in a random potential (no Anderson localization limit).

Quantum dimer models based on coverings of bipartite lattices (such as the square lattice) have an important property: they display conserved topological numbers, known as winding numbers, when considering systems with periodic boundary conditions. The definition of the two conserved winding numbers (denoted $W_{x}$ and $W_{y}$) is given in Fig. 1 (c). It is easy to see that they are conserved by the Hamiltonian Eq. 1, and by any local moves (i.e. that do not go through periodic boundary conditions), hence their topological nature. $W_{x}$ and $W_{y}$ both range from $-L/2$ to $L/2$ on a square lattice. The Hamiltonian is block-diagonal, with sectors identified by $(W_{x},W_{y})$. We note that the resulting abelian $U(1)$ symmetry of the model Eq. 1 is compatible with many-body localization Potter and Vasseur (2016). Without loss of generality, we will restrict to the largest sector $(W_{x}=0,W_{y}=0)$ in what follows.

II.2 Lattices and numerical methods

This work uses two-dimensional lattices (represented in Fig. 14 in the Appendix) with $N$ sites and $2N$ bonds of the following form: regular square lattices with $N=L^{2}$ sites (with $L=4,6,8$), rectangular lattices with $N=L_{x}\times L_{y}$ sites (with $L_{x,y}=4,6,8,10$) as well as tilted square lattices with $N=32,40$ and $52$ sites (see Appendix for details). Periodic boundary conditions are taken. All the lattices used in this work are presented in Tab. 1, together with the corresponding Hilbert space sizes in the $(W_{x}=0,W_{y}=0)$ sector.

We use different numerical methods to treat exactly the Hamiltonian Eq. 1. In Sec. III, we first use methods which allow to obtain exact eigenstates of $H$. As usually performed in the MBL context, we concentrate on eigenstates with eigenvalue $E$ located right in the middle of the many-body spectrum, that is at $\epsilon=0.5$ with $\epsilon=(E-E_{\rm min})/(E_{\rm max}-E_{\rm min})$ where $E_{\rm min/max}$ are the extremal energies of $H$. Full numerical diagonalization is possible up to $N=36$ sites. We also use massively parallel shift-invert diagonalization Pietracaprina et al. (2018) that allows to target eigenstates at $\epsilon=0.5$ on larger samples up to $N=52$. We typically obtain about $100$ eigenstates near $\epsilon=0.5$. Results are averaged over more than 500 (up to 10.000) realizations of disorder. To show the interest of working with constrained models, the Hilbert space of the largest cluster on which we can obtain exact eigenstates ($N=52$) is almost twice as large as the one for $N=24$ spins 1/2 (in the zero magnetization sector), which constitutes the state of the art for the analysis of MBL in spin chains Pietracaprina et al. (2018); Macé et al. (2019, 2019). Obtaining $10$ eigenstates in the middle of the spectrum for one disorder sample at $N=52$ requires about $20$ minutes on 10.000 CPU cores.

In Sec. V, we also consider dynamical properties of various observables after a quench. The knowledge of all eigenstates allows to consider properties at arbitrary time up to $N=36$. We also use Krylov expansion time evolution techniques Nauts and Wyatt (1983), which allow simulations up to time of the order of thousand plaquette flips on large clusters up to $N=64$. Note that simulations on the largest sample $N=64$ are extremely demanding (see sizes of Hilbert space in Tab. 1). In the study of dynamics, we average over between 30 and 1000 realizations of disorder for each simulation set.

III Properties of eigenstates

III.1 Spectral statistics

We first consider the simplest property of a many-body system, namely its energy spectrum. The nature of the distribution of energy levels can already discriminate between different regimes: ergodic systems display level repulsion while localized systems display a random distribution of level spacing characterized by Poisson statistics. In the many-body context, and in order to counter the fact that the density of states is not uniform, it is useful to consider the distribution of consecutive ratios as introduced in Ref. Oganesyan and Huse, 2007. We define gaps in the many-body spectrum as $s_{n}=E_{n}-E_{n-1}$ and consider the consecutive reduced gap ratio $r_{n}=\frac{\min{(s_{n},s_{n+1})}}{\max{(s_{n},s_{n+1})}}$, for which $0\leq r\leq 1$. The distribution $P(r)$ and average value of $\langle r\rangle=\int_{0}^{1}rP(r)dr$, averaged over eigenstates and disorder realizations, have been computed for random matrix ensembles (corresponding to ergodic systems) and for independent energy levels (Poisson statistics) Atas et al. (2013). For real Hamiltonians such as Eq. 1, the random matrix ensemble to be considered is the Gaussian Orthogonal Ensemble (GOE). Fig. 2 displays both the distribution $P(r)$ and its average $\langle r\rangle$ for various values of disorder and system sizes. We compare them to predictions for the GOE ensemble Atas et al. (2013) $P_{\rm GOE}(r)=\frac{27}{4}\frac{r+r^{2}}{(1+r+r^{2})^{5/2}}$, $\langle r_{\rm GOE}\rangle\simeq 0.5307$ and Poisson statistics $P_{\rm Poisson}(r)=\frac{2}{(1+r)^{2}}$, $\langle r_{\mathrm{P}oisson}\rangle\simeq 0.38629$. Considering first the distribution (Fig. 2b), we find that the agreement with $P_{\rm GOE}(r)$ for low values of disorder ($V=3$) and $P_{\rm Poisson}(r)$ for high values ($V=25,30$) is excellent, with no adjustable parameter. Intermediate values of disorder $V=12,17$ present distributions which display a crossover behavior between the two limit distributions for the system size considered $N=40$. The average gap ratio $\langle r\rangle$ also displays an interesting crossover between the two limiting cases, with different system sizes showing an apparent crossing point (see inset of Fig. 2a for square samples $N=32,36,40$) around $V\simeq 15-20$. We note that rectangular samples also show the same crossover between the GOE and Poisson limits, but do not cross exactly at the same value which can be attributed to their different aspect ratios 222Considering only the clusters with $L_{x}=4$ leads to a smaller critical disorder $V\sim 12$ for the one-dimensional geometry of this four-leg tube. Results on other observables are consistent with this estimate.. The critical value of the gap ratio $\langle r\rangle^{*}\simeq 0.392$ is smaller than for the 1d standard MBL model Luitz et al. (2015), indicating that the putative transition point looks even closer to the Poissonian localized limit. Note that for the largest disorder $V=30$ considered here, we also checked on the largest system available to full diagonalization ($N=36$) that no band-like structure is present in the energy spectrum, a feature which was found Papić et al. (2015) to mimic MBL on too small system sizes.

III.2 Eigenstate statistics

We next consider the similarity between eigenstates at the same energy density. It is well-known that eigenstates in an ergodic phase look ’similar’ (e.g. they present similar expectation values for local observables), while eigenstates in a MBL phase are very different (e.g. local observables can differ strongly from one eigenstate to the following one in the many-body spectrum). To quantify the degree of similarity, it is useful to consider Luitz et al. (2015) the Kullback-Leibler divergence between two eigenstates $|n\rangle$ and $|n^{\prime}\rangle$ located nearby in the spectrum (here at the same energy density $\epsilon=0.5$). It is defined as $\mathrm{KL}=\sum_{i}p_{i}\log(p_{i}/p^{\prime}_{i})$ where $p_{i}$ and $p^{\prime}_{i}$ are the participation coefficients, i.e. the square amplitudes $p_{i}=|\langle n|i\rangle|^{2}$ , $p^{\prime}_{i}=|\langle n^{\prime}|i\rangle|^{2}$ of eigenstates in a given basis (here the computational basis where dimer occupations are diagonal). The sum runs over the full Hilbert space, and the eigenstates are normalized such that $\sum_{i}p_{i}=\sum_{i}p^{\prime}_{i}=1$. Eigenstates at ‘infinite temperature’ ($\epsilon=0.5$) in an ergodic phase should be very close to eigenstates of GOE random matrices, for which one has $\langle\mathrm{KL}\rangle_{\rm GOE}=2$. On the other hand, for eigenstates which are many-body localized, we expect $\langle\mathrm{KL}\rangle$ to diverge with system size as the eigenstates are increasingly different. Fig. 3 precisely shows this behavior, with expectation values $\langle\mathrm{KL}\rangle$ converging to their limiting GOE value at low disorder, and diverging for large disorder. We also observe a crossing point at $V\simeq 15$, albeit slightly drifting, for different system sizes. This drifting behavior is also observed for 1d MBL in the random-field Heisenberg model Luitz et al. (2015). The value of $\langle\mathrm{KL}\rangle$ at the transition point is larger than for the 1d spin model Luitz et al. (2015), pointing again to a putative critical point closer to the localized limit.

III.3 Dimer occupation and columnar imbalance

*Dimer occupation — * A qualitative understanding of the behavior of eigenstates in the different regimes is provided by considering local observables, and in particular by lattice bond occupation by dimers. For each lattice bond $b$, we compute the dimer occupation $\langle n|n_{b}|n\rangle$ in the eigenstate $|n\rangle$ with $n_{b}|i\rangle=1$ if the bond $b$ is occupied by a dimer in the basis configuration $|i\rangle$, [math] otherwise.

Fig. 4 represents graphically the dimer occupation $\langle n|n_{b}|n\rangle$ in three different eigenstates, consecutive in the many-body spectrum, for three values of disorder strength. In the low-disorder phase ($V=1$), all dimer occupations are equally likely: $\langle n_{b}\rangle\simeq 1/4$ for all bonds and all eigenstates. The dimer occupations at large disorder ($V=30$) mainly overlap with a single dimer configuration, which differs totally from one eigenstate to the next one. Resonances between configurations which differ by one or two plaquette flips can be evidenced by almost equal occupation of parallel dimers around a plaquette in some eigenstates (e.g. left-middle plaquette for the left-most eigenstate at $V=30$). At intermediate disorder ($V=17$), while for the two first left eigenstates dimer occupations look similar (at a superficial level), the right eigenstate configuration is markedly different.

*Columnar imbalance — * We now consider expectation values for another local observable, which will be useful later when discussing dynamics in Sec. V. We define the imbalance with respect to some specific state, the columnar configuration represented in the inset of Fig. 5, as:

[TABLE]

The phase $\varphi_{p}$ carried out by each plaquette is simply chosen to be [math] or $\pi$ on a columnar arrangement (see inset of Fig. 5) such that the imbalance of this perfect columnar state is maximal and equal to ${\cal I}=1$. At thermal equilibrium and infinite temperature, the imbalance vanishes on average.

We display in Fig. 5 the probability distribution $P({\cal I})$ of $\langle n|{\cal I}|n\rangle$ taken over eigenstates $|n\rangle$ located at $\epsilon=0.5$ for different values of disorder and computed on the tilted square lattice with $N=40$ sites. For small values of disorder (left column), the probability distribution is gaussian around the mean value [math], as expected for an ETH phase. The typical width of the distribution (which is very small for the smallest disorder $V=1$) broadens with disorder. Close to $V\simeq 15$, peaks for specific values of the imbalance start to appear in the distribution which become clearly visible for all strong values of disorder $V\geq 20$ (right panels). The distribution of imbalance for pure dimer configurations $P_{\rm conf}({\cal I})$ (obtained by computing ${\cal I}$ for all dimer configurations in the Hilbert space) is represented in the $V=30$ panel for comparison (there are $1/N$ values possible for the imbalance ${\cal I}$). The location of the peaks of $P({\cal I})$ matches precisely those of $P_{\rm conf}({\cal I})$, with however different amplitudes. This indicates that at large $V$, eigenstates are close (but not exactly equal) to pure dimer configurations, in agreement with the qualitative view provided by the snapshots of Fig. 4.

III.4 Participation entropies

To understand the localization properties of the eigenfunctions in the Hilbert (configuration) space, we next consider the participation entropy

[TABLE]

where $p_{i}=|\langle n|i\rangle|^{2}$ and the sum runs over all basis states $|i\rangle$ of the computational basis. At low disorder, we expect eigenfunctions to be fully delocalized with a participation entropy scaling as $S_{p}=a\log({\cal H})-\dots$ (with $a=1$ for a fully-delocalized state and $\dots$ denoting finite size corrections). In the large disorder regime on the other hand, based on data for the MBL phase in the 1d random field Heisenberg Macé et al. (2019), we do not expect a true localization in configuration space but instead a volumic, multifractal, behavior where $S_{p}=a\log($H$)$ with $a$ strictly less than 1 (and increasingly small as disorder is increased). Fig. 6 presents the numerical data for the average participation entropy divided by the log of the size of the Hilbert space $\langle S_{p}\rangle/\log({\cal H})$. Two regimes can be clearly distinguished: one at low $V$ where data tend to their maximal value ($a=1$) with system size, and another at larger disorder where curves tend to group together towards sensibly smaller values of the volume law coefficient $a$. This is in accordance with the results for the MBL phase in one dimension Macé et al. (2019). We find that data for $S_{p}/\log{\cal H}$ for different sizes tend to regroup themselves around $V\simeq 15-20$, consistent with previous estimates of the possible transition point.

III.5 Entanglement entropy

Many-body eigenstates in ergodic or localized phases have strikingly different entanglement properties. Consider the von Neumann entanglement entropy of an eigenstate $|n\rangle$ for a system divided in two parts $A$ and $B$:

[TABLE]

with $\rho_{A}$ the reduced density matrix obtained by tracing out bond degrees of freedom in B $\rho_{A}=\mathrm{Tr}_{B}|n\rangle\langle n|$. For our samples, we choose equal parts $A$ and $B$ containing identical number of lattice bonds and sites (see Appendix for a graphical representation of the entanglement cut).

The entanglement entropy $S$, averaged over eigenstates and realizations of disorder, is represented in Fig. 7. In the ETH phase, a volume law of the entanglement entropy is expected with $S$ scaling with the number of degrees of freedom in $A$. This is what is found in the low-$V$ regime of Fig. 7, with the following caveat: due to the constrained nature of the Hilbert space, not all configurations of $A$ are compatible with those of $B$ (see a discussion in Ref. Morampudi et al., 2020 and in the Appendix) and the number of these ‘constraint sectors’ in $\rho_{A}$ depends in particular on the area between $A$ and $B$. This results in different samples with the same volume having different constraints on $\rho_{A}$ when bipartitioned. This is clearly seen when computing entanglement entropy on random states in the different samples denoted as $S^{\rm random}$ and represented by the star symbol in Fig. 7: for instance, the $8\times 4$ rectangular lattice has a much higher entropy than the $4\times 8$ lattice, even though both have the same volume (number of bonds and sites) for $A$. Fig. 7 shows that the low-$V$ entanglement entropies of all samples converge to the random-state values, confirming clearly that the low-$V$ phase fulfills the ETH. For larger values of $V$, the entanglement entropy drops down for all samples. The inset of Fig. 7 shows a reasonable scaling of $S$ with the area ${\cal A}$ at large enough $V$, in agreement with the expected area law scaling expected in a localized phase. Finite-size effects combined with the constraint effect on $S$ discussed above prevent us from pinpointing a transition/crossover point between the volume and area law scaling with great precision.

In previous studies of the behavior of entanglement entropy close to the MBL transition, the variance of entanglement entropy was shown Kjäll et al. (2014); Khemani et al. (2017a) to peak when approaching the transition, which can be rationalized as the result of the crossover between the volume (area) law of entanglement entropy in the ETH (MBL) phase. This is also seen in Fig. 8 where the standard deviation ${\rm std}(S)$ of the entanglement entropy (computed over all eigenstates and realizations of disorder), scaled by the entanglement entropy of a random state $S^{\rm random}$, is represented. Focusing on the square samples, a right-shift of the peak with system size towards larger values of $V$ is apparent: a similar shift was always found to be present in numerical studies of the MBL transition in 1d Kjäll et al. (2014); Khemani et al. (2017a). For the largest square (rectangular) sample, the peak is located at $V=10$ ($V=12$), providing a lower bound for the possible transition.

IV Machine learning analysis of entanglement spectra

As a complementary approach, we study the quantum dimer model Eq. 1 using machine learning techniques and follow the supervised scheme introduced in Ref. Carrasquilla and Melko, 2017. As inputs representative of the two phases, we provide entanglement spectra (i.e. the eigenvalues $\lambda_{i}$ of the reduced density matrix $\rho_{A}$) obtained deep in the ETH and the MBL phases and we train a neural network to classify them accordingly. The supervised learning has been used in multiple occasions in the context of delocalization-localization transitions Schindler et al. (2017); Li et al. (2017); Venderley et al. (2018); Hsu et al. (2018); Zhang et al. (2019); Huembeli et al. (2019) producing results in good qualitative agreement with more conventional analyses. In particular, Refs. Schindler et al., 2017; Venderley et al., 2018; Hsu et al., 2018 showed the interest of using labeled entanglement spectra.

For this scheme, we consider the $N=32,36,40$ square samples and the largest rectangular sample $N=6\times 8$. For each of them, we provide 10.000 entanglement spectra (2000 for $N=40$, $6\times 8$) including between 100 and 200 disorder realizations per disorder strength at $V=1$ for the ETH-labelled phase and $V=30$ for the MBL-labelled phase. The spectra being rather large (1972 values for $N=32$, up to 21286 for $N=6\times 8$), we found that sorting them allowed for both perfect training and test accuracies for all system sizes. Fig. 9 represents the form of the entanglement spectra used in these two limits for the sample $N=36$ for $~{}100$ different realizations of disorder. The very similar form for various disorder samples in the ETH phase clearly contrasts with the stronger dispersion observed in the MBL phase. In the later, superposing the entanglement spectra for various samples sizes (inset of Fig. 9) highlights that the larger values of the spectrum decay as a power-law, similar to what is found in 1D MBL Serbyn et al. (2016).

We used a fully-connected neural network consisting of one hidden layer of 32 neurons followed by two softmax output neurons. We follow a cross-validation procedure where we randomly selected half of the dataset to form the training dataset, the rest being assigned to the test set. This process is repeated multiple times, generating new training and test partitions each time. This allowed us to track whether the neural networks were overfitting depending on the training conditions. Namely, we checked that data from $V=1$ and $V=30$ give perfect training and test accuracies for each system size.

Fig. 10 displays features that are consistent with the previous analysis. The left panel displays the average output of the neural network for the different samples as a function of disorder strength. At low $V$, the machine learning analysis validates a fully-ETH phase (*i.e. * where all samples are classified ETH) that extends up to a value $V=V_{1}$, and at large $V$ a fully-MBL phase (with more than $99\%$ accuracy) at large for $V\geq V_{2}$ ($V_{1}\simeq 6$ and $V_{2}\gtrsim 20$ for the largest $N=6\times 8$ sample). Notice how these bounding values (in particular $V_{1}$) shift to higher values of disorder with system size. This reflects that ETH is easier to disrupt on a too small sample, in perfect agreement with the trend in all other observables discussed in Sec. III. We find no crossing point with system size in the current data, different from what is observed with the 2d Ising model Carrasquilla and Melko (2017).

The right panel of Fig. 10 shows the standard deviation of the neural network output as a function of disorder. The standard deviation is low in both limits where the phases are well distinguished (at low and large $V$), and peaks at an intermediate value of disorder. The location of the peak (which shifts with system size) is the point where the neural network has most difficulties to classify phases. It can be interpreted as a finite-size estimate of a possible transition point. Notice the similarity between the standard deviation of the neural network output and the standard deviation of the entanglement entropy (Fig. 8), in particular the positions of the peaks are almost the same for both quantities for the different sample sizes.

In conclusion of this section, we find that a neural network only fed with entanglement spectra is able to learn how to correctly distinguish the ETH and MBL phases for the quantum dimer model with random potential Eq. 1 as well as to provide finite-size estimates of the transition point between the two. This automated method gives results in very good agreement with the analysis based on more standard, feature-enginereed, estimators of the phases presented in the previous Sec. III. Our results also indicate that the machine learning analysis is also subject to finite-size effects, indicating that the use of numerical data coming from a single size may not be able to provide quantitative results for the study of phase transitions (see a discussion in Ref. Théveniaut and Alet, 2019). One noticeable interest of the neural network analysis (already pinpointed earlier Schindler et al. (2017); Venderley et al. (2018)) is that the required amount of data and overall computational cost is considerably lower than with more traditional observables to obtain approximately similar quality of prediction: for instance good statistics on the gap ratio (Fig. 2) requires $\simeq 100$ times more realizations of disorder than with the machine learning analysis (Fig. 10).

V Dynamical properties after a quench

We consider here dynamical properties after a quench of the model Eq. 1, in order to probe for the arrested transport dynamics or slow propagation of information expected in a MBL phase. We perform quenches from an initial product state (with no entanglement) and take here the columnar configuration depicted in the inset of Fig. 5 (a product state). This columnar configuration is the most flippable state and is well connected to other dimer configurations by the action of the Hamiltonian. We thus intuitively expect this state to rapidly thermalize if it does. The choice of a columnar configuration echoes the experimental protocols for which the system is often initalized in a charge density state Schreiber et al. (2015); Bordia et al. (2017), or numerical simulations which choose a Néel state in spin chains, both states sharing the same properties as the columnar configuration.

V.1 Dynamics of imbalance from the initial state

We start with the dynamics of the columnar imbalance defined in Eq. 2. By definition, ${\cal I}(t=0)=1$, whereas ${\cal I}=0$ when averaged over all eigenstates (see Fig. 5). In the long-time limit, ${\cal I}(t)$ is thus expected to vanish in a thermal phase, whereas a non-vanishing long-time value would signal a localized phase. Note that on finite systems, ${\cal I}(t\rightarrow\infty)$ never strictly vanishes. Fig. 11 displays the columnar imbalance as a function of time, with panels corresponding to different disorder strengths $V$. In the low-disorder limit, the imbalance quickly decays to a vanishing value on short time-scales (this is most readily seen at $V=1$). When disorder is increased, decays become slower and for intermediate disorder values $V=10,15,20$, the saturation value is not reached for the largest clusters within the time scale studied (for $N=32,36$ data not shown at longer time show that the saturation is reached within less than a percent for $t=1000$ for all disorders considered). We can nevertheless observe that the imbalance clearly decreases with system size for the largest time available, in a manner compatible with a vanishing value in the thermodynamic limit for $V\lesssim 15$. For larger values of disorder, one clearly observes that the saturation value is much larger and, within available time and cluster sizes, the imbalance does not tend to vanish: this is most speaking for $V=30$ where more than $80\%$ of the initial imbalance is retained even on the largest clusters. Pinpointing exactly a transition point with dynamics is difficult due to the time scale limitations for the largest clusters, but the present data appear to be compatible with a transition/crossover for values of $V$ between $15$ and $20$, similar to what is found with the analysis of eigenstate data in Sec. III.

V.2 Entanglement and participation entropies dynamics

The nature of the growth of entanglement entropy $S(t)$ after a quench is a hallmark feature to distinguish MBL in one-dimension, where it takes a logarithmic form with time Žnidarič et al. (2008); Bardarson et al. (2012), from an ETH regime (with ballistic growth $S(t)$ Kim and Huse (2013)) or an Anderson localized regime (where it saturates to a constant $O(1)$). In both MBL and ETH regimes, the saturation value at long time is expected to fulfill a volume law, with a prefactor equal to (smaller than) the corresponding thermal entropy for the ETH (MBL) case Serbyn et al. (2013b); Nanduri et al. (2014).

We represent in Fig. 12 the dynamics after a quench from the columnar configuration of the average entanglement entropy $\overline{S(t)}$ normalized by the entanglement entropy of a random state $S^{\mathrm{random}}$ (which is the thermal entropy up to $O(1)$ corrections) for different disorder strenghts. The ETH regime with a rapid growth of entanglement is readily identified at low disorder $V=1$, where $\overline{S(t)}$ reaches its saturation value (close to the maximum value reached by a random state) in less than $~{}10$ plaquette flips. In the top panels where the same scale is used for the vertical axis, the entanglement entropy is observed to grow slower towards its saturation value, with visible finite-size effects. Again the time and size limitations do not strictly allow to draw definitive conclusions, but it appears very likely that for large enough clusters, the entanglement entropy will reach its thermal value at long time for $V\lesssim 15$. On the bottom panels (disorders $V=20,25,30$) on the other hand, the entanglement entropy develops at an even much slower pace. The logarithmic horizontal scale chosen for the bottom panels highlights the striking feature that the growth of entanglement entropy appears to behave as $\log(t)$ for the largest clusters which are not hampered by finite-size saturation effects. Note that for the largest data at $V=15$, a logarithmic growth is also compatible with the data (not shown). Naively exporting the argument (based on the exponentially-decaying interactions between local integral of motions) that explains the logarithmic growth of entanglement in 1d MBL, we would also expect a logarithmic growth in 2d for the entanglement cut geometry chosen here, as entanglement can grow only in one direction.

We finally present the evolution after a quench of the participation entropy $\overline{S_{p}(t)}$ in Fig. 13. The dynamics of the participation entropy is almost identical to the one of the entanglement entropy discussed above: at small disorder, rapid saturation to the maximum attainable value $\log({\cal H})$, slower growth at intermediate values and logarithmic growth for the largest samples for strong disorder (note again the logarithmic scale on the bottom panels of Fig. 13). The only noticeable difference in the dynamics of both entropies comes from different finite-size effects at short time, with participation data showing an improved nesting of curves as system size is increased which is not observed for entanglement. From participation entropy, we conclude again on slow dynamics starting from $V\gtrsim 15-20$ consistent with the existence of a localized phase.

VI Discussion and perspectives

We presented a comprehensive and extensive large-scale exact numerical study of the localization properties of a two-dimensional constrained quantum many-body system with disorder. All data, both for eigenstate or dynamical properties, can be interpreted as consistent with the existence of two distinct phases: an ETH phase at low disorder and a many-body localized phase at strong disorder, which are separated by a transition located around disorder strength $V\simeq 15-20$.

Our evidence for a MBL transition in the 2d quantum dimer model comes from numerical simulations on finite lattices. Of course, finite-size simulations can always be argued to artificially detect a MBL phase even when only a ETH phase occurs in the thermodynamic limit. We would like to emphasize that the level of numerical evidence for a MBL transition in the 2d quantum dimer model is in our opinion and experience similar to the one obtained for the standard model of MBL in one dimension Pal and Huse (2010); Luitz et al. (2015), with similar or larger Hilbert space sizes and time scales probed in the current work. This is of course not a definite proof that a MBL transition occurs in the thermodynamic limit.

What is perhaps more important with respect to potential experiments is the fact that even if the behavior at large values of disorder (say $V\geq 25$) is ultimately ergodic, the time scales and / or system sizes needed to probe ergodicity would be extremely large. An experimental realization of Eq. 1 or of a similar constrained model in 2d would see localization for all practical purposes on the time scales available in the lab. We deliberately avoided to attempt to perform a data collapse analysis to extract e.g. critical exponents due to the two following limiting points: the moderate linear system sizes that we can reach as well as the non-uniform aspect ratio between available samples would likely provide biased estimates for critical exponents. For one-dimensional MBL, the larger linear length scales that can be reached have been argued (see e.g. Ref. Khemani et al., 2017b) not to be large enough to provide correct estimates of asymptotic critical behavior. The situation is likely to be the same here. With these numerical limitations in mind, we can nevertheless observe that critical values of the gap ratio or of the Kullback-Leibler divergence of eigenstates are closer to their Poisson than their ETH limits, indicating that the putative transition point is even less ergodic than for the one-dimensional MBL transition in the random-field Heisenberg spin chain Luitz et al. (2015).

Going back to the *l-*bits arguments which predict the asymptotic unstability of MBL in two dimensions De Roeck and Huveneers (2017); De Roeck and Imbrie (2017); Potirniche et al. (2019), we remark that we do not know how to simply construct l-bits for the quantum dimer model. The crucial point is that there is no tensor product structure for the Hilbert space of quantum dimer models, thus constructing a basis of commuting operators with local support appears ill-defined. As the bubble argument is explicitly based on *l-*bits, it cannot strictly speaking apply here.

Regarding experiments, the relative strong disorder needed to induce the MBL phase does not provide immediate relevance of our results in frustrated magnets, where QDMs naturally emerge. However, the artificial realizations of quantum dimer models (and other related 2d constrained models) have been extensively studied in cold-atomic and mesoscopic systems. Precise experimental protocols have been proposed using Ryberg atoms Glaetzle et al. (2014); Celi et al. (2020), trapped ions Nath et al. (2015), large-spin ultracold atoms Sundar et al. (2019), Josephson Junction arrays Ioffe et al. (2002) or superconducting quantum circuits Marcos et al. (2014). Further explorations are needed to understand which experimental platform is the most suited to implement random potential or interactions. Constrained models in 1d have already been implemented (with no randomness) exploiting Rydberg blockade in atomic chains Bernien et al. (2017), already exhibiting rich unusual quantum dynamics.

There are several perspectives opened up by our work. First, the roadmap to two-dimensional MBL can be exploited using the QDM on other two-dimensional lattices, allowing to test for universality and to search for other features. A recent investigation Pietracaprina and Alet (2020) considered the honeycomb lattice, which has an effective smaller local Hilbert space allowing to reach larger samples, and obtained results similar to ours. The QDM on the kagome lattice is also an interesting case worth pursuing at it possesses conserved $Z_{2}$ quantum numbers, allowing the exciting possibility of 2d topological order in MBL states Moessner and Sondhi (2001). The possibility of MBL in other quantum constrained models such as quantum ice or loop models also provide an interesting follow-up.

The large-scale numerical techniques that we use in the present work would be useful to investigate other aspects of the model Eq. 1 not addressed here, in particular to make connection with 1d models where MBL is best understood. The current investigation focused on the ’infinite temperature limit’, that is excited states in the middle of the spectrum ($\epsilon=0.5$). It would be interesting to construct a full energy-disorder phase diagram to see if a many-body mobility edge exists as in 1d spin chain models Luitz et al. (2015). Also the detailed dynamics occurring in the ETH phase remains to be investigated. The relatively slow dynamics exhibited around $V\sim 10,15$ in Figures 11, 12, 13 could be signs of a precursor subdiffusive regime, similar to what is observed in 1d systems Agarwal et al. (2015); Luitz et al. (2016); Žnidarič et al. (2016). Finally, it would be intriguing to see if the machine learning techniques used in Sec. IV would be able to distinguish 1d from 2d MBL, and if not, to use neural networks trained on 1d spin chain models to probe 2d MBL (transfer learning).

On the analytical side, the absence of an obvious scheme for constructing *l-*bits points towards the interest of building an alternative effective description for the many-body localized phase in such 2d constrained models. Recall as well that there is no free-particle limit to start with in the model Eq. 1. The height description Blote and Hilhorst (1982) of dimer coverings (and of constrained models in general) appears as a good candidate to develop such a framework. Finally from a purely classical perspective, we note that the strong disorder limit $V\rightarrow\infty$ constitutes an interesting classical statistical mechanics problem (with for instance non-trivial degeneracies), unsolved to the best of our knowledge and distinct from other interacting dimer models Alet et al. (2005).

Acknowledgements.

We gratefully thank Juan P. Garrahan and Stephen Powell for useful discussions at an initial stage of this project. We also thank Nicolas Laflorencie for useful comments on the manuscript and Nicolas Macé for an insightful remark on the scaling of entanglement growth. H.T. is supported by a grant from the Fondation CFM pour la Recherche. Z.L. acknowledges support from EPSRC Grant No. EP/M019691/1. This work benefited from the support of the project THERMOLOC ANR-16-CE30-0023-02 of the French National Research Agency (ANR) and by the French Programme Investissements d’Avenir under the program ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT. We acknowledge PRACE for awarding access to HLRS’s Hazel Hen computer based in Stuttgart, Germany under grant number 2016153659, as well as the use of HPC resources from CALMIP (grants 2017-P0677 and 2018-P0677) and GENCI (Grant 2018-A0030500225). The shift-invert Pietracaprina et al. (2018) and Krylov time evolution numerical calculations are based on the linear algebra libraries PETSc Balay et al. (1997, 2017), SLEPc Hernandez et al. (2005), MUMPS Amestoy et al. (2001, 2006) and STRUMPACK Ghysels et al. (2016, 2017). The neural network analysis used the library TensorFlow Abadi et al. (2015).

Appendix A Lattices

A.1 Lattice geometries

We use regular and tilted square lattices as well as rectangular lattices with $N$ sites and $2N$ bonds. Regular square lattices have $N=L^{2}$ sites, and we used samples with $L=4,6,8$. Rectangular lattices, which have $N=L_{x}\times L_{y}$ sites (with $L_{x,y}=4,6,8,10$), have a non-unity aspect ratio, which can affect the analysis of finite-size effects. Finally, tilted square lattices host $N=p^{2}+q^{2}$ sites, and we use the combination $p=4,q=4$ (with $N=32$ sites), $p=6,q=2$ ($N=40$), and $p=6,q=4$ ($N=52$). Tilted square lattices have a unity aspect ratio, but they do not have all symmetries of the full regular square: this is however irrelevant as the Hamiltonian Eq. 1 does not allow any lattice symmetries due to disorder. All these lattices, which are represented in Fig. 14, have in common that they host a $(W_{x},W_{y})=(0,0)$ winding sector: this is true when $L,L_{x},L_{y},p$ and $q$ are even. One can also easily define columnar state on these samples (see example for $N=52$ in Fig. 14).

A.2 Entanglement cut

When computing bipartite entanglement properties, there is some freedom in choosing where is located the entanglement cut, i.e. which degrees of freedom belong to $A$ or $B$ in Eq. 4. We choose the cuts displayed in Fig. 14, where all lattices are divided in two parts with equal number of bonds. With this cut, the area ${\cal A}$ between $A$ and $B$ is expected to scale with $L$ for regular square lattices, $L_{y}$ for rectangular lattices, and $\sqrt{p^{2}+q^{2}}$ for tilted square lattices.

This choice of cuts for the regular square, rectangular and $N=32$ tilted square lattice also induces the following property: one of the two windings numbers (the one associated to a line parallel to the cut, e.g. $W_{y}$ for the rectangular lattices) is conserved by the cut, but not the other one. The reduced density matrix on $A$ is block-diagonal with blocks labelled by winding numbers in $A$ ($W_{y}^{A}$ which is forced to be [math], and $W_{x}^{A}$), which we take advantage of in numerical computations. The existence of blocks in the reduced density matrices and its effect on entanglement has been discussed in other models with conservation laws Vidmar and Rigol (2017); Garrison and Grover (2018); Zhou and Luitz (2017); Morampudi et al. (2020).

We note also a supplementary specificity of constrained models for entanglement properties (see discussion in Ref. Morampudi et al., 2020). Due to the constraint and even within the conserved-number blocks (e.g. configurations of $A$ with $W_{x}^{A}$ and configurations of B with $W_{x}^{B}$ such that $W_{x}^{A}+W_{x}^{B}=W_{x}=0$), not all configurations of $A$ are compatible with those of $B$. This property has an important consequence, for instance when computing volume-law entanglement, as different samples with the same volume can have different constraints on configurations in $A$ and $B$ when bipartitioned. This is illustrated in the main text with the example of the $8\times 4$ and $4\times 8$ rectangular samples. This property can also cause difficulties in comparing scaling of entanglement with system size for different samples.