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

TL;DR

This paper demonstrates that a 1D Hubbard model with cavity-mediated long-range interactions can exhibit many-body localization, observable through dynamics of emitted light, expanding understanding of localization phenomena in long-range interacting systems.

Contribution

It introduces a model showing many-body localization with cavity-mediated long-range interactions and provides numerical evidence and observable signatures.

Findings

Localization occurs in a 1D Hubbard model with all-to-all coupling.

Time evolution shows features consistent with many-body localization.

Localization can be observed via light emitted from optical cavities.

Abstract

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.

Figures

Equations

$$H=H_{\mathrm{A}}+H_{\mathrm{C}},$$

$$H_{\mathrm{A}}=-J\sum _j^K\left(b_{j+1}^{†}{b}_j+H.c.\right)+\sum _j^K{E}_j{n}_j+H_{\mathrm{int}}^{\mathrm{F},\mathrm{B}}.$$

$$H_{\text{int }}^{\mathrm{F}}=U\sum _j^K{n}_j{n}_{j+1},$$

$$H_{\text{int }}^{\mathrm{B}}=U\sum _j^K{n}_j(n_j-1).$$

$$H_{\mathrm{C}}=-\frac{U_1}{K}{\left(\sum _j^K(-1{)}^j{n}_j\right)}^2=-\frac{U_1}{K}\sum _{i,j}^K(-1{)}^{i+j}{n}_i{n}_j,$$

$$I(t)=\sum _i(-1{)}^i{n}_i$$

$$r_n=\min ({\delta }_n,{\delta }_{n+1})/\max ({\delta }_n,{\delta }_{n+1})$$

$$W\to \left(W-W_C(U_1)\right)K^{1/\nu (U_1)},$$

$$C(t)=D\sum _{i=1}^K\left({\bar{n}}_i(t)-\bar{n}\right)\left({\bar{n}}_i(0)-\bar{n}\right),$$

$$S=-\sum _i{\rho }_{ii}(t)\log ({\rho }_{ii}(t)),$$

$$S(t)=-\sum _{n=0}^N{p}_n\log (p_n)-\sum _{n=0}^N{p}_n\sum _i{\rho }_{ii}^{(n)}\log ({\rho }_{ii}^{(n)})\equiv S_P(t)+S_C(t).$$

$$H=\sum _{i=1}^K{J}_i^{(1)}{\tau }_i^z+\sum _{i,j=1}^K{J}_{ij}^{(2)}{\tau }_i^z{\tau }_j^z+\sum _{i,j,k=1}^K{J}_{ijk}^{(3)}{\tau }_i^z{\tau }_j^z{\tau }_k^z+...,$$

$$J_{ij}^{(2)}=J_0{\mathrm{e}}^{-|i-j|/\xi },$$

$${\tilde{J}}_{ij}^{(2)}=J_0{\mathrm{e}}^{-|i-j|/\xi }+\frac{{\mathrm{J}}_1}{\mathrm{K}}(-1{)}^{\mathrm{i}+\mathrm{j}}{\mathrm{r}}_{\mathrm{ij}},$$

$$E_{\mathrm{out}}(t)\propto ⟨I(t)⟩,$$

$$E_j=W\cos (2\pi \beta j+\varphi ).$$

$$H=-J\sum _j^K\left(b_{j+1}^{†}{b}_j+H.c.\right)+\sum _j^K{E}_j{n}_j+U\sum _j^K{n}_j(n_j-1)-\frac{U_1}{K}\sum _{i,j}^K(-1{)}^{i+j}{n}_i{n}_j.$$

$${\hat{\mathcal{H}}}_A=\sum _{j=1,2}\int dx{\hat{\Psi }}_j^{†}(x){\hat{H}}_j(x){\hat{\Psi }}_j(x)+2U_{12}\int dx{\hat{\Psi }}_1^{†}(x){\hat{\Psi }}_2^{†}(x){\hat{\Psi }}_2(x){\hat{\Psi }}_1(x),$$

$${\hat{H}}_j(x)=-\frac{{\hslash }^2{\partial }_x^2}{2m}+V_{\mathrm{cl}}^{(j)}{\cos }^2(kx)+\frac{1}{2}\int d{x}^{\prime }{U}_{jj}(x-x^{\prime }){\hat{\Psi }}_j^{†}(x^{\prime }){\hat{\Psi }}_j(x^{\prime })-\hslash {\Delta }_a|2⟩⟨2|.$$

$${\hat{\mathcal{H}}}_C=-\hslash {\Delta }_C{\hat{a}}^{†}\hat{a},$$

$${\hat{\mathcal{H}}}_{\mathrm{int}}=\hslash \int dx{\hat{\Psi }}_2^{†}(x){\hat{\Psi }}_1(x)(g(x)\hat{a}+{\Omega }_z(x))+H.c.,$$

$${\hat{\Psi }}_2(x,t)\simeq$$

$$\dot{\hat{\Psi }}=$$

$$\dot{\hat{a}}=$$

$$\hat{n}(x)={\hat{\Psi }}^{†}(x)\hat{\Psi }(x).$$

$${\int }_{\bar{t}}^{\bar{t}+\Delta t}\hat{a}(\tau )d\tau /\Delta t={\hat{a}}_{\mathrm{st}}(\bar{t}),$$

$${\hat{a}}_{\mathrm{st}}(\bar{t})=\frac{\int dxS(x)\hat{n}(x,\bar{t})}{({\Delta }_C-\int dx{U}_0(x)\hat{n}(x,\bar{t}))+\mathrm{i}\kappa }+\frac{\mathrm{i}\sqrt{2\kappa }{\hat{\bar{a}}}_{\mathrm{in}}(\bar{t})}{({\Delta }_C-\int dx{U}_0(x)\hat{n}(x,\bar{t}))+\mathrm{i}\kappa },$$

