Haldane insulator in the 1D nearest-neighbor extended Bose-Hubbard model with cavity-mediated long-range interactions

In the one-dimensional Bose-Hubbard model with on-site and nearest-neighbor interactions, a gapped phase characterized by an exotic non-local order parameter emerges, the Haldane insulator. Bose-Hubbard models with cavity-mediated global range interactions display phase diagrams, which are very similar to those with nearest-neighbor repulsive interactions, but the Haldane phase remains elusive there. Here we study the one-dimensional Bose-Hubbard model with nearest-neighbor and cavity-mediated global-range interactions and scrutinize the existence of a Haldane Insulator phase. With the help of extensive quantum Monte-Carlo simulations we find that in the Bose-Hubbard model with only cavity-mediated global-range interactions no Haldane phase exists. For a combination of both interactions, the Haldane Insulator phase shrinks rapidly with increasing strength of the cavity-mediated global-range interactions. Thus, in spite of the otherwise very similar behavior the mean-field like cavity-mediated interactions strongly suppress the non-local order favored by nearest-neighbor repulsion in some regions of the phase diagram.


Introduction
For several decades, the Bose-Hubbard model (BHM) [1] has attracted continued interest. In its most simplistic form, it exhibits two phases in the ground state: a Mott insulator (MI) phase and a superfluid (SF) phase, depending if on-site repulsion or nearest-neighbor hopping dominates. Through the years, quantum Monte-Carlo (QMC) methods contributed greatly to the investigation of quantum critical phenomena. Here, one must especially emphasize path-integral Monte-Carlo [2,3], world-line QMC [4,5], worm-algorithm QMC [6][7][8][9] and, as a derivative, the stochastic Green's function algorithm [10,11].
An interesting question is, whether the HI phase can also occur in the BHM with cavity-mediated interactions, since on the mean-field level it is equivalent to the BHM with NN interactions [45]. In this paper we will address this question with the exact QMC worm-algorithm [6,7].
The paper is organized as follows: in Section 2 we introduce the BHM with extended NN and LR interactions and show that both additional interactions lead to similar terms in the mean-field approximation. Section 3 outlines the QMC worm-algorithm. The measured observables and results for different parameter settings are discussed in Section 4. We account chain lengths up to L = 192.

Model
We consider the one-dimensional extended Bose-Hubbard model with nearest-neighbor short-range and cavitymediated long-range interactions (NNLR-BHM) in the grand-canonical ensemble. It is defined by the Hamiltonian The first term describes the nearest-neighbor hopping between two sites with the hopping strength t.b † i (b i ) are the bosonic creation (annihilation) operators. The second term describes the on-site repulsion (U > 0) andn i is the number operator. The third term defines a nearestneighbor repulsion (V > 0). The fourth term contains the chemical potential µ. We use the grand-canonical ensemble, thus keeping µ fixed and let the total particle number N vary. The cavity-mediated interaction is introduced in the last term. U d is the misbalance parameter, L the chain length and e(o) the sum over all even (odd) sites. We apply periodic boundary conditions. The total misbalance is defined as where D ranges between −N (all bosons on odd sites) and +N (all bosons on even sites). In Dogra et al. [45] it was shown that the mean-field (MF) expressions of the LR and NN interactions are equivalent. In mean-field approximation the decoupled expressions for kinetic energy, NN and LR interactions readb † with the coordination number z and the order parameters ψ = b † i = b j , ϑ = n i − n j and θ = 2 D /L. Hence, on a mean-field level the LR-BHM term in (1) is equivalent to the NN-BHM term: In (4) the transformed LR interaction increases the NN interaction, as intuitively expected. Furthermore, the chemical potential is shifted, such that even for negative values of µ the system can localize in a DW or MI phase and the rescaling of the on-site potential leads to a narrower lobe width.
Ultimately, since the HI phase was already found in Hamiltonian (1) with U d = 0 [28,39,43,48] and since on the MF level NN and LR interactions are equivalent one is lead to ask whether the HI phase exists also for U d > 0. This is the question that we will answer in the following using a QMC algorithm.

