Many-body localization in presence of cavity mediated long-range interactions

We show that a one-dimensional Hubbard model with all-to-all coupling may exhibit many-body localization in the presence of local disorder. We numerically identify the parameter space where many-body localization occurs using exact diagonalization and finite-size scaling. The time evolution from a random initial state exhibits features consistent with the localization picture. The dynamics can be observed with quantum gases in optical cavities, localization can be revealed through the time-dependent dynamics of the light emitted by the resonator.


Introduction
Many-body localization (MBL) is the most robust manifestation of ergodicity breaking in interacting many-body systems. The interactions are generically considered as leading to ergodic dynamics as far as local observables are concerned. The standard formulation of this belief is the eigenvector thermalization hypothesis (ETH) [1,2]. Numerous studies over the last decade have shown that the many-body interacting systems do not thermalize in presence of a quenched disorder. In particular, for a sufficiently strong disorder MBL may occur, leading to a long-time memory of the initial state (for recent reviews see [3] as well as a topical issue of Annalen der Physik [4]).
While most studies so far considered short-range interactions, for instance in the XXZ spin Hamiltonian (which became a paradigmatic model of MBL [5]) or bosons and fermions in optical lattices in tight binding models [6][7][8][9][10], it is by no means clear whether MBL persists for genuinely long-range interactions, such as for Coulombic or dipolar potentials. In recent works [11] it was argued that MBL may appear in disordered long-range interacting systems, others suggest [12][13][14][15] lack of MBL for, e.g., dipole-dipole interactions in three dimensions (3D) [16,17]. Furthermore, MBL in presence of power-law decaying tunneling elements and interactions has been studied in [18][19][20] where algebraic decay of correlation functions as well as of algebraic growth of entanglement entropy was found.
In the present work we numerically analyse whether MBL occurs in a disordered Hubbard model with all-to-all interactions. This model is expected to describe the dynamics of atoms in an external lattice and interacting dispersively with a mode of a standing-wave optical cavity. [21] Long-range interactions appear naturally in this system -the mode of the cavity mediates a two-body interaction whose range is as large as the system size. Our investigations are based on exact diagonalization supplemented by numerical techniques for sparse Hamiltonian matrices for a system of fermions (or bosons). Our results show features which can be attributed to the occurence of MBL in the system. The paper is structured as follows. In Sec. 2 we describe an experimentally realizable system of atoms inside a resonant cavity. In Sec. 3 we determine the phase diagram by a finite-size scaling analysis assuming that the atoms are spinless fermions. In Sec. 4.1 we turn to dynamics of the system and provide the evidence of ergodicity breaking resulting from the interplay of disorder and interactions in the system. In Sec. 4.2 we show that the ergodicity breaking can be understood within an appropriately modified picture of local integrals of motion (LIOMs). We then argue that in cavity QED setups ergodicity breaking can be revealed via the light emitted at the cavity output. In Sec. 6 we provide arguments that the reported properties on the system are similar both for purely random disorder as well as for quasiperiodic potential. Finally, we discuss MBL for bosons loaded into the cavity system. Details of the time-propagation algorithm are to be found in the Appendix.

