Finite temperature mean-field theory with intrinsic non-hermitian structures for Bose gases in optical lattices

We reveal a divergent issue associated with the mean-field theory for Bose gases in optical lattices constructed by the widely used straightforward mean-field decoupling of the hopping term, where the corresponding mean-field Hamiltonian generally assumes no lower energy bound once the spatial dependence of the mean-field superfluid order parameter is taken into account. Via a systematic functional integral approach, we solve this issue by establishing a general finite temperature mean-field theory that can treat any possible spatial dependence of the order parameter without causing the divergent issue. Interestingly, we find the theory generally assumes an intrinsic non-hermitian structure that originates from the indefiniteness of the hopping matrix of the system. Within this theory, we develop an efficient approach for investigating the physics of the system at finite temperature, where properties of the system can be calculated via straightforward investigation on the saddle points of an effective potential function for the order parameter. We illustrate our approach by investigating the finite temperature superfluid transition of Bose gases in optical lattices. Since the underlying finite temperature mean-field theory is quite general, this approach can be straightforwardly applied to investigate the finite temperature properties of related systems with phases possessing complex spatial structures.


I. INTRODUCTION
Since the first experimental realization of the Bose-Hubbard model [1] with ultracold atoms in optical lattices [2], ultracold gases in optical lattices have become one of the most important experimental platforms [3] to investigate rich physics that are relevant for a wide range of different fields, ranging from solid-state physics, over condensed matter physics, to high energy physics. This is attributed to the tremendous experimental development in ultracold gases in optical lattices in the last two decades, including realizing, for instance, optical lattices with different dimensionality and geometries [4][5][6][7], different types of interactions [8][9][10], artificial gauge fields [11][12][13], etc, which at the same time gives rise to considerable types of systems theat support rich unconventional quantum phases with nontrivial spatial structures.
A case in point is ultracold Bose gases in optical lattices, where for this type of systems with, for instance, long-range interactions [8][9][10], artificial gauge fields [11][12][13][14], etc, experimental and theoretical investigations [10][11][12][15][16][17][18][19][20][21][22] have shown that they can support exotic phases with nontrivial spatial structures, such as chargedensity waves, spin-density waves, unconventional superfluid, etc. In investigations of these systems, straightforward mean-field decoupling approach [23,24] and bosonic Gutzwiller mean-field variational approach [25,26] are two of the most widely employed theoretical tools that have efficiently revealed considerable nontrivial physics of these systems at zero temperature. However, concerning finite temperature properties of these systems, which are naturally indispensable for relevant experiments and the thorough understanding of their physical behavior, these two approaches seem not as efficient to be employed as in the zero temperature investigations [21], not even to mention that the legitimacy of the straightforward meanfield decoupling is not self-evident in different application scenarios and often overlooked in the literature [1,23,24]. This thus poses the natural demand for a reliable meanfield theory that could efficiently investigate the finite temperature properties of these systems.
We address this demand for the type of mean-field theory where the mean-field treatment is performed on the hopping term of the system to investigate the superfluid transition. Using single-component Bose gases in optical lattices as the concrete type of systems, we show that the widely used straightforward mean-field decoupling of the hopping term [23,24] generally results in problematic mean-field Hamiltonians with no lower energy bounds once the spatial dependence of the meanfield superfluid order parameter is taken into account (cf. Eq. (4) and the discussion below). We solve this problem by establishing the proper finite temperature mean-field theory [cf. Eq. (18)] via a systematic functional integral approach. In particular, we find the proper mean-field Hamiltonian generally assumes an intrinsic non-hermitian structure [cf. Eqs. (19,20)] that originates from the indefiniteness of the hopping matrix of the sys-tem. Based on this non-hermitian mean-field Hamiltonian, we develop an efficient and versatile approach for investigating the physics of the system at finite temperatures, where properties of the system can be calculated via straightforward investigations on the saddle points of an effective potential function for the order parameter [cf. Eqs. (21,22) and Fig. 1]. We illustrated our approach by using both a homogeneous and an inhomogeneous two-sublattice finite temperature mean-field theory to investigate the finite temperature superfluid transition of Bose gases in optical lattices, the results of which at low temperature agree with the well-established results from, for instance, bosonic Gutzwiller variational ansatz [25,26]. And we directly map out the finite temperature phase diagram of the system (cf. Fig. 1) and show how Mott lobes "melt" in the presence of thermal fluctuations (cf. Fig. 2 and Fig. 3). Since the underlying finite temperature mean-field theory is quite general, this approach can be straightforwardly applied to efficiently investigate the finite temperature properties of related systems with phases possessing nontrivial spatial structures.

