Many-body localization in a fragmented Hilbert space

We study many-body localization (MBL) in a pair-hopping model exhibiting strong fragmentation of the Hilbert space. We show that several Krylov subspaces have both ergodic statistics in the thermodynamic limit and a dimension that scales much slower than the full Hilbert space, but still exponentially. Such a property allows us to study the MBL phase transition in systems including more than $50$ spins. The different Krylov spaces that we consider show clear signatures of a many-body localization transition, both in the Kullback-Leibler divergence of the distribution of their level spacing ratio and their entanglement properties. But they also present distinct scalings with system size. Depending on the subspace, the critical disorder strength can be nearly independent of the system size or conversely show an approximately linear increase with the number of spins.


I. INTRODUCTION
Many-body localization (MBL) and its transition [1][2][3][4] have been the subjects of numerous studies over the recent decades. They are directly related to core physical concepts and properties of the physics of closed quantum systems, namely thermalization, transport, and the effects of disorder. Interacting systems at weak disorder thermalize and present ergodic features seemingly following the so-called strong eigenstate thermalization hypothesis (ETH) [5][6][7] . It states that, at high energy, a generic closed quantum system has all its eigenstates display thermal values for all local observables. At strong disorder, on the other hand, theoretical arguments and numerical studies show a breakdown of the ETH in one dimensional systems, arising from emergent integrability and approximate integrals of motions [8][9][10][11][12][13][14][15][16] . The high-energy eigenstates are then characterized by low entanglement, following an area law instead of a volume law as in the ergodic phase [17][18][19][20][21][22] . While the existence of the MBL phase is now generally well accepted, despite a recent debate [23][24][25][26][27][28][29] , describing the transition from the ergodic phase to the localized phase remains a numerical and theoretical challenge. Indeed, the physics of MBL arises from the rich interplay of various phenomena: many-body interactions allowing for nonintegrability, (strong) disorder leading to localization, and high-energy physics. Many theoretical studies rely on phenomenological renormalization group arguments, based on various physical arguments on thermalization and predictions of random matrix theory [30][31][32][33][34][35][36][37][38][39] . Numerical studies are also limited by the complexity of the problems: matrix-product-states based approaches [40][41][42][43][44][45] perform well deep inside the MBL phase but become unreliable close to the transition due to rapidly increasing entanglement, leaving exact diagonalization and variants thereof, with its generally limited system sizes, as the main source of exact numerical resources 27,[46][47][48][49][50][51][52] . The existence of MBL, as a means to break the strong ETH beyond integrability, spurred the growth of interest in other phenomenas leading to such a breakdown. [53][54][55][56][57] Two majors archetypes have emerged: many-body quantum scars [58][59][60][61][62][63][64][65][66] and Krylov fragmentation [67][68][69][70][71][72][73] . Systems with quantum scars present a set of measure zero of highly excited non-thermal eigenstates, typically characterized by a sub-volume law entropy. The other eigenstates remain thermal. The presence of these states has especially strong consequences on nonequilibrium dynamics in such systems, with partially suppressed thermalization 58,61 . Depending on the initial state, time-evolution under a Hamiltonian presenting these scar states can typically present much slowler relaxation of observables towards the thermal equilibrium states, with slowly suppressed revivals at long times. Generic methods to embed such states into a thermal spectrum have been proposed 74,75 and scars have been proved to be resilient to the effect of disorder 66 . More relevant to this work is the concept of Krylov subspaces or Hilbert space fragmentation [67][68][69][70][71][72][73]76 . Due to the interplay between different U (1) symmetries such as charge and dipole conservation, each symmetry sector of the Hilbert space shatters into an exponential number of sectors or Krylov subspaces that are not connected by the Hamiltonian dynamics. Importantly, these subspaces are not fully labeled by quantum numbers. The exponential number of small disconnected sectors leads to anomalous and effectively localized dynamics 70,71,76 . Conversely, exponentially large Krylov subspaces have recently been the subject of several studies. 68,73 Remarkably, in the same model and in the absence of disorder, some of these Krylov subspaces follow the ETH, while other subspaces have completely integrable statistics.
A natural question for these systems is the effect of disorder on these Krylov subspaces, and in particular whether it can preserve the fragmented nature of the Hilbert space, and lead to a localization of the different arXiv:2011.04659v1 [cond-mat.dis-nn] 9 Nov 2020 ergodic subspaces. More importantly, we identify sets of ergodic subspaces whose dimension grows much slower than the total Hilbert space dimension, albeit still in an exponentional fashion. This gives us the possibility to investigate through exact (and full) diagonalization one-dimensional systems of unprecedented physical sizes. A similar approach was proposed in Ref. 77: in the conventional XXZ model, strong disorder permits to approximately separate the Hilbert space into quasiindependent subspaces. There, the separation is only a strong-disorder induced approximation. In our model, it is exact at all disorder strengths. Our approach also shares some similarities with the studies of frustrated models such as quantum dimers whose Hilbert space shows slow exponential scalings and signs of MBL even in two dimensions 78,79 .
The outline of our paper is as follows. In Section II, we introduce the pair-hopping model and its main properties and symmetries. Section III is dedicated to Krylov subspaces. After a formal definition and a discussion of the Krylov subspaces studied in Ref. 68, we introduce our slowly-growing ergodic subspaces. We study the level spacing ratio statistics in the pair-hopping model in Section IV. We discuss the distribution of the level spacing ratio of the full Hilbert space, showcasing the need to consider the individual Krylov subspaces. To probe the ETH-MBL transition within the Krylov subspaces, we rely on the Kullback-Leibler divergence of the level ratio distribution and the reference GOE and Poisson distributions 4,47,80 . We identify the critical disorder strength as the point of maximum confusion. We show that our Krylov subspaces present all signs of the MBL phase transitions. The different critical disorders nonetheless have radically different scalings with system size, ranging from quasi absent finite-size effects to an approximately linear shift (within the size we have access to). Our results underline the importance of the structure of the Hilbert space in the behavior of the MBL phase transition. We also briefly discuss the existence of a mobility edge 1,2,48 . We then turn to the von Neumann entanglement entropy (vNEE) of highly excited states in Section V. The critical disorder strengths, identified by the maximum of the standard deviation of the midchain entanglement entropy 20 , are in qualitative agreement with our previous results. We also carefully discuss the scaling of the vNEE with subsystem size which present unusual plateaus due to the strongly constrained nature of our model.

