Self-organized topological insulator due to cavity-mediated correlated tunneling

Topological materials have potential applications for quantum technologies. Non-interacting topological materials, such as e.g., topological insulators and superconductors, are classified by means of fundamental symmetry classes. It is instead only partially understood how interactions affect topological properties. Here, we discuss a model where topology emerges from the quantum interference between single-particle dynamics and global interactions. The system is composed by soft-core bosons that interact via global correlated hopping in a one-dimensional lattice. The onset of quantum interference leads to spontaneous breaking of the lattice translational symmetry, the corresponding phase resembles nontrivial states of the celebrated Su-Schriefer-Heeger model. Like the fermionic Peierls instability, the emerging quantum phase is a topological insulator and is found at half fillings. Originating from quantum interference, this topological phase is found in"exact"density-matrix renormalization group calculations and is entirely absent in the mean-field approach. We argue that these dynamics can be realized in existing experimental platforms, such as cavity quantum electrodynamics setups, where the topological features can be revealed in the light emitted by the resonator.


Introduction
Manifestation of topology in physics [1,2] created a revolution which is continuing for almost four decades. With the discovery of topological materials, condensed matter physics has gained a new terrain where quantum phases of matter are no longer controlled by local order parameters as in paradigmatic Landau theory of phase transitions but rather by the conservation of certain symmetries [3]. These new phases of matter, so called symmetry-protected topological (SPT) phases, display edge and surface states that can be robust against perturbations, making them genuine Titas Chanda: titas.chanda@uj.edu.pl candidates for quantum technologies [4].
To date, a detailed understanding of noninteracting topological materials, such as e.g., topological insulators and superconductors [5], has been obtained through a successful classification based on fundamental symmetry classes, the so-called "ten-fold way" [6][7][8]. On the other hand, interactions are almost unavoidable in many-body systems. It is therefore a central issue to understand whether SPT phases can survive the inter-particle interactions, or perhaps whether interactions themselves might stabilize SPT phases and even give rise to novel topological properties [9][10][11][12][13][14][15][16]. These questions are at the center of an active and growing area of research [17,18]. Recent works have argued that the range of interactions plays a crucial role on the existence of SPT phases. Specifically, in frustrated antiferromagnetic spin-1 chain with power-law decaying 1/r α interactions, topological phases are expected to survive at any positive value of α [19]. This behavior shares similarities with the topological properties of noninteracting Kitaev p-wave superconductors [20] that are robust against long-range tunneling and pairing couplings [21][22][23][24]. It was found that for infinite-range interactions the one-dimensional Kitaev chain can exhibit edge modes for appropriate boundary conditions [25]. On the other hand, in spin chains and Hubbard models, the infinite-range interactions suppress the onset of hidden order at the basis of the Haldane topological insulator [26,27].
In this work, we report a novel mechanism where nontrivial topology emerges from the quantum interference between infinite-range interactions and singleparticle dynamics. We consider bosons in a onedimensional (1D) lattice, where the global interactions have the form of correlated tunneling resulting from the coupling of the bosons with a harmonic oscillator. We identify the conditions under which this coupling spontaneously breaks the discrete lattice translational symmetry and leads to the emergence of non-trivial edge states at half filling. The resulting dynamics resembles the one described by the famous Su-Schrieffer-Heeger (SSH) model [28,29]. The topology we report shares analogies with the recent studies of symmetry breaking topological insulators κ Ω Figure 1: Interference-induced topological phases can be realized in a cavity quantum electrodynamics setup. The bosons are atoms (red spheres) confined by a one-dimensional optical lattice and dispersively interacting with a standingwave mode of the cavity. The cavity standing-wave field is parallel to the lattice, its wavelength is twice the lattice periodicity, and the atoms are trapped at the nodes of the cavity mode. Correlated hopping originates from coherent scattering of laser light into the cavity, the laser Rabi frequency Ω controls the strength of the interactions. The field at the cavity output is emitted with rate κ and provides information about the phase of the bosons. [14-16, 30, 31]. Nevertheless, differing from previous realizations, here the interference between quantum fluctuations and global interactions is essential for the onset of the topological phase and cannot be understood in terms of a mean-field model.
We argue that these dynamics can be realized, for instance, in many-body cavity quantum electrodynamics (CQED) setups [32][33][34][35][36][37], like the one illustrated in Fig. 1 highlighting the experimental feasibility of our proposal. In the following, we first give a brief description of the model in relation to the experimental setup. Then, we ascertain the phase diagram of the system at half filling, characterize the topological phase, and point out how to experimentally verify its existence.

