Coexistence of localized and extended phases: Many-body localization in a harmonic trap

We show that the presence of a harmonic trap may in itself lead to many-body localization for cold atoms confined in that trap in a quasi-one-dimensional geometry. Specifically, the coexistence of delocalized phase in the center of the trap with localized region closer to the edges is predicted with the borderline dependent on the curvature of the trap. The phenomenon, similar in its origin to Stark localization, should be directly observed with cold atomic species. We discuss both the spinless and the spinful fermions, for the latter we address Stark localization at the same time as it has not been analyzed up till now.

We show that the presence of a harmonic trap may in itself lead to many-body localization for cold atoms confined in that trap in a quasi-one-dimensional geometry. Specifically, the coexistence of delocalized phase in the center of the trap with localized region closer to the edges is predicted with the borderline dependent on the curvature of the trap. The phenomenon, similar in its origin to Stark localization, should be directly observed with cold atomic species. We discuss both the spinless and the spinful fermions, for the latter we address Stark localization at the same time as it has not been analyzed up till now.
For a long time, it has been believed that many-body systems tend to thermalize as expressed by eigenstate thermalization hypothesis [1,2]. The many-body localization (MBL) phenomenon (for reviews see [3][4][5][6])) is a direct counterexample -for a sufficiently strong disorder, the system preserves the memory of its initial state. However, recent results on the existence of many-body quantum scar states [7,8], Hilbert-space fragmentation [9][10][11], and lack-of-thermalization in gauge theories [12][13][14] show strong non-ergodic behaviors even in the absence of disorder. Another example considers Stark localizationwhere the presence of a static electric field resulting in a tilt in the many body system may lead to localization [15,16].
It seems that the more many-body physics is explored the less ergodic the many-body dynamics turns out to be. The present work provides another example of such a situation. While currently there exits a vivid debate about the very existence of MBL in the thermodynamic limit [17][18][19][20], we consider finite system sizes only. Such systems are directly amenable to experimental studies in cold atomic settings [21][22][23][24]. We consider one-dimensional (1D) chains with chemical potentials quadratically dependent on position. Such a situation is quite common in quasi-1D situations realized in optical lattices [25,26], where a tight confinement in directions perpendicular to a chosen one is due to illumination by strong laser beams with gaussian transverse profiles. Those profiles may be well approximated as a harmonic trap along the considered direction [25]. Thus we shall consider models with Hamiltonians being where A is the curvature of the harmonic trap and l is the site index (we assume unit spacing between sites of the chain). H 0 is the model Hamiltonian considered, which may represent the Heisenberg chain (equivalent to interacting spinless fermions), bosons represented by In the absence of external curvature A = 0 the crossover to localized regime occurs for W ≈ 2.2 for such a small system. Above A = 0.3, system seems localized for arbitrary W . Right: r as a function of scaled energy and curvature A at disorder strength W = 0.5. The blue and red lines are contours for r ≈ 0.46 and 0.42 respectively and serve as the guide to the eyes for distinguishing localized and ergodic phases. Bose-Hubbard model, or spinful fermions with Hubbard Hamiltonian. In contrast to the study of [16], where small quadratic potential has been considered on top of the dominant uniform linear potential, we consider the effect of harmonic trap alone, which as we shall show acts differently in the different parts of the system. Heisenberg spins or spinless fermions.-As the simplest possible model, first we consider the Heisenberg chain, where H 0 becomes with S l 's being spin-1/2 operators and h l is a diagonal disorder (a magnetic field along z-axis) drawn from random uniform distribution in [−W, W ] interval. We set J = 1 to be the unit of energy. The harmonic trapping potential in this case is given by H trap = Site-observed spin dynamics, as measured by the local expectation values 2 × S z x for L = 50 Heisenberg chain with no disorder. For small A the system "thermalizes" and initial spin-wave configuration is destroyed by interactions. For larger curvature, one clearly observes the coexistence of delocalized (in the center) and localized regions (at the edges). The black dashed lines give the border of localization as given by Stark localization prediction with Fc ≈ 2 [15]. chain of spinless fermions via Jordan-Wigner transformation. A typical random matrix theory (RMT) based measure is the gap ratio defined as a minimum of the ratio of consecutive level spacings, r n = min{ sn+1 sn , sn sn+1 } with s n = E n+1 − E n . The mean gap ratio is r ≈ 0.53 for delocalized system, well described by Gaussian orthogonal ensemble (GOE), while r ≈ 0.38 for Poisson spectra characteristic for integrable, localized cases [29]. We find that for sufficiently large curvature A the spectrum shows Poissonian features regardless of the disorder amplitude (see Fig. 1 for L = 16). The figure resembles that observed for Stark localization [15].
To get insights into the physics observed, let us consider the time dynamics. We prepare the chain in the separable state with every second spin being up and down respectively as |↑, ↓, ↑, .., ↓ and observe whether this spinwave arrangement is preserved in time evolution. For small disorder and small A the system thermalizes (upper row in Fig. 2). For larger A different picture emergeswhile at the center of the chain delocalization still occurs, at a sufficient distance from it we observe preservation of the initial spin texture -i.e. localization.
One can easily explain this phenomenon. For a given distance l 0 from the center the local static field can be expressed as F = ∂ ∂l0 A 2 l 2 0 = l 0 A. If this local field exceeds the border of Stark localization [15,16] -the part of the system localizes, while the region close to the center remains extended. Therefore, unlike the usual Stark localization, one can always find localized regions for any finite values of A for large enough systems under harmonic trapping potential. The dashed lines in Fig. 2 give the Stark localization border, as predicted in [15] to be F ≈ 2, which nicely fits numerical data obtained using time-dependent variational principle (TDVP) algorithm using matrix product states (MPS) ansatz [30][31][32][33] (see [34] for details). Instead of spin profiles, one may look at the entanglement entropy growth in time (see Fig. 3). We see that entropy grows rapidly in the central region, while remaining low in the localized parts. Quite surprisingly, the entanglement entropy can also show logarithmic growth in time, even in the delocalized regions when the effect of trapping potential becomes strong. For example, in case of A = 0.2, entropy the central bond grows rapidly in time and approaches the maximum allowed value by the MPS ansatz (ln(512) in this case, see [34] for details). On the other hand, the entropy in the bonds 30 and 32 shows a logarithmic growth, despite being on the delocalized side of the system. The entanglement growth in case of A = 0.4 is greatly modified by the harmonic trap even in the central bond. This is indeed very unusual dynamics where delocalized parts behave as systems showing MBL in terms of entropy growth. Spinful fermions: Stark localization.-Since the interactions between spinless fermions are hard to realize experimentally, we consider the spinful case, represented by the Hubbard model, as in experiments [21,22]. The curvature free part H 0 is with l ∈ [−L/2, L/2]. As before, we set J = 1 to be the unit of energy and consider U = 1 throughout the paper unless otherwise stated. Let us first consider the Stark localization problem under linear potential as it was only addressed for spinless fermions [15,16] until now. Thus, we add to the Hamiltonian a tilt term F l l(n l↑ +n l↓ ) and analyze the gap ratio statistics. The corresponding r statistics is shown in Fig. 4(a) for L = 12 quarter-filled chain (i.e., N ↑ = N ↓ = L/4). To obtain this plot we break the SU(2) symmetry of the Hamiltonian by adding a local magnetic field to the Hamiltonian via the term H break = B(n L/2↑ − n L/2↓ ) with B = 0.5 following the prescription and discussion in [35]. We observe that, in comparison to spinless fermions [15], the crossover seems quite broad, possibly due to the small system size taken. To get more precise critical value of F , we consider time evolution of staggered density-wave state |↑, 0, ↓, 0, ↑, 0, ↓, ... for larger system-sizes and measure the density correlation C(t) = D l (n l (t) − ρ) (n l (0) − ρ) , wheren l (t) = n l↑ +n l↓ , ρ is the average number of particles per site, and the constant D is chosen so that C(0) = 1. The illustration of such time dynamics for L = 32 sites system is shown in Fig. 5. We observe that for the high field value, Top row and bottom left panel: Time dynamics of spinful fermions for the initial staggered density-wave state without any disorder for three static field amplitude, F , across the crossover. The localization is complete for F = 3, while for intermediate field values a partial localization is observed with some fraction of particles accumulating near the bottom of the effective potential. Bottom right panel: Time evolved density correlator C(t) at large times as a function of F . We discard 4 sites from the boundaries to minimize their effects. We take the inflection point of the curve at Fc ≈ 2.8 as a critical field amplitude. All data are for L = 32 obtained using TDVP algorithm. e.g., F = 3, the localization is almost complete (bottom left panel). On the other hand, at lower fields, we observe a partial localization (as revealed also by standard local densitiesn l (t)), see top row of Fig. 5). The bottom right panel of Fig. 5 shows large time values of the density correlator C(t) and its derivative dC(t) dt . We approximate critical F c from the inflection point of C(t), as obtained from the maximal value of its derivative, to be F c ≈ 2.8 for disorderless scenario. Spinful fermions: Localization under harmonic trap.-Having established the estimate for the critical field amplitude corresponding to the crossover to localized phase, we may turn again to the harmonic confinement case. Thus we again consider (1) now for spinful fermions (3), with H trap = A 2 L/2 l=−L/2 l 2 (n l↑ +n l↓ ). The right panel of Fig. 4 shows level spacing statistics for quarter-filled L = 12 chain in such case and Fig. 6 depicts the time dynamics of the density profile for different curvatures of the harmonic potential with the staggered density-wave state being the initial one. We observe that, as in the spinless fermions case, while in the center of the trap apparent fast "thermalization" occurs, closer to the edges the effective electric field coming from the curvature of the trap leads to localization. The dashed lines give the estimate of the threshold assuming F = l 0 A condition with F c = 2.8 from the previous analysis.
We may also analyze the time evolved entanglement entropy in different regions. While in the center of the trap the entropy grows linearly with time and soon saturates due to insufficient bond dimension of the MPS ansatz rendering the results in this region not accurate, the entropy beyond the localization boundary given by l 0 > F/A (and symmetrically for negative l 0 ) seems to grow logarithmically providing a further evidence for many-body localization in the outer regions.
The corresponding spin dynamics is also interesting. As in the standard MBL case, we observe a subdiffusive decay of initial spin configuration for the initial staggered density-wave state, characteristic of the remaining SU(2) symmetry of the problem [36][37][38]. Visualization of this effect as well as the plots for larger system size L = 64 may be found in [34]. Conclusions.-We have shown that many-body localization behavior can be observed in the presence of the harmonic trap and in the absence of the disorder. The effect is due to a local static field that induces, for sufficient curvatures, Stark localization as recently shown for spinless fermions [15,16] and announced in the spinful case [39]. Since the effect has a lower bound on the curvature, the central region of the trap always remains delocalized. Thus a harmonic trap makes a possible realization of a very interesting situation -coexistence of delocalized and MBL phases in a single system. Let us note that these regions do not suffer from avalanche effects [40] at least on the experimental time scales considered. Finally, let us mention that the harmonic trap may play some role in all the experiments on MBL up to date (e.g. [21,22]) as the residual trap due to Gaussian beam profiles is most probably present there. While in these experiments disorder induced effects play a dominant role the residual harmonic-like trap may affect the details of the time dynamics for large systems. This aspect is a subject of a current research. Finally, let us note that the effect is not limited to harmonic trap but can be generalized to arbitrary potentials with non-vanishing first order derivatives. While in the main text we have shown the localization in a harmonic trap for L = 50, it is fruitful to consider smaller system sizes. Fig. 1 shows the effect of the harmonic trap for different curvatures A. For too small curvature no localization effect is observed. With increasing A the system starts to localized at the edges, when the local electric field F = lA (where l ∈ [−L/2, L/2] is the site number) exceeds the Stark localization value F c ≈ 2 [1]. In Fig. 2 we show the time-evolved pattern of entanglement entropy for the same. The entropy grows rapidly in the central region which is apparently delocalized, and grows slowly in the localized sites close to the edges of the system. For smaller values of A (e.g., 0.16), the entropy follows 'volume-law' growth with the size of the bipartition in the time-evolved asymptotic states, while for larger values of A it is concentrated only in the central region.