Description of the system
We consider an ensemble of atoms trapped in a quasi 1D geometry and tightly bound by an optical lattice. The atoms dispersively interact with an optical cavity in the regime, in which the cavity mode can be adiabatically eliminated from the dynamics. We assume that the cavity field is a perturbation to the optical lattice, so that the dynamics can be restricted to the lowest band of the external lattice and the cavity-mediated long-range interactions effectively describe all-to-all interactions between the lattice sites. The dynamics is governed by the effective Hamiltonian with H A being the standard Hubbard-like Hamiltonian for the dynamics of atoms in an optical l/2 k z x W z (x) Figure 1: (color online) The Hubbard model we consider describes the dynamics of atoms in an optical lattice and interacting with the standing-wave mode of a high-finesse optical cavity. In the limit in which the atom-photon interactions are dispersive and the cavity field can be adiabatically eliminated from the atomic dynamics, the resulting Hubbard model is characterized by all-to-all interactions mediated by the cavity photons. The trasverse arrow symbolizes an external laser pumping photons into the cavity via coherent atom scattering. The additional disorder shifts the energy of the lattice potential. The quantum state of the system can be inferred by measuring the light at the cavity output, by time-of-flight measurements, or by Bragg spectroscopy using a weak probe.
lattice and in presence of disorder and H C the energy of the interaction with the cavity field. In detail, the optical lattice is composed of K sites, we assume periodic boundary conditions and the atomic Hamiltonian reads where b j and b † j are the onsite annihilation and creation operators of a fermion or a boson at site j = 1, . . . , K, (with K + 1 identified with the first site), J is the tunneling coefficient scaling the nearest-neighbour hopping, n j = b † j b j is the occupation operator at site j, E j is the onsite energy at site j, and H F,B int is the interaction term, which takes different forms depending on the quantum statistics of the atomic gas. For spinless fermions and U > 0, whereas the first non-trivial interaction term in tight-binding expansion for bosons reads In turn, the cavity-mediated long-range interactions take the form [22][23][24]: with U 1 > 0. This interaction is here derived under the assumption that the wavelength of the cavity field equals the one of the electric field generating the optical lattice. This interactions favour density-wave (DW) ordering. We determine the existence of the MBL phase by exact diagonalization of the Hubbard Hamiltonian, focusing majorly on spinless fermion. This situation is in fact more accessible to numerical analysis since the local Hilbert space has only dimension 2. We then only briefly show that analogous properties of the fermionic case are also found for bosons. We finally note that the disorder in our model is in the onsite energy. Here we assume two cases. Throughout most of the work we make the theoretically elegant assumption that E j are uncorrelated random variables uniformly distributed in [−W, W ] interval, where 2W denotes the interval width. In Sec. 6 we then analyse the situation where E j is due to a quasi-periodic optical potential.
In the rest of this manuscript we report energies in units of J and time in units of 1/J.

Phase diagram
Energy level statistics encode an answer to the question of whether a disordered system is localized or ergodic and satisfies ETH. Level statistics of ergodic systems with (generalized) time reversal symmetry have properties akin to the Gaussian Orthogonal Ensembe (GOE) [25]. As the disorder strength increases and the system becomes localized, the level statistics becomes Poissonian [26,27] (an accurate model for level statistics across the localization transition was recently proposed in Refs. [28,29]). The level statistics can be characterized using the gap ratio. This is defined as with δ n = E n − E n−1 being the spacing between two consecutive eigenvalues [30]. Averaging over different energy levels within a certain interval as well as over disorder realizations results in the mean gap ratio, r, that may be used to characterize the spectra. The mean gap ratio changes from r ≈ 0.53 in the ergodic regime [30,31] to r ≈ 0.39 for a localized system and is thus a straightforward probe of the MBL transition especially as it does not require level unfolding, a tricky procedure for a many body system [32]. Figure 2 displays the contour plot of the mean gap ratio r in the W − U 1 phase diagram, namely, as a function of the disorder and of the long-range interaction strength. The colour code refers to the calculations performed for a gas of N = 8 fermions in a lattice with K = 16 sites, the gap ratio was first determined for 500 eigenvalues E n for which n = (E n − E min )/(E max − E min ) ≈ 0.5, where E max , E min are respectively the largest and the smallest eigenvalue for given disorder realization, and then averaged over 400 disorder realizations. One clearly identifies two regions: (i) the yellow region, corresponding to r ≈ 0.53 where the system has GOE level statistics and is thus ergodic, and (ii) the blue region with r ≈ 0.39, where the system is MBL. The white stripe separates the ETH from the MBL regimes and gives the disorder strength at which r = 0.45. This disorder strength depends on the system size and shifts to larger values as we increase the system size K. Nevertheless, a finite-size scaling analysis suggests that the ergodic-MBL boundary converges to the yellow line connecting the red markers. The red markers are obtained as follows. We first consider the scaling form of the disorder strength given by where the critical disorder strength W C (U 1 ) and the exponent ν(U 1 ) depend on the long-range interaction strength U 1 . We then consider the system sizes K = 16, 18, 20, which can be numerically simulated using the shift-invert technique, Ref. [33], implemented in PETSc/SLEPc setting, see Refs. [34,35]. The onset of Fig. 3 displays the mean gap ratior in the band center as a function of disorder strength for U 1 = 1 and U 1 = 4. The scaling (7) allows us to collapse the mean gap ratio for different system sizes r as a function of disorder strength W onto universal curves g U 1 [(W − W C (U 1 )) K 1/ν ] with good accuracy. From these curves we extract the critical disorder strengths W C (U 1 ) for U 1 ∈ {0, 0.5, 1, 2, 3, 4}, which correspond to the markers in Fig. 2. From this ansatz we also determine the exponent ν(U 1 ). This increases with U 1 from the value ν(U 1 = 1) = 1.0(1) to ν(U 1 = 4) = 1.5 (1), showing that the observed ETH-MBL transition gets increasingly sharp with the long range interaction strength U 1 .
We remark that the scaling ansatz (7) is analogous to the one used for the standard MBL system with short-range interactions [36]. The fact that the same scaling seems to hold even in presence of long-range interactions suggests that the underlying physics of our system is similar. The mean gap ratio r, in the center of the spectrum, is displayed as a function of the disorder amplitude W and for K = 16, 18, 20 lattice sites. The left panel corresponds to U 1 = 1 , the right one to U 1 = 4. The insets display the data rescaled according to Eq. (7). Error bars are smaller than the marker's size. The universal functions Let us note that the considered system possesses another non-ergodic phase at large values of the long range interaction strength U 1 10 (not shown in Fig. 2). This regime emerges when the all-to-all coupling term dominates in the Hamiltonian leading to non-ergodic dynamics associated with intrinsic non-ergodicity of systems with long-range interactions [37].

