Phase diagram of the extended Bose Hubbard model

By means of the Density Matrix Renormalization Group technique, we accurately determine the zero-temperature phase diagram of the one-dimensional extended Bose Hubbard model with on-site and nearest-neighbor interactions. We analyze the scaling of the charge and of the neutral ground-state energy gaps, as well as of various order parameters. In this way we come to an accurate location of the boundaries between the superfluid and the insulating phases. In this last region we are able to distinguish between the conventional Mott insulating and density-wave phases, and the Haldane Insulator phase displaying long-range string ordering, as originally predicted by E.G. Dalla Torre, E. Berg and E. Altman in Phys. Rev. Lett. 97, 260401 (2006).


Introduction
In the last decade, ultracold atomic gases loaded in optical lattices successfully established as an excellent setup to probe the equilibrium, as well as the out-ofequilibrium physics of strongly correlated quantum systems. The great advantages of these experimental setups are essentially related to two aspects. On one hand, they have an extremely high degree of flexibility: besides the ability to address different geometries, and to deal with Bosonic and with Fermionic species, they admit the possibility to manipulate the underlying Hamiltonian system parameters to a large extent. Moreover, the remarkably high degree of isolation from any environmental source of decoherence opened up entirely new scenarios in the observation of genuinely many-body quantum phenomena [1].
The paradigm model to describe cold Bosonic atoms trapped in an optical lattice is obtained by combining the kinetic energy in the lowest band with the on-site repulsion arising for sufficiently deep lattices. This leads to the celebrated Bose-Hubbard model (BHM) [2]. The rich physics of the BHM stems from the competition between the kinetic energy J, which is gained by delocalizing particles across the lattice in an extended Bloch state, and the repulsive on-site interaction U, which disfavors having more than one particle per site. When the kinetic energy term dominates, the system is in a coherent Superfluid (SF) phase; on the other hand, repulsive interactions tend to favour a Mott insulating (MI) phase [3].
The recent advances in manipulating magnetic atoms and molecules with large dipole momentum make it possible to achieve longer-range interactions, which can be accurately tuned as well, thus permitting to probe the interplay between strong correlations and charge-ordering effects [4]. Dipolar bosons confined in optical lattices are typically described by an extended version of the BHM, the so called Extended Bose Hubbard Model (EBHM), which also includes a two-body non local repulsive term typically decaying as r −3 with the distance r.
Interestingly the presence of long range interactions noticeably enriches the phase diagram of the BHM, for example leading to a stabilization of a peculiar insulating phase, named the Bosonic Haldane insulator (HI). This gapped phase presents some analogies with the well known Haldane phase in integer spin chains [5]; it does not break the translational symmetry of the lattice and is characterized by an underlying hidden order, that is a non-trivial ordering of the fluctuations which appear in alternating order separated by strings of equally populated sites of arbitrary length [6,7,8]. More recently a mean field analysis suggested the existence of various supersolid phases and a series of non-trivial density-wave states with an increased ground-state degeneracy, which can be stabilized in the strong coupling regime, but with a relatively weak onsite interaction [9]. It has been also shown that the presence of disorder in the on-site potential [10], as well as occupation-dependent hopping and pair tunneling terms arising from dipolar interactions [11], can drive important modifications in the phase diagram leading to the stabilization of a plethora of novel quantum phases, such as structured insulating states, Wigner crystals, or pair-supersolids.
The aim of this paper is to work out the zero-temperature phase diagram of the pure one-dimensional (1D) EBHM, by means of the numerical Density Matrix Renormalization Group (DMRG) technique [12], in order to accurately determine the phase diagram first proposed in Refs. [6,7] and extend their preliminary results. This method has been already revealed its great potentialities in the study of ground states of 1D strongly correlated bosons on a lattice, providing very accurate results for these kinds of systems (see e.g. Refs. [13,14,15,16,17,18] and references therein).
The paper is organized as follows. In Sec. 2 we introduce the model we are going to study, and qualitatively discuss the various emerging quantum phases. We then present, in Sec. 3, all the relevant quantities we will address: i) three kinds of local and non local correlators, which will enable us to locate all the different phases identifying suitable order parameters; ii) the ground state charge and neutral energy gaps, which will serve to discriminate critical phases as well as insulating QPT points. A phase diagram of the model is shown in Sec. 4, followed by a detailed discussion of the finite-size scaling of the various order parameters (SubSec. 4.1) and of the charge and energy gap (SubSec. 4.2), up to a few hundreds of sites. Finally, in Sec. 5 we draw our conclusions.