II. MODEL HAMILTONIAN AND LIMITATIONS OF THE STRAIGHTFORWARD MEAN-FIELD DECOUPLING APPROACH
A. Model and its straightforward mean-field decoupled Hamiltonian To proceed with the discussion on a concrete basis, we consider the simplest single-band Bose-Hubbard model [1] on a two-dimensional (2D) square lattice which is widely used in describing the physics of ultracold gases in optical lattice [3]. Its explicit form readŝ is the bosonic creation (annihilation) operator for the Wannier state on the site i in the lowest band,n i ≡b † ib i is the particle number operator. Here, t is the positive hopping amplitude of bosons between nearest neighboring lattice sites denoted by i, j , U is the strength of the on-site interaction and µ is the chemical potential.
To investigate the properties of the system, straightforward mean-field decoupling for the hopping term is frequently employed [23,24]. It is done by plugging the decompositionb i = ψ i + b i − ψ i into the hopping term ofĤ and keeping only up to the first-order terms in (b i − ψ i ). This gives rise to the straightforward meanfield decoupled (SMFD) HamiltonianĤ SMFD ({ψ * i , ψ i }), the explicit form of which readŝ Here, ψ i is a complex variable and assumes the physical interpretation of the local superfluid (SF) order parameter. At zero temperature, the value of ψ i is expected to be determined by minimizing the ground state energy of B. Limitations of the straightforward mean-field decoupling approach When utilizing the straightforward mean-field decoupled HamiltonianĤ SMFD ({ψ * i , ψ i }), a homogeneous ansatz for ψ i , i.e., ψ i = ψ for ∀i is usually further assumed (see for instance Refs. [23,24]). Under this homogeneous ansatz,Ĥ SMFD assumes the form with z ≡ j= i = 4 being the coordination number of the 2D square lattice and N s being the total number of the lattice sites. Direct calculations based on H SMFD (ψ * , ψ) can give reasonable ground state properties of the system [23,24] and are consistent with related results from other mean-field approaches, for instance, bosonic Gutzwiller variational wave function approach [25,26]. In this regard, one would naturally expect this straightforward mean-field decoupling approach with a generic ansatz for the mean-field ψ i that takes into account its possible spatial dependence, should also be able to give reasonable predictions.
Somewhat unexpected, one actually finds this natural expectation is not true in general. To see this point, let us first take a simple inhomogeneous ansatz for ψ i as an example, where ψ i assumes the same value on each of the two sub-lattices of the 2D square lattice and can assume different values on different sub-lattice, i.e., ψ i = ψ σ if i ∈σ with σ = e, o. Here,σ denotes the set of all lattice sites of the sub-lattice with the index σ, and we denote two sub-lattice indices as e and o. The straightforward mean-field decoupled Hamiltonian based on this ansatz can be directly obtained. Its explicit form readŝ We notice that the first term ofĤ SMFD ({ψ * σ , ψ σ }) is a quadrature with respect to {ψ * σ , ψ σ } that is not positive definite, i.e., clearly not being a positive definite matrix. This directly indicates thatĤ SMFD ({ψ * σ , ψ σ }) generally assumes no lower energy bound with respect to {ψ * σ , ψ σ }, hence making the mean-field theory Eq. (4) not reliable anymore.
More generally, we notice from Eq. (2) that the straightforward mean-field decoupled Hamiltonian is not guaranteed to assume a lower energy bound with respect to {ψ * i , ψ i }, since its first term is a quadrature of {ψ * i , ψ i } that is not guaranteed to be positive definite. Although this may not cause serious issues in investigating zero temperature properties in related systems [27][28][29][30] since iterative self-consistent approaches are usually employed, the lack of the lower energy bound is fatal in the finite temperature case due to the presence of the thermal fluctuations of the order parameter field. This thus indicates the reliability of the mean-field theory Eq. (2) is limited to the applications combined with the ansatzes for ψ i that give positive definite quadratures of {ψ * i , ψ i }. In this regard, it is desirable that a reliable construction for the mean-field theory beyond this limitation can be established. Indeed, as we shall present in the following, such a finite temperature mean-field theory that can reliably work together with generic ansatzes for ψ i can be established via formulating the exact partition function of the system in the functional integral form and taking the classical limit of the quantum order parameter field introduced by the standard Hubbard-Stratonovich transformation [31][32][33].

