# Superfluid stiffness for the attractive Hubbard model on a honeycomb optical lattice

###### Abstract

In addition to the conventional contribution that is directly controlled by the single-particle energy spectrum, the superfluid phase stiffness of a two-component Fermi gas has a geometric contribution that is governed by the quantum metric of the honeycomb’s band structure. Here, we take both contributions into account, and construct phase diagrams for the critical superfluid transition temperature as a function of the chemical potential, particle filling, onsite interaction and next-nearest-neighbor hopping. Our theoretical approach is based on a self-consistent solution of the BCS mean-field theory for the stationary Cooper pairs and the universal BKT relation for the phase fluctuations.

###### pacs:

## I Introduction

Following the pioneering works by Peotta and Törmä on the origins of superfluidity in topologically nontrivial flat bands torma15 (), the deeper connection between some of the superfluid (SF) properties of a two-component Fermi gas and the quantum geometry of its non-interacting Bloch bands came as a complete surprise in recent years torma16 (); torma17a (); iskin17 (); iskin18a (); iskin18b (). It has been found on general grounds that the SF weight of a multi-band SF with a uniform order parameter can be separated into two distinct parts, depending on the physical mechanisms involved. While the real intraband processes are attributed to the conventional contribution, the virtual interband processes are attributed to the geometric one. Alternatively, in contrast to the conventional contribution that is solely controlled by the derivatives of the energy dispersions, the geometric one is also associated with the derivatives of the underlying Bloch wave functions torma17a (). For instance, unless the geometric interband contribution vanishes, superfluidity prevails in a flat band thanks to the presence of other flat or dispersive bands torma15 (); torma16 (). More recently, the root cause of this deeper connection has been identified as a mass-renormalization mechanism for the SF carriers, i.e., the quantum geometry governs not only the SF weight but also some other SF properties through renormalizing the effective mass of the two-body bound states and of Cooper pairs in general iskin17 (); iskin18b ().

Furthermore, in the particular cases of flat-band and two-band systems, the geometric contribution to the SF weight is simply controlled by the so-called quantum metric torma15 (); torma16 (); torma17a (); iskin17 (); iskin18a (). The quantum metric corresponds to the real part of the quantum geometric tensor, and its geometrical importance reveals itself as a measure of the quantum distance between nearby Bloch states provost80 (); berry89 (); thouless98 (). Note that the imaginary part of the quantum geometric tensor corresponds to the so-called Berry curvature, which is a distinct but related quantity associated with the emergent gauge field in momentum space, i.e., characterizing its quantum topology provost80 (); berry89 (); thouless98 (). Some of the two-band SFs that have already been analyzed in this context are the Haldane-Hubbard model torma17a (), Kane-Mele-Hubbard model torma17a (), time-reversal-invariant Hofstadter-Hubbard model torma15 (); iskin17 (), and the spin-orbit coupled Fermi gases iskin18a (). These works show clear signs that understanding the quantum metric effects on any one of these models may eventually have far reaching implications for a wide class of two-band SFs.

Motivated by these theoretical proposals as well as ongoing experimental efforts utilizing cold Fermions on various forms of honeycomb optical structures tarruell12 (); uehlinger12 (); polini13 (); jotzu14 (); flaschner16 (), here we calculate the critical SF transition temperature of the attractive Hubbard model on a two-dimensional honeycomb lattice for a large window of model parameters. Our theoretical approach is based on a self-consistent solution of the BCS mean-field theory for the stationary Cooper pairs and the universal BKT relation for the phase fluctuations, and we have two main goals. In addition to constructing the phase diagrams for the critical SF transition temperature, we plan to uncover the critical role played by the quantum geometry of the underlying band structure. For instance, while the highest attainable critical temperature is found to be around for the nearest-neighbor-hopping model, it increases quite rapidly with the inclusion of next-nearest-neighbor hoppings. In addition, the relative weight of the quantum metric contribution to the SF phase stiffness is found to be a non-monotonous function of the interaction strength, and it may reach beyond depending on the parameters. Thus, these findings arguably suggest that a SF Fermi gas that is loaded on a honeycomb lattice is one of the ideal platforms for studying quantum geometric effects with cold atoms.

The rest of the paper is organized as follows. The theoretical framework is presented in Sec. II, where we first discuss the honeycomb’s band structure and highlight the presence of Dirac cones in Sec. II.1, then introduce the BCS mean-field theory and derive the order-parameter and number equations in Sec. II.2, and then review the BKT relation and the SF stiffness in Sec. II.3. Having a complete set of self-consistency equations for determining the critical SF transition temperature, we present its numerical analysis in Secs. III and IV, and conclude the paper with our final remarks in Sec. V.

