{'fraction_non_alphanumeric': 0.03540615512308062, 'fraction_numerical': 0.007186073433370465, 'mean_word_length': 4.6558369363804815, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical and numerical study of trapped strongly correlated bosons in two- and three-dimensional lattices

We study the ground-state properties of trapped inhomogeneous systems of hardcore bosons in two- and three-dimensional lattices. We obtain our results both numerically, using quantum Monte Carlo techniques, and via several analytical approximation schemes, such as the Gutzwiller-mean-field approach, a cluster-mean-field method and a spin-wave analysis which takes quantum fluctuations into account. We first study the homogeneous case, for which simple analytical expressions are obtained for all observables of interest, and compare the results with the numerical ones. We obtain the equation of state of the system along with other thermodynamic properties such as the free energy, kinetic energy, superfluid density, condensate fraction and compressibility. In the presence of a trap, superfluid and insulating domains coexist in the system. We show that the spin-wave-based method reproduces the quantum Monte-Carlo results for global as well as for local quantities with a high degree of accuracy. We also discuss the validity of the local density approximation in those systems. Our analysis can be used to describe bosons in optical lattices where the onsite interaction U is much larger than the hopping amplitude t.


I. INTRODUCTION
The successful realization of the superfluid-to-Mottinsulator transition in ultracold bosonic gases trapped in optical lattices in one [1], two [2,3], and three [4] dimensions has opened the way to a myriad of experimental and theoretical studies of strongly correlated lattice systems [5]. One of the main ideas driving these studies is that ultracold atoms in optical lattices can be used as analogue simulators of Hamiltonians of the Fermi-and Bose-Hubbard type. Intensive efforts are currently under way to validate this approach by comparing experimental and theoretical results for systems such as the Bose-Hubbard model which is amenable to both treatments [6,7].
In the Bose-Hubbard model, the basic Hamiltonian consists of a hopping term of amplitude t and an onsite two-body repulsion term with amplitude U . The phase diagram of this homogeneous model is known to consist of two phases: (i) a superfluid phase that is present for all incommensurate fillings and arbitrary values of the ratio U/t, and for commensurate fillings below some critical value (U/t) c , which depends on the dimensionality of the system and on the (integer) filling, and (ii) a Mott insulator which is present for commensurate fillings for U/t > (U/t) c [8][9][10][11][12].
In experiments with ultracold gases, a confining potential is always present. This confining potential is to a good approximation harmonic, and generates an inhomogeneous density profile in which superfluid and Mottinsulating phases coexist in spatially-separated domains [13][14][15][16][17][18]. In such systems, the appearance of Mott insulating domains in different regions of the trap depends not only on the total filling and the ratio U/t (as in the homogeneous case) but also on the curvature of the har-monic confining potential.
So far, an accurate characterization of the local properties of the system in the trap has only been achieved by means of quantum Monte Carlo (QMC) simulations [14,15,17,18], and density-matrix renormalization group (DMRG) techniques in one dimension [16]. Our goal in this paper is to introduce an analytical approach that provides an accurate prediction of the behavior of the system in two and three dimensions. In this study we shall consider the strongly correlated limit where the ratio U/t is very large and the local occupation of the lattice sites is lower than or equal to one. In that case, the system can be described, to a good approximation, by impenetrable (hardcore) bosons. As we will show, this in turn enables the implementation of an analytical treatment which goes beyond simple mean-field calculations, where within this approach, quantum fluctuations are taken into account by the addition of spin-wave corrections.
In what follows, we shall examine the extent to which the results of this method are valid by comparing them against quantum Monte Carlo simulations. As we will show, spin-wave-corrected results provide a big improvement over the usual Gutzwiller-mean-field approach and its cluster-mean-field extension. In most cases, for the system sizes considered here, they are also more accurate than the QMC-based local density approximation (LDA). In some regimes, and for some local observables, the latter is found to be less accurate than the simple Gutzwiller-mean-field approach.
The paper is organized as follows. In Sec. II we review the model at hand and discuss its properties. In Sec. III, we study homogeneous systems in two and three dimensions. We obtain the equation of state along with basic thermodynamic properties and examine the extent to which several analytical approaches compare against exact numerical QMC simulations. In Sec. IV, we analyze the global and local properties of hardcore bosons trapped in a confining harmonic potential in a two-dimensional setup. We also examine the validity of the local-density approximation for finite systems and compare the various analytical approximations against exact QMC simulations. In Sec. V, we analyze the behavior of three-dimensional harmonically trapped hardcore bosons using the various approximation schemes. Finally, in Sec. VI, we conclude with a discussion and summary of our results.

II. MODEL
The Hamiltonian for hardcore bosons confined in a harmonic potential in a d-dimensional hypercubic lattice with N = L d sites, can be written as: Here, ij denotes nearest neighbors andâ i (â † i ) destroys (creates) a hardcore boson on site i located at a distance r i = |r i | from the center of the trap, where the coordinates r i = (x 1i , . . . , x di ) are given here in units of the lattice spacing a, which we set to unity.
The operatorn i =â † iâ i is the local density operator, µ is the global chemical potential, and V is the harmonic potential strength. The hopping parameter t (which we shall fix at t = 1) sets the energy scale. The hardcore boson creation and annihilation operators satisfy the constraintsâ †2 i =â 2 i = 0, which prohibit double or higher occupancy of lattice sites, as dictated by the U → ∞ limit of the Bose-Hubbard model. For any two different sites i = j, the creation and annihilation operators obey the usual bosonic relations To gain a general understanding of the zerotemperature phases of this model, let us first analyze the atomic (t = 0) limit. In this limit, there is no kinetic (hopping) term, and the boson number operatorŝ n i commute individually with the Hamiltonian, so every lattice site is occupied by a fixed number of bosons. In this case, the Hamiltonian is diagonal in the Fock-states basis, and the ground-state wave-function has the form |GS = ⊗ i |ψ i and is solved for each site separately, giving where we have denoted µ i = µ − V r 2 i as the local chemical potential, and |0 and |1 denote vacant and occupied sites, respectively. In the special V = 0 case, where no trap is present, the model is translationally invariant, and the ground-state boson occupancy is the same throughout the lattice: For µ < 0 the minimal energy configuration is simply the vacuum, i.e., the completely-empty lattice, and for µ > 0 the minimal energy configuration is the completely-filled lattice. The ground-state energy of these phases is degenerate at µ = 0.
In the t = 0 case, the model has in general no analytic solution [19]. As it turns out, however, the phase boundaries separating the insulating phases from the superfluid one can easily be obtained analytically in the V = 0 case even for nonzero t. To see this, we use the fact that our Hamiltonian commutes with the total-number-of-bosons operatorN b = in i . This simply means that for any given µ and t, the ground-state wave-function is a linear combination of product states each having the same number of occupied (vacant) sites. In the completelyfilled phase (which we denote by F) the wave-function is simply with energy ε F = −µN . In the infinitesimally thin layer outside of the F phase, the state of the system is characterized by exactly one vacant site, and thus its wavefunction is the sum with energy ε F + = −2dt − µ(N − 1). The phase boundary separating the superfluid phase and the (insulating) completely filled phase is the curve along which the F state, Eq. (3), is no longer energetically favorable. This happens when its energy becomes degenerate with the energy of the defect state, Eq. (4). Matching the two, we obtain the phase boundary: By the same token, the boundary between the superfluid phase and the completely empty insulating phase (which we denote by E) is given by µ/(2dt) = −1. This result can also be obtained by repeating the above exercise with the substitution |0 ↔ |1 . Between the two insulating phases, in the region |µ/2dt| < 1, the system is superfluid. When a trap is introduced, the system will no longer be entirely superfluid or entirely insulating. Depending on the trap curvature and the total chemical potential, the system may be in a coexistent state of a superfluid in one part of the lattice and an insulator in another. Before discussing why this is so, let us first note that unlike the homogeneous case which is scale-invariant, in the presence of a trap the model has a length scale determined by the ratio of the trap curvature to the hopping amplitude, given by ξ = (V /t) −1/2 . As a result, our Hamiltonian can be rescaled by introducing the dimensionless scaled lengthr i ≡ r i /ξ through which the Hamiltonian can be rewritten in a way that eliminates the amplitude of the harmonic trap V : As a result of this rescaling, it is clear that systems with the same global chemical potential exhibit the same qualitative behavior for all positive trap curvatures. In addition, away from the center of the trap, the local onsite potential µ i always becomes very large and negative and hence outside of some radiusr out , the lattice must be empty ( n i = 0). On the other hand, in the center of the trap, large enough local chemical potentials will produce an insulating (circular) domain with n i = 1 inside some critical radiusr ≤r in . In the intermediate region, r in ≤r ≤r out , the system is in the so-called superfluid regime. The radial density profiles of the two possible scenarios, with and without the insulator in the center of the trap, are sketched in Fig. 1.
in ). Outside of this radius the system is superfluid up untilr =r (1) out , outside of which the lattice is empty. If the chemical potential is small (dashed line) a Mott insulator will not form in the center of the trap, and the only superfluid-insulator transition will take place atr =r (2) out .

