Interacting fermions on the honeycomb bilayer: from weak to strong coupling

Interacting fermions on the honeycomb bilayer: from weak to strong coupling

Oskar Vafek National High Magnetic Field Laboratory and Department of Physics,
Florida State University, Tallahasse, Florida 32306, USA
September 16, 2019

Many-body instabilities of the half-filled honeycomb bilayer are studied using weak coupling renormalization group as well as strong coupling expansion. For spinless fermions and assuming parabolic degeneracy, there are 4-independent four-fermion contact couplings. While the dominant instability depends on the microscopic values of the couplings, the broken symmetry state is typically a gapped insulator with either broken inversion symmetry or broken time reversal symmetry, with a quantized anomalous Hall effect. Under certain conditions, the dominant instability may appear in the particle-particle (pairing) channel. For some non-generic fine-tuned initial conditions, weak coupling RG trajectories flow into the non-interacting fixed point, although generally we find runaway flows which we associate with ordering tendencies. Additionally, a tight binding model with nearest neighbor hopping and nearest neighbor repulsion is studied in weak and strong couplings and in each regime a gapped phase with inversion symmetry breaking is found. In the strong coupling limit, the ground state wavefunction is constructed for vanishing in-plane hopping but finite inter-plane hopping, which explicitly displays the broken inversion symmetry and a finite difference between the number of particles on the two layers. Finally, we discuss the spin-1/2 case and use Fierz identities to show that the number of independent 4-fermion contact couplings is 9. The corresponding RG equations in the spin-1/2 case are also presented, and used to show that, just as in strong coupling, the most dominant weak coupling instability of the repulsive Hubbard model (at half-filling) is an anti-ferromagnet.

I Introduction

The problem of interacting fermions on the stacked honeycomb bilayer at half-filling has attracted attention due to a confluence of several factors. First, purely on theoretical grounds, in its simplest form with the nearest neighbor hoping only, the tight-binding approximation gives rise to a band structure with two bands touching quadratically at the Fermi levelMcCann and Fal’ko (2006); Castro Neto et al. (2009) near two non-equivalent points in the Brillouin zone, and . Even at the non-interacting level, such quadratic degeneracy gives rise to logarithmically divergent susceptibilitiesCastro Neto et al. (2009); Vafek and Yang (2010) in several channels as temperature, or frequency, are taken to zeroZhang et al. (2010); Nandkishore and Levitov (2010); Lemonik et al. (2010). As a result, some form of spontaneous symmetry breaking is expected at finite temperature upon inclusion of even weak interactionsMin et al. (2008); Castro Neto et al. (2009); Sun et al. (2009); Vafek and Yang (2010); Zhang et al. (2010); Nandkishore and Levitov (2010); Lemonik et al. (2010). And while fine-tuning is necessary to achieve such band-structure, in that (with the exception of square checkerboard and Kagome lattices studied in Ref.Sun et al. (2009)) inclusion of trigonal warping termsMcCann and Fal’ko (2006) eventually gives rise to four Dirac fermions at each -point, non-interacting susceptibilities may be sufficiently enhanced that many-body instabilities appear, albeit at finite coupling strength. In this sense, the A-B stacked honeycomb bilayer problem is another example of the observation that there are no generic weak coupling particle-hole instabilitiesRaghu et al. (2010). Rather, fine-tuning, in the form of nesting for example, is necessary to bring the strong coupling physics down to weak coupling. If we are interested in accessing the symmetry-broking phases in the particle-hole channel, as we are in this case, then fine tuning is a small price to pay for this access, made available within perturbative RG. Second, the isolation of graphene bilayers and the experimental ability to perform, for example, electricalNovoselov et al. (2006); Feldman et al. (2009); Zhao et al. (2010), angle resolved photoemissionOhta et al. (2006), Raman spectroscopyYan et al. (2008) or infra-redZhang et al. (2009) measurements, while controlling the gate voltage through the neutrality point, gives rise to the opportunity to test such theoretical expectations in a reasonably well controlled physical setting. In addition, the technological promise of this material fuels further need to understand its electronic structure and with it the many-body interactions. Finally, the problem of interacting fermions on the AB-stacked honeycomb bilayer may soon be realized in cold atom optical lattices, where the theory may also be tested.