II. PAIR HOPPING MODEL IN A TRANSVERSE FIELD
The pair hopping model 68,69,81 is an interacting model of spinless fermions on a one-dimensional lattice whose Hamiltonian is given by: where c j (c † j ) is the fermionic annihilation (creation) operator on site j, J j are site-dependent pair hopping terms and we fix the number of sites to L. We take J j uniformly sampled in [0.9, 1.1] in order to break inversion symmetry and translation symmetry. For convenience, we perform a Jordan-Wigner transform, and work with the spin-1 and consider either open boundary conditions (OBC) or periodic boundary conditions (PBC) in the spin basis. PBC in the spin basis are not equivalent to fermionic PBC due to the presence of the fermionic string, but the observables we consider in the remainder of this article are largely unaffected. This model preserves the total polarization (or charge in the fermionic language) and the dipole moment (or center-of-mass position) defined as Additionally, for L even with PBC and all L's with OBC, the pair hopping terms preserve the sublattice symmetry, i.e., P o − P e where P e (P o ) is the total charge of the even (odd) sites, and therefore these two charges are also conserved quantities. We denote with p e , p o , p tot and c the quantum numbers respectively associated to P e , P o , P tot and C. Similar pair-hopping terms appear naturally in different experimentally relevant set-ups, such as electrons in a Landau level 69,82,83 and in the Wannier-Stark problem [84][85][86] .
We introduce disorder in the form of a transverse field (a random on-site chemical potential in the fermionic language), resulting in the total hamiltonian where W j is taken uniformly in [−W, W ]. This disorder does not break any of the identified symmetries.
To better understand the physics of the hopping terms, it is convenient to represent the system in terms of pair of spins, using the notations introduced in Ref. 68. For convenience, we assume L even in the rest of the paper. Let the local Hilbert space be defined by The system is composed of N = L/2 pairs of spins, decomposed in the basis |↑ = |10 , |↓ = |01 , |− = |00 and |+ = |11 . (7) We denote |↑ and |↓ as pseudo-spins, |+ and |− as fractons and the combination |+− and |−+ as dipoles.
While σ z 's are still diagonal in this basis, the hopping terms of Eq. (1) lead to some complex algebra. The transformation rules are the following: Pseudo-spins exchange with each other (Eq. (8)), and dipoles can move each in one flavor of pseudo-spins (Eq. (9)). Conversely, well-chosen trio of fractons can transform into a fracton and a pair of pseudo-spins, and reversely the presence of a fracton may lead a pair of pseudo-spins to transform into a pair of fractons.

III. KRYLOV SUBSPACES
In this Section, we introduce the notion of Krylov subspaces, vector subspaces of the symmetry-resolved Hilbert space which are stable under the application of the Hamiltonian. In particular, we present several families of Krylov subspaces, ergodic at zero disorder despite their simple structure, whose dimension grows slowly with system size.

A. Definition
A natural approach to find the stable subspaces induced by the fragmentation of the Hilbert space is to work with Krylov subspaces 67,68,[70][71][72][73] . Generally, the Krylov space K 1 (|Ψ , H) generated by the state |Ψ and the Hamiltonian H is defined as follows: This definition has the drawbacks of being numerically unstable for generic Hamiltonians and requiring the orthogonalization of a set of large vectors. Working in the configuration basis it is more convenient to use the modified definition: K 2 (|Ψ , H) = Span({|w ∈ B|∃n ∈ N, w|H n |Ψ = 0}). (14) For simplicity, we consider in the following |Ψ ∈ B, i.e., a product state in the σ z 's basis. Then, seeing the Hamiltonian as a graph in the configuration space B, where each node is an element of the basis and each link a non-zero coefficient of the Hamiltonian, K 2 (|Ψ , H) is simply the connected subgraph including |Ψ . This definition is numerically stable and straightforward to compute. Note that contrary to K 1 , K 2 strongly depends on the basis B, which should be specifically crafted for the studied model. Our choice for B for the pair-hopping model ensures the following property: In fact, K 1 and K 2 here nearly always coincide 87 . We therefore only work with K 2 (|Ψ , H) and drop the subscript and explicit dependence on the Hamiltonian in what follows.

B. Ergodic Krylov subspace with a single pair of dipoles
In Refs. 68 and 73, the authors observed that exponentially large subspaces with either Poissonian or ergodic statistics might coexist in the absence of disorder. An especially convenient family of exponentially large ergodic subspaces were generated by the pair of dipoles |− + +− or |+ − −+ in a sea of pseudo spins |↑ and |↓ . In the absence of dipoles, the pair-hopping Hamiltonian acts in the sea of pseudo spins, with conservation of each flavour of pseudo-spins. Introducing a pair of dipoles breaks down integrability, as each can only move through a single flavour of pseudo-spins. The pair-hopping Hamiltonian acting on this subspace conserves the number of each flavor of pseudo-spin and fracton, leading to a remarkably simple structure of the corresponding Krylov subspace.
The first Krylov subspace we consider is therefore generated by the state where (w) n marks that we repeat n times the sequence w. The system is comprised of N = 4n + 4 pairs of spins and − + +− is placed exactly at the center of the chain. This Krylov space was shown to be ergodic in the absence of a transverse field 68 and its dimension can be readily computed. For simplicity, we first consider periodic boundary conditions. The dipoles |−+ and |+− can be seen as separating the pseudo-spins into two sequences: either in between |−+ and |+− or outside of them. We denote by N ↑ = 2n the conserved total number of pseudo-spins ↑, and by n ↑ (resp. n ↓ ) the number of ↑ (resp. of ↓) between |−+ and |+− . As a concrete example, the state belongs to K( Ψ 1 n , H PBC ) with n = 4, N ↑ = 8, n ↑ = 2 and n ↓ = 3. The dimension of the Krylov subspace is simply given by: using the Chu-Vandermonde identity for simplification, and the Stirling's approximation for the factorial. The dimension of this Krylov subspace therefore scales as 2 N = √ 2 L , i.e., much slower than the full Hilbert space's dimension, which scales as 2 L .
For OBC, the two dipoles separate the chain in three. There is a fixed number N L ↑ = n of pseudo-spins ↑ (↓) to the left (right) of the dipoles. The n ↓ ↓-pseudo-spins in between |−+ and |+− originally came from the left of |− + +− via acting by H OBC ; similarly the n ↑ ↑ in between pseudo-spins came from the right. Hence the state represented in Eq. (17) also belongs to K( Ψ 1 n ) OBC with N L ↑ = 4, n ↑ = 2 and n ↓ = 3. The dimension of K( Ψ 1 n ) OBC is thus given by Using twice the Chu-Vandermonde equality (see App. A 1), we obtain The exponential scaling is similar, albeit with a more favorable prefactor.
For both boundary conditions, the reduced Hilbert space scaling is not due to an extremal choice of quantum numbers (such as the linear scaling of the one particle sector of a particle number conserving U (1) model). It is a simple consequence of the presence of an extensive number (∝ N ) of freely exchanging pseudo-spins (composed of two real spins) that make most of the degrees of freedom. A table summarizing the properties of the Krylov subspace for numerically relevant values of N can be found in App. A 1. Due to the significantly smaller Krylov space's dimension, we focus on OBC for this family.

