Nematicity with a twist: Rotational symmetry breaking in a moiré superlattice

See allHide authors and affiliations

Science Advances  05 Aug 2020:
Vol. 6, no. 32, eaba8834
DOI: 10.1126/sciadv.aba8834


Motivated by recent reports of nematic order in twisted bilayer graphene (TBG), we investigate the impact of the triangular moiré superlattice degrees of freedom on nematicity. In TBG, the nematic order parameter is not Ising like, as in tetragonal crystals, but has a three-state Potts character related to the threefold rotational symmetry (C3z) of the moiré superlattice. We find that, even in the presence of static strain that explicitly breaks the C3z symmetry, the system can still undergo a nematic-flop phase transition that spontaneously breaks in-plane twofold rotations. Moreover, elastic fluctuations, manifested as acoustic phonons, mediate a nemato-orbital coupling that ties the nematic director orientation to certain soft directions in momentum space, rendering the Potts-nematic transition mean field and first order. In contrast to the case of rigid crystals, the Fermi surface hot spots associated with these soft directions are maximally coupled to low-energy nematic fluctuations in the moiré superlattice case.


Recent experiments in twisted bilayer graphene (TBG) have unveiled a rich phase diagram displaying phenomena often observed in strongly correlated systems, such as superconductivity, metal-insulator transition, nematicity, ferromagnetism, and the anomalous quantum Hall effect (17). In contrast to correlated materials, however, in TBG, these phases can be precisely tuned with gating and twist angles, rather than doping and pressure. This offers a compelling venue to elucidate correlation-driven electronic orders (821), which avoids the typical complications related to bulk compounds.

In correlated materials, nematic order—a correlation-driven lowering of the point group symmetry of a crystal—is often associated with unconventional superconductivity (22, 23). In TBG, scanning tunneling microscopy (STM) (2426) and transport measurements (27) reported evidence that the threefold rotational symmetry of the moiré superlattice, denoted by C3z, is broken in different regions of the TBG phase diagram, including close to the superconducting dome. Moreover, spontaneous C3z symmetry breaking has been invoked to explain the observed Landau level degeneracy at charge neutrality (28, 29). While a number of theoretical proposals have explored the possibility of such an electronic nematic phase (3034), experimentally, it remains a difficult task to distinguish spontaneous nematic order from an explicit broken symmetry caused by strain, whose presence is ubiquitous in TBG (26, 3537).

An important question is whether nematicity in TBG is a straightforward generalization of nematicity in bulk correlated systems or whether new effects arise from the unique properties of twisted systems. A defining property of nematic order is its strong coupling to the lattice. Because the moiré superlattice in TBG is triangular (see Fig. 1A), nematic order is not described by an Ising order parameter, as it is the case in tetragonal pnictides and cuprates (22), but rather by a three-state Potts model order parameter, associated with the three C3z-related lattice bonds (see Fig. 1B). Importantly, the moiré superlattice is an emergent lattice, whose elastic properties are fundamentally distinct than the rigid crystals of correlated materials (38, 39).

Fig. 1 Nematic order in a moiré superlattice.

(A) The triangular moiré superlattice of TBG (blue bonds), formed by the AA stacking regions (black dots). (B) In the presence of nematic order, one of the superlattice bonds becomes different (red bond), while the other two remain equivalent. (C) The symmetries of the moiré superlattice involve twofold rotations with respect to orthogonal in-plane axes (C2x and C2y) and sixfold rotations with respect to the z axis (C6z). (D) Nematic director n̂=(cosθ,sinθ): black (red) corresponds to γ < 0 (γ > 0) in the action in Eq. 1. Note that n̂ and n̂ (dashed arrows) are identified. (E) Bond order pattern and lattice distortion pattern associated with each nematic director.

