{'fraction_non_alphanumeric': 0.029275178494711646, 'fraction_numerical': 0.011435936253349972, 'mean_word_length': 4.316150081566069, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The one-dimensional Bose-Hubbard Model with nearest-neighbor interaction

We study the one-dimensional Bose-Hubbard model using the Density-Matrix Renormalization Group (DMRG).For the cases of on-site interactions and additional nearest-neighbor interactions the phase boundaries of the Mott-insulators and charge density wave phases are determined. We find a direct phase transition between the charge density wave phase and the superfluid phase, and no supersolid or normal phases. In the presence of nearest-neighbor interaction the charge density wave phase is completely surrounded by a region in which the effective interactions in the superfluid phase are repulsive. It is known from Luttinger liquid theory that a single impurity causes the system to be insulating if the effective interactions are repulsive, and that an even bigger region of the superfluid phase is driven into a Bose-glass phase by any finite quenched disorder. We determine the boundaries of both regions in the phase diagram. The ac-conductivity in the superfluid phase in the attractive and the repulsive region is calculated, and a big superfluid stiffness is found in the attractive as well as the repulsive region.


I. INTRODUCTION
At zero temperature superfluid-to-insulator phase transitions can be found in bosonic lattice systems 1 . Experimental realizations of such systems are superconducting islands or grains connected by Josephson junctions, in which the relevant particles are the bosonic Cooper pairs. The transitions at zero temperature belong to the class of quantum phase transitions that are not driven by thermal, but by quantum fluctuations. These quantum fluctuations are controlled by system parameters like the charging energy of the superconducting islands and the Josephson coupling between them. Depending on such parameters, the system can assume different forms of long-range order. In two dimensions superfluid-toinsulator phase transitions at zero temperature were observed in thin granular films [2][3][4] and in fabricated Josephson junction arrays 5,6 .
Recently experiments were carried out in one dimensional systems of fabricated Josephson junctions. In chains one junction wide and 63, 127 and 255 junctions long with a tunable Josephson coupling a superfluidinsulator transition was observed 7 . In another experiment with fabricated Josephson junctions, long and narrow arrays formed an effectively one-dimensional lattice for lattice fluxes formed by Cooper pairs 8,9 . The density of these fluxes was controlled by an external magnetic field perpendicular to the array, and for small ratios of Josephson coupling to charging energy, insulating charge density wave phases with densities 1/3, 1/2, 2/3, 1.. were found.
In all of these systems the relevant particles, the Cooper pairs or the lattice fluxes, are, at least approximately, bosonic. In this paper we will study the phase diagram of strongly interacting bosons on one-dimensional lattices. In both of the above experiments on one dimensional Josephson junction arrays the range of the interac-tions were reported to be several sites long. To study the effect of longer-ranged interactions we will take nearestneighbor interactions into account.
In the presence of on-site interactions only, Mottinsulators are found at integer densities, surrounded by a superfluid phase 1 . Additional nearest-neighbor interaction leads to charge density wave phases at half-integer densities. At the transition from the charge density wave phase to the superfluid phase two forms of long range order are involved: charge density wave order and superfluid order. If there is a direct phase transition from the charge density wave phase to the superfluid phase, one type of order appears at the same point where the other vanishes. In two and higher dimensions an intermediate phase, called the supersolid phase, that has both types of order, was found between the charge density wave phase and the superfluid phase in theoretical models 10,11 . For experimental systems the existence of supersolids is still controversial 12 . In one-dimensional bosonic models supersolids have not been found so far 13 .
The existence of a normal or metallic phase that has neither superfluid nor charge density wave order was claimed in one 14 and two 15 dimensional bosonic models in the high density limit. While normal phases at zero temperature have not been ruled out rigorously, it has been argued that they should not exist 16 . We will determine whether normal or supersolid phases exist in the one dimensional case.
Effects of disorder and impurities can be especially strong in one-dimension. If the effective interactions are repulsive, a single impurity can make some regions of the superfluid phase insulating 17 , while an even larger area of the superfluid phase is turned into a lattice glass by any finite quenched disorder 18,1 . The regions in which impurities and disorder become relevant are known 18,1,17 , and can be determined in the phase diagram without actually adding impurities or disorder to the model. Since disorder and impurities are important in experiments, we will determine these regions in the phase diagram. We will also compare the conductivity and Drude weight for the attractive and repulsive regions of the superfluid phase in pure systems without impurities or disorder.
The outline of this paper is as follows: in Section II the basic phase diagram and the possible phase transitions of the Bose-Hubbard model are discussed. Some aspects of the density-matrix renormalization group are discussed in Section. III. The calculation of the phase boundaries is presented in Section IV. The correlation functions in the different phases are shown in Section V. In Section VI the phase diagram with on-site interaction is presented. The possible existence of normal or supersolid phases is discussed in Section VII, and the phase diagram with nearest-neighbor interaction is determined. In Section VIII the ac-conductivity and the superfluid stiffness is calculated for the Mott-insulator and different regions of the superfluid phase. Conclusions are given in Section IX.

II. THE BOSE-HUBBARD MODEL
The basic physics of interacting bosons on a lattice is contained in the Bose-Hubbard model 1 . We use an extended version which includes nearest-neighbor repulsion: where the b i are the annihilation operators of bosons on site i, n i = b † i b i is the number of particles on site i, and t is the hopping matrix element. U is the on-site Coulomb repulsion and V is the nearest-neighbor repulsion. The energy scale is set by choosing U = 1.
For small interactions or large t the bosons are completely delocalized, the system is in a superfluid phase. If the density is commensurate with the lattice, and there is an interaction with the corresponding wavevector, the bosons become localized at small t. In the presence of on-site interaction only (V = 0), Mott-insulating regions with integer density are found. Fig. 1 is a sketch of the phase diagram showing Mott-insulating regions surrounded by the superfluid phase.
On most of the phase boundaries between the insulating phases and the superfluid phase the density of the system changes as the phase boundary is crossed from the incompressible insulator to the compressible superfluid. The location of this commensurate to incommensurate density transition can be directly determined as the energy it takes to add a particle or hole to the insulator: Here E 0 is the energy of the insulator groundstate, E p is the energy of a state with the density of the groundstate and an additional particle and E h of that with an additional hole. These energies can be calculated using DMRG, which will be discussed in Section IV. Note that the chemical potentials µ p c and µ h c are not equal to each other, the compressibility of the insulator is κ = ∂ρ/∂µ is zero. The superfluid phase is compressible, and for states in the superfluid phase µ p c = µ h c . The values of µ p c and µ h c at t = 0 shown in Fig. 1 can be easily calculated analytically.
At the phase transitions from the insulator to the superfluid phase where the density remains an integer, the model is in the universality class of the xy-model, and there is a Kosterlitz-Thouless 19-22 phase transition. This transition is purely driven by phase fluctuations that are determined by t. The particle-hole excitation gap at the Kosterlitz-Thouless transition closes as: giving the insulating regions a very pointed shape. The commensurate-incommensurate phase boundaries can be determined directly by calculating the particle and hole excitation energies (Eq. (2) and Eq. (3)), and in principle the Kosterlitz-Thouless transition could also be found by locating the t at which E g is zero. But since the energy gap closes very slowly (Eq. (4)), small errors in the energies lead to a big error in the location t c of the critical point. Instead, we will study the decay of the correlation functions to find the critical point. The superfluid phase of interacting bosons in one dimension has a linear dispersion relation for small wavevectors q and no excitation gap. The low energy physics of this phase is that of a Luttinger Liquid 1,23,24 with the basic Hamiltonian Here Π(x) are density fluctuations and Φ(x) phase fluc- , v is the second sound velocity and K determines the decay of the correlation for r ≫ ρ. The interactions and the lattice introduce an extra term: where ∂ x Θ(x) = π[ρ 0 + Π(x)], ρ is the density of the system, ρ 0 is the density of the insulator (e.g. ρ 0 = 1 for the Mott insulator), and n is the denominator of the density of the insulator: ρ 0 = m/n. For K > K c this term becomes relevant and drives the system into an insulating phase. At the Kosterlitz-Thouless transition K c = n 2 2 and K c = n 2 at the commensurateincommensurate transition [25][26][27] .
At the Mott-insulators with integer densities, the denominator of the density is n = 1. At the sides of the insulator the Luttinger-Liquid parameter is K = 1, and at the Kosterlitz-Thouless transition at the tip it is K = 1/2. The parameter K can be determined from the correlation functions, and we will locate the t c of the Kosterlitz-Thouless transition by finding the t at which K = 1/2.
If additional nearest-neighbor interactions are included in the model, charge density wave phases are found at half integer densities. They also have a Kosterlitz-Thouless transition at the tips, giving them a similar shape as the Mott-insulators (Fig. 2). Since the denominator of their density is n = 2, the parameter K = 4 at the sides of the phase boundary of the charge density wave, and K = 2 at the Kosterlitz-Thouless transition at the tip.
The possible existence of an intermediate phase, supersolid or normal, between the charge density wave phase and the superfluid phase will be addressed in Section VII.

III. DMRG
To determine the energies and correlation functions we use the Density-Matrix Renormalization Group (DMRG) 28,29 , a numerical method capable of delivering precise results for groundstate properties of low dimensional strongly interacting system. We use the finite-size version of the DMRG algorithm, in which the system is built up to a certain size, and the basis of the system is then optimized to represent the chosen target states by sweeping through the system repeatedly until the basis is converged.
The density matrix weight of the states discarded in a DMRG step is a measure of the numerical errors caused by the truncation. We found this truncation error to depend on the correlation length in the system. At a fixed number of states kept, we find very small truncation errors in the insulating phases, that grow as the phase transition to the superfluid phase is approached, and are biggest in the superfluid phase. Note that in one dimension the whole superfluid phase is critical with a diverging correlation length, but the correlation length is always finite in finite systems.
In each DMRG calculation for a given set of model parameters we first use the groundstate, the state with an additional particle and the state with an additional hole as target states. To obtain adequate numerical accuracy in all cases, we require the density matrix weight of the truncated states ∆ < 5 × 10 −6 (see Appendix B). The energies of these states are used to calculate the chemical potentials (Eq. (2) and Eq. (3)).
For further sweeps only the groundstate is used as a target state. We require the weight of the truncated states ∆ < 10 −9 , and the number of states kept is increased if necessary. After the basis is converged, which usually takes two sweeps, the groundstate correlation functions are calculated.
At the same parameters and number of states, the truncation error in a system with periodic boundary con- With this additional term, a particle on the boundary on average has the same potential energy as in the rest of the system.

IV. PHASE BOUNDARIES
As pointed out in Section II, the phase boundaries with the exception of the Kosterlitz-Thouless transition can be determined by calculating the particle and hole excitation energies (Eq. (2) and (3)). Using DMRG we calculate these energies in finite systems. We expect quadratic system size dependence of the energies of the insulator groundstates, and linear system size dependence for the states with additional particles and holes. Fig. 3(a) and Fig. 4 show that the leading term in the scaling of µ p c and µ h c in the insulator and the superfluid phases is 1/L. Fig.  3(b) shows µ p c of the Mott-insulator at ρ = 1 with V = 0 and t = 0.1. In this case the system size dependence is very weak, and on the scale of Fig. 3(b) the quadratic part of the scaling can be seen. Since the quadratic and higher parts contribute only very weakly to the scaling, we ignore them and use linear extrapolations from the finite system sizes to determine µ p c and µ h c in the thermodynamic limit. In the insulator phase µ p c = µ h c , since  there is a finite gap E g = µ p c − µ h c (Eq. (4)). In the superfluid phase the extrapolations for µ p c and µ h c should result in the same µ since E g = 0 and the system is compressible. In Fig. 4 small deviations from this can be seen, which we will ignore in our analysis.

A. The local density
Since open boundary conditions are used, special care has to be taken to reduce boundary effects. The most obvious form of these are local density oscillations. In the superfluid phase they show the power-law decay away from the edge of the system characteristic for the Luttinger Liquid. If the density of the bosons is given as a rational number ρ = n/m, the wavelength of the oscillations is the denominator m of the density -the same wavelength as in the density-density correlation functions (Eq. (7)). Fig. 5 shows the local density in a system at density ρ = 3/4 with only on-site interactions, and ρ = 1/4 with additional nearest-neighbor interactions, both in the superfluid phase.
In the insulators the boundary effects decay exponentially. In addition to the boundary induced density oscillations, in the charge density wave phase at density ρ = 1/2 ( Fig. 6) there is real long range order in n i in the form of a charge density wave with n i = ρ + S π (−1) i , where S π is the structure factor. In an infinite system the long range order is due to spontaneous symmetry breaking -in finite systems with even numbers of sites we allow this by adding a small symmetry breaking term to the chemical potential on the left boundary. Without this symmetry breaking term reflection symmetry would cause the groundstate to be the linear combination of two charge density wave phases with a phase difference of π, canceling out the long ranged oscillations in n i . Breaking the symmetry in the finite system reduces the Hilbert space necessary to represent the groundstate by half and also leads to a better convergence of the DMRG. In the superfluid phase the symmetry breaking term just modifies the density oscillations at the boundary.  b † r b0 correlation function in superfluid systems with density ρ = 1/4 at t = 0.1, U = 1 and V = 0.4 and different system sizes L. As the system size is increased the boundary effects become weaker, and the curves for different system sizes look similar over increasing regions. The inset shows the decay parameter K from power-law fits b † r b0 ∼ r −K/2 to 16 ≤ r ≤ 32 for the different system sizes (L = 32, 64, 96, 128 and 256). As the system size grows, K converges to the value for infinite systems. The dashed lines show the extrapolation to the upper and lower limit of K. The mean value of the upper and lower limit is taken as K, and the difference gives an estimate of the error. In the shown case it is K = 1.40 ± 0.03.

B. The hopping correlation function
In the superfluid phase the hopping correlation function Γ(r) = b † r b 0 decays with the power-law behavior given in Eq. (6), which can be used to determine the Luttinger Liquid parameter K. As discussed above, there are local density fluctuations n i in the finite systems. The creation operator can be represented by a density n i and a phase φ i part: b † i = n i exp (iφ i ). As discussed above, there are local density fluctuations n i in the finite systems, and they will affect the correlation function b † i b j . Since the local density oscillations, with the exception of the charge density wave phase, are boundary induced, and we are interested in the properties of an infinite system, we reduce the effect of the local density oscillations averaging over pairs of To minimize boundary effects we place i and j symmetrically around the center. Fig. 7 shows the power-law behavior in the superfluid phase for small r, which is modified to a faster decay closer to the system boundaries. The bigger the systems are, the bigger is the region in which the correlation functions fit the algebraic decay, and also the region in which the correlation functions look the same for different system sizes. In all cases the two biggest systems we calculate are at least 128 and 256 sites long. To estimate K, we fit a·r −K/2 to the numerical data for 16 ≤ r ≤ 32. For systems with 128 and 256 sites the boundary effects in this region are small, while the distance is also big enough to avoid short-ranged (non-Luttinger Liquid) effects.
With increasing system size the boundary effects get weaker, resulting in a decreasing K that asymptotically approaches the infinite size value. To find a simple estimate of this we use the K determined in the biggest system as an upper limit K u , and the linear extrapolation from the values in the two biggest systems as a lower limit K l . We take the mean value K = (K u + K l )/2 and estimate the error as ∆K = (K u − K l )/2.

C. The density-density correlation function
The density-density correlation functions are calculated in the same way as the hopping correlation function. However, in this case it is necessary to subtract the static expectation values, measuring n i n j − n i n j , instead of just taking n i n j . Fig. 8 shows the density-density correlation function at density ρ = 1/4. A fit with Eq. (7) works fairly well, but the first term 2 K (2πρr) −2 could not be observed. Instead of the correlation function being bigger for small r, we find it to be smaller. Eq. (7) only necessarily holds at large distances, and the short range behavior we see is dominated by the repulsive interaction between the particles. In fitting A(ρr) −2/K cos 2πρr to the data, a cut-off at small distances has to be made. While this works well enough to confirm that the correlation functions decay with a power-law behavior, the uncertainties in the fit are too high to determine K.

VI. ON-SITE INTERACTIONS
Using the methods described above we determine the phase diagram in the presence of on-site interactions only. Fig. 9 shows the Mott-insulator with density ρ = 1, surrounded by the superfluid phase. In Fig. 10 the tip of the insulator is shown on an expanded scale. The very pointed tip of the insulating region reflects the closing of the energy gap given by Eq. (4). Fig. 9 also shows results for the commensurate-incommensurate phase transition from twelfth order perturbation theory. The excellent agreement with DMRG confirms the high accuracy

achieved.
To find the critical point of the Kosterlitz-Thouless transition at the tip of the Mott-insulator, we determine the t at which K = 1/2. As described in Section V B, the K are determined by fitting power-law behavior to the decay of the hopping correlation function Γ(r) = b † r b 0 . Due to finite-size effects, the result can depend on the interval of r that is used for this fit. While we found 16 ≤ r ≤ 32 to be reasonable in general, there are logarithmic finite-size effects at the Kosterlitz-Thouless transition that require a more detailed inspection. Fig. 11 shows K's that were determined with different fitting intervals plotted versus t.
For the t = 0.29, t = 0.3 and t = 0.31 systems with L = 512 and L = 1024 sites were calculated to keep finite-size effects small. Table I shows the t c found with different fitting intervals. Due to the finite-size effects the t c go to higher t as the fitting interval is shifted to bigger r. While using a fitting interval of 16 ≤ r ≤ 32 is a good compromise between avoiding finite-size effects by using small r, and finding an asymptotic value for the decay by using big r, the error determined by this fit alone underestimates the true error. With an error estimated from the effect of the choice of the fitting interval, we  34 . In a combination of exactdiagonalization for system with up to L = 12 sites and renormalization group Kashurnikov and Svistunov found t c = 0.304 ± 0.002 35 , and together with Kravasin found t c = 0.300 ± 0.005 32 in a quantum Monte Carlo study. An exact diagonalization approach reported the critical point to be at t c = 0.275 ± 0.005 36 , and a study using 12th order strong coupling expansions 30 located it at t c = 0.26 ± 0.01.
There were two previous studies using DMRG to locate the Kosterlitz-Thouless transition in the Bose-Hubbard model. In the first study 37 a particle cut-off of n = 4 per site was used, and the critical point was determined by using a phase twist φ = π and periodic boundaries to calculate the superfluid stiffness. The critical point was found at t c = 0.298. In another DMRG study 33 , the critical point was found at t c = 0.277 ± 0.01 by determining the t at which K = 1/2, similar to the procedure used in the present work. In both of these DMRG studies the numerical accuracy was limited due to the use of the infinite size version of the DMRG algorithm, periodic boundary conditions and a small number of states in the DMRG basis. The biggest system sizes were smaller than L = 80.
The range of these results demonstrates the difficulty involved in determining the critical point of the Kosterlitz-Thouless transition, which is mostly due to logarithmic finite-size effects close to the critical point. The large system sizes used in this paper should compensate for this within the given error bars, and yield a reliable result.
The phase diagram shown in Fig. 9 has a very interesting feature. Imagine moving on a line of constant chemical potential µ, for example µ = 0.15, and starting at small t, moving toward bigger t. The particle density along this line is illustrated in Fig. 12. For small t the system is in the Mott-insulator phase. At t ≈ 0.1 there is a phase transition to the superfluid phase, with densities ρ < 1. The density decreases up to a minimum, then it start increasing again. At t ≈ 0.26 the density goes up to ρ = 1 again, and there is another phase transition, this time reentering the Mott-insulating phase from the superfluid phase. Increasing t further leads to another phase transition from the Mott-insulator to the superfluid phase, this time with ρ > 1, and the density increasing further with increasing t.
To gain more insight into this, we compare to results from a mean field approach 38 . Fig. 13 shows the phase diagram with the Mott-insulator at density ρ = 1, surrounded by the superfluid phase. In the superfluid phase, the lines of constant density slope downward as t is increased. This is not only found in one dimension, but in all dimensions. The limit of t → ∞ corresponds to keeping t constant and setting the interactions to zero. If the interactions are zero the system goes from a superfluid phase to a Bose-Einstein condensate, in which every particle has an energy of −2t. If the chemical potential is smaller than −2t, the system is empty, because it costs energy to put a particle in, and for chemical potentials bigger than −2t the number of particles goes to infinity, because every additional particle reduces the total energy of the system. Going back to the picture of constant interactions and changing t, this means that the density of the system always goes to infinity as t is increased.
In dimensions two and higher the superfluid-insulator transition on the line of constant density is a second order transition, and the tip of the insulating region is round.  as t is increased, and then increases again.
In one dimension the tip of the insulating region is very long and narrow due to the Kosterlitz-Thouless transition, and it is possible to reenter into the insulator at ρ = 1.

VII. NEAREST-NEIGHBOR INTERACTION
Longer ranger interactions have been found to be important in experiments [7][8][9] . We now include nearestneighbor interactions by setting V = 0.4. Due to the nearest-neighbor interactions a new insulator phase appears at half integer densities. It is a charge density wave phase (CDW) with a wavelength of two sites, and like the Mott-insulator at integer density it has an excitation gap and is incompressible. The crystalline order is characterized by a non-zero structure factor In Fig. 6 the local density oscillations in the charge density wave phase are shown. A small boundary effect can be seen, but the main feature are long-range density oscillations throughout the system that do not decay. An order parameter n i − ρ can be defined to describe this charge density wave, even in one dimension. In one dimension the superfluid phase is signalled by a diverging correlation length and a non-zero superfluid stiffness which is proportional to the Drude weight D ∼ ρ s . In two and higher dimensions there is also an order parameter b † i = 0. In one dimension the whole superfluid phase is critical, and there is no order parameter.
At the transition from the charge density wave phase to the superfluid phase both types of order are involved: the crystalline order in the charge density wave phase and the superfluid order in the superfluid phase. In addition to a direct phase transition from the charge density wave to the superfluid at which the crystalline order vanishes at the same point where the superfluid order appears, there is the possibility of an intermediate phase. Tab. II shows the possible phases close to density ρ = 1/2 in a bosonic system with on-site and nearest-neighbor interaction in two or higher dimensions. In addition to the charge density wave and the superfluid phase, supersolids that have both forms of order were found in two dimensional models 10,11 . Baltin and Wagenblast 14 found a region that has neither superfluid stiffness nor charge density wave in a one dimensional bosonic model in the high density. The possible existence of such a phase was also recently predicted for a two dimensional bosonic model in the high density limit by Das and Doniach 15 , who call it a Bose-metal.
In Appendix C strong coupling expansions are used to illustrate the difference between the commensurateincommensurate phase transition at ρ = 1/2 in one and two dimensions. Strong coupling expansions can be used to study the insulator, but not the superfluid. To study the low-energy behavior of the superfluid phase the Luttinger liquid can be used. In addition to the basic Luttinger liquid Hamiltonian (Eq. (5)), the lattice and the interactions introduce scattering terms (Eq. (8)). These only contribute at ρ = 1/2, where they can drive the system into a different phase, but not at nearby densities. At incommensurate densities close to ρ = 1/2, the wavefunction is incommensurate with the lattice and hence cannot be pinned to the lattice to form an insulator. Of course this would be changed if there were impurities or disorder, but a pure system in one dimension is in the Luttinger liquid phase unless it is at a density commensurate with the lattice and the interactions.
At density ρ = 1/2 DMRG can be used to determine if there is an intermediate phase or a direct phase transition from the charge density wave phase to the superfluid phase. To do this, we investigate the relationship between the superfluid and crystalline order at the phase transition. The onset of superfluidity is signaled by a  diverging correlation length ξ (Eq. (11)), charge density order is meassured by S π > 0 (Eq. (10)).
The model is in the universality class of the xy-model at the phase transition on the line of constant density at ρ = 1/2. Due to the Kosterlitz-Thouless transition expected on this line, it is difficult to determine exactly when the structure factor or the inverse correlation length go to zero. But it is possible to study the dependence of the structure factor on the correlation length. Fig. 15 shows that for small values the structure factor depends on the correlation length by a power-law: To keep the effect of the boundaries small for the calculation of both the structure factor S π and the correlation length ξ, only sites that were at least a quarter of the system size away from the boundaries were taken into account. In Fig. 16 the extrapolated inverse correlation length 1/ξ(S π = 0) at zero structure factor is plotted against the inverse system size. The linear fit to the three biggest system sizes shows that 1/ξ(S π = 0) goes to zero for infinite systems. From this and the power-law behavior in Fig. 15 we conclude that there is a power-law dependence of the structure factor on the correlation length, and that in infinite systems the correlation length diverges at the same point at which the structure factor goes to zero. This means that there is a direct phase transition from the charge density wave phase to the superfluid phase, and no supersolid or normal phase in between.
The phase boundaries of the charge density wave phase can be found in the same way as those of the Mottinsulator, and we use the methods used for the on-site only interaction case to calculate the phase diagram. where systems with up to L = 512 sites are used, are much more accurate.
At the phase boundaries of the charge density wave phase the Luttinger Liquid parameter is K = 4 except for the Kosterlitz-Thouless transition at the tip, where it is K = 2. The charge density wave phase is surrounded by a region where K > 1. The effective interactions in the Luttinger-Liquid are attractive for K < 1 and repulsive for K > 1. Kane and Fisher showed that in the repulsive region a single impurity turns the system into an insulator 17,27 . An even larger region of the Luttinger liquid is driven into a glass phase if K > 2/3 for any finite quenched disorder 18,1 .
In the phase diagram these regions (K > 1 and K > 2/3) are determined by doing calculations for different densities and t. For example, to determine the t for which K(t) = 1 at a given density ρ, we do calculations for different t until we find a pair of t 1 and t 2 that are close to each other, and K(t 1 ) > 1 > K(t 2 ). We then determine t(K = 1) by linear interpolation. To find the boundaries of the repulsive Luttinger Liquid region, we calculate t(K = 1) for various densities.
The lines with K = 1 and K = 2/3 are shown in Fig.  18. The repulsive Luttinger Liquid (K > 1) region completely surrounds the charge density wave phase. Instead of going to t = 0 as the Mott-insulator is approached, the K = 1 line ends in the side of the Mott-insulating region, where K = 1 for all of the commensurateincommensurate phase boundary. Fig. 17 illustrates how the two lines with K = 1 meet at the phase boundary of the Mott-insulator. Although we could not obtain more detailed results for densities closer to ρ = 1, we argue that the lines with K > 1 bend towards t = 0 as the density gets closer to one, while those with 1/2 < K < 1 bend towards the tip of the Mott-insulator.
Lines of constant K with 1/2 < K < 1 are discontinuous at ρ = 1, where the system is an insulator for K > 1/2, with the tipped shape reflecting the Kosterlitz-Thouless behavior. Lines with K ≤ 1/2 are round at ρ = 1, and do not reflect the Kosterlitz-Thouless behavior. Analogous to ρ = 1, we also observe that the lines of constant K are round for K ≤ 2 and ρ = 1/2, where the charge density wave phase ends in a Kosterlitz-Thouless transition at K = 2.

VIII. AC-CONDUCTIVITY
The repulsive region of the Luttinger liquid is turned into an insulator by a single impurity 17 . This raises the question if there is a qualitative difference between the conductivity in the repulsive region and the attractive region of the Luttinger liquid in the pure system. The regular part of the conductivity is given by: and the current operator is Recent developments with DMRG make the calculation of dynamical correlation functions like the acconductivity possible 39 . The conductivity at a frequency z = ω + iη can be calculated as the direct product of the current operator applied to the groundstate and a correction vector By using these two states and the groundstate |0 as DMRG target states, the conductivity σ reg 1 (z) can be calculated very precisely. To calculate the conductivity over an interval of width η ranging from ω 1 to ω 2 , we use correction vectors | x(ω 1 + iη) and | x(ω 2 + iη) as target states. At the end of the DMRG calculation, when the DMRG basis is optimized to represent these states, we calculate the conductivity from ω 1 to ω 2 . Repeating this procedure for neighboring intervals, we piece together the conductivity for a whole range of frequencies.
The finite broadening η in the correction vectors is only used to obtain appropriate DMRG target states. To calculate the spectrum within the DMRG basis, we use a Lanczos method that yields approximate eigenstates of the Hamiltonian. The broadening used in our plots is then applied to the discrete peaks found with the Lanczos method, and is only used for better visualization.
DMRG calculations work best with open boundary conditions. How the current operator is applied in a system with open boundary conditions is discussed in Appendix D.
The conductivity in the Mott-insulator phase is shown in Fig. 19. There is an energy gap of ∆ω = 0.54, with a big peak after that and only small excitations at higher TABLE III. The Drude weight D for different system sizes in the Mott-insulator with density ρ = 1, t = 0.05, U = 1,V = 0.4. T is the kinetic energy, σ reg 1 (ω)dω the integral over the ac-conductivity. Also shown is the number of states m, the broadening η of the correction vectors, and the truncation error ∆. energies. Since it is an insulator phase, we expect to find no Drude weight. With the kinetic energy defined as the Drude weight is given by Note that the Drude weight is proportional to the superfluid stiffness ρ s given in Eq. (12).
Since the kinetic energy T can be calculated directly with DMRG, and we expect D = 0, this is an opportunity to verify the consistency of the calculation. In Tab. III the Drude weight is shown for various system sizes. A small finite-size effect can be seen in the data. From the differences in the individual values we estimate the error of the Drude weight ∆D = 0.02t, or 2% of − T /L.
In the superfluid phase we find precursor peaks at small frequencies in the conductivity. They are due to the finite    20 shows these precursor peaks in the repulsive Luttinger liquid for different system sizes. As the system size is increased the precursor peaks move towards ω = 0. For the calculation of the Drude weight these peaks should not contribute, and we use low-energy cut-offs to ignore them. Fig. 21 and Fig. 22 show the conductivity in the superfluid phase in the repulsive and attractive Luttinger liquid regions. In the Luttinger liquid the conductivity was predicted to increase with a power-law for small frequencies, and decay exponentially for big frequencies 24,25 . The conductivity in the attractive region shown in Fig.  22 is in good qualitative agreement with this. In the repulsive case (Fig. 21) there are too few peaks to clearly identify this behavior. Bigger systems would have to be studied to determine if the overall shape is qualitatively different from the attractive region.  The Drude weight in the repulsive and attractive regime of the superfluid phase is shown in Tab. IV and Tab. V. For some system sizes data with different numbers of states m and broadening η is shown. The numerical accuracy depends on these parameters, with bigger m and smaller η for higher accuracy. The data in Tab. IV and Tab. V shows that the impact of m and η on the Drude weight is small. In both the attractive and repulsive case we find big non-zero values that are close to the  V. The Drude weight D for different system sizes in the attractive region of the superfluid phase at density ρ = 3/4, t = 0.5, U = 1,V = 0.4. Same notation as in Tab.III. kinetic energy per site in the systems. The differences in the Drude weight in different system sizes, with the exception of the smallest systems, are rather due to numerical errors that grow with the system size, than due to finite-size effects.

IX. CONCLUSIONS
In summary, we have studied the phase diagram of the one-dimensional Bose-Hubbard model with on-site only interactions and with additional nearest-neighbor interaction. The density matrix renormalization group (DMRG) was used to calculate chemical potentials for given densities and model parameters, and by doing this for sets of parameters the phase boundaries of the Mottinsulators and the charge density wave phase were determined.
The low-energy behavior of the superfluid phase of onedimensional bosonic systems is that of a Luttinger liquid. We determined the Luttinger liquid parameter K from the decay of the hopping correlation functions. Since the value of K is known for insulator-superfluid transitions, we could use it to locate the Kosterlitz-Thouless transitions at the tips of the Mott-insulators and the charge density wave phase in the µ-t phase diagram.
In the charge density wave phase we found that close to the phase transition the structure factor depends on a power-law of the superfluid correlation length. From this we conclude that there is a direct phase transition from the charge density wave phase to the supersolid, and no intermediate phase like a supersolid or normal phase.
The charge density wave phase is surrounded by a region of the superfluid phase where K > 1, which corresponds to a Luttinger liquid with repulsive effective interactions. Kane and Fisher have shown that this region will be turned into an insulator by a single impurity. We determined the boundary of the repulsive region by finding the line where K = 1 in the phase diagram. We found that this boundary does not go to t = 0 as the Mott-insulator is approached, but ends in the side of the Mott-insulating region, where the Luttinger liquid parameter also is K = 1.
We calculated the ac-conductivity in the Mottinsulator and the superfluid phase. In the Mott-insulator and the attractive region of the superfluid phase the ac-conductivity has the expected shape. In the repulsive region of the Luttinger liquid we found a different shape, but could not determine if this is due to the finite system sizes. The Drude weight or superfluid stiffness was found to be big in both the attractive and the repulsive region. To verify the effect of this truncation on the correlation function Γ(r) = b † r b 0 , we calculate systems with different n max . Fig. 23 shows Γ(r) in the Mott-insulator. Due to the small particle-hole excitations, the correlation functions are almost identical for n max ≥ 3. In the superfluid phase there are more particle-hole excitations which are affected by the truncation. Fig. 24 shows  that the correlation functions are independent of n max for n max ≥ 4. By choosing n max = 5 for all calculations the effect of the truncation should be small enough not to affect the results presented in this work.

APPENDIX B: TRUNCATION OF THE DMRG BASIS
In every DMRG step the basis is truncated, and only the eigenstates of the density matrix with the biggest eigenvalues are kept. The density matrix weight of the discarded states ∆ is a measure of the error caused by these truncations. To verify to which extent the truncation errors affect the results, we calculate the correlation  function Γ(r) = b † r b 0 with different numbers of states kept in the DMRG basis. Fig. 25 shows Γ(r) with different truncation errors in the Mott-insulator. Even for very small numbers of states the discarded weight is very small, and the dependence on the weight of the discarded states is weak. Note that the discrepancies are mostly apparent due to the logarithmic scale.
The correlation function in the superfluid phase is shown in Fig. 26. We find that the discarded weight with the same number of states is bigger than it is in the insulator. At small distances r the correlation functions are very similar for all numbers of states, with increasing differences as r is increased. If the discarded weights ∆ are smaller than 8 × 10 −8 , the correlation functions coin-  cide even at the boundaries of the system. By requiring the discarded weight to be smaller than ∆ ≤ 10 −9 for the calculation of the correlation functions, accuracy should be high enough in all cases.
The chemical potentials are calculated from the energies it takes to add a particle or hole. Figures 27 and 28 show chemical potentials calculated with different numbers of states kept versus the discarded weight ∆. The error bars correspond to the changes in the energies during a DMRG sweep. Differences in the chemical potentials are small for ∆ < 10 −5 . We require the discarded weight to be smaller than ∆ < 5 × 10 −6 for the calculations of the chemical potentials, and to improve the results further, we extrapolate linearly from the two values with the lowest ∆.

APPENDIX C: STRONG COUPLING EXPANSION AT THE COMMENSURATE-INCOMMENSURATE TRANSITION
The fundamental difference between the commensurate-incommensurate phase transition at ρ = 1/2 in one and two dimensions can be illustrated with the help of a strong coupling expansion. In the strong coupling limit the kinetic energy is zero. The zero order states are the groundstates of the Hamiltonian only including the particle-particle repulsion: The series expansion is made in terms of the kinetic energy term: Strong coupling expansions of this type have been successfully used to study the phase diagram with on-site only interaction 30,40,41 . To determine the phase boundaries of the Mott-insulator at ρ = 1, first the groundstate of Eq. (C1) has to be found. In this state there is simply one boson sitting on every site. Higher terms of the perturbation series introduce local particle-hole excitations. The chemical potentials on the boundaries can be determined from the energy it costs to add a particle (Eq. (2)) or a hole (Eq. (3)). But in these cases, the zeroth order ground state is degenerate, since the additional particle or hole can sit on any site. This degeneracy is lifted in first order perturbation theory. In first order the problem is reduced to the additional single particle moving on a uniform background of completely localized particles. Since the extra particle gains energy by hopping from site to site, it becomes completely delocalized. Although this behavior can be modified in higher order of the perturbation series, and there are limitations due to the radius of convergence, it is interesting to note the difference between the perturbation series in the insulator and with an additional particle or hole. At integer density the series starts out with a completely localized state, while it starts with a completely delocalized state if there is an additional particle or hole. This is in good agreement with the fact that there is a Mott-insulator for small t at integer density, and a direct phase transition to a superfluid (delocalized) phase if the density is changed.
A similar strong coupling expansion can also be used at the charge density wave phase at ρ = 1/2. The charge density wave phase at t = 0 is a state with alternating particle numbers, in one as well as two dimensions ( Fig.29 and Fig.30). Higher order terms in the perturbation series introduce local particle hopping without destroying the charge density wave order. At the commensurate-incommensurate transition, an additional particle (or hole) enters the system. For V < U/2 the energy is smallest if the additional particle goes to one of the empty sites. Fig. 31 shows how the additional particle fits into the two-dimensional charge density wave. In higher orders of perturbation theory the additional particle, as well as particle-hole excitations, hop on the charge-density background without destroying it. From this we cannot infer if the true (non-perturbation theory) groundstate is superfluid or not, and if the charge density wave order is destroyed by the particle hopping. Nevertheless, it is interesting to note that for this case supersolids have been found 10,11 in two dimensions. Close to the charge density wave phase at ρ = 1/2, the charge density order survives at small doping.
In contrast to this, the one-dimensional case looks quite different. The additional particle also goes to an unoccupied site. If the charge density wave remains unchanged, the additional energy is ∆E = 2V . With the structure factor S π , the charge density wave is given by an order parameter n l = ρ + S π exp (iπl + iφ 0 ). An additional particle or hole can also be added by shifting the phase φ 0 by π over a region with an odd number of sites. Fig.32 shows an example of such a state. In the center there is a domain with a π phase shift, and the number of particles compared to the charge density wave (Fig.29) is increased by one. And the additional energy is also ∆E = 2V . To lift the degeneracy between all these states in first order perturbation theory, it can be seen that a particle hopping next to the domain boundary is equivalent to the domain wall moving. Since energy is gained by this, the domain walls are completely delocalized in first order perturbation theory. Unlike the two dimensional case, where the charge density wave order survives in all orders, in one dimension it is destroyed in first order. While a perturbation series does not necessarily converge, this striking feature illustrates the fundamental difference between the one and two dimensional case.

APPENDIX D: CURRENT OPERATOR WITH OPEN BOUNDARIES
DMRG calculations work best with open boundary conditions. To calculate the conductivity with DMRG, the current operator has to be implemented. The current operator as it is given in Eq. (15) can be used directly with periodic boundary conditions, but to apply the current operator with open boundary conditions we modify it with a filter function: The filter function P (x n /M ) used here is a Parzen filter, x n is the distance of site n from the middle of the system, and M = L/2 is half the system size. The Parzen filter looks very similar to a Gauss function, but goes smoothly to zero at the system boundaries. It is given as: A prefactor a is chosen to provide results with the same amplitude as those found in systems with periodic boundary conditions, with a chosen so that P (x n /M ) 2 = 1. and one with open boundaries, for otherwise identical system parameters. Fig. 33 and Fig. 34 show the conductivity in the Mott-insulator and the superfluid phase with open and with periodic boundary conditions. Even in the small systems with L = 32 the curves are quite similar. In the superfluid phase a precursor peak at small frequencies can be seen in the system with open boundary conditions. Fig. 20 shows how the precursor peaks move to smaller frequencies as the system size is increased. They are an artifact of the open boundary conditions, and we use frequency cut-offs to exclude them from the calculation of the Drude weight. Tab. VI shows the Drude weights in the insulator and the superfluid. The values for open and periodic boundary conditions compare quite well, and we estimate an error of ∆D/t ≈ 0.02.