C. Slowly-growing ergodic Krylov subspaces
It is possible to construct a series of ergodic subspaces with even more favorable scaling with system size, which cannot be mapped to any simple quasi-particle picture. Working with sets of |↑↓ and pairs of dipoles |− + +− allows us to keep a simple analytical structure of the Krylov subspace while working in the sector (p o , p e , c) = (0, 0, 0). If we fix the number of dipoles to be constant when increasing system size, the dimension of the resulting Krylov space ultimately grows as 2 N . A natural way to go beyond this limit is to alternate between a finite number of pairs of pseudo-spins and the set of two dipoles. The less pseudo-spins, the slower the growth of the Hilbert space. The pair-hopping Hamiltonian in Eq. (1) cancels the state |− + + − − + + − ... , which therefore forms a Krylov subspace of dimension 1. We define the state |Φ m n for periodic boundary conditions as: where (w) m again marks that we repeat m times the sequence w. Hence, |Φ m n is therefore a state of size N = (2n + 4)m, and can be seen as several dipoles oscillating in a sea of pseudo-spins. Note that Ψ 1 n also belongs to this family when considering periodic boundary conditions. The dimensions and scaling of the Krylov spaces with total system sizes at fixed n can also be computed analytically using a transfer matrix approach, for both open and periodic boundary conditions. The dimension asymptotically scales as t m + where t + is the largest eigenvalue of the transfer matrix T (n) whose entries are given by We summarize in Table I  particular, for n = 1, we show that t + = 3 + 2 √ 2 and that the dimension of the Krylov space scales as In this example, the slow growth of the Krylov space cannot be understood from a simple quasi-particle picture. Indeed we prove in App. A 2 that there exists no p ∈ N * such that t p + is rational. The subspaces also always scale slower than 2 N as can be seen from Table I (see also App. A 2). Actually, 2 N matches the Krylov spaces scaling when n → +∞, keeping m fixed. In the rest of the paper, we focus on the n = 1 and n = 2 families, i.e., the families with the two slowest scalings. We will show in Secs. IV B and V that these two families have indeed ergodic statistics at zero and low disorders.

IV. LEVEL SPACING RATIO STATISTICS
Level spacing statistics are a convenient tool to determine whether a system is integrable or ergodic 4,47,80 . Random matrices without conservation laws, i.e., describing non-integrable models, have level-repulsion: the probability of having two eigenstates with the same energy is vanishing. Integrable models, on the other hand, are characterized by the presence of an extensive number of conserved quantities. Each sector then behaves as an independent random matrix and therefore there is no level repulsion between different sectors. Additionally, given a symmetry sector, directly studying the level spacing statistics requires unfolding the spectrum. Indeed, in order to obtain universal signatures, we are required to work with a uniform density of states. Several unfolding procedures exist, but finite-size effects may lead to different physical interpretations depending on the exact choice of method. [88][89][90][91] Instead, an efficient way to characterize quantitatively the level repulsion is to look at the level spacing ratio defined as follow 92 . The study of this quantity does not require flattening the density of states. Let {e n } be the ordered eigenspectrum of the Hamiltonian. We denote by r n the level spacing ratio Its probability distribution P (r) distinguishes between ergodic and integrable models. For an integrable model, P (r) is the Poisson distribution P Poi (r) = 1 (1+r) 2 , while for non-integrable systems, it depends on the symmetries of the Hamiltonian and is well-approximated by functionals of the form 93,94 : The real Hermitian Hamiltonians we consider fall into the Gaussian Orthogonal Ensemble (GOE) 88 with Z β = 8 27 and β = 1. In practice, it is more convenient to studỹ r n = min(r n , 1 r n ) (27) which is bounded between 0 and 1 and therefore has no heavy tails. For the classes we are interested in, P (r) = 2P (r)θ(1 − r). In the following, references to level ratio are references tor.
Finally, we remind the reader of the definition of the Kullback-Leibler (KL-)divergence: where p and q are the probability densities associated to the distributions P and Q. It trivially satisfies D KL (P, P ) = 0. The KL−divergence is asymmetric in (P , Q). It corresponds to the relative entropy from Q to P , that is to say the amount of additional information required to model P starting from the prior Q. Hence, when D KL (P num (r), P Poi (r)) < D KL (P num (r), P GOE (r)), the numerical distribution P num is better modelled by the Poisson distribution.

