Finite-temperature degenerate perturbation theory for bosons in optical lattices

TL;DR

This paper develops a finite-temperature degenerate perturbation theory for bosons in optical lattices to accurately describe the Mott-insulator to superfluid phase transition, overcoming limitations of previous non-degenerate methods.

Contribution

It introduces a degenerate perturbation theory approach based on projection operators to resolve degeneracy issues near the phase transition in the Bose-Hubbard model.

Findings

Provides physically consistent order parameters near the phase transition.

Addresses degeneracy problems in mean-field approximations.

Improves theoretical understanding of phase behavior in optical lattice bosons.

Abstract

Bosonic atoms confined in optical lattices can exist in two different phases, Mott-insulator and superfluid, depending on the strength of the system parameters, such as the on-site interaction between particles and the hopping parameter. This work is motivated by the fact that non-degenerate perturbation theory applied to the mean-field approximation of the Bose-Hubbard Hamiltonian at zero and finite temperature fails to give consistent results in the vicinity of the Mott-insulator-superfluid phase transition, e.g., the order parameter calculated via non-degenerate perturbation theory reveals an unphysical behavior between neighboring Mott lobes, which is an explicit consequence of degeneracy problems that artificially arise from such a treatment. Therefore, in order to fix this problem, we propose a finite-temperature degenerate perturbation theory approach based on a projection operator formalism which ends up solving such degeneracy problems in order to obtain physically consistent results for the order parameter near the phase transition.

Equations

$${\hat{H}}_{BH}=\frac{U}{2}\sum _i{\hat{a}}_i^{†}{\hat{a}}_i^{†}{\hat{a}}_i{\hat{a}}_i-t\sum _{⟨i,j⟩}{\hat{a}}_i^{†}{\hat{a}}_j-\mu \sum _i{\hat{a}}_i^{†}{\hat{a}}_i.$$

$${\hat{H}}_{MF}=\frac{U}{2}\sum _i\left({\hat{n}}_i^2-{\hat{n}}_i\right)-\sum _i\mu {\hat{n}}_i-tz\sum _i\left({\Psi }^{\ast }{\hat{a}}_i+\Psi {\hat{a}}_i^{†}-{\Psi }^{\ast }\Psi \right),$$

$$\hat{H}=\frac{U}{2}\left({\hat{n}}^2-\hat{n}\right)-\mu \hat{n}-tz\left({\Psi }^{\ast }\hat{a}+\Psi {\hat{a}}^{†}-{\Psi }^{\ast }\Psi \right).$$

$$\mathcal{F}\left({\Psi }^{\ast },\Psi \right)=a_0+a_2|\Psi {|}^2+a_4|\Psi {|}^4+\cdots .$$

$${\hat{H}}_0=\frac{U}{2}\left({\hat{n}}^2-\hat{n}\right)-\mu \hat{n}+tz{\Psi }^{\ast }\Psi ,$$

$$\hat{V}=-tz\left({\Psi }^{\ast }\hat{a}+\Psi {\hat{a}}^{†}\right).$$

$$E_n=\frac{U}{2}\left(n^2-n\right)-\mu n+tz|\Psi {|}^2,$$

$$\mathcal{Z}=\mathrm{Tr}\left[{\mathrm{e}}^{-\beta \hat{H}}\right],$$

$$\hat{U}={\mathrm{e}}^{-\beta {\hat{H}}_0}{\hat{U}}_{\mathrm{I}}(\beta ),$$

$$\frac{d{\hat{U}}_{\mathrm{I}}(\tau )}{d\tau }=-{\hat{V}}_{\mathrm{I}}(\tau ){\hat{U}}_{\mathrm{I}}(\tau ),$$

$${\hat{V}}_{\mathrm{I}}(\tau )={\mathrm{e}}^{\tau {\hat{H}}_0}\hat{V}{\mathrm{e}}^{-\tau {\hat{H}}_0}$$

$${\hat{U}}_{\mathrm{I}}(\beta )$$

$$+{\int }_0^{\beta }d{\tau }_1{\int }_0^{{\tau }_1}d{\tau }_2{\int }_0^{{\tau }_2}d{\tau }_3{\int }_0^{{\tau }_3}d{\tau }_4{\hat{V}}_{\mathrm{I}}({\tau }_1){\hat{V}}_{\mathrm{I}}({\tau }_2){\hat{V}}_{\mathrm{I}}({\tau }_3){\hat{V}}_{\mathrm{I}}({\tau }_4).$$