The issue of band-structure fine-tuning notwithstanding, the type of leading instability in a graphene bilayer (with spin fermions) has been a subject of debate as well. A mean-field approach has been used to argue for an insulating state with broken inversion symmetryMin et al. (2008). A similar approach has also been argued to lead to trivial gapped insulating phasesNandkishore and Levitov (2010) as well as to an anomalous quantum Hall phaseNandkishore and Levitov (2010). On the other hand, the leading weak coupling instability can be analyzed without resorting to uncontrolled approximations by using weak coupling renormalization group. This approach was used in Ref.Vafek and Yang (2010) where a nematic phase was found to be the dominant instability within the model studied. Such instability was subsequently also argued for in Ref. Lemonik et al. (2010). On the other hand, an inversion symmetry breaking insulating phase has been claimed in Ref. Zhang et al. (2010).

To determine what type of broken symmetry state is preferred in the case of spinless fermions, we perform weak coupling RG analysis by studying the flow of independent symmetry allowed short-range interactions. We find that generically, depending on the initial values of the 4-fermion contact couplings, the system flows into a gapped phase with either broken inversion symmetry and a finite difference between the total number of particles on the two layers, or broken time reversal symmetry. The former state was not found to be preferred in the model for spin- fermions studied in Ref.Vafek and Yang (2010) (where the nematic state was found to dominate), but an example of the latter state corresponded to one of the fixed points found therein. In particular, for the spinless case studied here, we find that a gapped state with anomalous (zero B-field) quantum Hall conductivity has the most divergent susceptibility for a range of initial couplings as determined by the (right) sink of the RG trajectories shown in Fig.(3). While non-generic, we also specify special conditions under which the interacting model flows back to the non-interacting fixed point.

In addition, we analyze the specific microscopic model with nearest neighbor hopping(s) (and ) and nearest neighbor repulsion in both the weak coupling RG and in strong coupling. In both regimes we find the (trivial) insulating phase with broken inversion symmetry to dominate. As discussed in more detail below, in weak coupling the RG flow tends to the left sink shown in Fig.(3), with a susceptibility that dominates over other broken symmetry states mainly due to subdominant terms. In strong coupling, we construct a ground state wavefunction for , , , but , which shows explicitly the broken layer inversion. Since in this model, the same symmetry appears to be broken in the limit of both weak and strong coupling, it is reasonable to assume that such a broken symmetry state appears at any .

A similar analysis is presented in the spin- case with short range interactions. For the repulsive Hubbard model, we find that the most dominant weak coupling instability is towards an anti-ferromagnetic state. Since the same ordering tendency happens in the strong coupling, it is reasonable to assume that in this model, the Neel ordering appears at any .

This paper is organized as follows: in Section II we write down the (non-interacting) bilayer Hamiltonian first in the tight-binding approximation and then within perturbation theory. In Section III we construct the low energy effective theory at the neutrality point by fine-tuning the trigonal warping terms to zero. The rest of that section deals with identifying microscopic-symmetry-allowed 4-fermion contact interaction terms using the method of Herbut, Juricic and RoyHerbut et al. (2009) used for the same purpose in single-layer graphene. Before the reduction due to Fierz identities, there are such couplings which further reduce to 4 once Fierz identities are taken into account. The weak coupling RG is presented in Section IV, along with the flow diagram in the space of coupling constant ratios and the analysis of the susceptibility growth. The model with weak and strong coupling limits is studied in Section V. In Section VI, the spin- case is revisited. Symmetry is used to construct an eighteen-dimensional Fierz vector along with the Fierz matrix to show that there are independent couplings in this case. Their RG equations are determined and while more general, they are shown to reduce to the ones studied in Ref.Vafek and Yang (2010) under conditions outlined therein. In Section VII we study the Hubbard model in weak and strong coupling. Section VIII is devoted to conclusions. Details of the derivation are presented in the Appendices.

