Superfluid to Mott-insulator transition of hardcore bosons in a superlattice

# Superfluid to Mott-insulator transition of hardcore bosons in a superlattice

Itay Hen    Marcos Rigol Department of Physics, Georgetown University, Washington, DC 20057, USA
July 4, 2019
###### Abstract

We study the superfluid to Mott-insulator (SF-MI) transition of hardcore bosons in commensurate superlattices in two and three dimensions. We focus on the special case where the superlattice has period two and the system is at half filling. We obtain numerical results by using the stochastic series expansion (SSE) algorithm, and compute various properties of the system, such as the ground-state energy, the density of bosons in the zero-momentum mode, the superfluid density, and the compressibility. We employ finite-size scaling to extrapolate the thermodynamic limit, and find the critical points of the phase transition. We also explore the extent to which several approximate solutions such as mean-field theory, with and without spin-wave corrections, can help one gain analytical insight into the behavior of the system in the vicinity of the phase transition.

superfluidity, Mott-insulator, hardcore bosons
###### pacs:
64.70.Tg, 03.75.Lm, 02.70.Ss, 67.85.-d

## I Introduction

Recent developments in the field of ultracold Bose gases have opened a new promising avenue of theoretical and experimental research in the study of the phases of quantum matter. A gas of bosonic atoms in an optical trap has been reversibly tuned from a Bose-Einstein condensate to a state composed of localized atoms as the strength of a periodic optical potential was varied.gr1 (); gr2 () This is an example of a quantum phase transition; a phase transition generated by quantum fluctuations and correlations rather than by a competition between the energy of a system and the entropy of its thermal fluctuations.sachdev () Understanding this phenomenon has emerged as one of the most challenging and interesting tasks of condensed matter physics. Theoretically, it is generally accepted that it can be studied using the Bose Hubbard model, where the transition is thought to be from a superfluid phase to a Mott-insulator (SF-MI), as examined in the seminal paper by Fisher et al.,r4 () with an application to He absorbed in porous media in mind. The relevance of the Bose-Hubbard model to gases of alkali-metal atoms in optical lattices was realized in Ref. r5, , and recent developments have been reviewed in Refs. r6, and bloch08, .

Interestingly, the Bose-Hubbard model is nonintegrable even in one dimension (as opposed to, say, its fermionic counterpartfhb ()). Gaining analytical insight into the SF-MI phase transition thus normally requires resorting to numerical and variational methods such as strong-coupling expansion,sce1 (); sce2 () coarse graining,cg () mean-field theories,mft () field-theoretical approachesftapp () or other perturbative methodsper1 () for a better understanding of this phenomenon. Within the variational approach, the phase transition is taken to be the point at which the variational ansatz has lower energy than a delocalized Bogoliubov state (where a fixed particle number at each lattice site is constrained).

In a recent paper, Aizenman et al.Aizenman () considered an alternative model for the study of the SF-MI phase transition. They studied the half-filled Bose-Hubbard model in the limit of infinite on-site repulsion (i.e., the case of hardcore bosons), in the presence of an alternating on-site chemical potential (a superlattice with period two). They showed that this model exhibits all the salient properties apparent in the Bose-Hubbard model, while also being more ‘treatable’ analytically. Specifically, they were able to rigorously prove the existence of superfluid and Mott-insulating phases in three dimensions. In addition, it is also known that this very same model is exactly solvable in one-dimension through a mapping to noninteracting fermions. In this case, the half-filled system is insulating for any nonzero alternating potential.Val () The off-diagonal one-particle correlations and the momentum distribution function of this model, as well as its nonequilibrium dynamics, were computed by exact meansrigol04 () in Ref. rigol06, .

Motivated by the aforementioned results, here we study the SF-MI phase transition of hardcore bosons in the presence of an alternating potential in two and three dimensions. We focus on the case where the system is at half-filling, in which case the transition between the superfluid state and the insulating state occurs at fixed density. Our first goal is to accurately determine the critical values of the alternating potential strength at which the phase transition takes place. As a next step, we analyze the results of different approximate solutions, such as mean-field theory with and without the addition of spin-wave corrections, as these allow for an analytical treatment of the problem.