$$Z=\sum _{n=0}^{\infty }{\mathrm{e}}^{-\beta E_n}⟨n|{\hat{U}}_{\mathrm{I}}(\beta )|n⟩,$$

$${\mathcal{Z}}^{(0)}=\sum _{n=0}^{\infty }{\mathrm{e}}^{-\beta E_n}.$$

$$\mathcal{F}=-\frac{1}{\beta }\ln \mathcal{Z}.$$

$$\mathcal{F}\approx -\frac{1}{\beta }\left[\ln {\mathcal{Z}}^{(0)}+\frac{{\mathcal{Z}}^{(2)}}{{\mathcal{Z}}^{(0)}}+\frac{{\mathcal{Z}}^{(4)}}{{\mathcal{Z}}^{(0)}}-\frac{1}{2}{\left(\frac{{\mathcal{Z}}^{(2)}}{{\mathcal{Z}}^{(0)}}\right)}^2\right].$$

$$a_0$$

$$a_2$$

$$a_4$$

$$a_0=E_n-tz|\Psi {|}^2,$$

$$a_2=tz+(tz{)}^2\left(\frac{n+1}{{\Delta }_{n,n+1}}+\frac{n}{{\Delta }_{n,n-1}}\right),$$

$$a_4={\left(tz\right)}^4\left[\frac{n\left(n-1\right)}{{\Delta }_{n,n-1}^2{\Delta }_{n,n-2}}+\frac{\left(n+1\right)\left(n+2\right)}{{\Delta }_{n+1,n}^2{\Delta }_{n,n+2}}+\frac{n{}^2}{{\Delta }_{n-1,n}^3}+\frac{{\left(n+1\right)}^2}{{\Delta }_{n+1,n}^3}+\frac{n\left(n+1\right)}{{\Delta }_{n+1,n}^2{\Delta }_{n-1,n}}+\frac{n\left(n+1\right)}{{\Delta }_{n,n-1}^2{\Delta }_{n+1,n}}\right],$$

$$\frac{\partial \mathcal{F}}{\partial |\Psi {|}^2}=a_2+2a_4|\Psi {|}^2=0,$$

$$|\Psi {|}^2=-\frac{a_2}{2a_4}.$$

$$\hat{P}=|n⟩⟨n|+|n+1⟩⟨n+1|,$$

$$\hat{Q}=\sum _{m\notin P}|m⟩⟨m|.$$

$$\hat{H}={\hat{H}}_0+\left(\hat{P}+\hat{Q}\right)\hat{V}\left(\hat{P}+\hat{Q}\right)$$

$${\hat{\mathcal{H}}}_0$$

$$\hat{\mathcal{V}}$$

$${\hat{\mathcal{H}}}_0^{(\mathrm{nd})}=\left(\begin{array}{cc}{E}_n& -tz\Psi \sqrt{n+1}\\ -tz{\Psi }^{\ast }\sqrt{n+1}& E_{n+1}\end{array}\right).$$

$${\mathcal{E}}_{\pm }$$

$$|{\Phi }_{\pm }⟩$$

$$\hat{U}={\mathrm{e}}^{-\beta {\hat{\mathcal{H}}}_0}{\hat{\mathcal{U}}}_{\mathrm{I}}.$$

Full text

I Introduction

Optical lattices are laser arrangements which enable a spatially periodic trapping of atoms due to the interaction between the external electric field and the induced dipole moment of the atoms jaksch ; bloch ; pethick ; pitaevskii . A gas composed of bosonic atoms in an optical lattice can be described by the Bose-Hubbard model ueda ; lewenstein , which has three main parameters: the on-site interaction parameter, the hopping parameter, and the chemical potential. Depending on the magnitude of the parameters, the system can realize two different phases, the Mott insulator or the superfluid phase fisher ; greiner1 ; greiner2 ; widera ; folling ; gunter ; ospelkaus ; lewenstein2 ; gerbier . If the on-site interaction parameter is much larger than the hopping parameter, the system is in the Mott insulator (MI) phase. This phase is characterized by a strong localization of the atoms. By decreasing the amplitude of the periodic potential, so that the hopping parameter becomes much larger than the atom-atom interaction parameter, the system undergoes a phase transition to a superfluid (SF) phase, where the atoms are delocalized. Such differences in the localization of bosons make it possible to measure the phase the system is currently in through time-of-flight experiments greiner1 ; hoff . The MI-SF transition can happen even at zero-temperature, driven by quantum mechanical fluctuations, thus characterizing a quantum phase transition sachdev .