III. FINITE TEMPERATURE MEAN-FIELD THEORY WITH INTRINSIC NON-HERMITIAN STRUCTURE
To solve the issue of the mean-field theory constructed by the straightforward mean-field decoupling approach, we first write down the exact partition function of the system Z = tr[e −βĤ ] in the standard coherent state functional integral formulation, the explicit form of which reads with the action S[{b * i (τ ), b i (τ )}] assuming the explicit form ) is the complex field that corresponds to the bosonic operatorb † i (b i ). Via the standard Hubbard-Stratonovich transformation (HST) [31][32][33], the fluctuating quantum superfluid order parameter field ψ i (τ ) can be directly introduced into the partition function to decouple the hopping term in the action as what has been done routinely in literature (see, for instance, Refs. [1,24]). However, we remark here that one should pay special attention to the convergence of the complex Gaussian integral involved in the HST, which is frequently overlooked in literature, since the "hopping matrix" (to be defined explicitly below) is generically an indefinite matrix.

A. Hubbard-Stratonovich transformation in the presence of the indefinite hopping matrix
To perform the HST in the presence of the indefinite hopping matrix, let us first reformulate the hopping term in the action, i.e., −t where T is the "hopping matrix" with its matrix elements denoted by T ij , B is an N s dimensional column vector that collects all b i , i.e., B ≡ (· · · , b i , · · · ) T . Since T is indefinite in general, the quadrature in Eq. (7) is indefinite too. We separate the positive definite part of the quadrature from its negative definite part by first diagonalizing T with a unitary matrix U , i.e., U and can be further written as a sum of the positive and the negative definite part, i.e., with I (±) being an N ± dimensional identity matrix and 0 (±) being an N ± × N ∓ zero matrix.
The HST of the whole hopping term is facilitated by performing two standard HST for the positive definite part and the negative definite part, separately. For the positive definite part, the HST reads while for the negative definite part, the corresponding HST reads In particular, we notice from Eq. (9) that the imaginary unit i appears as an overall prefactor of the coupling term betweenΨ (+) andB (+) in the HST for the positive definite part, while this is not the case for the negative definite one.
With the two separate HST in Eqs. (9, 10), we can straightforwardly express the hopping term in the partition function as an integral with respect to the N s di- where In particular, we notice that T H is a positive definite hermitian matrix, while T NH is a non-hermitian matrix in general, and gives rise to the non-hermitian structure of the mean-field Hamiltonian as we shall see in the following.