Spinful fermions
In the main text, we have considered Stark localization in case of spinful fermions, as described by the Hubbard Hamiltonian in presence of the tilt term F l l(n l↑ + n l↓ ). We have analyzed the system under quarter filling with equal number of spin up and down particles, and only examined the charge (or density) sector during the time evolution via local densitiesn l (t) and the density correlator C(t). We now also consider the spin sector via the local spin operator, 2 × S z l =n l↑ −n l↓ , of the same system starting from the same staggered densitywave state |↑, 0, ↓, 0, ↑, 0, ↓, ... . Fig. 3  32 chain at F = 3 where the localization of the charge sector has been shown to be completed. However, we can easily appreciate the subdiffusion of spin sector (c.f. [2][3][4]) where the memory of the initial spin configuration starts to wear off after t ≈ 200. Let us now move to spinful fermions in the harmonic trap, H trap = A 2 L/2 l=−L/2 l 2 (n l↑ +n l↓ ). We immediately go to larger system sizes and present data for L = 64. The supplement allows to present them in a full glory. Fig. 4 shows the both the density profile (top row) as well as the spin profile (bottom row) for two different values of curvature A with the initial state being |↑, 0, ↓, 0, ↑, 0, ↓, ... . As in the case of L = 32 reported in the main text, the predicted borders between the localized and the delocalized regions obtained from the critical field strength F c ≈ 2.8 match very well with the density plots for L = 64 too. On the other hand, a slight spreading of the delocalized region in the spin sector at later times is an indicator of the subdiffusive spin dynamics [2][3][4]. Over all the qualitative features of MBL in a harmonic trap are very similar to that in Stark localization as well as in the standard disorder induced MBL.

Details about TDVP simulations
To simulate the time dynamics, mentioned in the paper, we employ time-dependent variational principle (TDVP) [5][6][7][8] using matrix-product states (MPS) ansatz with open boundary condition. Specifically, we use a hybrid variation of the TDVP scheme mentioned in [8][9][10][11], where we first use two-site version of TDVP to dynamically grow the bond dimension upto an allowed value, say χ max . When the bond dimension in the MPS is saturated to the value χ max , we move to the one-site version to avoid any error due to SVD truncation. In the paper, the results are produced with χ max = 512, so that the maximum allowed value of the entanglement entropy at any given bond in the bulk of the system is ln χ max = ln(512).
In the delocalized center the entropy grows fast, and therefore, the simulations in this region may not be accurate. However, as we move towards the boundaries of the system, the results with χ max = 512 becomes 'exact', even within the delocalized region. This we confirm by performing the same simulations with χ max = 384.