The main difficulty in solving the Bose-Hubbard Hamiltonian is the nonlocality of the hopping term. Thus, a common path for obtaining a first approximation of the MI-SF quantum phase diagram is a mean-field calculation, which approximates the Bose-Hubbard Hamiltonian by a sum of local Hamiltonians fisher . Following this simplification, Rayleigh-Schrödinger perturbation theory (RSPT) is typically used for obtaining the mean-field phase diagram at zero temperature sachdev . However, there are problems that arise from RSPT, since it does not properly deal with the degeneracies that occur between two consecutive Mott lobes. One of such RSPT problems concerns the calculation of the condensate order parameter which falsely vanishes between consecutive Mott lobes melo ; martin .

Also, other methods have been suggested in order to improve the mean-field quantum phase diagram for bosons in optical lattices, such as in Ref. santos , which uses a variational method and the field-theoretic concept of the effective potential. Also, the MI-SF phase transition at arbitrary temperature was investigated in Ref. bradlyn using an effective action approach. Furthermore, in Refs. grass1 ; grass2 an effective action approach was derived for the Bose-Hubbard model within the Schwinger-Keldysh formalism in order to handle time-dependent problems at finite temperature. Likewise, melo implemented a nearly degenerate perturbation theory for the zero-temperature case, which led to better results for the order parameter (OP) when compared to those from the RSPT calculations. More recently, Brillouin-Wigner perturbation theory was applied in order to correct such degeneracy-generated unphysical results at zero-temperature martin . It turns out that nondegenerate finite-temperature perturbation theory, as applied in Refs. bradlyn ; ednilson , also presents degeneracy problems similar to RSPT. Indeed, RSPT is equivalent to the usual finite-temperature perturbation theory in the zero-temperature limit. Therefore, degeneracy-related problems are also expected to appear at low enough temperatures.

The present work is concerned with correcting the degeneracy problem that artificially arises from such perturbative approaches regarding a system composed of bosonic atoms confined in an optical lattice at finite temperature. Starting from the mean-field approximation for the Bose-Hubbard Hamiltonian and considering the Landau expansion for the order parameter in the vicinity of the MI-SF phase transition, we perform a perturbation theory in imaginary time. In addition, in order to fix the degeneracy problem, we introduce a projection operator formalism for the finite-temperature system. The main idea of this degenerate approach is to separate the Hilbert subspace in which the degeneracies are contained from its complement. This system is then exactly diagonalized inside the degenerate subspace, while the effects of the interaction between the two subspaces are taken into account perturbatively. Such a procedure corrects the degeneracy problem and leads to physically consistent results for the condensate density.

In Sec. II, we introduce the mean-field approximation for the Bose-Hubbard Hamiltonian in order to get rid of its nonlocality, which transforms the originial Bose-Hubbard Hamiltonian into a sum of local Hamiltonians, thus allowing us to work with separated lattice sites. In Sec. II.1, we briefly discuss the Landau theory for second-order phase transitions, which enables us to expand the system free energy as a power series of the OP in the vicinity of the MI-SF phase transition. Furthermore, the Landau theory gives us the equation for the phase boundary. In Sec. II.2, we apply the imaginary-time-dependent nondegenerate perturbation theory (NDPT) considering the system in the vicinity of the phase transition. In this way, we treat the OP perturbatively in order to get expressions for the Landau coefficients and, consequently, obtain the phase boundary as well as the OP, close to the phase boundary. Following the calculation of the Landau coefficients, Sec. II.3 exposes the unphysical behavior of the OP between two consecutive Mott lobes, which is an explicit consequence of the degeneracies that are not adequately handled within the framework of NDPT at low enough temperatures. This inconsistency in the theory is the motivation for our proposed degenerate approach, which is worked out in detail in Sec. III. In Sec. III.1, we evaluate the condensate densities for different temperatures and values of system parameters making use of our proposed finite-temperature degenerate perturbation theory (FTDPT). In Sec. III.2, we turn our attention to a region between two consecutive Mott lobes in the phase diagram, where the superfluid clearly dominates and also a region where the NDPT fails at very low temperatures. We compare the results of the NDPT and the FTDPT to conclude that our degenerate method corrects all inconsistencies. Finally, we calculate the particle densities in Sec. III.3 for different temperatures and describe the existence of a melting of the wedding cake like structure.