Our approach is to first perform high-precision numerical simulations using the stochastic series expansion (SSE) algorithm SSE1 (); SSE2 () in order to find the critical points of the superfluid to Mott-insulator phase transition in the various dimensions. The quantities we calculate are the free energy , the density of bosons in the zero-momentum mode , 111For homogeneous systems, the number of bosons in the zero-momentum mode coincides with the condensate occupation. However, this need not be the case in general, as the condensate occupation is defined as the largest eigenvalue of the one-particle density matrix. The latter quantity will in general be different from the occupation of the zero-momentum state if the system is inhomogeneous.rigol06 () the superfluid density and the compressibility . The latter three quantities signify the transition from a superfluid to an insulator by dropping to zero at this point (while having nonzero values in the superfluid regime). We then employ mean-field and spin-wave analyses, which allow for some analytical insight into the behavior of our observables of interest and the location of the critical point. Our use of these approximation methods is partly motivated by results previously reported by Bernardet et al.,hcb () who studied the homogeneous version of the model in two dimensions. There, the mean-field approximation alone was shown to provide a fairly good qualitative description of the model, and remarkably enough, when spin-wave corrections were added, quantities such as the superfluid density and the condensate fraction were found to be virtually indistinguishable from their exact-numerical counterparts.

The paper is organized as follows. In Sec. II we briefly review the model at hand. In Sec. III, we present the exact numerical solutions obtained using the stochastic series expansion (SSE) algorithm. We compute the various physical quantities at zero temperature, and find the critical values of the SF-MI phase transition. In Sec. IV we proceed to study several approximation schemes, namely mean-field approaches and spin-wave corrections, comparing the critical values obtained using these methods, with the SSE results. In Sec. V we conclude with a few comments.

## Ii The model

The Hamiltonian for hardcore bosons in a period-two superlattice in -dimensions, with sites and periodic boundary conditions, can be written as:

 ^H=−t∑⟨ij⟩(^a†i^aj+^a†j^ai)−A∑i(−1)σ(i)^ni−μ∑i^ni.

Here, denotes nearest neighbors, () destroys (creates) a hardcore boson on site , is the local density operator, is the global chemical potential, and is an alternating local potential with on the even sublattice and on the odd sublattice. The hopping parameter sets the energy scale.

The hardcore boson creation and annihilation operators satisfy the constraints

 ^a†2i=^a2i=0, {^ai,^a†i}=1, (2)

which prohibit double or higher occupancy of lattice sites, as dictated by the limit of the Bose-Hubbard model. For any two different sites , the creation and annihilation operators obey the usual bosonic relations

 [^ai,^aj]=[^a†i,^a†j]=[^ai,^a†j]=0. (3)

The expected phase diagram of the model, in dimensions higher than one, is sketched in Fig. 1. Our model has two (trivial) insulating regimes corresponding to a completely filled lattice (with particle density ), obtained for large and positive chemical potential values, and a second insulating regime which corresponds to an empty lattice, formed in the case where the chemical potential is large and negative. These two regimes are also present in the absence of the alternating potential. The alternating one-body potential creates another insulating phase, one for which the density of particles is . In this case, the alternating potential, will in some cases (depending on its strength) create a gap in the energy spectrum, generating a superfluid to Mott-insulator transition. As the latter regime is the one which is of interest to us, we shall henceforth set the global chemical potential to . In this case, the model has particle-hole symmetry which in turn fixes the density at as desired.

Before moving on, a remark is in order. The insulating phase of the model at hand is a consequence of a counterbalance between strong on-site interactions (which in our model are in fact infinite) and an alternating potential. The resulting local density will thus be different on the odd sublattice than on the even sublattice. While this state is sometimes referred to as a charge density wave,Val () in what follows, we shall address this phase as a Mott-insulator, in the spirit of Ref. Aizenman, .

## Iii Numerical results

We obtain numerically-exact results for the model at hand by performing numerical simulations based on the stochastic series expansion (SSE) algorithm.SSE1 (); SSE2 () As our main objective is to find the critical points of the SF-MI phase transition in the various dimensions, simulations are performed for a range of values of the ratio and for various system sizes. Since we are interested in the zero-temperature properties of the system, simulations are performed with high inverse-temperature (in our units, ), where in most cases we will find it sufficient to have in order to obtain virtually zero-temperature results. (The effects of increasing beyond this value are indiscernible.)

Finite size effects are corrected by repeating the simulations with different system sizes. The thermodynamic-limit value of the phase transition is then extrapolated by performing finite size scaling of the results in the vicinity of the phase transition: Around the critical point, most physical quantities (which we denote here by ) scale according to the general rule:

 XLξ/ν=F(|A−Ac|L1/ν), (4)

where is a universal scaling function, is the shifted control parameter ( being the control parameter, and – the critical value), is the correlation length critical exponent and is the critical exponent belonging to the observable . The values of these exponents are determined by the universality class the transition belongs to. In our case (and in the Bose-Hubbard model for integer filling as well), it is the () dimensional universality class.r4 (); Bruder () We note that the above universality class characterizes only the fixed-density transition (the dashed line in Fig. 1). The transition driven by changing the density belongs to the mean-field universality class and is characterized by different critical exponents.