Worm-algorithm
To determine the ground state properties of the Hamiltonian (1) we use the quantum Monte-Carlo (QMC) worm-algorithm [6,7,51]. It relies on the Dyson series where the d-dimensional quantum system is mapped onto a (d + 1)-dimensional classical one. The partition function reads with the inverse temperature β andV ninj = n i |V |n j . The off-diagonal part of the Hamiltonian (1) isV = t i,j (b † ib j + h.c.), |n i > are the Fock states of the diagonal Hamiltonian and i are the diagonal energy values of the respective Fock states.
In the worm-algorithm the configuration space is expanded by a creator and annihilator pair in the form b i (τ ) = (e τ ab i e −τ a ) and vice versa forb † i (τ ). One operator is assigned as the worm head and the other one as the worm tail. The head moves through the given state and can create, delete or relink vertices, where bosons hop from one site to a neighboring one. When it reaches the tail the worm gets deleted.
Between a head and tail, the total particle number will be increased or decreased in comparison to the initial state. Thus the worm-algorithm performs grand-canonical update steps. When the update procedure is complete and the worm deleted, one can calculate canonical observables by importance sampling with C denoting states with fixed total particle numbers. We consider chain lengths up to L = 192 and an onsite repulsion of U = 1. The inverse temperature was set to β = 128 as comparisons with lower temperatures showed no sufficient difference for the obtained order parameters.

Measured observables
In the NNLR-BHM various observables are of interest. The particle density is ρ = i n i /L and the superfluid density can be calculated via [2,3] where W is the winding number, which is the difference of bosons crossing over one side of the periodic boundary conditions minus the crossing over the other side. The density-density correlation and structure factor are defined as For k = π the structure factor gives the misbalance order parameter θ = 2 D /L. These order parameters would be sufficient to distinguish between the Mott insulator (MI) phase, superfluid (SF) phase, density wave (DW) phase and supersolid (SS) phase.
For a fixed density ρ = 1 and large U and V values, the site occupation is practically restricted to 0, 1 and 2. As a result, the difference between particle number and density δn i =n i − ρ has the same Eigenvalues as theŜ z spin operator from the Spin-1 Heisenberg chain (−1, 0, 1) and thus both models are similar.
(12) while the numbers represent δn i with an undetermined amount of 0 between each +1 and −1. On the other hand, the MI is a dilute gas of particle-hole pairs, thus no global ordering emerges but local fluctuations which can result Table 1. Order parameters at ρ = 1 for different phases. (*) In finite chain lengths, the superfluid density can attain non-zero values. ρs S(π), θ Os(L/2) Op(L/2) where two consecutive −1 or +1 emerge and counteract the global ordering [48]. Therefore in finite systems, it is possible that the particles (+1) and holes (−1) wind around the chain, which results in a non-zero superfluid density. However, in the limit L → ∞, ρ s disappears for the HI and MI phases.
We consider three different parameter regimes of the NNLR-BHM. In the first part, we set the NN interaction to V = 0.75 and the LR interaction to U d = 0 and compare our results with [28]. In the second part we fix V = 0 and U d = 0.6. Finally, we include NN as well as LR interactions and compare the behavior of obtained phases to the previous cases. Therefore, we set V = 0.75 and increase U d in size.
For all cases, we focus on the ρ = 1 lobe of the µ/U − t/U diagram and determine the phases for commensurate filling to see whether a HI phase occurs. So, we tune the chemical potential carefully such that the particle density ρ equals 1 and then perform measurements of the order parameters.