II Mean-Field Approximation

The description of spinless bosonic atoms confined in an optical lattice is given by the Bose-Hubbard Hamiltonian,

[TABLE]

The respective parameters are the following: $t$ represents the hopping parameter, $U$ stands for the on-site interaction parameter describing the interaction between particles, and $\mu$ denotes the chemical potential within a grand-canonical description. Furthermore, $\hat{a}^{\dagger}_{i}$ and $\hat{a}_{i}$ are the usual bosonic creation and annihilation operators at site $i$, respectively. Note that in this model only nearest neighbor hopping is allowed and this restriction is depicted by $\langle i,j\rangle$.

Due to the nonlocal character of the hopping parameter, a standard mean-field approximation is usually the simplest way to solve a problem of this kind. The fundamental concept behind such an approach is to approximate the nonlocal hopping term by a local one. This procedure results in the mean-field Hamiltonian melo ; ednilson ; santos ; fisher ,

[TABLE]

where $z$ denotes the number of nearest neighbors, $\Psi\equiv\langle\hat{a}_{i}\rangle$ and $\hat{n}_{i}\equiv\hat{a}^{\dagger}_{i}\hat{a}_{i}$ is the number operator. Since (2) is a sum of local Hamiltonians, we restrict ourselves in the following to the one lattice site Hamiltonian,

[TABLE]

II.1 Landau theory

Landau argued that the free energy can be written as a polynomial function of the order parameter in the vicinity of a phase transition landau ,

[TABLE]

Since $\mathcal{F}$ is considered to be an analytic function of $\Psi$ and since the Bose-Hubbard Hamiltonian described by (1) possesses a global $U(1)$ phase invariance, $a_{n}$ will not vanish only for even values of $n$. Therefore, for small values of $|\Psi|$, we can consider only the lowest-order terms in (4), i.e., $a_{0}$, $a_{2}$, and $a_{4}$ as nonvanishing. For $a_{4}>0$, a second-order phase transition may occur. This originates from the fact that, if $a_{2}>0$, the only solution of $\partial\mathcal{F}/\partial\Psi=0$ is $\Psi=0$, thus corresponding to the MI symmetric phase, while if $a_{2}<0$, the effective potential $\mathcal{F}$ has infinitely many minima with $|\Psi|\neq 0$ which differ only in the phase of $\Psi$ and corresponds to the SF phase ednilson . Thus we conclude that the condition $a_{2}=0$ defines the boundary between the two phases.

II.2 Nondegenerate perturbation theory

As mentioned before, the transition from Mott insulator to superfluid is followed by a symmetry breaking and can be characterized by a change of the order parameter from zero to a nonzero value. Since we are considering our system in the vicinity of a phase transition, $\left|\Psi\right|$ has a small value and hence we treat the hopping term in (3) as a perturbation. Thus, (3) decomposes according to $\hat{H}=\hat{H}_{0}+\hat{V}$ into an unperturbed Hamiltonian

[TABLE]

and the perturbation

[TABLE]

The unperturbed eigenenergies are

[TABLE]

where the quantum number $n=0,1,2,\cdots$ corresponds to the number of bosons per site.

At this point we are interested in evaluating how the perturbation changes the free energy of the system. For this purpose, we must work out the partition function,

[TABLE]

in order to obtain the free energy of the system. The quantum-mechanical evolution operator with imaginary time, i.e., $\hat{U}=\mathrm{e}^{-\beta\hat{H}}$, can be factorized according to

[TABLE]

where $\hat{U}_{\rm{I}}(\beta)$ is the interaction picture imaginary-time evolution operator. The equation for the time evolution operator in the interaction picture is sakurai

[TABLE]

with

[TABLE]