Equation (4) above will help us find the critical point, as it tells us that (a) the quantity should be independent of the system-size at the phase transition, and (b) when plotting against the resulting curve should be independent of the system-size as well.

The quantity we shall be using to that end is the superfluid density, which has the critical exponent (see Ref. r4, for details) where is the dimension, and is the dynamical critical exponent, which in our case is .Bruder () The correlation length exponent is dimension-dependent and takes the values , and in one, two and three dimensions, respectively.

### iii.1 One dimension

In one dimension, our model has an analytic solution.Val () This is due to the Jordan-Wigner transformation which enables the mapping of the hardcore bosons Hamiltonian to that of noninteracting spinless fermions.Val () The latter Hamiltonian may be diagonalized to produce exact analytical results. In this case, the SF-MI phase transition occurs at , i.e., the system is superfluid only when the alternating potential is absent, in which case it exhibits off-diagonal quasi-long-range order (a power-law decay of the one-particle correlations). In that sense, one may say that the system exhibits quasi-condensation when .Val (); rigol04 (); rigol06 ()

Simulations in one dimension were thus performed only as a check on our computational method. No discrepancies between the analytical solution and the numerical one were found: In Fig. 2, the superfluid density is plotted against for different system sizes (here, ). In the figure, all curves intersect at the critical point , indicating the location of the phase-transition in the thermodynamic limit, in agreement with the analytic results. The inset shows the scaled superfluid density as a function of the scaled control parameter, in which case all curves should be, and in fact are, on top of each other. The numerical value for the superfluid density at the transition coincides with the expected value of given by the analytic solution.Val ()

As superlattices such as the one we study here have already been realized in experiments with ultracold bosons in optical lattices,peil03 (); strabley06 (); strabley07 (); lee07 () and the observable usually measured in those kind of experiments is the momentum distribution function , we plot it in Fig. 3 for two different values of . Due to the quasi-long-range decay of one-particle correlations in the superfluid phase, the momentum distribution function has a peak at [Fig. 3(a)]. On the other hand, in the insulating phase, the one-particle correlations decay exponentially, yielding a broad momentum distribution [Fig. 3(b)]. This leads to the following observation: As one increases the size of the lattice (while keeping the density fixed), one finds that in the superfluid phase increases for small values of [Fig. 3(a)], while for the insulating phase [Fig. 3(b)] this does not happen. Exact results for (using the approach described in Ref. rigol04, ), are also presented in Fig. 3. As expected, the SSE results are right on top of the exact ones.

### iii.2 Two dimensions

In dimensions higher than one, no analytic solution to the model exists, so accurate results are obtainable only numerically. Here, we have applied the SSE algorithm to systems of sizes ranging from to , with inverse-temperature . In Fig. 4, the scaled superfluid density is plotted against for the different system sizes (the errors are on the order of magnitude of the symbol sizes). All curves intersect at , signifying the phase transition. The inset shows the scaled superfluid density as a function of the scaled control parameter. As in the one-dimensional case, all data points fall into a single curve. The value for the critical point we obtained here agrees with the value recently obtained by Priyadarshee et al.2dcrit ()

The momentum distribution function in the superfluid and insulating regimes are shown in Figs. 5(a) and 5(b), respectively. In two dimensions, the superfluid regime exhibits true off-diagonal long-range order, which means that the peak is sharper that in one dimension, which exhibits only quasi-long-range order. This can be seen in Fig. 5(a). The Mott-insulating regime is once again characterized by an exponential decay of one-particle correlations. The corresponding momentum distribution function has a broad peak around as shown in Fig. 5(b).

### iii.3 Three dimensions

In three dimensions, we have performed simulations for system sizes ranging from to and an inverse temperature of . Here, the SF-MI phase transition is found at , as indicated by the scaled superfluid density plotted as a function of in Fig. 6, for the different system sizes. The inset in Fig. 6 depicts the scaled superfluid density as a function of the scaled control parameter, exhibiting the collapse of all data points into a single curve, as in one and two dimensions. The momentum distribution function in three dimensions is qualitatively similar to its two-dimensional counterpart, both in the superfluid phase and in the insulating phase, and thus will not be presented here.

## Iv Approximation schemes

Having obtained the critical values via quantum Monte Carlo techniques, we now turn to look for approximation schemes that would provide analytical descriptions of the phase transition. We start this investigation with the Gutzwiller mean-field approach. Before doing so however, we recall that the model at hand can also be viewed as the model of a spin-1/2 system.mats () We shall make use of this correspondence, utilizing the exact mapping between bosonic operators and generators, namely,

 ^a†i ↔ S+i, (5) ^ai ↔ S−i, ^a†i^ai ↔ Szi+1/2.

