Ergodicity breaking with long range cavity induced quasiperiodic interactions
Piotr Kubala, Piotr Sierant, Giovanna Morigi, Jakub Zakrzewski
aa r X i v : . [ c ond - m a t . d i s - nn ] D ec Ergodicity breaking with long range cavity induced quasiperiodic interactions
Piotr Kubala ∗ Institute of Theoretical Physics, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland
Piotr Sierant
The Abdus Salam International Center for Theoretical Physics, Strada Costiera 11, 34151, Trieste, Italy andInstitute of Theoretical Physics, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland
Giovanna Morigi
Theoretical Physics, Saarland University, Campus E2.6, D–66123 Saarbrücken, Germany
Jakub Zakrzewski
Institute of Theoretical Physics, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland andMark Kac Complex Systems Research Center, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland. † Many-body localization (MBL) behavior is analyzed in an extended Bose-Hubbard model withquasiperiodic infinite-range interactions. No additional disorder is present. Examining level statisticsand entanglement entropy of eigenstates we show that a significant fraction of eigenstates of thesystem is localized in the presence of strong interactions. In spite of this, our results suggest that thesystem becomes ergodic in the standard thermodynamic limit in which the energy of the system isextensive. At the same time, the MBL regime seems to be stable if one allows for a super-extensivescaling of the energy. We show that our findings can be experimentally verified by studies of timedynamics in many-body cavity quantum electrodynamics setups. The “quench spectroscopy” is aparticularly effective tool that allows us to systematically study energy dependence of time dynamicsand to investigate a mobility edge in our system.
I. INTRODUCTION
Many-body localization (MBL) is a robust way of er-godicity breaking in interacting many-body systems. Theseminal early works [1–3] show that such a behavior is as-sociated with the presence of strong disorder in the sys-tem. Hundreds of studies in the last 15 years addressedvarious aspects of MBL so we refer the reader to reviewsof the problem [4–8] Despite the numerous works, eventhe issue of the existence of MBL in thermodynamic limitis not settled yet. On one side there is a proof of exis-tence MBL in certain one-dimensional (1D) spin systems[9, 10] under assumptions regarding properties of energylevels. On the other side recent numerical studies putinto question [11] the very existence of MBL as a sta-ble dynamical phase of matter or point out difficulties inits unambiguous verification [12, 13]. Subsequent worksprovided further arguments against the stability of MBLphase [14, 15], or supported the scenario of stable MBLphase [16] indicating that earlier predictions [17] of crit-ical disorder strengths may underestimate the real value[18] (see also [19, 20]). Those works concentrate mainlyon the model of spin − / systems with short range in-teractions and a uniform random disorder present.Instead, the experimental attempts to observe signa-tures of MBL often employ a quasiperiodic potential [21–23] realized by superimposing two incommensurate opti-cal lattices, the primary one holding the atoms (allowing ∗ [email protected] † [email protected] for the tight binding approximate description) and thesecondary one producing the disorder by weakly perturb-ing local minima of the primary lattice. The quasiperi-odic potential leads to MBL similarly to the uniform ran-dom disorder. It is, however, claimed that the univer-sality class of the transition from delocalized, extendedphase to MBL phase is dependent on the disorder type[24, 25]. Similarly, the range of interactions plays a roleas well as the dimension of the problem. We shall restrictourselves to 1D case, still here, the existence of MBL seemto be dependent on the interaction range. In particularfor power-law decaying interactions V ∝ /r α where r is the distance between interacting particles it is claimedthat MBL exists for α > [26, 27]. This condition is nat-urally broken in a system with cavity mediated all-to-allinteractions [28]. The results presented by two of us inthis work, while valid for a finite system may neverthelessbe fragile in the thermodynamic limit as pointed out by[29].The aim of this work is to investigate ergodicity break-ing in the presence of all-to-all cavity mediated interac-tions in the quasi-periodic disorder resulting from themismatch between the cavity mode wavelength and thewavelength of the optical lattice parallel to the cavityaxis. Hence, contrary to [28] we do not assume any ad-ditional random or quasi-periodic on-site disorder. Weshall show that while an MBL regime may be seen fora finite size system, the thermodynamic limit is subtleand the results depend on the details of the model im-plementation. The paper is organized as follows. Themodel and its Hamiltonian are discussed in the next Sec-tion followed by the presentation of spectral propertiesof the system. Then we discuss the time dynamics andconsider different possible quenches in the system. As weshall see some care must be taken in order to excite thesystem in the controlled way. II. HAMILTONIAN
Figure 1. Quasi-periodic infinite-range interactions can berealised by a quantum gas of bosons in an optical lattice (solid,grey line) and strongly coupling to a standing-wave mode of ahigh-finesse optical resonator (dashed blue line). The dipolartransitions couple dispersively with the cavity mode and witha transverse laser beam. As a result the multiply scatteredphotons effectively mediate atom-atom interactions with theperiodicity of the cavity mode wavelength. We assume thatthe ratio between the periodicity of the optical lattice and ofthe cavity mode is incommensurate.
We consider N bosons in a one-dimensional lattice with K sites and open boundary conditions. Their dynamicsis described by the extended Bose-Hubbard model withinfinite-range interactions, whose Hamiltonian reads: ˆ H = ˆ H BH + ˆ H cav . (1)Here, the first term ˆ H BH describes the standard Bose-Hubbard dynamics with nearest-neighbour hopping atrate J > and repulsive onsite interaction with ampli-tude U : ˆ H BH = − J K − X i =1 (ˆ b † i ˆ b i +1 + c.c. ) + U L X i ˆ n i (ˆ n i − , (2)where ˆ b i and ˆ b † i are the annihilation and creation opera-tors for a boson on site i , and ˆ n i = ˆ b † i ˆ b i is the number ofparticle at site i . The infinite-range interactions is givenby the Hamiltonian [30]: ˆ H cav = − U K K X i =1 Z β ( i, φ )ˆ n i ! . (3)where U scales the interaction amplitude, while Z β is afunction of the site i and of the real parameters β and φ , that is bounded to unity, |Z β | ≤ . This functionhas periodicity β in units of the lattice periodicity. In the Appendix we provide its specific dependence on theparameters of the setup of Fig. 1. φ is an arbitraryphase factor corresponding to different realisations of theeffective “disorder”. In this work we choose β to be anirrational number, and in particular β = ( √ − / . Asa result Hamiltonian (1) is quasiperiodic.In the rest of this manuscript we assume that U > .For this choice, the interaction term favours quasiperi-odic density distributions which maximize the expecta-tion value of Eq. (3) and the quantum ground state ex-hibits the features of a Bose glass [30, 31]. In these worksit was shown that Hamiltonian (1) can be realised ina cavity quantum electrodynamics setup by overlappingthe cavity mode to an optical lattice, tightly confining theatoms along the cavity axis. The setup is illustrated inFig. 1. In this configuration parameter β is proportionalto the ratio between the cavity wave number k and theoptical lattice wave number k l , namely, β = k/ (2 k l ) . It isknown [32] that the value of β affects the distribution ofon-site energies between neighboring sites so the resultsquantitatively depend on the β choice.We finally remark that the single-particle ground stateand spectrum of these dynamics have been studied inRefs. [33, 34]. Here, it was shown that the ground statecan be localized for sufficiently large U . For these valuesthe spectrum exhibits mobility edges.In the rest of this manuscript we analyse the dynamicsfor J = U and report the interaction strength U in unitsof J . III. SPECTRAL PROPERTIES
The localisation properties of our system are studiedusing exact diagonalisation (ED) of the Hamiltonian of asystem with a unit filling for different number of sites K ranging from 7 to 9. These modest system sizes are de-termined by necessity of including the possibility to haveseveral bosons at single sites and correspond to dimen-sions of Hilbert space between 1716 and 24310. All quan-tities obtained are also averaged over 500-30000 random φ phases - corresponding to independent realizations ofquasi-periodic disorder.We start our analysis with the density of states (DOS)and its dependence on the system size, K and cav-ity mediated interaction strength U , shown in Fig. 2.We rescale eigenenergies E to a unit interval as ǫ =( E − E min ) / ( E max − E min ) with E min , E max denotinga minimal and maximal energy from each disorder re-alization, respectively. For moderate value of all-to-allinteractions, U = 5 , the maximal DOS occurs near themiddle of the spectrum. For the higher value U = 25 the maximum of DOS is shifted to higher energies, withlow and middle energy states laying in the tail of DOS.When the system size is increased the shift of DOS toslightly lower energies is observed. The strong asym-metry of DOS reflects itself in localization properties ofeigenstates in different parts of the spectrum as we show . . . . . . ǫ . . . . . . . D O S K = 7 , U = 5 K = 7 , U = 25 K = 9 , U = 5 K = 9 , U = 25 Figure 2. The density of states (DOS) for K = 7 , and U = 5 , as indicated in the figure. below. U . . ¯ r a) K = 7 K = 8 K = 9 0 10 20 U . . ¯ r b) Figure 3. The dependence of the mean gap ratio ¯ r on theall-to-all interaction strength U for different energy ranges:a) ǫ ∈ [0 . , . , b) ǫ ∈ [0 . , . . To study the localization properties of eigenstates weutilize the mean gap ratio [3]. It is defined as ¯ r = (cid:28) min (cid:26) ǫ i +1 − ǫ i ǫ i − ǫ i − , ǫ i − ǫ i − ǫ i +1 − ǫ i (cid:27)(cid:29) , (4)where the averaging is performed in a band of eigenener-gies around a given rescaled energy ǫ and over disorderrealizations. For ergodic system, with level statistics welldescribed by Gaussian Orthogonal Ensemble of randommatrices, ¯ r ≈ . while ¯ r ≈ . indicates localisationand Poissonian level statistics [35]. The dependence ¯ r on the interaction strength U for different parts of spec-trum and system sizes is shown in Fig.3. The ¯ r plotsshow a couple of interesting properties of our system.Firstly, for a fixed system size and energy range, the sys-tem undergoes a crossover between ergodic and localizedregimes as we increase U , meaning that stronger inter-actions actually prohibit thermalization. Thus, the inter-action induced randomness to some extent plays a role ofa (quasi)random on-site potential in the interacting sys-tem. A similar phenomenon occurs for bosons interactingby random contact interactions [36–38]. Here, we keepthe uniform on-site interactions fixed at low U = 1 valueand the randomness comes from infinite-ranged cavitymediated interactions only. A second observation is that the higher the energy is, the stronger U is necessary toinduce the crossover to localized phase. Therefore, for afixed U , we observe a so-called mobility edge, namelythe lower energetically part of the spectrum is localized,while the higher is thermal (delocalized). However, Fig. 3b) shows that the value of U at which the system entersinto the localized regime increases with system size K ,implying that the mobility edge shifts to lower energieswith increasing K . This suggests that the mobility edgemay not survive in the thermodynamic limit and in thislimit all eigenstates are extended. V . . ¯ r a) K = 7 K = 8 K = 9 0 1 2 3 V . . ¯ r b) Figure 4. The dependence of ¯ r on V = U /K for variousenergy ranges. a) ǫ ∈ [0 . , . , b) ǫ ∈ [0 . , . . Observe, however, that we can choose a different scal-ing the coupling strength in ˆ H (0) cav term, (3). The factor infront of the sum, U /K , has been chosen, in accordancewith the customary approach [30, 39] that ensures theappropriate extensive scaling of the Hamiltonian with in-creasing K . In fact, it may be justified by assuming thatthe increase of the system size K , for a fixed wavenumber k , implies the increase in the cavity length and thus theappropriate scaling of the cavity mode volume. However,another approach is possible - to keep the cavity lengthfixed and change the lattice spacing. Then, instead ofthe coefficient U /K in front of the all-to-all interactionsterm, we would have a constant V that remains fixedwhen K is increased. Such a scaling implies a nonstan-dard thermodynamic limit as the energy becomes a su-per extensive quantity. A similar procedure was adoptedin [40] to study MBL in presence of power-law interac-tions. We show the mean gap ratio ¯ r as a function of V in Fig. 4. For a sufficiently large V a crossover to lo-calized regime is observed. However, contrary to fixed U case, the resulting r ( W ) curves deviate from the er-godic value r ≈ . at a system size independent valueof V and the crossover becomes steeper with increasingsystem size. This suggests that the observed crossoverbecomes a sharp transition to an MBL phase at a fixedcritical value of V in the considered non-standard ther-modynamic limit. We note that a similar behavior atthe ergodic-MBL crossover was observed for disordered U (1) lattice gauge theory [41], in which the eliminationof gauge field leads to long-range interactions.Let us note that identifying the transition for V = const. limit implies the delocalized limit for U = const. for arbitrary U . On the other hand all experiments areperformed for a finite K value that may reach severaltens or hundreds but necessarily remains finite so ourobservations for the standard U = const scenario areexperimentally relevant. . . . ǫ . . . S / S a) K = 7 K = 8 K = 9 0 . . . ǫ . . . S / S b)0 . . . ǫ . . . S / S c) 0 . . . ǫ . . . S / S d) Figure 5. The ǫ dependence of von Neumann biparite entan-glement entropy S/S of eigenstates normalized by a randomGaussian state value S . Subsequent panels correspond to dif-ferent values of V = U /K : a) V = 0 . , b) V = 1 , c) V = 2 ,d) V = 3 . The crucial difference between extended and localizedstates can be seen in their entanglement properties. Forextended states quantum entanglement between subsys-tems which we call A and B resulting from splitting ofthe entire system into two, follows a volume law, whilefor localized states area law is expected. To quantify theentanglement, we calculate the von Neumann bipariteentanglement entropy defined as S = − Tr B [ˆ ρ B ln ˆ ρ B ] , (5)where ˆ ρ B = Tr A [ˆ ρ AB ] is a (usually mixed) state of thesubsystem B obtained by tracing the state ˆ ρ AB = | ψ i h ψ | over degrees of freedom of complementary subsystem A . If subsystems A and B are not entangled, ˆ ρ B is apure state and S = 0 . Fig. 5 presents the entangle-ment entropy of eigenstates ( ˆ ρ AB = | ψ n i h ψ n | , where ˆ H | ψ i = E n ψ n ) averaged over eigenstates with energiescorresponding to the rescaled energy ǫ for different val-ues of V and K , normalized by the entropy S of randomGaussian states [42]. In this case, A and B subsystemsare two equal K/ parts.Again we observe two main features. Firstly, for low ly-ing energy states the normalized entropy is close to 0 andindependent of system size K , indicating localized states.However for higher energies ǫ , the entanglement entropy S approaches S as expected for thermal states. Inter-estingly, for very high energies we observe the decreasein S , but this appears at the edge of the spectrum wherethe density of states is again very low. Secondly, increas-ing V shifts the the thermal region to higher energies,which is consistent with level statistics. What is impor-tant, for constant V and changing K we again observe the crossing of curves corresponding to different systemsizes. Interestingly, the crossing is visible on both sidesof S maximum, suggesting that the localized eigenstates,fulfilling the area law, found at the highest energies ǫ canalso survive in the non-standard thermodynamic limit. − t − . . . . . . I U = 25 U = 15 U = 10 U = 5 U = 1 Figure 6. Population imbalance I ( t ) during time evolution ofa density-wave like state for K = 8 and a several of U values:25, 15, 10, 5, 1 respectively, starting from the top. IV. TIME DYNAMICS
In the previous Section we analysed spectral proper-ties discussing various aspects of the crossover betweenergodic and MBL regimes as well as the presence of themobility edge. Here, we address more experimentally rel-evant issue - whether the localization properties may beobserved in the time dynamics starting from a judiciouslyprepared initial state. Following the usual experimentalstrategy [21, 22] we might consider the population imbal-ance, defined as I ( t ) = * K K X i =1 ( − i ˆ n i ( t ) + . (6)An initial state with an imprinted density pattern, forwhich the imbalance I ( t ) is non-zero at t = 0 can be usedto demonstrate ergodicity breaking. For unit density as-sumed by us such a state can be chosen as | . . . i .The system evolves governed by ˆ H and the I ( t ) is mon-itored to reveal the system properties. The initial unitimbalance decays to 0 for systems obeying the eigenstatethermalization hypothesis [43, 44]. In such a case theinformation about the initial configuration is lost dur-ing the course of time evolution. In contrast, a non-zero long-time value of I suggests that the informationis partially preserved and is indicative of the localizedphase. The evolution of the population imbalance for aseveral values of U is shown in Fig. 6. The results areobtained using Chebyshev propagation [45, 46] relying onsparse matrix operations, which allows us to reach systemsizes up to K = 12 with reasonable computational timeand resources. Surprisingly, although the level statistics(Figs. 3,4) indicate localization in some ranges of U , forall parameters studied here, the imbalance for the ini-tial density wave state rapidly decays to zero suggestingthermalization. It is easily understood after inspectingthe energy of the density wave state | ... i (energy isconserved during propagation) which, for all reasonablechoices of the parameters of the problem, corresponds toregions identified as ergodic by the mean gap ratio. −
20 0 20( U ) . . . ǫ a) −
20 0 20( U ) . . . ǫ b) .
625 0 .
650 0 . β . . . ǫ c) 0 .
625 0 .
650 0 . β . . . ǫ d) . . . . φ [rad]0 . . . ǫ e) 0 . . . . φ [rad]0 . . . ǫ f) Figure 7. The dependence of average energy ǫ of the quenchedstate w.r.t. the final Hamiltonian on several quenched param-eters for K = 8 . On both panels in each row the same energyvalues are presented, but with different types of errors rep-resented by vertical lines. In the left column they representthe quantum error rD ˆ H final E − D ˆ H final E , while in the rightcolumn – the standard deviation of the random phase distri-bution of mean energies. One can imagine a preparation of a different initialproduct state, e.g. a uniformly occupied state | . . . i and use the transport distance [47] (see also below) tolook for signatures of localization. Our attempts in thatdirection failed (except for highly unusual states whichdid not allow us to study the breaking of ergodicity inour system in a systematic manner). The reason for thatis simple. To reach the MBL regime in our system weneed a sufficiently strong disorder which appears solelyin the cavity mediated interaction term. For large U this term dominates the energy of product states. Thus,the occupation numbers of an initial product state mustbe adapted in such a way that the all-to-all interactionterm is possibly small in order to reach localized regime.To inspect the low energy regime in a more systematic manner we employ an alternative approach introducedby [48] and termed the quantum quench spectroscopy(QQS). In the first step of QQS the system is preparedin the ground state of some initial Hamiltonian ˆ H ini . Anabrupt change of one of its parameters to the desired ˆ H follows. In such a procedure the initial state is notchanged but now becomes a non-stationary initial su-perposition of eigenstates of ˆ H . One may hope that anappropriate choice of the initial Hamiltonian ˆ H ini mayallow to probe different parts of the spectrum of the finalHamiltonian ˆ H . It was shown in [48] that this proce-dure, when ramping the disorder strength, may allow oneto probe the mobility edge in the system of interactingspinless fermions in a quasiperiodic potential.It is important, however to realize two possible limita-tions. As the initial state is not an eigenstate of ˆ H itsenergy suffers from quantum uncertainty ∆ E = rD ˆ H E − D ˆ H E , (7)the fate shared, in fact, also by the initial product states.Moreover, as a target mean energy depends on a random φ phase in (3), it is also a random variable, whose disper-sion can be characterized by a sample standard deviation.This second uncertainty is known to decrease with an in-creasing system size [48], so it is not a major concern inthe thermodynamic limit. On the other hand it may bean obstacle for relatively small system sizes treated byED in this work.Comparison of mean energies obtained and both uncer-tainties is shown in Fig. 7 for different possible quenchescenarios. We fix β = ( √ − / ≈ . and U = 25 asthe final values. In the first approach we quench the dis-order amplitude from some initial ( U ) value to U - thepath similar to [48]. The results presented in the upperrow of Fig. 7 show low energy selectivity. For small dif-ference between the initial and final disorder amplitudelow excitations are reached only. The rapid jump in rela-tive energy ǫ appears for low initial interaction strength.When ( U ) has a different sign than U a relatively large ǫ values may be reached but with quite significant quan-tum uncertainty. Similarly, the region of small ( U ) where a rapid change of final energies is observed suffersfrom that uncertainty.Situation is entirely different if, instead of the interac-tion strength, the ratio of cavity and laser wavenumbersis changed, i.e. β as shown in the middle row of Fig. 7.The quantum uncertainties are small and a broad rangeof final energies may be achieved, however φ averaginguncertainty is very large for the system sizes under con-sideration. Moreover, the maximal energy achieved isaround ǫ ≈ . , which is below the mobility edge esti-mated for K = 8 from gap ratio as ǫ ≈ . − . notenabling us to probe the crossover between ergodic andMBL regimes.Finally, the bottom row in Fig. 7 shows yet anotherpossible quench. For each (random) φ we choose a dif-ferent φ (0) and define the quench as a rapid change ofthe phase from φ (0) to φ . We denote the magnitude ofthe phase change by ∆ φ . As observed the quantum andstatistical uncertainties are relatively small. Moreover,the final energy approaches 1 for ∆ φ = π/ . Thus, weperform further the QQS simulations using this type ofquench. − − t − − − ∆ x a) − − t . . . . F b) 0 . . . . . . . . ǫ . . . . V x c) K = 7 K = 8 K = 9 K = 10 K = 12 . . . ǫ . . V F d) . . . ǫ h ∆ x i e) . . . ǫ h F i f) Figure 8. a) Time evolution of transport distance for K = 8 and a several ∆ φ values. b) A similar evolution of fluc-tuations F . c) Maximal transport pseudo-velocity depen-dence on initial state energy for a several K . d) Fluctuationspseudo-velocity energy dependence. e) Average value of ∆ x for t ∈ [100 , energy dependence. f) Average value of ∆ x for the same time interval energy dependence. After a change of the cavity wavenumber by ∆ φ weconsider time evolution of the system. Instead of theImbalance we calculate, following [47] the transport dis-tance, defined as ∆ x = 2 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K − X d =1 d × D G (2) c ( i, i + d ) E i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (8)where G (2) c ( i, i + d ) = h ˆ n i ˆ n i + d i−h ˆ n i i h ˆ n i + d i and averagingover sites h . . . i i is done before averaging over realizations ( . . . ) . Note that the measurement of ∆ x requires accessto two-site correlations available in state of the art ex-periments [47]. The second quantity we calculate is anedge fluctuation defined as F = (cid:16)(cid:10) ˆ n (cid:11) − h ˆ n i (cid:17) + (cid:16)(cid:10) ˆ n K (cid:11) − h ˆ n K i (cid:17) . (9)It requires a site dependent resolution at the edges ofthe system but does not involve two-point correlations. Both quantities were analyzed previously with the aim ofidentifying the mobility edge in a Bose-Hubbard modelin a tilted lattice [49]. Time evolution of both quan-tities for K = 8 , U = 25 and different initial stateenergies are depicted in the first row of Fig. 8. Ini-tial states were prepared in quantum quench of φ , with ∆ φ = 0 . , . , . , . and . , where two first val-ues yield ǫ ≈ . and . in the localized regime, thirdgives ǫ ≈ . in the vicinity of mobility edge, fourthgives ǫ ≈ . which is well in extended regime and fifth– ǫ ≈ . , which is near the high end of spectrum, wherethe states localize again. Observe that ∆ x ( t ) evolutionstarts with an abrupt, ballistic-like transient behaviour,but after a time of the order of inverse tunneling energy /J = 1 one can observe a slower, subdiffusive growth.When the Heisenberg time scale (which is proportionalto the inverse of the mean level spacing) is reached, thetransport distance value saturates. For energies corre-sponding to the localized part of the spectrum, ∆ x satu-rates at low values below a single site unit distance whilefor those in the thermal domain ∆ x is of the order of sys-tem size K = 8 . Here, in fact, the saturation is reachedwell before the Heisenberg time when the system size isreached. In the vicinity of the crossover between ergodicand MBL regimes the saturation of ∆ x occurs only laterand at an intermediate value (of the order of unity). No-tice, that ∆ x for the highest energy considered . satu-rates at lower value then for . confirming the presenceof the localized regime at the high end of the spectrum.We observe a similar situation for the edge fluctuation, F . The initial transient behavior is followed by a typical,logarithmic growth followed by a saturation of F ( t ) .It was proposed in [49] to detect the mobility edge us-ing the notion of pseudovelocity. The authors computedthe ratio V O = O ( t max ) − O ( t min ) t max − t min (10)where O ( t ) = ∆ x or F for fixed times t max and t min for evolution of states with different energies and plottedthem against ǫ . They noticed a peak of this value nearmobility edge, providing an experimental way for its ob-servation. However, in our case, the mobility edge isrelatively high in the rescaled energy ǫ hindering the ob-servation of a dip of pseudovelocity after the initial rise.Also, for large K and maximal simulation time t = 10 some observables are not yet saturated. In effect thepseudovelocities obtained had large fluctuations. In or-der to mitigate that issue, we modified the approach –instead of fixing t max and t min , we scanned ∆ x ( t ) and F ( t ) dependencies with ∆ t = 30 window starting from t = 10 , computed a linear fit slope in this window andchose its maximal value as a pseudovelocity. This allowedus to make V O ( ǫ ) dependencies much smoother (Fig. 8c,d). For both V x and V F , on the localized side of spectrumthe values are near 0 and start to visibly rise around thetransition energy ǫ ≈ . . They also correctly recognizethe shift of transition energy to lower values for largersystem sizes. Experimental accessibility of ∆ x is lowerthan F , because it requires measuring all the correla-tions, while the latter one concerns the sites on the edgesonly. On the other hand V F is more affected by fluc-tuations than V x . Instead of pseudovelocities one may,however, consider simply the average values of observ-ables in a fixed time interval as shown in the bottom rowof Fig. 8. Their behavior is much smoother, providing asimilar observation – a rapid increase of both the meantransport distance and mean fluctuation at the crossoverfrom localized to extended regime. − t . . . S / S a) 0 . . . . .
50 10 − t . . . S / S b)10 − t . . . S / S c) . . . ǫ . . . S m a x / S d) K = 8 K = 10 K = 12 Figure 9. a)-c) Time evolution of von Neumann biparite en-tanglment entropy for ∆ φ -quenched states for, respectively, K = 8 , , . d) Maximal value (for t = 10 ) of S ( t ) /S as afunction of energy ǫ for those 3 system sizes. We may also consider time evolution of the entangle-ment entropy using the same quenched initial states. Forinitially separable state the hallmark of MBL is the log-arithmic entanglement entropy growth [50–52]. This hasbeen typically tested for separable, product-like initialstates. It is interesting, therefore, to observe that en-tropy growth also for our initial entangled states - seeFig. 9. For localized initial states S grows very slowly.In the vicinity of the crossover we observe a delayed sat-uration, while in the extended regime S saturates at thevalue of the order of the average entanglement entropyof a random Gaussian state S . However, S is neverreached for system sizes and times under consideration.One may notice, that, once again, the saturation valuefor highest energy ǫ = 0 . is lower that for ǫ = 0 . ,which is consistent with previous results. We also plot,in Fig. 9(d), the dependence of maximal S max againstinitial state energy. In this case, the visible rise of S max starts around ǫ ≈ . , below the transition energy. Nev-ertheless, it provides additional evidence supporting themobility edge, at least for system sizes studied in thispaper. V. SUMMARY
We considered signatures of many body localization inthe cavity QED model combined with the optical lat-tice whose wavelength is incommensurate with the cav-ity mode wavelength. In such a case the effective long-range interactions between atoms (bosons) in the systembecome dependent on the incommensurability ratio re-sulting in the effective quasi-random potential. Impor-tantly, we do not assume any other type of disorder tobe present.We analyse both spectral properties of the system(mean gap ratio as well as entanglement entropy of eigen-states) as well as the time dynamics from specially pre-pared intitial states. For finite (small, amenable to exactdiagonalization) systems we observe clear signatures ofMBL with mean gap ratio reaching the Poisson level char-acteristic for a localized regime. However, the interactionstrength leading to localization seems to increase with thesystem size. In effect, in the standard thermodynamiclimit, we predict that the system will remain delocalizedfor almost all the states. One may define, however, anonstandard thermodynamic limit, in which the interac-tion strength is not scaled with system size. While thismakes the energy of the system a super-extensive quan-tity, it also gives rise to a clear mobility edge betweenlocalized and delocalized regimes.These findings, based on spectral statistics, were con-firmed and verified by time dynamics in which we ad-dressed quantities that are within the reach of currentexperiments. Interestingly, meaningful studies requiredquantum quench spectroscopy with initial states ob-tained by a parameter quench. We have shown that agreat care should be taken in the choice of the quenchedparameter, in particular for necessarily small system sizesconsidered. We have also shown that in the localizedregime the entanglement entropy growth for such initialweakly entangled states follows a logarithmic law, a hallmark of MBL observed typically for initial product states.Our results indicate that many-body cavity quantumelectrodynamics setups are well suited to studies of quan-tum ergodicity and its breaking. The setup considered inthis work is considerably simpler than the system investi-gated in [28]. Here, the cavity mediated interactions acteffectively as a source of randomness leading to a non-trivial interplay between ergodic and MBL eigenstates atdifferent energies.
ACKNOWLEDGEMENTS
The support by National Science Centre (Poland) un-der project OPUS 2019/35/B/ST2/00034 (J.Z.) is ac-knowledged. P.K. would like to acknowledge the supportof Ministry of Science and Higher Education (Poland)grant no. 0108/DIA/2020/49. P.S. acknowledges thesupport of Foundation for Polish Science (FNP) throughscholarship START. G. M. acknowledges the support bythe German Research Foundation (Project-ID 429529648– TRR 306 QuCoLiMa "Quantum Cooperativity of Lightand Matter) and by theGerman Ministry of Educationand Research (BMBF), QuantERA project NAQUAS.Project NAQUAS has received funding from the Quan-tERA ERA-NET Co-fund in Quantum Technologiesimplemented within theEuropean Union?s Horizon2020program. Numerical simulations were carried out withthe support of the Interdisciplinary Center for Mathe-matical and Computational Modeling (ICM) at the Uni-versity of Warsaw under grant no. GB76-1. For linear al-gebra computations, Armadillo C++ language bindings[53, 54] with Intel® Math Kernel Library as BLAS andLAPACK backend were used.
Appendix A: Details on the Bose-Hubbard model
We consider the setup of Fig. 1 in the regime wherethe cavity mode dynamics is characterized by a fastertime scale than the atomic motion. In this limit one canapply a coarse graining to the time evolution [28, 30].Eliminating the cavity degrees of freedom leads to an ef-fective infinite-range interaction term [55]. By assumingthat the atoms are tightly bound in the lowest band of an optical lattice with wave number k l , the interactionterm can be cast in the form of Eq. (3), with [28, 30, 39] Z β ( i, φ ) = Z d x cos( kx + φ ) w i ( x ) w j ( x ) , (A1)where w i ( x ) is lowest band Wannier function centred at i th site. The other quantities are the wave vector of thecavity mode k , the interaction strength U , and an ar-bitrary phase α . For a sufficiently deep optical latticepotential, we approximate w i ( x ) ≈ δ ( x − π/k l ( i + 1 / and obtain: Z β ( i, φ ) ≈ cos(2 πβi + φ ) , (A2)with β = k/ (2 k l ) .This approach assumes that the cavity dynamics adi-abatically follows the atomic motion. In particular thecoarse-grained cavity field operator is now proportionalto the atomic observables ˆ a ∝ K X i =1 Z β ( i, φ )ˆ n i . (A3)This quantity can be revealed by homodyne detection ofthe field at the cavity output [56, 57]. [1] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov,Phys. Rev. Lett. , 206603 (2005).[2] D. Basko, I. Aleiner, and B. Altschuler, Ann. Phys. (NY) , 1126 (2006).[3] V. Oganesyan and D. A. Huse,Phys. Rev. B , 155111 (2007).[4] R. Nandkishore and D. A. Huse,Annual Review of Condensed Matter Physics , 15 (2015),https://doi.org/10.1146/annurev-conmatphys-031214-014726.[5] D. J. Luitz and Y. B. Lev,Annalen der Physik , 1600350 (2017).[6] F. Alet and N. Laflorencie,Comptes Rendus Physique , 498 (2018), quantumsimulation / Simulation quantique.[7] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn,Rev. Mod. Phys. , 021001 (2019).[8] S. Gopalakrishnan and S. Parameswaran,Physics Reports , 1 (2020), dynamics and transportat the threshold of many-body localization.[9] J. Z. Imbrie, Phys. Rev. Lett. , 027201 (2016).[10] J. Z. Imbrie, Journal of Statistical Physics , 998 (2016).[11] J. Šuntajs, J. Bonča, T. Prosen, and L. Vid-mar, arXiv e-prints , arXiv:1905.06345 (2019),arXiv:1905.06345 [cond-mat.str-el].[12] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor,and M. Žnidarič, Europhysics Letters , 67003 (2020).[13] P. Sierant, D. Delande, and J. Zakrzewski,Phys. Rev. Lett. , 186601 (2020).[14] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer,and J. Sirker, Phys. Rev. Lett. , 243601 (2020).[15] D. Sels and A. Polkovnikov, “Dynamical obstruction to localization in a disordered spin chain,” (2020),arXiv:2009.04501 [quant-ph].[16] D. J. Luitz and Y. B. Lev,Phys. Rev. B , 100202 (2020).[17] D. J. Luitz, N. Laflorencie, and F. Alet,Phys. Rev. B , 081103 (2015).[18] P. Sierant, M. Lewenstein, and J. Zakrzewski,“Polynomially filtered exact diagonalization ap-proach to many-body localization,” (2020),arXiv:2005.09534 [cond-mat.dis-nn].[19] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D.Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi,Phys. Rev. B , 174202 (2018).[20] T. Chanda, P. Sierant, and J. Zakrzewski,Phys. Rev. B , 035148 (2020).[21] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch,and C. Gross, Science , 1547 (2016).[22] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet,E. Altman, U. Schneider, and I. Bloch,Phys. Rev. Lett. , 260401 (2017).[23] H. P. Lüschen, S. Scherg, T. Kohlert, M. Schreiber,P. Bordia, X. Li, S. Das Sarma, and I. Bloch,Phys. Rev. Lett. , 160404 (2018).[24] V. Khemani, D. N. Sheng, and D. A. Huse,Phys. Rev. Lett. , 075702 (2017).[25] A. Maksymov, P. Sierant, and J. Zakrzewski,“Many-body localization in one dimensional op-tical lattice with speckle disorder,” (2020),arXiv:2008.00219 [cond-mat.dis-nn].[26] A. L. Burin, Phys. Rev. B , 094202 (2015). [27] A. L. Burin, Phys. Rev. B , 104428 (2015).[28] P. Sierant, K. Biedroń, G. Morigi, and J. Zakrzewski,SciPost Phys. , 8 (2019).[29] A. O. Maksymov and A. L. Burin,Phys. Rev. B , 024201 (2020).[30] H. Habibian, A. Winter, S. Paganelli, H. Rieger, andG. Morigi, Phys. Rev. Lett. , 075304 (2013).[31] H. Habibian, A. Winter, S. Paganelli, H. Rieger, andG. Morigi, Phys. Rev. A , 043618 (2013).[32] E. V. H. Doggen and A. D. Mirlin,Phys. Rev. B , 104203 (2019).[33] K. Rojan, R. Kraus, T. Fogarty, H. Habibian, A. Min-guzzi, and G. Morigi, Phys. Rev. A , 013839 (2016).[34] J. Major, G. Morigi, and J. Zakrzewski,Phys. Rev. A , 053633 (2018).[35] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Phys. Rev. Lett. , 084101 (2013).[36] P. Sierant, D. Delande, and J. Zakrzewski,Phys. Rev. A , 021601 (2017).[37] P. Sierant, D. Delande, and J. Zakrzewski,Acta Phys. Polon. A , 1707 (2017).[38] P. Sierant and J. Zakrzewski,New Journal of Physics , 043032 (2018).[39] N. Dogra, F. Brennecke, S. D. Huber, and T. Donner,Phys. Rev. A , 023632 (2016).[40] S. Gopalakrishnan and D. A. Huse,Phys. Rev. B , 134305 (2019).[41] G. Giudici, F. M. Surace, J. E. Ebot, A. Scardicchio, andM. Dalmonte, Phys. Rev. Research , 032034 (2020).[42] L. Vidmar and M. Rigol,Phys. Rev. Lett. , 220603 (2017).[43] J. M. Deutsch, Phys. Rev. A , 2046 (1991). [44] M. Srednicki, Phys. Rev. E , 888 (1994).[45] H. Tal-Ezer and R. Kosloff,J. Chem. Phys. , 3967 (1984),https://doi.org/10.1063/1.448136.[46] H. Fehske and R. Schneider, Computational many-particle physics (Springer, Ger-many, 2008).[47] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai,J. Léonard, and M. Greiner, Nature , 385 (2019).[48] P. Naldesi, E. Ercolessi, and T. Roscilde,SciPost Phys. , 010 (2016).[49] R. Yao and J. Zakrzewski, “Many-body localization ofbosons in optical lattice: Dynamics in disorder-free po-tentials,” (2020), arXiv:arXiv:2007.04745.[50] M. Žnidarič, T. Prosen, and P. Prelovšek,Phys. Rev. B , 064426 (2008).[51] J. H. Bardarson, F. Pollmann, and J. E. Moore,Phys. Rev. Lett. , 017202 (2012).[52] M. Serbyn, Z. Papić, and D. A. Abanin,Phys. Rev. Lett. , 260601 (2013).[53] C. Sanderson and R. Curtin, Journal of Open SourceSoftware , 26 (2016).[54] C. Sanderson and R. Curtin, Lecture Notes in ComputerScience (LNCS) , 422 (2018).[55] S. Schütz, S. B. Jäger, and G. Morigi,Phys. Rev. A , 063808 (2015).[56] A. T. Black, H. W. Chan, and V. Vuletić,Phys. Rev. Lett. , 203001 (2003).[57] K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger,Phys. Rev. Lett.107