Bose-Hubbard model with global correlated tunneling
The model we consider describes the motion of bosons in a 1D lattice with N s sites and open boundaries. The specific experimental situation realizing this dynamics is sketched in Fig. 1. The dynamics is governed by the HamiltonianĤ in terms of the field operatorsb j andb † j , which annihilate and create, respectively, a boson at site j in the lowest lattice band. The bosons, tightly confined in the lowest band of the lattice, strongly couple with a cavity standing-wave mode (the harmonic oscillator), which is parallel to the lattice. When the atoms are transversely driven by a laser, the Hamiltonian takes the detailed form of a Bose-Hubbard (BH) model with the additional optomechanical coupling with the oscillator [38]: where we have set = 1. Here, the first two terms on the right-hand side (RHS) of Eq. (1) are the nearestneighbor hopping with amplitude t and onsite repulsion with magnitude U . The third term on the RHS is the bosons-cavity coupling, where operatorsâ and a † annihilate and create, respectively, a cavity photon, while the operatorB acts on the bosonic Hilbert space. The last term of (1) gives the cavity oscillator energy in the reference frame rotating at the frequency of the transverse laser with ∆ c the detuning between cavity and laser frequency. Finally, the bosons-cavity coupling is scaled by the coefficients S and y that can be independently tuned. The first one, in fact, is proportional to the amplitude Ω of the transverse laser field, c.f. Fig. 1. The coefficient y, instead, depends on the localization of the scattering atoms at the lattice sites (see Appendix A). The specific form of operatorB depends on the spatial dependence of the cavity-bosons coupling [39][40][41][42]. In the present case, we consider The staggered coupling, (−1) i+1 b † ib i+1 + H.c , originates from the spatial mode function of the cavity mode along the lattice when the lattice sites are localized at the nodes of the cavity mode, and when the periodicity of cavity mode standing wave is twice the periodicity of the optical lattice.
Hamiltonian (1) with (2) is reminiscent of the phonon-electron coupling of the SSH model [43]. Differing from the ionic lattice of the SSH model, the bosons couple to a single oscillator -the cavity mode. Here, the instability is associated with a finite stationary value of the field quadraturex =â +â † : for x = 0 the bonds connecting even-odd and odd-even sites differ by a quantity proportional to x . Sincex is a dynamical variable, this process is accompanied by a spontaneous breaking of the lattice translational symmetry. In a cavity, for instance,x is the electric field that is scattered by the atoms and depends on the atomic mobility. The resulting bosonic dynamics is thus intrinsically nonlinear.
We analyze the quantum phases of the system in the limit in which the cavity field (oscillator) can be eliminated from the equations of the bosonic variables assuming that |∆ c | is the largest frequency scale of the dynamics. In this limit, the time-averaged field iŝ ε(τ ) = 1 ∆t τ +∆t τ dtâ(t) ≈ −SB(τ )/∆ c , where τ is the coarse-grained time in the grid of step ∆t [39,44]. The effective Bose-Hubbard Hamiltonian takes the form where U 1 = S 2 N s /∆ c and the explicit dependence on N s warrants the extensivity of the Hamiltonian in the thermodynamic limit when U 1 is fixed by scaling S ∝ 1/ √ N s . The term proportional toB 2 emerges from the back-action of the system through the global coupling with the oscillator. It describes a global correlation between pair tunnelings.
Before discussing the emerging phases, we note that the model in Eq. (1) has extensively been employed for describing ultracold atoms in optical lattices and optomechanically coupled to a cavity mode [32,33,39,41,45]. The specific form of the coupling with operator (2) has been discussed in Ref. [41] and can be realized by suitably choosing the sign of the detuning between laser and atomic transition as in Ref. [34]. The scaling 1/N s of the nonlinear term corresponds to scaling the quantization volume with the lattice size [46]. In the time-scale separation ansatz, the bosons dynamics is coherent provided that |∆ c | is also larger than the cavity decay rate κ: in this regime shot-noise fluctuations are averaged out. Correspondingly, there is no measurement back-action, since the emitted fieldε(τ ) is the statistical average over the time grid ∆t [47]. We note that topological phases can also be realized in plethora of platforms where the range of interactions can be tuned and the geometry controlled, prominent examples are optomechanical arrays [48], photonic systems [49], trapped ions [50], and ultracold atoms in optical lattices [32,[51][52][53][54][55][56]. Specifically, this model could be also realized in the experimental setups as discussed in [57][58][59]. In this work, we shall keep on referring to a CQED setup, where the dynamics predicted by Hamiltonian (1) has been extensively studied and can be realized [32][33][34].

Phase diagram at half filling
In the rest of this work, we consider the system at half filling, i.e., density ρ = 1/2. We determine the ground state and its properties using the density matrix renormalization group (DMRG) method [60,61] based on matrix product states (MPS) ansatz [62,63] -for details see Appendix B. The phase diagram is determined as a function of t/U and U 1 /U . The phases and the corresponding observables are summarized in Table 1, the observables are detailed below. We characterize off-diagonal long-range order by the Fourier transform of the single particle correlations, i.e., the single particle structure factor which can be experimentally revealed by means of time-of-flight measurements [64,65]. The coupling with the cavity induces off-diagonal long-range order, that is signaled by the bond-wave order parameter This quantity is essentially the cavity field in the coarse graining time scale and is directly measurable by heterodyne detection of the electric field emitted by the cavity [66]. We also consider the density-wave order parameter, O D = 1 Ns j (−1) jn j , which signals the onset of density-wave order and typically characterizes phases when the lattice sites are at the antinodes of the cavity field (see Appendix A). Moreover, we analyze the behavior of the string and parity order parameters: These order parameters depend nonlocally on onsite fluctuations δn j =n j − ρ from the mean density ρ. In our calculations of O S and O p , we take i = 10 and j = N s − 11 to neglect the boundary effects. We first remark that, for U 1 = 0, hence in the absence of global interactions, the phase is superfluid (SF). We also expect that for U 1 > 0 the quantum phase of Hamiltonian in Eq. (3) is SF, since the formation of a finite cavity field costs energy. Instead, we expect that correlated hopping becomes relevant for U 1 < 0. We have first performed a standard Gutzwiller mean-field analysis of the model assuming two-site translational invariance. At sufficient large values of |U 1 | mean-field predicts the formation of a SF phase of the even or odd bonds accompanied by a finite value of bond-wave order parameter O B . This phase maximizes the cavity field amplitude and we denote it by Bond Superfluid phase (BSF). We point out that mean-field predicts that the ground state exhibits off-diagonal long-range order for any value of U 1 .
We now discuss the quantum phase obtained from DMRG calculations. The phase diagram is reported as a function of the ratio t/U , which scales the strength of the single-particle hopping in units of the onsite repulsion, and of the ratio U 1 /U , which scales the strength of the correlated hopping. We sweep U 1 from positive to negative values. We note that in a cavity the sign of U 1 is tuned by means of the sign of the detuning ∆ c . The effective strength, in particular, shall be here scaled by the parameter y 2 , depending on the particle localization. Here, y is constant across the diagram, since we keep in fact the optical lattice depth constant and tune the ratio t/U by changing the onsite repulsion U (experimentally, this is realized by means of a Feshbach resonance). where OB has been determined using DMRG. The mean-field Gutzwiller approach for bond-wave order parameter agrees well with DMRG in SF and BSF regimes, but misses the appearance of BI phase, which is further verified by computing the mean-field results for the maximum of single particle structure factor max M M F the SF from the BSF phase, where the effective tunneling amplitudes b † ib i+1 + H.c. attain a staggered pattern characterized by a finite value of bond-wave order parameter O B . The long range coherence of the BSF phase is manifested by narrow peaks of M 1 (k) centered at k = ±π/2 (see Fig. 2). Remarkably, we observe a reentrant insulating phase separating the SF and the BSF. The insulator is signaled by vanishing off-diagonal long-range order and therefore by vanishing structure factor M 1 (k). It is characterized by the non-zero (zero) values of the string order parameter and by vanishing (non-vanishing) parity order parameter depending on the boundary sites of these non-local parameters (see the next section for details). We denote this phase as a Bond Insulator (BI). This phase is separated from the SF by a continuous phase transition. The transition BI-BSF is also continuous.
We note that the bond insulator phase is entirely absent in the standard Gutzwiller mean-field approach [67][68][69][70][71], where the bosonic operators are decomposed asb with Φ j = b j being the superfluid order parameter in the Gutzwiller analysis. By performing such meanfield analysis with two-site unit cell, we can obtain the bond-wave order parameter O M F B at the mean-field level. Figure 3 displays the mean-field bond-wave order O M F B (top left panel) and the deviation of it from the DMRG result (top right panel). While the exact borders between various phases quantitatively differ, the mean field Gutzwiller approach agrees well with DMRG in SF and BSF regimes, but misses the appearance of BI phase. To be sure about this, we check the mean-field value of the maximum of single particle structure factor M 1 (k), which is to be redefined as The max M M F 1 (k) is presented in the bottom panel of Fig. 3. It has possesses high values in the entire (t/U, U 1 /U )-plane confirming that the mean-field analysis cannot capture the insulating BI phase where M 1 (k) must be vanishingly small. Moreover, the mean-field value of the string-order parameter O S remains exactly zero in the entire phase diagram. It is to be noted that max M M F 1 (k) does not match the DMRG results presented in Fig. 2, as in reality the offdiagonal correlations b † ib j decay algebraically with the distance |i − j| in the superfluid phases, while at the mean-field level Φ i Φ * j does not. Indeed the novel BI phase is entirely due to the interplay between the long-range coupling induced by the cavity and the single-particle tunneling. Due to this quantum interference, the insulating BI phase reveals itself in quasi-exact calculations like DMRG, while Gutzwiller analysis can only capture two superfluid phases. We further note that studies of the ground state of (3), based on exact diagonalization for small system sizes, did not report the existence of the BI phase [41]. In the following section, we characterize the topology associated with the BI phase and argue that it is stable in the thermodynamic limit.

Emergent topology associated with the BI phase
By means of excited-state DMRG, we reveal that the BI phase has triply degenerate ground state (quasidegenerate for finite N s ) separated by a finite gap from the other excited states. The site distribution is visualized in Fig. 4 which shows that the absolute ground state has a uniform mean half-filling, while the other two states possess edge excitations, namely, fractional 1/2 particle-hole excitations with respect to the mean half-filling (bottom two rows of Fig. 4). Such edge excitations are characterized by the bondwave order parameter with opposite sign than the trivial phase. They suggest that the BI phase is a symmetry protected topological (SPT) phase. Similar topological edge states have been reported e.g., for noninteracting system [56] or in superlattice BH model [72], where the superlattice induces a tunneling structure resembling that of the SSH model [28,29]. In our case, instead, the effective tunnelings are spontaneously generated by the creation of a cavity field that breaks discrete Z 2 translational symmetry of the system. However, bond centered inversion symmetry still remains intact -it protects this SPT phase. This ib i+1 + H.c. as a function of the bonds (i, i + 1). Orange (teal) bars denote the even (odd) bond. Right panels: Density ni as a function of the lattice site. The dashed line is a guide for the eye. Observe the characteristic alternate weak/strong pattern in the bonds with weak bonds occurring at the edges for topological states that reveal topological particle-hole edge excitations. These excitations are due to exactly fractional ±1/2 particle localizations on the edges. Here, we set U1/U = −10, t/U = 0.05, and y = −0.0658 (see Appendix A) to obtain the states using DMRG algorithm.
On a further inspection, it is found that the string order O S and parity order O P can be non-zero (zero) depending on the location of the two separated sites that sit at the boundaries of the non-local operators (sites i and j in Eqs. (6)). We illustrate it in Fig. 5(a). We find that O S = 0 and O P = 0, as reported in Fig. 2, when the non-local operators start at the second site of a strong bond (i.e., the bond with larger tunneling element) and end at the first site of a weak bond (the bond with smaller tunneling element) further away. That is why we have chosen i = 10 (even site) and j = N s −11 (odd site) to calculate these nonlocal operators for the trivial state in Fig. 2. These are unusual properties when compared to, say, topological Haldane phases of extended BH models at unit filling [27,74].
To check whether we are indeed dealing with topological states, we calculate the entanglement spectrum of the system. For this purpose we partition the chain into a right (R) and left (L) subsystem as |ψ GS = n λ n |ψ L ⊗ |ψ R where λ n are the corresponding Schmidt coefficients for the specific bipartition. The entanglement spectrum is then defined as the set of all the Schmidt coefficients in logarithmic scale ε n = −2 log λ n and is degenerate for phases with topological properties in one dimension [75]. We find that ε n are degenerate near the chain center when the bipartition is drawn across a strong bond, while it is non-degenerate at the weak bonds. In Fig. 5(b), we display the entanglement spectrum for trivial and topological ground states of the BI phase for N s = 60, when ε n are measured across the bipartition at the chain's center. The entanglement spectrum, together with the density pattern and the behavior of string and parity order parameters, provide convincing proof of the topological character of the BI phase. Furthermore, to show that the BI phase is stable in the thermodynamic limit, we consider the entanglement gap, ∆ε = ε 1 − ε 0 , for different system sizes. Fig. 5(c) presents the variation of ∆ε across BSF-BI-SF phases for fixed t/U = 0.05 -confirming the stability of SPT BI phase in the thermodynamic limit. We perform similar analysis with different values of t/U too. Such a finite-size analysis with the entanglement gap also confirms that the phase boundaries predicted for N s = 60 in Fig. 2 remain stable with increasing system-size.
In order to reveal the bulk-edge correspondence, we determine the many-body Berry phase [76] and show that it is Z 2 -quantized in the BI phase. We first note that, because of the strong interactions, the winding number or the Zak phase [56,77,78] is not a good topological indicator in our case. Therefore, we follow the original proposition of Hatsugai [79] and determine the local many-body Berry phase, which is a topological invariant playing the role of the local "order parameter" for an interacting case [79]. For this purpose, we introduce a local twist t → te iθn in the Hamiltonian (3), such that the system still remains gapped in the BI phase. Then the many-body Berry phase is defined as where ψ θn 's are the ground states with θ 0 , θ 1 , ..., θ K = θ 0 on a loop in [0, 2π]. Here, we consider the local Berry phase corresponding to a bond by giving the local twist in tunneling strength t → te iθn only on that particular bond, and take K = 10. Since the ground state manifold is (quasi-)degenerate, we need to put small local terms to distinguish between the different states in the manifold in order to compute local many-body Berry phases. Specifically, we add ∓0.02 b † 1b 2 +b † Ns−1b Ns + H.c. respectively for the trivial and topological states to the twisted Hamiltonians. However, in case of topological states, two edge states (bottom two rows of Fig. 4) are exactly degenerate and cannot be distinguished by the local term mentioned above. In that case, to calculate the many-body Berry phase we put one extra particle, i.e., N s /2 + 1 bosons in total, so that we have only one edge state where both edges are localized with an extra +1/2 particle. Similarly, we could have reduced one particle (N s /2−1 bosons) where the unique edge state would have an extra −1/2 particle localization on both the edges. The local Berry phases γ (K) 's are displayed in Fig. 5(d) for the system size N s = 60. Similar to the entanglement spectrum, we find γ (K) = π for the strong bonds, while γ (K) = 0 on the weak bonds.

Discussion
The BI phase of this model is a reentrant phase. It separates the SF phase, where correlated hopping is suppressed by quantum fluctuations, from the BSF phase, where correlated tunneling is dominant and single-particle tunneling establishes correlations between the bonds. We have provided numerical evidence that the emerging topology is essentially characterized by the interplay between quantum fluctuations and correlated tunneling. Interactions are here, therefore, essential for the onset of topology. Their global nature is at the basis of the sponta-neous symmetry breaking that accompanies the onset of this phase and which induces an asymmetry between bonds. In this respect, it is reminiscent of the Peierls instability of fermions in resonators [31], where the topology is associated with the spontaneous breaking of Z 2 symmetry. Differing from that case, where photon scattering gives rise to a self-organized superlattice trapping the atoms, in our model photon scattering interferes with quantum fluctuations. Like in [31], gap and edge states can be measured in the emitted light using pump-probe spectroscopy. The single-particle structure factor may be directly accessible by the time-of-flight momenta distributions [64,65] enabling the detection of insulator-superfluid phase transition. The two combined measurements of the cavity output and of the structure factor shall provide a clear distinction between the BI, SF and BSF phases.
We observe that the global long-range interaction of this model inhibits the formation of solitons. For other choice of periodicity, and thus of the form of op-eratorB, one could expect glassiness associated with the formation of defects [39], whose nature is expected to be intrinsically different from the one characterizing short-range interacting structures.
To conclude, we have presented a new paradigm of topological states formation via interference between single particle dynamics and interaction induced hopping.
A Coefficients of the extended Bose-Hubbard model To fix the notation, let us consider N atoms of mass m confined within an optical cavity in a quasi-onedimensional configuration (almost) collinear with a one-dimensional optical lattice created by light with wave number k L = 2π/λ, which may be different from k -the wave number of the cavity field. The optical lattice is created by the trapping potential, V 0 making the atomic motion effectively onedimensional. In our calculations we take V = 50E R , where E R = 2 π 2 /2ma 2 is the recoil energy and a denotes the periodicity of the optical lattice. The Wannier basis in the lowest band is W i (x, y, z) = w i (x)Φ 0 (y, z) with w i (x) being the standard onedimensional Wannier function centered at site i and Φ 0 (y, z) -a two-dimensional Wannier function for the transverse deep lattice. For our purposes, Φ 0 (y, z) is essentially equivalent to a Gaussian with width σ = a 2 π 2 E R /V . After adiabatically eliminating the cavity field one gets an effective Hamiltonian describing atomic motion in terms of atomic creation/annihilation operatorsb † i /b i for atoms at site i is decomposed into the sumĤ =Ĥ BH +Ĥ Cav [39,47], where the motion in the lattice and atom-atom interactions are described by the standard Bose-Hubbard Hamiltonian (we assume contact interactions and neglect densitydependent tunneling terms known to be small for such interactions [80]) are given by the integrals The rescattering due to the cavity mode leads to a long range "all-to-all" interaction terms expressible aŝ (13) We rewrite the integrals above as Here, β = k/k L is the ratio between k, the wavenumber of the cavity mode, and k L , the wavenumber corresponding to the optical lattice, and φ is a phase in the mode function. In our work we consider β = 1, and note that arbitrary values of β would lead to a quasiperiodic Hamiltonian [81,82]. For β = 1 and for i = j the magnitude of the integral becomes independent of i, that can be written as Please note here that we have defined y ij (Eq. (14)) coming from Eq. (13) with a minus sign to fix y 11 = z to be positive, as in our convention the lattice index starts from i = 1. For non-diagonal y ij , we observe that due to localization of Wannier functions i = j ± 1 terms may be only significant ones. For our choice of β = 1, the magnitude of the integral again becomes independent of i, and can be written as ThenĤ Cav may be put in the form, Typically, one assumes no phase difference (φ = 0) between the optical cavity and the external optical lattice, i.e., the lattice sites are located at the antinodes of the cavity mode. Then |z| = 0 and |y| ≈ 0, and the terms proportional to y are dropped off leading to the standard case considered in the past [83]. The importance of additional terms was noticed already in [41] where an identical Hamiltonian is obtained for a slightly different arrangement of the cavity and external optical lattice. Here, we have focused on a immediate vicinity of φ = π/2. This corresponds to the atoms being trapped at the nodes of the cavity mode, where z vanishes (see left panel of Fig. 6). However, as φ starts to deviate from π/2 the y term rapidly decreases and z term becomes significant (see right panel of Fig. 6). Note that the quadratic form ofĤ Cav is responsible for the long-range character of the couplings. SquaringD leads then to all-to-all density-density interactions, responsible for a spontaneous formation of density wave phase for sufficient U 1 [84]. For the case considered by us, z ≈ 0 and B 2 term leads to the all-to-all long-range correlated tunnelings alternating in sign. Throughout the paper, we have fixed V = 4E R resulting in y = −0.0658 and z = 0 for φ = π/2.
In Fig. 7, we plot the order parameters and the phase diagram in the vicinity of φ = π/2 for t/U = 0.05. As φ starts to deviate from π/2, y starts to diminish and z becomes increasingly larger (see right panel of Fig. 6). As a result, BI phase is replaced by a more standard density wave (DW) phase [84], when φ becomes sufficiently different from π/2. In the DW phase O DW as well as O S and O P are both non-zero, while the structure factor vanishes.

B Numerical implementation
We use standard matrix product states (MPS) [62,63] based density matrix renormalization group (DMRG) [60,61] method to find the ground state and low-lying excited states of the system with open boundary condition, where we employ the global U (1) symmetry corresponding to the conservation of the total number of particles. For that purpose, we use ITensor C++ library (https://itensor.org) where the MPO for the all connected long-range Hamiltonian can be constructed exactly [85,86] using AutoMPO class. The maximum number bosons (n 0 ) per site has been truncated to 6, which is justified as we only consider average density to be ρ = 1/2.
We consider random entangled states, |ψ ini = are random product states with density ρ = 1/2, as our initial states for DMRG algorithm. The maximum bond dimension of MPS during standard two-site DMRG sweeps has been restricted to χ max = 200. We verify the convergence of the DMRG algorithm by checking the deviations in energy in successive DMRG sweeps. When the energy deviation falls below 10 −12 , we conclude that the resulting MPS is the ground state of the system.
To obtain low-lying excited states, we first shift the Hamiltonian by a suitable weight factor multiplied with the projector on the previously found state. To be precise, for finding the n th excited state |ψ n , we search for the ground state of the shifted Hamiltonian, where W should be guessed to be sufficiently larger than E n − E 0 .