Dynamics of the system
Imagine we prepare the system in a well defined separable state |ψ 0 . To probe the dynamics of the system we consider the time-dependent density correlation function wheren is the average number of particles andn i (t) = ψ(t)|n i |ψ(t) is evaluated over the evolved state |ψ(t) = exp(−iHt)|ψ 0 . Here, the constant D warrants that C(0) = 1. According to ETH, an ergodic system loses the memory of the initial state and the C(t) correlation decays to zero. Conversely, in the MBL phase the density correlation function C(t) reaches a nonzero asymptotic value after a transient of the order of few J −1 [36]. The second quantity with which we probe the dynamics is the bipartite entanglement entropy S(t). This is obtained after splitting the lattice into two subsystems A and B and calculating the density matrix of the subsystem A ρ(t) = Tr B |ψ(t) ψ(t)|, where Tr B denotes the trace over subsystem B's degrees of freedom. The entanglement entropy is then defined as where ρ ii are Schmidt basis coefficients squared with i ρ ii = 1 (see e.g. [38]). In systems with short-range interactions the logarithmic growth of the entanglement entropy S(t) during time evolution of the system is a hallmark of MBL [39,40] and can be understood within the picture of LIOMs [41,42]. Systems with strong long-range interactions, on the other hand, manifest dynamical properties typical of ergodicity breaking, such as the logarithmic growth of the entanglement entropy after quenches even in absence of disorder [37,[43][44][45]. In order to single out the onset of localization and the effect of long-range interactions on the localization properties, in the following subsections we explore the transition non-ergodic to ergodic as a function of U 1 and for constant disorder amplitude W = 8. We then turn to the regimewhere the long-range interactions are a perturbation to the dynamics.