Here, we show that nematicity in twisted systems can be fundamentally different than in bulk correlated systems. While static strain completely smears the nematic transition in tetragonal lattices, it allows the moiré superlattice to still undergo an Ising-like nematic-flop transition, in which in-plane twofold rotational symmetries are spontaneously broken. Finite-momentum strain fluctuations, manifested as acoustic phonons, mediate a nonanalytic nemato-orbital coupling in the moiré superlattice. The latter makes certain directions in momentum space, which are tied to the nematic director’s orientation, softer than others across the nematic transition. This is expected to render the three-state Potts-nematic transition mean field and first order and also constrains the electronic states that can exchange low-energy nematic fluctuations to a discrete set of Fermi surface hot spots. Because the moiré superlattice is not a rigid crystal (38, 39), the nematic form factor is maximum at these hot spots. This contrasts with rigid lattices, where the form factor vanishes at the hot spots, effectively decoupling the electronic system from low-energy nematic fluctuations. Thus, the maximum coupling between hot spots and nematic fluctuations makes moiré superlattices promising systems to elucidate the impact of nematicity on electronic properties, such as transport and superconductivity.


Electronic nematic order is described by a traceless symmetric tensor, which, in two dimensions, has two independent components Φ1 and Φ2 corresponding to the charge quadrupole moments with dx2−y2 and dxy symmetries, respectively. In systems with tetragonal symmetry, these two d-waves have distinct symmetries and must thus be treated as two independent Ising order parameters. This is markedly different in hexagonal systems, such as TBG with point group D6: The two nematic components belong to a single irreducible representation of D6 and transform as partners under its symmetries, defining a two-component order parameter Φ = (Φ1, Φ2). It is natural to parametrize it as Φ = Φ(cos2θ, sin2θ), where the angle θ can be identified with the orientation of the nematic director n̂=(cosθ,sinθ) (see Fig. 1D); note that Φ(θ) = Φ(θ + π), as expected. Note that all of this remains true if one considers a TBG model with D3 point group symmetry instead.

Although this parametrization might suggest that Φ is an XY order parameter, the lattice symmetries of TBG introduce crystal anisotropy effects that pin the nematic director to a discrete set of high-symmetry directions. The Landau-type action Snem[Φ] is [see also (30, 4042)]Snem[Φ]=S0[Φ]+γ6x(Φ+3+Φ3)(1)where x = (r, τ) denotes spatial coordinate r and imaginary time τ, and Φ± ≡ Φ1 ± iΦ2. The first term, S0[Φ]=12rΦΦ2+14uΦΦ4, is a standard Φ4 action with U(1) symmetry. The cubic term reflects the crystalline anisotropy of the hexagonal lattice and is expressed as 13γΦ3cos6θ, which is minimized by θ = 2nπ/6 for γ < 0 and θ = (2n + 1)π/6 for γ > 0. These solutions correspond to sets of threefold degenerate nematic directors, as shown in Fig. 1D [recall that angles differing by π (dashed arrows) must be identified], and manifested as bond orders in real space (Fig. 1E). Equation 1 is the continuum version of the three-state Potts model, where the Z3 symmetry is identified with the threefold rotation about the z axis, C3z. Below the nematic transition temperature Tnem, the sixfold rotation symmetry C6z about the z axis is broken down to a twofold symmetry C2z, while the twofold rotation symmetries C2x and C2y with respect to the in-plane x and y axes (or their symmetry-related equivalents) are preserved (see Fig. 1C). Despite the presence of a cubic term in (Eq. 1), the three-state Potts transition is continuous in two dimensions.

Static strain

According to the Landau theory of phase transitions, once threefold rotation symmetry C3z is broken because of condensation of the nematic order parameter Φ, all order parameters that also break C3z symmetry and linearly couple to Φ will become nonzero—irrespective of their microscopic origin (be it electronic or structural). For the triangular moiré superlattice, nematic order implies a unit cell–preserving redistribution of charge between nearest-neighbor bonds, making one bond inequivalent from the remaining two. This is shown in Fig. 1E and is referred to as bond order. We emphasize that, due to the three-state Potts model nature of this order, two of the three bonds remain equivalent. Depending on the sign of the nematoelastic coupling, the distinct bond, shown in red in Fig. 1E, is either stronger or weaker compared to the other two bonds, as we discuss below.

