## Abstract

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 (*C*_{3z}) of the moiré superlattice. We find that, even in the presence of static strain that explicitly breaks the *C*_{3z} 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.

## INTRODUCTION

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 (*1*–*7*). 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 (*8*–*21*), 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) (*24*–*26*) and transport measurements (*27*) reported evidence that the threefold rotational symmetry of the moiré superlattice, denoted by *C*_{3z}, is broken in different regions of the TBG phase diagram, including close to the superconducting dome. Moreover, spontaneous *C*_{3z} 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 (*30*–*34*), 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*, *35*–*37*).

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 *C*_{3z}-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*).

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.

## RESULTS

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 *d*_{x2−y2} and *d _{xy}* 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

*D*

_{6}: The two nematic components belong to a single irreducible representation of

*D*

_{6}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

**Φ**(θ) =

**Φ**(θ + π), as expected. Note that all of this remains true if one considers a TBG model with

*D*

_{3}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 *S*_{nem}[**Φ**] is [see also (*30*, *40*–*42*)]*x* = (**r**, τ) denotes spatial coordinate **r** and imaginary time τ, and Φ_{±} ≡ Φ_{1} ± *i*Φ_{2}. The first term, ^{4} action with *U*(1) symmetry. The cubic term reflects the crystalline anisotropy of the hexagonal lattice and is expressed as *n*π/6 for γ < 0 and θ = (2*n* + 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 Z_{3} symmetry is identified with the threefold rotation about the *z* axis, *C*_{3z}. Below the nematic transition temperature *T*_{nem}, the sixfold rotation symmetry *C*_{6z} about the *z* axis is broken down to a twofold symmetry *C*_{2z}, while the twofold rotation symmetries *C*_{2x} and *C*_{2y} 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 *C*_{3z} is broken because of condensation of the nematic order parameter **Φ**, all order parameters that also break *C*_{3z} 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 **u**. In the present case, the electronic nematic order parameter **Φ** may refer to any fermionic bilinear that breaks *C*_{3z} symmetry, which could be in the charge, valley, or spin channel. On the other hand, ε* _{ij}* and ω

*refer only to the lattice geometric degrees of freedom. The elasto-nematic action is given by*

_{ij}We consider first the effect of static strain applied in TBG. For compressive (tensile) uniaxial strain ε < 0 (ε > 0) applied parallel to an arbitrary direction *S*′ = − λ∫* _{x}*ε Φ cos(2α − 2θ), where

*T*≫

*T*

_{nem}, where

*T*

_{nem}is the transition temperature of the unstrained system, we can approximate

_{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

*T*≪

*T*

_{nem}and set ∣

**Φ**∣ = Φ

_{0}as approximately constant. Expanding around the high-temperature director, θ = θ

_{0}+ δθ, gives

*n*π/6). It follows that

*a*

_{θ0}= 0 and that the system remains symmetric under twofold rotations with respect to in-plane axes,

*C*

_{2x}and

*C*

_{2y}. 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

*S*

_{nem}. In this case, once Φ

_{0}reaches the critical value

_{0}to

*C*

_{2x}and

*C*

_{2y}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.

From the perspective of bond order, the breaking of the twofold *C*_{2x}/*C*_{2y} symmetries (in addition to *C*_{3z}), makes all three bonds connecting nearest-neighbor sites inequivalent. Recall that, in the nematic phase, where only the *C*_{3z} 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

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 *C*_{3z} but also *C*_{2x}/*C*_{2y} symmetry. Consequently, no true phase transition can happen. Nevertheless, near

### Fluctuating strain

Apart from static strain, finite-momentum elastic fluctuations are expected to strongly affect the nematic transition (*44*–*46*). For a system with *D*_{6} symmetry, diagonalization of the harmonic elastic action *T*) and a longitudinal (*L*) one with sound velocities *v*_{L, T}*q* = (**q**, ω* _{n}*), with ω

*the (bosonic) Matsubara frequency. The displacement field*

_{n}_{q}= arctan (

*q*/

_{y}*q*). According to (

_{x}*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*,

*47*–

*49*). Integrating out the acoustic phonons leads to an additional contribution to the nematic action,