A. Level spacing ratio statistics of the full Hilbert space
Before we turn to the study of the individual Krylov subspaces themselves in Section IV B, we point out that it is crucial to decompose the symmetry resolved Hilbert space into its fractured components, in order to study any thermalization properties and transition.
We study a system of length N = 16 (L = 32 spins) with OBC in the symmetry sector p e = p o = c = 0 (see Eqs. (3) and (4)). The dimension of this symmetry sector is 4.8 × 10 6 , beyond the reach of full diagonalization. It fractures into approximately 2.5 × 10 5 Krylov subspaces, whose dimension varies from 1 to 12870. To emphasize the role of the dipole conservation, note that a system of 32 spins with only the two U (1) sublattice symmetries has already a symmetry sector (p e = 0, p o = 0) that includes 165M states.
We compute the exact full spectrum taking advantage of the decomposition into Krylov subspaces and identify whether each Krylov subspace has GOE or Poissonian statistics by computing the KL-divergences defined in Eq. (28), using the theoretical distributions as prior.
We fix W = 0.01 in order to avoid accidental degeneracies and limit finite-size effects, and average over 100 disorder realizations. We consider the Krylov subspace to have Poissonian statistics if D KL (P num (r), P Poi (r)) < 0.03, (29) and to have GOE statistics if Otherwise, we do not assign a label as the subspace is either afflicted by finite-size effects or presents signs of criticality. For comparison, the KL-divergences between the Poisson and GOE distributions are given by In practice, as we do not compare the distributions directly but histograms with 50 bins between 0 and 1, the effective divergence D KL (P Poi (r), P GOE (r)) is slightly reduced to 0.305 [D KL (P GOE (r), P Poi (r)) is almost unaffected]. Our choice of cut-off comes from the following observation: distributions that are maximally confusing with the KL-divergence verify as will be discussed in Sec. IV B. Choosing a threshold lower than 0.05 allows us to only select distributions that are convincingly Poissonian or GOE. For the sake of simplicity, we also focus only on Krylov spaces of dimension larger than 50 to minimize the number of samples to average over. This removes around 2.3 × 10 5 Krylov spaces associated to 1.8 × 10 6 states (41% of the symmetry sector), including 5 × 10 4 dark states, i.e., Krylov spaces consisting of a single state. The fractions of dark states decreases with system size. Fig. 1 summarizes our results and the nature of the Krylov spaces. Of the approximately 1.9 × 10 4 remaining spaces, a significant fraction present intermediate statistics (around 1.1 × 10 4 spaces, comprising 1.1 × 10 6 states, i.e., 23% of the symmetry sector). 4.2 × 10 3 Krylov spaces present clear Poissonian statistics and the remaining 2.6 × 10 3 have GOE statistics. They nonetheless represent a significant proportion of the total symmetry sector, approximately 7.9 × 10 5 states (16% of the total symmetry sector) and 9.9×10 5 states (20%) respectively.
We now turn towards the study of the level spacing ratio statistics in this symmetry sector, without resolving the Krylov spaces. As shown in Fig. 1d, the statistics are essentially undistiguishable from Poisson. In a given symmetry sector, the occupancies of each Krylov subspace act as an exponential number of additional good  . 1. a), b) and c) : histograms of the dimensions of the Krylov subspaces for N = 16 in the symmetry sector pe = po = c = 0 depending on their level spacing ratio distribution with W = 0.01. We only represent subspaces with dimension larger than 50. In a), we represent the Krylov subspaces that present ergodic level spacing ratio statistics, in b) Poissonian. c) provides the Krylov spaces that cannot be classified by our KL-divergence criteria defined in Eqs. (29) and (30). In each panel, we also give the number of Krylov spaces in that category, and the total dimension of these Krylov spaces. d): distribution of the level ratio for the same system if we do not resolve the Krylov spaces. Even if we discard the smallest Krylov subspaces, the level spacing ratio distribution is virtually undistiguishable from the Poisson distribution. quantum numbers. Therefore there is no apparent level repulsion. Theoretically, analytical formulas have been recently derived 95,96 to predict the distribution of the level ratios for matrices decomposing in several independent blocks. These studies computed the level spacing ratio distribution obtained from considering a small number (up to 12, but easily generalizable) of independent ergodic blocks as a single matrix. In Ref. 96, it was numerically shown that the mean level spacing ratio obtained from M ergodic blocks converges toward the Poissonian statistics approximately as M −2 . This means that, already for M = 12, the two average values differ only by 10 −3 . With the exponentially large number of blocks, and the additional scrambling induced by our Poissonian blocks, the precision required to differentiate our numerically obtained distribution from the true Poissonian distribution goes well-beyond any numerically achievable sampling. Indeed, we numerically obtain that the full numerical distribution P full , including all the Krylov subspaces, has a KL-divergence with respect to P Poi of