$${\hat{\bar{a}}}_{\mathrm{out}}(t)=\sqrt{2\kappa }{\hat{a}}_{\mathrm{st}}-{\hat{\bar{a}}}_{\mathrm{in}},$$

$$\dot{\hat{\Psi }}=$$

$$H_{\mathrm{CQED}}=\hslash \int dx\frac{{\Omega }_z(x{)}^2}{{\Delta }_a}\hat{n}(x)+\frac{1}{{\Delta }_C}\hslash {\left(\int dxS(x)\hat{n}(x)\right)}^2,$$

Full text

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

Piotr Sierant1, Krzysztof Biedroń1, Giovanna Morigi2, and Jakub Zakrzewski1,3

1 Instytut Fizyki imienia Mariana Smoluchowskiego, Uniwersytet Jagielloński, Łojasiewicza 11, 30-348 Kraków, Poland

2 Theoretische Physik, Universität des Saarlandes, D-66123 Saarbrücken, Germany

3 Mark Kac Complex Systems Research Center, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland

• piotr.sierant@uj.edu.pl

Abstract

**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. **

Contents

1. 1 Introduction

2. 2 Description of the system

3. 3 Phase diagram

4. 4 Dynamics of the system

5. 4.1 Density correlations and entanglement entropy

6. 4.2 Weak long-range interactions

7. 5 Detection of ergodicity breaking: Light at the cavity-output

8. 6 Quasi-random disorder

9. 7 Bosons in the cavity

10. 8 Conclusions

11. A Derivation of the Hubbard Hamiltonian

12. B Nonergodic regime for strong all-to-all interactions

13. C Time evolution with Chebyshev expansion technique

