# Phase diagrams of the Bose-Fermi-Hubbard model at finite temperature

###### Abstract

The phase transitions at finite temperatures in the systems described by the Bose-Fermi-Hubbard model are investigated in this work in the framework of the selfconsistent random phase approximation. The case of the hard-core bosons is considered and the pseudospin formalism is used. The density-density correlator is calculated in the random phase approximation and the possibilities of transitions from superfluid to supersolid phases are investigated. It is shown that the transitions between uniform and charge ordered phases can be of the second or the first order, depending on the system parameters.

###### pacs:

37.10.Jk, 67.85.Pq, 71.10.Fd^{†}

^{†}: J. Phys.: Condens. Matter

## 1 Introduction

Properties of the systems of ultracold atomic gases confined in optical lattices are intensively studied in the last years both theoretically and experimentally [1, 2, 3, 4, 5]. Special attention is paid to the mixture of bosons and spin polarised fermions (e.g., Li- Li, K- Rb, Li- Rb atoms). Such systems can be well described by the Bose-Fermi-Hubbard model (BFHM) [5], which is an extension of the Bose-Hubbard model. The BFHM can also be considered as a generalisation of the fermionic Hubbard model. It is known for the case of the Bose-Hubbard model that competition between two terms, one is connected with the on-site energy and another describes the nearest-neighbour hopping with the tunnelling parameter , defines the state of the system (when the kinetics energy dominates the ground state of the system is superfluid, in the opposite case the ground state is a Mott insulator) [6, 7]. For the case of the BFHM phase diagrams are more complicated because due to the presence of fermions the effective interaction between bosons is generated.