Density correlations and entanglement entropy
By inspecting the phase diagram ( Fig. 2), one can see that for the disorder strength W = 8 the system is deep in the MBL phase for U 1 = 0. In fact, at U 1 = 0 Eq. (1) for fermions reduces to the XXZ Heisenberg spin chain and undergoes the ETH-MBL transition in the vicinity of W C (U 1 = 0) = 3.7 [36]. This transition is accompanied by the appearance of non-vanishing values of the correlation function at the asymptotics, C(t → ∞) = 0, as well as by the logarithmic growth of the entanglement entropy S(t). We now consider a nonvanishing value of the long-range interaction strength, and in particular analyse the cases (i) U 1 = 1 that correspond to a MBL phase for finite K as well as in the thermodynamic limit (as obtained by the finite-size scaling, see Fig. 2) (ii) U 1 = 3, and (iii) U 1 = 5. The latter two cases are both in the localized regime for K = 16, and yet they are delocalized in the thermodynamic limit according to the finite size scaling in Fig. 2.
The time evolution of the correlation function C(t) as well as of the entanglement entropy S(t) is presented in Fig. 4. We consider system sizes K = 16, 18, 20 and refer the reader to Appendix A for details on the numerical methods we use. For U 1 = 1 we observe the features characteristic of ergodicity breaking: the correlation function C(t) acquires a stationary value which very weakly depends on the system size, as visible in Fig. 4(a). The entanglement entropy S(t) shows an increase with time which is sublinear, and indeed seems weaker than logarithmic, see Fig. 4(b). As the strength of long-range interactions U 1 increases there appears a slow decay of the correlations C(t) towards zero which becomes more pronounced as the system size K is increased. This result suggest that the correlations C(t) decay to zero in the thermodynamic limit, which is in agreement with the results of finite-size scaling. On the other hand the entanglement entropy S(t) for U 1 = 3, 5 clearly grows logarithmically in time. Such a behavior is consistent with the picture of LIOMs and is believed to be a feature of MBL system, which seems to lead to an apparent paradox: In fact, while the dynamics of C(t) suggest that large systems would be ergodic, at the same time, the entanglement entropy growth S(t) shows no signs of delocalized behaviour. Yet the behaviour of S(t) could also originate from the long-range nature of the interactions. We observe, in particular, that the slope of S(t) increases with U 1 and with the system size.  S P (t) In order to gain insight we analyse in detail the behaviour of the entanglement entropy. Following Ref. [46] we express the entanglement entropy as the sum of two contributions S(t) = S P (t) + S C (t), where S P (t) stems from particle number fluctuations between subsystems A and B and S C (t) is the entanglement entropy of different configurations of particles within the two subsystems. Let us denote by p n the probability of populating the n-particle sector in subsystem A and by ρ (n) the corresponding block of the density matrix ρ for subsystem A, such that ρ = n p n ρ (n) . A simple manipulation shows that Eq. (9) can be rewritten as The resulting behaviors of S P (t) and S C (t) are shown in Fig. 5. We first notice that exchange of particles between subsystems A and B occurs due to tunneling. As visible in subplot (a), S P (t) grows significantly at a time scale of few tunneling times, independently of the value of U 1 and of the system size. After this transient, its behaviour depends significantly on U 1 and on the system size. In particular, for U 1 = 1 it grows very slowly with time and weakly depends on the system size, hinting towards a strong suppression of particle number fluctuations. For U 1 = 5, instead, it has a clear logarithmic growth in time and a significant dependence on the system size. The former case is a standard MBL behavior [46]: the logarithmic growth of S(t) is mainly due to the increase in the configuration entropy S C (t). The latter behavior, in which S P (t) grows logarithmically in time enabling also faster and faster growth of S C (t) leads eventually to thermalization. The dynamics at large U 1 thus points towards ergodicity for larger sizes, in agreement with the finite-size scaling analysis. Yet, it is so slow that the decay time of the correlation function C(t) is much slower than in the ergodic regime at small W and U 1 = 0. We remark that deeply in the MBL phase of standard models with short-range interactions the time evolution can be efficiently simulated to large times (≈ 10 3 ) and for large systems sizes (K ≈ 10 3 ) [47] with algorithms based on tensor network approach. Such an approach is ruled out by the infinite range of interactions in our model.