## Ii Theoretical Approach

The honeycomb lattice is a two-dimensional crystal structure with a hexagonal Bravais lattice and a two-site basis. In this paper, we denote its lattice spacing by , and choose and as the primitive lattice vectors for its Bravais lattice as shown in Fig. 1. The corresponding reciprocal lattice vectors and also form a hexagonal lattice in reciprocal space, leading to a first Brillouin zone that has the shape of a hexagon with side-length . Due to its two-site basis on a Bravais lattice, a honeycomb lattice gives rise to a two-band structure with important features for the single-particle problem. For the sake of completeness, let us first discuss its band structure and highlight the presence of Dirac cones, as they turn out to play critical roles in the many-body problem as well.

### ii.1 Band Structure

Within the tight-binding approximation, the single-particle Hamiltonian can be written as where the pseudo-spin denotes the two components of a Fermi gas, the index refers to a site in the hexagonal sublattice , the hopping parameter characterizes the tunneling amplitude from site to , and the operator () creates (annihilates) a particle on . In this paper, we set the nearest-neighbor (i.e., inter-sublattice) hopping parameter as our energy scale, and vary the next-nearest-neighbor (i.e., intra-sublattice) one accordingly.

Using the Fourier expansion of the creation and annihilation operators in the reciprocal () space, e.g., where is the number of sites in the hexagonal sublattice, a compact way to express this Hamiltonian is Here, the creation operator is in the form of a two-component sublattice spinor with the Hamiltonian density where is a unit matrix and is a vector of Pauli matrices. Note that while the diagonal element of is due solely to the next-nearest-neighbor hoppings, the off-diagonal element with and is due solely to the nearest-neighbor ones.

Thus, the single-particle energy eigenvalues are simply determined by where denotes the upper/lower band and reduces to . It can be shown that both energy bands exhibit a total of two Dirac cones that are equally distributed among the six corners of the first Brillouin zone. Since and at the tips of the cones, the density of single-particle states vanishes at , i.e., where the upper and lower cones touch each other. This analysis suggests that the low-temperature behavior of a Fermi gas that is loaded on a two-dimensional honeycomb lattice is directly controlled by the presence of these Dirac cones in the band structure. For instance, the Fermi gas shows a semi-metallic behavior that persists up to a critical interaction threshold depending on the temperature zhao06 (); lee09 ().

Our main goals in this work are twofold. We would like not only to construct the phase diagrams for the critical SF transition temperature of the attractive Hubbard model on a honeycomb lattice, but also to uncover the critical role played by the quantum geometry of the underlying band structure. These are achieved by adapting a self-consistent approach that is based on the simultaneous solution of the BCS mean-field theory for the stationary Cooper pairs and the universal BKT relation for the phase fluctuations.

### ii.2 BCS Mean-field Theory

For a two-component Fermi gas that is considered in this work, the attractive Hubbard model can be written as , where with takes the onsite intercomponent interactions into account. We treat the interaction term within the BCS mean-field approximation for pairing, and characterize various SF phases through the complex order parameter where denotes the thermal average. However, thanks to the time-reversal symmetry of , turns out to be uniform for a given sublattice, i.e., Furthermore, in order to fix the total number of particles in the thermal state, we include an additional term in , where is the chemical potential.

Similar to the single-particle problem, a compact way to express the mean-field Hamiltonian is

(1) |

where is a constant with denoting the trace over sublattices, is a four-component spinor operator, and with the Kronecker-delta is diagonal in the sublattice sector. Since turns out to be uniform for the entire lattice thanks to the inversion symmetry of the and sublattices, we take as real without losing generality. Combining this expression with the number equation we eventually obtain a set of self-consistency equations

(2) | ||||

(3) |

for and . Here, is the number of lattice sites, is the shifted dispersion, is a thermal factor with the Boltzmann constant and the temperature, is the quasi-particle energy spectrum, and is the total particle filling. Thus, we use Eqs. (2) and (3) to determine and for any given set of , , and parameters.

While the critical BCS transition temperature is simply determined by setting in Eqs. (2) and (3), the mean-field theory is known to give qualitatively reliable results for only. This is because growing phase fluctuations eventually break the mean-field approximation down in the limit, for which characterizes the pair formation temperature nsr85 (); randeria92 (). In particular to two dimensions, the SF phase coherence temperature is determined by the universal BKT relation, leading to a much lower result.