With this mapping, the hardcore bosons Hamiltonian, Eq. (II), becomes that of the antiferromagnet with an alternating magnetic field applied along the direction:

 ^H= − t∑⟨ij⟩(S+iS−j+S+jS−i) (6) − ∑i[μ+A(−1)σ(i)](Szi+12).

### iv.1 Mean-field approach

We start our mean-field calculation with the following product state as an initial ansatz:

 |0⟩{\scriptsize MF}=⊗∏j(sinθj2|↓⟩+cosθj2eiφj|↑⟩), (7)

where specify the orientation of the -th spin. Obviously, we expect every other site to be described by the same wave function, due to the symmetry of the problem. This is schematically shown in Fig. 7.

As we are using the grand-canonical scheme, the orientations of the spins will be determined by minimizing the grand-canonical potential (per site)

 Ω{\scriptsize MF} = {\scriptsize MF}⟨0|^H|0⟩{\scriptsize MF}=−t2N∑⟨ij⟩sinθisinθjcos(ϕi−ϕj) (8) −12N∑i[μ+A(−1)σ(i)](1+cosθi).

with respect to these angles. For the azimuthal angles, this simply implies a constant (yet arbitrary) value , while for the polar angles, extremization yields

 cosθ1 =μ1√1+μ221+μ12, (9a) cosθ2 =μ2√1+μ121+μ22, (9b)

where . These values correspond to a minimal configuration only in the region . Outside this region, the system is saturated, and the solution is one where all spins are aligned – pointing either all up or all down. In bosonic language, these latter configurations correspond to the completely full/empty insulating configurations.

At this point we can calculate the following quantities. First, the density of particles is:

 ρ{\scriptsize MF} = 1N∑i{\scriptsize MF}⟨0|^a†i^ai|0⟩{\scriptsize MF}=12+12N∑icosθi (10) = 12+14(cosθ1+cosθ2).

Next, the free energy becomes

 Ω{\scriptsize MF} = {\scriptsize MF}⟨0|^H|0⟩{\scriptsize MF}=−dt2sinθ1sinθ2−μ2 (11) − 14(μ+A)cosθ1−14(μ−A)cosθ2,

and the density of bosons in the zero-momentum mode is calculated as:

 ρ0,{\scriptsize MF} = 1N{\scriptsize MF}⟨0|^a†k=0^ak=0|0⟩{\scriptsize MF} = 14N2∑i,jsinθisinθj=116(sinθ1+sinθ2)2.

The superfluid density requires a special treatment of the boundary conditions. As is well known,rhoS () one can relate the superfluid density to the “spin stiffness”. To accomplish this, one needs to compare (the free energy) of the system under periodic conditions with the free energy under a “twist” in the boundary conditions along one of the linear directions (say, the direction). In the periodic case, which we treated above, the azimuthal angles were all identical. To implement a twist, we take this angle to be site-dependent and with a constant gradient such that the total twist across the system in the direction is , namely . Within the mean-field treatment, one can show that addition of this gradient is equivalent to substituting . Now, the square of the gradient twist is related to the superfluid density via,rhoS (); hcb ()

 Ωtwisted−Ω=tρsδφ2, (13)

which in turn yields the simple expression

 ρs=−12d∂Ω∂t. (14)

Setting , this expression for the superfluid density coincides with that of :

 (15)

giving the critical value for the phase transition . Figures 8 and 9 show: (a) the free energy, (b) the superfluid density, (c) the density of bosons in the zero-momentum mode, and (d) the compressibility of the system as a function of in two and three dimensions. The dashed and solid curves represent the mean-field and SSE results, respectively. As one can immediately see, the critical values obtained within the mean-field approximation do not agree with the exact-numerical results. In two dimensions the error is and in three dimensions, it is . The very large errors here merely reflect the fact that the mean-field approach used here is not fit to describe the model at hand, especially in the vicinity of the SF-MI phase transition.

As pointed out earlier, the addition of spin-wave corrections yields virtually exact results in the homogeneous case in two dimensions.hcb () For the reader’s convenience, we review the mean-field calculations of the homogeneous () case and its spin-wave corrections in Appendix A (thereby also correcting some misprints that appeared in the original manuscript examining this case, Ref. hcb, ). Let us see how the mean-field results are modified by the addition of spin-wave corrections in our case. To include these, we proceed in the usual way.SW1 (); SW2 (); SW3 (); SW4 () We first introduce a set of local rotations that align the direction of each of the spins with its mean field orientation. This is accomplished by switching to new spin operators defined by

 (16)

