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

Twisted 2D materials realize a unique electronic nematic state, quite different from its counterpart in bulk quantum materials.

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)(25)(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)(31)(32)(33)(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)(36)(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 nematicflop 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 lowenergy 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 x 2 −y 2 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 ˆ n = (cos, sin) (see Fig. 1D); note that () = ( + ), as expected. Note that all of 1 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)(41)(42) where x = (r, ) denotes spatial coordinate r and imaginary time , and  ± ≡  1 ± i 2 . The first term, S 0 [  ] = 1 _ 2 r  || 2 + 1 _ 4 u  || 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 1 _ 3   3 cos6 , 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 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 cellpreserving 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  ij ≡ 1 _ 2 ( ∂ i u j + ∂ j u i ) and the rotation tensor  ij ≡ 1 _ 2 ( ∂ i u j − ∂ j u i ) , 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 C 3z 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 is the elastic free energy and 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 x . At high temperatures T ≫ T nem , where T nem is the transition temperature of the unstrained system, we can approximate S nem ≈ 1 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 T ≪ T nem and set || =  0 as approximately constant. Expanding around the high-temperature director,  =  0 + , gives . 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, 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 ¯ , resulting in an Ising-like transition that spontaneously breaks the 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 T nem flop , 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 nematicflop 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 T nem flop , 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 (44)(45)(46). For a system with D 6 symmetry, diagonalization of the harmonic elastic action S el [ ˆ  , ˆ  ] leads to two acoustic phonon modes, a transverse (T) and a longitudinal (L) one with sound velocities v L, T where q = (q,  n ), with  n the (bosonic) Matsubara frequency. The displacement field u = ∑  u ̃  ˆ e  has been decomposed into its longitudinal and transverse components u ̃  with ˆ e L = (cos  q , sin  q ) and ˆ e T = (− sin  q , cos  q ) and  q = arctan (q y /q x ). According to (38,39), for the dominant acoustic phonons that act on the moiré superlattice scale, u corresponds to the relative displacement of the marks an Ising-like nematic-flop transition in which the twofold rotational symmetries C 2x and C 2y are spontaneously broken. The sixfold rotation C 6z is explicitly broken to C 2z 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 C 2x and C 2y . two graphene sheets. These and other phonon modes have been proposed to be linked to superconductivity in TBG (18,(47)(48)(49). Integrating out the acoustic phonons leads to an additional contribution to the nematic action,  S nem = − 1 _ 2 ∑ ij ∫ q  i,q ˆ  ij (q )  j,−q . In the static limit,  n = 0, we find where ˆ is the identity matrix and  ≡ 1 − v T 2 / v L 2 . The first term of ˆ  gives an overall enhancement of 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 nematoorbital coupling where is the momentum space (i.e., orbital) quadrupolar form factor. Recasting the nemato-orbital coupling term as  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  nem (q → 0, ˆ q ) is largest near T nem 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 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 threestate 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 S elec = ∫ k,q g(k)  q c k−q/2 † c k+q/2 , with form factor g(k) = g 0 cos(2 − 2 k ), where g 0 is a constant and  k = arctan (k y /k x ). 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" k hs where the Fermi surface's tangent is parallel to ˆ q (0) , i.e., ˆ q (0) ·  k hs = 0 . The issue is how strong these fermions are coupled to the nematic fluctuations, i.e., what the magnitude of g(k hs ) is. For a circular Fermi surface, the hot spots are located at  k hs =  q (0) +  _ 2 , and thus g( k hs ) = − g 0 cos(2 − 2  q (0) ) . As we saw above, if  > 0, then the soft directions are  q (0) =  ±  / 4 , yielding 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  q (0) = ,  ±  / 2 , implying that |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 depends only on the strain ˆ  , since global rotations do not cost energy. In a triangular lattice, there are only two independent elastic constants, C 11 ≡ C xxxx and C 12 ≡ C xxyy , yielding v L 2 = C 11 and v T 2 = ( C 11 − C 12 ) / 2 . Lattice stability requires 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,  S el [ ˆ  ] = 1 _ 2 ∫ x K  xy 2 (39). This term contributes only to the transverse velocity and when K > 2(C 11 + C 12 ), v T becomes larger than v L (i.e.,  < 0), implying that |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 Fig. 3. Phonon-mediated nemato-orbital coupling. Momentum directional dependence of the nematic susceptibility  nem (q → 0, ˆ q ) caused by the nematoorbital 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, v T < v L , the soft direction is rotated by ±/4 with respect to the nematic director n (red arrow), for TBG, v T > v L , the rotation is 0, /2. Note that the director can point in any of the directions  of Fig. 1B. 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 weakcoupling 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)(55)(56)(57)(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 strongcoupling 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)(48)(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 Here, H p z and H p ± are the subblock Hamiltonians in the p z and p ± subspaces, and H  describes the coupling between the kagome lattice sites. The subblocks C XY describe the couplings between the 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 where ± labels the K and K′ valleys. The full Hamiltonian H k is then given by where H k is the Hamiltonian of Eq. 8 and U ≡ U C 2z 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 symmetrybreaking 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 H =  ∑ k  k+ †  k− + H . c. breaks the 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 where the phases  lm are defined as  lm = e −ik · (la 1 + ma 2 ) [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 H p z to the Hamiltonian H p z of the p z orbital appearing in Eq. 8 given by 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 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 7 of 8 Eq. 10. More precisely, the Hamiltonian of Eq. 10 is modified according to where H  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 where 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)(12)(13)(14)21).