and $\hbar=1$.

Equation (10) has to be solved with the initial value $\hat{U}_{\rm{I}}(0)=1$. This can be done iteratively, thus allowing the construction of a perturbative expansion. Up to fourth order in the interaction we have ednilson

[TABLE]

It turns out that for the perturbative Hamiltonian in (6) all odd-order terms in (II.2) vanish. Therefore, we can restrict ourselves to the calculation of the zeroth, second, and fourth-order terms in (II.2).

Making use of the time-evolution operator in the interaction picture $\mathcal{Z}=\mathrm{Tr}\left[\mathrm{e}^{-\beta\hat{H}_{0}}\hat{U}_{\rm{I}}(\beta)\right]$, we calculate the partition function,

[TABLE]

with the single-site eigenstates $|n\rangle$ corresponding to the eigenvalues in Eq. (7). The zeroth-order term is

[TABLE]

The second and fourth-order terms $\mathcal{Z}^{(2)}$ and $\mathcal{Z}^{(4)}$ are calculated in detail in App. A.

From the partition function, we then evaluate the free energy,

[TABLE]

Up to fourth order we get

[TABLE]

Therefore, by comparing (4) and (16) we read off the Landau expansion coefficients

[TABLE]

At zero temperature, we obtain results which are equivalent to RSPT. In particular, the Landau expansion coefficients reduce to

[TABLE]

where $\Delta_{i,j}\equiv E_{i}-E_{j}$.

The explicit solution for $a_{2}=0$ gives the phase boundaries in Fig. 1, as in ednilson . From Fig. 1 we read off that thermal fluctuations mainly affect the phase boundary between two Mott lobes due to fluctuations in the number of bosons per site.

II.3 Nondegenerate perturbation theory inconsistency

As already pointed out, NDPT is expected to exhibit degeneracy-related problems. Indeed, by directly observing the coefficient denominators in (18b) and (18c) we clearly identify such degeneracy problems. Whenever $\mu/U$ becomes an integer $n$, there is an equality between two consecutive energy values, for instance $E_{n}$ and $E_{n+1}$, thus characterizing a divergence in these expressions.

According to (4), we can consider the Landau expansion up to fourth order for the free energy in the vicinity of a phase transition. Extremizing (4) with respect to the order parameter leads to

[TABLE]

with the solution in the superfluid phase

[TABLE]

The plot of $|\Psi|^{2}$ as a function of $\mu/U$ making use of (17b) and (17c) is interesting for our purposes since it reveals an unphysical behavior, which is a consequence of NDPT: the order parameter approaches zero at a point where no phase transition occurs. Fig. 2 shows equation (20) for $tz/U=0.2$ for a varying chemical potential. We observe that, indeed, the OP is well-behaved in most parts of the diagram. However, it also shows an inconsistency: at integer values of $\mu/U$ the order parameter for the zero-temperature system goes to zero, while for $T>0$ it mimics the zero-temperature behavior by decreasing its values despite of not vanishing.

Since, for finite temperatures, NDPT also shows a similar unphysical behavior typical of RSPT, in the following section we demonstrate how such problems can be fixed at finite temperature.

III Degenerate Approach

In this section, we introduce our method which consists of a degenerate perturbative calculation making use of projection operators. As we aim at describing the region between two neighboring Mott lobes, we define a subspace of the Hilbert space which is composed of two degenerate states with $n$ and $n+1$ particles via the projection operator

[TABLE]

this way the corresponding complementary operator is given by

[TABLE]

We begin our analysis by considering the one-site mean-field Hamiltonian (3) and regard, as in Sec. II.2, the hopping term (6) as a perturbation for (5). We multiply both sides of the perturbation by the identity operator, $\hat{\mathbb{I}}=\hat{P}+\hat{Q}$, and get

[TABLE]

Then we define the new unperturbed Hamiltonian and the new perturbation as

[TABLE]

The Hamiltonian in equation (24a), written in the basis of the unperturbed eigenstates, is a block diagonal matrix, whose only nondiagonal block is

[TABLE]

Its eigenvalues and eigenstates are given by

[TABLE]

As pointed out in Sec. II.2, we must evaluate the partition function (8) in order to calculate the free energy (15). The only difference is that now we are working with the new unperturbed Hamiltonian (24a) and the new perturbation (24b). With this the time evolution operator now reads