Weak long-range interactions
We now analyse the role of the long-range interactions by considering the limit in which U 1 is sufficiently weak with respect to the tunneling rate. We first focus on the case U = 0, when the particles solely interact via the long-range interactions. For U 1 J we expect that the dynamics is first dominated by the hopping, and is then affected by U 1 on a longer time scale. Figure 6(a) displays the entanglement entropy S(t) for different values of U 1 , ranging from 1.6 × 10 −3 till 0.8. The entanglement entropy first rapidly grows during an initial transient, which is of the order of 1/J and is the same for all considered values of U 1 . After this transient S(t) saturates to a value for a time interval, up to a time scale T 1 of the order of T 1 ≈ 1/U 1 . We understand this behaviour as the system being in the Anderson localization regime, since for this time scale the dynamics is the one of non-interacting particles. After T 1 the long-range interactions start to significantly affect the dynamics and the entanglement entropy S(t) increases approximately linearly with time, till it reaches its saturation value. This saturation value depends on the given system size, on the strength of disorder W , and on the long-range interaction strength U 1 . Figure 6(b) displays the behaviour of the entanglement entropy S(t) but now for U = 1, thus when the interaction U is equal to the hopping rate. Thanks to the time scale separation the growth of entanglement entropy can be now understood as the combination of two effects: the growth which goes logarithmic in time, as for a standard MBL system, and the rapid transient of S(t) at the time scale T 1 due to the all-to-all coupling. The saturation values of the entanglement entropy seems to be only weakly dependent on the value of U . This leads us to conjecture that the origin of the nonergodic dynamics here is also MBL, and it arises from a quasi-degenerate manifold of states coupled by the long-range interactions.  In order to test this conjecture, in Fig. 7 we separately plot the particle and the configuration entanglement entropy as a function of time and corresponding to the curves with U 1 = 0.016 and U = 0, 1 of Fig. 6. For both values of U the particle entanglement entropy S P (t) increases rapidly as the initial occupation of lattice sites spreads due to tunnelling. This transient dynamics occurs on a time scale of the order of 1/J, after which the particle number fluctuations change only marginally: observe that the associated S P (t) is approximately independent of the considered system sizes. Instead, the configuration entanglement entropy S C (t) shows different behaviours as a function of U and of the system size. The configuration entanglement entropy S C (t) shows also a dynamics that can be separated by the time scale T 1 : in the first part for U = 0 it is roughly constant, while for U = 1 it grows logarithmically with time. After T 1 it rapidly grows and then saturates to a value which increases with system size. We understand this increase as the number of accessible configurations grows with K. On the basis of this analysis we conclude that the dynamics for U = 0 and small U 1 is a textbook case of Anderson localization perturbed by weak long-range interactions that couple only states that are closely spaced in energy. Such states are localized in different regions of space. Thus long-range interactions leads to a spread of the initial state among relatively few localized eigenstates, the corresponding dynamics is strongly nonergodic. The analogies with the dynamics observed in the case U = 1 and small U 1 leads us to conclude that the dynamics is MBL, where the long-range interactions strongly couple the quasi-degenerate manifold of localized states, while the dynamics remains nonergodic.
We explore the properties of this peculiar MBL phase by using a LIOM picture. We first recall that the Hamiltonian of a generic (fully) many-body localized system can be expressed as [41,42] where τ z i are quasi-local operators knows as LIOMs or l-bits. They can be thought of as dressed occupation number operators as τ z i = U † n i U where U is a quasi local unitary transformation. The couplings J ij fall off exponentially with distance between interacting l-bits (similar holds for higher order couplings J ijk , ...), ξ is the localization length in the system. It has been analytically shown [48] that the l-bit model (11) leads to a logarithmic growth of the Renyi-2 entropy S 2 (t) = − log Trρ 2 (which for large times and large system sizes behaves essentially the same as the von Neumann entanglement entropy S(t)) assuming that the initial state is an equal superposition of all Fock states. Moreover, to observe the logarithmic growth of the Renyi entropy S 2 (t) is suffices to introduce only the J (2) ij couplings (12) keeping the higher order couplings J ijk , ... equal to zero. In Fig. 8 we show that the Renyi entropy of the l-bit Hamiltonian (11) reproduces the behaviour of the configuration entropy S C (t) for the extended Hubbard model with spinless fermions. Specifically, assuming the exponential decay of the coupling terms J (2) ij , Eq. (12), we reproduce the logarithmic growth of entanglement entropy characteristic of standard MBL. In order to reproduce the rapid growth of S C (t), we introduce long-range couplings between LIOMs according toJ where the term J 1 /K(−1) i+j r ij mimics the coupling experienced by l-bits caused by the longrange interaction term H C (5) and r ij ∈ [0, 1] is a random variable which models the overlap of τ z i and n i . If we set J 0 = J 1 = 0 in Eq. (13), then there is no growth of entanglement entropy in the l-bit model which corresponds to no growth of configuration entropy -a situation of Anderson localization. Introducing non-zero J 0 in the model (11) we obtain the logarithmic growth of S 2 (t) -the hallmark of MBL. Setting then J 1 to a finite, non-vanishing value one gets a rapid growth of entanglement entropy which starts at a certain time scale set by J 1 : After this growth S 2 (t) saturates at the same value as in the case J 1 = 0. All of those feature are in qualitative agreement with the growth of configuration entropy S C (t) for our system at U 1 = 0.016 and U = 1. The analysis of this Section strictly applies for very small values of U 1 , such that the time scale T 1 = 1/U 1 is much larger than the time scale J −1 set by the tunneling rate. As visible in Fig. 6, this separation of time scales takes place as long as U 1 0.16. The saturation value of the entanglement entropy S(t) is already slightly larger for U 1 = 0.16 than for smaller values of U 1 meaning that a slight modification of structure of LIOMs (possibly an increase of the support of τ z i ) happened. For larger U 1 = 0.8, the saturation value of S(t) is significantly larger. On the basis of the discussion so far, in particular of the studies in Sec. 3, we conclude that the system is still MBL, however, the properties of l-bits τ z i are significantly affected by the long-range interactions.
We have shown, therefore, that the cavity mediated long-range interactions lead to nonergodic behavior of the system, which in presence of strong disorder can be understood as MBL. The observed ergodicity breaking exhibits novel features, such as the rapid growth of entanglement entropy due to the long-range interactions. When the effect of the long-range interactions can be separated from the short-range interactions, then the entanglement entropy increases in good approximation linearly in time. Interestingly the entanglement entropy is still bounded by a constant which only moderately changes with U 1 as long as the system is in the localized phase. The observed phenomenon is different from the non-ergodic phase observed for systems with single particle mobility edge [49][50][51][52] where long-range interactions are not involved.