B. Non-hermitian structure of the mean-field Hamiltonian
With the HST in Eq. (11) for the hopping term, we can straightforwardly reformulate the partition function Z as a functional integral with respect to both B(τ ) and Ψ(τ ), the explicit form of which reads Up to now, Eqs. (14,15) essentially correspond to an exact reformulation of the quantum partition function of the system, based on which one can formulate different perturbative approaches following the similar way as the ones presented in, for instance, Refs. [1,24,34]. Here, to investigate the finite temperature properties of the system, we shall formulate a different approach based on analyzing the minimum of the effective potential function [cf. Eq. (22)] of the system constructed from its proper mean-field Hamiltonian [cf. Eq. (19)]. More specifically, to proceed further, we make the mean-field (classical) approximation that neglect the quantum fluctuations of the order parameter field Ψ(τ ), i.e., assume the superfluid order parameter field Ψ(τ ) does not depend on τ , i.e., Ψ(τ ) = Ψ. This gives rise to the mean-field partition function Z MF with the explicit form where the mean-field action S MF reads Noticing that the τ dependence of the second term of the mean-field action S MF only comes from {b * i (τ ), b i (τ )}, we can reformulate Z MF as (see Appendix A for more details) where the mean-field HamiltonianĤ MF (Ψ † , Ψ) assumes the explicit form withĤ (i) SS (Ψ † , Ψ) being a single site Hamiltonian that only involves operators on site i, the explicit form of which readsĤ We noticeĤ MF (Ψ † , Ψ) is generally a non-hermitian Hamiltonian due to the general non-hermiticity of T NH . Such a non-hermitian structure is reminiscent of nonhermitian Hamiltonians widely used to describe physics of open quantum systems (see Ref. [35] for a recent review), where the non-hermiticity originates from dissipations. In contrast, here, it originates from the hermitian hopping term with an indefinite hopping matrix. It is worth mentioning that although the non-hermitian structure of the mean-field Hamiltonian manifests itself explicitly in concrete calculations via giving rise to complex eigenvalues, this non-hermiticity does not give rise to any physical effect characteristic of open quantum manybody systems (cf. Sec. III D 2 for a detailed and concrete discussion). Moreover, we would like to remark that the above derivation of the mean-field Hamiltonian is completely general and can be directly applied to systems with complex hopping terms, for instance, systems with spin-orbit coupling [13,[18][19][20][21][22], effective magnetic flux [11,12,14], etc.
C. Finite temperature mean-field theory in terms of the potential function for the order parameter field To make the mean-field theory Eq. (18) a useful and efficient tool for investigating both zero temperature and finite temperature properties of the system, we further reformulate the integrand in the expression for Z MF in Eq. (18) as an exponential of a potential function of the superfluid order parameter field alone, i.e., where the potential function Ω t,U,µ,β (Ψ † , Ψ) reads .
The most appealing feature of this reformulation is that we can analyze the properties of the system by simply investigating the saddle points of the potential function Ω t,U,µ,β (Ψ † , Ψ). This is due to the fact that in the thermodynamic limit (N s → +∞), the partition function Z MF is exactly determined by the saddle point value of Ω t,U,µ,β (Ψ † , Ψ), hence the value of the SF order pa-rameterΨ at temperature T is determined by the value of Ψ that minimize Ω t,U,µ,β (Ψ † , Ψ). Although the explicit form of Ω t,U,µ,β (Ψ † , Ψ) can not be obtained analytically due to the trace in its expression, its value at given (Ψ † , Ψ) can be calculated very efficiently at sufficiently high accuracy by employing a large enough cut-off n max in the dimension of the local Hilbert space for the site i. For the temperature regime of most interests, where thermal fluctuations compete strongly with the two other energy scales in the system (the ones associated with the on-site interaction and the hopping), i.e., k B T ∼ U ∼ zt, a cut-off n max of O(10) is already large enough.
D. Finite temperature mean-field theory combined with the homogeneous and the inhomogeneous two-sublattice ansatz for ψ i To illustrate concretely how the finite temperature mean-field theory is applied, we used it together with the homogeneous and an inhomogeneous two-sublattice ansatz to investigate the superfluid transition of the system at finite temperature. To this end, let us first diagonalize the hopping term that appears in the action in Eq. (6) by a straightforward lattice Fourier transformation, i.e.,
In the upper panels of Fig. 1, we show two typical landscapes of the potential function Ω t,U,µ,β (ψ * , ψ) at different temperatures (other system parameters are kept the same with 2zt/U = 0.24 and the filling factor ρ ≡ N −1 s i n i = 1). We see that at the low temperature (k B T /U = 0.1 in this case) the potential function assumes a typical Mexican hat structure with its minimums located on a ring where ψ assumes non-zero modulus, while at the high temperature (k B T /U = 0.3 in this case) the potential function assumes a bowl structure with its unique minimum located at ψ = 0. This kind of change in the structure of the potential function directly indicates the temperature-driven phase transition from the superfluid to the normal phase is a second-order transition. Indeed, from the lower-left panel of Fig. 1, we see that the modulus of the superfluid order parameter |ψ| of the system continuously decreases to zero as T increases. The complete finite temperature phase diagram of the system at unit filling (ρ = 1) is shown in the lowerright panel of Fig. 1, where we notice that the critical temperature T C for the superfluid to normal phase transition only show obvious increasing behavior with respect to the hopping amplitude once the two competing energy scales, i.e., 2zt and k B T , are comparable with each other.
Moreover, away from unit filling, we calculated the superfluid order parameter |ψ| as a function of the chemical potential µ and hopping amplitude t at different temperatures as shown in Fig. 2. We see that at nearly zero temperature (k B T /U = 10 −2 ), the distribution of |ψ| on the µ-t plane still manifests a clear Mott-lobe structure, which agrees well with zero temperature results from other mean-field type theories [1,23,24,26] (cf. the black curve in each plot of Fig. 2, which corresponds to the zero temperature Mott insulator-superfluid transition boundary obtained by the mean-field theory employed in Refs. [1,23,24]). As T increases the superfluid region with relatively small hopping amplitude t vanishes first, making the Mott-lobe structure less and less apparent.
Finally, we remark that by comparing Eq. (25) with Eq. (3), we notice that under the homogenous ansatz for ψ i , the mean-field Hamiltonian constructed by the straightforward mean-field decoupling assumes exactly the same form as the one obtained via the systematic functional integral construction. On the one hand, this provides a solid theoretical ground for the straightforward mean-field decoupling approach combined with the homogenous ansatz. On the other hand, we should emphasize that such a formal consistency is purely due to the coincidence caused by the form of the homogeneous ansatz, under which only the negative definite part of the quadrature associated with the hopping term [cf. Eq. (8)] eventually appears in the mean-field Hamiltonian. Indeed, as we shall see explicitly in the following concrete example with an inhomogeneous two-sublattice ansatz, the functional integral construction gives rise to a meanfield Hamiltonian with a lower energy bound with respect to the superfluid order parameter field, which is in sharp contrast to the one constructed by the straightforward mean-field decoupling shown in Eq. (4).
FIG. 2. Superfluid order parameter |ψ| as a function of the chemical potential µ and hopping amplitude t at different temperatures. The calculations are performed within the homogeneous finite temperature mean-field theory. The black curve in each plot corresponds to the zero temperature Mott insulator-superfluid transition boundary obtained by the mean-field theory employed in Refs. [1,23,24]. At nearly zero temperature (kBT /U = 10 −2 ), the distribution of |ψ| on the µ-t plane still manifests a clear Mott-lobe structure. As the temperature T increases the superfluid region with relatively small hopping amplitude t vanishes first, making the Mott-lobe structure less and less apparent. See text for more details.

Two-sublattice finite temperature mean-field theory
Now, let us proceed to construct the mean-field theory by using a more general ansatz for ψ i , which assumes a two-sublattice structure, i.e., ψ i = ψ σ if i ∈σ with σ = e, o, andσ denoting the set of all lattice sites of the sub-lattice with the index σ. Without losing any generality, the explicit form of ψ i under this ansatz reads with K ≡ (π, π) T . Directly plugging this ansatz into Eq. (24), we can directly obtain (see Appendix B for derivation details) wherê The corresponding mean-field partition function Z MF reads where where η σ ≡ +1, −1 for σ = e, o, respectively. From the explicit form ofĤ SS is the hermitian conjugate ofĤ SS ] · tr[e −βĤ (o) SS ] is guaranteed to assume positive real values. Therefore, the potential function Ω t,U,µ,β ({ψ * σ , ψ σ }) is still a real-valued function despiteĤ MF ({ψ * σ , ψ σ }) being a non-hermitian Hamiltonian. Moreover, following the same line of derivation, one can also straightforwardly show that the potential function Ω t,U,µ,β is real-valued for a class of ansatzes for ψ i assuming the form ψ i = ϕ + δϕ cos(K · i) with K = (π/p, π/q) and p(q) being a generic positive integer, since H The concrete calculations based on the two-sublattice finite temperature mean-field theory Eq. (32) can proceed in a similar way as the one for the homogeneous meanfield theory. At different temperatures, the superfluid order parameter on the even and odd sublattice, i.e., |ψ e | and |ψ o |, as a function of the chemical potential µ and hopping amplitude t are shown in Fig. 3. We directly observe that |ψ e | and |ψ o | show exactly the same behavior, indicating the superfluid order parameter being homogeneous over the whole system as expected. Moreover, compared with results from the homogeneous finite temperature mean-field theory shown in Fig. 2, we can easily see that |ψ e | and |ψ o | manifest the same behavior as |ψ|, indicating the two-sublattice finite temperature mean-field theory is consistent with the homogeneous one. This is also consistent with the natural physical expectation that equilibrium phases of the system should be homogeneous. The calculations are performed within the two-sublattice finite temperature mean-field theory. The black curve in each plot corresponds to the zero temperature Mott insulatorsuperfluid transition boundary obtained by the mean-field theory employed in Refs. [1,23,24]. The distribution of |ψe| on the µ-t plane is exactly the same as the one of |ψo|, indicating the superfluid order parameter is always homogeneous over the whole system as expected. Moreover, the |ψe| and the |ψo| distribution on the µ-t plane are exactly the same as the |ψ| distribution obtained within the homogeneous mean-field theory shown in Fig. 2. See text for more details.