### ii.3 BKT Temperature and Superfluid Stiffness

Going beyond the mean-field theory, and including the phase fluctuations, the critical SF transition temperature is determined by the universal BKT relation through an analogy with the XY model b (); kt (); nk (). This approach has long been applied to the single-band SFs with great success denteneer93 (), and it has recently been generalized to the case of multi-band SFs with uniform order parameters torma16 (); torma17a (); iskin17 (). In the case of two-band SFs, one finds torma17a (); iskin18a ()

(4) | ||||

(5) | ||||

(6) |

where is the SF phase stiffness, is the area of the system, and is a thermal factor. Here, while the conventional contribution is of the usual single-band form accounting for real intraband processes, the geometric contribution is due to virtual interband processes that are directly controlled by the total quantum metric of the single-particle bands, i.e., with a unit vector.

In contrast to the standard BCS mean-field theory, it turns out that a self-consistent solution of Eqs. (2)-(6) for , and provides a qualitatively reliable description of the critical SF transition temperature for both and limits denteneer93 (). Even though this approach is still far from being quantitatively accurate in comparison to the numerically-exact ones lee09 (), its much simpler analytical construction provides considerable insight into the main features. For instance, Eq. (4) puts as the upper bound on in such a way that from below when , and that when . The latter limit can be shown by noting that and leading to a diagonal SF stiffness with Thus, the BKT relation (4) gives showing that, independently of its sign, increases for a given in the limit. Except for the and limits, the self-consistency equations are not analytically tractable in general, and we resort to numerical methods instead.

## Iii Numerical Results

Having introduced the theoretical framework, here we implement an iterative numerical approach to find fully self-consistent solutions for , and that satisfy all Eqs. (2)-(6) simultaneously for a given set of model parameters. This allows us not only to construct phase diagrams for the critical SF transition temperature, but also to uncover the critical role played by the quantum geometry of the underlying band structure. For this purpose, next we choose a set of exemplary ratios, and present the self-consistent results for and as a function of , and . Note that thanks to the apparent symmetry between the parameter sets and , our phase diagrams cover the entire parameter regime of the model Hamiltonian.

Let us first set and analyze the non-interacting limit. Since the lower single-particle band lies within the energy interval and the upper one lies within , the Fermi gas first fills the lower band as a function of increasing its Fermi energy up until . For this reason while region is denoted as a vacuum of particles with and region is denoted as a vacuum of holes (i.e., a band insulator) with , corresponds precisely to the half filling with . Thus, except for the symmetric point where the single-particle density of states vanishes like a semi-metal, the ground state of a non-interacting Fermi gas is normal when . All of these regions are clearly seen in Figs. 2 and 3.

Once the interactions are turned on, the BCS theory suggests that the normal region immediately transits into a SF with coinciding with in the limit. While such an exponential growth is clearly seen in the darker region in Figs. 2 and 3, our convergence scheme eventually fails in the white region when , where becomes comparable to our relative numerical accuracy (i.e., ) between two consecutive iterations. In Fig. 2, we note that the periphery of the white region nicely follows the general structure of , including its van Hove singularities at and . On the other hand, both the particle and hole vacuums as well as the semi-metal phase transit into a SF at finite interaction thresholds, beyond which grows as nearby the former regions and as nearby . For instance, we find that the critical interaction threshold decreases with for the semi-metal to SF phase transition at half filling. These are consistent with the known results in the literature where for zhao06 (); cichy18 (). We emphasize that, in contrast to the normal to SF transition boundary, our vacuum, insulator and semi-metal to SF transition boundaries are very accurate as vanish quite rapidly near .

In addition, for , Fig. 3 shows that a maximal value of is attainable at half filling when with and . However, we also find that higher critical temperatures may be achieved, respectively, with at fillings when . This is in agreement with our expectation that increases with independently of its sign in the limit. Indeed, our numerical results benchmark very well with this analytic expression in the regime of its validity. Note that our result is considerably higher than the maximal value reported in the literature for zhao06 ().

The discrepancy between our findings and the literature may well be caused by the use of a non-fully-self-consistent approach that is based on the numerical extraction of at zhao06 (). In addition, suspecting that the novel geometric contribution to the SF stiffness may partly be responsible for the apparent disagreement, we also present the relative weight of this contribution for the same parameters. For instance, for , reaches its maximum value at when , and it is at the location of the maximal . Thus, the geometric contribution is a non-monotonous function of , and it accounts for a sizeable fraction of the SF stiffness in general reaching beyond . In particular, Fig. 3 shows that its role becomes more and more (less and less) critical at lower (higher) fillings with decreasing .