where is the rotation matrix

 (17)

The corresponding new annihilation and creation operators and describe low-energy fluctuations about the mean-field ground state – these are spin waves. They too obey hardcore bosons commutation relations. Substituting these expressions into our Hamiltonian, and ignoring cubic and quartic terms in these bosonic operators (thus assuming a dilute gas of spin waves), the new Hamiltonian reads

 ^H{\scriptsize SW} = ^H{\scriptsize MF}+D∑i^b†i^bi+C∑i(−1)σ(i)^b†i^bi + B∑⟨ij⟩(^b†i^b†j+^bi^bj)−A2∑⟨ij⟩(^b†i^bj+^bi^b†j),

where the coefficients are

 A =t(1+cosθ1cosθ2), (19a) B =t/2(1−cosθ1cosθ2), (19b) C =dt(μ1cosθ1−μ2cosθ2), (19c) D =dt(2sinθ1sinθ2+μ1cosθ1+μ2cosθ2). (19d)

This quadratic Hamiltonian can be diagonalized by first going to Fourier space, using . This in turn yields the Hamiltonian:

 ^H{\scriptsize SW}=^H{\scriptsize MF% } + ∑k(D−Aγk)^b†k^bk+C∑k^b†k^bk+L/2 (20) + B∑kγk(^b†k^b†L−k+^bk^bL−k),

where, , and are the components of the momentum vector in each of the directions. We note that the Fourier-space operators and no longer obey the hardcore bosons commutation relations. These field operators are only excitations about the ground state, and are treated as soft-core bosons.hcb (); SW1 (); SW2 (); SW3 (); SW4 () At this point, our Hamiltonian may be diagonalized in a straightforward manner (we review the diagonalization process in Appendix B). Once diagonalized, the Hamiltonian takes the form

 ^H{\scriptsize SW}=^H{\scriptsize MF% }+∑kΛk^η†k^ηk+E0, (21)

where the ’s are energy levels and is the correction to the ground-state energy of the system, given by:

The operators and in Eq. (21) are modified spin-wave creation and annihilation operators, respectively, and are each a linear combination of and their adjoints. The coefficients of these linear combinations are fixed during the diagonalization process, and using them, all physical observables can be calculated in a straightforward manner (we elaborate on this matter in Appendix B).

The results of the spin-wave analysis are indicated by the dotted lines in Figs. 8 (two dimensions) and 9 (three dimensions). They show: (a) the free energy, (b) the superfluid density, (c) the density of bosons in the zero-momentum mode, and (d) the compressibility, after the addition of spin-wave corrections, as a function of .

As one can see in those figures, in the superfluid phase, the spin-wave corrected values for the free energy are almost on top of the exact-numerical ones; and more so in the three-dimensional case than in the two-dimensional one. As for the other measured observables, the spin-wave corrections are clearly an improvement over the mean-field results, especially for small values of where the spin-wave corrections yield virtually exact results. Unfortunately however, as one approaches the phase transition itself, the spin-wave corrections lose their accuracy, eventually leaving the phase-transition at its mean-field value, namely at .

Another issue worth noting here is the behavior of the spin-wave corrected superfluid density [Figs. 8(b) and 9(b)] in the vicinity of the predicted phase transition, . On the superfluid side of the transition the superfluid density becomes negative, indicating the breakdown of the spin-wave approximation for that quantity. The transition point is still signaled by a discontinuity in . However, the overall behavior of the superfluid density around the transition point is clearly an artifact of the spin-wave approximation and should not be considered further.

### iv.3 Improved mean-field approach

Having seen that spin-wave corrections, albeit accurate in the weak-potential regime, do not modify the critical point predicted by the mean-field solution, we have devised an improved mean-field approach. As we show now, this method provides a significant improvement over the mean-field results (and the spin-wave corrections) discussed previously, particularly in the context of the location of the critical point.

We start with a variational ansatz which, as before, is a product state. However, this time we do not choose a product of single-site wave-functions. The new ansatz is a product of wave-functions each describing the state of a ‘block’ of sites, such that with this block as the basic cell, the model turns homogeneous. In two dimensions a block consists of cells (as shown in Fig. 10) each of which is described by the general wave function

 |0⟩{\scriptsize IMF}=⊗∏blocks⎛⎝∑i,j,k,l∈{↓,↑}cijkl|ijkl⟩⎞⎠, (23)