[TABLE]

The equation for the imaginary-time-evolution operator in the interaction picture is

[TABLE]

with

[TABLE]

The solution for equation (28) with the initial condition $\hat{\mathcal{U}}_{\rm{I}}(0)=1$ up to second order is given by

[TABLE]

Evaluating the partition function $\mathcal{Z}=\mathrm{Tr}\left[\mathrm{e}^{-\beta\hat{\mathcal{H}}_{0}}\hat{\mathcal{U}}_{\rm{I}}(\beta)\right]$, we have

[TABLE]

The zeroth-order term in (30) yields in (31)

[TABLE]

Furthermore, we read off from (6), (29) and (24b) that the first-order contribution in (31) must vanish.

Finally, the second-order term, which is calculated in detail in Appendix B, gives

[TABLE]

where we have introduced the abbreviation $\Delta_{i,\pm}\equiv E_{i}-\mathcal{E}_{\pm}$.

From Eq. (III), we observe that the difference between the degenerate energies $E_{n}$ and $E_{n+1}$ will no longer appear in the denominator of the free energy as it did in the NDPT treatment, thus solving the degeneracy-related problems discussed above.

III.1 Condensate Density

Now we turn our attention to the calculation of the condensate density, which turns out to coincide with the superfluid density in the mean-field approximation martin . Our degenerate approach, up to second order, results in the partition function given by $\mathcal{Z}=\mathcal{Z}^{(0)}+\mathcal{Z}^{(2)}$ with (32) and (III), which is free from any divergence despite of the degeneracies. From the partition function, we obtain for the system free energy (15)

[TABLE]

Hence, we evaluate the condensate density $|\Psi|^{2}$ by evaluating

[TABLE]

We apply this procedure by considering different temperatures between the Mott lobes $n=0$ and $1$ in Fig. 3, and between $n=1$ and $2$ in Fig. 4.

In order to check the fidelity of the calculated condensate densities we must observe the phase boundary evaluated by FTDPT, which is given by

[TABLE]

This procedure leads to the same phase diagram evaluated by NDPT. From Fig. 1 we read off that for small values of $tz/U$ there are bigger portions of values of $\mu/U$ where the condensate density can be evaluated, since we regard the Landau expansion of the order parameter being valid in the vicinity of the phase transition, i.e. the smaller the hopping, the bigger the region of the calculated condensate density. Therefore, we conclude that we are able to reliably calculate $\left|\Psi\right|^{2}$ via FTDPT near the phase boundary in Fig. 3. Also, we observe that for $\mu/U=0$ the condensate densities no longer vanish or approach zero as they do when calculated from NDPT. Regarding the condensate densities calculated between the first and the second Mott lobes, Fig. 4, we also find that the decreasing behavior characteristic of the NDPT in the degeneracy point $\mu/U=1$ is absent.

III.2 Comparison between NDPT and FTDPT

Now we turn our attention to the point between two consecutive Mott lobes in order to analyze the differences between the condensate densities calculated via NDPT and FTDPT between the Mott lobes $n=0$ and 1, and $n=1$ and 2, as shown in Fig. 5. We observe that the NDPT gives condensate densities that approach zero or have a decreasing behavior at the degeneracy point, which corresponds to $\mu/U=0$ for the region between $n=0$ and $n=1$ and is depicted in Figs. 5(a) and 5(b); while for the region between the first and the second Mott lobes, i.e., $n=1$ and 2, the degeneracy occurs at $\mu/U=1$ and is depicted in Fig. 5. Such behavior indicates an inaccuracy of the theory, since it mimics the unphysical vanishing of the OP typical of RSPT, which is a direct consequence of not taking into account the degeneracies that happen in between two consecutive Mott lobes. While NDPT presents such unphysical behavior due to the incorrect treatment of degeneracies, FTDPT gives consistent results for the condensate density between two consecutive Mott lobes.

We observe from Fig. 5 that the condensate densities calculated via FTDPT, which are represented by the solid lines, do not present any decreasing behavior in the vicinity of the degeneracy, concluding that they are consistent in all considered regions of the phase diagram. In particular, at $\mu/U=0$ and 1 the condensate densities no longer vanish or present a decreasing behavior as they do when calculated from NDPT. The decreasing behavior presented by the condensate densities calculated via NDPT can clearly be observed by the solid lines in Fig. 5. Such decreasing behavior is a direct consequence of the incorrect treatment of degeneracies by NDPT, which happens to occur between two consecutive Mott lobes.