III. HOMOGENEOUS (V = 0) CASE IN TWO AND THREE DIMENSIONS
Having discussed the general properties of the model in the previous section, we now turn to analyze it in a quantitative manner. We shall try to keep the discussion as general as possible so that it applies to arbitrary dimensions and confining potentials. Only after the general formulation has been introduced, we will obtain analytical expressions for homogeneous (i.e., V = 0) systems in two and three dimensions. (We note that some of the observables of the homogeneous two-dimensional case we analyze here have been studied in Ref. [20].) This general formulation will be useful later on when we address the problem of harmonically-confined bosons.
Our main objective in this section is to obtain the equation of state of the model, i.e., the average density as a function of the chemical potential, along with the basic thermodynamic properties of the system. Since in dimensions higher than one the model has no analytic solution, exact results for finite systems (up to the relevant statistical errors) are obtained numerically. Here we use the stochastic series expansion (SSE) algorithm [21,22] and perform simulations over a wide range of chemical potentials. Zero-temperature properties of the system are obtained by choosing large inverse-temperatures β = 1/T (in our units, k B = 1), where we have found it sufficient to have β ≥ L, as the effects of increasing β beyond this value are indiscernible. In two dimensions, we simulate systems with 48 × 48 sites and 60 × 60 sites, with an inverse-temperature of β = 60, and in three dimensions, simulations are performed on a 16 × 16 × 16 lattice and β = 20. After obtaining the QMC results, we examine the validity of several approximation schemes by comparing them against the results of the SSE simulations.