As a result, each one of the three allowed electronic nematic directors (black or red arrows in Fig. 1D) corresponds to selecting one of the three bonds in real space. In the nematic state, the “red” bonds form unidirectional stripes along the moiré superlattice, as shown in Fig. 1B. Such a pattern was recently observed in STM experiments in TBG (26). Despite the appearance of stripes, this type of charge redistribution pattern is described by a zero-momentum order parameter, since there is no additional translational symmetry breaking.

To quantitatively account for the coupling between the electronic and lattice degrees of freedom, we include the latter via the strain tensor εij12(iuj+jui) and the rotation tensor ωij12(iujjui), which are defined in terms of the moiré superlattice displacement vector u. In the present case, the electronic nematic order parameter Φ may refer to any fermionic bilinear that breaks C3z symmetry, which could be in the charge, valley, or spin channel. On the other hand, εij and ωij refer only to the lattice geometric degrees of freedom. The elasto-nematic action is given by Selnem[Φ,ε^,ω^]=Sel[ε^,ω^]+S[Φ,ε^], where Sel[ε^,ω^] is the elastic free energy andS[Φ,ε^]=λx[(εxxεyy)Φ1+2εxyΦ2](2)describes the nematoelastic coupling with coupling constant λ.

We consider first the effect of static strain applied in TBG. For compressive (tensile) uniaxial strain ε < 0 (ε > 0) applied parallel to an arbitrary direction d^, the action above becomes S′ = − λ∫xε Φ cos(2α − 2θ), where cosα=d^·x^. At high temperatures TTnem, where Tnem is the transition temperature of the unstrained system, we can approximate Snem12xχnem1Φ2. Thus, strain not only triggers a finite nematic order parameter Φ ∝ χnem∣ε∣ but also pins the nematic director parallel or perpendicular to the strain direction; i.e., θ0 = α or θ0 = α + π/2, depending on whether λε > 0 or λε < 0, respectively. To understand what happens as temperature is lowered, we consider TTnem and set ∣Φ∣ = Φ0 as approximately constant. Expanding around the high-temperature director, θ = θ0 + δθ, givesSnem+S=x[aθ0(δθ)+bθ0(δθ)2](3)with coefficients aθ0=2γΦ03sin6θ0 and bθ0=2Φ0(λε3γΦ02cos6θ0). Let us focus first on the case when strain is applied along a high-symmetry direction (α = nπ/6). It follows that aθ0 = 0 and that the system remains symmetric under twofold rotations with respect to in-plane axes, C2x and C2y. These symmetries can nevertheless be broken spontaneously, which is signaled by bθ0 < 0 in Eq. 3. This can only happen if θ0 coincides with the maxima, but not the minima, of the cubic term in Eq. 1—in other words, if the strain term S′ is minimized by a director that is maximally penalized by the cubic term of Snem. In this case, once Φ0 reaches the critical value Φ¯0=λε3γ, usually at a temperature Tnemflop>Tnem, the minimum changes from θ0 to θ0±θ¯0, with θ¯0=12arccos(121+3Φ¯02Φ02), resulting in an Ising-like transition that spontaneously breaks the C2x and C2y symmetries. Because of its resemblance to a spin-flop transition, we call the reorientation of the nematic director under an external field a nematic-flop transition. Thus, as illustrated in the schematic phase diagram of Fig. 2A, a nematic-related phase transition can still occur in a strained triangular lattice (43), in contrast to the case of a strained tetragonal lattice, where only a crossover exists.

Fig. 2 Nematic transition in the presence of static strain.

(A) Schematic temperature versus strain phase diagram with strain applied along the y axis (α = π/2) and λ < 0, γ > 0. For compressive strain (ε < 0), because the director is fixed at θ0 = π/2, which is a minimum of the cubic term, no phase transition occurs, and only a crossover temperature Tnem* survives. For tensile strain (ε > 0), the director is at θ0 = 0, which is a maximum of the cubic term, for T>Tnemflop, and at ±θ¯0 for T<Tnemflop. Thus, Tnemflop marks an Ising-like nematic-flop transition in which the twofold rotational symmetries C2x and C2y are spontaneously broken. The sixfold rotation C6z is explicitly broken to C2z everywhere for ε ≠ 0. (B) Evolution of the bond order patterns as temperature is lowered in the cases of compressive and tensile strain. Only in the latter case, the three bonds become inequivalent at low temperatures, signaling the breaking of C2x and C2y.