where the generalization to three dimensions, in which case the basic block is a cell, is straightforward (note that the coefficients for each of the blocks are the same). As before, we minimize the free energy with respect to the coefficients of the wave function (this time we do so numerically). Obtaining the various observables in terms of the wave function given in Eq. (23) is straightforward, and was performed in much the same way as the usual mean-field approach discussed in Sec. IV.1. The results of this approximation are given by the dash-dotted lines in Figs. 8 and 9. They depict: (a) the free energy, (b) the superfluid density, (c) the density of bosons in the zero-momentum mode, and (d) the compressibility, as a function of .

As the figures indicate, in most instances, the results of this method are more accurate than those of the previous approximation schemes, in particular, for the location of the phase transition. The critical values given by this approximation are in two dimensions ( error) and in three dimensions ( error). Also, we note that while the spin-wave corrected values for the various thermodynamic quantities are a better approximation in the weak potential [small ] regime, as one moves away from this region, the improved mean-field technique proves to be a better estimator for all quantities but the free energy. It is clear still that the improved-mean-field method presented here is far from being very accurate.

## V Conclusions

We have studied the superfluid to Mott-insulator phase transition of hardcore bosons in a period-two superlattice in two and three dimensions. We focused on the case where the system is at half filling, for which the quantum phase transition belongs to the dimensional universality class.

Using quantum Monte Carlo simulations and finite size scaling, we have determined the critical value of the alternating potential parameter at which the SF-MI phase transition occurs. In two dimensions, our results agree with previous calculations.2dcrit ()

We have also compared our numerical results against several approximation schemes, some of which have been used successfully in the two-dimensional homogeneous version of the model. We have seen that employing a mean-field approach using the usual Gutzwiller ansatz works very poorly ( error in two dimensions and error in three dimensions). This is a clear indication of the fact that this mean-field approach is not suitable for describing this model, especially in the vicinity of the phase transition, as it breaks down in the strong coupling regime.

The spin-wave corrections to the mean-field solution turned out to be very useful, especially in the superfluid phase, where the spin-wave corrected estimation of the free-energy is very close to the exact values, and also reproduced the exact results for all observables for small values of . However, as one moves away from the weak potential regime, the spin-wave corrections become more and more inaccurate, and their predictions of the critical points eventually coincide with those of the mean-field approach, therefore indicating their unusefulness in that region.

The improved mean-field approximation scheme we have devised here, which was based on the underlying homogeneity of the problem, has proved to be an improvement over the previous methods, albeit still far from being accurate. This approach provides an analytical description of the superfluid-to-Mott-insulator transition and gives an estimate of the critical value for the transition with (around) one half the error of the usual Gutzwiller ansatz, i.e., it is an improvement in terms of the location of the critical point.

###### Acknowledgements.
This work was supported by the US Office of Naval Research Award No. N000140910966 and by startup funds from Georgetown University. We thank Tommaso Roscilde and Valery G. Rousseau for useful discussions.

## Appendix A Mean-field results and spin-wave corrections in the homogeneous case

In what follows, we briefly review the results of the mean-field approximation of Sec. IV.1 and its spin-wave corrections in the homogeneous () case for arbitrary values of .

Starting with the ansatz given in Eq. (7), minimization of the free energy Eq. (8) with respect to the spin orientation angles yields

 cosθj=μ2dt, (24)

where the azimuthal angle takes on, once again, a constant yet arbitrary value . The density of particles becomes

 ρ{\scriptsize MF} = 1N∑i{\scriptsize MF}⟨0|^a†i^ai|0⟩{\scriptsize MF}=12+12N∑icosθi (25) = 12(1+μ2dt),

and the free energy is

 Ω{\scriptsize MF} = {\scriptsize MF}⟨0|^H|0⟩{\scriptsize MF}=−dt2sin2θ−12μ(1+cosθ) (26) = −12dt(1+μ2dt)2.

The density of bosons in the zero-momentum mode turns out to be

 ρ0,{\scriptsize MF} = 1N{\scriptsize MF}⟨0|^a†k=0^ak=0|0⟩{\scriptsize MF} = 14N2∑i,jsinθisinθj=14[1−(μ2dt)2].

Using Eq. (14), it can be easily shown that the expression for the superfluid density in the homogeneous case coincides with the expression obtained for above (as is the case with the alternating potential).

The addition of spin-wave corrections to the mean-field results is carried out in exactly the same manner as with the staggered potential. The Hamiltonian in this case has the same form as the one given in Eq. (IV.2) but with coefficients

 A =t[1+(μ2dt)2], (28a) B =t2[1−(μ2dt)2], (28b) C =0, (28c) D =2dt. (28d)