The model
In the following we will focus on the 1D extended Bose Hubbard model (EBHM), with on-site and nearest-neighbor repulsive interactions. The sharp cutoff of the r −3 interaction range, although presenting quantitative differences with the standard longrange EHBM, does not qualitatively alter its physics: here b † j , b j are creation and annihilation operators of bosons on site j, J denotes the hopping strength, while U and V are the on-site and the nearest-neighbor interaction strengths. Hereafter we set the energy scale by taking J = 1, and work in units of = 1.
We will concentrate on the model in Eq. (1) at zero temperature and at integer fillingn, which is known to exhibit a quite rich phase diagram with various quantum phases, ranging from the superfluid (SF) for low U, V interaction strengths, to insulating phases in the opposite regime. Depending on the relative strength of U and V , a Mott insulator (MI) or a density-wave state (DW) can form, the first one establishing with large values of U, while the second one with large V . In between the two conventional insulating phases, a peculiar gapped phase emerges: the Haldane insulator (HI).
The phase diagram and the phase transitions of such kind of model, describing bosons loaded in optical lattices with short range interactions, have been addressed since long time [3]. Emphasis was put on the experimentally accessible SF-MI transition [19]. In the 1D case and with a constant integer filling, the transition is of the Berezinskii-Kosterlitz-Thouless type and was extensively analyzed with DMRG in Refs. [13,14]. The same authors also demonstrated that, in presence of nearest-neighbor interactions, there is no normal or supersolid phase, but a direct phase transition from the DW to the SF phase occurs. Only a few years ago, by means of a mean field theory describing lowenergy excitations in the continuum limit, in Refs. [6,7] it was realized that extended interactions can stabilize a further HI phase, possessing hidden order that is revealed by non-local string correlations. In those papers the DMRG was employed to extract some preliminary quantitative results and draw the resulting phase diagram in 1D. Further DMRG calculations have been also used to address closely related effective spin-1 chains, including finite-size scaling studies [20,21], the effects of trapping potentials [22], and to study the entanglement spectrum in order to probe the degeneracies in the HI phase [23].
Here we numerically probe the EBHM over all its parameters space, corroborating the investigations previously done in the literature and providing an accurate description of its complete phase diagram at integer filling. In order to locate the different quantum phases of Eq. (1) we resort to the DMRG algorithm with open boundary conditions ‡. Our code is based on a variational ansatz using Matrix Product States (MPS), in the restricted subspace of integer filling [12]. ‡ The choice of open boundary conditions is popularly used in the DMRG community, and it is dictated by the fact that it has been proven that they perform more efficiently than periodic boundaries [12]. We will take special precautions for border effects emerging in our simulated finite-size systems.

