Two-state Bose-Hubbard model in the hard-core boson limit

Phase transition into the phase with Bose-Einstein (BE) condensate in the two-band Bose-Hubbard model with the particle hopping in the excited band only is investigated. Instability connected with such a transition (which appears at excitation energies $\delta<\lvert t_0' \rvert$, where $\lvert t_0' \rvert$ is the particle hopping parameter) is considered. The re-entrant behaviour of spinodales is revealed in the hard-core boson limit in the region of positive values of chemical potential. It is found that the order of the phase transition undergoes a change in this case and becomes the first one; the re-entrant transition into the normal phase does not take place in reality. First order phase transitions also exist at negative values of $\delta$ (under the condition $\delta>\delta_{\mathrm{crit}}\approx-0.12\lvert t_0' \rvert$). At $\mu<0$ the phase transition mostly remains to be of the second order. The behaviour of the BE-condensate order parameter is analyzed, the $(\Theta,\mu)$ and $(\lvert t_0' \rvert,\mu)$ phase diagrams are built and localizations of tricritical points are established. The conditions are found at which the separation on the normal phase and the phase with the BE condensate takes place.


Introduction
During the recent years Bose-Hubbard model (BHM) is proved to be a valuable tool in the theory of systems of strongly correlated particles. The model achieves a wide recognition due to a successful description of thermodynamics and dynamics of ultracold Bose atoms in optical lattices where a phase transition to the phase with the Bose-Einstein (BE) condensate (so-called Mott insulator (MI) -superfluid state (SF) transition) occurs at very low temperatures. Experimental evidence of BE condensation in optical lattices was found for the first time in works [1,2] while theoretical predictions of such an effect were given earlier [3]. Starting from the 90-ies of the past century, a series of papers were devoted to the theory of this phenomenon. Among the first key articles on the subject one should mention the work [4] where BHM was studied in the mean field approximation. The calculated therein phase diagrams demonstrate that in the simplest case (i.e., hopping of Bose particles in the presence of a single-site Hubbard repulsion) the MI-SF transition is of the second order. Moreover, it is supposed that particles reside in the ground state of local potential wells in the lattice. Forthcoming theoretical investigations in this field were performed with the use of various techniques, e.g., the random phase approximation (RPA) in the Green function method [5,6], a strong-coupling perturbation theory [7,8], the dynamical mean field theory (Bose-DMFT) [9,10] as well as quantum Monte-Carlo calculations [11,12] and other numerical methods.
The Bose-Hubbard model is also intensively used for a theoretical description of a wide range of phenomena: quantum delocalization of hydrogen atoms adsorbed on the surface of transition metals [13,14], quantum diffusion of light particles on the surface or in the bulk [15,16], thermodynamics of the impurity ion intercalation into semiconductors [17,18].
In the last mentioned applications, there is usually a restriction on the position occupation number (n i 1), which corresponds to the limit of an infinite Hubbard repulsion for the considered model. Such a model of hard-core ions (where particles are described by the Pauli statistics) is also known as the fundamental one for the investigation of a wide range of problems, e.g., superconductivity due to a local electron pairing [19] or ionic hopping in ionic (superionic) conductors [20,21].
The study of a quantum delocalization or diffusion reveals an important role of excited vibrational states of particles (ions) in localized (interstitial) positions with a much higher probability of ion hopping between them [15,22,23]. A similar issue of a possible BE condensation in the excited bands in optical lattices is also considered but the condition of their sufficient occupation due to the optical pumping (see, e.g. [24]) is imposed. An orbital degeneration of the excited p-state is accompanied by anisotropy of hopping parameters and causes the appearance of variously polarized bands in the one-particle spectrum. Such bands correspond by convention to different sorts (so-called "flavours") of bosons and their number correlates with the lattice dimensionality. In the framework of the necessary generalization of the Bose-Hubbard model, a possibility of the MI-SF transition to the phase with BE condensate in the pumping-induced quasi-equilibrium long-living state of the system has been established [25].
In the equilibrium case, the issue of BE condensation involving the excited states in the framework of ordinary Bose-Hubbard model was not considered in practice. The exception is the system of spin-1 bosons [26,27] where a hyperfine splitting gives rise to multiplets of local states resulting in closely-spaced excited levels. As demonstrated in [28,29], the MI-SF phase transition could be of the first order when a single-site spin interaction is of the "antiferromagnetic" type. A similar change of the phase transition order also takes place for multicomponent Bose systems in the optical lattices [30].
In the present work we consider an equilibrium thermodynamics of the Bose-Hubbard model taking into account only one nondegenerated excited state on the lattice site besides the ground one. On the one hand, such a model corresponds to 1D or strongly anisotropic (quasi-1D) optical lattice, and on the other hand, it is close to a situation that is characteristic of a system of light particles adsorbed on the metal surface. For example, the excited states of hydrogen atoms on the Ni(111) surface are sufficiently distant [22] so only the lowest one could be taken into account. We shall investigate a condition of instability of a normal state of the Bose system with respect to BE condensation considering a criterion of divergence of the susceptibility (χ ∼ c l |c + p ω | q=0,ω=0 ) characterizing the system response with respect to the field related to a spontaneous creation or annihilation of particles. We shall also study the behaviour of the order parameter c 0 ( c + 0 ) as well as the grand canonical potential in the region of the MI-SF transition and shall build relevant phase diagrams. Special attention will be paid to a change of the phase transition order and localization of tricritical points at different values of excitation energy, particle hopping parameter and temperature.
We shall limit ourselves to the hard-core boson (HCB) limit where a limitation on occupation numbers is present: no more than one particle per site regardless of the state (excited or ground) occupied by it. Thus, the single-site problem is a three-level one (contrary to the two-level ordinary HCB case). For this reason, it is convenient to use the formalism of Hubbard operators [31] (standard basis operators [32]).
where b i and b + i (c i and c + i ) are Bose operators of annihilation and creation of particles in the ground (excited) state, ε and ε ′ are respective energies of state and µ is the chemical potential of particles. Such a Hamiltonian includes the single-site Hubbard repulsions with energies U b , U c and U bc as well as the particle hopping between ground (t b ), excited (t c ) and different (t bc ) states. Hereinafter we assume U b = U c = U bc for simplicity.
Let us define a single-site basis |n b i , n c i ≡ |i; n b i , n c i (which is formed by particle occupation numbers in the ground and in the excited states, i.e., eigenvalues of operators In the X-operator representation, the single-site part of Hamiltonian (2.1) can be written aŝ Terms describing an inter-site transfer in Hamiltonian (2.1) are transformed in a similar way. Our primary goal is to calculate the two-time temperature boson Green's functions b|b + and c|c + , which describe an excitation spectrum and make it possible to investigate the conditions of the system's instability with respect to the spontaneous symmetry breaking and the appearance of a BE condensate. As follows from definitions (2.3) The latter are originated from the commutation of an initial X-operator with the inter-site transfer terms of the Hamiltonian, thus producing the higher-order Green's functions where X ··· l stands for operators on the right-hand side of expressions (2.10a)-(2.10d). Decoupling of functions (2.11) in the random phase approximation (RPA) is performed in the following way: In the case of the normal phase (which will be studied herein and equation (2.8) can be rewritten as 14) The following notations are introduced for the occupation difference of adjacent levels and the related transition energies when the number of Bose particles in the ground state (with the energy ε) on the site increases by one.
Proceeding from X-operators in equation (2.14) to the Bose operators b and b + according to definition (2.3) we obtain where the function has the meaning of the unperturbed Green's function for bosons residing in the single-site ground state. Equations of motion for "mixed" Green's functions c|b + are obtained in the way similar to the above described scheme. Using decoupling (2.12) one can write which results in the equation Here, similarly to (2.15) and (2.17) (2.20) and the function g ′ 0 (ω) is the unperturbed Green's function for bosons residing in the excited state. By means of the Fourier transform one can proceed to the momentum representation obtaining a system of equations where t q , t ′ q and t ′′ q stand for the Fourier transforms of hopping parameters. A pair of equations for other Green's functions are obtained in a similar way Solutions of equations (2.22) and (2.23) are as follows: The equation D q (ω) = 0 gives the excitation spectrum which is obtained here in the RPA. On the other hand, the divergence of boson Green's functions (2.24) at the zero values of wave vector and frequency is the criterion of instability with respect to BE condensation [5,33], thus giving the following condition which can be rewritten in the explicit form and δ = ε ′ − ε is the excitation energy.

13004-5
We should point out that divergence of the b|b + 0,0 function correlates with the appearance of the BE condensate in the ground state while at the divergence of the c|c + 0,0 function, BE condensation takes place in the excited state. In general, both condensates appear simultaneously except the case t ′′ q = 0 (e.g. due to symmetry reasons) where these effects become independent and only the one type of condensate arises in the instability point.
Equation (2.27), mutually relating the chemical potential, hopping parameters and temperature, allows us to construct spinodal surfaces (or lines) in the above mentioned coordinates and to find the temperature of the phase transition to the phase with BE condensate (so-called SF phase) where such a transition is of the second order. Below, this problem (especially the issue of the phase transition order) will be investigated more in detail.

NO phase instability in HCB limit
Let us consider now a simple special case of the HCB limit when occupation numbers in the |n, m state are restricted by a condition n + m 1. In the framework of the model, it formally means U → ∞.
In this case, the model becomes a three-level one with the local energies and the following transition energies Thus, equation (2.27) can be rewritten in the form where Q 00 = 1 − e βµ 1 + e βµ + e β(µ−δ) , Q ′ 00 = 1 − e β(µ−δ) 1 + e βµ + e β(µ−δ) (3.4) in the zero approximation with respect to hopping.  1. Lines of the NO phase instability (spinodals) with respect to the appearance of BE condensate in the (Θ, µ) plane in the HCB limit at various excitation energies (t0 = 0, We take into account (according to estimations made in [15,25]) that boson wave functions in adjacent potential wells overlap in greater extent in the excited states compared to the ground ones. Accordingly, we shall put here t 0 = 0. For a centrosymmetric lattice and in the case of different parity of wave functions of ground and excited states we have also t ′′ 0 = 0. Finally, we follow a usual convention of the BH model for optical lattices taking t ′ 0 < 0. In this way equation (3.3) can be reduced to (3.5) Its solutions determine the stability region boundaries of the normal (NO) phase. Respective lines of spinodals are numerically calculated and presented in figure 1 (here and below the energy quantities are given in units of |t ′ 0 |). As illustrated in figure 1, at δ < |t ′ 0 | spinodals surround an asymmetric area in the 13004-6 (Θ, µ) plane which is located between the points µ = δ − |t ′ 0 | and µ = 0 of the abscissa axis. In this region, the NO phase is unstable; this is connected with the appearance of BE condensate. At δ < |t ′ 0 |/2 and µ > 0 the backward path of spinodal is observed and a lower temperature of the NO phase instability appears, thus suggesting a possibility of the SF phase existence in the intermediate temperature range (so-called "re-entrant transition"). However, as will be shown further, in the mentioned region a real thermodynamic behaviour is even more complicated. The order of the NO-SF transition can change to the first one and the SF-phase remains stable up to the zero temperature.

Phase diagrams in MFA
For a more detailed treatment of the NO-SF transition issue, let us study the thermodynamics of the considered system in the HCB limit, thus reducing the problem to a three-state model with the Hamiltonian where the shorthand notations are used Possibility of BE condensation will be studied in the MFA. Average values of creation (annihilation) operators for Bose particles in the ground or excited local state play the role of order parameters for the SF-phase. Hence, the mean-field Hamiltonian is as follows: Limiting our consideration to the case of particle hopping only through excited states (t ′ 0 = 0, t 0 = t ′′ 0 = 0) we can diagonalize Hamiltonian (4.4) by a rotation transformation  where and λ 2 − λ 0 = δ − µ. In terms of operators X rs = |r s| New energies of single-site states arẽ In the new basis which yields after averaging the order parameter ξ: (4.12) Solution ξ = 0 corresponds to the NO phase. A nonzero solution describing the BE condensate is obtained from the equation (4.13) In the limit ξ → 0 this equation determines the line where the order parameter for the SF phase tends to zero. One can readily see that it coincides with spinodal equation (3.5) thus defining the line of the second order NO-SF phase transition (when just the transition of such an order takes place).
Numerical solutions of equation (4.13) make it possible to study the behavior of the order parameter ξ depending on chemical potential µ at various temperatures as illustrated in figure 2. In the main, at negative values of chemical potential the parameter ξ changes smoothly and the phase transition to the SF phase is of the second order. But at µ 0 and low enough temperatures, the  ξ(µ) dependence has an S-like bend. In this case, the first order phase transition with an abrupt change of the parameter ξ takes place. This phase transition occurs at a certain value of the chemical potential which could be calculated using the Maxwell rule or considering the minimum of the grand canonical potential Ω(µ) as a function of the chemical potential (see below). Obviously, the point of ξ nullification does not anymore correspond here to the phase transition. Similar behaviour of the parameter ξ holds even at zero excitation energy (δ = 0) where the first order phase transition remains for nonzero temperatures whereas at T = 0 its order changes to the second one ( figure 3). At negative values of δ (which corresponds to inversion of ε and ε ′ levels and to hopping between ground states) the second order of the transition is preserved in the low-temperature region close to T = 0 transforming to the first order transition at the temperature increase and recovering henceforth (  13004-9 significant (figure 5). In the case of δ < 0, two critical points appear at a certain distance; the latter tends to zero at δ = δ crit ≈ −0.12|t ′ 0 | and the first order phase transitions at a further increase of δ (figure 6) is suppressed.
Phase diagrams in the (|t ′ 0 |, µ) plane at various temperatures for δ > 0 are depicted in figure 7 with indication of tricritical points. In distinction to the standard two-level HCB model [34] (where the SF phase transition is of the second order) the diagrams are asymmetric. In the limit T → 0 for µ > 0 the first order transition occurs at µ = ( √ δ − |t ′ 0 |) 2 (see the next section) whereas for µ < 0 they are of the second order on the line µ = δ − |t ′ 0 |.

Phase separation at fixed boson concentration
Let us consider now the thermodynamics of the model at a fixed concentration of Bose particles. We will utilize a connection between the concentration and the chemical potential of bosons which can be established using its definition in such a form or basing on the relationship In the first case similarly to equality (4.10) one can obtain a relation which results in

(5.4)
In the second case, taking into account that and differentiating with respect to µ, one can come to the same expression as (5.4).
There are different relationships between n and µ in NO and SF phases; in the last case, a nonzero value of ξ (a solution of equation (4.13)) should be substituted into expression (5.4). Order parameter ξ has a jump at the first order phase transition, so a stepwise change of concentration n takes place. In the n = const regime (at the value of n in the region of step) it means a phase separation into two phases with different concentrations: the NO phase (ξ = 0 and a larger concentration of bosons) and the SF phase (ξ = 0 and their smaller concentration).
The above described situation is illustrated in figure 8, where the numerically calculated (Θ, n) phase diagrams are presented. At δ > 0, phase separation region spans up to tricritical temperatures. When δ goes to zero and finally reverses its sign, the shape of the separation region changes in a peculiar way moving off abscissa axis (figure 8). Now the phase separation begins at nonzero temperatures and vanishes at δ < δ crit ; the line of the second order phase transition remains only. At the further increase of |δ| (in the µ < 0 region) the (Θ, n) diagram becomes more and more symmetric, approaching by its shape the diagram known for the usual HCB model [35] (see also [36]).
Phase diagrams in the (|t ′ 0 |, n) coordinates are given in figure 9 where the regions of NO, SF and separated phases are shown at various temperatures. The case of the zero temperature can be studied more in detail in a pure analytic way. In this limit there are three branches of order parameter ξ as a function of the chemical potential (see figure 2): After elimination of ξ parameter, one can obtain the grand canonical potential Ω as follows: : : Ω/N = 0, µ < 0, −µ, µ > 0.

13004-11
Differentiating expressions (5.7) with respect to µ we have At the first order phase transition from the SF phase to the NO phase, the order parameter ξ jumps from branch (1) to branch (3). This occurs at the µ = µ * ≡ ( |t ′ 0 | − √ δ) 2 value given by equality of respective grand canonical potentials Ω (1) = Ω (3) . Then boson system separates into SF and NO phases with concentrations of bosons:

Discussion and conclusions
As was shown in this work, the transition to the SF phase (the phase with BE condensate) in the Bose-Hubbard model with two local states (the ground and excited ones) on the lattice site can be of the first order in the case, when the particle hopping takes place only in the excited band. Calculations and estimates for optical lattices give evidence of significant distinction between hopping parameters t 0 and t 1 in the ground and excited bands, respectively. It follows from estimates [25] that t 1 /t 0 ≈ 30 − 50 depending on depth V 0 of local potential wells (one can produce effect on V 0 changing the intensity of laser beams which create an optical lattice). Similar results are obtained in the studies of quantum delocalization of the adsorbed hydrogen atoms. One can see from calculations [22,23] of energy spectrum of the H-atom subsystem on the nickel surface that the ground-state band has a negligible bandwidth. At the same time, for excited bands, the bandwidth varies in the range from 15 to 45 meV (depending on the excited state symmetry and on the crystallographic orientation of metal surface), being mostly of the order of half the corresponding excitation energy ∆ε α = ε α − ε 0 . There are, however, the cases of strong delocalization (e.g. H on the Ni(110) surface) where the excited bands overlap, and the width of the lowest one is of the same order as ∆ε α [22].
The values of hopping parameters greatly increase at the decrease of V 0 ; the distance between the local energy levels becomes smaller in this case (see [37,38]). It is one of the possible ways of changing the relation between the hopping parameters and excitation energy (|t ′ 0 | and δ in our model). Another possibility (discussed in [39]) is connected with an essential reduction of the energy gap between local s-and p-levels due to sufficiently strong interspecies Feshbach resonance in the presence of Fermi atoms added to the Bose system in optical lattice.
Along with investigation of BE condensation in the excited band (or bands p x , p y (p x , p y , p z ) in two-(three-) dimensional case) on condition that certain concentration of Bose-atoms has been created in the band by optical pumping [25,38], an attempt was made in [40] to study the effect of excited bands on the physics of BE condensation in the lowest (s-) band (when the s-band hopping is taken into account). The case of finite values of the one-site interaction U was considered. The possibility of the re-entrant behaviour of the MI-SF transition was claimed. However, the order of phase transition was not investigated; the consideration was restricted to the case of zero temperature. As we show in this work, re-entrant type dependence on T or µ takes place only for spinodals and the return to the initial MI phase from the SF phase could be possible only in the case of the second order phase transitions. In reality, the order of phase transition changes to the first order in this region. In the HCB limit (no more than one particle per lattice site), it takes place mainly at positive values of chemical potential of particles; at µ < 0, the transition remains, for the most part, of the second order. The region of existence of SF phase is restricted, as a whole, to the interval −|t ′ 0 | < µ < |t ′ 0 |, while excitation energy should obey the inequality δ < |t ′ 0 |. We have constructed the corresponding phase diagrams and established localization of tricritical points, where the order of phase transition changes. The separation on SF and NO phases at the fixed particle concentration is investigated; the conditions of the appearance of phase-separated state are analyzed.
It should be mentioned that phase diagrams in figures 2-9 are close by their shape to the diagrams obtained in the framework of Bose-Hubbard model for Bose atoms with spin S = 1 in optical lattices [29]. The excited levels are formed in that case by the higher spin single-site states and corresponding interactions of the "ferromagnetic" or "antiferromagnetic" type (the Hund-rulelike splitting), while the hopping parameter is taken the same for all bands. The similarity of the mentioned diagrams points out to the fact that the role of the excited states in the change of the phase transition order in going to the phase with the BE condensate is the same in both cases. Distinction, however, consists in another genesis of the single-site spectrum. In our model, in the limiting case of HCB there are no effects connected with the level splitting due to interaction; the excited single-particle states are taken by us into account instead.
The consideration developed in this work can be extended to the systems with the close or degenerate excited local levels. Generalization of the model by adding inter-site interactions is also important. It could even make it possible to take into consideration other phases (densitymodulated or supersolid) besides NO and SF ones.
We finally emphasize that the hopping parameter t ′ ij in the excited band can be positive; in particular, this concerns the p-bands [39]. In such a situation, the condensation takes place into states with wave vector Q on the boundary of the Brillouin zone, while the order parameters c Q , c + Q describe the modulated condensate. Since t ′ Q = −t ′ 0 , the results obtained in this work are also valid (with |t ′ Q | in place of |t ′ 0 |) in that case.