The Bose-Fermi mixtures in optical lattices have been studied using a variety of methods [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. In [8], it was demonstrated that a two-dimensional mixture of bosons and fermions develops a supersolid phase (this phase is characterized by the simultaneous presence of a density wave and phase order in condensate). The case of small fermion hopping was investigated in [9] in the framework of composite fermion approach and composite fermions were formed by a fermion and one or several bosons (bosonic holes) for attractive (repulsive) Bose-Fermi-interactions. In [10, 11], inhomogeneous (due the the presence of the trapping potential) mixtures of bosons and fermions were studied. Enhancement of the superfluidity due to the presence of fermions was predicted in [12]. The existence of the supersolid phase was confirmed in [15] using quantum Monte-Carlo simulations. A mixture of the mean field approximation for a bosonic part and the dynamical mean field theory for a fermionic part of the Hamiltonian was applied in [16] and the presence of a supersolid phase at weak Bose-Fermi interaction was established. The case of 1D Bose-Fermi-Hubbard model (BFHM) in the limit of large fermion hopping was investigated in [17] (the case of half filling was considered only and they did not observe the supersolid state).

It should be noted that the Bose-Fermi-Hubbard-type model can also be applied for the description of intercalation of ions in crystals (for example, lithium intercalation in TiO crystals). Theoretical investigation of such process in most cases were restricted to the numerical ab-initio and density-functional calculations [18, 19, 20]. It was shown that Li is almost fully ionized once intercalated and reconstruction of electron spectrum at intercalation takes place. Thus, ion-electron interaction can play a significant role in such systems. Another interesting feature of such crystals is a displacement of the chemical potential at intercalation into the conduction band. As a result, these crystals have metallic conductivity. At intercalation of lithium in TiO, phase separation into Li-poor and Li- rich phases occurrs and this two-phase behaviour leads to a constant value of the electrochemical potential [21, 22] (this fact is used when constructing batteries).

In our previous works [23, 24] we have formulated the pseudospin-electron model of intercalation. We have revealed that the effective interaction between bose-atoms (ions) can change its character depending on fermionic band filling, that leads to the charge-ordered phase or phase separation into the uniform phases with different concentrations of bosons and fermions. The ion-electron interaction was also considered in [25] at the investigation of thermodynamics of the spin- model of intercalation (the model was similar to the known Blume-Emery-Griffiths model), but the electron as well as ion transfer was not taken into account. Models of pseudospin-electron model type are widely used in physics of the strongly-correlated electron systems. Application of such models to high-temperature superconductors allows one to describe thermodynamics of anharmonic oxygen ion subsystem and explain the appearance of inhomogeneous states and the bistability phenomena (see [26]). Models of a lattice gas are also used at the description of ionic conductors and at the calculation of their conductivity starting from works of Mahan [27] and others [28, 29].

In this work we consider the hard-core limit (infinite on-site boson-boson interaction) of the BFHM at finite temperature (most previous investigations considered the zero-temperature case). Our paper is organized as follows. In section 2 we present the description of the model and give a self-consistent scheme for calculation of the density-density correlator (susceptibility) in the random phase approximation (RPA). In section 3 we present phase diagrams for different values of the model parameters. Special attention is paid to the influence of the temperature change on the phase transitions. We present our conclusions in section 4.

## 2 Model and Method

We consider the BFHM in the hard core limit. Using the pseudospin formalism, the Hamiltonian of the model is written in the following form

(1) |

The pseudospin variable takes two values ( when boson is present in a site and in the opposite case), while and are fermionic creation and annihilation operators, respectively. The first and the second terms in equation (2) are responsible for the nearest neighbour boson and fermion hopping, respectively; -term accounts for the boson-fermion interaction energy. To control the number of bosons and fermions we introduce the bosonic and fermionic chemical potentials and , respectively.

The unperturbed Hamiltonian in the mean feald approximation (MFA) is obtained using the following simplification:

(2) | |||

(3) |

The Hamiltonain becomes

(4) | |||

(5) | |||

(6) |

It is worth noting that application of the MFA to the strongly correlated systems in the limit of a weak on-site correlation (when there is no correlational splitting of the fermionic band) allows one to satisfactorily describe their properties.

To diagonalize the Hamiltonian we pass to -representation and perform the unitary transformation in the pseudospin subspace:

(7) | |||

(8) | |||

(9) | |||

(10) | |||

(11) |

where is the number of lattice sites.

To calculate the density-density correlator , we perform an expansion in powers of

(12) | |||

(13) | |||

(14) | |||

(15) |

the averaging is performed over the distribution with , where is the imaginary time ordering operator and is the inverse temperature.

To calculate the average values of the -products of the pseudospin and fermion operators, we utilize the diagrammatic technique based on Wick’s theorem for the spin operators [30] (besides the usual procedure for the Fermi operators). After elimination in this way of the nondiagonal operators we perform the semi-invariant expansion in order to calculate the mean values of the remaining products of the operators. At the summation of diagrams we restrict ourselves to the diagrams having a structure of multi-loop chains in the spirit of the random phase approximation (see [31]). The junctions between bosonic (pseudospin) Green’s functions and semi-invariants are realised by bosonic hopping and the fermionic loop . It is useful to introduce unperturbed bosonic and fermionic Green’s functions and , respectively, and semi-invariant .

Let us consider the Green’s function (with ). Typical RPA diargams for this Green’s function in the frequency representation are shown in figure 1.

We used the notations for the unperturbed bosonic Green’s function

(16) |

the fermionic loop

(17) |

semi-invariant and average value of the pseudospin variable .

The equation for this Green’s function is

(18) |

where (with ) and is the bosonic Matsubara frequency. These matrix equations (18) form three independent sets of equations of the third order and can be separately solved. For the case of the Green’s functions , , the matrices , and the unperturbed Green’s functions are

(19) | |||

(20) | |||

(21) | |||

(22) | |||

(23) | |||

(24) | |||

(25) | |||

(26) | |||

(27) | |||

(28) | |||

(29) |

Similar matrix equations can be written for the Green’s functions , , and , , with the corresponding matrices and (we do not present here these matrices). As a result, we can solve these three sets of equations of the third order and after some tedious algebra we derive the expression for the density-density correlator

(30) |

(31) | |||

(32) |

If we use the equation of motion method developed for the two-time Green’s function and decoupling in the spirit of Tyablikov approximation we can obtain the expression for the correlator which is similar to the equation (2) but differs from it due to the absence of the terms proportional to . These terms have appeared in the diagrammatic technique due to the presence of the semi-invariants and are important when we investigate the static limit . The equation of motion method does not allow us to take into account these terms and because of this we should use the diagrammatic technique. It should be noted that such an peculiarity was also pointed out in [32] at the investigation of the Bose-Hubbard model in the hard-core case.

## 3 Phase Diagrams

Lines of the instability with respect to the transition into the phase with charge ordering can be obtained using the condition of divergence of the static density-density correlator . We consider two cases: i) the transition from a normal (NR) nonsuperfluid uniform to nonsuperfluid charge-density-wave (CDW) phase ii) the transition from a superfluid phase to superfluid phase with long-range ordering (a supersolid phase). The equations for averages , , are obtained in the mean field approximation. Let us introduce two sublattices: , is a sublattice index, is an elementary cell index. Using the Hamiltonian , we can obtain the following equations for averages [24]