The spin-wave field operators which diagonalize the Hamiltonian are given by the simple relation:hcb ()

 ^bk=coshϕk^ηk−sinhϕk^η†L−k, (29)

with obeying

 sinh2ϕk =12(D−Aγk√(D−Aγk)2−(2Bγk)2−1), (30a) cosh2ϕk =12(D−Aγk√(D−Aγk)2−(2Bγk)2+1), (30b)

and so, the various spin-wave corrected quantities may be written explicitly: the corrected density of particles is

 ρ{\scriptsize SW}=ρ{\scriptsize MF}−1Nμ2dt∑k≠0sinh2ϕk, (31)

and the free energy becomes

 Ω{\scriptsize SW} = Ω{\scriptsize MF} + 12∑k≠0[√(D−Aγk)2−(2Bγk)2−(D−Aγk)].

Using Eq. (14), the superfluid density immediately follows. Finally, the density of bosons in the zero-momentum mode turns out to be:

 ρ0,{\scriptsize SW}=ρ0,{\scriptsize MF% }−1N[1−(μ2dt)2]∑k≠0sinh2ϕk.

## Appendix B Diagonalization of quadratic bosonic Hamiltonians

Following the prescription given in Ref. trans, for the diagonalization of fermionic quadratic Hamiltonians, we provide here the analogous prescription for the diagonalization of bosonic quadratic Hamiltonians of the general form

 ^H=∑k,m(Akm^b†k^bm+12Bkm(^b†k^b†m+^bk^bm)), (34)

where and are bosonic annihilation and creation operators, respectively, and and are real-valued and symmetric. For the spin-wave Hamiltonian of Eq. (20), the coefficients are

 Akm =(D−Aγk)δkm+Cδk,m+L/2, (35a) Bkm =2Bγkδk,L−m. (35b)

The diagonalization process starts by defining the following linear transformation

 ^ηk =∑m(gkm^bm+hkm^b†m), (36)

where and are real-valued and we require and be bosonic operators. This is enforced by the constraint

 δkm=[ηk,η†m]=∑l(gklgml−hklhml). (37)

The coefficients and are determined in such a way that the transformed Hamiltonian takes the diagonal form

 H=∑kΛk^η†k^ηk+E0, (38)

once the ’s are substituted for the ’s, and . As the new Hamiltonian is already in diagonal form, the new field operators obey the eigenvalue equation

 [^ηk,^H]=Λk^ηk. (39)

Plugging in the transformations given in Eqs. (36), we obtain the relations

 Λkgkm =∑l(gklAlm−hklBml), (40a) Λkhkm =∑l(gklBml−hklAlm). (40b)

These relations may be further simplified by defining the new coefficients

 ϕkm =gkm+hkm, (41a) ψkm =gkm−hkm, (41b)

for which, the constraint (37) translates to

 12∑l(ϕklψml+ψklϕml)=δkm. (42)

With the above definitions, Eqs. (40) may be cast in vector notation:

 \boldmathϕk(\boldmathA−% \boldmathB) =Λk\boldmathψk, (43a) \boldmathψk(\boldmathA+% \boldmathB) =Λk\boldmathϕk. (43b)

These can be solved by simply plugging each of these equations into the other, resulting in the eigenvalue equations

 \boldmathψk(\boldmathA+% \boldmathB)(\boldmathA−\boldmathB) =Λ2k\boldmathψk, (44a) \boldmathϕk(\boldmathA−% \boldmathB)(\boldmathA+\boldmathB) =Λ2k\boldmathϕk. (44b)

These equations are to be solved by standard techniques. Once the ’s, ’s and ’s are found, all physical observables can be readily calculated: First, the observable of interest should be expressed in terms of normal-ordered ’s. This may be accomplished by using the inverse of the transformation given in Eq. (36):

As a next step, one should use the fact that as excitations, the ’s obey . This leads to

 {\scriptsize MF}⟨0|^b†k^bm|0⟩{\scriptsize MF}=14∑l(ϕ−1kl−ψ−1kl)(ϕ−1ml−ψ−1ml).

As an example, consider the spin-wave corrected density of particles in our model. It is calculated as

 ρ{\scriptsize SW} = 1N∑i{\scriptsize MF}⟨0|^a†i^ai|0⟩{\scriptsize MF}=ρ{\scriptsize MF}−1N∑i{\scriptsize MF}⟨0|^b†i^bicosθi|0⟩{% \scriptsize MF} = ρ{\scriptsize MF}−12N(cosθ1+cosθ2)∑k{\scriptsize MF}⟨0|^b†k^bk|0⟩{\scriptsize MF}−12N(cosθ1−cosθ2)∑k{\scriptsize MF% }⟨0|^b†k^bk+L/2|0⟩{% \scriptsize MF} = ρ{\scriptsize MF}−18N((cosθ1+cosθ2)∑mk(ϕ−1km−ψ−1km)2+(cosθ1−cosθ2)∑mk(ϕ−1km−ψ−1km)(ϕ−1(k+L/2),m−ψ−1(k+L/2),m)).

