Emergent Spin Liquids in the Hubbard Model on the Anisotropic Honeycomb Lattice
We study the repulsive Hubbard model on an anisotropic honeycomb lattice within a mean-field and a slave-rotor treatment. In addition to the known semi-metallic and band-insulating phases, obtained for very weak interactions, and the anti-ferromagnetic phase at large couplings, various insulating spin-liquid phases develop at intermediate couplings. Whereas some of these spin liquids have gapless spinon excitations, a gapped one occupies a large region of the phase diagram and becomes the predominant phase for large hopping anisotropies. This phase can be understood in terms of weakly-coupled strongly dimerized states.
Metal-insulator transitions and other electronic transitions Lattice fermion models Atoms in optical lattices
One of the most salient features of graphene, an atomically thin graphite sheet with carbon atoms arranged in a honeycomb lattice, is certainly the relativistic character of its low-energy electronic excitations . An important issue of present research efforts is a deeper understanding of the role of interactions in this novel two-dimensional material and of the persistance of the semi-metallic (SM) phase when such interactions become relevant. As a consequence of the poor screening properties of electrons in the SM phase, the Coulomb interaction remains long-range and may be described in terms of an effective fine-structure constant that turns out to be density-independent, (for recent reviews on electronic interactions in graphene, see Refs.  and ). Here is the Fermi velocity in terms of the speed of light , and is the dielectric function of the surrounding medium. Whereas lattice-gauge theories predict a flow to strong coupling above some critical value , even suspended graphene seems to be weakly correlated, with a stable SM phase . A reason for this effective flow to weak coupling may be an intrinsic dielectric constant in graphene due to virtual interband excitations . Indeed, recent renormalization-group studies confirm this picture of weakly-interacting electrons in graphene .
A perhaps more promising system for the study of the interplay between strong (short-range) correlations and the relativistic character of Dirac fermions may well be a gas of cold fermionic atoms trapped in an optical honeycomb potential [8, 9]. As compared to graphene, such a system has several advantages. First, neutral atoms residing on the same site exhibit short-range interactions which can often be tuned over orders of magnitude and turned repulsive or attractive by using Feshbach resonances . Indeed, recent experiments have proven the feasibility of implementing optical honeycomb lattices and probing interaction physics in the context of bosonic Rb atoms . Second, the hopping rates between neighboring sites are easily controlled by the laser configuration and beam intensities as both determine the depth and position of the optical potential wells. Rather moderate changes in the laser intensities can significantly imbalance the tunneling rates, leading to situations ranging from weakly-coupled zig-zag linear chains to weakly-coupled dimers . This situation needs to be contrasted to graphene, where unphysically large lattice distortions are required to obtain the limit , where novel physical phenomena are expected .
Starting from the repulsive fermionic Hubbard model (RFHM) for spin-1/2 particles with onsite interactions and identical nearest-neighbor hopping rates , mean-field calculations for a half-filled lattice and zero temperature predict
a SM phase at small values of and an anti-ferromagnetic (AF) phase developing above [14, 15]. Quantum Monte Carlo (QMC) calculations confirm this SM-AF transition but at a higher transition point [14, 16] while dynamical-mean-field estimates yield even larger values, . Most saliently, recent QMC investigations have revealed an intriguing insulating spin-liquid (SL) phase, with localized charges but no spin ordering, that emerges below the AF transition  and that may be related to the exotic Mott insulator identified in slave-rotor studies [18, 19, 20]. The transition points derived from the slave-rotor theory  also occur at globally smaller values of than in QMC calculations.
Here, we study the RFHM on a half-filled honeycomb lattice where one of the nearest-neighbor hopping parameters is larger than the other two . Such a lattice reveals an astonishingly rich phase diagram (see Fig. 1), in which SL phases predominate at large values of . At , this system develops a topological phase transition between a semi-metal and a band insulator (BI). Indeed, at , the two Dirac points responsible for the SM phase merge and eventually disappear while a band gap opens [8, 22]. These topological properties are prominently reflected in the SL phase, where the spin excitations acquire a gap as a function of renormalized hopping parameters. Finally, we show that the SL phase finds a compelling interpretation in terms of weakly-coupled dimer states in the limit .
As the honeycomb lattice is made of two shifted identical triangular sublattices and , the kinetic term of the RFHM, , reads
where and annihilate (create) a fermion with spin on the sites of the and sublattices, respectively. It also distinguishes the hopping amplitude , chosen along neighboring sites linked by the vector , which differs from the other two chosen along neighboring sites linked by or by , being the lattice constant [see Fig. 2(a)]. At half-filling, the onsite interaction term reads
with , the summation running over both sublattices and or being the corresponding number operators. Here we consider a balanced population between spin- and spin-. When , the Hamiltonian is readily diagonalized in reciprocal space, and one obtains the two energy bands , in terms of the weighted sum of phase factors with and . When , one recovers the usual Dirac points located at K and K, as depicted in Fig. 2(b). The Dirac points move towards the point M as is increased, where they finally merge when [8, 22]. At half-filling, this topological phase transition separates a SM phase (for ), with two Dirac points, from a BI (for ), the insulating gap opening at the point M.
For this half-filled lattice, one expects a Mott insulating (MI) phase at large in which each lattice site is exactly occupied by one atom. As a consequence of the superexchange interaction , the spins of the particles are further frozen in an AF Néel state that breaks the lattice inversion symmetry, with e.g. spin- atoms on the sublattice and spin- on the sublattice. This state is described by the staggered mean-field order parameter, . The mean-field approximation in the interaction term leads to a quadratic Hamiltonian which is diagonalized through a Bogoliubov-Valatin transformation. The corresponding excitation spectrum is . The Mott gap is found to be and is therefore intimitely related to the AF order parameter, i.e. it opens as soon as the AF order sets in. In the thermodynamic limit, the self-consistent gap equation at zero temperature then reads
where is the area of the first Brillouin zone (FBZ). The transition line delineating the AF phase (blue solid line in Fig. 1) is obtained from Eq. (3) by setting . Notice that the value , obtained at , agrees with former mean-field calculations [14, 15]. The critical value is however shifted to larger values when increases. This shift may be understood from a simple scaling argument when considering the full band width of the non-interacting case. Since all particles are localized in the AF phase, one needs to consider states with energies up to . The critical value for the isotropic case may now be expressed as . If we consider this value to remain constant when varying , one obtains , a relation that describes to great accuracy the transition line for (Fig. 1). At larger values of , the mean-field Mott gap competes with the band insulator gap , and the transition line may be understood qualitatively as the line where both gaps are equal. This yields the asymptotic behavior in the large -limit (dashed line in Fig. 1).
In order to decouple the Mott transition from the AF phase and to study possible intermediate phases, more sophisticated methods than a simple mean-field treatment are required. Apart from QMC calculations , an intermediate MI spin-liquid phase has been identified in the honeycomb lattice within a slave-rotor treatment [18, 19, 20], where the fermion operators and are viewed as products of two auxiliary degrees of freedom, . We adopt, here, a U(1) slave-rotor treatment that is expected to provide qualitatively correct results within the mean-field approximation. However, if one aims at an effective low-energy theory for the spinon degrees of freedom, one needs to take into account the coupling to an SU(2) gauge field as discussed in Ref. . The bosonic “rotor” field is conjugate to the total charge at site , described by the angular momentum , and the fermion operator carries the spin (“spinon”). As this procedure artificially enlarges the Hilbert space, double counting needs to be cured by imposing the constraint
at each site . Whereas the rotor and spinon fields are coupled via the kinetic Hamiltonian , the interaction term is described solely in terms of the rotor degrees of freedom, , the particle-hole-symmetric form of which also renders the chemical potential . Following Ref. , the rotor and spinon degrees of freedom may now be decoupled within a mean-field treatment by defining the averages and , where and and are nearest neighbors connected by for the primed averages and by otherwise. Notice that the mean-field parameters along are assumed to be equal, hence the corresponding symmetries are not broken here. The decoupled mean-field Hamiltonian may then be written as , with
where is a local Lagrange multiplier ensuring the constraint (4). At half-filling, particle-hole symmetry imposes that . The Mott transition may now be interpreted in terms of rotor condensation. Indeed, in the rotor-condensed phase, the phase is fixed and the number of particles (or the angular momentum ) therefore fluctuates on the lattice sites. This corresponds to the SM phase for or the BI for . In the MI phase, however, there is no spin ordering since the spinon Hamiltonian has no interaction term.
Most saliently, can be mapped onto through , such that the low-energy spinon excitations are governed by the same topological properties as the non-interacting Dirac fermions. In particular, similarly to , we find the weighted sums of phase factors
Therefore, there exists a critical value
that separates a SL with gapless spinon excitations, around two distinct Dirac points at the Fermi level (for ), from a gapped SL phase (for ). The spinon dispersion calculated from the Hamiltonian is depicted in Fig. 3 for two characteristic sets of parameters.
To obtain the Mott transition line, one needs to self-consistently solve for , , , and . This is achieved with the help of the imaginary-time action
where is the inverse temperature and the Boltzmann constant. We next rewrite (9) in terms of the fields and constrain the normalization via a Lagrange multiplier . The resulting Green’s functions for the fields and are
where and are the fermionic and bosonic Matsubara frequencies, respectively. A change of is performed in equation (10a), in order to preserve the correct atomic limit . Based on the form of the rotor Green’s function, one may define the charge gap as
where is the maximum of over the FBZ. The normalization of the -field yields the equation 
We consider the rotor-disordered and rotor-condensed phases separately. In the former case, we perform the sum over Matsubara frequencies in Eq. (12) and consider the zero temperature limit. We find
which, together with the equations for the s, completes the set of self-consistency equations for the s and for the rotor-disordered (MI) phase. For the rotor-condensed (SM/BI) phase, as in the case of a normal Bose-Einstein condensate, a macroscopic fraction, namely per lattice site, of the particles occupies the state, which corresponds to the lowest energy . The chemical potential () is equal to this lowest energy. The corresponding equation becomes
which, instead of Eq. (13), completes the set of self-consistent equations for the s and for the rotor-condensed (SM/BI) phase.
In the rotor-disordered (MI) phase, one has () and , while in the rotor-condensed (SM/BI) phase, one has and . A second-order transition line is defined by and . As we show below, part of the phase transition curve that we obtained is of first order, indicated by a jump of or from some finite values to zero. For the rotor-disordered (MI) phase, one finds three distinct solutions with different types of spinon excitations:
The chain solution has vanishing mean-field parameters on the horizontal bonds, i.e. , and thus describes a system composed of decoupled vertical zig-zag chains [corresponding to the bonds with a hopping in Fig. 2(a)]. A self-consistent calculation yields , while and are functions of . The associated SL phase in this case is gapless, since .
The dimer solution has , and hence describes a system composed of decoupled dimers on the horizontal bonds of the honeycomb lattice [see Fig. 2(a)]. Self-consistency requires that and . One may easily verify that for this solution, independent of . This relation also fixes the value of . This SL is gapped in the spinon channel, since .
A third (more general) solution in the rotor-disordered SL phase is continuously connected to the rotor-condensed SM or BI phase, and it therefore separated from the latter by a second-order phase transition. This solution may be obtained by starting in the rotor-condensed phase with a sufficiently small value of , where the mean-field parameters are obtained numerically. These parameters serve as the starting point for the next calculation step (i.e. for , with an incremental energy step ). One thus obtains a continuous line of solution in parameter space. The condensate fraction decreases with increasing on-site repulsion , and the rotor-disordered MI phase is obtained when . Furthermore, the charge gap increases from zero continuously when is increased, as one expects for a second-order phase transition.
In order to identify which of the above-mentioned solutions is chosen by the system for a particular set of parameters (), we have performed a free-energy analysis of the three solutions. The free energy per lattice site reads in terms of the partition function and the number of lattice sites , and one has, in the zero- limit, in which the free energy becomes the internal energy
|is the free energy of the fermionic part, and|
that of the bosonic part. Note that in the condensed phase, . The last term in equation (15a) accounts for the correction to the free energy due to the dynamically unimportant constants, which have been ignored so far.
For a particular pair of , the true ground state is the one with the smallest . The transition curves resulting from this free energy analysis are plotted in Fig. 1 (lower thicker green and dotted line), and Fig. 4 for the vicinity of the isotropic limit (). One finds that, for , the system undergoes a second-order phase transition (within the third solution) from the SM phase to a gapless SL phase (GLSL in Fig. 1) upon increase of ; for , the system first experiences the above-mentioned second-order phase transition from a SM to the GLSL (third solution), and then a first-order phase transition to another gapless one of the chain solution [GLSL (Chain)]; for , the second-order phase transition disappears and the system jumps directly from a SM phase to the GLSL (Chain) via a first order phase transition. Finally, for , the system undergoes a first-order phase transition from the SM/BI phase to a gapped spin liquid of the dimer solution [GDSL (Dimer)]. The AF transition, which cuts the GLSL-GLSL (Chain) transition at , but which is otherwise well above the SM/BI-SL transitions described here, is not treated within the current slave-rotor description. In principle, this AF transition may be described within a slave-rotor theory by reintroducing spin-charge correlations into the Hamiltonians (5b) .
For the isotropic case , our result of the second-order transition point agrees with Refs. [18, 20]. In this limit, the free-energy analysis indicates that the gapless SL of the third solution is the lowest while the gapped SL of the dimer solution is the highest in energy among the three solutions. This is in contrast to the QMC calculations , where a gapped SL was identified as the true ground state. However, only an upper bound of the free energy of the dimer solution is provided by our free energy analysis. This is especially the case in the isotropic limit, where a discrimination among the three hopping parameters is unphysical. Indeed, possible kinetic dimer terms (e.g. resonating dimer moves around a hexagon) could lower the energy of the gapped SL of the dimer solution . At a general value of , free energies are modified by these terms such that they are adiabatically connected to the (also modified) isotropic case. Similarly, the gapless chain solution may be unstable to an infinitesimal interchain coupling that, in the context of the square lattice, is known to open a spin gap . An analysis of kinetic dimer terms and the stability of the chain solution is, however, beyond the scope of this paper.
The predominating SL phase in the limit is described by the dimer solution, where both and are zero. This in effect is equivalent to setting , such that only sites connected by horizontal bonds are coupled by tunneling (Fig. 2) and the whole system reduces to a set of decoupled dimers. It is therefore sufficient to solve the Hubbard Hamiltonian for two interacting fermions occupying two lattice sites coupled by the hopping amplitude . The two-particle Hilbert space obviously reduces to six states as the Pauli principle forbids fermions with identical spin to occupy a same site. By the same token, states where fermions with identical spin occupy different sites are obviously zero-energy eigenstates as hopping is Pauli-blocked. The relevant two-particle Hilbert space is thus spanned by , , , and , where the first entry describes the occupancy of the site and the second one that of the site. The diagonalization of the Hubbard Hamitonian then yields the ground state
where is the dimer singlet and with . The first excited state is the triplet state and is separated from by the gap . The ground state is therefore essentially the singlet state , with an admixture of states with double occupancy the weight of which () vanishes in the large- limit, where the gap is dominated by the exchange energy . Because remains the ground state over the whole range of values of , there is just one single phase that continuously connects the BI to the GDSL phase. Notice that, within the mean-field U(1) slave-rotor treatment, there is not such a continuous connection because can only take the values (for the BI) and (for the GDSL), double occupancy of a single site being ruled out in the MI phase. Furthermore, only at the asymptotic point, , do the singlet and the triplet states become degenerate, and the AF state, which requires a degeneracy (and a superposition) of and , may be formed. In the intermediate region, the gap protects the ground state which is only marginally perturbed by a small inter-dimer coupling mediated by . Within this simple dimer picture, one may therefore understand the gapped SL [GDSL (Dimer)] in Fig. 1 as being adiabatically connected to the state .
In conclusion, we have investigated the repulsive fermionic Hubbard model on the anisotropic honeycomb lattice. This system could be experimentally realized by loading ultracold fermions at half-filling in a honeycomb optical lattice. Beside the SM, BI, and AF phases, which are readily obtained at the mean-field level, we have used a slave-rotor description to show the emergence of various SL phases. Two gapless SL phases may be found from a self-consistent mean-field treatment of the slave-rotor theory: a chain solution that consists of essentially decoupled zig-zag chains and a phase that is connected continuously to the small- SM/BI phase. The latter is found in an intermediate coupling regime for the isotropic case, . Most saliently, a gapped dimer SL phase may be stabilized at rather low hopping anisotropies, , and is adiabatically connected to the large- limit, where the system may be described in terms of essentially decoupled dimer states. One may speculate that the dimer picture, if modified by an additional kinetic term that allows for dimer moves, yields insight also into the physical properties of the system in the isotropic limit .
Acknowledgements.We thank A. H. Castro Neto, F. Crépin, B.-G. Englert, J.-N. Fuchs, M. Gabay, N. Laflorencie, K. Le Hur, G. Montambaux, F. Piéchon, and M. Rozenberg for fruitful discussions. We acknowledge support from the France-Singapore Merlion program (CNOUS grant 200960 and FermiCold 2.01.09) and the CNRS PICS 4159 (France). The Centre for Quantum Technologies is a Research Centre of Excellence funded by the Ministry of Education and the National Research Foundation of Singapore.
- A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- M. O. Goerbig, (preprint) arXiv:1004.3396.
- V. N. Kotov, B. Uchoa, V. M. Pereira, A. H. Castro Neto, F. Guinea, (preprint) arXiv:1012.3484.
- J. E. Drut and T. A. Lähde, Phys. Rev. Lett. 102, 026802 (2009); Phys. Rev. B 79, 165425 (2009).
- X. Du, I. Skachko, A. Barker, Eva Y. Andrei, Nature Nanotechnology 3, 491-495 (2008).
- J. González, F. Guinea, and M. A. H. Vozmediano, Phys. Rev. B 59, R2474 (1999).
- I. F. Herbut, V. Jurićić, and B. Roy, Phys. Rev. B 79, 085116 (2009).
- S.-L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev. Lett. 98, 260402 (2007).
- K. L. Lee, B. Grémaud, R. Han, B.-G. Englert, and Ch. Miniatura, Phys. Rev. A 80, 043411 (2009).
- T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, Phys. Rev. Lett. 93, 050401 (2004).
- P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nat. Phys., doi:10.1038/nphys1916.
- C. Lee, X. Wei, J. K. Kysar, and J. Hone, Science 321, 385 (2008); M. O. Goerbig, J.-N. Fuchs, G. Montambaux, and F. Piéchon, Phys. Rev. B 78, 045415 (2008); V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres, Phys. Rev. B 80, 045401 (2009).
- P. Dietl, F. Piéchon, and G. Montambaux, Phys. Rev. Lett 98, 236405 (2008); B. Wunsch, F. Guinea, and F. Sols, New J. Phys. 10 103027 (2008).
- S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992).
- N. M. R. Peres, M. A. N. Araújo, and Daniel Bozi, Phys. Rev. B 70, 195122 (2004).
- Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).
- S. A. Jafari, Eur. Phys. J. B 68, 537 (2009); W. Wu, Y.-H. Chen, H.-S. Tao, N.-H. Tong and W.-M. Liu, arXiv:1005.2043v2.
- S.-S. Lee and P. A. Lee, Phys. Rev. Lett 95, 036403 (2005).
- M. Hermele, Phys. Rev. B 76, 035125 (2007).
- S. Rachel and K. Le Hur, Phys. Rev. B 82, 075106 (2010).
- S. Florens and A. Georges, Phys. Rev. B 66, 165111 (2002), ibid. 70, 035114 (2004); for a pedagogical review of this technique, see E. Zhao and A. Paramekanti, Phys. Rev. B 76, 195101 (2007).
- G. Montambaux, F. Piéchon, J.-N. Fuchs, and M. O. Goerbig, Phys. Rev. B 80, 153412 (2009); Europhys. J. B, 72, 509 (2009).
- See e.g. E. Zhao and A. Paramekanti, Phys. Rev. B 76, 195101 (2007).
- R. Moessner, S. Sondhi, and P. Chandra, Phys. Rev. B 64, 144416 (2001).
- For a review, see I. Affleck, J. Phys.:Condens. Matt. 1, 3047 (1989).