Ii Bilayer Hamiltonian

In this section we will define the non-interacting model by using two different approximation methods. First, the well known tight binding approximation Castro Neto et al. (2009) will be used and then the -method, or equivalently the method of invariantsBir and Pikus (1974); McCann and Fal’ko (2006); Aleiner et al. (2007); Lemonik et al. (2010). Both methods lead to the same form of the low energy Hamiltonian and it is ultimately a question of convenience which one should be adopted.

Figure 1: (i) Schematic representation of the A-B stacked bilayer. The low energy wavefunction near is also sketched, with . The primitive lattice vectors are and . The area of the unit cell is . (ii) Schematic representation of the (reciprocal) -space.

ii.1 Tight-binding approximation

The non-interacting Hamiltonian in the tight-binding approximation can be written as




In the case of bilayer graphene, the values of the hopping integrals , , were extracted experimentally in Ref.Zhang et al. (2008). If we define the Fourier transform of a Fermi field as where or , or , and is the number of unit cells. Next, we let to write the non-interacting Hamiltonian (1) as


In the above, the wavevector dependent function where the sum runs over , and . Near , . Near , . The low energy spectrum of this (well-known) HamiltonianMcCann and Fal’ko (2006); Nilsson et al. (2008), which is easily diagonalized, will be discussed in the next section.

ii.2 approach

Instead of resorting to the tight-binding approximation, we can also arrive at the low energy Hamiltonian by analyzing the symmetry of the bilayer potential alone. This is a standard technique when dealing with semiconductorsBir and Pikus (1974) and one which has also been applied to grapheneAleiner et al. (2007). For the sake of self-inclusiveness, we present this method as well to show that one arrives at the same general form of the Hamiltonian as in the tight-binding approximation, although in practice the coefficients of various symmetry-allowed terms must be determined from experiment. We start with the Schrodinger equation for a particle moving in potential due to the atoms in layers and separated by




The low energy field theory is written in terms of the eight-component Fermi fields (two layers, and , two valleys, and , and two sublattices and as sketched in Fig.(1)):


The rapidly-varying Bloch functions at and at are related by complex conjugation, , irrespective of the layer or sublattice index. Moreover, the Bloch functions and transform irreducibly under point group operations of the lattice (see Fig.1). For the sake of concreteness, within the nearly free electron approximation for electron wavefunctions confined to layers and respectively we have


where .


i.e. the interlayer hopping arises from the mixing of the sublattices and . The matrix elements of the in-plane momentum operator are also dictated by symmetry to be


Defining , gives us the effective Hamiltonian near to read


This is equivalent to what we found in the tight-binding approximation.

The spectra of the and the tight-binding Hamiltonians are well known and have been discussed extensively in the literature (See e.g. McCann and Fal’ko (2006); Zhang et al. (2008); Castro Neto et al. (2009)). In the vicinity of each -point, there are four Dirac points: one isotropic at and three anisotropic ones arranged in accordance with 3-fold lattice symmetry around the isotropic one. When we neglect trigonal warping terms, by setting , or set the higher order hopping terms , the four Dirac points merge into a parabolic degeneracy.

Iii Low energy effective theory

In the weak coupling limit, the kinetic energy dictates which modes are important to determine the behavior of the system at low energies. Clearly, at we have two degenerate levels and two levels at . Since we wish to work with a theory for the low energy modes only, we need to project out the bands which originate from the two ”split-off” bands. We can do so in several equivalent ways. The method used here implements the path integral formalism, where we integrate out the Fermi fields associated with and modes (sites), and arrive at an effective action with an effective ”Hamiltonian” for the low energy modes. In addition to the wave vector dependence, this ”Hamiltonian” is frequency dependent as well. Near the -point, the effective quadratic action after integrating out the -modes is