Results for the NN-BHM
In this subsection, we neglect the LR interaction, i.e. (U d = 0) and consider only the NN extended BHM with V = 0.75. Since zV > U the DW(2,0) phase appears instead of the MI(1) phase. Furthermore, calculations for the t = 0 case yield that the DW(2,0) phase occurs for 1 < µ/U < 2.
The phase diagram for the ρ = 1 lobe is depicted in Figure 1. There are four different phases: the DW(2,0) phase in the red lobe, surrounded by a SS phase up to the blue line where the transition to the SF phase occurs. Between the green lines, the HI phase is present and transits into the DW phase at around t = 0.22. Our results are in good agreement with [28], where the 1D phase diagram with NN potential was calculated with density matrix renormalization group and stochastic Green's function methods.
Next, we fix µ to a linear equation depending on t, like depicted in Figure 1 and increase t from 0 to 0.35, thus traversing the whole ρ = 1 lobe to the tip and into the HI phase. The results for the order parameters are shown in Figure 2. Here, the DW phase, in which ρ s is zero and all other parameters are non-zero, persists up to t ≈ 0.22  where the transition to the HI phase occurs. There, S(π) and O p drop to zero while O s decays but stays finite. The superfluid density is non-zero but size dependent and vanishes for larger system sizes until the transition to the SF phase occurs. For the SF phase all order parameters are zero except the superfluid density.
Next, we study the finite size behavior of the order parameters in the different phases. In Figure 3  For the DW phase all order parameters stay constant. As expected S(π), O s and O p attain finite values while ρ s = 0 for all chain lengths. In the HI phase S(π) and O p are finite for small sizes and become zero for infinite L. The string order parameter persists for larger sizes and approaches a finite value. The superfluid density decreases for increasing lengths and vanishes completely for a large (but finite) L.
This behavior underlines the breaking of the hidden Z 2 symmetry. For small sizes the particles and holes can potentially wind around the chain, leading to a finite superfluid density value. This fluctuation is dependent on the amount of particles and holes in the system, meaning when there are a lot -like close to the DW phase -the fluctuations get small and when there are only a few, the fluctuations get high. When L is smaller than this fluctuation length, winding happens and the superfluid density is greater than zero. Otherwise, there exists a finite chain length where no winding and thus no superfluid density exists any more.
The order parameters in the SF phase behave opposite to the DW phase. The superfluid density is non-zero and does not vanish for infinite sizes while S(π), O s and O p are very small for tiny lengths and become zero for L → ∞.
We can visualize the finite size effects in the HI phase also by the Green's function The worm algorithm can directly calculate the Green's function, which is a degree of spatial movement in the system and thus correlated to the winding number [3].  Figure 4 depicts the Green's function G(r) for various chain lengths in the different phases. When the worm moves through the extended configuration space the distance between worm head and tail varies with every Monte Carlo move. In the MI and DW phase, this movement is rather restricted and head and tail stay close to each other, which lead to an exponentially decaying Green's function. In the SF phase bosons become delocalized implying that both worm ends can be arbitrarily far from each other. The Green's function in the SF phase is expected to decay algebraically in one dimension, but in a finite system G(r) has a minimum at G(L/2) due to the periodic boundary conditions, as is visible in Figure 4. In the HI phase G(r) behaves similar to the SF phase for small system sizes, but for larger system sizes the winding of the worm becomes unlikely and -as in the MI and DW phase -the decay of the Green's function approaches an exponential form.
The Fourier transformation of the Green's function yields the momentum distribution and the k = 0 mode gives the condensate fraction, which is experimentally accessible [53]. Note that due to the algebraic decay of the Greens function in one dimension the condensate fraction vanishes in the thermodynamic limit, in contrast to the superfluid density.

Results for the LR-BHM
Now set V = 0 and U d = 0.6. For t = 0 it is straightforward to determine the different phases of equation (1). First, since 2U d > U , all phases are DW(X,0) phases with X being any integer number. That means every second site is occupied by X particles, while all other sites are empty. The transition from vacuum to the DW(1,0) phase is at µ 0,1 = −0.3 and at µ 1,2 = 0.1 the DW(2,0) phase starts. Every DW phase has a width of ∆µ = 0.4, so the next transitions are at µ 2,3 = 0.5, µ 3,4 = 0.9 and so on. On the basis of the MF analysis resulting in equation (4) we can obtain a first guess of the approximate shape of the phase diagram of the LR-BHM and introduce a set of rescaled parameters:Ũ = 0.4,μ = µ + 0.3 and V = 0.3, thus the ratioṼ /Ũ = 0.75 is identical to the NN case. The rescaled on-site repulsion accounts for the same lobe width ∆µ =Ũ and the rescaled chemical potential is shifted by the equal amount as discussed above for the t = 0 case. Then the DW(2,0) phase exists in the interval µ/Ũ ∈ [1,2], the same range as in the NN-BHM (Fig. 1).
Therefore, we present the results for the LR-BHM in the ratio ofμ/Ũ and compare it directly to the NN-BHM case. In Figure 5 the DW(2,0) lobe is depicted. In comparison with the phase diagram of the NN-BHM, Figure 1, we see that no HI phase emerges at the tip of the DW lobe, but a MI phase instead. Otherwise, the DW phase is broader and its tip shifted downwards.
Next, analogous to the NN-BHM case, we keep ρ = 1 constant by fixingμ to a linear function depended on t and increase t from 0 to 0.25. Our results are shown in Figure 6. In the DW phase the structure factor, string and parity order parameters are non-zero, while the superfluid density vanishes. Approaching the transition point at around t = 0.23 the structure factor and string order parameter drop to zero, while the parity order parameter persists. Also, the superfluid density increases but is strongly dependent on the system size as for larger sizes the superfluid density tends to zero.
The behavior of the ρ = 1 phase transition in the LR-BHM is similar to the NN-BHM (Fig. 2), whereby the MI phase replaces the HI phase. Here, not the string order parameter persists, but the parity.
We show the finite size scaling of the LR-BHM with commensurate filling in Figure 7. The order parameters in the DW phase behave as in the NN-BHM case. Otherwise, in the MI phase only the parity operator is non-zero, Fig. 6. Order parameters as a function of t/U forμ(t) = −1.5 t + 2.125Ũ and U d = 0.6 in the LR-BHM at ρ = 1. Depicted are the DW(2,0) phase left to t = 0.23 and the MI phase right to it. For the DW phase the structure factor, string order parameter and parity order parameter are non-zero, while the superfluid density is zero. In the MI phase the parity order parameter remains greater than zero and the superfluid density increases but approaches zero for larger chain lengths. where structure factor, string order parameter and parity order parameter have finite values, while the superfluid density remains zero. Bottom: MI phase for t = 0.24, where the structure factor and string order parameter tend to zero, whereas the parity operator has a finite value and the superfluid density vanishes for a large, but finite chain length.
while the superfluid density becomes zero at a large, but finite L.
The existence of a MI phase at the tip of the DW phase was not found in two (or more) dimensional systems. The emergence of a MI phase can be understood in the following way: for the DW(2,0) phase a large long-range interaction exists, preventing site occupation fluctuations from the underlying checkerboard structure. On the other hand in the MI phase for small occupation imbalances, the resulting long-range interaction is proportional to 1/L, thus negligible for large L. This effect is based upon the description of the HI and MI phases as explained above [48].