Correlation functions, order parameters and energy gaps
The nature of the different quantum phases emerging in the EBHM, as well as all the Quantum Phase Transition (QPT) points determining the various phase boundaries, are revealed by the behavior of certain two-point correlation functions. The following three correlators are able to univocally distinguish between all the emerging quantum phases: C string (r) = δn j e iπ j≤k<j+r δn k δn j+r , where δn j = n j −n denotes the boson number fluctuations from the average filling (hereafter we will always consider the casen = 1). The long range off-diagonal order, typical of superfluid states, manifests itself into a power-law decay to zero of C SF (r); on the other hand, in the insulating phases the absence of such ordering is characterized by an exponential suppression of C SF (r) with r. The large-r limit of the second correlator, C DW (r), identifies a staggered diagonal order which naturally emerges in the DW phase: atn = 1, large values of V in Eq. (1) tend to stabilize such ordering, by inducing a pattern . . . − 2 − 0 − 2 − 0 − . . . for the occupation number. Such phases are then characterized by a finite DW order parameter O DW ≡ lim r→∞ C DW (r). Finally, the two-point correlator C string (r) contains a non-local string correlation function accounting for the number of bosons that are present in between the two points. For large values of r, this has been shown to define an order parameter for the HI, namely O string ≡ lim r→∞ C string (r), which distinguishes from the usual MI though such highly non-local ordering. On the contrary, in the incoherent MI phase all the correlation functions in Eqs. (2)-(3)-(4) vanish exponentially.
In addition to the previously defined correlation functions, we also analyze the energy gap. The nature of the different quantum phases, as well as the QPT points are indeed sensitive to it. In particular, we will focus on the neutral (∆E n ) and on the charge (∆E c ) gaps, respectively defined as: where E (k) m indicates the energy of the k-th excited state (k = 0 labelling the ground state, k = 1 the first excited state, and so on), for a system with a given number of m bosons in L lattice sites §. In practice we run our DMRG simulations in the canonical ensemble. The charge gap is evaluated by running three instances of DMRG, each of them separately targeting the ground state of n = {L, L ± 1} particle number § The case m = L corresponds exactly to the integer filling sectorn = 1, in which the neutral gap is evaluated. The charge gap in Eq. (6) expresses the difference in the energy costs ∆E (0) ± to add and to subtract a single particle from the ground state at unitary filling, i.e. ∆E c = ∆E (0) sectors [13,14]. The neutral gap is extracted by running only one instance of DMRG, in which both the ground state and the first excited state are targeted .
Unfortunately, the charge gap ∆E c only detects particle and hole excitations, and it is able to locate the superfluid-to-insulator transition, as well as the MI-HI QPT point. On the other hand, the HI phase can display another type of excitation which in energy goes below the particle-hole excitation, and is revealed only by the neutral gap ∆E n . It turns out that the presence of such low-energy neutral modes becomes crucial for the identification of the HI-DW transition, as we will explicitly show later.
In the next Section we will discuss in details our results. We first show the phase diagram of the EBHM, distinguishing between the various quantum phases of the model. Then we will present our data on the order parameters and the ground-state energy gaps.

Phase diagram
We start the discussion of our results by presenting the phase diagram of Eq. (1) in the (U, V ) plane. Our DMRG simulations were performed for systems with open boundary conditions up to L = 400 sites, keeping at most m = 500 states, and working in a canonical ensemble with a fixed number of bosons. Due to computational reasons, for a fillingn = 1 we admitted a maximum of n max = 3 bosons per site. We checked that the location of the phase boundaries is not qualitatively affected by the values of the various cutoffs ¶. Figure 1 summarizes the results of our investigation. As we stated before, model (1) has four different quantum phases. For large U values, the system is always insulating. In order to distinguish between MI, HI and DW phases, we carefully analyzed the string and the density-wave order parameter, respectively O string and O DW . An accurate finitesize scaling of the O string data allowed us to draw the black line (data are shown as circles) in Fig. 1. A similar analysis was also performed on O DW , producing the red curve in the figure (red squares). For further details in the specific situation with U = 5 we refer to Sec. 4.1. We performed a cross-check of such obtained phase boundaries by looking at the charge and the neutral energy gaps (see Sec. 4.2), obtaining minimal values for the gaps at the QPT points.
For small values of U and not too large V , the system enters a critical SF phase, characterized by long-range off diagonal order. We located the SF by looking at the charge energy gap (blue diamonds). While for small V values the determination of the SF-insulator boundary is quite accurate at relatively small sizes L ≈ 300 (where, by increasing U, at some critical point the charge energy gap undergoes a sudden increase), this becomes harder for larger V . The dashed line in the figure has been plotted as a In general, fixing the system size and the bond-link dimension, the evaluation of ∆E c is much faster and more accurate than that of ∆E n , since only the ground state has to be addressed. ¶ We also performed simulations fixing n max = 2, thus mapping the problem into a spin-1 Heisenberg Hamiltonian with a single-site uniaxial anisotropy, and explicitly found a phase diagram with the same qualitative features as those depicted Fig. 1. We however noticed important quantitative changes, with the HI phase being shifted at lower values of V , and the DW phase enlarged. guide to the eye, and is an interpolation between the blue data in the small-V region and the red DW boundary at large values of V . As a matter of fact, despite the relatively low numerical accuracy that we were able to achieve in that region, our data are indeed compatible with the presence of a SF phase below the HI, for small U values and 1.5 V 3, in agreement with the findings of Ref. [10].