Since the integral is Gaussian, we can easily perform it and find that up to an additive constant





Within the theory, the parameters and should be determined from experiment. To make contact with the notation in literature, Ref.Nilsson et al. (2008) have and (see their Eqs. 6 and 15).

If we are interested in the modes near the Fermi level of an unbiased bilayer, we can simply set in the effective action (19-22). As will be obvious from the discussion in the next section, terms arising from the corrections are perturbatively irrelevant near the Gaussian fixed point in the sense discussed in a different context in Ref.Shankar (1994).

In what follows we will also set to fine-tune the system to quadratic degeneracy. Such a situation arises if in the tight-binding formulation we consider only the nearest neighbor hopping integrals, and . Otherwise, as mentioned in the introduction, the ultimate low energy dispersion involves four (one isotropic and three anisotropic) Dirac conesMcCann and Fal’ko (2006); Nilsson et al. (2008); Lemonik et al. (2010). While such fine-tuning appears artificial, it is an example of the maximRaghu et al. (2010) that there are no generic weak coupling particle-hole instabilities. Rather, fine-tuning, in the form of nesting for example, is necessary to bring the strong coupling physics down to weak coupling. If we are interested in accessing the symmetry-breaking phases in the particle-hole channel, as we are in this case, then fine tuning is a small price to pay for this access made available within perturbative RGRaghu et al. (2010).

Putting back the point, the low energy degrees of freedom can now be expressed in terms of a four component Fermi field

i.e. the electronic degrees of freedom are expanded as


The non-interacting low energy (imaginary time ) Lagrangian, which includes both and valleys, and which will serve as our (gaussian) fixed point of departure, can therefore be written as


where we defined the vector function and the matrices as


The effective mass parameter entering the above equations is (In the tight binding approximation ). The four component Fermi objects appearing in Eq. (24) were defined as the envelope Fermi fields in Eq.(23). In the above, the first Pauli matrix acts in the valley -space and the second in the layer -space. To make contact with the literature we also use Dirac -matrices which we represent as


The action is invariant under the scale transformation


This means that the ”dynamical critical exponent” for the gaussian theory, which will be our point of departure when analyzing weak coupling instabilities.

iii.1 Short range interactions

From the above discussion of the gaussian fixed point, it is evident that the short range interactions, when projected onto our low energy modes, will contain among other (perturbatively irrelevant) terms, contact four-fermion terms which are marginal by power counting. The rest of this section deals with identifying such symmetry-allowed interaction terms. The method used here follows almost verbatim the method used by Herbut, Juricic and RoyHerbut et al. (2009) in their analysis of the short range interactions in single layer graphene. In addition to the lattice symmetries used in Ref.Herbut et al. (2009), we also include the three-fold rotational symmetryAleiner et al. (2007), which reduces the number of independent four-fermion couplings to .

We can therefore start by writing the general Lagrangian


where was introduced in Eq.(24) and


where, at this point, the sum over includes all sixteen independent four-by four matrices (generators of ) and so does the sum over . Naively, we have couplings to consider. Just as in the case of the single-layer grapheneHerbut et al. (2009), this number will be dramatically reduced first by using the discrete symmetries of the lattice and second by using Fierz identities.

The key role in this reduction is played by the behavior of the Bloch functions under symmetry operations, which dictates the transformation properties of the four component, slowly varying, envelope Fermi fields Bir and Pikus (1974); Aleiner et al. (2007). The dimer centered rotation by , mirror reflection about the yz-plane and about the xz-axis followed by xy-plane respectively give


The time reversal symmetry and translational symmetry give


In the above, where and . And since , , where . The lattice translational symmetry therefore corresponds to the discrete analog of the chiral generated by .

iii.1.1 Symmetry reduction

Following Herbut et al. (2009), we split the sixteen linearly independent four-by-four matrices and into four sets