Results for the NNLR-BHM
In this section, we analyze nearest-neighbor and longrange interactions simultaneously. As discussed in the subsections before, the HI phase exists in the NN-BHM while it is absent in the LR-BHM. Hence, we start with the NN-BHM and add increasing long-range interactions to see if the HI phase gets destructed. Figure 8 shows the evolution of the HI phase with the inclusion of LR interaction. Already for LR parameters U d one order of magnitude smaller than the NN parameter V the HI phase disappears quickly. We see that the phase transition point of t = 0.22 for the NN-BHM moves to t = 0.25 for U d = 0.02 and to t = 0.275 for U d = 0.04. The phase transition between HI and SF phases is not affected by the increasing LR parameter.
The inclusion of a system wide LR coupling has a strong influence on the DW phase, while it is negligible for the HI and SF phases. Thus, the DW phase expands and supersedes the HI phase. For strong enough LR parameter strength the HI phase disappears completely as expected from the LR-BHM results.
We have not checked the behavior of the MI phase in the LR-BHM when increasing the NN interaction, but we assume this phase to vanish analogously to the HI phase above. Since an additional NN interaction does not increase the energy of the DW(2,0) phase but the energies of the MI and SF phases, we expect the DW(2,0) phase to expand and the MI to shrink for increasing NN interactions.

Conclusions
In this paper, we investigated the extended 1D Bose-Hubbard model with a nearest-neighbor interaction, a cavity-mediated long-range interaction and both combined. For the NN-BHM we confirmed earlier results that a Haldane insulator phase exists. For the LR-BHM we found the absence of a HI phase at the tip of the DW lobe and its replacement by a MI phase. The reason is the global long-range interaction preventing site occupation fluctuations in the DW phase while it is possible in the MI phase. For the NNLR-BHM we increased the longrange parameter gradually and observed the HI phase to shrink as it is replaced by a growing DW phase. The LR parameter was one order of magnitude smaller than the NN parameter, showing the instability of the HI phase against a global ordering via cavity-mediated interactions.
Furthermore, we conclude that the NN interaction prevents the creation of a MI phase at the tip of the DW lobe due to the commensurate filling of all sites, while the LR interaction suppresses a HI phase at the tip since the specific global order, necessary to form a HI phase, becomes dominated by the long-range ordering effect of the cavity-mediated interactions. Hence, neither HI nor MI phases exist in the NNLR-BHM with strong LR and NN couplings.
Open access funding provided by Projekt DEAL. We like to thank Benjamin Bogner for his comments on the paper and many beneficial discussions and Chao Zhang for her support with the worm-algorithm and helpful suggestions.

Author contribution statement
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.