From the perspective of bond order, the breaking of the twofold C2x/C2y symmetries (in addition to C3z), makes all three bonds connecting nearest-neighbor sites inequivalent. Recall that, in the nematic phase, where only the C3z symmetry is broken, two of the three bonds remain equivalent. Therefore, referring to the phase diagram of Fig. 2A, the sequence of bond order patterns observed as temperature is lowered would be the one illustrated in Fig. 2B. Without applied strain, the three bonds are equivalent. Application of compressive strain along one of the bonds would make that bond distinct and the other two equivalent. This is an explicit symmetry breaking induced by the strain. Upon lowering the temperature, no additional symmetry breaking occurs. However, for tensile strain, once the temperature becomes smaller than Tnemflop, one of the two remaining equivalent bonds would become different. This is a spontaneous Ising-like symmetry breaking, since one among two bonds is chosen.

Since STM has been able to identify bond order patterns in the nematic phase (26), it provides a potential tool to probe such a nematic-flop transition. An appealing experiment, capable of distinguishing between spontaneous nematicity and strain-induced anisotropies, would consist of comparing the types of bond order that appear in the presence of tensile and compressive strain applied along one of the crystalline directions. If in both cases one sees unidirectional bond order patterns, then this would imply that the system is in the nematically disordered phase and that strain is causing the anisotropies. However, if in one case the type of bond order evolves as a function of temperature according to what is shown in the lower half of Fig. 2B, then this would be evidence for long-range nematic order.

Let us now consider the case when strain is not applied along one of the high-symmetry directions, i.e., α ≠ nπ/6; then, aθ0 ≠ 0 in Eq. 3. As a result, the applied strain breaks not only C3z but also C2x/C2y symmetry. Consequently, no true phase transition can happen. Nevertheless, near Tnemflop, the nematic director θ will display strong changes as temperature is lowered, which is manifested as a sharp temperature dependence of the bond order pattern. If the system does not have a nematic ground state, however, then θ will change only weakly as a function of temperature.

Fluctuating strain

Apart from static strain, finite-momentum elastic fluctuations are expected to strongly affect the nematic transition (4446). For a system with D6 symmetry, diagonalization of the harmonic elastic action Sel[ε^,ω^] leads to two acoustic phonon modes, a transverse (T) and a longitudinal (L) one with sound velocities vL, TSel=12μ=L,Tqũq,μ(ωn2+vμ2q2)ũq,μ(4)where q = (q, ωn), with ωn the (bosonic) Matsubara frequency. The displacement field u=μũμe^μ has been decomposed into its longitudinal and transverse components ũμ with e^L=(cosζq,sinζq) and e^T=(sinζq,cosζq) and ζq = arctan (qy/qx). According to (38, 39), for the dominant acoustic phonons that act on the moiré superlattice scale, u corresponds to the relative displacement of the two graphene sheets. These and other phonon modes have been proposed to be linked to superconductivity in TBG (18, 4749). Integrating out the acoustic phonons leads to an additional contribution to the nematic action, δSnem=12ijq Φi,qΠ^ij(q)Φj,q. In the static limit, ωn = 0, we findΠ^=λ2vT2[I^ηP^],P^=(cos22ζq12sin4ζq12sin4ζqsin22ζq)(5)where I^ is the identity matrix and η1vT2/vL2. The first term of Π^ gives an overall enhancement of Tnem. The second term couples the two nematic components Φ1 and Φ2 in a way that depends on the direction but not on the magnitude of q. Such a nonanalytic term typically appears when order parameters couple linearly to an elastic mode (50) and was previously studied for Ising nematic order in tetragonal lattices (44, 45). Here, it is manifested as a nemato-orbital couplingSnem(eff)[Φ]=S0[Φ]+γ6x(Φ+3+Φ3)+λ2vT2[xΦ2+ηq(Φ·D^)2](6)where D^=(cos2ζq,sin2ζq)=(q^x2q^y2,2q^xq^y) is the momentum space (i.e., orbital) quadrupolar form factor. Recasting the nemato-orbital coupling term as Φ2cos2(2θ − 2ζq), we see that it makes only certain directions of momentum space to become soft; i.e., the static nematic susceptibility χnem(q0,q^) is largest near Tnem only along special directions q^ [see also (45)]. While the cubic term in Eq. 6 forces the director n^ to point along one of three directions (see Fig. 1D), the nemato-orbital coupling makes only two momentum space directions ζq soft, namely, the ones that make a relative angle of 0 and π2 (for η < 0) or ±π4 (for η > 0) with respect to n^ (see Fig. 3). The reduction of the soft-direction phase space from continuous to discrete is known to effectively enhance the dimensionality of the Φ4 action S0[Φ] from d to d + 1. Thus, one expects that the nematic transition in the moiré superlattice will be the same as a three-dimensional three-state Potts model transition, which is mean field and first order.