In the limit, we may relate to the density and effective mass of the SF pairs through the identity where with the filling of condensed pairs iskin17 (). This leads to as the filling of SF pairs whose effective mass increases with but decreases with in the limit. For completeness, the fraction of condensed particles is also shown in Figs. 2 and 3 for the parameter regimes of interest. In comparison to the half-filling case where half of the pairs or holes may at most be condensed with in the limit, all of the particle (hole) pairs are condensed with () in the low particle (hole) filling () limit.

## Iv Discussion

In order to provide a better contextualization of our results, here we first calculate in a fully-self-consistent manner both with and without the geometric contribution , and benchmark these results directly with those of the literature. For this purpose, we set in Fig. 4, and present the resultant phase diagrams for precisely the same parameter window as the one that is shown in Ref. zhao06 ().

This figure clearly illustrates that exclusion of the contribution from the universal BKT relation leads to a substantial reduction of for .This is in accordance with the discussion given above in Sec. II.3, where is determined by the vanishing () of the BCS mean-field for . In addition, since the latter phase diagram turns out to be quantitatively similar to that of Ref. zhao06 (), we speculate that the geometric contribution may not be taken fully into account by their approach. This is because while the self-consistent equations that are solved numerically for the pairing order parameter and particle filling are exactly the same in both works, there is one important difference in the way the SF stiffness is calculated. That is our analytic expression for the SF stiffness is derived under the linear response theory, but Ref. zhao06 () extracts it numerically from the dispersion of the Goldstone mode which in turn is derived by considering the Gaussian fluctuations of the order parameter on top of the BCS ground state. Given our detailed benchmark, we believe these two approaches are equivalent for the region but not away from it for the honeycomb lattice. As both approaches are routinely used to evaluate the SF stiffness and the related critical SF transition temperature, further work is needed to assess which method is more reliable and accurate, and whether they could be reconciled to some extent. In addition, a more detailed comparison with the state of the art unbiased numerical method is also highly desirable.

## V Conclusions

In summary, here we calculated the critical SF transition temperature of the attractive Hubbard model on a two-dimensional honeycomb lattice via a theoretical approach that is based on a self-consistent solution of the BCS mean-field theory for the stationary Cooper pairs and the universal BKT relation for the phase fluctuations. For instance, we found that the highest attainable is around for the nearest-neighbor-hopping model, and that it increases quite rapidly with the inclusion of next-nearest-neighbor hoppings. In addition to the construction of the phase diagrams for a large window of model parameters, we also uncovered the critical role played by the quantum geometry of the underlying band structure. In particular, we found that the relative weight of the quantum metric contribution to the SF phase stiffness is a non-monotonous function of the interaction strength, and that it is generally far from being negligible reaching beyond . These findings arguably suggest that a SF Fermi gas that is loaded on a honeycomb lattice tarruell12 (); uehlinger12 (); polini13 (); jotzu14 (); flaschner16 () is one of the ideal platforms for studying quantum geometric effects with cold atoms. The possible outcomes of such a realization would clearly have a broader impact in solid-state, condensed-matter and some other physics communities, given the modern surge of interest in the quantum topological and/or quantum geometrical concepts in general.

This line of work offers many extensions for future research. For instance, since the SF stiffness is determined by the ratio of the density and the mass of the SF carriers, we expect sizeable geometric contributions for those observables (e.g., sound velocity) that have explicit dependence on the mass of the two-body bound states or of the Cooper pairs on general grounds iskin18b (). Given our findings for the geometric SF stiffness, we expect these observables to have a non-monotonous dependence on the interaction strength as well. This distinguishes the honeycomb lattice from the square-like Bravais lattices, for which the corresponding ground-state observables are known to evolve monotonously in the usual BCS-BEC crossover problem nsr85 (); randeria92 (). In addition to the honeycomb system, we are aware of other two-band systems with a non-trivial quantum geometry, exhibiting similar geometric effects. For instance, the Haldane-Hubbard model torma17a (), Kane-Mele-Hubbard model torma17a (), time-reversal-invariant Hofstadter-Hubbard model torma15 (); iskin17 (), and the spin-orbit coupled Fermi gases iskin18a () are described by the single-particle and mean-field Hamiltonians that are of exactly the same form as the ones considered in this paper. Thus, there is no doubt that understanding the quantum metric effects on a honeycomb lattice may eventually have far reaching implications for a wide-class of two-band SFs.