IV. CONCLUSIONS
The proper mean-field Hamiltonian can generally manifest distinct intrinsic structure from its original one, as the finite temperature mean-field theory for the Bose gases in optical lattices shows: due to the general indefiniteness of the hopping matrix, the proper mean-field treatment of the hopping term directly gives rise to a mean-field Hamiltonian that is non-hermitian. This nonhermitian structure poses no hindrance to the subsequent calculations. On the contrary, it facilitates the application of the mean-field theory in combination with generic space-dependent ansatzes for the order parameter field, based on which an efficient and versatile approach for calculating finite temperature properties of the system can be developed, as illustrated in the investigation of the finite temperature superfluid transition of Bose gases in optical lattices. As a direct application, our approach can be employed to study finite temperature properties of lattice Bose gases in optical cavities [10] whose order parameters can be inhomogeneous in certain sys-tem parameter region. Also, we believe our work will stimulate further efforts to investigate finite temperature properties of the unconventional superfluids in ultracold atom systems with spin-orbit coupling, effective magnetic flux, etc, and also other distinct intrinsic structures of the mean-field Hamiltonians that could appear in various mean-field treatments. Moreover, it is also intriguing to promote the finite-temperature mean-field theory developed here to systematically include corrections from quantum fluctuations similar to what has been done in the zero-temperature case [36].
We thank Lijun Lang for helpful discussions. This work was supported by NSFC ( Therefore, we can rewrite Z MF as the form shown in Eq. (18), i.e., In this appendix, we present the relevant derivation details for the explicit form ofĤ MF (Ψ † , Ψ) in the general case, and also its form under the homogeneous and the inhomogeneous two-sub-lattice ansatz for ψ i . From Eq. (23) we see that the hopping term is directly diagonalized by the lattice Fourier transformation. With the explicit form of the lattice Fourier transformation, one can directly write down the following relations concerning the corresponding unitary matrix U that diagonalizes the hopping matrix T , i.e., From Eq. (19) we see that in order to obtain the explicit form ofĤ MF (Ψ † , Ψ), one needs to get the explicit form for the quadratures, Ψ † T H Ψ, i (Ψ † T NH ) ibi , and ib † i (T NH Ψ) i . This can be done by plugging the explicit forms of T H and T NH in Eqs. (12,13), and the above relations for the unitary matrix U into the explicit form of the quadratures, i.e.
With the above explicit forms of the quadratures, we directly can get the explicit form forĤ MF (Ψ † , Ψ) shown in Eq. (24).
2. Explicit form ofĤMF(ψ * , ψ) in the homogeneous finite temperature mean-field theory To obtain the explicit form ofĤ MF (ψ * , ψ), we directly plug the homogeneous ansatz ψ i = ψ into the explicit forms of the three quadratures shown in Eqs. (B5, B6, B7), and can obtain where the explicit expression for E(k), i.e., E(k) = −4t (cos k x + cos k y ), is also used in the above derivations. With the above explicit forms of the quadratures, we directly get the explicit form ofĤ MF (ψ * , ψ) shown in Eq. (25).