Fig. 3 Phonon-mediated nemato-orbital coupling.

Momentum directional dependence of the nematic susceptibility χnem(q0,q^) caused by the nemato-orbital coupling, with q^=(cosζq,sinζq). Light blue (dark blue) denotes softer (harder) directions, corresponding to higher (lower) susceptibility. While for a rigid crystal, vT < vL, the soft direction is rotated by ±π/4 with respect to the nematic director n (red arrow), for TBG, vT > vL, the rotation is 0, π/2. Note that the director can point in any of the directions θ of Fig. 1B.

Electronic degrees of freedom

If the first-order character of the nematic transition discussed above is weak, then nematic fluctuations are still expected to affect the electronic degrees of freedom. For a single-band system with fermionic operator ck, the electronic nematic coupling is Selec=k,qg(k)Φqckq/2ck+q/2, with form factor g(k) = g0 cos(2θ − 2θk), where g0 is a constant and θk = arctan (ky/kx). The electronic states that exchange low-energy nematic fluctuations are at the Fermi surface and separated by the small momentum q of the nematic mode. Since the nematic fluctuations are the softest (albeit nondiverging) along the special directions qˆ(0)=(cosζq(0),sinζq(0)) discussed above, the relevant pairs of fermions are located around the “hot spots” khs where the Fermi surface’s tangent is parallel to q^(0), i.e., q^(0)·ξkhs=0. The issue is how strong these fermions are coupled to the nematic fluctuations, i.e., what the magnitude of g(khs) is. For a circular Fermi surface, the hot spots are located at θkhs=ζq(0)+π2, and thus g(khs)=g0cos(2θ2ζq(0)). As we saw above, if η > 0, then the soft directions are ζq(0)=θ±π/4, yielding g(khs) = 0. Thus, in this case, the hot spots effectively decouple from the softest nematic fluctuations, similarly to what was obtained for an Ising nematic tetragonal lattice (45). On the other hand, if η < 0, then the soft directions are ζq(0)=θ,θ±π/2, implying that ∣g(khs)∣ = ∣g0∣; i.e., the hot spots are maximally coupled to the soft nematic fluctuations. For a generic noncircular Fermi surface respecting D6 symmetry g(khs) remains maximum for η < 0 but is expected to be nonzero albeit small for η > 0.

The sign of η1vT2/vL2 is determined by the elastic action Sel[ε^,ω^]. For a rigid crystal, Sel[ε^,ω^]=12x[(τu)2+Cijklεijεkl] depends only on the strain ε^, since global rotations do not cost energy. In a triangular lattice, there are only two independent elastic constants, C11Cxxxx and C12Cxxyy, yielding vL2=C11 and vT2=(C11C12)/2. Lattice stability requires C11 > ∣C12∣, which makes η > 0, implying that g(khs) is small. However, the moiré superlattice is not a rigid crystalline structure for small twist angles, as lattice relaxation leads to sharp domain walls separating the regions with AB and BA stacking. Because of this, arbitrary rotations of the moiré superlattice cost energy, and the elastic free energy acquires an extra term, δSel[ω^]=12xK ωxy2 (39). This term contributes only to the transverse velocity and when K > 2(C11 + C12), vT becomes larger than vL (i.e., η < 0), implying that ∣g(khs)∣ is maximum. Recent calculations of the acoustic phonon spectrum of TBG found that this condition is satisfied for small twist angles (38, 39), making TBG a rather unique system in which the Fermi surface hot spots are maximally coupled to the nematic fluctuations.