A. Quantum Monte Carlo results
We first study the equation of state, which is very relevant for experiments as it provides the dependence of the number of bosons in the system (equivalently, the average density) on the global chemical potential. It is depicted in Figs. 2(a) and 3(a) for two-and three-dimensional lattices, respectively. The corresponding compressibilities defined by κ = ∂ρ/∂µ are also given, in Figs. 2(b) and 3(b). These provide the response of the average density to a change in the chemical potential. We have computed them in two independent manners: (i) as the numerical derivative of the average density with respect to the chemical potential and (ii) within the SSE algorithm compressibility κ = ∂ρ/∂µ for hardcore bosons in a homogeneous two-dimensional lattice (60×60 sites). The SSE results are denoted by circles, whereas the dotted, dashed, and solid lines, are the results of the Gutzwiller-mean-field, the cluster-meanfield, and the spin-wave-corrected Gutzwiller-mean-field approximations, respectively. using the formula: where Z is the partition function and we have taken advantage of the fact that [N b ,Ĥ] = 0. As expected, the results of both methods were found to coincide. In the figures, the SSE results are denoted by circles, whereas the various lines indicate the predictions obtained by several approximation schemes that will be introduced later.
As the figures indicate, at µ/(2dt) = 0 the average density is 1/2. This is because of the particle-hole symmetry of the model. Upon increasing the chemical potential, the boson density increases too until it reaches ρ = 1 at µ/(2dt) = 1, as predicted by the analysis given in the previous section. At ρ = 1 the system becomes insulating and the compressibilities drop sharply to zero. Note that due to the particle-hole symmetry, the results for negative µ can be immediately read off the figures by the transformations µ → −µ and ρ → 1 − ρ. ] further indicate that for ρ < 1, when the systems are compressible, they are also superfluid and exhibit Bose-Einstein condensation. In the presence of the strong correlations generated by the infinite onsite repulsion, one can see that the superfluid density is always greater than the condensate density. This is similar to what happens in liquid helium where strong correlations deplete the condensate density to a very small value while the superfluid density remains very large. Interestingly, as the dimensionality of the system increases one can see that the difference between the superfluid density and the condensate density decreases. This point will be touched upon later, when we discuss the Gutzwiller-mean-field method which corresponds to the exact solution in infinite dimensions. Also evident from the figures is the fact that the superfluid density and the condensate density become equal in the low density limit.

B. Approximation schemes
Having analyzed the homogeneous case via the SSE algorithm, we now proceed to analyze the model with a number of approximation schemes. We start this investigation with the Gutzwiller mean-field approach. We note here that this approach is a particular case of a more general treatment initially developed for soft-core bosons [8,13]. , whereas the dotted, dashed, and solid lines indicate the Gutzwiller-mean-field, the cluster-mean-field, and the spinwave-corrected Gutzwiller-mean-field approximations, respectively.