Detection of ergodicity breaking: Light at the cavity-output
The dynamics of the system can be revealed by experimentally measuring the correlation function C(t) via a site-resolved measurement of n i , which can be performed in cold atoms experiments [7]. In this section we argue that the cavity setup can allow one to measure the breakdown of ergodicity via the light at the cavity output. This measurement is nondestructive, the light which is emitted is scattered by the atoms and thus contains the information of their density distribution, the corresponding observable is the intensity I(t) which reads [24,53,54] where I(t) = i (−1) i n i is the population imbalance between even and odd sites of the cavity standing-wave mode, namely, the sites x i = ia where the standing wave mode with wave number k takes the value cos kx i = +1 and −1, respectively. Recall that a similar imbalance I(t) was employed in Ref. [7] to demonstrate MBL for a system of fermions in the optical disordered lattice. Below we discuss the time evolution of I(t) for two possible states in which the system can be initially prepared. We first consider a density-wave like state |DW 10 = |101010.. , with odd sites occupied and even sites empty. This state maximizes the population imbalance I for fractional densityn = 1/2 and minimizes the short-range interactions. In the absence of disorder, it is an eigenstate of the Hamiltonian in the atomic limit [21]. In the presence of strong disorder and finite hopping we expect that I(t) decays to zero in the ergodic phase, while it can saturate to some constant finite value in the MBL phase. In this analysis one shall keep in mind that the state |DW 10 becomes the ground state for U 1 W , and its energy thus gets closer to the ground state as the ratio U 1 /W is increased. Figure 9 diplays the time evolution of I 2 (t) for the system initially prepared in |DW 10 state. The dynamics is shown for W = 0.5 (left panel) and W = 10 (right panel), corresponding to the ergodic and the MBL regimes, respectively, for the considered values of U 1 . For comparison we show the dynamics of the correlation function C(t) for the same parameters, but when the initial state is a random Fock state with the same density. In the ergodic regime (W = 0.5), the correlation function C(t) decays to zero signifying total relaxation of the initial density profile regardless of the value of U 1 . The dynamics of the population imbalance for the initial state |DW 10 depends strongly on the ratio U 1 /W . For W = 0.5 and U 1 = 0.16, the intensity I 2 (t) decays to zero and is a valid probe of the ergodic properties of the system. For the same disorder amplitude and larger U 1 = 1.6, instead, the intensity I 2 (t) saturates at a constant value, from which we infer that the state |DW 10 has already significant overlap with the ground state of the system. For strong disorder, W = 10, the information about localization properties of the system provided by the correlation function C(t) becomes also visible in time evolution of I 2 (t), even for U 1 = 6.4. This behavior is due to the fact that the energy spectrum of the system is much broader, such that even for values of the interaction strength such as U 1 = 6.4 the state |DW 10 is still in a region of the spectrum with a relatively large density of states. We now assume that the system has been prepared in the state |DW 2 = |11001100... , where the initial site occupations form a pattern of intertwining doubly occupied and double empty lattice sites. State |DW 2 is an eigenstate of the imbalance operator I to the eigenvalue 0, therefore it has energy lying around the band center. Figure 10 displays the time evolution of the light intensity at the cavity output starting from |DW 2 = |11001100... at t = 0 for three different values of the disorder amplitude, W = 1, 3, 7. The solid lines display for comparison I(t) for U 1 = 0, where W = 1 is in the ergodic case and W = 3, 7 are MBL. The dashed-dotted lines give the corresponding behaviour of the population imbalance for U 1 = 2, where according to the phase diagram, Fig. 2, the phase is expected to be ETH for W = 1, 3 and MBL for W = 7 at K = 16. When the system is ergodic, after a short transient the density profile relaxes to uniform density, which correspond to I(t) 2 → 0 and is observed for W = 1 and for both values U 1 = 0 and U 1 = 2. On the other hand, in the situation of many-body localized system, the way in which the imbalance I(t) evolves in time from the |DW 2 state depends strongly on the disorder realization. This means that the square of the imbalance I 2 (t) is expected to saturate at a non-zero value, which is what we observe in Fig. 10 for W = 7.
We have therefore shown that the measurement of light emitted by the cavity can allow one to determine the phase of the system by starting from well defined states. This probe of of ergodicity breaking is an appealing alternative to the band mapping technique of Refs. [55,56] used in standard population imbalance measurements.