To apply these results specifically to TBG, we use the six-band model in (51), generalizing it to the case where different types of nematic order are present (for details, see Materials and Methods). As shown in Fig. 4A, there are two Fermi surfaces associated with the two valley degrees of freedom and thus related by a C2z rotation. Because the two pairs of hot spots for a given nematic director θ are related by π/2 rotations, they correspond to different valley symmetries. Setting θ = 0 for concreteness, we find that the pair of hot spots located at θkhs = 0, π is associated with intravalley nematicity (Fig. 4B), whereas the pair located at θkhs = ± π/2 is associated with intervalley nematicity (Fig. 4C).

Fig. 4 Fermi surface distortion and nematic hot spots.

(A) Fermi surface of the six-band model in (51); red and black correspond to the two valleys. The two pairs of hot spots are marked by open and full symbols. (B) Distortion of the Fermi surface in the presence of intravalley Potts-nematic order, with nematic director n along the x axis. (C) Same as (B), but for intervalley nematic order. In (B) and (C), the undistorted Fermi surface is shown by the dashed lines.


The manifestations and consequences of electronic nematicity in TBG that we have analyzed and discussed do not depend on the microscopic nature of nematic order; i.e., they do not depend on whether nematicity originates from charge, valley, or spin degrees of freedom. Nevertheless, it is important to discuss the possible microscopic mechanisms for Potts-nematic order in TBG. In weak-coupling approaches, a Pomeranchuk instability breaking the C3z rotational symmetry of the Fermi surface can be favored by van Hove singularities (52, 53). In strong-coupling approaches, where charge degrees of freedom are quenched, a widely used effective Hamiltonian is described in terms of an SU(4) “super spin” associated with spin and orbital variables (21, 30, 5458). These orbital variables are associated with eigenstates of the C3z symmetry operation and, in certain formulations, are nearly fully valley polarized (21). In this SU(4) formulation, nematicity is then described by an ordering of the orbital variables, i.e., ordering in the SU(2) orbital sector, which breaks spatial rotational symmetry. Whether the ground state of the effective SU(4) Hamiltonian is a nematic phase is an interesting open question. A third possible mechanism is a nematic phase that is a vestigial order of a primary electronic ordered state that breaks C3z and some additional symmetry (23), such as p + p-wave/d + d-wave superconductivity (30, 34, 40) or charge/spin stripe density waves (41).

Thus, while there are several other possible ground states in TBG, a number of them are fully compatible with electronic nematic order, as they break C3z symmetry (among other symmetries). The study of the detailed microscopic mechanisms responsible for the various forms of correlated behavior in TBG is currently a very active and rapidly progressing area of research. Although a detailed microscopic description is not the focus of this work, which instead focuses on the manifestations and consequences of nematic order, motivated by experimental evidence, we can nonetheless provide a rough estimate of the nematic transition temperature. We start by noting that the energy scales of the narrow band bandwidth and of the screened Coulomb repulsion are both of the order of 10 meV. Therefore, since weak-coupling calculations often find instabilities at temperatures of a few percent of the bandwidth, whereas strong-coupling calculations usually give instabilities at temperatures of a few percent of the Coulomb repulsion, a reasonable range for the nematic transition temperature would be between 1 and 10 K.

To conclude, we showed that the Potts-like character of the nematic order parameter in triangular moiré superlattices leads to unique nematic behaviors seen neither in tetragonal systems nor in rigid triangular crystals. Notably, a nematic-flop phase transition that spontaneously breaks the in-plane twofold rotational symmetries can still take place even when C3z symmetry-breaking strain is applied. This can only happen, however, for one sign of strain (i.e., either compressive or tensile) and if the ground state of the unstrained system is nematic. Thus, because the breaking of the twofold in-plane rotational symmetries is manifested as a peculiar type of bond order in the triangular lattice, this makes it possible to unambiguously detect long-range nematic order in TBG even in the presence of strain.