(33) | |||

(34) | |||

(35) |

with

(36) | |||

(37) | |||

(38) |

The grand canonical potential can be written as [24]

(39) |

The doubling of the unit cell leads to the splitting in the fermionic spectrum with the gap . The differences , , play the role of the order parameter for the modulated phase ( in the superfluid phase and in the supersolid phase). Coming from the set of equations (3), (34) and (35), we can write the equations for , , and separate the contributions of the first order. As a result, we obtain the condition of the appearance of nonzero solutions for , and . It can be shown that this condition coincides with the condition when the static density-density correlator diverges. Therefore our scheme for calculation of the density-density correlator in the RPA and the corresponding averages , and in the MFA is a self-consistent scheme.

At the numerical calculations of the density-density correlator we consider a three-dimensional case with a lattice constant and in our calculations we choose a half width of the fermionic band to be our energy scale (). First we investigate the uniform phase (). From the set of equations (3), (34) and (35) it follows that the solutions of these equations with can be realised when . Therefore at finite temperature we can consider the transition from the normal uniform nonsuperfluid phase (at low temperatures this is a Mott insulating phase) to the CDW phase for small values of the bosonic hopping parameter (). In figure 2(a), we plot lines of the instability with respect to the transition into the charge ordered phase at the fixed fermionic chemical potential (the case corresponds to the half-filling case) for the case of the nonsuperfluid phase (the bosonic concentration ). As seen in figure 2(a), the highest temperature of the instability is realised for the case of the chess-board phase with the wave vector . In figure 2(b) the lines of the instability for the case (when the system goes away from the half-filling case) are plotted. From figure 2(b) we observe that the incommensurate charge ordered phase with the wave vector has the highest temperature of the instability and the system undergoes the transition to the incommensurate modulated phase. It should be noted that the condition of the divergence of the static density-density correlator allows one to investigate the phase transitions of the second order only.

Now let us consider the transition to the supersolid phase. In figure 3 lines of instability of the superfluid phase with respect to the transition into the supersolid phase for the half-filling case , are depicted. As shown in figure 3, the transition to the supersolid phase with modulation wave vector is realised. We revealed that when the system goes away from the half-filling case and the supersolid phase with the modulation wave vector also has the highest temperature of the transition.

It should be emphasised that appearance of the CDW and supersolid phases is connected with the presence of the effective interaction between bosons which is formed due to the boson-fermion correlation. This interaction depends on the filling of the fermionic band.

Now we want to investigate the case of the chess-board phase in more detail. We use the equations for averages (3), (34), (35) and the expression for the grand canonical potential (3) to find thermodynamically stable states (in this part of our numerical calculations we use the semielliptic density of states, ). The phase transition lines and particle concentrations as functions of the bosonic chemical potential are shown in figure 4 and figure 5. The phase transition from the normal uniform nonsuperfluid to chess-board phase can be of the second or first order, see figure 4a. The existence of the phase transition of the first order leads to phase separation in the regime of the fixed concentrations into the NR and CDW phases. As shown in figure 5, similar picture is obtained for the transition from the superfluid to the supersolid phase and the transition from the superfluid to the supersolid phase can be of the first or second order depending on the system parameters.

In figure 6, we show the phase diagrams in the plane () at low temperature. As temperature increases, the regions of the existence of the CDW phase and the supersolid phase are possible for smaller parameter space and the first order phase transitions from the normal uniform nonsuperfluid (superfluid) into the CDW (SS) phases transforms into the second one. It should be noted that similar diagrams at were obtained in [17], but they did not reveal the possibility of the transition to the supersolid phase.

## 4 Conclusions