1 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 the tight–binding limit [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 other studies [21, 22, 23, 24] a model of a single spin coupled to all sites of otherwise short-range interacting systems was analysed showing that the presence of additional coupling may strongly modify MBL properties.

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. Such a model has been extensively studied in the past concentrating mainly, however, on ground state properties [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37] (for a review with extensive list of references to earlier works see [38]. The certainly incomplete list of recent experimental works in that area includes [39, 40, 41, 42]). 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. When the atoms are tightly confined by an external optical lattice, the cavity-mediated long-range interactions tends to order the atoms in structures maximizing the intracavity field intensity. We investigate MBL in this extended Hubbard model and with local disorder using exact diagonalization supplemented by numerical techniques for sparse Hamiltonian matrices for a gas of fermions (or bosons). Our results show features which can be attributed to the occurence of MBL in the system. We argue that these features can be revealed in the light emitted by the resonator.

The paper is structured as follows. In Sec. 2 we describe an experimentally realizable system of atoms inside a resonant cavity. We summarize the details of the derivation of the effective model in Appendix A, as described e.g. in [28, 33]. 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. Appendix B contains a short discussion of the nonergodicity in the specific limit of very strong cavity-mediated coupling which is outside the MBL regime; Appendix C provides the details of the time-propagation algorithm used.

2 Description of the system

We consider an ensemble of atoms trapped in a quasi one-dimensional 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 (for the detailed derivation see Appendix A):

[TABLE]

with $H_{\rm A}$ being the standard Hubbard-like Hamiltonian for the dynamics of atoms in an optical lattice and in presence of disorder and $H_{\rm 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

[TABLE]

where $b_{j}$ and $b_{j}^{\dagger}$ are the onsite annihilation and creation operators of a fermion or a boson at site $j=1,\ldots,K$, (with $K+1$ identified with the first site), $J$ is the tunneling coefficient scaling the nearest-neighbour hopping, $n_{j}=b_{j}^{\dagger}b_{j}$ is the occupation operator at site $j$, $E_{j}$ is the onsite energy at site $j$, and $H^{\rm F,B}_{\rm int}$ is the interaction term, which takes different forms depending on the quantum statistics of the atomic gas. For spinless fermions

[TABLE]

and $U>0$, whereas the first non–trivial interaction term in tight-binding expansion for bosons reads

[TABLE]

In turn, the cavity-mediated long-range interactions take the form [28, 33, 37]:

[TABLE]

with $U_{1}>0$. This Hamiltonian term is here derived under the assumption that the wavelength of the cavity field equals the one of the electric field generating the optical lattice, see Appendix A. This interaction is proportional to the squared population imbalance, $H_{\rm C}\propto-I(t)^{2}$, where

[TABLE]

and the sign is due to the fact that the standing wave cavity mode with wave number $k$ takes value $\cos(kia)=(-1)^{i}$ at the optical lattice site centered at $x_{i}=ia$. Its expectation value is thus positive (negative) when the even (odd) sites of the cavity standing-wave mode are prevailingly occupied and it favours density-wave (DW) ordering [28, 33, 34, 35]. We determine the existence of the MBL phase by an exact diagonalization of the Hubbard Hamiltonian, Eq. (1), for spinless fermions. 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$.

3 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) [43]. As the disorder strength increases and the system becomes localized, the level statistics becomes Poissonian [44, 45] (an accurate model for level statistics across the localization transition was recently proposed in Refs. [46, 47]). The level statistics can be characterized using the gap ratio. This is defined as

[TABLE]

with $\delta_{n}=E_{n}-E_{n-1}$ being the spacing between two consecutive eigenvalues [48]. Averaging over different energy levels within a certain interval as well as over disorder realizations results in the mean gap ratio, $\overline{r}$, that may be used to characterize the spectra. The mean gap ratio changes from $\overline{r}\approx 0.53$ in the ergodic regime [48, 49] to $\overline{r}\approx 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 [50].

Figure 2 displays the contour plot of the mean gap ratio $\overline{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 $\epsilon_{n}=(E_{n}-E_{\rm min})/(E_{\rm max}-E_{\rm min})\approx 0.5$, where $E_{\rm max},$ $E_{\rm min}$ are respectively the largest and the smallest eigenvalue for given disorder realization, and then averaged over $400$ disorder realizations. The statistics is sufficient to determine the mean gap ratio with an accuracy below 1% of its value. One clearly identifies two regions: (i) the yellow region, corresponding to $\overline{r}\approx 0.53$ where the system has GOE level statistics and is thus ergodic, and (ii) the blue region with $\overline{r}\approx 0.39$, where the system is MBL. The white stripe separates the ETH from the MBL regimes and gives the disorder strength at which $\overline{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

[TABLE]

where the critical disorder strength $W_{C}(U_{1})$ and the exponent $\nu(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. [51], implemented in Portable, Extensible Toolkit for Scientific Computation (PETSc) in Scalable Library for Eigenvalue Problem computations (SLEPc) setting, see Refs. [52, 53]. The onset of Fig. 3 displays the mean gap ratio $\bar{r}$ in the band center as a function of disorder strength for $U_{1}=1$ and $U_{1}=4$. The scaling (8) allows us to collapse the mean gap ratio for different system sizes $\overline{r}$ as a function of disorder strength $W$ onto universal curves $g_{U_{1}}[\left(W-W_{C}(U_{1})\right)K^{1/\nu}]$ with good accuracy. From these curves we extract the critical disorder strengths $W_{C}(U_{1})$ for $U_{1}\in\{0,0.5,1,2,3,4\}$, which correspond to the markers in Fig. 2. From this ansatz we also determine the exponent $\nu(U_{1})$. This increases with $U_{1}$ from the value $\nu(U_{1}=1)=1.3(1)$ to $\nu(U_{1}=4)=1.8(1)$. We remark that the scaling ansatz (8) is analogous to the one used for the standard MBL system with short-range interactions [54]. 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.

So far the gap ratio analysis together with the finite-size scaling indicates the existence of a boundary between ergodic (ETH-like) and MBL phase. The mean gap ratio reaches the value characteristic for poissonian ensemble i.e. for integrable systems. This occurs for sufficiently large disorder values, strongly dependent on the long range interaction strength $U_{1}$ which realizes all-to-all couplings.

For completeness, let us note that the considered system possesses another non-ergodic phase at large values of the long range interaction strength $U_{1}\gtrsim 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 due to global interactions [55]. We discuss this regime in Appendix B considering here the regime of small $U_{1}$.

4 Dynamics of the system

Imagine we prepare the system in a well defined separable state $|\psi_{0}\rangle$. To probe the dynamics of the system we consider the time-dependent density correlation function

[TABLE]

where $\bar{n}$ is the average number of particles and $\bar{n}_{i}(t)=\langle\psi(t)|n_{i}|\psi(t)\rangle$ is evaluated over the evolved state $|\psi(t)\rangle=\exp(-iHt)|\psi_{0}\rangle$. Here, the constant $D$ warrants that $C(0)=1$. According to ETH, an ergodic system loses the memory of the initial state and the correlation $C(t)$ decays to zero. Conversely, in the MBL phase the density correlation function $C(t)$ reaches a nonzero asymptotic value after a transient time of the order of few $J^{-1}$ (which is here set to unity) [54].

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 $\rho(t)$ of the subsystem $A$: $\rho(t)=\mathrm{Tr}_{B}\{|\psi(t)\rangle\langle\psi(t)|\}$, where $\mathrm{Tr}_{B}$ denotes the trace over subsystem B’s degrees of freedom. The entanglement entropy is then defined as

[TABLE]

where $\rho_{ii}$ are Schmidt basis coefficients squared with $\sum_{i}\rho_{ii}=1$ (see e.g. [56]). In systems with short-range interactions the logarithmic growth of the entanglement entropy $S(t)$ during the time evolution of the system is a hallmark of MBL [57, 58] and can be understood within the picture of LIOMs [59, 60]. 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 [61, 62, 63, 55]. 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 between non-ergodic and ergodic regimes as a function of $U_{1}$ and for constant disorder amplitude $W=8$. We then turn to the regime where the long-range interactions are a weak perturbation to the dynamics.

4.1 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$ [54]. This transition is accompanied by the appearance of non-vanishing values of the correlation function at the asymptotics, $C(t\rightarrow\infty)\neq 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$, (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 case $U_{1}=1$, instead, corresponds to a MBL phase for all system sizes we consider as well as in the thermodynamic limit (as predicted by the finite-size scaling, see Fig. 2).

The time evolution of the correlation function $C(t)$ as well as of the entanglement entropy $S(t)$ is displayed in Fig. 4.

We analyze system sizes $K=16,18,20$ using the Chebyshev expansion technique for the time evolution (see Appendix C for details). 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 (compare with the left panel in Fig. 4). The entanglement entropy $S(t)$ (right panel) shows an increase with time which is sublinear, and indeed seems weaker than logarithmic. 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$ and $U_{1}=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. This seems to lead to an apparent paradox: In fact, while the dynamics of $C(t)$ suggests 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 [55]. We observe, in particular, that the slope of $S(t)$ increases with $U_{1}$ and with the system size.

In order to gain insight we analyse in detail the behaviour of the entanglement entropy. Following Ref. [64] 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 $\rho^{(n)}$ the corresponding block of the density matrix $\rho$ for subsystem $A$, such that $\rho=\sum_{n}p_{n}\rho^{(n)}$. A simple manipulation shows that Eq. (10) can be rewritten as

[TABLE]

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 left panel of Fig. 5 $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 [64]: 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 in standard models with short-range interactions deeply in the MBL phase, the time evolution can be efficiently simulated to large times ($\approx 10^{3}$) and for large systems sizes ($K\approx 10^{3}$) [65] with time-dependent density matrix renormalization group related approaches [66]. Such approaches are ruled out by the infinite interaction range of our model. We note that algorithms based on time dependent variation principle [67, 68] could in principle tackle large system sizes. We leave this task for future work and here consider small sizes, amenable to diagonalization-like treatments.

4.2 Weak long-range interactions

We now analyse the role of the long-range interactions on MBL 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}\ll 1$ we expect that the dynamics is first dominated by the hopping, and only on a much longer time scale it is going to be visibly affected by $U_{1}$. The left panel of Fig. 6 displays the entanglement entropy $S(t)$ for different values of $U_{1}$, ranging from $1.6\times 10^{-3}$ till $0.8$. The entanglement entropy first rapidly grows during an initial transient, which is the same for all considered values of $U_{1}$ and is of the order of $1/J$ (which is here the unit of time). After this transient $S(t)$ saturates to a value for a time interval, up to a time scale $T_{1}\approx 1/U_{1}$. We understand this behavior 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 saturates. The corresponding saturation value depends on the given system size, on the strength of disorder $W$, and on the long-range interaction strength $U_{1}$, being, however, much lower, than the maximal (ergodic) value for a given system size.

The behaviour of the entanglement entropy $S(t)$ for $U=1$ (thus when $U=J$) is shown in the right panel of Fig. 6. The time scale separation allows us to identify two behaviours characterizing the entanglement growth: 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 interaction strength $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 of Fig. 6 with $U_{1}=0.016$ and $U=0$ and $U=1$. For both values of $U$ the particle entanglement entropy $S_{P}(t)$ increases rapidly as the initial occupation of the lattice sites spreads due to tunneling. 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 $K$. Instead, the configuration entanglement entropy $S_{C}(t)$ shows different behaviors as a function of $U$ and of $K$. The dynamics of $S_{C}(t)$ can be characterized by the time scale $T_{1}\sim 1/U_{1}$: for $t\lesssim T_{1}$ and for $U=0$ it is roughly constant, while for $U=1$ it grows logarithmically with time. After $T_{1}$ the configuration entanglement entropy rapidly grows and then saturates to a value which increases with the 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. The long-range interactions couple only states that are closely spaced in energy and that are localized in different regions of space. Thus long-range interactions lead to a spread of the initial state among relatively few localized eigenstates, the corresponding dynamics is strongly nonergodic. The analogies shared between the dynamics at $U=0$ and $U=1$ for small $U_{1}$ suggest that for $U\neq 0$ the dynamics is MBL-like. MBL is perturbed by the long-range interactions, which strongly couple the quasi-degenerate manifold of localized states and at the same time 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 [59, 60]

[TABLE]

where $\tau^{z}_{i}$ are quasi-local operators knows as LIOMs or l-bits. They can be thought of as dressed occupation number operators as $\tau^{z}_{i}=U^{{\dagger}}n_{i}U$ where $U$ is a quasi local unitary transformation. The couplings $J_{ij}$ fall off exponentially with the distance between interacting l-bits as

[TABLE]

where $\xi$ is the localization length (a similar relation holds for higher order couplings $J^{(3)}_{ijk},...$). It has been analytically shown [69] that the l-bit model (12) leads to a logarithmic growth of the Renyi-2 entropy $S_{2}(t)=-\log\rm Tr\rho^{2}$ (which for large times and large system sizes behaves essentially 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)$ it suffices to keep only the coupling coefficients $J^{(2)}_{ij}$, Eq. (13), setting the higher-order couplings $J^{(3)}_{ijk},...$ equal to zero.

In Fig. 8 we show that the Renyi entropy of the l-bit Hamiltonian (12) 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. (13), 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 to

[TABLE]

where the term $J_{1}/K(-1)^{i+j}r_{ij}$ mimics the coupling experienced by l-bits caused by the long-range interaction term $H_{C}$ (5) and $r_{ij}\in[0,1]$ is a random variable which models the overlap of $\tau^{z}_{i}$ and $n_{i}$. If we set $J_{0}=J_{1}=0$ in Eq. (14), then there is no growth of entanglement entropy in the l-bit model which corresponds to no growth of configuration entropy – a situation characteristic of Anderson localization. Introducing non-zero $J_{0}$ in the model (12) 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}\lesssim 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 the LIOMs (possibly an increase of the support of $\tau^{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 $\tau^{z}_{i}$ are significantly affected by the long-range interactions.

In summary, we have shown that the cavity mediated long-range interactions lead to non-ergodic behavior of the system, which in presence of strong disorder may be interpreted as MBL. However, the logarithmic growth of entanglement entropy, a hallmark of standard MBL is missing in our model, at least for very small $U_{1}$. The observed ergodicity breaking exhibits novel features, such as the rapid growth of entanglement entropy which we attribute to the long-range interactions. When the effect of the long-range interactions can be separated from the short-range interactions in the dynamics, 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. On the other hand, the dynamics and in particular the growth of the configuration entropy may be explained using the language of LIOMs as described above, the system is integrable (as the gap ratio analysis reveals poissonian statistics) and the entropy growth saturates to smaller values than those found for ergodic dynamics. For these reasons we consider the behavior observed as a nonstandard variant of MBL. Let us also note that the observed phenomenon is different from the non-ergodic phase observed for systems with single particle mobility edge [70, 71, 72, 73] where long-range interactions are not involved.

5 Detection of ergodicity breaking: Light at the cavity-output

The dynamics of the system can be experimentally revealed via a site-resolved measurement of $n_{i}$, which can be performed in cold atoms experiments [7] and which allows one to reconstruct the correlation function $C(t)$. In this section we argue that the cavity setup can allow one to measure the breakdown of ergodicity by photo-detection of the light at the cavity output. The emitted light, in fact, is scattered by the atoms and thus contains the information on their density distribution. In particular, the electric field amplitude, which we denote by $E_{\rm out}$, is directly proportional to the expectation value of population imbalance $I(t)=\sum_{i}(-1)^{i}n_{i}$, Eq. (6) in the dispersive cavity regime (see Appendix A and [74, 31, 32, 37, 75]):

[TABLE]

which can be measured via heterodyne detection [74] (Recall that an observable similar to $I(t)$ was employed in Ref. [7] to demonstrate MBL for a system of fermions in the optical disordered lattice). This measurement introduces projection noise that affects the atomic dynamics. Nevertheless, as we detail in Appendix A, this noise is negligible in the limit we consider, where the cavity field dynamics occurs on a much faster time scale than the atomic motion. The information extracted from the light at the cavity output has been used in this fashion to extract information about ground (or metastable) state phases of similar systems [31, 37, 75].

Below we discuss the time evolution of the imbalance $\langle I\rangle$ and of its square $\langle I^{2}\rangle$ (corresponding to the light intensity) for two states in which the system can be initially prepared. We first consider a density-wave like state $|DW_{10}\rangle=|101010..\rangle$, with odd sites occupied and even sites empty. This state maximizes $\langle I\rangle$ for the fractional density $\bar{n}=1/2$ and minimizes the energies of both short-range and long-range interactions. In the absence of disorder $|DW_{10}\rangle$ is an eigenstate of the Hamiltonian in the atomic limit [34]. Starting from this state, in the following we investigate to what extent the experimentally accessible time evolution of the population imbalance $\langle I\rangle$ allows one to probe ergodicity and its breaking in the system. In this analysis one shall keep in mind that for $U_{1}\gtrsim W$ the state $|DW_{10}\rangle$ has a significant overlap with the ground state, and its energy thus gets closer to the ground state energy as the ratio $U_{1}/W$ is increased. Figure 9 diplays the time evolution of $\langle I\rangle$ for the system initially prepared in $|DW_{10}\rangle$ state. The dynamics is shown for $W=0.5$ (upper panels) and $W=10$ (lower panels). These two disorder strengths $W$ correspond to the ergodic and the MBL regimes, respectively, for the considered values of $U_{1}$. For comparison, we also display the dynamics of the correlation function $C$ 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$ 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}\rangle$ depends strongly on the ratio $U_{1}/W$. For $W=0.5$ and $U_{1}=0.16$, the imbalance $\langle I\rangle$ decays to zero and is a valid probe of the ergodic properties of the system. For the same disorder amplitude and larger values of $U_{1}$, instead, $\langle I\rangle$ saturates at a constant value, from which we infer that the state $|DW_{10}\rangle$ has already significant overlap with the ground state of the system. In this regime, the time evolution of $\langle I\rangle$ is not a valid probe of ergodicity of a typical state of the system.

For strong disorder, $W=10$, on the other hand, the information about localization properties of the system provided by the correlation function $C$ is also visible in the time evolution of the imbalance $\langle I(t)\rangle$, even for $U_{1}=3.2$. This behaviour is due to the fact that the energy spectrum of the system is much broader, such that even for values of the interaction strength as large as $U_{1}=3.2$ the state $|DW_{10}\rangle$ is still in a region of the spectrum with a relatively large density of states.

We now consider the initial state $|DW_{2}\rangle=|11001100...\rangle$, where the initial site occupations form a pattern of singly-occupied and empty doublets of lattice sites. State $|DW_{2}\rangle$ is an eigenstate of the imbalance operator $I$ at the eigenvalue [math], therefore its energy lies about the band center. Figure 10 displays the time evolution of squared population imbalance $\langle I^{2}\rangle$ when the initial state is $|DW_{2}\rangle=|11001100...\rangle$ and for different values of the disorder amplitude: $W=1,3,7,10$. The upper panels report, for comparison, the correlation function $C$ for the corresponding values of the disorder amplitude but when the initial state is a random Fock state. The dynamics of the density-density correlation function $C$ shows that the system is ergodic for $W=1,3$ with correlations eventually decaying to zero. For $W=7,10$, instead, $C$ saturates to non-zero values at large times signalling MBL. The behaviours for vanishing long-range interactions and for $U_{1}=2$ are similar, see Fig. 10 a) and b). If the system is ergodic, the “intensity” $\langle I^{2}\rangle$ increases from zero over a time scale of few $J^{-1}$ till it reaches a nonvanishing saturation value which depends on the disorder strength $W$. As visible in Fig. 10 c), the saturation value is maximal in the ergodic phase at $W=1$ and is smaller in the MBL phase at $W=7,10$. The saturation value of $\langle I^{2}\rangle$ thus reflects the ergodicity of the system dynamics. The saturation values of $\langle I^{2}\rangle$ still point towards non-ergodicity for $U_{1}=2$ at $W=7,10$, consistently with the value of mean gap ratio $\overline{r}$ for system size $K=16$, see Fig. 2. Note that the saturation values of $\langle I^{2}\rangle$ at $U_{1}=2$ are larger than the ones reached for $U_{1}=0$ at the same disorder strength. This behaviour is consistent with the intuition that long-range interactions tend to delocalize. In Fig. 10(d) one observes that the saturation value of $\langle I^{2}\rangle$ is lower for $W=1$ than for $W=3$. This suggests that the subspace of the Hilbert space accessible during the time evolution from $|DW_{2}\rangle$ state is constrained by long-range interactions at small disorder strength $W=1$, and that this constraint becomes weaker at $W=3$.

In summary, we have argued that the measurement of light emitted by the cavity can allow one to determine the localization properties 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. [76, 77] used in standard population imbalance measurements.

6 Quasi-random disorder

Up till now we have considered a random on-site disorder. Such a situation may be realized by an off-resonant pumping laser with a random intensity distribution. This laser shall drive an atomic transition that does not scatter into the cavity field, and thus generate a random a.c. Stark shift of the onsite energy. The disorder may. also be realized with the setups of Ref. [37, 78] by introducing additional weak laser beams creating optical lattices, whose periodicity is incommensurate with the cavity lattice periodicity. This setup creates quasi-random disorder analogously to the experiment of Ref. [7]. We now analyse localization in such a case i.e. when the onsite energy is due to the contribution of an incommensurate periodic potential, namely,

[TABLE]

Here, $W$ plays a role of disorder amplitude, $\beta$ is the ratio of the two lattice periodicities (we take $\beta=(\sqrt{(}5)-1)/2$) and the value of the phase $\phi\in[0,2\pi]$ determines the disorder realization. Quasiperiodic potentials such as the one of Eq. (16) have been employed in a number of experiments investigating MBL in optical lattices [7, 79, 80, 73].

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 [81] and of the gap ratio [47]. 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 similar features are found in our extended Hubbard model, we here discuss the dynamics of spinless fermions when the onsite energy in Eq. (2) is given by Eq. (16).

The time evolution of the correlation function $C(t)$ and of the entanglement entropy $S(t)$ are displayed 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 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 of spinless fermions with quasiperiodic disorder and in the presence of cavity-mediated all-to-all interactions are analogous to the ones found when the disorder is instead random.

7 Bosons in the cavity

Many-body localization in bosonic systems was studied numerically in [9, 82, 10] as well as in experimental realizations [64, 83]. 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 [37, 78] are performed with bosons, here we briefly discuss signatures of the MBL in the extended Bose-Hubbard Hamiltonian:

[TABLE]

Figure  13 displays the phase diagram for a lattice consisting of $K=8$ sites at unit filling. The diagram reports the mean gap ratio $\bar{r}$, obtained from the gap ratio $r$ averaged over the states at the center of the spectrum $\epsilon\approx 0.5$. It 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 similar to the phase diagram of spinless fermions in Fig. 2. We note however, that the values of disorder required to induce many-body localization in the bosonic system of Eq. (17) are larger than the ones of the fermionic counterpart.

To provide further insight into the physics of the MBL phase in all-connected bosonic system we calculate the bipartite entanglement entropy during the course of time evolution of the system – c.f. Fig. 14. The features are overall similar to the ones found for spinless fermions and thus seem to not be limited to small Hilbert spaces per site as for fermions.

8 Conclusions

In this work we have analysed the occurrence of many-body localization in a system of particles with all-to-all interactions. The dynamics is described by an extended Hubbard model, 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 that the transition to the ergodic phase occurs at disorder amplitudes whose critical value increases with the long-range interaction strength. The all-to-all interactions affect the spreading of the entanglement during the time evolution: the growth of entanglement entropy is faster than logarithmic in the MBL phase. It is instead linear when the tunneling coefficient is much larger than long-range interaction strength. Nevertheless, the saturation value of entanglement entropy is constrained by particle number fluctuations between partitions of the system. A fixed distribution of particle number between subsystems is attained during the course of time evolution 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 the 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. Similar 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, photo-detection of the light at the cavity output can allow one to probe the gas phase and thus provide evidence of ergodicity breaking.

We finally remark that MBL phase in presence of global coupling to $d$ dimensional system was found in [24], particular scaling of $d$ with system size leads in high frequency limit to Hamiltonian similar to the one considered in this work.

Acknowledgements

This work was performed with the support of EU via Horizon2020 FET project QUIC (nr. 641122). Numerical results were obtained with the help of PL-Grid Infrastructure. We acknowledge support of the National Science Centre (PL) via project No.2015/19/B/ST2/01028 (K.B. and P.S.), 2018/28/T/ST2/00401 (Etiuda scholarship – P.S.) and the QuantERA QTFLAG programme No. 2017/25/Z/ST2/03029 (J.Z.). G.M. acknowledges support by the German Research Foundation (DFG, Priority Program No. 1929 GiRyd) and by the German Ministry of Education and Research (BMBF) via the QuantERA project NAQUAS is acknowledged. Projects NAQUAS and QTFLAG have received funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme.

Appendix A Derivation of the Hubbard Hamiltonian

Here we present the summary of the derivation of Hamiltonian ((1)-(5)). The details in different variants of the model are discussed in e.g. Refs. [26, 28, 31, 33, 34, 35].

We consider $N$ atoms of mass $m$, they are confined within an optical cavity in a quasi-one-dimensional configuration colinear with a one-dimensional optical lattice. The lattice oriented along the cavity axis has the same wave number $k=2\pi/\lambda$ as the cavity field. The atoms are prepared in the electronic ground state $|1\rangle$, the electric dipole transition $|1\rangle\to|2\rangle$ couples dispersively to the cavity and to the transverse laser with Rabi frequency $\Omega_{z}(x)$, see Fig. 1. Atoms and cavity field are treated in second quantization: $\hat{a}$ and $\hat{a}^{\dagger}$ are the annihilation and creation operators of a cavity photon, respectively, and obey the commutation relation $[\hat{a},\hat{a}^{\dagger}]=1$; the atomic field operator $\hat{\Psi}_{j}(x,t)$ destroys an atom in the internal state $\left|j=1,2\right>$ at position $x$ and time $t$, the commutation relations are $[\hat{\Psi}_{i}(x,t),\hat{\Psi}^{\dagger}_{j}(x^{\prime},t)]_{\pm}=\delta_{ij}\,\delta(x-x^{\prime})$ (where $\pm$ indicate the anti- and commutation relations, depending on whether the atoms are fermions or bosons, respectively). Cavity and atomic operators commute at equal times. The Hamiltonian governing the dynamics can be decomposed into the sum of atoms, electromagnetic field, and atom-field interactions. The atomic Hamiltonian reads

[TABLE]

where

[TABLE]

Here, $\Delta_{a}=\omega_{L}-\omega_{0}$ is the detuning between the laser frequency $\omega_{L}$ and the atomic transition frequency in the frame rotating with the laser frequency $\omega_{L}$, $V_{\rm cl}^{(j)}$ is the state-dependent depth of the optical lattice tightly binding atoms to its minima. The optical lattice potential is due to light shifts by a laser driving a different atomic transition. The other parameters are the collision rate $U_{jj}$, which is also state dependent, and the collision rate $U_{12}$ between atoms in different electronic states. Since the laser is here described by a classical field, the only quantum field is the cavity mode at frequency $\omega_{C}=ck$, with $c$ the speed of light. The field Hamiltonian in the frame rotating with the laser frequency takes the form

[TABLE]

where $\Delta_{C}=\omega_{L}-\omega_{C}$. Finally, the Hamiltonian describing the interaction between the atomic dipoles and the electric fields is given by

[TABLE]

where $g(x)=g_{0}\cos(kx)$ is the cavity vacuum Rabi frequency and $\Omega_{z}(x)$ is the position-dependent (real-valued) Rabi frequency of a laser propagating along the $z$ direction and coupling to the atomic transition (compare Fig. 1).

We now assume that the photon scattering processes are elastic, which is valid assuming that the largest frequency is the detuning $|\Delta_{a}|$. In particular, $|\Delta_{a}|\gg\gamma$, where $\gamma$ the radiative linewidth of the excited state, and $|\Delta_{a}|\gg\Omega,g_{0}\sqrt{n_{\rm cav}},|\Delta_{C}|$, namely, the detuning is much larger than the strength of the coupling between the ground and excited state, where $n_{\rm cav}=\langle\hat{a}^{\dagger}\hat{a}\rangle$ is the mean intracavity photon number. In this regime we can eliminate the excited state approximating [27, 28]

[TABLE]

which is valid to lowest order in the expansion in $1/|\Delta_{a}|$. Using Eq. (22) in the Heisenberg equation of motion for the field operator $\hat{\Psi}_{1}(x,t)$ results in the equation (we now drop the subscript from the field operator: $\hat{\Psi}_{1}\equiv\hat{\Psi}$)

[TABLE]

where we have taken care to keep the ordering between cavity and atomic operators. This equation is now solely coupled to the Heisenberg-Langevin equations for the cavity field:

[TABLE]

which depends on the atomic operators through the atomic density

[TABLE]

Equation (24) includes the quantum noise due to the cavity losses, with $\kappa$ the cavity loss rate and $\hat{a}_{\rm in}(t)$ is the input noise operator, with $\langle\hat{a}_{\rm in}(t)\rangle=0$ and $\langle\hat{a}_{\rm in}(t)\hat{a}_{\rm in}^{\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$ [84]. The other parameters are the frequency $U_{0}(x)=g(x)^{2}/\Delta_{a}$, which scales the depth of the intracavity potential generated by a single photon, and the scattering amplitude $S(x)=g(x)\Omega_{z}(x)/\Delta_{a}$ [26, 27].

We now eliminate the cavity degrees of freedom from the atomic dynamics by identifying the time-scale $\Delta t$ over which the atomic motion does not significantly evolve while the cavity field has relaxed to a state which depends on the atomic density at the given interval of time. This is verified when $\Delta t\gg T_{c}$, where $T_{c}=1/|\Delta_{C}+{\rm i}\kappa|$, and $\sqrt{\omega_{R}E_{\rm kin}}\ll\hbar/\Delta t$, with $\omega_{R}=\hbar k^{2}/(2m)$ the recoil frequency and $E_{\rm kin}$ the mean kinetic energy [36]. We also require that the coupling strengths between atoms and fields, which determine the time scale of the evolution due to the mechanical effects of the interaction with the light, are much smaller than $1/\Delta t$. In this limit, we replace the field operator with its coarse-grained average $\hat{a}_{\rm st}(\bar{t})$, defined as

[TABLE]

such that $\int_{\bar{t}}^{\bar{t}+\Delta t}\dot{\hat{a}}_{\rm st}(\tau)d\tau=0$, with $\dot{\hat{a}}$ given in Eq. (24). The ”stationary” cavity field is a function of the atomic operators at the same (coarse-grained) time, and in particular it takes the form

[TABLE]

with $\hat{\bar{a}}_{\rm in}$ the input noise averaged over $\Delta t$, such that in the coarse-grained time scale $\langle\hat{\bar{a}}_{\rm in}(\bar{t})\rangle=0$ and $\langle\hat{\bar{a}}_{\rm in}(\bar{t})\hat{\bar{a}}_{\rm in}^{\dagger}(\bar{t}^{\prime})\rangle=\delta(\bar{t}-\bar{t}^{\prime})$. Note that the commutation relations of the new operators are modified, and in particular the commutator between $\hat{a}_{\rm st}$ and $\hat{a}_{\rm st}^{\dagger}$ scales as $T_{c}/\Delta t$, The quantum noise term is an effective term which provides the same averages as the corresponding quantum-noise operator one gets by formally integrating Eq. (24). It can be neglected when the mean intracavity photon number is larger than its fluctuations or by integrating for sufficiently long time. In this limit the dynamics is coherent and the field at the cavity output,

[TABLE]

allows one to monitoring the state of the atoms [84, 25, 85, 33].

Using Eq. (26) in place of the field $\hat{a}$ in Eq. (23), and discarding the noise term, leads to an equation of motion for the atomic field operator which depends solely on the atomic variables[27, 28]:

[TABLE]

where we have assumed $|\Delta_{C}|\gg|U_{0}|n,\kappa$ and thus discarded the terms scaling with $U_{0}$ in the denominators of Eq. (26) (see [28] for a discussion). The term of the right-hand side of Eq. (28) can be cast in terms of the commutator between $\hat{\Psi}$ and the effective Hamiltonian $H=H_{A}+H_{\rm CQED}$, where [28]:

[TABLE]

which contains an infinitely ranged density-density interaction. This interaction is attractive for $\Delta_{C}<0$, which is the regime we assume in this work.

The Hubbard model is obtained when the external optical lattice is sufficiently deep to tightly bind the atoms in the lowest band and when the cavity interactions are a sufficiently small perturbation. We denote by $K$ the number of lattice sites and perform the Wannier decomposition of the atomic field operator $\hat{\Psi}(x)=\sum_{i}w_{i}(x)\hat{b}_{i}$ [86, 87]. Here, $w_{i}(x)$ denotes the Wannier function of the classical optical lattice that is centered at lattice sites $x_{i}=ia$, with $a=\lambda_{0}/2$ the lattice periodicity, while $\hat{b}_{i}$ annihilates a particle at the corresponding lattice site. We further assume the scaling $g(x)=\tilde{g}(x)/\sqrt{K}$ (thus $S(x)=\tilde{S}(x)/\sqrt{K}$), which is equivalent to assuming that the cavity mode volume scales linearly with the size of the lattice [88, 28]. This scaling is equivalent to Kac’s scaling and warrants that the energy is extensive despite the all-to-all interactions (even though it remains non-additive). The Hubbard term due to the cavity global interactions takes the form (recall $\Delta_{C}<0$) [33, 34]:

[TABLE]

with $\hat{n}_{i}=\hat{b}_{i}^{\dagger}\hat{b}_{i}$ the occupation of site $i$. For $\tilde{S}(x)=S_{0}\cos(kx)$, corresponding to a homogeneous transverse pump $\Omega_{z}(x)=\Omega$, then $S_{0}=\tilde{g}_{0}\Omega/\Delta_{a}$ and

[TABLE]

Thus, one obtains Eq. (5) with the coefficient

[TABLE]

while the remaining terms in the tight binding approximation yield the Hubbard model. Let us note that the (coarse-grained) field $\hat{\bar{a}}_{\rm out}(t)$ may be expressed as (neglecting the noise term in the coarse-grained limit)

[TABLE]

Hence the mean electric field at the cavity output is proportional to the expectation value of the population imbalance,

[TABLE]

where we used the definition in Eq. (6). We remark that, for the same reason that the time-scale separation allows one to derive a Hamiltonian dynamics for the atomic motion, thus discarding the commutation relation between the cavity field operators, then projection noise due to cavity leakage and photo-detection can be neglected because it is of the same order as the non-adiabatic corrections.

Appendix B Nonergodic regime for strong all-to-all interactions

Even in the absence of disorder strong long-range interactions may lead to nonergodic behavior manifesting itself e.g. in the logarithmic or sublinear growth of the entanglement entropy for a sudden quench [89, 90, 55]. This behavior seems to be distinct from MBL. An appearance of such a regime in the model studied here may be visualized considering mean gap ratio for large $U_{1}$ values as shown in Fig. 15. For large $U_{1}$ system tends towards gap ratio $\overline{r}$ corresponding to Poisson statistics.

Similar conclusions may be obtained from the time dynamics of the density-density correlation function $C(t)$ and of bipartite entanglement entropy $S(t)$ starting from initial random Fock state as shown in Fig. 16. Large $U_{1}$ leads to slower decay of correlations accompanied by a slower sublogarithmic even growth of the entanglement entropy. While for $K>16$ we cannot follow the dynamics for sufficiently long times the results indicate a clear saturation of $C(t)$ for $U_{1}=25$ in agreement with mean gap ratio values.

Appendix C 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 [91, 92]

[TABLE]

where $a=(E_{\rm max}-E_{\rm min})/2$, $b=(E_{\rm max}+E_{\rm min})/2$, the Hamiltonian is rescaled $\mathcal{H}=\frac{1}{a}(H-b)$ so that spectrum of $\mathcal{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 (35) to time $t_{\rm max}$ is $N\approx 2at_{\rm max}$ [93].

The time-evolution of the initial state $|\psi_{0}\rangle$ is given by

[TABLE]

and reduces to matrix-vector multiplications

[TABLE]

where the recursion relation satisfied by Chebyshev polynomials was used. In order to get $|\psi(t)\rangle$ we generate iteratively a sequence of $N$ vectors $|\psi_{0}\rangle,T_{1}|\psi_{0}\rangle,...,\,T_{N}|\psi_{0}\rangle$. To reach long times of time evolution $t_{\rm max}\approx 10^{3}$ one needs relatively large $N$ which increases memory consumption. Therefore we split the time interval $[0,t_{\rm max}]$ into parts $[0,\Delta t],[\Delta t,2\Delta t],...$ in such a way that $|\psi\left((n+1)\Delta t\right)\rangle$ can be calculated from the state $|\psi\left(n\Delta t\right)\rangle$ with the expansion (36) involving only a limited number of terms e.g. – $N\approx 1000$ which allows us to obtain time evolution for the system size $K=20$ with memory consumption smaller than $5$GB (performing the matrix-vector multiplications in PETSc).