III.3 Particle density

We calculate the particle density,

[TABLE]

making use of our developed FTDPT. We consider different temperatures and different hopping values for the purpose of analyzing their effects on the density of particles. We plot the resulting equation of state for two different values of the hopping parameter and four different values of the temperature, thus observing the melting of the structure as in Refs. bloch2 ; gerbier2 , as shown in Fig. 6.

We observe the effects that the change of both the temperature and the hopping have upon the particle density in Fig. 6. First, we conclude that increasing the temperature makes the particle density to vary more smoothly when compared to those particle densities with lower temperatures. This fact is due to thermal fluctuations, which make the system more feasible to exist in the superfluid phase. Also, by comparing the left panel to the right one we observe the melting of the Mott lobes due to an increased hopping, which is also very intuitive: the particles, having more kinetic energy, are more likely to hop from one site to another, which is characteristic for the SF phase. Another factor responsible for making the curves smoother is the increase of the chemical potential, $\mu/U$. The reason for this relies on the fact that the bigger $\mu/U$ becomes the smaller the Mott lobes are, as can be seen in Fig. 1. Thus, the system is more likely to exist in the superfluid phase for bigger values of $\mu/U$.

Now we must turn our attention to the points of the figures where the degeneracies happen, which correspond to $\mu/U=1$ and $2$. We observe that our calculations lead to no decreasing behavior happening at those regions, meaning that our developed FTDPT method possess no inconsistency in the calculation of the equation of state for the mean-field approximation of bosonic atoms confined in optical lattices. As the NDPT leads to a vanishing behavior of the condensate densities, or at least a decreasing one, in the vicinity of the degeneracies, i.e., $\mu/U=1$ and $2$ in Fig. 6, we should also expect the same behavior in the evaluation of the particle density. Finally, we also conclude that FTDPT gives reliable results for the particle density since there is no decreasing behavior in the vicinity of the degeneracies in Fig. 6.

IV Conclusions

By using a projection operator formalism we were able to generalize the usual nondegenerate perturbation theory for the mean-field approximation of the Bose-Hubbard Hamiltonian at finite temperatures. With this we have solved the degeneracy problems which are typical for nondegenerate perturbation theories. We introduced the mean-field considerations applied to the Bose-Hubbard Hamiltonian, followed by the Landau theory regarding second-order phase transitions. Also, we showed that NDPT results in an inconsistent behavior for the order parameter: it predicts a phase boundary in a region of the phase diagram where there should be none. Subsequently, we developed a degenerate perturbative method based on a projection operator formalism that corrects all such contradictions which arise from NDPT due to degeneracies that occur between two adjacent Mott lobes. Our approach was able to eliminate all the problematic divergences in the thermodynamic potential, which allowed us to accurately evaluate the condensate densities and the particle densities in the vicinity of the MI-SF phase transition for different temperatures and different hopping values. Further, we drew the finite-temperature phase diagrams in order to check the consistency of the previously calculated condensate densities.

It must be noted that for the zero-temperature regime, which is depicted in Fig. 4(d), the results for the condensate densities are similar to those obtained in martin , which uses a Brillouin-Wigner treatment for the perturbation expansion followed by a proper diagonalization in order to calculate the system free energy. That approach differs from the one used in this work. While the Brillouin-Wigner approach in martin is also able to correct the degeneracy problems from NDPT, it can only be applied to the zero-temperature case. On the other hand, the theory presented in this paper corrects degeneracy problems for both zero and finite temperatures, thus providing a relatively simple method for calculating the condensate density in a wide range of optical-lattice systems.

V Acknowledgments