Moreover, we showed that the emergence of a nemato-orbital coupling mediated by acoustic phonons affects not only the character of the Potts-nematic transition, which becomes mean field and first order, but also the impact of the low-energy nematic fluctuations on the electronic properties, which is maximized because of the nonrigid nature of the moiré superlattice. Such a maximal coupling between nematic fluctuations and electrons, which is, in turn, mediated by phonons, could, in principle, lead to unusual temperature dependencies of thermodynamic and transport quantities in an intermediate temperature regime, provided that the resulting first-order character of the Potts-nematic transition is weak. In this context, it is interesting to note the recent observation of linear-in-T resistivity in TBG, whose origin remains under debate (59, 60). It would be desirable to investigate whether this phonon-assisted electronic nematic coupling can promote nontrivial properties of response functions. Similarly, it will be important to assess the importance of the coupling between phonons and nematic fluctuations in phonon-based pairing mechanisms that have been proposed to explain superconductivity in TBG (18, 4749).


The results discussed in the section “Electronic degrees of freedom” of the main text for the coupling between nematic fluctuations and the low-energy electronic states are quite general, as they rely on the fact that the system has sixfold rotational symmetry. To make the connection with TBG more transparent, in the main text, we also considered a specific six-band model for TBG in the presence of nematic order. Here, we provide the details of implementing nematic order in this tight-binding model, introduced in (51). We will generally follow the notation in (51), with a few exceptions.

Six-band tight-binding model

The six-band tight-binding model in (51) is defined by (pz, p+, p) orbitals on the sites of a triangular lattice and s orbitals on the sites of a kagome lattice. The purpose of the six-band model is to reproduce the low-energy flat bands of TBG in such a way that all symmetries manifestly present in the continuum description are respected (and are implemented locally).

Note that the six-band model describes the two low-energy flat bands originating from a single valley of the two graphene sheets forming the bilayer system. This implies that a model that includes the full set of flat low-energy bands (apart from spin) must have 12 bands: two copies of the 6-band model related by the symmetries that exchange valleys. Here, we sketch how such a model is constructed and show how it can be supplemented with the appropriate symmetry-breaking terms to account for nematic order.

We begin by recalling the definition of the six-band model. The orbital degrees of freedom can be represented by a fermion operator ψk given byψk=(pkz,pk+,pk,ak,bk,ck)(7)

Note that this notation is a slight departure from the one in (51). In terms of these degrees of freedom, the (six-band) Hamiltonian for a single valley is given byHk=(Hpz+μpzCp±pz0Cp±pzHp±+μp±Cκp±0Cκp±Hκ+μκ)(8)