Gutzwiller mean-field
Generalizing the discussion in Refs. [23] and [24] to inhomogeneous systems, we start our mean-field calculation with the following product state as our ansatz for the ground state wave-function: where the unknowns (θ j , ϕ j ) are to be determined via the minimization of the grand-canonical potential (per site) which evaluates to In the homogeneous case, the wave-functions of each of the lattice sites are identical. However, in the following, we consider general local potentials, as later we shall apply this method to the case of the harmonic potential. For the azimuthal angles ϕ j , it is sufficient to require that all sites share the same value φ j = Φ, where Φ can be chosen arbitrarily. The conditions for the polar angles are obtained by differentiating Ω MF with respect to θ j , giving: Here, S i = j sin θ j δ i,n(j) where δ i,n(j) is unity if i and j are neighbors and is zero otherwise. This set of equations is to be solved numerically in the general case.
In the limit where the lattice spacing is very small compared to the average interparticle distance, i.e., in the low density limit, the θ j 's may be assumed to vary smoothly over the lattice and the above set of equations may be approximated by its continuous version, which turns out to be the ordinary differential equation which also has to solved numerically in the general case.
In the homogeneous case, where V = 0 and the system is translationally invariant, Eq. (10) can be solved exactly, as we can write θ j = θ from which it follows that Now that the ground state wave-function has been obtained, all physical properties of interest can easily be calculated by computing the expectation value of the appropriate operator. The average density of particles, for example, becomes where in the homogeneous case, this expression simplifies to The density of bosons is plotted in Figs. 2 and 3 for two and three dimensions, respectively, along with the corresponding compressibilities, as a function of the chemical potential. As can be seen in Figs. 2(a) and 3(a), when compared against the quantum Monte Carlo results, the Gutzwiller mean-field approach provides a good approximation for the density close to half filling and again at µ/(2dt) = 1 where the lattice becomes completely filled, but deviates from the QMC results in between. The behavior of the compressibility, as predicted by the Gutzwiller mean-field approach, is on the other hand qualitatively incorrect for all densities in the superfluid regime. Next, the free energy, Eq. (9), becomes in the homogeneous case and the density of bosons in the zero-momentum mode (the condensate density) turns out to be whereâ † k=0 (â k=0 ) creates (destroys) a particle with momentum k. The superfluid density, ρ s,MF , requires a special treatment of the boundary conditions. As is well known [25], the superfluid density is related to the "spin stiffness" of the system: 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 x direction). In the homogeneous case we consider here, the azimuthal angles ϕ j are 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 x direction is π, namely δϕ = ϕ j+x − ϕ j = π/L. Within the mean-field treatment, one can show that addition of this gradient is equivalent to substituting t → t/d [(d − 1) + cos δϕ]. Now, the square of the gradient twist is related to the superfluid density via the relation Ω twisted − Ω = tρ s δϕ 2 which in turn yields the simple expression At the mean-field level, this expression evaluates to that of the condensate density. Since mean-field theory becomes exact in infinite dimensions, the differences between both observables should decrease as the dimension of the systems increases. This explains why the differences between the superfluid density and the condensate density (as they were computed with the QMC simulations in Sec. III A) were smaller in three dimensions than they were in two dimensions. The quantities discussed in Eqs. (15)- (17) are plotted in Figs. 4 and 5 for two and three dimensions, respectively. As one can immediately see, the Gutzwiller-meanfield predictions agree rather poorly with the QMC data with respect to all quantities but the free energy, particularly close to half-filling. Hence, the mean-field approach described in this section does not provide an accurate picture for most observables of interest close to half filling. Nonetheless, it is clear that the agreement with the QMC data is better in three dimensions: As already noted, the mean-field approximation becomes exact in the limit where d → ∞.

Cluster mean-field
The cluster-mean-field approach is an improvement over the simple Gutzwiller-mean-field approach, which is also applicable to harmonically trapped hardcore bosons. It was introduced in Refs. [23] and [24] to improve the Gutzwiller-mean-field prediction for the phase diagram of hardcore bosons when a period-two superlattice was added to the homogeneous model. Within this approach, one starts with a variational ansatz which is also a product state, but not of single-site wave-functions. The new ansatz is a product of wave-functions each describing the state of a 'block' of 2 d sites. In two dimensions, a block consists of 2 × 2 square cells, so the ground state wavefunction has the form: Analogously, in three dimensions, the basic block is a 2 × 2 × 2 cubic cell [23].
As with the Gutzwiller-mean-field case, we minimize the free energy Ω CMF = N −1 CMF GS|Ĥ|GS CMF with respect to the coefficients c b ijkl of the wave-function. Obtaining the various observables in terms of the wave-function given in Eq. (18) is straightforward, and is performed in much the same way as with the usual mean-field approach discussed in Sec. III B 1.
The two-and three-dimensional equations of state, as determined by the cluster mean-field approximation, are depicted by the dashed lines in Figs. 2 and 3, respectively. The predictions for the other observables are presented in Figs. 4 and 5. As the figures indicate, the clustermean-field predictions are an improvement over the plain Gutzwiller-mean-field approach but are still far from accurately reproducing the results of the QMC simulations.

Addition of spin-wave corrections
In what follows, we show that the results of the meanfield approaches discussed above can be significantly improved by the addition of spin-wave corrections which take into account quantum fluctuations [26][27][28][29]. This method was also discussed in the context of homogeneous two-dimensional systems in Ref. [20]. In a later study [23], it was shown that in the presence of a superlattice, spin wave corrections do not modify the predictions of the simple Gutzwiller-mean-field theory for the phase diagram. Here, we generalize this method to arbitrary dimensions and general inhomogeneous systems.
The spin-wave analysis begins with introducing a set of local rotations to the boson creation and annihilation operators. This is accomplished by switching to new field operators via the transformation: The new annihilation and creation operatorsb j andb † j 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, ignoring cubic and quartic terms in these bosonic operators (thus assuming a dilute gas of spin waves), the new Hamiltonian becomes a sum of a constant termĤ 0 , a linear termĤ 1 and a quadratic termĤ 2 , wherê and the matrix coefficients inĤ 2 are: The spin-waves analysis continues with the determination of the θ i 's via the requirement that the linear term H 1 , Eq. (20), vanishes. Equating each term in the sum to zero, one ends up with the set of equations: which is analogous to the set of equations obtained previously when the Gutzwiller-mean-field approach was discussed, and which similarly has to be solved numerically in the general case. In the homogeneous case where µ i = µ, we have θ i = θ and therefore also S i = 2d sin θ. This leads to the solution: Substituting the solution into the zeroth order Hamilto-nianĤ 0 , Eq. (20), the zeroth order free energy yields, as it should, the free energy of the Gutzwiller mean-field approach: In the homogeneous case, the above expression evaluates to:Ĥ which is precisely the mean-field result, Eq. (15). In order to obtain the spin-wave corrections, one must diagonalize the remaining quadratic termĤ 2 in Eq. (20), for which one assumes that the spin waves are very dilute and hence that the hardcore constraint can be ignored. The diagonalization process is a lengthy but straightforward process. Here, we shall only review the basic steps. The interested reader is referred to Appendix B of Ref. [23] for a more detailed account. The diagonalization process consists of the following steps: (i) Diagonalize the matrix (A + B)(A − B) and obtain its (non-negative) eigenvalues Λ 2 k along with the unnormalized eigenstates φ k = {φ k1 , . . . , φ kl , . . . , φ kN } where k = 1 . . . N (here, A and B are matrices whose elements are A ij and B ij , respectively). (ii) Evaluate the set of vectors ψ k = {ψ k1 , . . . , ψ kl , . . . , ψ kN } through the relation (A − B)φ k = Λ k ψ k . (iii) Normalize φ k and ψ k according to k φ k · ψ k = 1. (iv) Obtain the basic twoparticle operators through the inverses of φ k and ψ k : (v) Express any desired physical quantity in terms of the new field operators using the transformation, Eq. (19), neglecting cubic and quartic terms, and use the expressions for the two-particle operators, Eqs. (26), to evaluate the resulting expression.
As an example, consider the local density of particles. Within the spin-wave analysis, it evaluates to following which the average density becomes: The potential energy (per site) is then: The off-diagonal terms are more cumbersome to calculate but can still be obtained in a straightforward manner. They evaluate to: Once these are calculated, one can obtain the kinetic energy (per site): and also the momentum distribution function defined by: In the homogeneous case, the above equations simplify substantially, and one can obtain analytic expressions for the various physical observables [23]. The spin-wave corrected density of particles is given by: where we have defined The spin-wave corrected density is plotted in Figs. 2 and 3 for the two-and three-dimensional cases, respectively. As is evident from the two figures, the spin-wave corrected densities provide a decisive improvement over the plain Gutzwiller-mean-field results and they are also closer to the exact results than the cluster mean-field calculations for the most part. Interestingly, we find that for the compressibility, the results of the cluster meanfield approach are slightly better than those provided by the spin-wave corrections in three dimensions. The spin-wave corrections also yield the following expression for the free energy: where the superfluid density is obtained by evaluating the expression ρ s,SW = −(2d) −1 ∂Ω SW /∂t. Finally, the density of bosons in the zero-momentum mode is found to be: The above quantities are plotted in Figs. 4 (two dimensions) and 5 (three dimensions). As is clear from the figures, not only are the spin-wave corrected predictions a major improvement over the mean-field approaches for the thermodynamic quantities, but they are essentially on top of the QMC results for the kinetic energy and the superfluid density, and very close to the QMC results for the free energy and the condensate density.

IV. TRAPPED HARDCORE BOSON IN TWO DIMENSIONS
Having discussed the properties of the homogeneous case, we are now ready to analyze the more experimentally-relevant case of lattice bosons in the presence of a harmonic trap. We shall limit the discussion to the two-dimensional case in this section and address the three-dimensional case in the next section.
As discussed in Sec. II, in the presence of a harmonic potential, the trapping amplitude V can be scaled out by using the dimensionless parameter ξ = (V /t) −1/2 . In order to verify that this is indeed the case, we have performed simulations over a range of chemical potentials for two different trapping amplitudes, V = 0.01 and V = 0.02. We have verified that once properly normalized, results for local quantities for the two trapping potentials are virtually the same.
In much the same way we analyzed the homogeneous case in the previous section, here we obtain the equation of state of the model along with the various thermodynamic properties in the trapped system. We examine the accuracy of the analytical approaches discussed in the previous section in this inhomogeneous setup and compare them against QMC simulations in a trap. (For the trapped bosons, however, the equations derived in the previous section cannot be reduced to simple analytical expressions as in the homogeneous case; they must be computed numerically.) In addition, we explore the regimes of validity of another approximate method, namely, the LDA which is based on the QMC results of the homogeneous system.

A. Equation of state
In a trapped system, where local densities change with position, the meaning of 'equation of state' may become somewhat unclear. This is because in an inhomogeneous trapped system the local phases are determined not only by the onsite interaction and total filling of the system but also through the strength of the confining potential. Fortunately, it turns out that the total filling and the curvature of the confining potential come into play only as a combination of the two: This is the characteristic density. It is the trap-equivalent of the overall density in homogeneous systems.
The characteristic density (ρ) was introduced for lattice fermions in one dimension in Refs. [31,32], where it was shown by means of QMC simulations, that the scaled dimensionless variableρ = ρN (a/ξ) d can be used to generate a state diagram of the coexisting phases in the trap in the planeρ vs U/t in a way that is independent of the specific values of the number of particles and the amplitude of the confining potential. Further studies validated this approach for the one-dimensional fermionic case within Bethe-Ansatz schemes [33][34][35][36], and state diagrams have been obtained for fermions in higher dimensions by means of dynamical mean-field theory (DMFT) [37]. The same idea has been shown to work with equal success in bosonic systems [18], in the framework of which the concept of characteristic density has been further discussed within the local density approximation [38] and finite-size scaling [39,40]. It should be noted however that the concept of characteristic density remains a meaningful quantity even in cases where the LDA fails [18].
As follows from the above discussion, one can then generate equations of state for trapped systems based on the characteristic density. This has recently been done by Roscilde in Ref. [30] for soft-core bosons in various superlattice potentials. The equation of state for two-dimensional hardcore bosons is given in Fig. 6(a), which shows the dependence of the characteristic density in our two-dimensional trapped systems in terms of the global chemical potential. The SSE results are denoted in the figure by circles, whereas the various approximation schemes are represented by the different lines. The inset shows the deviations of the various approximation schemes from the QMC data. It is clear from the figure that, at least to some extent, all approximation schemes without exception provide a fairly accurate description of the equation of state. The spin-wave corrected results however provide a better match than the cluster-meanfield approximation which in turn shows only negligible improvement over the Gutzwiller-mean-field results. Fig-ure 6(b) shows the density in the center of the trap as a function of the characteristic density, a relation that is very useful to experimentalists for controlling their density profiles by only changing the total filling in their systems once the value of the trapping potential is known. The circles indicate the QMC results whereas the dotted, dashed, solid, and dot-dashed lines are the Gutzwillermean-field, the cluster-mean-field, the Gutzwiller-mean-field with spin-wave corrections, and the LDA results, respectively. The inset in panel (a) shows the deviations of the various approximation schemes from the SSE data. Figure 6 also shows another useful approximation method that is usually applied to trapped systems. This is the LDA, which provides predictions of local quantities based on results obtained from matching homogeneous systems. Within the LDA, local observables of the confined system are approximated by their values in the corresponding homogeneous system, where the chemical potential of the homogeneous case is taken to be the corresponding local chemical potential in the trap. In general, the LDA has been shown to give reasonably accurate descriptions of other confined systems, but is known to fail in regions where the local potential changes rapidly when compared with the correlation length of the relevant phase in the homogeneous case. The LDA results here were generated from the QMC results of the homogeneous systems studied in the previous section. As the figure illustrates, the LDA provides a very good approximation of the QMC results for the density in the center of the trap, although, as expected, it underestimates the characteristic density required to form a Mott insulator in the center of the trap [ Fig. 6(b)].

B. Local densities and compressibilities
In recent experiments, it has become possible to image local density profiles in optical lattices in-situ [41][42][43]. This now allows experimentalists to explore the behavior of local quantities in a trap. Two quantities that provide insight into the properties of the trapped system are the local densities and compressibilities [17]. In Fig. 7, we show the local radial densities ρ -averaged over the azimuthal angle -and the local compressibilities ∂ρ/∂µ i both as a function ofr i . The local compressibility was computed as a numerical derivative of the radial density when viewed as a function of the local chemical potential, in much the same way as it would be derived in an actual experiment. This is done for two values of chemical potential: The top panels correspond toρ match to the QMC results, although as is evident from the insets of the density profile panels, in the vicinity of the transition between the superfluid and the insulating phase, the LDA method shows larger deviations. This is also reflected in the compressibilities.

C. Local phases
The concept of characteristic density can also provide a general picture of the spatially-separated local phases inside the trap as the total filling and/or trapping potential are varied. To show that, we have computed the inner and outer radiir in andr out which enclose the superfluid region of the bosonic cloud. This was done by inspection of the local density profiles obtained via the QMC simulations for the different chemical potentials. In practice, the inner radiusr in has been determined as the radius inside of which the local density is greater than 0.999 whereasr out was defined as the radius outside of which the local density is lower than 0.001. We have verified that our results are robust against small deviations from the above values. Figure 8 shows the resulting diagram: Both radii are plotted for two different systems (two different trapping potentials) but with the same characteristic density at each point in the diagram. The circles indicate the SSE results with V = 0.01 while squares correspond to V = 0.02. The lines in the figure indicate the predictions of the different approximation schemes. This figure too illustrates that the concept of characteristic density is justified even when the LDA is not accurate, as the two systems corresponding to the two values of trapping amplitudes trace virtually the same curve. Also, remarkably the various approximation schemes but the LDA yield very accurate predictions of both radii.
Within the LDA, the two radiir in andr out , at which the local densities are n i = 1 and n i = 0, respectively, can be obtained in a straightforward manner: From the homogeneous results discussed in the previous sections, we know that these densities are reached at µ/(2dt) = ±1. Translating back to the harmonic trap case, this condition becomes µ/(2dt) −r 2 i /(2d) = ±1. Solving for r i for the two cases, we arrive at: The two radii are indicated by the dot-dashed lines plotted in Fig. 8. As one can see, the LDA yields a rather poor prediction of the two radii. This is expected however since these radii mark the transition from the superfluid phase to the insulating phase in the periodic system where diverging correlations are present.

D. Thermodynamic properties and the momentum distribution function
To conclude the analysis of the two dimensional case, here we present measurements of the kinetic and total (free) energies of the system. These are shown in Fig. 9 as a function of the characteristic density. The SSE results are denoted by circles, whereas the dotted, dashed and solid lines indicate the Gutzwiller-mean-field, the clustermean-field, and the spin-wave-corrected approximations. Both the kinetic energy and the inset in the free-energy panel, indicate that spin-wave-corrected results provide a far better match than the two mean-field methods.  Next, we analyze the momentum distribution of the model n(k) as a function of the radial-momentum coordinate k r = k 2 x + k 2 y (where, as with the radial densities, the azimuthal angle is averaged over). This quantity can be directly probed in experiments with ultracold atomic gases via absorption imaging after time-of-flight [1][2][3][4][5]. In the homogeneous case, due to the long-range decay of one-particle correlations in the compressible phase, the momentum distribution function has a delta peak singularity at k r = 0. In general, in the insulating phase, the off-diagonal one-particle correlations decay exponentially (they are zero in our hardcore model), yielding a broad momentum distribution (completely flat in our case). When a harmonic trap is present, the compressible and incompressible insulating phases coexist, so the zero-momentum peak is a smoother function of k r . for which such a phase exists. Interestingly, it is the latter value for which the zero-momentum peak is higher, despite the presence of the insulating phase in the center of the trap. This is because the system corresponding to the higher characteristic density has a higher number of bosons in the superfluid domain that surrounds the Mott insulator.
Looking at the Gutzwiller-and cluster-mean-field approximations, they both provide rather accurate predictions away from zero momentum but overestimate the peak by ≈ 16%. The spin-wave corrected prediction is much better in both cases, giving only a ≈ 5% error in the worst case.

