Many-body localization on finite generation fractal lattices

We study many-body localization in a hardcore boson model in the presence of random disorder on finite generation fractal lattices with different Hausdorff dimensions and different local lattice structures. In particular, we consider the Vicsek, T-shaped, Sierpinski gasket, and modified Koch-curve fractal lattices. In the single-particle case, these systems display Anderson localization for arbitrary disorder strength if they are large enough. In the many-body case, the systems available to exact diagonalization exhibit a transition between a delocalized and localized regime, visible in the spectral and entanglement properties of these systems. The position of this transition depends on the Hausdorff dimension of the given fractal, as well as on its local structure.


Introduction
The interplay between disorder and interactions in closed quantum many-body systems gives rise to transitions between ergodicity and many-body localization [1].Manybody localized systems are nonergodic and do not obey the eigenstate thermalization hypothesis.This phenomenon goes beyond the ground state properties and involves the physics of highly excited states.The transitions have been studied theoretically [2][3][4][5][6][7] and experimentally [8][9][10][11].The transition between ergodicity and many-body localization has also exposed rich structures in the entanglement entropy [12][13][14][15].
The study of many-body localization has been extensively carried out in onedimensional systems [2][3][4][5][6][7], and the prospect of many-body localization has also been studied in two-dimensional systems [10,16,17].The existence of a many-body localized phase in the thermodynamic limit is debated [7,[18][19][20].While the thermodynamic limit is important for condensed matter systems, the progress of quantum simulator technologies has made it possible to build small and highly controllable systems that can behave differently from the same system in the thermodynamic limit [21].For these few-particle systems, relevant for experiments with cold atoms, the existence of a manybody localized regime is well established experimentally and theoretically both in one dimension [2,5,8] and two dimensions [10,22].
Fractal lattices are aperiodic structures which exhibit fractional dimensions.Such lattices are playgrounds to explore properties, which depend on the dimension of space and how the particles can move.They provide an intermediate regime between integer dimensions (e.g. between one and two, or two and three).The fractals also display peculiar properties which are absent in integer dimensions, such as the lack of distinction between the bulk and edge.Several physical properties were studied in these systems [23][24][25][26][27][28][29][30][31][32].In particular, single-particle (Anderson) localization phenomena were investigated in fractal and bifractal setups [33][34][35][36].The possibility of manybody localization in systems other than periodic lattices, was in general rarely studied (an example is [37], where the formalism was developed for models defined on any graph).The studies of fractal lattices are also motivated by experimental developments to generate fractals [38,39], including the possibility to realize arbitrary lattices with Rydberg atoms [40] and the progress in techniques of constructing arbitrary lattice geometries and addressing single sites in optical lattices [41,42].
Here, we investigate models of hardcore bosons on finite generation fractal lattices in the presence of disorder and compare the results to one-and two-dimensional systems.The considered fractal lattices have different local structures and different Hausdorff dimensions, and we generally find that the systems transition from an ergodic to a many-body localized regime as the disorder strength is increased.Our results suggest that the strength at which the transition happens is influenced by several factors, but typically increases with increasing Hausdorff dimension.The transitions are not sharp phase transitions, as the considered lattices are far from the thermodynamic limit, and we stress that our results do not say whether the many-body localized phase exists in the thermodynamic limit or not.Extrapolating to the thermodynamic limit is challenging due to the rapid increase of the Hilbert space dimension and is not realistic with the exact diagonalization tools used here.Our results instead show how small, experimentally relevant systems are affected by disorder.We also briefly discuss the single-particle case, which shows Anderson localization.
The paper is organized as follows.We describe the fractal lattices and the hardcore boson model in section 2. In section 3, we study the single-particle localization in these systems.In section 4, we discuss the relation between the many-particle hardcore boson model and fermionic models.The spectral and entanglement properties of the manybody systems, showing a transition from ergodic to many-body localized behavior, are investigated in sections 5 and 6, respectively.We conclude the paper in section 7.

Fractal lattices and the model
In the following, we consider fractal lattices with different Hausdorff dimensions and different local lattice structures, and we start with a description of these lattices.The fractal lattices are shown in figure 1.
The Vicsek fractal lattice is constructed from a square as follows.We first divide the square into 3 × 3 squares.Then the four squares at the corners and the square in the middle are kept and the other squares are removed.Hence the remaining squares form a " × " shape.We repeat the procedure recursively for each of the five remaining squares to obtain different generations.Once the desired generation has been reached, a lattice site is placed on each of the remaining squares to obtain the fractal lattice.The Vicsek fractal lattice has a fractal dimension D = ln(5)/ ln(3) ≈ 1.465.A second generation lattice with N = 25 sites is shown in figure 1(a).
The construction of the T-shaped fractal lattice is the same as for the Vicsek fractal lattice, except that we keep the squares forming a "T" shape.The T-shaped fractal lattice has fractal dimension D = ln(5)/ ln(3) ≈ 1.465.A second generation lattice with N = 25 sites is shown in figure 1(b).
The Sierpinski gasket fractal lattice is constructed from an equilateral triangle by repeated removal of triangular subsets as follows.We begin with an equilateral triangle, which we subdivide into four smaller congruent equilateral triangles.Then we remove the central triangle.We repeat this step with each of the remaining smaller triangles to obtain different generations, and finally we obtain the fractal lattice by putting one lattice site on the center of each of the remaining triangles.The Sierpinski gasket fractal lattice has fractal dimension D = ln(3)/ ln(2) ≈ 1.585.A third generation lattice with N = 27 sites is shown in figure 1(c).
The modified Koch-curve fractal lattice is constructed as follows.We begin with a straight line and divide the line segment into three segments of equal length.We draw an equilateral triangle that has the middle segment from the previous step as its base.We place one lattice site at each of the vertices and we repeat the procedure to have different generations.The modified Koch-curve fractal lattice has a fractal dimension D = ln(5)/ ln(3) ≈ 1.465.A second generation lattice with N = 25 sites is shown in figure 1(d).
The fractals in figure 1(a), (b), and (d) all have the same fractal dimension but have different local lattice structures.For comparison, we also consider a one-dimensional chain and a two-dimensional square lattice with open boundary conditions.
On the fractal lattices, as well as on the one-and two-dimensional lattices, we define a hardcore boson model.Its Hamiltonian is where c i is the annihilation operator of a hardcore boson, acting on the ith lattice site, and The first sum is over the nearest neighbor sites of the lattice that are connected with bonds in figure 1 (or, in the case of one-and two-dimensional lattices, over nearest-neighbor sites).Note that each bond contributes two terms making the Hamiltonian Hermitian.We take h i ∈ [−h, h] to be random variables, which are uniformly distributed in the interval from −h to h, where h is the disorder strength.This model is motivated by recent experiments [43] in quasi-one-dimensional optical lattices.

Single-particle properties
Before delving into the physics of many-particle systems, we study the single-particle localization problem on the fractal lattices.That is, we consider the Hamiltonian (1), but fix the number of particles to be M = 1.In this case, the interaction enforced by the hardcore condition, as well as the statistics of the particles, have no effect, and the problem is equivalent to the single-particle free fermion problem.
In infinite one-and two-dimensional systems, arbitrarily weak disorder is enough to localize the wavefunctions [44].In finite systems, the localization is seen when the system size is large compared to the localization length.Here, we consider the fractals from section 2 to see if their behavior is similar.To probe the localization of the singleparticle wavefunctions, we study the level spacing statistics and inverse participation ratio.We average over 1000 disorder realizations.
Level spacing statistics.-In the ergodic phase of systems with time-reversal symmetry, the statistical distribution of the energy level spacings obeys the Wigner-Dyson surmise of the Gaussian orthogonal ensemble as the energy levels exhibit level repulsion.On the contrary the Poisson distribution is expected in the localized phase.Therefore to discriminate between these two behaviors we consider the ratio of consecutive level spacings [2,5,[45][46][47] defined as where the energies E n of the eigenstates are labeled in increasing order with the index n ∈ {1, 2, . . ., dimH} and dimH is the dimension of the Hilbert space (for M = 1, dimH = N ).The value r = ⟨r n ⟩ represents the average of rn over disorder realizations, and the value rn represents the average of r n over the middle 1/3 portion of the energy spectrum.The value of r approaches r = 4 − 2 √ 3 ≈ 0.536 in ergodic systems with time-reversal symmetry and r = 2 ln(2) − 1 ≈ 0.386 in localized systems [48].
The results are shown in figure 2. All the fractals seem to display a behavior qualitatively similar to the behavior of the one-dimensional system.For small systems and weak disorder, r has a value significantly higher than 0.386, but not necessarily close to 0.536, suggesting a delocalized regime.With increasing disorder strength, r gets close to 0.386, which suggests localization.As we increase the system size, the maximum of r decreases, and the entire curve gets closer and closer to 0.386, suggesting that in the thermodynamic limit, the system will be localized for any disorder strength.For the two-dimensional system, we observe a delocalized regime with r ≈ 0.536 at weak disorder.The position of the transition shifts toward lower h, however, as the system size increases, which is consistent with the theoretical result that in the infinite two-dimensional system localization occurs for an arbitrarily weak disorder.Inverse participation ratio.-For the nth eigenstate of the single-particle Hamiltonian, the inverse participation ratio is defined as where ψ i is the coefficient of the eigenstate at site i.The inverse participation ratio measures the degree of localization.If a wavefunction is spread evenly over the sites, we have P 2 (n) = 1/N , while if it is fully localized on a single site, P 2 (n) = 1 (see [45] for an example of how the inverse participation ratio is used in an Anderson localization problem).Thus, in the ergodic system, P 2 (n) is inversely proportional to the number of sites, while in the localized system it is constant.Similarly to r, we define P 2 = ⟨P 2 (n)⟩ as the value of P 2 (n) averaged over the middle 1/3 of the spectrum and over disorder realizations.
In figure 3, we plot ln(N P 2 ) as a function of ln(N ).For strong disorder, the dependence is seen to be almost linear with slope close to 1, which signifies localization.For weak disorder, the dependence becomes less linear, but otherwise roughly similar, except of two cases: the one-and two-dimensional systems at weak disorder.In these cases, for h = 0.1 (the weakest investigated disorder) we observe that ln(N P 2 ) is almost constant.
The results from figures 2 and 3 suggest that the studied fractals indeed behave as an interpolation between one and two dimensions, in the sense that they seem to display localization for arbitrarily weak disorder, provided the system size is large enough.Moreover, it seems that the fractals do not display the P 2 ∼ 1/N regime even at h = 0.1.This suggests that the obstructions to localization in small systems may be less pronounced in the fractal lattices.An alternative explanation is that the fractal structures are not uniform, so there is no reason for the wavefunction to spread evenly over the sites.We note that very pronounced inhomogenities were observed on the Cayley trees (finite Bethe lattice), which have a large fraction of edge sites.In such a case, it was possible for localized states to occur near the edges while delocalized states exist near the center of the lattice in finite systems [49,50].We do not rule out similar effects in our tree-like fractals (Vicsek, T-shaped), but we note that our systems are more complicated: there are more than two types of sites with different local environments.Thus, we focus on the quantities averaged over 1/3 of the spectrum.We also note that there are significant differences between these fractals and Cayley trees.

Jordan-Wigner mapping for many-particle systems
Since the M = 1 case is equivalent to the problem of free fermions, it is instructive to map our hardcore bosons to fermions using the Jordan-Wigner transformation.In order to do it, we have to number the lattice sites from 1 to N , which can also be seen as a mapping to a one-dimensional chain with hoppings that are generally non-local.The Jordan-Wigner transformation transforms the hopping −c † i c j + h.c. between sites i and j, with i < j, into Thus, the nearest-neighbor hoppings of bosons turn into nearest-neighbor hoppings of fermions, but lth-neighbor hoppings (l = i − j) turn into interaction terms involving up to l particles.The simplest case is the one-dimensional system, which transforms to a chain of free fermions.The two-dimensional case is more complicated, but it has a natural way of ordering the sites when mapping to a one-dimensional chain, which leads to interactions involving up to five particles for the 5 × 5 system considered in this work.For fractal lattices, there is no obvious way of ordering the sites (and e.g. for 25 sites there are more than 10 25 permutations).One example for a T-shaped fractal is shown in figure 4. It can be seen that the system maps into several free fermion chains of different lengths (the sites connected by nearest-neighbor hoppings in figure 4), connected by interaction terms involving several particles (more distant hoppings in figure 4).We expect a similarly complicated result for other fractals as well.Therefore, the Jordan-Wigner transformation simplifies the problem only in the one-dimensional case.The one-dimensional system, being equivalent to free fermions, is integrable for any value of the disorder and is Anderson localized for the disorder strong enough that the localization length is smaller than the system size (see figure 3(a)).For the two-dimensional case we can rely on previous experiments and calculations, showing that finite systems of hardcore bosons exhibit a transition between ergodic and manybody localized regimes [10,16,17].The fractal lattices will be studied numerically in the following sections.

Many-body spectral properties
We now use exact diagonalization to determine the properties of model ( 1) on the finite fractal lattices for M > 1. Due to the rapid increase of the Hilbert space dimension with system size, we limit the study to the system sizes listed in table 1.We choose the particle number in such a way that the filling factor M/N is nearly the same in all the lattices.We differentiate between thermal and many-body localized behaviors by studying the level spacing statistics and the many-body mobility edge with tools from random matrix theory.
Level spacing statistics.-Thelevel spacing statistics are investigated using the parameter r, defined in section 3.In figure 5(a), we plot r as a function of the disorder strength h.The value of r is obtained by averaging over the middle 1/3 of the states in 250 disorder realizations (except from the Koch lattice, for which we study 2500 realizations, because for this case we consider M = 3 for which the Hilbert space is much smaller).
For the one-dimensional lattice, we observe that r is close to 0.386 for all h.This is related to the fact that the one-dimensional case is equivalent to a free fermion system The adjacent gap ratio r as a function of the disorder strength h for the considered lattices in the many-particle case.The dotted horizontal black lines at r ≈ 0.536 and at r ≈ 0.386 are the expected values for ergodic and many-body localized systems, respectively.For all the cases, there is a transition from ergodic behavior to many-body localized behavior as the disorder strength is increased.(b) The adjacent gap ratio r as a function of the disorder strength h for the same systems as in (a), but with one particle only (M = 1).
and thus integrable.For all the other lattices, we find a transition from ergodic behavior to a many-body localized behavior.This confirms that the fractals (even the tree-like ones) cannot be mapped to a free fermion model.The three lattices with Hausdorff dimension 1.465 all have the transition at comparable disorder strengths.The transition happens at a slightly larger disorder strength for the lattice with Hausdorff dimension 1.585, and for the two-dimensional lattice the transition happens at even stronger disorder.The fractal lattices hence show features intermediate between the one-dimensional and the two-dimensional case, and the Hausdorff dimension seems to be an important factor regarding the location of the transition.
A comparison between the single-and many-body results can be performed by replotting all the blue curves from figure 2 in one picture, which we do in figure 5(b).We note that the many-body results display a clearer transition between a localized and delocalized regime, as r for fractal lattices in figure 5(b) never reaches the Wigner-Dyson value.Also, the Hausdorff dimension seems somewhat less important in the single-particle than in the many-particle case -note the similarity between the r versus h curves for the Koch fractal (D ≈ 1.465) and Sierpinski triangle (D ≈ 1.585), as well as between the one-dimensional case and the Vicsek and T-shaped fractals (D ≈ 1.465), in figure 5(b).The influence of the Hausdorff dimension may be obscured by the variations in the system size (N = 20 for Koch fractal and N = 27 for Sierpinski lattice -the latter has 1.35 times more sites).
For a given energy E, the corresponding energy density is defined as where E min (E max ) is the lowest (highest) energy eigenvalue in the energy spectrum.We plot the energy density resolved gap spacings r(ϵ) (i.e.r n averaged over disorder realizations as well as over the eigenstates falling into a certain small window of ϵ) in the ϵ versus h plane in figure 6.For both the fractal lattices and the two-dimensional lattice a D-shape appears in the plots.In the region with the mobility edge, there are localized states at low and high ϵ and thermal states for intermediate ϵ.Also in this plot it is again seen that the transition from the thermal to the many-body localized behavior happens at lower disorder strength for the models on the fractal lattices than for the model on the two-dimensional lattice.

Many-body entanglement entropy
The entanglement entropy is another quantity to identify the transition between ergodicity and many-body localization [53,55].The variance of the entanglement entropy has a peak at the transition point.We study this quantity for the considered fractal lattices.(d,i) the Sierpinski gasket fractal lattice, and (e,j) the two-dimensional square lattice.We take two choices of the subsystems A and B for each case as shown in the insets, and we find similar behaviors for both.The peaks in var(S) show the transition points from thermal behavior at weak disorder to the many-body localized behavior at strong disorder.The mean entropy ⟨ S⟩ is high in the thermal region and low in the many-body localized region.
1.5 1.6 1.7 1.8 1.9 2.0 3 We partition the system into two parts, which we denote A and B. The von Neumann entanglement entropy of part A of the state |Ψ⟩ is where is the reduced density matrix of part A and the traces are over part A and B, respectively.We consider 20 energy eigenstates closest to the middle of the spectrum ϵ = 0.5.The mean entropy ⟨ S⟩ is obtained by averaging ( 6) over all these states and all the disorder realizations.Because we focus only on a small number of states, we can obtain them using the shift-invert approach, which is faster than full matrix diagonalization and allows us to study a larger number of disorder realizations: 10000 for the Koch lattice and 2500 for the other lattices.
The variance var( S) of the entanglement entropy is computed in the following way: first, we obtain the average S of the entropy over 20 states in the middle of the energy spectrum in each disorder realization, and then compute the variance var( S) of that average among the disorder realization.We plot the variance var( S) and the mean ⟨ S⟩ of the entanglement entropy as a function of the disorder strength h in figure 7 for all the lattices.We note that the choices of A and B are not unique on fractal lattices and hence we make two different choices in each case for comparison.The values of var( S) are small both deep into the ergodic region and deep into the many-body localized region but are large in the vicinity of the transition point.The mean value of the entropy ⟨ S⟩ is high in the thermal region and low in the many-body localized region.
For a quantitative comparison of the transition in different lattices, we study the position of the maximum of var( S).To estimate the error bar on the maximum, we apply the following procedure.We divide the disorder realization into blocks of size 100.For each block we compute the variance var( S) block as a function of h.Then, we find the value h max corresponding to the maximum of var( S) block .Finally, we average h max over all the blocks, obtaining the mean ⟨h max ⟩ and the standard deviation σ max of the mean.We use max{σ max , 0.1} as an error bar, because 0.1 is the separation between the data points in the first axes of the plots in figure 7 (note that this separation is 0.5 in most parts of the plots, but we use a finer grid near the maxima of var( S)).
The resulting transition points ⟨h max ⟩ are shown in figure 8 as a function of the Hausdorff dimension.Both choices of bipartition yield the same results, and thus we present the plot only for one of them.A probable reason is the correlation between the results, caused by the fact that we use the same disorder realizations for both bipartitions (i.e. for a given disorder realization, we diagonalize the Hamiltonian and use its eigenstates to compute entropies for both bipartitions).
The transition points ⟨h max ⟩ from figure 8 fulfil the following condition: if we take two lattices with fractal dimensions D 1 and D 2 and transition points ⟨h max ⟩ 1 and ⟨h max ⟩ 2 , respectively, then if D 1 < D 2 , then ⟨h max ⟩ 1 < ⟨h max ⟩ 2 , suggesting the importance of Hausdorff dimension in the description of many-body localization in fractals.On the other hand, if D 1 = D 2 , it is not necessarily the case that ⟨h max ⟩ 1 = ⟨h max ⟩ 2 , suggesting the importance of other factors, such as local lattice structure.

Conclusions
We have studied the transition between ergodicity and many-body localization in a model of hardcore bosons in the presence of disorder on finite fractal lattices with Hausdorff dimensions between one and two and different local lattice structures.In particular, we have studied models on the Vicsek fractal lattice, the T-shaped fractal lattice, the Sierpinski gasket fractal lattice, and the modified Koch-curve fractal lattice.We considered a one-dimensional chain and a two-dimensional square lattice for comparison.We also compared the results to the single-particle case, where much larger system sizes are computationally available.
For the single-particle case, the results suggest Anderson localization for all disorder values, as is known for one and two dimensions.In the many-body case, a transition from an ergodic regime to a many-body localized regime was found on all the fractal lattices by studying the level spacing statistics, the many-body mobility edge, and the von Neumann entropy of highly excited states.The analysis of the variance of the entanglement entropy suggests that the disorder strength at the transition point generally increases with the Hausdorff dimension of the lattice, but is also affected by the local lattice structure.
The fractal lattices studied here provide further examples of systems beyond one dimension, which exhibit a many-body localized regime for finite size.At the same time, the number of sites scales more slowly with linear system size than in two dimensions.The fractal lattices considered here are bigger in the sense of linear size than twodimensional lattices with the same number of sites.
In recent years, there has been developments in the field of many-body localization in two dimensions.An avalanche theory for many-body localization [56][57][58] has been proposed, within which the many-body localization in two dimensions has been analyzed both analytically and numerically.It is concluded that a transition exists between ergodicity and many-body localization in a finite two-dimensional system.The disorder strength, at which the transition happens, increases, however, indefinitely with increasing system size.Hence we expect the transition disorder strength of our results in two-dimensions to be a growing function of the system size.Whether similar conclusions hold for fractal lattices remains an interesting open question.

Acknowledgments
This work has been supported by Danmarks Frie Forskningsfond under Grant No. 8049-00074B and by Carlsbergfondet under Grant No. CF20-0658.S.M. thanks Weizmann Institute of Science, Israel Deans fellowship through Feinberg Graduate School for financial support.

Figure 1 .
Figure 1.We consider different fractal lattices, namely (a) the Vicsek fractal lattice, (b) the T-shaped fractal lattice, (c) the Sierpinski gasket fractal lattice, and (d) a modified Koch-curve fractal lattice.To the left we show the first generation of the fractals, and the lattice sites are illustrated as grey circles.To the right we show either the second or the third generation of the lattices, and the allowed nearest-neighbor hopping terms are shown with black, solid lines.Note that in (d) we only keep those nearest-neighbor hoppings, which maintain the self-similarity of the fractal lattice.The number of lattice sites is (a-b) N = 25, (c) N = 27, and (d) N = 20, respectively.

Figure 2 .
Figure 2. The average level spacing ratio r for the single-particle systems as a function of disorder strength h.The dashed horizontal lines are located at r ≈ 0.386 and r ≈ 0.536, corresponding to Poisson and Wigner-Dyson statistics, respectively.

Figure 3 .
Figure3.The average inverse participation ratio for the single-particle systems, for different disorder strengths, as a function of system size.

Figure 4 .
Figure 4.One of the possible ways of folding the T-shaped fractal (figure 1(b))into a one-dimensional chain, with black lines denoting the hoppings.After the Jordan-Wigner transformation, the nearest-neighbor hoppings become single-particle fermionic hoppings, while the more distant hoppings transform into hoppings with phases depending on the occupation of all intermediate sites.

Figure 5 .
Figure 5. (a) The adjacent gap ratio r as a function of the disorder strength h for the considered lattices in the many-particle case.The dotted horizontal black lines at r ≈ 0.536 and at r ≈ 0.386 are the expected values for ergodic and many-body localized systems, respectively.For all the cases, there is a transition from ergodic behavior to many-body localized behavior as the disorder strength is increased.(b) The adjacent gap ratio r as a function of the disorder strength h for the same systems as in (a), but with one particle only (M = 1).

Figure 6 .
Figure 6.The energy density resolved adjacent gap ratios r(ϵ) as a function of the disorder strength h for (a) the Vicsek fractal lattice, (b) the T-shaped fractal lattice, (c) the modified Koch-curve fractal lattice, (d) the Sierpinski gasket fractal lattice, and (e) the two-dimensional square lattice.

Figure 7 .
Figure 7.The variance var(S) and mean ⟨ S⟩ of the entanglement entropy as a function of the disorder strength h for (a,f) the Vicsek fractal lattice, (b,g) the Tshaped fractal lattice, (c,h) the modified Koch-curve fractal lattice, (d,i) the Sierpinski gasket fractal lattice, and (e,j) the two-dimensional square lattice.We take two choices of the subsystems A and B for each case as shown in the insets, and we find similar behaviors for both.The peaks in var(S) show the transition points from thermal behavior at weak disorder to the many-body localized behavior at strong disorder.The mean entropy ⟨ S⟩ is high in the thermal region and low in the many-body localized region.

Figure 8 .
Figure 8.The position of the maximum of the variance of the entanglement entropy in the studied lattices plotted as a function of the Hausdorff dimension.

Table 1 .
The Hausdorff dimension D, the number of lattice sites N , the number of particles M , the lattice filling factor M/N , and the corresponding Hilbert space dimension dimH for the considered many-body systems.