*= 0, we find*

_{n}*T*

_{nem}. 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 coupling

^{2}cos

^{2}(2θ − 2ζ

_{q}), we see that it makes only certain directions of momentum space to become soft; i.e., the static nematic susceptibility

*T*

_{nem}only along special directions

*45*)]. While the cubic term in Eq. 6 forces the director

_{q}soft, namely, the ones that make a relative angle of 0 and

^{4}action

*S*

_{0}[

**Φ**] 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.

### 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 *c*_{k}, the electronic nematic coupling is *g*(**k**) = *g*_{0} cos(2θ − 2θ_{k}), where *g*_{0} is a constant and θ_{k} = arctan (*k _{y}*/

*k*). The electronic states that exchange low-energy nematic fluctuations are at the Fermi surface and separated by the small momentum

_{x}**q**of the nematic mode. Since the nematic fluctuations are the softest (albeit nondiverging) along the special directions

**k**

_{hs}where the Fermi surface’s tangent is parallel to

*g*(

**k**

_{hs}) is. For a circular Fermi surface, the hot spots are located at

*g*(

**k**

_{hs}) = 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

*g*(

**k**

_{hs})∣ = ∣

*g*

_{0}∣; i.e., the hot spots are maximally coupled to the soft nematic fluctuations. For a generic noncircular Fermi surface respecting

*D*

_{6}symmetry

*g*(

**k**

_{hs}) remains maximum for η < 0 but is expected to be nonzero albeit small for η > 0.

The sign of *C*_{11} ≡ *C _{xxxx}* and

*C*

_{12}≡

*C*, yielding

_{xxyy}*C*

_{11}> ∣

*C*

_{12}∣, which makes η > 0, implying that

*g*(

**k**

_{hs}) 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,

*39*). This term contributes only to the transverse velocity and when

*K*> 2(

*C*

_{11}+

*C*

_{12}),

*v*becomes larger than

_{T}*v*(i.e., η < 0), implying that ∣

_{L}*g*(

**k**

_{hs})∣ 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 *C*_{2z} 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).

## DISCUSSION

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 *C*_{3z} 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*, *54*–*58*). These orbital variables are associated with eigenstates of the *C*_{3z} 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 *C*_{3z} 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 *C*_{3z} 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 *C*_{3z} 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*, *47*–*49*).

## MATERIALS AND METHODS

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 (*p _{z}*,

*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

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 by

Here, *H _{pz}* and

*H*

_{p±}are the subblock Hamiltonians in the

*p*and

_{z}*p*

_{±}subspaces, and

*H*

_{κ}describes the coupling between the kagome lattice sites. The subblocks

*C*describe the couplings between the

_{XY}*X*and

*Y*sectors (where

*X*,

*Y*=

*p*,

_{z}*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* and *K*′ valleys. The full Hamiltonian *H*_{k} is the Hamiltonian of Eq. 8 and *U* ≡ *U*_{C2z} is the matrix representation of the twofold rotation *C*_{2z}. Note that (*51*) uses a gauge for which *H*_{k + G} = *H*_{k} 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 *U _{v}*(1) valley conservation symmetry. Various symmetry-breaking terms can be considered, and here, we are specifically interested in terms that break

*C*

_{3z}symmetry but preserve

*C*

_{2z}symmetry, as required by quadrupolar nematic order.

Note that a coupling of the valleys of the form *U _{v}*(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 factors* _{lm}* are defined as ϕ

*=*

_{lm}*e*

^{−ik · (la1 + ma2)}[see (

*51*)]. Here, we use the notation

*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 δ*H _{pz}* to the Hamiltonian

*H*of the

_{pz}*p*orbital appearing in Eq. 8 given by

_{z}_{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

*H*

_{p±}, 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 δ*H*_{p±}, which achieves this, is given by

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 *C*_{2z}. Specifically, the perturbation δ*H*_{κ} to the kagome lattice Hamiltonian *H*_{κ} is given by

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 to_{Φ} collects all terms that describe nematic distortions and takes the form

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*H _{pz}* and δ

*H*

_{p±}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 (*11*–*14*, *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.

## REFERENCES AND NOTES

**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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).