Scaling of the order parameters
In order to determine the HI-DW and the MI-HI phase boundaries, we respectively looked at the correlation functions in Eqs. (3)-(4). In particular we computed the string and the density-wave order parameters by extracting the large-r limit of such correlators. To minimize border effects due to open boundaries, we computed the average expectation values of all the two-point correlators between spins in the middle part of the chain, discarding the outer L/4 sites on both sides. We also studied the dependence on the system size by performing a finite size scaling up to L = 400 system sites.
An example of such analysis is shown in Fig. 2, where we fixed the on-site interaction strength U = 5, and varied V . Both the HI order parameter (filled symbols, on the left In the inset we show how the two order parameters approach the zero value very close to the QPT. Namely, we plotted O string (black circles) and O DW (red squares) respectively at V MI−HI = 2.95 and V HI−DW = 3.525. In both cases we obtained an exponential decay of O with the system size L.

Scaling of the energy gaps
Further analysis of the QPT points has been performed by looking at the neutral and at the charge energy gaps, ∆E c and ∆E n [see Eqs. (5), (6)]. In addition to corroborating the results obtained from the order parameters scaling analysis as we will show in Fig. 4 and Fig. 6, the closure of the energy gaps discriminates between the SF and the insulating phases. This is a method to locate the SF boundaries, which reveals more accurate than directly looking at the off diagonal ordering. As a representative example, in Fig. 3 we show the behavior of the charge energy gap in proximity of the SF-MI transition of model (1), as long as U is increased upon a critical value (see blue diamonds in Fig. 1). Namely, we plotted ∆E c as a function of U, for a fixed value of V and for a system of L = 400 sites. The SF-MI QPT is revealed by a sudden growth of ∆E c at a critical value U ⋆ , which is estimated to be U ⋆ ≈ 3.33 at V = 0 and U ⋆ ≈ 2.3 at V = 1 + . We + The value U ⋆ ≈ 3.33 corresponds to the quantum critical point for the SF-MI transition in the 1D BH model at filling ρ = 1, which has been thoroughly investigated in the past literature (see, for example, Ref. [14]).  remark that, at the sizes that are typically studied, it is easy to discriminate between a power-law decay of C SF (r) deep in the SF phase and its exponential decay deep in the MI, as it is clearly shown in the two insets of Fig. 3. Nonetheless, this becomes much trickier in proximity of the phase boundaries and would require very large sizes to distinguish between the two behaviors [14]. We point out that we are always interested in the bulk energy gaps, even if our DMRG method preferably works with open boundary conditions. Unfortunately the appearance of edge states in the HI phase can interfere with the determination of the bulk gaps. For this reason, we lifted the excitations at the borders by forcedly requiring to have zero particles on the leftmost site of the chain, and two particles on the rightmost site. This corresponds to apply some strong field at the edges. Such a bias also breaks the ground state degeneracy of the DW phase, thus stabilizing the numerical algorithm at large values of V .

Ground-state charge energy gap.
In Fig. 4 we show the charge energy gap as a function of V and for different number of sites L, after fixing the on-site interaction strength U = 5 (same value as in Fig. 2). The non monotonic behavior, with a minimum gap that closes for L → ∞, indicates the presence of an insulating QPT between a MI phase on the left side, and a HI phase on the right side.
A finite-size scaling of the obtained results is performed in Fig. 5. In the left panel we plot the corresponding value of V where the gap reaches its minimum, as a function of the system size. By linearly fitting our numerical data as a function of 1/L, we are able to extrapolate a value at the thermodynamic limit that is equal to V L→∞ min ≈ 2.98. Under numerical accuracies, this is accordance with the value V ≈ 2.95 of the MI-HI QPT point that was found from the analysis of O string at U = 5. The right panel clearly shows that the minimum charge energy gap ∆E c , which is reached at the MI-HI transition, drops to zero as a power law with the size, following the behavior ∆E c ∼ 1/L. This kind of behavior at the QPT transition can be explained within the conformal field theory approach, where it has been proven that finite-size corrections to the energies of low-lying gapless excitations behave as 1/L [24,25].

Ground-state neutral energy gap.
Let us now concentrate on the neutral gap ∆E n , which is plotted in Fig. 6 for the same Hamiltonian parameters of Fig. 4. Besides the QPT between the MI and the HI located at V ≈ 3, another critical point at V ≈ 3.5 emerges, separating the HI and the DW phase (see also Fig. 1). In the MI phase and in part of the HI phase, the particle and hole excitations (detected by ∆E c ) are the lowest-energy excitations. As it is seen from the figure, this is no longer the case for V > V ⋆ ≈ 3.3, starting from the middle of the Haldane phase. The presence of a marked cusp at V ⋆ indeed indicates a modification in the nature of the lowest excitation, which becomes of a different kind. The charge gap is not able to detect such excitation, and the transition to the DW phase, occurring   at V > V ⋆ . On the contrary, the non monotonic behavior of the neutral gap at V * and its dependence on the size, clearly indicate it through a gap closure.
In the inset of Fig. 6 we analyze the behavior of the minimum neutral gap at the HI-DW transition, showing that, similarly to the minimum charge gap at the MI-HI transition, it drops to zero with the size as 1/L [24,25].
Finally we remark that, for the finite sizes we considered, the neutral gap does not perfectly coincide with the charge gap at V < V ⋆ . This slight discrepancy is barely visible from a direct comparison of Fig. 4 and Fig. 6, and comes from boundary effects which are more important for small sizes, and split the difference in the energies of the first charge and neutral excitations. Differences are also induced by the bias added to lift the border excitations.

Conclusions
In this paper we worked out an accurate phase diagram of the extended Bose Hubbard model in one dimension, using the DMRG algorithm with open boundary conditions. We located the superfluid, as well as the Mott insulating, the density wave and the Haldane insulating phases, by performing an up-to-date finite-size scaling with systems of up to L = 400 sites. Our analysis involved the order parameters for the HI and the DW, and the study of the ground state charge and neutral gaps. We point out that it would be tempting to analyze also the bipartite fluctuations of the boson number in a sub-block of the whole lattice system; this method can be employed within standard DMRG simulations and without any additional computational effort. Being able to clearly distinguish between critical and gapped phases, it could provide more accurate data for the location of the superfluid region at large values of V [18,26].
Besides giving rigorous predictions for future investigations of strongly correlated quantum phases in cold gases of dipolar atoms, our quantitative study is also relevant in the context of quantized transport through adiabatic pumping [27]. The model we discussed here indeed supports two distinct insulating phases separated by a critical point: the Mott and the Haldane insulator. One can open a gap at such point by breaking the inversion symmetry of the system (for example, via a correlated tight binding hopping for the bosons). Therefore, by suitably tuning the various parameters of the Hamiltonian, it is possible to adiabatically encircle the MI-HI critical point. As a consequence, the non trivial topology of the isolated gapless point inside the loop would induce the transport of one single boson through the chain [28].