The phase transitions in the Bose-Fermi-Hubbard model at finite temperature has been considered in this work. We studied the hard-core limit and used pseudospin formalism. The thermodynamics of the model was investigated in the case of the weak boson-fermion interaction. The analytical expression for the density-density correlator has been obtained in the framework of the self-consistent scheme of the random phase approximation. The effective boson-boson interaction is formed due to the boson interaction with fermions, this effective interaction depends on the filling of the fermionic band. It is revealed that at small values of the bosonic tunnelling amplitude the system undergoes the phase transition from the uniform nonsuperfluid phase to the chess-board phase (the case of half filling of the fermion band) or to the incommensurate phase (when the system goes away from the half filling) at the lowering of the temperature. At increase of the bosonic hopping parameter the phase transition from the superfluid phase to the supersolid phase with a doubly modulated lattice period takes place (it should be noted that the presence of the supersolid phase at the weak boson-fermion interaction and zero-temperature was also established in [16] in the framework of a generalized dynamical mean field theory). The transition from the uniform to modulated phase can be of the first or second order, depending on the model parameters and temperature. The presence of the first order phase transition means that in the regime of the fixed fermionic concentrations the phase separation into the uniform and modulated phases is possible.

## References

## References

- [1] Jaksch D, Bruder C, Cirac J I, Gardiner C W and Zoller P 1998 Phys. Rev. Lett.81 3108
- [2] Greiner M, Mandel O, Esslinger T, Hänsch T W and Bloch I 2002 Nature 415 39
- [3] Günther K, Stöferle T, Moritz H, Köhl M and Esslinger T 2006 Phys. Rev. Lett.96 180402
- [4] Ospelkaus S, Ospelkaus C, Wille O, Succo M, Ernst P, Sengstock K and Bongs K 2006 Phys. Rev. Lett.96 180403
- [5] Albus A, Illuminati F and Eisert J 2003 Phys. Rev.A 68 023606
- [6] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 Phys. Rev.B 40 546
- [7] Sheshadri K, Krishnamurthy H R, Pandit R and Ramakrishnan T V 1993 Europhys. Lett. 22 257
- [8] Büchler H P and Blatter G 2003 Phys. Rev. Lett.91 130404
- [9] Fehrmann H, Baranov M, Lewenstein M and Santos L 2004 Optics Express 12 55
- [10] Cramer M, Eisert J and Illuminati F 2004 Phys. Rev. Lett.93 190405
- [11] Varney C N, Rousseau V G and Scalettar R T 2008 Phys. Rev.A 77 041608(R)
- [12] Lutchyn R M, Tewari S and Das Sarma S 2008 Phys. Rev.B 78 220504(R)
- [13] Refael G and Demler E 2008 Phys. Rev.B 77 144511
- [14] Iskin M and Freericks J K 2009 Phys. Rev.A 80 053623
- [15] Hébert F, Batrouni G G, Roy X and Rousseau V G 2008 Phys. Rev.B 78 184505
- [16] Titvinidze I, Snoek M and Hofstetter W 2008 Phys. Rev. Lett.100 100401
- [17] Mering A and Fleischhauer M 2010 Phys. Rev.A 81 011603(R)
- [18] Stashans A, Lunell S, Bergström R, Hagfeldt A and Lindquist S-E 1996 Phys. Rev.B 53 159
- [19] Koudriachova M V, Harrison N M and de Leeuw S W 2001 Phys. Rev. Lett.86 1275
- [20] Koudriachova M V, Harrison N M and de Leeuw S W 2002 Phys. Rev.B 65 235423
- [21] Wagemaker M, van de Krol R, Kentgens A P M, van Well A A and Mulder F M 2001 J. Am. Chem. Soc. 123 11454
- [22] Wagemaker M, Kearley G J, van Well A A, Mutka H and Mulder F M 2003 J. Am. Chem. Soc. 125 840
- [23] Mysakovych T S and Stasyuk I V 2007 J. Phys. Studies 11 195
- [24] Mysakovych T S, Krasnov V O and Stasyuk I V 2008 Condens. Matter Phys. 11 663
- [25] Stasyuk I V and Dublenych Yu I 2005 Phys. Rev.B 72 224209
- [26] Stasyuk I V 2007 Phase transitions in the pseudospin-electron model (Order, disorder and criticality. Advanced problems of phase transition theory vol 2) ed Yu Holovatch (Singapore: World Scientific)
- [27] Mahan G D 1976 Phys. Rev.B 14 780
- [28] Tomoyose T 1997 J. Phys. Soc. Japan66 2383
- [29] Stasyuk I V and Dulepa I R 2007 Condens. Matter Phys. 10 259
- [30] Izyumov Yu A and Skryabin Yu N 1989 Statistical Mechanics of Magnetically Ordered Systems (New York: Consultants Bureau)
- [31] Stasyuk I V and Mysakovych T S 2002 Condens. Matter Phys. 5 473
- [32] Stasyuk I V and Menchyshyn O 2009 J. Phys. Studies 13 4707