Quasi-random disorder
The model we discuss can be realized with the setups of Ref. [24,57] by introducing additional additional weak beams with wavelength incommensurate with the cavity lattice periodicity, which create quasi-random disorder analogously to the experiment of Ref. [7]. A second possibility would be to add a second, off-resonant pumping laser with a random intensity distribution. This laser shall drive to an atomic transition that does not scatter into the cavity field, and thus generate a random a.c. Stark shift of the onsite energy. We now analyse the predictions on MBL when the onsite energy is due to the contribution of an incommensurate periodic potential, namely, Here, W plays a role of disorder amplitude, β is the ratio of the two lattice periodicities (we take β = ( (5) − 1)/2) and the value of the phase φ ∈ [0, 2π] determines the disorder realization. Even though there are long-range correlations between on-site energies, the quasiperiodic potential (15) was employed in a number of experiments with MBL in optical lattices [7,52,58,59]. It is well known that the properties of the MBL transition for purely random disorder differ substantially from the ones of quasi-periodic disorder, as it was shown by the analysis of the entanglement entropy [60] and of the gap ratio [29]. Yet, while the important aspects of the transition itself are different, the ergodic and MBL phases are similar in the two settings, which is the reason why the seminal observation of MBL [7] was feasible in a setup with quasiperiodic disorder. In order to show that the similar property is shared in our extended Hubbard model, we here discuss the dynamics of spinless fermions when the onsite energy in Eq. (2) is given by Eq. (15). The time evolution of the correlation function C(t) and of the entanglement entropy S(t) are shown in Fig. 11. Two distinct regimes can be here identified: (i) the ergodic phase at W = 1 in which the correlations C(t) rapidly decays to zero and entanglement spreads ballistically and (ii) the MBL phase at W = 8 characterized by an asymptotic non-vanishing value of the correlation function C(t) as well as by the logarithmic growth of entanglement entropy. The intermediate disorder strength W = 4 corresponds to the regime in which the   Figure 12: Configurational entanglement entropy S C (t) and particle entanglement entropy S P (t) for the system of N = 7 spinless fermions on K = 14 sites with quasirandom disorder. The parameters are U = 0 and W = 8. The results averaged over the phase φ in the interval [0, 2π], the initial state is a Fock state with random site occupation and N = 7.

S(t)
localized-to-ergodic transition takes place as the strength of all-to-all coupling U 1 is increased.
To provide further evidence that the physics is similar for both random and quasiperiodic disorders, we calculate the time evolution of the configurational entanglement entropy S C (t) and of the particle entanglement entropy S P (t). The results are presented in Fig. 12. Similarly to the case of random disorder, small values of U 1 lead to rapid growth of S C (t) at the time scale T 1 = 1/U 1 . For U 1 = 0.16, S C (t) reaches a larger asymptotic value. Further increase of the long-range interaction strength U 1 leads to the logarithmic growth of S C (t) with time, the entanglement entropy saturates at much larger value. We thus predict that the properties of the MBL phase in the system of spinless fermions in the cavity with quasiperiodic disorder are analogous to the ones discussed in presence of random disorder.