###### Acknowledgements.

The author acknowledges financial support from TÜBİTAK.## References

- (1) S. Peotta and P. Törmä, “Superfluidity in topologically nontrivial flat bands”, Nat. Com. 6, 8944 (2015).
- (2) A. Julku, S. Peotta, T. I. Vanhala, D.-H. Kim, and P. Törmä, “Geometric Origin of Superfluidity in the Lieb-Lattice Flat Band”, Phys. Rev. Lett. 117, 045303 (2016).
- (3) L. Liang, T. I. Vanhala, S. Peotta, T. Siro, A. Harju, and P. Törmä, “Band geometry, Berry curvature, and superfluid weight”, Phys. Rev. B 95, 024515 (2017).
- (4) M. Iskin, “Berezinskii-Kosterlitz-Thouless transition in the time-reversal-symmetric Hofstadter-Hubbard model”, Phys. Rev. A 97, 013618 (2018).
- (5) M. Iskin, “Exposing the quantum geometry of spin-orbit coupled Fermi superfluids”, Phys. Rev. A 97, 063625 (2018).
- (6) M. Iskin, “Quantum metric contribution to the pair mass in spin-orbit-coupled Fermi superfluids”, Phys. Rev. A 97, 033625 (2018).
- (7) J. P. Provost and G. Vallee, “Riemannian structure on manifolds of quantum states”, Commun. Math. Phys. 76, 289 (1980).
- (8) M. V. Berry, “The quantum phase, five years after in Geometric Phases in Physics”, edited by A. Shapere and F. Wilczek (World Scientific, Singapore, 1989).
- (9) D. J. Thouless, “Topological Quantum Numbers in Nonrelativistic Physics”, (World Scientific, Singapore,1998)
- (10) L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, “Creating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb lattice”, Nature (London) 483, 302 (2012).
- (11) T. Uehlinger, G. Jotzu, M. Messer, D. Greif, W. Hofstetter, U. Bissbort, and T. Esslinger, “Artificial Graphene with Tunable Interactions”, Phys. Rev. Lett. 111, 185307 (2013).
- (12) M. Polini, F. Guinea, M. Lewenstein, H. C. Manoharan, and V. Pellegrini, “Artificial honeycomb lattices for electrons, atoms and photons”, Nat. Nano. 8, 625 (2013).
- (13) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, “Experimental realization of the topological Haldane model with ultracold fermions”, Nature (London) 515, 237 (2014).
- (14) N. Fläschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. Lühmann, K. Sengstock, and C. Weitenberg, “Experimental reconstruction of the Berry curvature in a Floquet Bloch band”, Science 352, 1091 (2016).
- (15) E. Zhao and A. Paramekanti, “BCS-BEC crossover on the two-dimensional honeycomb lattice”, Phys. Rev. Lett. 97, 230404 (2006).
- (16) K. L. Lee, K. Bouadim, G. G. Batrouni, F. Hébert, R. T. Scalettar, C. Miniatura, and B. Grémaud, “Attractive Hubbard model on a honeycomb lattice: Quantum Monte Carlo study”, Phys. Rev. B 80, 245118 (2009).
- (17) P. Noziéres and S. Schmitt-Rink, “Bose condensation in an attractive fermion gas: From weak to strong coupling superconductivity”, J. Low Temp. Phys. 59, 195 (1985).
- (18) L. Belkhir and M. Randeria, “Collective excitations and the crossover from Cooper pairs to composite bosons in the attractive Hubbard model”, Phys. Rev. B 45, 5087(R) (1992).
- (19) V. L. Berezinskii, “Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group I. classical systems”, JETP 32, 493 (1971).
- (20) J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transitions in two-dimensional systems”, J. Phys. C: Solid State Phys. 6, 1181 (1973).
- (21) D. R. Nelson and J. M. Kosterlitz, “Universal Jump in the Superfluid Density of Two-Dimensional Superfluids”, Phys. Rev. Lett. 39, 1201 (1977).
- (22) P. J. H. Denteneer, G. An, and J. M. J. van Leeuwen, “Helicity modulus in the two-dimensional Hubbard model”, Phys. Rev. B 47, 6256 (1993).
- (23) A. Cichy and A. Ptok, “Reentrant Fulde-Ferrell-Larkin-Ovchinnikov superfluidity in the honeycomb lattice”, Phys. Rev. A 97, 053619 (2018).