We acknowledge financial support from the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) within the Collaborative Research Center SFB/TR 49 ”Condensed Matter Systems with Variable Many-Body Interactions” and SFB/TR 185 ”Open System Control of Atomic and Photonic Matter” (OSCAR Project No. 277625399), and from the binational project between CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, Improvement Coordination of Higher Level Personnel) and DAAD (Deutscher Akademischer Austauschdienst, German Academic Exchange Service) Probral (Programa Brasil-Alemanha, Brazil-Germany Program) No. 488/2018 Grant No. 88881.143936/2017-01. F.T.S. acknowledges CAPES for financial support. We also thank M. Bonkhoff, S. Eggert, M. Kübler, and C. Sá de Melo for helpful discussions. Support from Optics and Photonics Research Center (CePOF, Centro de Pesquisas em Óptica e Fotônica) through Grant No. 2013/07276-1 is acknowledged. F.E.A.d.S. acknowledges CNPq (Conselho Nacional de Denvolvimento Científico e Tecnológico, National Council for Scientific and Technological Development) for support through Grant No. 305586/2017-3.

Appendix A Nondegenerate calculation of $\mathcal{Z}^{(2)}$ and

$\mathcal{Z}^{(4)}$

In this appendix we are concerned with the detailed calculation of $\mathcal{Z}^{(2)}$ and $\mathcal{Z}^{(4)}$ via the nondegenerate approach. The second-order term reads

[TABLE]

Inserting (11) in (38) we have

[TABLE]

As $|n\rangle$ are eigenstates of $\hat{H}_{0}$ (39) reduces to

[TABLE]

According to (6) we have

[TABLE]

yielding

[TABLE]

The scalar products reduce (A) to

[TABLE]

Finally, the integrations yield

[TABLE]

where we have used the abbreviation $\Delta_{i,j}\equiv E_{i}-E_{j}$ for differences between two eigenvalues (7).

For the fourth-order term we have

[TABLE]

Inserting (5) and (11) in (45) gives

[TABLE]

According to (6) we have

[TABLE]

Using again (5) and (11) in (A) results in

[TABLE]

which, from the scalar products, reduces to

[TABLE]

The integrations result in

[TABLE]

Appendix B Degenerate calculation of $\mathcal{Z}^{(2)}$

This appendix is devoted to the evaluation of (31) for the second order of (30)

[TABLE]

We shall perform the calculation of each term separately and identify them as $\mathcal{Z}^{(2)}=\mathcal{Z}^{(2)}_{+}+\mathcal{Z}^{(2)}_{-}+\mathcal{Z}^{(2)}_{m}$.

As the evaluation of $\mathcal{Z}^{(2)}_{+}$ and $\mathcal{Z}^{(2)}_{-}$ are completely equivalent we perform a generic calculation for both contributions. Inserting the expression for the perturbation in the interaction picture (29) in the first term we have

[TABLE]

As $|\Phi_{\pm}\rangle$ are eigenstates of $\hat{\mathcal{H}}_{0}$ we get

[TABLE]

As we have $\hat{Q}|\Phi_{\pm}\rangle=0$, $\hat{P}|\Phi_{\pm}\rangle=|\Phi_{\pm}\rangle$ as well as $\hat{Q}$ and $\hat{P}$ represent hermitian operators (53) reduces to

[TABLE]

From (6) and (26b) and using the scalar products

[TABLE]

we have

[TABLE]

The evaluation of (B) leads to

[TABLE]

Evaluating the integrations in (57) yields finally

[TABLE]

The last term to be calculated is $\mathcal{Z}_{m}^{(2)}$. The first steps of this calculation are similar to those from the evaluation of $\mathcal{Z}_{\pm}^{(2)}$. Therefore, we have

[TABLE]

Applying the exponential operators to the eigenstates we are left with

[TABLE]

When we evaluate the multiplication among the terms between brackets, we must be aware of the fact that the cross terms, i.e., those that contain $\Psi^{2}$ or $\Psi^{*2}$ give zero since they contain the products $\langle m-1|\Phi_{\pm}\rangle$ and $\langle m+1|\Phi_{\pm}\rangle$, which cannot be both nonzero because it is not possible for $m+1$ and $m-1$ be equal to $n$ or $n+1$ at the same time. So, we are left with

[TABLE]

The integrations lead finally to

[TABLE]

Combining the contributions (58) and (B) the second-order term of the partition function reads

[TABLE]

Taking into account that the scalar products $\langle m-1|m^{\prime}\rangle$ and $\langle m+1|m^{\prime\prime}\rangle$ lead to one further restriction each in the summations, thus we finally obtain equation (III).