Bosons in the cavity
Many-body localization in bosonic systems was studied numerically in [9,10,61] as well as in experimental realizations [46,62]. The essential features of the MBL phase in bosonic models are analogous to the system of spinless fermions (or spins). However, the possibility that the lattice site occupations exceeds unity leads to the natural appearance of many-body mobility edges at higher energy densities. Motivated by the fact that several experiments with cold atoms in optical resonators [24,57] are performed with bosons, here we briefly discuss signatures of the MBL in the extended Bose-Hubbard Hamiltonian: Figure 13 displays the phase diagram for a lattice consisting of K = 8 sites at unity filling. The diagram is obtained from the gap ratio r averaged over states at the center of the spectrum ≈ 0.5 and shows that for a given value of the all-to-all coupling U 1 there exists a disorder strength sufficient to localize the system. The result is analogous to the phase diagram of spinless fermions in Fig. 2. We note however, that the values of disorder needed to induce many-body localization in the bosonic system (16) are larger than the ones required for the fermionic counterpart.  To provide further insight into physics of the MBL phase in bosonic system in the presence of long-range interactions we calculate the bipartite entanglement entropy during the course of time evolution of the system -c.f. Fig. 14. The features are overall analogous to the spinless fermions and thus seem to be not limited to small Hilbert spaces per site as for fermions.

Conclusion
In this work we have analysed the occurrence of many-body localization in a system of particles with all-to-all interactions, whose dynamics is described by an extended Hubbard models, where the onsite energy follows a random distribution. Our study is numerical and is based on exact diagonalization as well as on methods for sparse Hamiltonian matrices. By means of finite-size scaling we have shown that the MBL phase is indeed present in the system and the transition to ergodic phase occurs at disorder amplitudes whose critical value increases with the long-range interaction strength. The all-to-all interactions affect spreading of the entanglement in the system. Entanglement entropy growth due to long-range interactions in MBL phase is faster than logarithmic -it is linear in cases where the tunneling coefficient is much larger than long-range interaction strength. However, the saturation value of entanglement entropy is constrained by particle number fluctuations between parts of the system. This vanishes even in presence of all-to-all interactions and hence, the dynamics is non-ergodic. We have shown that the features of the MBL phase can be qualitatively understood within picture of LIOMs with long-range couplings. The interplay between disorder and long-range interactions results in transition from the MBL to an ergodic phase. The characteristic features of the observed ETH-MBL transition seem to be independent of the quantum statistics, as the numerical analysis for finite systems of spinless fermions and of bosons show. Analogous features are observed both for true random disorder as well as for quasi-periodic potential. This dynamics can be observed experimentally with quantum gases in an optical resonator. For these systems, indeed, the light intensity at the cavity output can probe the gas phase and provide evidence of ergodicity breaking.
We finally remark that MBL phase in presence of global coupling to d dimensional system was found in [63], particular scaling of d with system size leads in high frequency limit to Hamiltonian similar to the one considered in this work.

A Appendix: Time evolution with Chebyshev expansion technique
Time evolution of fermionic systems at half filling with K = 16 (and smaller) can be obtained easily by full exact diagonalization of the Hamiltonian matrix followed by exact calculation of the evolution operator U (t) for arbitrary time t.
To deal with larger system sizes we employ the expansion of the evolution operator into series involving Chebyshev polynomials [64,65] where a = (E max − E min )/2, b = (E max + E min )/2, the Hamiltonian is rescaled H = 1 a (H − b) so that spectrum of H belongs to the [−1, 1] interval, J k (t) is the Bessel function of the order k and T k (x) is the Chebyshev polynomial of order k. The number of terms N needed to assure convergence of the expansion (17) to time t max is N ≈ 2at max [66].
The time-evolution of the initial state |ψ 0 is given by and reduces to matrix-vector multiplications where the recursion relation satisfied by Chebyshev polynomials was used. In order to get |ψ(t) we generate iteratively a sequence of N vectors |ψ 0 , T 1 |ψ 0 , ..., T N |ψ 0 . To reach long times of time evolution t max ≈ 10 3 one needs relatively large N which increases memory consumption. Therefore we split the time interval [0, t max ] into parts [0, ∆t], [∆t, 2∆t], ... in such a way that |ψ ((n + 1)∆t) can be calculated from the state |ψ (n∆t) with the expansion (18) involving only a limited number of terms e.g. -N ≈ 1000 which allows us to obtain time evolution for the system size K = 20 with memory consumption smaller than 5GB (performing the matrix-vector multiplications in PETSc).