Self-consistent approach for Bose-condensed atoms in optical lattices

Bose atoms in optical lattices are considered at low temperatures and weak interactions, when Bose-Einstein condensate is formed. A self-consistent approach, based on the use of a representative statistical ensemble, is employed, ensuring a gapless spectrum of collective excitations and the validity of conservation laws. In order to show that the approach is applicable to both weak and tight binding, the problem is treated in the Bloch as well as in the Wannier representations. Both these ways result in similar expressions that are compared for the self-consistent Hartree-Fock-Bogolubov approximation. A convenient general formula for the superfluid fraction of atoms in an optical lattice is derived.


Introduction
Systems with Bose-Einstein condensate are interesting objects from both theoretical and experimental points of view. That is why they have been intensively studied in recent years. Vast literature on this problem can be found in the books [1][2][3][4] and review articles [5][6][7][8][9][10][11][12][13][14][15]. Creation of optical lattices has made it possible to achieve a new dimension in the physics of cold atoms, providing an opportunity for numerous novel applications and for modeling many effects typical of condensed matter [16][17][18][19].
The occurrence of Bose-Einstein condensate is intimately related to the global gauge symmetry breaking [2,11] that is a necessary and sufficient condition for Bose-Einstein condensation. In the theory of Bose-condensed systems, there exists an old problem, formulated by Hohenberg and Martin [20], who showed that, as soon as gauge symmetry is broken, the description of such a system suffers from one of the defects, either yielding unphysical spectrum of collective excitations or resulting in broken conservation laws and incorrect thermodynamics. Any of these deficiencies implies that the description is not self-consistent, corresponding to an unstable system. This problem has been solved by employing representative statistical ensembles [21][22][23] to systems with a broken gauge symmetry [24][25][26][27]. This approach was shown to be completely self-consistent and gapless, with the Hartree-Fock-Bogolubov (HFB) approximation [28,29] providing an accurate description for uniform Bose systems [27,[30][31][32], as well as for these systems in random external potentials [33,34].
In the present paper, this self-consistent approach is applied to Bose-condensed atoms in optical lattices. Sections 2 and 3, contain the main definitions related to optical lattices and Bose-condensed atoms, respectively. In section 4, the Bloch representation is used, which can be more suitable for weak binding, while in section 5, the Wannier representation is employed, which is more convenient for tight binding. Both these cases are treated in the HFB approximation leading to similar results. However, the Wannier representation, yielding the Hubbard Hamiltonian, is a bit simpler. Some thermodynamic characteristics are considered in section 6, where a general and convenient formula for superfluid fraction is derived. Section 7 concludes.
Throughout the paper, the system of units is used, where the Planck and Boltzmann constants are set to be one.

Optical lattices
Optical lattices are created by laser beams forming standing waves, which corresponds to the formation of a periodic lattice potential with a being a lattice vector with the components a α = λ α /2, where λ α is a laser wavelength and α = 1, 2, . . . , d enumerates spatial components in a d-dimensional space. The standard form of the lattice potential is with the laser wave vector The lattice depth is defined by the parameter Another important quantity, characterizing an optical lattice, is the recoil energy where m is atomic mass. The ratio E R /V 0 characterizes the relative lattice depth.

Bose atoms
The lattice is loaded with Bose atoms, whose interactions are measured by means of the scattering length a s entering the effective interaction strength The energy operator is given by the Hamiltonian in which U = U (r) is a trapping potential, if any, and V L = V L (r) is a lattice potential. The atom field operatorsψ(r) satisfy the Bose commutation relations. The existence of Bose-Einstein condensate necessarily requires that global gauge symmetry should be broken [2,11]. The most straightforward way of the gauge symmetry breaking is by means of the Bogolubov shift of the field operatorψ (r) = η(r) + ψ 1 (r) . Here, the first term is the condensate wave function normalized to the number of condensed atoms N 0 = |η(r)| 2 dr .

23002-2
The second term is the field operator of uncondensed atoms, whose number is given by the statistical average (3.5) of the number-of-particle operatorN 1 . The uncondensed atoms are normal in the sense that the average of their field operator is zero, To avoid double counting of the degrees of freedom, the orthogonality condition η * (r)ψ 1 (r) dr = 0 is required. This condition is a direct consequence of orthogonality of wave functions serving as a basis for the expansion of the field operatorψ(r) [28,29]. The number of atoms per lattice site is called a filling factor that is defined as the ratio in which a is a mean interatomic distance and ρ is the average atomic density, The representative ensemble for a system with a broken gauge symmetry is characterized [24][25][26][27] by the grand Hamiltonian It is worth stressing that the introduction of two Lagrange multipliers, µ 0 and µ 1 is necessary due to the presence of two independent variables in the Bogolubov shift (3.3) and the related two normalization conditions (3.4) and (3.5). It is a general mathematical fact that the number of Lagrange multipliers should be equal to the number of imposed constraints, such as the normalization conditions. The theory can become non-self-consistent if the number of Lagrange multipliers is smaller than that of the imposed constraints. Introducing two Lagrange multipliers does not exclude that in particular cases, these multipliers could become equal, as it happens in the Bogolubov approximation [28,29]. The physical meaning of using two Lagrange multipliers has been thoroughly explained in the previous papers [11,14,19,[23][24][25][26][27][30][31][32][33].

Bloch representation
One usually considers optical lattices by reducing the problem to a Hubbard Hamiltonian by means of the Wannier representation which is convenient in the case of a tight binding. Here, we show that it is equivalently possible to employ the Bloch representation that can be more appropriate for weak binding and leads to the results similar to those in the Wannier representation to be considered in the following section. Below, we assume that there is no trapping potential, so that the system is ideally periodic.
Let {ϕ nk (r)} be the basis of Bloch functions labeled by the zone index n and quasi-momentum multiindex k. Then, the field operators of uncondensed atoms can be expanded over this basis, ψ 1 (r) = nk a nk ϕ nk (r) .

23002-3
The basis should be chosen so that the Bloch functions are natural orbitals [35], that is, the eigenfunctions of the density matrix Then, the density matrix enjoys a diagonal expansion In other words, the use of natural orbitals simplifies the consideration due to the following properties Here, the first term contains only a condensate wave function, but no field operators of uncondensed atoms. The term, linear in ψ 1 , is canceled by the Lagrange termΛ. In the following expressions, the pair {n, k}, for brevity, will be denoted as k, while the set {n, −k}, as −k. Then, the term, containing the products of two operators of uncondensed atoms, reads as The term of third order, with respect to the products of the field operators of uncondensed atoms, is And the fourth-order term is In the Hartree-Fock-Bogolubov (HFB) approximation, the third-order term H (3) yields expressions linear in ψ 1 , which should be canceled by the Lagrange cancelerΛ. The fourth-order part takes the form in which the notations for the so-called normal n k ≡ 〈a † k a k 〉 , (4.11) and anomalous averages are used. The normal average (4.11) is the distribution of atoms, while the absolute value |σ k | of the anomalous average (4.12) is the distribution of the correlated atomic pairs [19,25,29]. Let us introduce the notation (4.14) Then, the grand Hamiltonian (4.5) in the HFB approximation can be written as where the first term is the nonoperator quantity The quadratic Hamiltonian (4.15) can be diagonalized and all observables calculated. However, the resulting expressions are rather complicated. In order to simplify the calculations, it is possible to assume that the main contribution in the above formulas comes from diagonal terms, since the Bloch functions are mutually orthogonal. This can be referred to as the diagonal approximation, when expressions (4.13) and (4.14) take the form The use of the diagonal approximation is not compulsory and it is possible to diagonalize the quadratic form (4.15) without it. This approximation, however, essentially simplifies the formulas. Justification of this approximation is based on the fact that the expansion functions ϕ k are mutually orthogonal, which makes it reasonable to assume that the matrix elements over these functions are such that their diagonal elements are larger than off-diagonal.
In the diagonal approximation, Hamiltonian (4.15) reduces to This form is much simpler to diagonalize using the Bogolubov canonical transformation [28,29].

23002-5
Following a standard procedure by diagonalizing Hamiltonian (4.20), we find the Bogolubov spectrum of elementary excitations The condition of condensate existence [14,19] requires that the spectrum should be gapless, This condition is equivalent to the Hugenholtz-Pines theorem [36]. Hence, we get where the notations are used for the condensate density density of uncondensed atoms 25) and the anomalous average The equation for the condensate wave function, in the case of an equilibrium system, is defined by the variational condition δH δη * (r) = 0, (4.27) which yields the equation (4.28) The latter gives the condensate chemical potential Comparing expressions (4.23) and (4.29), we see that they are connected by the relation   This leads to the expression µ = µ 0 n 0 + µ 1 n 1 , (4.33) in which the condensate fraction n 0 and the fraction of uncondensed atoms, n 1 , are introduced, Invoking equation (4.30), we get Sometimes, one requires that µ should be equal to µ 0 and µ 1 , which forces us to assume that the anomalous average σ 1 should be zero. Such a requirement has no physical reason. In addition, it can be shown by direct calculations [14,19,37] that the anomalous average is always comparable with or larger than either the density of uncondensed atoms or that of condensed atoms. Therefore, there is no such a region of parameters, where it could be admissible to neglect the anomalous average, but to keep the normal density and the density of condensed atoms. The sole possibility could be at temperatures close to zero and asymptotically weak interactions, when, though the anomalous average is three times larger than the normal density, both of them are much smaller than the condensate density. Then, it could be possible to omit both the anomalous average and the normal density, keeping only the condensate density. But neglecting one of them, though keeping another one, is mathematically wrong. Moreover, neglecting the anomalous average is not merely mathematically incorrect, but it is qualitatively deficient, making thermodynamics non-self-consistent, disturbing the condensate transition to the first order, and resulting in unphysical divergences of compressibility and structure factor [38].
Since this section is based on the Bloch representation, it is necessary to briefly describe how the Bloch functions could be defined. Formally, as has been mentioned above, the basis of Bloch functions should be chosen as a set of natural orbitals [35], since this gives a diagonal expansion for the density matrix (4.3). However, the problem is that the density matrix (4.2) is not known explicitly. Hence, it is impossible to find its exact eigenfunctions representing the natural orbitals. A standard way is to define the Bloch functions as solutions to the equation − ∇ 2 2m + V L (r) ϕ nk (r) = E nk ϕ nk (r) . It is also possible to define Bloch functions as eigenfunctions of the nonlinear Schrödinger equation [19], including the interaction term into equation (4.35). Then, calculations become essentially more complicated. In addition, there arises a problem of nonorthogonality of eigenfunctions of the nonlinear equation. Thus, the simplest way is to use the solutions to the linear equation (4.35) as a basis, complimenting it by conservation conditions (4.4).

Wannier representation
The field operator of atoms can be expanded over the basis of Wannier functions, ψ(r) = n jĉ n j w n (r − a j ) ,

23002-7
here, the operatorsĉ j satisfy the Bose commutation relations. The parameters entering the Hubbard Hamiltonian (5.2) can be calculated in the tight-binding approximation. A detailed demonstration of this calculation can be found in reference [19]. For a threedimensional space in this approximation, we find the expressions The explanation of the notations for V 0 , E R , and k 0 are given in section 2.
The single-band Hamiltonian (5.2) is called the boson Hubbard model. It is possible to generalize this model by taking into account two or more bands [39,40]. Here, we consider the single-band case, when the system displays Bose-Einstein condensation, though.
Employing the Bogolubov shift  In terms of the operators c j , the Bogolubov shift reads as follows: where the number of the nearest neighbors is denoted as (5.10) The first-order term is canceled by the linear cancelerΛ. The second-order term is The third-order term reads as follows: The fourth-order terms is The fraction of uncondensed atoms takes the form (5.14) where the lattice ideality is used. For the dimensionless anomalous average, we have The necessary condition of the system stability The operators c j can be expanded over the Fourier basis, (5.18) where k runs over the Brillouin zone.
Let us consider a cubic lattice. Then, the second-order term (5.11) becomes The third-order and fourth-order terms are (5.20) and, respectively, In the HFB approximation, the third-order term is zero, due to condition (5.8). And the fourth-order term in the HFB approximation reads as follows: ∆ ≡ νU (n 0 + σ) (5.24) 23002-9 for the grand Hamiltonian (4.5), we obtain The condensate chemical potential (5.17) in the HFB approximation becomes µ 0 = −z 0 J + νU (1 + n 1 + σ) .  (5.27) in which and the Bogolubov spectrum is The condition of the condensate existence (4.22) yields As is seen, µ 0 does not coincide with µ 1 , by analogy with relation (4.30). The anomalous average cannot be neglected, as is explained in section 4. For the quasi-momentum atomic distribution and for the quasi-momentum representation of the anomalous average, respectively, we find 2T . (5.34) This shows that the normal and anomalous averages are connected by the relation σ 2 k = n k (1 + n k ) − 1 4sinh 2 (ε k /2T ) .
The condensate fraction reads as with the integration over the Brillouin zone. Let us emphasize again that the anomalous average cannot be neglected for the principal reason. As is evident form the above formulas, the anomalous average can be zero only when there is no condensate, n 0 = 0. Hence, there is no gauge symmetry breaking. However, as soon as there appears Bose-Einstein condensate, the gauge symmetry becomes broken, and the anomalous average is never zero. It is always comparable with or larger than either the density of uncondensed atoms or that of condensed atoms.

(6.2)
For the ground-state energy we have is the operator of the total number of atoms. Since the first term N 0 is a non-operator number, one has var(N ) = var(N 1 ).

23002-11
The number-of-atom operator variance defines the isothermic compressibility The atomic fluctuations are, of course, normal and the compressibility is finite everywhere below T c . The compressibility can diverge only at the critical point T c .
Bose-Einstein condensation is a second-order phase transition occurring at a temperature T c , where n 0 = 0 and σ = 0. At this point, the atomic density is Solving this equation in the Debye approximation, we obtain the critical temperature (6.11) This tells us that T c is not defined for d = 1 and T c = 0 for d = 2. In three dimensions, we have The general equation for the superfluid fraction [14,19] can be written in the form n s = 1 − Q Q 0 , (6.13) with the classical dissipated heat where d is spatial dimensionality, and Q = var(P ) 2mN (6.15) is the actual dissipated heat, expressed through the variance of the momentum operator P ≡ ψ † 1 (r)(−i ∇)ψ 1 (r) dr . (6.16) In an equilibrium system, this variance is var(P) = 〈P 2 〉 . (6.17) Note that the condensed fraction does not contribute to the operator of momentum (6.16) due to the lattice periodicity [19].
For a three-dimensional cubic lattice, with a lattice spacing a, we obtain is used, derived in the tight-binding approximation. Here, the notation means an effective localization length. For a three-dimensional cubic lattice, the relations are valid, which yield the ratio a 2 Then, equation (6.19) can be written as follows: Therefore, the dissipated heat (6.18) is written as follows: In this way, the self-consistent mean-field approximation allows us to calculate any thermodynamic characteristic.

Conclusion
A self-consistent approach, based on the use of a representative statistical ensemble, developed earlier for uniform Bose-condensed systems, is extended to Bose atoms in optical lattices. The approach ensures a gapless spectrum of collective excitations, the validity of conservation laws, and self-consistent thermodynamics. It is shown that the approach can be applied to the lattices with a weak binding as well as with tight binding. For the former case, the Bloch representation is more appropriate, while for the latter case, the Wannier representation is more suitable. Both the Bloch and the Wannier representations lead to a similar description. The results are compared for the self-consistent Hartree-Fock-Bogolubov approximation. A convenient general formula for the superfluid fraction of atoms in an optical lattice is derived.
The HFB approximation, used here, is based on the assumption of the condensate existence, which is taken into account by means of the Bogolubov shift, explicitly breaking the global gauge symmetry of the system. This approximation, therefore, is assumed to provide good description, when the Bose condensate is present, and may be inappropriate when the system passes to an insulating state. This implies that the HFB approximation for optical lattices can provide an accurate description for spatial dimensions larger than one (d > 1) and nonzero temperatures below the Bose-Einstein condensation temperature, 0 < T < T c .
The case of zero temperature requires a special consideration. Cubic optical lattices at zero temperature and unity filling factor ν = 1 have been extensively studied, mainly from the viewpoint of an insulating state, with the purpose of defining the stability boundary of this state, corresponding to the critical transition to the superfluid state. The dimensionless parameter u ≡ U z 0 J 23002-13 has been varied. For a cubic lattice, the number of nearest neighbors is z 0 = 2d. This zero-temperature problem has been treated in the Gutzwiller approximation [41,42], dynamical mean-field approximation [43], direct numerical diagonalization [44], density-matrix renormalization group [45], strong-coupling perturbation theory [46,47], and Monte Carlo simulations [48][49][50]. The critical values of the above parameter were found for d = 1 as u c = 1.8, for d = 2, as u c = 4.2, and for d = 3, as u c = 4.9. The HFB approximation underestimates quantum fluctuations at zero temperature. That is why it is applicable only for nonzero temperatures, when thermal fluctuations become more important. The advantage of using the developed approach for Bose-condensed atoms in optical lattices at finite temperatures is its relative simplicity, correct gapless spectrum, the validity of conservation laws, and self-consistent thermodynamics.