The matrices which belong to the set are even under both reflection operations (38) and (39). The matrices in the set are odd under -reflections (38) and even under ””-reflections (39). The matrices in the set are even under -reflections (38) and odd under ””-reflections (39). And finally, matrices belonging to the set are odd under both (38) and (39). This means that only quartic terms combining matrices from the same set are allowed by symmetry. Each such set contains such terms and that leaves couplings.

Eight matrices , , and are left invariant under the spatial translation operation (41). These give rise to couplings, eight direct ( and ), as well as four mixed . In addition, there are four sets of pairs which transform as vectors under (41): , , and . These give rise to additional couplings. Schematically, four of them are and two mixed ones are and . Altogether, after inclusion of the translation symmetry, we are left with couplings.

The unitary part of the time reversal operations , Eq. (40), happens to correspond to the mirror reflection about , (Eq.38), which has already been taken into account. However, complex conjugation, further restricts the number of couplings. Specifically, mixed terms with one purely real and one purely imaginary matrix cannot appear, therefore . This leaves couplingsHerbut et al. (2009).

The lattice symmetries considered by Herbut et.alHerbut et al. (2009) did not contain site- or plaquette- centered rotationAleiner et al. (2007) by . As stated in Eq.(37), this symmetry is generated by . Including this symmetry requires that the cross-terms . Moreover, it requires that , , and .

This leaves us with the following 9 terms


The 9 terms can be further reduced to 4 independent ones by using Fierz identities.

iii.1.2 Fierz identities

Figure 2: Diagrams appearing at 1-loop RG. Each vertex can be represented by either a 44 matrix (spinless) or 88 matrix (spin-) case.

We set to continue with the notation of Ref.Herbut et al. (2009). We use the method employed therein to write down Fierz identitiesHerbut et al. (2009); Itzykson and Zuber (2005) which, due to the Grassman nature of the Fermi fields, relate various seemingly unrelated couplings.

The starting point is the algebraic identity (see Eq.(A4) of Ref.Herbut et al. (2009))


which leads to


The minus sign comes from and being anti-commuting (four component) Grassman fields. For contact terms and the above equation (III.1.2) constitutes a set of linear relations between different terms of our symmetry reduced interaction Lagrangian (46).

If we arrange the quartic terms into a vector


then the Fierz identities lead to the linear constraint


A straightforward, though somewhat laborious, application of (III.1.2) leads to the explicit form of the Fierz matrix in the case of spinless fermions


The matrix has four zero eigenvalues and as a result there are four independent couplingsHerbut et al. (2009).

In order to make a connection with the previous workVafek and Yang (2010), we choose to eliminate




in favor of the remaining four terms. These equations will be used in deriving our RG equations, since elimination of fast modes will generate terms such as, for example, . The above equations show that such a term does not correspond to a new coupling in a renormalized action, but rather is a linear combination of terms already present.

Finally, we arrive at our interaction Lagrangian


Above is the most general four-fermion contact interaction Lagrangian for spinless fermions allowed by the symmetry of the A-B stacked honeycomb bilayer. In the next section, we study the weak coupling RG flow of the four couplings and . The first three couplings appeared in our previous workVafek and Yang (2010) where we called them , , and . The fourth coupling, , did not appear there since the starting point assumed only finite and, as we will see later, is not generated if its starting value is zero.

Figure 3: Flow diagram in the coupling constant ratio space assuming that (generic behavior). There are two sinks given by Eqs.(73) and (75) in the text and two mixed fixed ratios (74) and (76). For the flow is reversed, and generically, the divergent coupling constant ratios simply mean that has shrunk and crossed . After this point becomes negative and the directionality shown here is restored.

Iv Renormalization group analysis

Clearly, ’s are marginal by power-counting and the question is how they flow. The RG procedure employed here follows Ref.Shankar (1994) and consists of integrating out the fermionic modes in a thin shell between the initial cutoff and , while the integral over extends from to . Since