Here, Hpz and Hp± are the subblock Hamiltonians in the pz and p± subspaces, and Hκ describes the coupling between the kagome lattice sites. The subblocks CXY describe the couplings between the X and Y sectors (where X, Y = pz, p±, κ). The form of all these subblocks is given in (51), including the parameter set that we use here [see table VI in (51).

To promote the six-band model to a full twelve-band model, we take two copies and introduce a valley degree of freedom asΨk=(ψk+,ψk)(9)where ± labels the K and K′ valleys. The full Hamiltonian Hk is then given byHk=(HkUHkU)(10)where Hk is the Hamiltonian of Eq. 8 and UUC2z is the matrix representation of the twofold rotation C2z. Note that (51) uses a gauge for which Hk + G = Hk holds, where G is a reciprocal lattice vector. Equation 10 implicitly assumes a tight-binding gauge, in which case matrix representations of symmetries, in particular U, are momentum independent.

The Hamiltonian of Eq. 10 with the parameters specified in (51) defines a tight-binding model for TBG that respects all symmetries, including a Uv(1) valley conservation symmetry. Various symmetry-breaking terms can be considered, and here, we are specifically interested in terms that break C3z symmetry but preserve C2z symmetry, as required by quadrupolar nematic order.

Note that a coupling of the valleys of the form δH=Δkψk+ψk+H.c. breaks the Uv(1) valley but respects all lattice symmetries as well as time-reversal symmetry. Added to the Hamiltonian of Eq. 10 it enters as an off-diagonal block.

Rotational symmetry breaking

Consider next the rotational symmetry-breaking terms. Since the model is built from multiple degrees of freedom, there are a number of different ways in which rotation symmetry breaking can be implemented. We first focus on the triangular lattice sector of the model. Within this sector, there are two possibilities: Nematic order can occur as a result of hopping anisotropy or because of a lifting of the orbital degeneracy. To model the first possibility, we introduce the two d-wave form factorsdk1=ϕ01+Re ω* ϕ1¯1¯+Re ω ϕ10,+c.c.dk2=Im ω* ϕ1¯1¯+Im ω ϕ10+c.c.(11)where the phases ϕlm are defined as ϕlm = eik · (la1 + ma2) [see (51)]. Here, we use the notation l¯l and ω = exp (2πi/3). These form factors have precisely the same symmetry as (Φ1, Φ2) introduced in the main text.

The triangular lattice d-wave form factors can then be used to introduce a symmetry-breaking perturbation in the triangular lattice (p-orbital) sector. For instance, we can add a perturbation δHpz to the Hamiltonian Hpz of the pz orbital appearing in Eq. 8 given byδHpz=Φ1dk1+Φ2dk2(12)where Φ1,2 are the nematic order parameters as defined in the main text. This term gives the Fermi surface distortion of Fig. 4B of the main text. Clearly, the same perturbation (but proportional to the appropriate identity matrix) can be added to Hp±, which describes the p± orbitals.

In the p±-orbital sector, the nematic order parameter couples to another symmetry-breaking perturbation, which is independent of momentum. Making the two orbitals inequivalent lifts their degeneracy and necessarily breaks threefold rotation symmetry. In particular, the perturbation δHp±, which achieves this, is given byδHp±=(0Φ1iΦ2Φ1+iΦ20)(13)

Note that the diagonal terms are zero, since time-reversal symmetry must be preserved. (An overall energy can be absorbed in μp±.) The Fermi surface distortion due to Eq. 13 is qualitatively similar to a distortion originating from d-wave form factors in the kinetic terms (Fig. 4B of the main text).

Consider next the kagome lattice sector of the model. The kagome sector does not have an orbital degree of freedom, but it does have multiple sites in the unit cell. The simplest coupling to the nematic order parameter is given by a charge ordering perturbation within the unit cell, which breaks threefold rotations but preserves the twofold rotation C2z. Specifically, the perturbation δHκ to the kagome lattice Hamiltonian Hκ is given byδHκ=(Φ1iΦ2)(1ωω*)+H.c.(14)

The rotational symmetry-breaking perturbations introduced so far are all intravalley perturbations; they should be considered as perturbations to Eq. 8, with the full Hamiltonian given by the prescription of Eq. 10. One may, however, also consider intervalley nematic coupling terms, which enter the off-diagonal blocks in Eq. 10. More precisely, the Hamiltonian of Eq. 10 is modified according toHk(HkUHkU)+δHΦ(15)where δHΦ collects all terms that describe nematic distortions and takes the formδHΦ=(ΔΦΔΦ)(16)

The form of ΔΦ depends on the choice of nematic coupling; as in the case of intravalley nematic coupling, in principle, many possibilities of intervalley nematic coupling exist. One simple type of nematic coupling is given byΔΦ=δHpzδHp±(17)where δHpz and δHp± are given by Eqs. 12 and 13.

We have based our microscopic discussion of rotation symmetry breaking on the model introduced in (51). This was motivated by the natural implementation of all relevant symmetries in this model. It is worth stressing that an analysis of nematic order in TBG similar to the one presented here can also be obtained from different microscopic (tight-binding) models proposed for TBG (1114, 21).

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We thank A. Chubukov, L. Fu, P. Jarillo-Herrero, J. Kang, H. Ochoa, H. C. Po, J. Schmalian, T. Senthil, and O. Vafek for fruitful discussions. Funding: R.M.F. was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award no. DE-SC0020045. Author contributions: R.M.F. and J.W.F.V. conceived the problem, performed the analysis and calculations, and contributed to writing the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article