All other observables may be calculated in the same manner.

## References

• (1) M. Greiner, O. Mandel, T. Esslinger, T. E. Hänsch, and I. Bloch, Nature (London) 415, 39 (2002).
• (2) M. Greiner, O. Mandel, T. E. Hänsch, and I. Bloch, Nature (London) 419, 51 (2002).
• (3) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, England, 1999).
• (4) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
• (5) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
• (6) W. Zwerger, J. Opt. B: Quantum Semiclassical, Opt. 5, 9 (2003).
• (7) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
• (8) E. H. Lieb and F. Y. Wu, Physica A 321, 1 (2003).
• (9) J. K. Freericks and H. Monien, Phys. Rev. B 53, 2691 (1996).
• (10) N. Elstner and H. Monien, Phys. Rev. B 59, 12184 (1999).
• (11) A. P. Kampf and G. T. Zimanyi, Phys. Rev. B 47, 279 (1993).
• (12) K. Sheshadri, H. R. Krishnamurthy, R. Pandit, and T. V. Ramakrishnan, Europhys. Lett. 22, 257 (1993).
• (13) T. P. Polak and T. K. Kopeć, Phys. Rev. B 76, 094503 (2007).
• (14) N. Teichmann, D. Hinrichs, M. Holthaus, and A. Eckardt, Phys. Rev. B 79, 100503(R) (2009).
• (15) M. Aizenman, E. H. Lieb, R. Seiringer, J. P. Solovej, and J. Yngvason, Phys. Rev. A 70, 023612 (2004).
• (16) V. G. Rousseau, D. P. Arovas, M. Rigol, F. Hébert, G. G. Batrouni, and R. T. Scalettar, Phys. Rev. B 73, 174516 (2006).
• (17) M. Rigol and A. Muramatsu, Phys. Rev. A 70, 031603(R) (2004); Phys. Rev. A 72, 013604 (2005).
• (18) M. Rigol, A. Muramatsu, and M. Olshanii, Phys. Rev. A 74, 053616 (2006).
• (19) A. W. Sandvik, Phys. Rev. B 59, R14157 (1999).
• (20) A. Dorneich and M. Troyer, Phys. Rev. E 64, 066701 (2001).
• (21) K. Bernardet, G. G. Batrouni, J.-L. Meunier, G. Schmid, M. Troyer and A. Dorneich, Phys. Rev. B 65, 104519 (2002).
• (22) C. Bruder, R. Fazio and G. Schön, Annalen der Physik 14, 566 (2005).
• (23) S. Peil, J. V. Porto, B. Laburthe Tolra, J. M. Obrecht, B. E. King, M. Subbotin, S. L. Rolston, and W. D. Phillips, Phys. Rev. A 67, 051603(R) (2003).
• (24) J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. Porto, Phys. Rev. A 73, 033605 (2006).
• (25) J. Sebby-Strabley, B. L. Brown, M. Anderlini, P. J. Lee, W. D. Phillips, J. V. Porto, and P. R. Johnson, Phys. Rev. Lett. 98, 200405 (2007).
• (26) P. J. Lee, M. Anderlini, B. L. Brown, J. Sebby-Strabley, W. D. Phillips, and J. V. Porto, Phys. Rev. Lett. 99, 020402 (2007).
• (27) A. Priyadarshee, S. Chandrasekharan, J.-W. Lee, and H. U. Baranger, Phys. Rev. Lett. 97, 115703 (2006).
• (28) T. Matsubara and H. Matsuda, Prog. Theor. Phys. 16, 569 (1956).
• (29) M. E. Fisher, M. N. Barber, and D. Jasnow, Phys. Rev. A 8, 1111 (1973).
• (30) K. S. Liu and M. E. Fisher, J. Low Temp. Phys. 10, 655 (1973).
• (31) Yi-Chen Cheng, Phys. Rev. B 23, 157 (1981).
• (32) R. T. Scalettar, G. G. Batrouni, A. P. Kampf, and G. T. Zimanyi, Phys. Rev. B 51, 8467 (1995).
• (33) G. Murthy, D. Arovas, and A. Auerbach, Phys. Rev. B 55, 3104 (1997).
• (34) E. Lieb, T. Schultz and D. Mattis, Annals of Physics 16, 407 (1961).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters