The two-state Bose-Hubbard model in the hard-core boson limit: Non-ergodicity and the Bose-Einstein condensation

The Bose-Einstein condensation in the hard-core boson limit (HCB) of the Bose-Hubbard model with two local states and the particle hopping in the excited band only is investigated. For the purpose of considering the non-ergodicity, a single-particle spectral density is calculated in the random phase approximation by means of the temperature boson Green functions. The non-ergodic contribution to the momentum distribution function of particles (connected with the static density fluctuations) increases significantly and becomes comparable with the ergodic contribution in the superfluid phase near the tricritical point.


Introduction
The Bose-Hubbard model (BHM) [1] is widely used in describing the thermodynamics and the dynamics of Bose atoms in optical lattices [2,3]. Such systems demonstrate the phase transition into the phase with Bose-Einstein (BE) condensate at very low temperatures. Thus, the system can be in the normal (NO) phase (the state of the so-called Mott insulator, MI), or in the phase with the BE condensate (the superfluid state, SF). The BHM can be used for description of various phenomena such as the quantum delocalization of hydrogen atoms adsorbed on the surface of transition metals [4], quantum diffusion of light particles on the surface or in the bulk [5], thermodynamics of the impurity ion intercalation into semiconductors [6].
The BHM has been intensively studied for the two recent decades starting with the pioneering work [7], where the thermodynamics of the model was considered in the mean field approximation (MFA) for a particle hopping with an exact treatment of the Hubbard repulsion. Such a description revealed the first order NO-SF phase transition for particles residing in the ground state of local wells. A prospective generalization of the BHM consists in taking into account the excited vibrational states of bosons in local potential minima of lattice sites. For example, intersite particle hopping through the excited states is much easier in the case of quantum delocalization or diffusion [8,9] as well as in the case of optical lattices [10]. However, a possible BE condensation in the excited states was considered for the case of optical pumping in order to increase their occupations [11]. In this case, the orbital degeneration of the excited p-state can lead to the appearance of a special type of the BE condensate.
A thorough study of the phase transition into the phase with BE condensate in the BHM with two local states and the boson hopping only in the excited band was carried out in works [12,13]. In the MFA and the hard-core boson (HCB) limit, the instability related to the NO-SF transition was studied which occurs at excitation energies lower than the absolute value of the hopping parameter. Conditions of the change of the above mentioned phase transition from the second order to the first one were studied. The respective phase diagrams (Θ, µ) and (|t ′ 0 |, µ) were analyzed. A possible phase separation into NO and SF phases at a fixed concentration was revealed. Two-time boson Green's functions (GF) and a single-particle spectral density were calculated in the random phase approximation (RPA). The structure of excitation spectra of the particle and hole types in the region of the NO-SF phase transition was investigated. Some physical many-particle systems do not hold the ergodicity hypotesis. It means that there are regions of the phase space of a system which cannot be reached by the trajectory of the system evolution. From practical viewpoint, non-ergodicity leads to a distinction between isothermal and isolated responses (susceptibilities) of a system caused by the zero-frequency term present in the isothermal response only. In this case, the use of the "isolated" (Kubo) response leads to erroneous values of such physical characteristics as compressibility, specific heat and susceptibility for the systems described by means of the Ising model [14,15] (a zero Kubo response), the pseudospin-electron model [16] and the Falicov-Kimball model (the non-ergodic contribution describes the singularity related to the phase transition [17,18]). The BHM in the HCB limit, as will be shown below, also belongs to the non-ergodic class. The mentioned problems can be solved using the temperature (Matsubara's) GF [19] (see, e.g., a scheme for calculation of many-time correlation functions [20]).
In the present article we study the non-ergodic contribution to the particle momentum distribution function by the example of the two-state Bose-Hubbard model in the HCB limit. A detailed investigation is performed in the region close to the tricritical point where the order of the phase transition to the SF phase changes to the first one.

Thermodynamics of the model in the hard-core boson limit
The Hamiltonian of the BHM takes into account the tunnel hopping of particles between the nearest neighbour lattice sites and the mutual repulsion of the particles located in the same potential well [1]: Here, t i j characterizes the hopping energy, U is the energy of the Hubbard intrasite interaction (U > 0), µ is the chemical potential of the boson particles, b + i , b j are the Bose operators of particle creation and annihilation on the i -th site in the ground vibrational state. Taking into account the first excited state and considering the particle hopping only between these excited states, one can generalize the Hamiltonian (2.1) as follows:Ĥ where c + i , c i are the Bose operators of the particles in the excited state, ε (ε ′ ) is the energy of the particle in the ground (excited) state, U b ,U c ,U bc are the parameters of the Hubbard repulsion.
Starting from the basis of the single-site states |i ; n b i , n c i 〉, formed by the particle occupation numbers (eigenvalues of the operators Let us restrict our study to the HCB case (when n + m 1). In the considered model, such a situation is reached at U b ,U c ,U bc → ∞. In this limit, the single-site Hamiltonian becomes three-level with the respective energies λ 0 = 0, λ 1 = −µ, λ 2 = δ − µ (the following shortened designations for single-site states are used: |0〉 ≡ |00〉, |1〉 ≡ |10〉, |2〉 ≡ |01〉). Here δ = ε ′ − ε is the energy of transition to the local excited vibrational state. Respectively,

33002-2
The BE condensation takes place in the band that originated from the hopping between the excited states of the neighbouring wells (see below). Thus, the order parameter is equal to the average of boson creation or annihilation operators ξ = 〈X 20 Let us separate a MFA part of the Hamiltonian (2.4) where t ′ 0 is the Fourier transform of the hopping energy t ′ i j at q = 0 (hereinafter we consider t ′ 0 < 0, i.e., assuming symmetrical excited states, while t ′ 0 > 0 for antisymmetrical ones [22]). The ξ parameter is obtained from the self-consistency equation This approximation corresponds to the exact solution for the case of an infinite-range hopping. A strict derivation of the above statement for the standard BHM has been given in the work [23] (see also [24]) based on the approach in the style of the Bogolyubov-Tyablikov variational method.
Finally, a complete Hamiltonian looks as follows: The MF HamiltonianĤ MF can be reduced to a diagonal form by the transformation (2.8) In terms of operators X µν = |μ〉〈ν|Ĥ Energies of single-site statesλ µ in the phase with a broken symmetry (ξ 0) are equal toλ 0,2 = 1 2 (δ−µ∓E ), It is convenient to introduce the following linear combinations Thus, the X -operators are represented on a new basis as follows:  The Hamiltonian (2.6) expressed in terms of operators (2.10) is a three-state generalization of the standard Hamiltonian of the HCB model [25]. The presence of the third level |1〉 does not effect the pseudospin dynamics but modifies the occupations of levels |0〉 and |2〉, changing this way the thermodynamic behaviour of the model. Moreover, the application of pseudospin formalism allows one to utilize the usual scheme of the random phase approximation (RPA) for Green's functions of spin models, originated from the well-known Tyablikov decoupling for the Heisenberg model [26].
Let us make a short revision of the thermodynamics of a three-state model in the MFA. Using the Gibbs distribution with the Hamiltonian (2.9) and taking into account that the average values of transverse pseudospin projections 〈σ ± i 〉 are equal to zero, one can obtain the equation for the order parameter ξ = 〈X 02 i 〉: (2.13) The case of ξ = 0 corresponds to the NO phase. One can obtain nonzero solutions (indicating the appearance of the BE condensate) from the following equation (2.14) It has been shown in work [12] that in the region µ < 0, nonzero values of ξ gradually appear at the second order transition while in the region µ > 0 at a low enough temperature, the curve ξ(µ) is Unlike the standard (two-level) HCB model (where the Bose-particles remain in the local ground state, the excited local states are neglected and phase transitions between NO and SF phases are always of the second order), the phase diagrams in figure 1 are asymmetrical.

Temperature Green's functions and the non-ergodic contribution
The spectrum of boson particles and the single-site spectral density for this model are calculated in the RPA using the two-time GF in work [13]. Herein below we will perform similar calculations in the RPA but applying the temperature GF.
Assuming the hopping between the excited states only, the dynamics of the model corresponds to the HCB model and we consider only GF constructed on the operators which describe the transition between the states |0〉 and |2〉 T τ X 20 X 02 q,ω n = βξ 2 δ(q)δ(ω n ) + T τ (X 20 − ξ)(X 02 − ξ) q,ω n . (3.1) Using the above mentioned "pseudospin" representation (2.10), one can express (3.1) in the form T τ (X 20 − ξ)(X 02 − ξ) q,ω n = G zz sin 2 2ϑ + 1 4 G ++ (cos 2 2ϑ − 1) Temperature GF in the RPA are derived by means of the Larkin equation where unperturbed GF are expressed through the single-site GF and the average values of the pseudospin operator and the parameters of the effective interaction look as follows: The equation of the type (3.4) for the ordinary HCB model was derived in work [27]. The momentum distribution of particles in the excited state is obtained as the sum over Matsubara's frequencies from the temperature GF An explicit form of the spectrum depends on the phase where the system resides. For example, for NO phase (3.11) while there are two branches ±ε (SF) q for the SF phase This result formally coincides with its analogue for ordinary (two-state) HCB model (which can be reached in the limit δ → −∞). However, an important difference lies in the fact that the values of 〈σ z 〉 and ϑ are defined by an equation for the three-state system (2.12).
There are three parts in the momentum distribution: − cos 2ϑ single-particle excitations non-ergodic part . (3.13) The first and the second terms correspond to the results obtained in work [13] while the non-ergodic part (related to the static density fluctuations) can be obtained only by means of temperature GF: (3.14) The factor sin 2ϑ is non-zero in the SF phase only, so the existence of the non-ergodic part is also limited to this phase: Due to a sophisticated structure of the nonergodic part, a thorough analysis of its behavior is rather complicated. Dependencies of 〈σ z 〉 and b ′ parameters, entering the formula (3.14), on chemical potential and temperature are presented in figure 2. As one can conclude, the non-ergodic contribution is "frozen out" at low temperatures (because the hopping is possible between the excited states only). The momentum distribution, including according (3.13) all three parts, rapidly increases approaching the centre of the Brillouin zone (BE-condensation). The relative role of the non-ergodic part in the total momentum distribution is illustrated in figures 3-5. In the SF phase E = 2|t ′ 0 |〈σ z 〉 0 , therefore

33002-7
latter are characterized by the derivative ∂n/∂µ, that reduces to the b ′ parameter in the simplest approximation). It is confirmed by the behaviour of the n(µ) dependence at different temperatures obtained and illustrated in [13]. The total occupation of the excited state can be calculated by the summing over the wave vector (here we substitute the summation by the integration with the model density of states g 0 (z) for a simple cubic lattice [13]) n c RPA ≡ n = 1 N q n q = n z g 0 (z) dz. (3.15) It is interesting to compare partial contributions of the BE-condensate, the single-particle excitations and the non-ergodic part in the total occupation of the excited state: The most pronounced non-ergodic contribution is visible near the line of the second order phase transition in the vicinity of the tricritical point ( figure 6).

Conclusions
The non-ergodic part of the momentum distribution (related to the static density fluctuations) in the excited state for the BHM in the HCB limit is calculated by means of the temperature GF. This part becomes nonzero only in the SF phase. Its value rapidly grows and becomes on par with the ergodic part in the SF phase near the tricritical point (where ∂n/∂µ → ∞).
The questions remain, is the non-ergodicity (and the presence of corresponding contributions to the particle momentum distribution) a characteristic feature of the hard-core boson limit only and what could be at the transition to the standard BHM. Such a transition takes place when a constraint of occupation number n i is removed and the on-site repulsion energy U becomes finite. When the non-ergodicity would disappear in this case, it could mean that the effect is caused by the change of the particle statistics (due to the transition from Pauli-to Bose-operators). This problem needs a separate and more detailed consideration.
Dynamical characteristics of the model obtained in the present study and in the work [13] are generally based on the average values calculated in the MFA. There is a slight discrepancy between these characteristics and the respective results in RPA. For a complete self-consistency, one should use the appropriate expressions for averages, starting with the grand canonical potential obtained by summation of diagram series consisting closed cycles.