B. Level spacing ratio statistics in a single Krylov subspace
To study the effects of disorder, we focus on the families of Krylov subspaces defined in Sec. III. In particular, we specifically do not consider the largest Krylov subspace. Indeed, for OBC, the Hamiltonian restricted to this largest Krylov space is equivalent to a random XX model in a transverse field for all system sizes we considered. It is integrable and localizes at arbitrarily low disorder. Let the reduced energy of an eigenstate of energy E be with E min (E max ) the lowest (highest) energy of the reduced Hamiltonian in the Krylov subspace. In the rest of this article, we focus on states in the bulk of the spectrum with ε ∈ [0.4, 0.6].
We determine the level spacing ratio distribution by averaging over a large number of realizations, ranging from several thousands (for N = 8) down to 500 for the larger systems (for N = 30). We represent in Fig. 2a-c the KL divergence of the distribution of the level spacing ratio in the three families of Krylov subspaces defined in Secs III B and III C, using GOE and Poisson distributions as prior. For all families, at low disorder, we observe a quick convergence towards the GOE distribution of the level spacing ratio distribution, when increasing the system size N . The three families of Krylov subspaces appear indeed ergodic in the thermodynamic limit. Crossing of the KL-divergence universally occurs for D KL (P num (r), P Poi (r)) ≈ 0.05. This implies that we are maximally confused about which theoretical distribution better approximates the numerical one at this value of the KL-divergence. Thus, we take this crossing as a marker of the phase transition.
We first turn to the Krylov spaces generated by |Φ m 1 defined in Eq. (21), working with periodic boundary conditions due to the favorable scaling. As shown in Fig. 2a, for m > 2, the crossing point of the KL divergences shows very small finite-size effects at W 1 c ≈ 0.75 despite the small Hilbert spaces. We have performed the same analysis by studying the evolution of the mean level spacing ratio (see Appendix C 1). For the second family generated by |Φ m 2 , also with PBC, we observe similar results in Fig. 2b. Due to the faster growth of the Krylov subspaces, we are effectively limited to smaller systems plagued by stronger finite-size effects. For each m, we observe a transition from an ergodic phase to a localized phase, albeit at a significantly larger disorder strength, despite the similar Hilbert space dimensions and structures. The effective critical disorder strengths do not display any simple convergence behavior when increasing m, at least within the accessible system sizes. Finally, the family Ψ 1 n -studied with OBC due to the slower scaling -also exhibits signs of an MBL phase transition, as shown in Fig. 2c. Interestingly enough, the crossing point admits an approximately linear drift with increasing system sizes (see inset in Fig. 2c). Additionally, we observe in all Krylov subspaces that the critical disorder strength strongly depends on the relative energies of the eigenstates. Mobility edges 1,2,48 are therefore also present in these constrained systems. More details can be found in App. C 2.
Note that the three Krylov spaces considered here originate from the same initial Hamiltonian (up to boundary conditions) and therefore for a given W and N have the same disorder and hopping amplitudes in configuration space. Still, the corresponding critical values, as predicted by the level spacing ratio distributions, and scaling behavior are radically different. K( Φ 3 2 ) and K( Φ 4 1 ) both correspond to N = 24, K( Ψ 1 3 ) and K( Φ 2 2 ) to N = 16, and K( Ψ 1 4 ) and K( Φ 2 1 ) to N = 12 and yet admit different transition points. Conversely, the Krylov space dimension alone is also not a good indicator of the critical disorder: both K( Φ 2 2 ) (N = 16) and K( Ψ 1 4 ) (N = 20) have a dimension close to 2.2 × 10 4 . K( Φ 2 2 ) appears to localize at a larger disorder than K( Ψ 1 4 ) even though the disorder in the Fock basis is averaged over less sites.
The family K( Ψ 1 n ) shows a strong drift of the critical disorder towards larger values with increasing system sizes. This could be a sign of an absence of a transition for these subspaces in the thermodynamic limit. Paradoxically, this family also has a structure very close to an integrable one. Indeed, the Hamiltonian acting on the sea of pseudo-spins |↑ and |↓ reduces to a non-interacting XX Hamiltonian. The pair of dipoles breaks integrability  (31)). For very weak disorder and small Krylov space dimensions, the level spacing ratio distributions present strong finite-size effects due to the sparsity of the model and quasidegeneracies. Nonetheless, when increasing system sizes, all three families convincingly have ergodic statistics. For all systems, we observe an effective transition from GOE to Poisson statistics. The critical disorder is identified with the crossing points of the two divergences. Note the different effective critical disorder (Wc ≈ 0.75 for |Φ m 1 , Wc ≈ 1.15 for |Φ m 2 ) for the first two families with PBC. For the subspaces generated by Ψ 1 n (OBC), on the other hand, we observe an approximately linear shift of the effective critical disorder with increasing system sizes, as is shown in inset. Error bars (generally too small to be seen) are obtained through subsampling of our data.
by stitching together a set of triplets-the sea of pseudospins to the left, in-between, and to the right of the pair of dipoles-of integrable spaces. Yet, while the XX Hamiltonian is localized at arbitrarily low-disorder, with an effective critical disorder strength decreasing with system size, we observe the exact opposite for Ψ 1 n . The different behavior observed in our Krylov subspaces reinforces the need to distinguish between the Krylov subspaces if we want to study the MBL phase transition and the effect of disorder. In App. C 3, we show some additional numerical results showing the level spacing ratio statistics obtained when mixing the subspaces generated by Φ 3 2 and Φ 4 1 . We observe a significant difference between the distribution at low-disorder and P GOE , and a smoother crossover when studying the KL-divergences.