V. TRAPPED HARDCORE BOSONS IN THREE DIMENSIONS
In the previous section, we showed that the Gutzwillermean-field solution already provides a fairly good description of the hardcore Bose-Hubbard model in the presence of a harmonic potential, where the deviations from the QMC data were shown to be rather small in most cases. This was true both for global observables such as the energies shown in Fig. 9, and also for local observables such as local densities and compressibilities shown in Fig. 7. When spin-wave corrections were taken into account, the results were improved to a point where most errors dropped to below 1%, yielding virtually exact results.
The encouraging results of the previous sections suggest that in three dimensions the spin-wave corrected results should be an even better approximation of the exact solution. In light of this, in what follows we explore the harmonically trapped three-dimensional hardcore Bose-Hubbard model using the Gutzwiller-meanfield approach with and without the addition of spinwave corrections. The three-dimensional harmonicallyconfined system is a good example of a system that becomes very demanding computationally if one wants to study the ground state for lattice sizes that are on the same order of magnitude as the ones realized experimentally. Having a large linear system size is critical when probing local observables that change with distance as is the case of the harmonic trap.
Below we present the results of the three-dimensional analysis: The Gutzwiller-mean-field (dotted lines), the spin-wave-corrected (solid lines) and the LDA (dotdashed lines) predictions of the equation of state (Fig.  11), the boundaries of the local phases in the trap (Fig.  12) and the energies (Fig. 13) of harmonically trapped hardcore bosons in three dimensions (for a 24 × 24 × 24 lattice). These results can be used in current experiments to understand density profiles and other properties of the   As all three figures indicate, in three dimensions the spin-wave corrections are less pronounced than in two dimensions. However, they still provide a discernible improvement over the mean-field results, particularly for the kinetic energy.
In Fig. 14, we examine some of the local properties of the three-dimensional system, specifically, the radial density profiles ρ(r) and the local compressibilities, as they are predicted by the Gutzwiller-mean-field, the spinwave-corrected and the LDA methods. The figure shows the results for two values of chemical potential: The top panels correspond toρ ≈ 23.90 [µ/(2dt) = 0.3] for which there is no Mott-insulating region in the center of the trap, while the bottom panels correspond toρ ≈ 118.12 [µ/(2dt) = 1.5], in which case there is an insulating region in the center, as indicated by the plateau in the density profile whenr is small. As is evident from the figures, in the vicinity of the transition between the superfluid phase and the surrounding insulating phases which are rather abrupt, the LDA results detach from the meanfield and spin-wave ones. This is also reflected in the compressibilities. there is an insulating region in the center. The Gutzwillermean-field results are indicated by the dotted lines, whereas the dashed, solid, and dot-dashed lines correspond to the Gutzwiller-mean-field, the cluster-mean-field, the Gutzwillermean-field with spin-wave corrections, and the LDA results, respectively.

VI. SUMMARY AND CONCLUSIONS
In this paper, we have studied the ground-state properties of homogeneous and trapped inhomogeneous systems of hardcore bosons in two-and three-dimensional lattices. We have obtained results both by means of quantum Monte Carlo (SSE) simulations, and via several analytical approximation schemes, such as the Gutzwillermean-field approach, a cluster-mean-field method, and a spin-wave analysis which takes quantum fluctuations into account. We studied the equation of state and various quantities of interest such as densities and compressibilities, condensate and superfluid densities, as well as the momentum distribution function and other basic thermodynamic properties.
In the homogeneous case, we showed that the equation of state can be predicted reasonably-well by all approximation schemes. However, for other quantities such as the kinetic energy, and the superfluid and condensate densities, the Gutzwiller and cluster mean-field predictions were found to deviate considerably from the quantum Monte Carlo results, particularly close to half filling. Interestingly, the addition of quantum fluctuations through spin-wave corrections was shown to yield virtually exact results for all quantities but the compressibility.
In the inhomogeneous case, we found that for the equation of state and some quantities, such as the onsite densities, even the simplest method, i.e., the Gutzwillermean-field approach, is capable of providing a fair analytical description. For other quantities, on the other hand, Gutzwiller-mean-field theory is not a very good approximation. For example, for the momentum distribution function we found that it largely overestimates the very low momenta occupations. A cluster-mean-field approach, based on larger cell sizes, did not prove to be much of an improvement over the Gutzwiller-mean-field, especially considering the large amount of additional degrees of freedom in the ansatz.
The addition of spin-wave corrections turned out to be an excellent approximation for all quantities, yielding virtually exact results when compared against the QMC data. This, in turn, allowed us to study large systems of three-dimensional harmonically-trapped hardcore bosons of sizes that would be very demanding for QMC techniques to simulate, with a fraction of the effort. Hence, this analytical approach can be used to study trapped systems in which the onsite repulsion is much larger than the bandwidth and the filling-per-site that is always lower or equal to one.
Finally, we examined the validity of the local density approximation in these harmonically trapped systems. We showed that it provides a fairly good description of local quantities as long as the measurement takes place away from the superfluid-insulator boundaries, which correspond to a phase transition in the homogeneous case.