V. ENTANGLEMENT ENTROPY IN A CONSTRAINED MODEL
In the previous Section, we have seen that the different Krylov subspaces appear to undergo an MBL phase transition at different critical disorder strengths, according to their level spacing ratio distributions. We now turn to the study of the von-Neumann entanglement entropy (vNEE) of the many-body eigenstates as another complementary probe of this transition. For a subsystem A, the vNEE of the pure state |Ψ is given by: where Tr A marks the trace on the degrees of freedom not in A. We denote by H A (resp. H A ) the Hilbert space of ρ A (resp. ρ A ). In terms of the entanglement entropy, the MBL phase transition can be seen as a transition from thermal volume-law to an area-law [19][20][21][22] . In one dimension, the volume law is to be understood as with l A the number of sites (degrees of freedom) in A. s th takes the value log 2 in the thermal phase for conventional spin-1 2 systems. The strongly constrained model we study sees very irregular growth of the Hilbert space H A with subsystem size. Instead, we consider the entanglement entropy to be ergodic if it verifies: where the Page entropy 97 S Page is the average entanglement entropy of uniformly distributed random states. In the absence of symmetries or of Hilbert space fragmentation, the Page entropy S Page (A) is given by The Page entropy trivially satisfies the volume law as log m is roughly proportional to the number of degrees of freedom in A. Due to the presence of the multiple U (1) symmetries, we have to take into account the splitting of the wave functions down to submatrices in different symmetry subsectors. Correspondingly, the Hilbert space (or Krylov subspace) can be split into: where the subspaces H A,j and H A,j have dimension such that 98 dimH = j m j M j . The Page entropy is then given by The area-law remains here defined as We compute the entanglement entropy in the different Krylov spaces introduced in Secs. III B and III C, and average over all states with ε ∈ [0.4, 0.6] and over a large number of disorder realizations (see Sec. IV B). We work in the original spin basis (|0 , |1 ). We assume PBC for Φ 1 m and Φ 2 m , and OBC for Ψ 1 n due to the favorable scalings. In Fig. 3, we show the scaling of the vNEE S(l A ) as a function of the subsystem size, where A is the segment made of the l A consecutive spins [[1, l A ]] for different disorder values for Φ 5 1 , Φ 3 2 and Ψ 1 4 . The entropy we obtain therefore matches the one we would obtain in the pseudo-spins basis when l A is even. At low disorder values, the vNEE remains roughly proportional to the Page entropy, following the aforementioned volume law. The entropy varies only weakly with the disorder strength. At stronger disorder, we observe a crossover towards an area law where the entanglement remains (nearly) constant over several decades. This area law is typical of the predicted MBL phase and shows no signs of increasing again at larger scales.
The exact pattern followed by the vNEE depends on the Krylov subspaces, and can be very irregular. In particular, the family |Φ m 1 , for the cut we chose, has S(3l + 1) = S(3l + 2) = S(3l + 3) for l ≥ 1. It is not a consequence of any effective three-spin quasi-particles but a non-trivial interplay between the pair-hopping terms and the chosen starting state (see App. B).
Additionally, as can be seen in Fig. 3a, in the ergodic phase, the growth of the entanglement entropy from one plateau to the next alternates between large and small jumps. This irregular growth pattern comes from the lack of translation invariance at the single spin level in the starting generating state (while it remains invariant by translation of 12 spins). The dimension of the reduced density matrix grows faster when l A goes through a higher entropy jump. This irregular growth also affects the MBL phase. A larger growth of the reduced Hilbert space translates into more states connected by a pair-hopping term going through the entanglement cut. As the entanglement entropy at large disorder mainly arises from local resonant pairs, this structure leads to the observed alternating high and low plateaus. More details on the growth of the dimension of the reduced density matrix and the pairing structure can be found in App. B).
To pinpoint the transition, it is convenient to study the standard deviation of the entanglement entropy (typically at the midchain point) 20 . The transition point is taken to be at its maxima: the system can there be either in a thermal state with high volume-law entanglement or in a localized states with low entanglement. In Fig. 4, we show the standard deviation of the entanglement entropies obtained for all states with ε ∈ [0.4, 0.6] and for different disorder realisations for the Krylov subspace we considered. For all families, the larger the system, the more peaked the standard deviation is. For the family generated by |Φ m 1 , the peaks clearly concentrate around the critical disorder value W c ≈ 0.75. For |Φ m 2 there is no clear tendency emerging. Finally, for Ψ 1 n , the effective critical disorder values increase quasi-linearly with system size, preventing pinpointing any phase transition. The obtained values are in qualitative agreement with those obtained considering the level spacing ratio. Due to the limited number of sizes available in each family, we cannot perform a reliable scaling analysis.

VI. DISCUSSIONS AND CONCLUSIONS
In this paper, we have provided numerical evidence of a many-body localisation type transition within the ergodic Krylov subspaces of constrained models presenting a strong fragmentation of the Hilbert space. Due to the slow scaling of the Hilbert space dimensions, we have been able to study systems comprised of up to 60 spins using exact diagonalization. We observe the transition from an approximate linear scaling of the entanglement entropy at low-disorder to a clear area law over significantly larger scales than conventionally studied. We see no signs of a general breakdown of the manybody localization phenomenon in these large systems, though the Krylov spaces' dimensions remain comparable to other models which have been studied. Within the same constrained model, the different Krylov sub- spaces see a transition occuring at wildly different disorder strengths. This reinforces the importance of considering separately each Krylov space to study the localization properties of systems that see such a fragmentation of the full Hilbert space: without doing so, any sign of the transition will be blurred towards Poissonian statistics. More importantly, we see no significant correlations between effective critical disorder strengths and Krylov space dimensions or system sizes, as illustrated in Fig. 5.
The role of the structure of the Krylov space is therefore key to explaining the MBL phase transition in these models, and a detailed study is left for future works. The subspaces generated by the families |Φ m 1 appear to present a stable MBL phase transition in the thermodynamic limit, whether we consider the level spacing ratio distributions or the entropy properties. On the other hand, the subspaces generated by Ψ 1 n show an approx-imately linear scaling of the critical disorder strengths with system size, within the sizes and the numerical precision we have access to. It raises the question whether this subspace actually always thermalizes in the thermodynamic limit. This is especially remarkable given that this subspace and the action of the constrained Hamiltonian on it appear the closest to an effective integrable XX model given the presence of a single pair of dipoles in a sea of integrable spins. , it is harder to extract a tendency within the system size available, due to a strong finite size effect. Finally, for Ψ 1 n (c), we observe an approximately linear increase of the effective critical disorder with system sizes. In all cases, the predictions qualitatively agree with the ones obtained from the level spacing ratio statistics. Error bars are too small to be seen.  We discuss in this Appendix the simplification of the formulas given in Eqs. (18) and (19). We start with periodic boundary conditions, where the dimension of the Krylov subspace is given by: (A1) We reorganize the double summation introducing s = n ↑ + n ↓ , (A2) The bounds of the sum on n ↑ can be simplified as either one of the binomial coefficient is 0 for n ↑ < max(0, s−N ↑ ) or n ↑ > min(N ↑ , s). Namely, we get From here, application of the Chu-Vandermonde identity k j=0 m j  TABLE II. We summarize the properties of the Krylov space for the family generated by Ψ 1 n . The columns 3 to 5 (6 to 8) are for the OBC (PBC) case. The third and sixth columns list the dimensions of the Krylov subspaces according to Eqs. (18) and (19). The dimension of the Krylov spaces approximately grows as 2 N . The fourth and seventh columns show the connectivity of the Hamiltonian C (here defined as the ratio of the number of non-zero non-diagonal terms in the Hamiltonian over the Hilbert space dimension). Finally, the fifth and eighth columns list the dimensions dim ρ OBC 1 2 of the reduced density matrix for a cut exactly in the middle of the chain, i.e., separating the two +s in the generating state for different system sizes.
We now turn to open boundary conditions. The di-mension of the Krylov space (denoted here d OBC for convenience) is given by Using the Chu-Vandermonde identity and the fact that the term in the second sum is zero for s ≥ 2N L ↑ + 1, we obtain Similarly, we can instead introducet = 2N L ↑ + n ↑ + 1 to This leaves us with where we used Eq. (A10) a second time and obtain Eq. (20) in the main text. Note that as far as we know, this special identity for the triple sum of binomials is not registered in conventional tables. It can be generalized where we used, following Eq. (A11), the two changes of variables t x = 2X − x and t y = 2Y + y + 1, and applied Eq. (A10).

Krylov subspaces generated by |Φ m n
We now turn to the asymptotic dimension scaling of the Krylov spaces generated by |Φ m n . Let us first con-sider periodic boundary conditions for simplicity. ↑'s can move in between the two dipoles to their left (but not beyond), and similarly for ↓'s to their right. We denote x j = 0, ..., n and y j = 0, ..., n the number of pseudo-spins of the j th sequence of pseudo spins that have moved to their left and to their right. For example, we consider the Krylov subspace generated by Φ 3 2 , i.e., A typical configuration connected to this initial state looks like (A22) It satisfies (x 1 , y 1 ) = (0, 1) as only one of the ↓ pseudospins of the first subsequence of ↑↓↑↓ has moved the right of the first (leftmost) dipole −+, and none to the left of the last (rightmost) dipole +−. As can be straightforwardly observed, it also satisfies (x 2 , y 2 ) = (2, 1) and (x 3 , y 3 ) = (1, 2).
where the (n+1)×(n+1) transfer matrix T (n) has entries given by and x, y = 0, ..., n. Defining the matrices F x,y = f (x, y) and G 0 x,y = x+y x , we obtain the simple expression: Thus, the dimension of the Krylov subspace scales as (t (n) ) m with t (n) the largest eigenvalue of T (n) , as it corresponds to the translation by a single motif (↑↓) n − + + −. Using the relation N = (2n + 4)m, we readily obtain that the dimension of K(|Φ m n ) PBC scales as (t (n) ) N/(4+2n) = (t (n) ) L/(8+4n) .
Let us consider n = 1 as a concrete example. The matrix T (1) is given by Its eigenvalues are t (1) ± = 3 ± 2 √ 2, and therefore The dimension of the Hilbert space cannot be understood from a simple quasi-particle picture. Indeed, there exists no p ∈ N * such that (t (1) + 1 12 ) p is rational. The proof goes as follows: if such a p exists, then there existsp ∈ N * such that (T (1) )p has rational eigenvalues. Eigenvalues of (T (1) )p are roots of the polynomial det((T (1) )p − λI), which has integer coefficients, and a leading coefficient of 1. The eigenvalues are therefore real irrational integers, whose intersection with Q are integers only. They are given by The first parenthesis is trivially an integer, while the second term is of the form √ A 2 − 1, with A ∈ N and A > 1. This second term is therefore never an integer, and (T (1) )p can have neither integer nor rational eigenvalues.
For open boundary conditions, the Krylov subspace dimension can be obtained from a similar formula: The dimension of the OBC Krylov subspace therefore scales as in the periodic case.
In Tables III and IV, we summarize the dimensions of the Krylov subspace and of the reduced density matrix for a half-chain cut. In the Krylov subspace built from |Φ m 1 , the reduced density matrix follows an interesting simple pattern. As can be seen in Fig. 3a, each cut l = 3j + 1, l = 3j + 2 and l = 3j + 3 has the same entropy for j ≥ 1 and PBC, i.e.,

S(l
For OBC, this property is true already at j = 0. In this Appendix, we show that the reduced density matrices obtained for these cuts are actually identical. It can be proven by rewriting |Φ m 1 in terms of states combining three consecutive spins. We use the compact notation |ν 1 ν 2 ν 3 = |4ν 1 + 2ν 2 + ν 3 where ν j = 0 or 1. In this notation, |Φ m 1 can be written as |(4474) m . The relevant transformation rules induced by the pair-hopping terms now read |44 ↔ |30 and |47 ↔ |33 (B2) All other configurations of |0 , |3 , |4 and |7 are cancelled by the four fermions hopping terms and preserved by the transverse field. No state containing |1 , |2 , |5 , |6 and |8 is therefore connected to |Φ m 1 . The Krylov subspace only contains (a subset of the) states built from the triplets |0 , |3 , |4 and |7 .
For clarity, we start with the OBC case, and discuss the PBC case later. We will consider separately the case j = 0, j = 1 and j = 2 by direct inspection. Then we will address the generic j value. For any given starting state, we can write the single-site density matrix of the left most state explicitely as: where α, β are positive real coefficients and γ is a complex constant. Using the transformation rules in Eq. (B2), the configurations |44 can only be mapped to |30 or |43 . This means the configuration of the first triplet is either |3 or |4 . Once we fix the first spin, the next two are determined. Using Eq. (B3), we obtain the following expressions for the reduced density matrices at l A = 2 and l A = 3.
The coefficients of the density matrices, and therefore the vNEE, remain the same whether we cut after the first, second or third spin. Now, we turn to a cut through the second triplet (corresponding to j = 1 in Eq. (B1)). The of the reduced density matrix for a cut exactly in the middle of the chain.
(B6) Therefore, fixing the first 4 spins, i.e., the first triplet and the first spin of the second triplet again entirely determines the states of the second and third spin. Eq. (B1) is therefore also valid for j = 1. A similar reasoning can be applied to the third triplet and j = 2, with the sequences: To straightforwardly extend the results to the rest of the system, it is enough to consider all possible four triplets sequences in the Krylov subspace. There are only 44 such sequences (out of 2 12 = 4096 possible spin configuration and 4 4 = 256 combination of triplets), given in Tab. V. They can be obtained by brute force for m = 3 (m = 2 for PBC), and the limited propagation of the dipoles through the pseudo-spins ensures that no other configurations arise for larger systems. For all these sequences, fixing the first three triplets and the first spin of the fourth triplet is enough to determine the whole sequence. Eq. (B1) is therefore valid for any j.
For PBC, the same analysis can be performed, leading to the same property and pattern observed in Fig. 3a Additionally, the growth of the reduced Hilbert space from one three-site plateau to the next is very irregular, as shown in Fig. 3a. The alternating small and large jumps in entropy observed in the ergodic phase translates into alternating low and high entanglement plateaus in the MBL phase at strong disorder. Both phenomena can be explained by simple perturbative expansion arguments in the triplet language.
At strong disorder, the dominant energy terms are the disorder terms, and we can assume that the eigenstates are generally close to product states in the z spin-basis, and therefore in the triplet basis. Non-zero contributions to the entropy mainly come from local resonances, with two nearest neighbours triplet forming a pair due to the corresponding hopping term. In Table VI, we summarize how many pairs of states the hopping terms can generate depending on where they are applied on the Krylov spaces generated by |Φ m 1 .
Hopping term j ↔ j + 1 0 − 1 1 − 2 2 − 3 3 − 4 Number of pairs of states in Φ 1 for different values of m, taking PBC. The remaining terms can be obtained due to the invariance by translation of 12 sites. Due to the pattern of the generating state, we alternate between large and small numbers of connected states. This pattern explains the alternating large and small increase of entropies in the ergodic phase, and the alternating high and low entropy links in the MBL phase seen in Fig. 3a.
We treat as an example the case m = 1. With periodic boundary conditions, the Krylov subspace consists of only 6 states: |i ≡ |4474 , |ii ≡ |3074 , |iii ≡ |4334 , |iv ≡ |0473 , |v ≡ |0333 and |vi ≡ |0347 . The pairhopping terms linking the first two triplets only transform the state |i into |ii (and vice versa). Similarly, those connecting the 3rd and 4th triplets only transform |v into |vi . On the other hand, the pair-hopping terms connecting the 2nd and the 3rd triplet map |i into |iii and |iv into |vi . The alternating low and high number of pairs perfectly explain the observed entropy patterns. If this number is small, the average entropy in the MBL phase is lower as a limited number of local resonant pairs can exist. Conversely, at low disorder, the number of pairs reflect the number of connections in the configuration basis, that is to say the growth of the reduced Hilbert space when increasing system size. A lower number of pairs implies that the subspace grows less and therefore that the entropy increases less in the ergodic phase.

Appendix C: Additional numerical data
In this Appendix, we present briefly some additional numerical results mentioned in the main text.

Mean level spacing
We first turn towards the computation of the mean value of the energy level ratio. The mean level ratio is simply defined as the average of the distribution introduced in Eq. (27). It is a good indicator of the MBL phase transition as it crosses from 0.5307 for a GOE distribution to 0.38629 for Poisson distribution 92 . On the other hand, it generally shows a significant shift with system size and captures only partially the behavior of the distributions through the phase transition. In Fig. 6, we represent the mean level ratio for the three Krylov subspaces we consider. We observe in all cases, a crossover from GOE statistics to Poisson statistics. The behavior in each family is qualitatively different, with a sharper transition for |Φ m 1 , significant changes with system size for |Φ m 2 (compared to the other two Krylov subspaces), and a slower transition for Ψ 1 n . We also observe a significant shift of the transition point with system size for the family generated by Ψ 1 n . These results are consistent with those obtained by studying the KL-divergence in the main text. and Ψ 1 n for different values of m and n. We generically observe a crossover from the GOE value towards the Poisson value. The observed behavior is compatible with the results based on the KL-divergence in the main text.

Mobility edge
As was observed 1,2,20,48 in similar models, a MBL transition typically presents a mobility edge, that is to say that the critical disorder strength depends on the energy of the eigenstates. The same mobility edge can also be observed in our constrained model, as illustrated in Fig. 7. We compute the level spacing ratio distribution of states of normalized energies ε ∈ [ε t − 0.025, ε t + 0.025] for ε t varying from 0.075 to 0.925. We then estimate the critical disorder strengths as shown in Sec. IV B: it corresponds to the disorder strengths where the KLdivergences of the numerical distribution with the Poisson distribution or GOE one are equal. For all three spaces, we observe that the critical disorder strength significantly varies with the energy target, over a range and shape comparable to those observed in other models 48  for different values of m and n, as a function of the normalized energy. We compute the distribution of level ratio of states with normalized energies ε ∈ [εt − 0.025, εt + 0.025] for a range of εt. We estimate the critical disorder strength using the KL-divergence as shown in Sec. IV B. For all three families, we observe a significant variation of the transition point with energy level. Error bars are obtained through subsampling of our data.

Level spacing ratio for two mixed Krylov subspaces
Finally we study the level spacing ratio statistics obtained from mixing two Krylov subspaces seeing a transition at different disorder strengths. More precisely here, we consider the two Krylov subspaces generated by Φ 4 1 (of dimension 1154) and Φ 3 2 (of dimension 21873), taken with periodic boundary conditions with N = 24. We compute the level spacing ratio distributions by mixing the spectra obtained in the two spaces for the same disorder realisations. In Fig. 8, we show our results for eigenstates with ε ∈ [0.4, 0.6], averaging over 200 disorder realizations. The distribution of the level spacing ratio shows a significant departure from the GOE distri-bution at low disorder and converges towards the Poisson distribution at higher disorder. The behavior of the KLdivergence between the numerical distributions and our two reference distributions also shows a qualitative and quantitative difference with the distributions studied in Sec. IV B (see Fig. 2a-b). The difference at low disorder is better seen when studying D KL (P num , P Poi ), with a sat- 1 and Φ 3 2 , corresponding to N = 24. We observe a significant divergence from the ergodic GOE distribution at zero disoder, which converges slowly towards the Poisson distribution at higher disorder. b) KL-divergence of that distribution with the two reference distributions. The red lines mark the divergences between Pnum and PPoi). The discrepancy at low disorder is better seen in DKL(Pnum, PPoi) than in DKL(Pnum, PGOE). The phase transition is not as sharp than in Sec. IV B (see Fig. 2). Error bars are too small to be seen. uration value significantly lower than D KL (P GOE , P Poi ). We also observe a cross-over to the Poisson distribution, with a critical disorder strength in between the ones obtained in both subspaces. Unsurprisingly, the transition also appears much less sharp than in our study of the isolated Krylov spaces.