# Magic angle conditions for twisted 3D topological insulators

Aaron Dunbrack  Department of Physics and Astronomy, Stony Brook University,
Stony Brook, New York 11974, USA  Jennifer Cano  Department of Physics and
Astronomy, Stony Brook University, Stony Brook, New York 11974, USA  Center
for Computational Quantum Physics, Flatiron Institute, New York, New York
10010, USA

###### Abstract

We derive a general low-energy theory for twisted moiré heterostructures
comprised of Dirac materials. We apply our theory to heterostructures on the
surface of a three dimensional topological insulator (3D TI). First, we
consider the interface between two 3D TIs arranged with a relative twist
angle. We prove that if the two TIs are identical, then a necessary condition
for a perturbative magic angle where the Dirac cone velocity vanishes is to
have an interlayer spin-flipping hopping term. Without this term, the Dirac
cone velocities can still be significantly renormalized, decreasing to less
than half of their original values, but they will not vanish. Second, we
consider graphene on the surface of a TI arranged with a small twist angle.
Again, a magic angle is unachievable with only a spin-flipping hopping term;
under such conditions, the Dirac cone is renormalized, but only moderately. In
both cases, our perturbative results are verified by computing the band
structure of the continuum model. The enhanced density of states that results
from decreasing the surface Dirac cone velocity provides a tunable route to
realizing interacting topological phases.

##  I Introduction

Magic angle twisted bilayer graphene exhibits a variety of interacting phases
such as superconductivity and the quantum anomalous Hall effect [1, 2, 3, 4].
The relative twist angle between the two layers acts as a tuning knob that can
dramatically change their physical properties. Specifically, at specific
“magic” twist angles [5, 6, 7], the renormalized Dirac cone velocity vanishes,
leaving behind gapped flatbands that are susceptible to interacting
instabilities.

In this manuscript, we ask whether the physics observed in twisted graphene
heterostructures can be realized in other Dirac materials. We are particularly
motivated by renormalizing the velocity of the Dirac cone on the surface of a
topological insulator (TI), where an interacting gap will yield either a
quantum anomalous Hall insulator or a topological insulator [8]. Recently, the
effect of a slowly spatially-varying “moiré superlattice” potential on the
surface of a topological insulator has been studied [9, 10, 11]. Such a
potential will increase the surface density of states at the Dirac point by
renormalizing the Dirac cone. It will also generate van Hove singularities at
an energy slightly above and below the Dirac point. Both of these effects will
enhance the instabilities towards interacting topological phases [12, 13, 14,
15, 16, 17, 18, 19, 10, 11]. A potential that breaks time-reversal symmetry
may produce flat Chern bands [20]. The theory can be generalized to modulate
an interacting gap on the edge of a 2D TI [21].

Here, we consider a different route to manipulating the surface Dirac cone
velocity, which is to incorporate a 3D TI into a twisted heterostructure with
another Dirac material. The interlayer coupling between the two Dirac cones
offers a new tuning knob compared to the superlattice potential and, as we
show, can lead to perfectly flat bands under certain ideal (but not fine-
tuned) conditions.

To this end, we derive a general theory of a twisted moiré heterostructure
between two Dirac materials. The two layers need not be identical, but we
require them to feature Dirac cones which, upon being arranged with a small
twist angle, are separated by a small momentum difference. We derive an
analytical expression for the low-energy theory as a function of the
interlayer coupling and discuss the effect of rotation, time-reversal, and
interlayer mirror symmetries.

We then apply our theory to two types of twisted heterostructures on the
surface of a topological insulator. First, we consider the two-dimensional
(2D) interface between two three-dimensional (3D) TIs arranged with a small
relative twist angle. To meet our criterion that the two Dirac cones must be
separated by small momentum after twisting, we require that the Dirac cones
are not at the center of the Brillouin zone (BZ), but rather at the
$(\pi,\pi)$ point. We then derive conditions for gapless flat bands. Our most
significant result is that if the two 3D TIs are identical, a gapless flat
band requires an interlayer spin-flipping hopping term. We prove that this
condition is mathematically equivalent to the known condition in twisted
bilayer graphene that an interlayer sublattice-swapping is necessary to
achieve a magic angle [22]. However, while the latter is natural in twisted
bilayer graphene, the origin of a spin-flipping term is not clear, although it
is symmetry-allowed. Nonetheless, if the spin-flipping term is not present,
the Dirac cone velocity can still be significantly renormalized, reaching a
minimum around half of its original value. We back our analytical results with
a numerical simulation.

The second heterostructure we examine in detail is graphene on the surface of
a topological insulator, arranged with a small relative twist angle. Graphene
on a 3D TI surface has been studied extensively without a twist angle, both
theoretically [23, 24, 25, 26, 27, 28, 29, 30, 31, 32] and experimentally [33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. A twist angle was considered in
Refs. 26, 42, but flat bands were not considered. The primary goal of the
previous literature was to enhance spin-orbit coupling in graphene. Our
motivation is in the opposite direction: we hope the interlayer coupling from
graphene will renormalize the Dirac cone on the 3D TI. Similar to the
interface between two topological insulators, we find that a magic angle
condition is not achieved without spin-flipping interlayer hopping, although
the 3D TI Dirac cone can still be slightly renormalized. We verify our
analytical results with a numerical model.

##  II Magic Angle Condition for Dirac Cones at TRIM Points

We are interested in twisted bilayer heterostructures of two Dirac materials.
Although our theory is applicable to any Dirac material, we focus on the case
where one layer has a single time-reversal-invariant Dirac cone at a TRIM
point, which implies it is the surface of a strong 3D TI. We write down the
most general interlayer coupling Hamiltonian and compute the self-energy of
each Dirac cone perturbatively in the coupling. From the self-energy, we
extract both an energy shift and a renormalized Fermi velocity for the Dirac
cone, the latter of which gives a magic angle condition. We then discuss the
physical conditions necessary to achieve the magic angle. Although we use the
language specific to a twisted moiré heterostructure, the theory presented
applies equally well to a moiré heterostructure created by stacking two layers
with a small lattice mismatch.

###  II.1 Model for Dirac materials

A Dirac Hamiltonian for an isotropic gapless Dirac cone at charge neutrality
can be written in the form:

$H_{D}(k)=vk\cdot\sigma,$ (1)

where $\sigma$ is some here-unspecified degree of freedom such as spin,
sublattice, or orbital. (The proof that Eq. (1) is fully general is given in
Appendix A.) When the Dirac cone is not at the origin, we will denote the
momentum difference from the Dirac point at $k_{0}$ by $\bar{k}\equiv
k-k_{0}$. In each layer, there may be multiple Dirac cones of the form of Eq.
(1), which are labelled by valley, spin, or other indices. Unless symmetry-
related, each Dirac cone will have a different Fermi velocity and energy, as
well as anisotropy; for simplicity here we do not vary these parameters.

The Hamiltonian for each layer $L=1,2$ is written as
$H(\bar{k})=H_{D}(\bar{k})\delta_{ij}$ (where $i$, $j$ run over the different
Dirac cones, e.g., in graphene, over spin and valley):

$\hat{H}_{0}=\int
d^{2}\bar{k}\psi_{L,\bar{k}}^{\dagger}H(\bar{k})\psi_{L,\bar{k}}$ (2)

Interlayer coupling is then included as:

$\hat{T}=\sum_{Q}\hat{T}_{Q}=\sum_{Q}\int
d^{2}\bar{k}\hat{\psi}_{1,\bar{k}}^{\dagger}T_{Q}\hat{\psi}_{2,\bar{k}+Q}+h.c.$
(3)

Here, $\hat{\psi}_{L,\bar{k}}$ are the electron field operators on layer $L$
at crystal momentum $k_{0}+\bar{k}$ with spin, orbital, valley, etc. indices
implicit; note $k_{0}$ may also implicitly depend on valley and layer.

The vectors $Q$ are the displacements in momentum space of the Dirac cones
that result from twisting (see, e.g., Fig. 1). The collection of vectors
$\\{Q_{i}-Q_{j}\\}$ from each valley form the reciprocal lattice of the moirè
pattern and thereby define the moirè BZ. A derivation of which $Q$ enter and
how to compute the corresponding interlayer coupling matrix $T_{Q}$ is given
in Appendix B. Under the assumption that the lattices are nearly matched and
arranged with a small twist angle, and, in addition, that the interlayer
spacing is appreciably larger than the lattice constant of either material,
$T_{Q}$ is only non-negligible for a finite set of $Q$.

The low-energy continuum model in Eqs. (2) and (3) is valid when $\bar{k}$ is
small (i.e., close enough to a Dirac cone that the linearization is valid) and
when the valleys are well-separated (i.e., the distance between distinct
valleys is much larger than $Q$).

To study the effect of the interlayer coupling term in Eq. (3) on the Dirac
Hamiltonian in Eq. (2), we assume that the interlayer coupling is small, in
the sense that the energy scale of interlayer coupling, $t$, is much less than
$v|Q|$ for the smallest value of $|Q|$. (This implicitly assumes that there is
no $Q=0$ term; the possibility of such a term is discussed further in Sec.
III.) We then work perturbatively in the parameter $t/v|Q|$.

We further assume the presence of both time-reversal symmetry, $\mathcal{T}$,
with $\mathcal{T}^{2}=-1$, and an in-plane rotational symmetry by $2\pi/n$,
denoted $C_{n,z}$, with $n>2$. The rotation symmetry ensures that the Dirac
cone is isotropic, which is a convenient and physically relevant
simplification, but is not necessary. If the two layers are the same, there
may be additional symmetries that exchange the layers. The action of the
symmetry operators is derived in Appendix A.

###  II.2 Self-energy of a single Dirac cone

We now treat the interlayer coupling in Eq. (3) as a perturbation to the Dirac
Hamiltonian in Eq. (2) and compute the self-energy to lowest order in
$t/v|Q|$. From the self-energy, we will extract the renormalized velocity and
energy shift of the original Dirac cone due to interlayer coupling. We compute
the self-energy of a single Dirac cone located at a time-reversal-invariant-
momentum (TRIM); in Appendix C we discuss the generalization for a system with
multiple Dirac cones not necessarily at TRIM. Without loss of generality, we
take the layer in which we are computing the self-energy to be layer $1$, with
Fermi velocity $v_{1}$, while layer 2 has Fermi velocity $v_{2}$.

In this single Dirac cone case, the combination of time-reversal and
rotational symmetry constrain the self-energy to take the form:

$\Sigma^{\text{single-
cone}}(\omega,\bar{k})=-[A\omega+B+C(\bar{k}\cdot\sigma)+D(\bar{k}\crossproduct\sigma)]$
(4)

The parameter $A$ determines the quasiparticle weight via $Z^{-1}=1+A$, while
the parameters $C$ and $D$ are direct corrections to the renormalized Fermi
velocity, which takes the form:

$v^{*}_{1}=\frac{\sqrt{(v_{1}-C)^{2}+D^{2}}}{1+A},$ (5)

where $v_{1}$ is the original velocity before the interlayer coupling. The
parameter $B$ in Eq. (5) determines an overall energy shift of the
renormalized Dirac cone, given by:

$\Delta E=-\frac{B}{1+A}$ (6)

Eq. (5) is an important result: it imposes a condition to find a “magic angle”
where the renormalized Fermi velocity vanishes, specifically, $v^{*}_{1}=0$
when $D=0$ and $C=v_{1}$. In contrast, previous work [9] showed that a
superlattice potential on the surface of a 3D TI can renormalize the Dirac
cone velocity, but does not result in a magic angle condition; in the language
of Eq. (5), it yields $C=D=0$.

We now derive the physical conditions that must be satisfied to realize the
magic angle by expressing the abstract parameters $A,B,C,D$ in terms of the
physical interlayer coupling terms, $T_{Q}$, introduced in Eq. (3).

For simplicity, we first assume the original Dirac cone is coupled to a single
Dirac cone in an adjacent layer, so that $T_{Q}$ is a $2\times 2$ matrix. For
each $Q$, we decompose $T_{Q}$ as follows:

$T_{Q}=t_{0,Q}\mathbb{I}+t_{z,Q}\sigma_{z}+t_{||,Q}(\hat{Q}\cdot\sigma)+t_{\perp,Q}(\hat{Q}\times\sigma),$
(7)

where the parameters $t_{0,z,||,\perp}$ are complex numbers. Note that since
$\hat{Q}$ is ill-defined at $Q=0$, $t_{||}$ and $t_{\perp}$ are not continuous
through $\theta=0$; specifically, the behavior of $\hat{Q}$ implies that
$t_{||,\perp}\sim\text{sgn}(\theta)$, to lowest order in $\theta$.

In terms of these coefficients,

$\displaystyle A$
$\displaystyle=\sum_{Q}\frac{1}{v_{2}^{2}Q^{2}}(|t_{0,Q}|^{2}+|t_{z,Q}|^{2}+|t_{||,Q}|^{2}+|t_{\perp,Q}|^{2})$
(8a) $\displaystyle B$
$\displaystyle=\sum_{Q}\frac{2}{v_{2}|Q|}(\real[t_{0,Q}t_{||,Q}^{*}]+\imaginary[t_{\perp,Q}t_{z,Q}^{*}])$
(8b) $\displaystyle C$
$\displaystyle=\sum_{Q}\frac{1}{v_{2}Q^{2}}(|t_{\perp,Q}|^{2}-|t_{||,Q}|^{2})$
(8c) $\displaystyle D$
$\displaystyle=\sum_{Q}\frac{-2}{v_{2}Q^{2}}\real[t_{||,Q}t_{\perp,Q}^{*}]$
(8d)

These equations follow directly from the more general results presented in
Appendix C.1.

More generally, if the original Dirac cone is coupled to multiple Dirac cones
in the second layer, labelled by another (valley, spin, etc.) index, then the
decomposition in Eq. (7) should be applied to each valley, and the parameters
$A,B,C,D$ are found by summing Eq. (8) over all valleys.

In each of the sums in Eq. (8), all symmetry-related $Q$ values will
contribute the same summand, due to the definition of the coefficients in Eq.
(7) and the symmetry properties discused in Appendix B.3. Accordingly, if only
the smallest $Q$ values are non-negligible, it suffices to consider the
couplings for one choice of $Q$.

##  III Interface between two 3D topological insulators

The first case we will consider is the interface between two 3D TIs that are
stacked on top of each other with a small relative twist angle. This is also
the simplest case, since each surface hosts only a single time-reversal-
invariant Dirac cone.

If both layers have a Dirac cone at $\Gamma$, there will be a zero-momentum
($Q=0$) interlayer coupling term. Physically, this term appears because even
after twisting, the Dirac cones at $\Gamma$ have no momentum separation.
(Mathematically, the Q=0 term arises from our derivation in Appendix B.) The
presence of a $Q=0$ term invalidates our results in the previous section,
which were derived by expanding perturbatively in $t/v|Q|$. Instead of a
renormalized gapless Dirac cone, a $Q=0$ term will open a gap at the
interface, similar to the interface between two untwisted 3D TIs, but
modulated by a superlattice potential. Therefore, to realize a gapless Dirac
cone at an interface between two twisted 3D TIs, we require that at least one
layer has a Dirac cone at a TRIM that is not $\Gamma$.

This motivates us to consider the hypothetical situation of an interface
between two $C_{4}$-symmetric 3D TIs with the same lattice constant and Dirac
cones at the surface $M\equiv(\pi,\pi)$ point. When the two materials are
twisted, the surface BZ is illustrated in Fig. 1(a). Since the Dirac cones are
at $(\pi,\pi)$, they are $C_{4}$-invariant and thus isotropic to linear order
(any rotation of order greater than two will yield an isotropic Dirac cone.)
Without such a symmetry, the Dirac cones can flatten anisotropically; if the
velocity vanishes in only one direction, an effective 1D system results [45,
46]. The other TRIM on the surface of a four-fold-symmetric TI are not
$C_{4}$-symmetric: the $(\pi,0)$ and $(0,\pi)$ points mix under a four-fold
rotation. On the surface of a three- or six-fold symmetric TI, the TRIM also
mix under the three- or six-fold rotational symmetry, such that there are no
TRIM points (besides $\Gamma$) that are three- or six-fold invariant. Thus, if
we seek a Dirac cone at a TRIM that will be renormalized isotropically at a
twisted interface between two 3D TIs, our only option is to consider a Dirac
cone at $(\pi,\pi)$ on the surface of a four-fold symmetric topological
insulator.

We may also consider interfaces between two 3D TIs with different lattice
constants. For example, consider an interface between one 3D TI with a Dirac
cone at $M$ and a second 3D TI with a Dirac cone at $\Gamma$, such that the
two materials have a 2:1 lattice vector mismatch (their unit cell areas differ
by a factor of four). In this setup, the cone at $M$ from the first layer is
folded onto $\Gamma$ from the second when the layers are aligned. When the two
materials are twisted, the surface BZ is illustrated in Fig. 1(b). The low-
energy theory of the interface is similar to the case where the two layers
have the same lattice constant and both have Dirac cones at $M$, but subject
to different symmetries. Further, the global topology, i.e., where the Dirac
points are located on the surface BZ, may differ, which gives rise to
different mechanisms for the topological protection of gapless surface states.
One can generalize to construct interfaces with more complicated supercells
describing any integer ratio of lattice constants.

(a) Surface Dirac cones at $M$ in both layers, with identical lattice
constants.

(b) Dirac cone at $M$ on one surface and $\Gamma$ on the other, with lattice
constants that differ by a factor of two.

Figure 1: BZ at the interface between two 3D TIs arranged with a small
relative twist for two scenarios described in Sec. III. Filled circles
indicate Dirac cones; empty circles indicate other TRIM in the original BZ.
Vectors $Q$ indicate the momentum difference between nearest Dirac cones,
which is exactly the effective quasimomentum transfer in the continuum model,
per Eq. (3). Figure 2: Moiré BZ. Blue and red filled circles are the locations
of Dirac cones from each layer. Empty circles indicate other points along the
BZ path (which may be TRIM points depending on the configuration.)

###  III.1 Identical TI surfaces

We start by considering the first scenario, where both TIs are identical and
have Dirac cones at $M$. Then $C_{2}\mathcal{T}$ symmetry leaves each $Q$
invariant (as we derive in Eqs (45) and (46), it acts as
$-\sigma_{x}\mathcal{K}$, implying $\sigma_{x}T_{Q}^{*}\sigma_{x}=T_{Q}$.)
This constrains the coefficients $t_{0}$, $t_{||}$ and $t_{\perp}$ to be real,
whereas $t_{z}$ is purely imaginary. Accordingly, let us define
$t_{z}=i\hat{t}_{z}$.

We assume that only the smallest symmetry-related set of four $Q$ contribute
nonnegligibly. Without loss of generality, we may then choose these $Q$s along
the $x$ and $y$-axes, as illustrated in Fig. 1(a). This choice of $Q$s is
still valid in the presence of a small lattice mismatch, although the
decomposition into $t_{||}$ and $t_{\perp}$ will change. Eq. (8) simplifies
to:

$\displaystyle A$
$\displaystyle=\frac{4(t_{0,Q}^{2}+\hat{t}_{z,Q}^{2}+t_{||,Q}^{2}+t_{\perp,Q}^{2})}{v^{2}Q^{2}}$
(9a) $\displaystyle B$
$\displaystyle=\frac{8(t_{0,Q}t_{||,Q}+t_{\perp,Q}\hat{t}_{z,Q})}{v|Q|}$ (9b)
$\displaystyle C$
$\displaystyle=\frac{4(t_{\perp,Q}^{2}-t_{||,Q}^{2})}{vQ^{2}}$ (9c)
$\displaystyle D$ $\displaystyle=\frac{-8t_{||,Q}t_{\perp,Q}}{vQ^{2}}$ (9d)

One of the most important consequences of these formulas is that if $t_{||}$
and $t_{\perp}$ vanish, then the only non-zero term is $A$. However, according
to the formula for the renormalized velocity in Eq. (5), a magic angle
requires $C>0$. Thus, to lowest order in perturbation theory, a magic angle in
this system requires $t_{\perp}$ to be non-zero. Since, by its definition in
Eq. (7), $t_{\perp}$ is a coefficient of off-diagonal Pauli matrices, we
conclude that magic angles are only possible in the presence of inter-layer
spin-flipping hopping terms.

This result can be understood by analogy to twisted bilayer graphene, where
instead of spin, the $\sigma$ matrices act on the sublattice degree of
freedom. In Ref. [22], the authors showed that magic angles do not appear in
the “second chiral limit” of TBLG where inter-sublattice interlayer hopping
vanishes. When sublattice is replaced by spin, this is exactly our result that
magic angles do not appear at the interface between two identical 3D TIs
without spin-flipping interlayer hopping. Although in TBLG the Dirac cones are
in sublattice space and protected by spinless $C_{2}\mathcal{T}$, while on the
surface of a 3D TI, the Dirac cone is in spin space and protected by
$\mathcal{T}$, we show in Appendix D.1 that this analogy is mathematically
rigorous.

We now consider the magic angle constraint in more detail. If the two 3D TI
layers are identical, they have an additional layer-interchanging symmetry
$C_{2y}$ (which, in combination with $C_{4}$, generates three other layer-
interchanging symmetries). Using the orientation of $Q$s defined in Fig. 1
(along with the action of the rotation operator defined in Eq. 48), this
symmetry requires:

$\sigma_{y}T_{Q_{1}}^{\dagger}\sigma_{y}=T_{Q_{1}}$ (10)

Combining the above equation with the decomposition of $T_{Q}$ given in Eq.
(7) (and the aforementioned $C_{2}\mathcal{T}$ constraints) yields the single
additional constraint $t_{||}=0$, which, from Eq. (9), implies $D=0$,
$C=4t_{\perp,Q}^{2}/vQ^{2}$. According to the expression for the renormalized
Fermi velocity in Eq. (5), a magic angle will result when $D=0$, $C=v_{1}$.
Since varying the twist angle changes the magnitude $Q^{2}$ that appears in
the denominator of $C$, we conclude that as long as $t_{\perp,Q}\neq 0$, there
will exist some magic twist angle where $C=v_{1}$. Importantly, this magic
twist angle requires no fine-tuning.

On the other hand, as mentioned below Eq. (9), $t_{\perp,Q}$ corresponds to a
spin-flipping interlayer hopping term. While this term is symmetry-allowed, it
is not obvious how such a term would naturally arise, and we might expect it
to be small. Since the magic angle occurs when $C=v_{1}$, which implies
$t_{\perp,Q}\sim vQ$, then if $t_{\perp,Q}$ is small, the magic angle requires
that $Q$ also be small, which corresponds to a small twist angle (see Fig.
1(a)). At such small twist angles, there are two places where our perturbative
results may fail. First: if spin-preserving couplings dominate, then we
require small $|Q|$ for $t_{\perp}/v|Q|$ to be sufficiently large. However,
since our expansion in $t$ includes a leading order approximation in
$t_{0}/v|Q|$ as well, the smallness of $|Q|$ may demand working to a higher
order in perturbation theory in $t_{0}$. Second: physically, our assumption of
rigid layers may break down at small twist angle if disorder or lattice
relaxation effects dominate.

In the absence of a layer-interchanging symmetry, generically $t_{||}\neq 0$.
Since this implies $D\neq 0$ whenever $C>0$, by Eq. (5), lowest-order
perturbation theory indicates magic angles will not arise. To higher order in
perturbation theory, fine-tuned scenarios may arise where a magic angle
exists, but we do not expect it to be generic.

In Sec. III.3 we will compare our perturbative predictions for the magic angle
to a global band structure calculation.

(a) $H=vk\cdot\sigma$

(b) $H=vk\times\sigma$

(c) $H=vk\cdot\sigma^{*}$

(d) $H=vk\times\sigma^{*}$

(e) $H=-vk\cdot\sigma$

(f) $H=-vk\times\sigma$

(g) $H=-vk\cdot\sigma^{*}$

(h) $H=-vk\times\sigma^{*}$

Figure 3: Spin textures of the conduction bands of different Dirac cones. For
Dirac cones where $\sigma$ appears in the Hamiltonian, spin and momentum
rotate the same way, whereas for Dirac cones where $\sigma^{*}$ appears, spin
and momentum rotate in opposite directions. This winding number also results
in different (unitarily equivalent) angular momentum representations.

###  III.2 Spin vs. pseudospin

The low-energy Hamiltonian $v(\bar{k}\cdot\sigma)$ is not the unique
description of a Dirac cone invariant under time-reversal and rotation
symmetries: specifically, as derived in Appendix A.1, any linear combination
of $\bar{k}\cdot\sigma$ and $\bar{k}\times\sigma$ is invariant. Such linear
combinations can be brought into the form $v(\bar{k}\cdot\sigma)$ by a unitary
transformation, although the transformation also changes the meaning of the
spin basis. As a first example, if the Hamiltonian in one layer is described
by $v(\bar{k}\cdot\sigma)$ and the Hamiltonian in the other layer is described
by $v(\bar{k}\times\sigma)$ or $-v(\bar{k}\cdot\sigma)$, then we reach the
same conclusion as in the previous section that magic angles exist, but they
require a different mixture of $t_{\perp}$ and $t_{||}$. Importantly, spin-
flipping hopping is still required.

However, we get a different result if one layer has the low-energy Hamiltonian
$-\bar{k}\cdot\sigma^{*}$. The unitary transformation that transforms this
Hamiltonian to the form $\bar{k}\cdot\tilde{\sigma}$ exchanges spin-flipping
and spin-preserving interactions (following the procedure described in
Appendix A.1). Thus, the interlayer hopping term that results in a magic angle
to leading order in perturbation theory is now a spin-preserving term. For
other $\sigma^{*}$-type Hamiltonians, the hopping required to achieve a magic
angle preserves the projection of spin in the $z$-direction, but not in all
directions, i.e., the requirement to achieve a magic angle is a linear
combination of $t_{0}$ and $t_{z}$.

From a physical perspective, this set-up provides a physically reasonable
(spin-preserving) interlayer hopping term to achieve flat bands at the
interface between two topological insulators. However, this physical
motivation for the interlayer hopping term comes at a cost: it requires two
different topological insulators, one with effective Hamiltonian
$\bar{k}\cdot\sigma$ and the other with $-\bar{k}\cdot\sigma^{*}$.

Physically, materials with a $\bar{k}\cdot\sigma^{*}$ Hamiltonian can be
distinguished from materials with a $\bar{k}\cdot\sigma$ Hamiltonian by their
spin texture: their spins rotate in opposite directions around the Dirac cone,
as illustrated in Fig. 3. Specifically, writing the linearized Hamiltonian for
layer $L$ as $k\cdot M_{L}\cdot\sigma$, the spin-winding of the Dirac cone is
classified by the parameter:

$\chi_{L}=\text{sgn}(\det(M_{L}))=\pm 1,$ (11)

where $\chi_{L}(k\cdot\sigma)=1$ and $\chi_{L}(-k\cdot\sigma^{*})=-1$. To
determine whether spin-flipping inter-layer hopping is necessary to realize a
magic angle, the relevant quantity is $\chi:=\chi_{1}\chi_{2}$, the product of
the $\chi$ of each layer: if the product is $+1$, $z$-spin-preserving
interlayer hopping is required, whereas if the product is $-1$, then $z$-spin-
flipping interlayer hopping is required.

Note that $\chi_{L}$ does not completely classify the Hamiltonian, and as a
consequence, the mixture of $t_{||}$ and $t_{\perp}$ required to achieve a
magic angle depends on the low-energy Hamiltonian. For example, the ratio of
$t_{\parallel}$ to $t_{\perp}$ required for a magic angle at the interface
between two 3D TIs with $k\cdot\sigma$ Dirac cones will be different than that
required at the interface between two 3D TIs with $k\times\sigma$ Dirac cones,
even though both cases have $\chi_{1}=\chi_{2}=+1$.

###  III.3 Global band structure

So far we have derived the existence/absence of a magic angle using the
renormalized Dirac cone velocity derived in Eq. (5) to leading order. We now
study the validity of our conclusions to higher order in both $k$ and $t_{0}$
by computing the full band structure. We assume spin-preserving coupling and
consider two cases: one in which both TIs have effective Hamiltonians of the
form $\bar{k}\cdot\sigma$ (as considered in Sec. III.1), and one in which the
second TI instead has an effective Hamiltonian of the form
$-\bar{k}\cdot\sigma^{*}$ (as discussed in Sec. III.2), where $\sigma$ always
references the spin degree of freedom, rather than a rotated pseudospin degree
of freedom. We will call these the $\chi=+1$ and $\chi=-1$ cases. Note that
the band structure for the $\chi=-1$ case with spin-preserving hopping is
identical to the $\chi=+1$ case with spin-flipping hopping, i.e.,
$t_{0}\rightarrow t_{\perp}$.

Figure 4: Moiré band structures for the interface between two 3D TIs. The path
ABB’A refers to the moiré BZ shown in Fig. 2. In the top row, the low-energy
Hamiltonian for both materials is $vk\cdot\sigma$ ($\chi=+1$), while in the
bottom row, the low-energy Hamiltonian for one TI is $vk\cdot\sigma$ and for
the other $-vk\cdot\sigma^{*}$ ($\chi=-1$). In all cases, only a single spin-
preserving interlayer coupling term is included. The band structures are
identical to those that result from a spin-flipping term by making the
substitutions $t_{0}\rightarrow t_{\perp}$ and $\chi\rightarrow-\chi$.

Fig. 4 shows the spectrum plotted along the ABB’A slice of the moiré BZ, which
is labelled in Fig. 2. The top row of plots are spectra for the $\chi=+1$
case, and the bottom row of plots are spectra for the $\chi=-1$ case, with
varying $t_{0}/v|Q|$ in different columns.

We expect in the $\chi=-1$ case (from Eqs (9c) and (5), taking
$t_{0}\rightarrow t_{\perp}$ according to the basis change described in the
previous section and derived in Appendix A.1) that a magic angle will arise
for $t_{0}/v|Q|=0.5$. Examining the $t_{0}/v|Q|=0.5$ plot, we find bands that
are appreciably flatter, but not yet “magic.” However, we observe that at
higher coupling, $t_{0}/v|Q|=0.75$ (Fig. 4, bottom-right) produces appreciably
flatter bands near a magic angle. We attribute this mismatch between our
perturbative and numerical results to the fact that at the predicted magic
angle, our perturbative parameter $t_{0}/v|Q|=0.5$ is not very small. In Fig.
5, we compare the numerical and perturbative calculations of the Fermi
velocity and see that they begin to deviate around $t_{0}/v|Q|\sim 0.2$.

Nonetheless, the band structures plotted in Fig. 4 confirm the qualitative
intuition (discussed in Sec. III.2) that, while both the $\chi=+1$ and
$\chi=-1$ cases produce some band flattening (compare the leftmost column, in
which the layers are uncoupled, to the other columns), the $\chi=-1$ case
produces much flatter bands with spin-preserving interactions than the
$\chi=+1$ case. This can be seen by comparing the spectra in the two rows of
Fig. 4: the bandwidth along the AB line, and in particular at D, is
consistently narrower in the $\chi=-1$ case.

The same trend appears in Fig. 5, where the Fermi velocity is plotted against
$t_{0}/v|Q|$. The two cases are notably closer in Fermi velocity than first-
order perturbation theory would predict for intermediate values of
$t_{0}/v|Q|$ (between $0.2$ and $0.5$), but the $\chi=-1$ case nevertheless
consistently yields lower Fermi velocities, and for large $t_{0}/v|Q|\gtrsim
0.6$, the $\chi=-1$ Fermi velocity becomes smaller by an order of magnitude
than the $\chi=+1$ Fermi velocity. In addition, while the Fermi velocity
oscillates in amplitude in the $\chi=+1$ case, the Fermi velocity in the
$\chi=-1$ case decays monotonically, well past the validity of the
perturbation theory.

Figure 5: Fermi velocity vs. $t_{0}$ for spin-preserving couplings in the
$\chi=+1$ (black) and $\chi=-1$ (red) cases. Dashed lines are first-order
perturbation theory, solid lines are extracted from numerically-computed band
structures. The $\chi=-1$ case consistently produces lower Fermi velocities.

In Fig. 6, we illustrate the extremely flat low-energy bands that result in
the $\chi=-1$ case with $t_{0}/v|Q|=1.3$. The figure shows that not only do
the lowest energy bands become flat, but the adjacent bands also collapse onto
the flat bands, creating a larger density of states than would result from the
flat bands alone. Both the accumulation of flat bands and the exponential
suppression of Fermi velocity were also discussed for twisted square lattices
in Ref. 47.

Figure 6: Spectrum for spin-preserving coupling in the $\chi=-1$ case with
$t_{0}/v|Q|=1.3$. The original Dirac cone becomes extremely flat and nearby
bands collapse in energy, creating a large density of states.

Returning to the spectra in Fig. 4, it is also notable that the spectra are
consistently gapless at all energies. This is not required by the strong
$\mathbb{Z}_{2}$ topological invariant: time-reversal symmetry allows gaps to
open at the interface between two 3D TIs. Instead, the gapless surface states
are protected as a weak TI: the cones at $A$ and $B$ are separated in the
moiré BZ (as illustrated in Fig. 2), and consequently are protected by the
approximate translation symmetry of the moiré lattice [48], which also
protects the cones at $C$ and $C^{\prime}$. By “moiré BZ” we refer to the BZ
defined below Eq. (3), not the BZ of any particular commensurate lattice;
hence, the translation symmetries are not sensitive to any particular
commensurate structure. However, the protection of Dirac cones by moiré
translation symmetry implies that disorder is sufficiently small on the moiré
length scale. (The simplicity of this argument is in contrast to the analogous
circumstances in the second chiral limit of TBLG [49] and in twisted square
lattices [47] in which a more nuanced argument is required.)

Dirac cones also appear at $D$ and $D^{\prime}$ in the top row in Fig. 4.
These crossings are not topologically protected: there exist non-spin-
conserving interlayer couplings that are symmetry-allowed and would open a gap
at $D$. In particular, in the lower row in Fig. 4, the Dirac cones have the
opposite spin-momentum locking; hence, the spins are aligned where the Dirac
cones overlap at $D$ and a large gap opens.

Instead, these crossings’ gaplessness is a consequence of the specific
parameters used in our model. In particular: since $D$ lies on the $A-B$ line
(see Fig. 2), when the Dirac cones have the same spin-momentum locking (i.e.,
both are of the form $\bar{k}\cdot\sigma$), they have opposite spin at $D$.
Since the spectra in Fig. 4 only include spin-conserving hopping, there is no
term to open a gap at $D$. When the two Dirac cones have different velocities
or chemical potentials, the crossing still remains gapless, but moves along
the $A-B$ line.

Figure 7: Band structures for the interface between two identical twisted 3D
TIs with spin-preserving coupling (BZ path defined in Fig. 2). At
$t_{0}/v|Q|=0.75$, the bands nearest the Fermi level become nearly degenerate
along AB.

There is one additional interesting feature to examine. In the $\chi=+1$ case,
there are two nearly-degenerate bands at $t_{0}/v|Q|=0.75$. Fig. 7 shows a
finer range of hopping parameters, which reveals the lowest conduction bands
descend to produce additional zero-energy states at a critical value of
$t_{0}/v|Q|=0.75$, then separate again as the coupling increases. A similar
critical point where the lowest conduction band dips down and touches the flat
band also appears in TBLG, both in the usual Bistritzer-MacDonald model (Fig.
3A of [5]) and in the “second chiral model” with vanishing AB coupling (Fig.
2c of [22]), although the nearly-degenerate bands that result at this crossing
are specific to our model.

The touching point appears to be distinct in the two cases: in the TI-TI case,
the touching point is at the Dirac cones, whereas in TBLG the touching point
is at $\Gamma$ (away from the Dirac cones). However, this can be understood by
interpreting the touching point as occurring $-Q$ away from a Dirac point,
which in TBLG is the $\Gamma$ point and in our TI-TI model is at a Dirac cone.

##  IV TI on Graphene

We have shown that stacking two 3D TIs with a small twist angle can
significantly renormalize the Dirac cones at the interface, with nearly flat
bands appearing at a specific magic twist angle under certain conditions. In
the simplest case, the surface Dirac cone of each TI must not be at the center
of the surface BZ. However, none of the 3D TIs in the Bi2Se3 family meet this
condition: their surface Dirac cones are always at $\Gamma$ [50, 51].

This motivates us to consider a second example of a twisted heterostructure
with multiple Dirac cones: the interface between a topological insulator and
graphene. Such interfaces have been considered in previous work [23, 24, 25,
26, 27, 28, 29, 30, 31, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 32],
but the effect of a small relative twist angle was not considered.

We consider a 3D TI with a single Dirac cone at $\Gamma$, such as in Bi2Se3.
The lattice constants of the TI and graphene differ by a factor of $\sqrt{3}$
up to an approximate $3\%$ lattice mismatch[26]; Fig. 8 shows the lattice
matching ignoring the lattice mismatch. The aligned superstructure then folds
the $K$ and $K^{\prime}$ points in graphene onto $\Gamma$ in the TI BZ, with
the resulting twisted structure shown in Fig. 9. We assume the TI has a
sixfold rotational symmetry, as well as time-reversal symmetry. (The sixfold
symmetry is not exact: while each layer of Bi2Se3 has a sixfold rotation axis,
the layers are stacked such that the axes do not coincide [50]. Thus, surface
states mostly localized on the top layer will exhibit an approximate sixfold
rotational symmetry.) The TI Hamiltonian is given by $v_{TI}k\cdot\sigma$, as
in Sec. III.1. Graphene’s Hamiltonian consists of four Dirac cones including
spin and valley; the matrix form is given in Appendix A.3.

Figure 8: Graphene on Bi2Se3. Blue circles indicate Se atoms on the Bi2Se3
surface; red circles indicate carbon atoms in graphene. The superlattice
matching is accurate to $\sim 3$% [26].

Figure 9: The BZ of graphene (red) arranged on top of the 3D TI BZ (blue),
shown both with the lattice mismatch (top) and then with a small relative
twist angle while neglecting mismatch (bottom). The momentum differences
$Q_{i}$ are indicated in each case.

In the subsequent sections, we will derive a perturbative analytic calculation
of the renormalization of the TI Dirac cone that arises from its coupling to
graphene. The analogous calculation for graphene is more complicated due to
its spin degeneracy and is derived in Appendix C.4; when the interlayer
hopping is spin-preserving, almost all corrections vanish, with the
nonvanishing coefficients only rescaling the spectrum.

To go beyond our perturbative results, we compute the surface band structure,
specifying to the most physically interesting case where spin is conserved.
Our results both validate the first-order perturbative calculation for the TI
Dirac cone renormalization and provide insight into the higher-order
corrections to the graphene Dirac cone.

###  IV.1 Corrections to TI Dirac cone

We begin by computing the corrections to the TI Dirac cone to first order in
perturbation theory. In so doing, we make use of the corrections to a single
Dirac cone derived in Eqs. (4)-(8). We decompose the interlayer coupling terms
$T_{Q}$, defined in Eq. (3), into $T_{Q,k_{0},s}$, where $T_{Q,k_{0},s}$
couples each graphene Dirac cone to the TI Dirac cone, with the $k_{0}$
running over valley and $s$ over spin. At any particular $Q$, only one valley
will couple to the TI; without loss of generality, we consider $Q$
corresponding to the $K$ point, so that $T_{Q,K^{\prime}}=0$. The interlayer
coupling terms for other values of $Q$ can then be determined by time-reversal
symmetry.

When the lattices are aligned, there is a moiré pattern that arises from the
lattice mismatch. In this case, we can choose an orientation such that
$Q=Q_{x}\hat{x}$. On the other hand, if we ignore the lattice mismatch and
consider arranging the graphene layer on the TI surface with a small twist
angle, we can take $Q=Q_{y}\hat{y}$. In the most general case of a lattice
mismatch and a small twist angle, $Q$ could point in an arbitrary direction;
our formalism covers this case, but we here focus on the two limits.

Regardless of the combination of twist and mismatch, for a spin-preserving
coupling, the corresponding $T_{Q,K,s}$ will take one of two forms (depending
on $s$), with the rows corresponding to spins in the TI layer and the columns
corresponding to sublattice in the graphene layer:

$\displaystyle T_{Q,K,\uparrow}=\begin{bmatrix}t&t\\\ 0&0\end{bmatrix}$ (12a)
or $\displaystyle T_{Q,K,\downarrow}=\begin{bmatrix}0&0\\\ t&t\end{bmatrix}$
(12b)

For either of these $T_{Q,K}$, $t_{||}=\pm it_{\perp}$, where $t_{||}$ and
$t_{\perp}$ are defined in Eq. (7). Therefore, the direct corrections $C$ and
$D$ to the Fermi velocity (computed according to Eq. (8) and inputted into the
Fermi velocity in Eq. (5)) both vanish. Consequently, magic angles from spin-
preserving hoppings do not exist for any small twist or strain.

By the same logic, purely $z$-spin-flipping hoppings will also not produce
magic angles in this case. Unlike the interface between two 3D TIs (Sec. III),
we need both $z$-spin-preserving and $z$-spin-flipping hoppings to produce
magic angles.

Another change compared to the interface between two 3D TIs is that here, due
to the small lattice mismatch, the coupling to graphene shifts the TI’s
surface Dirac cone up or down in energy, while the Dirac cones at the
interface between two 3D TIs experience no such shift when the interlayer
coupling is spin-preserving. (The lack of an energy shift in that case is a
consequence of the parameter $B$ in Eq. (9b) vanishing for spin-preserving
coupling; here, it vanishes for a small twist, but not for lattice mismatch.)
In this case, the energy shift follows the formula:

$\Delta E=-\frac{12t^{2}v_{G}|Q|}{v_{G}^{2}Q^{2}+12t^{2}}$ (13)

The sign of the energy shift depends on the direction of mismatch, i.e.,
whether the graphene unit cell or the TI unit cell is larger compared to the
perfect supercell matching illustrated in Fig. 8.

###  IV.2 Global band structure

We now plot the band structure to see how our perturbation theory holds up. We
plot all spectra along the path shown in Fig. 10, which is the same path used
in [5].

Figure 10: Path along which spectra in Fig. 11 are plotted in the case of
lattice mismatch. For a twist, the figure is rotated 90 degrees
counterclockwise to account for the different values for $Q$.

We first examine spin-preserving, sublattice-independent coupling (which does
not produce magic angles). We present one representative spectrum each for the
cases of twist and lattice mismatch in Fig. 11, computed with
$t/v_{TI}|Q|=0.25$; other values of $t$ look qualitatively similar. The band
structure near the TI Dirac cone at A agrees with our perturbative
calculations: in the case of a lattice mismatch, it shows an energy shift
whose magnitude and direction matches Eq. (13); there is no energy shift in
the twist case; and in either case there is only modest velocity
renormalization. Fig. 12 illustrates the behavior of Fermi velocity and energy
shift, as extracted from band structure, as a function of twist angle. It
shows that the Fermi velocity does not experience nearly the renormalization
achieved at the interface between two 3D TIs (Fig. 5). It also shows that at
small twist angles and lattice mismatches (large couplings), twisting
suppresses Fermi velocity more than mismatch for the same moiré lattice size.

Figure 11: Spectrum of TI-graphene heterostructure with zero interlayer
hopping (left) and spin-preserving sublattice-independent interlayer hopping
$t/v_{TI}|Q|=0.25$ for a uniform lattice mismatch without twist (center) and
twist without mismatch (right). We assume $v_{G}=2v_{TI}$ [23].

Figure 12: Renormalized TI Dirac cone velocity $v_{TI}^{*}/v_{TI}$ (top) and
energy shift $E_{0}/v_{TI}|Q|$ (bottom) as a function of the spin-preserving
and sublattice-independent hopping parameter $t/v_{TI}|Q|$ extracted from band
structure calculations (solid lines) and computed perturbatively (dashed line,
drawn only in the one case where it does not overlap with a solid line), for
the case of a small twist angle (red) and small lattice mismatch (black). In
the twist case, the energy and Fermi velocity perfectly match the perturbative
calculations, so those lines are not drawn. In the lattice mismatch case, the
Fermi velocity has the same perturbative theory as the twist case (i.e.,
following the solid red line).

The graphene cones at B, C, and D resemble the quadratic bands of untwisted
Bernal-stacked bilayer graphene. This is because the TI couples graphene’s
spin-up and spin-down cones at $K$ (since propagation in the TI layer flips
the $z$-component of spin). These cones are not twisted relative to each
other, so they couple in much the same way as the Dirac cones in the different
layers of bilayer graphene, but with an effective interlayer hopping parameter
proportional to $\frac{t^{2}}{v_{TI}|Q|}$ (instead of $t_{BLG}$).

This correction to the graphene Dirac cones is fourth order in $t$, which
agrees with our analytic calculation showing that the only nonvanishing
correction to the graphene cone to quadratic order in $t$ is an overall
rescaling of the spectrum (details in Appendix C.4).

Like the TI-TI interface, the spectra in Fig. 11 remain gapless at all
energies. This is topologically required because the Dirac cone on the 3D TI
surface is protected by time reversal symmetry; unlike the TI-TI case
discussed in Sec. III.3, here the topological protection is strong and not
weak (i.e., not dependent on the approximate translation symmetries of the
moiré lattice) since there is only one TI surface state.

In summary, our calculations predict that arranging graphene on the surface of
a 3D TI with a small twist angle is not a promising platform for creating
flatter bands, both because no gaps are topologically allowed to open in the
spectrum and because under the most physically intuitive conditions, the Dirac
cone is only marginally renormalized.

##  V Conclusions

We have studied moiré heterostructures on the surface of a 3D TI. We derived
analytic expressions for the leading order corrections to the velocity of the
surface Dirac cone induced by coupling to another lattice-matched Dirac
material arranged with a small twist angle. We applied our results to two
types of heterostructures: an interface between two TIs and graphene on the
surface of a TI. We derived conditions for achieving a “magic” angle at which
the Fermi velocity vanishes.

One of our main results is that at the interface between two identical 3D TIs
arranged with a small twist angle, a spin-flipping interlayer hopping term is
a necessary ingredient to achieve a magic angle. The same is true for graphene
on a 3D TI. While such a term is symmetry-allowed (enforcing time-reversal and
rotational symmetry), it is not clear what physical mechanism would give rise
to it. In future work, it would be interesting to perform first-principles
calculations for different material combinations (such as graphene on Bi2Se3
or Bi2Te3, or the interface between two 3D TIs) to determine under what
conditions the spin-flipping interlayer hopping terms arise.

On the other hand, we found that magic angles are achieved without spin-
flipping interlayer hopping at twisted interfaces between two lattice-matched
TIs whose Dirac cones have opposite winding number (see Fig. 3). At the magic
angle, such an interface realizes locally flat bands. Unlike in twisted
bilayer graphene, the flat bands are not gapped from the conduction/valence
bands due to topological constraints. Nonetheless, such locally flat bands
will enhance the density of states at the Dirac point, which is favorable for
realizing instabilities to superconducting or quantum anomalous Hall states
[12, 13, 14, 15, 16, 17, 18, 19, 10].

This analysis therefore provides two routes for inducing magic angles on the
surface of a 3D TI. There are also several potential routes to creating flat
bands in such systems that go beyond the setups considered in this paper.

One route would be to consider an interface between two materials with Dirac
cones at $\Gamma$. We did not consider this case here because such interfaces
would generically gap, and therefore exhibit qualitatively different behavior
than other twisted interfaces, although this case has been studied recently
for thin slabs [52]. Whether the resulting bands contain interesting features
(e.g., van Hove singularities [10]) remains an open question.

Another route would be to reduce the rotational symmetry. We limited ourselves
to Dirac cones that would remain isotropic to linear order in the twisted
heterostructure, i.e., those that remain a center of three- or four-fold
rotation. If one considers Dirac cones with less rotational symmetry, then
even where a vanishing Fermi velocity is impossible, it may be possible to
achieve a significantly enhanced density of states via anisotropic band
flattening [45, 46]. In such a setup, one may have “partial magic angles”
where Dirac cone velocity vanishes in one direction. The interacting
instabilities of a 3D TI with such a Dirac cone provide a potential scenario
for further study.

##  VI Acknowledgments

We are grateful for useful conversations with Shiang Fang, Jed Pixley, Jie
Wang, and Justin Wilson and for a careful reading of the manuscript by Sayed
Ali Akbar Ghorashi. This work was partially supported by the Air Force Office
of Scientific Research under Grant No. FA9550-20-1-0260. J.C. acknowledges the
support of the Flatiron Institute, a division of the Simons Foundation.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/content/108/30/12233.full.pdf .
  * Lopes dos Santos  _et al._ [2012a] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012a).
  * Suárez Morell  _et al._ [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, Phys. Rev. B 82, 121407(R) (2010).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: topological insulators, Reviews of modern physics 82, 3045 (2010).
  * Cano  _et al._ [2021] J. Cano, S. Fang, J. H. Pixley, and J. H. Wilson, Moiré superlattice on the surface of a topological insulator, Phys. Rev. B 103, 155157 (2021).
  * Wang  _et al._ [2021] T. Wang, N. F. Q. Yuan, and L. Fu, Moiré surface states and enhanced superconductivity in topological insulators, Phys. Rev. X 11, 021024 (2021).
  * Guerci  _et al._ [2022] D. Guerci, J. Wang, J. Pixley, and J. Cano, Designer meron lattice on the surface of a topological insulator, arXiv preprint arXiv:2203.04986 (2022).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Superconducting proximity effect and majorana fermions at the surface of a topological insulator, Phys. Rev. Lett. 100, 096407 (2008).
  * Santos  _et al._ [2010] L. Santos, T. Neupert, C. Chamon, and C. Mudry, Superconductivity on the surface of topological insulators and in two-dimensional noncentrosymmetric materials, Phys. Rev. B 81, 184502 (2010).
  * Baum and Stern [2012a] Y. Baum and A. Stern, Magnetic instability on the surface of topological insulators, Phys. Rev. B 85, 121105(R) (2012a).
  * Baum and Stern [2012b] Y. Baum and A. Stern, Density-waves instability and a skyrmion lattice on the surface of strong topological insulators, Phys. Rev. B 86, 195116 (2012b).
  * Marchand and Franz [2012] D. J. J. Marchand and M. Franz, Lattice model for the surface states of a topological insulator with applications to magnetic and exciton instabilities, Phys. Rev. B 86, 155146 (2012).
  * Schmidt [2012] M. J. Schmidt, Strong correlations at topological insulator surfaces and the breakdown of the bulk-boundary correspondence, Phys. Rev. B 86, 161110 (2012).
  * Sitte  _et al._ [2013] M. Sitte, A. Rosch, and L. Fritz, Interaction effects on almost flat surface bands in topological insulators, Phys. Rev. B 88, 205107 (2013).
  * Mendler  _et al._ [2015] D. Mendler, P. Kotetes, and G. Schön, Magnetic order on a topological insulator surface with warping and proximity-induced superconductivity, Phys. Rev. B 91, 155405 (2015).
  * Liu  _et al._ [2021] Z. Liu, H. Wang, and J. Wang, Magnetic moiré surface states and flat chern band in topological insulators, arXiv preprint arXiv:2106.01630 (2021).
  * Chou  _et al._ [2021] Y.-Z. Chou, J. Cano, and J. Pixley, Band manipulation and spin texture in interacting moiré helical edges, arXiv preprint arXiv:2106.05270 (2021).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. III. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021).
  * Rossi and Triola [2019] E. Rossi and C. Triola, Van der waals heterostructures with spin‐orbit coupling, Annalen der Physik 532, 1900344 (2019).
  * Jin  _et al._ [2012] H. Jin, J. Im, J.-H. Song, and A. J. Freeman, Multiple Dirac fermions from a topological insulator and graphene superlattice, Phys. Rev. B 85, 045307 (2012).
  * Jin and Jhi [2013] K.-H. Jin and S.-H. Jhi, Proximity-induced giant spin-orbit interaction in epitaxial graphene on a topological insulator, Phys. Rev. B 87, 075442 (2013).
  * Zhang  _et al._ [2014] J. Zhang, C. Triola, and E. Rossi, Proximity effect in graphene–topological-insulator heterostructures, Phys. Rev. Lett. 112, 096802 (2014).
  * Cao  _et al._ [2016] W. Cao, R.-X. Zhang, P. Tang, G. Yang, J. Sofo, W. Duan, and C.-X. Liu, Heavy Dirac fermions in a graphene/topological insulator hetero-junction, 2D Materials 3, 034006 (2016).
  * De Beule  _et al._ [2017] C. De Beule, M. Zarenia, and B. Partoens, Transmission in graphene–topological insulator heterostructures, Phys. Rev. B 95, 115424 (2017).
  * Rodriguez-Vega  _et al._ [2017] M. Rodriguez-Vega, G. Schwiete, J. Sinova, and E. Rossi, Giant edelstein effect in topological-insulator–graphene heterostructures, Physical Review B 96, 235419 (2017).
  * Song  _et al._ [2018] K. Song, D. Soriano, A. W. Cummings, R. Robles, P. Ordejón, and S. Roche, Spin proximity effects in graphene/topological insulator heterostructures, Nano letters 18, 2033 (2018).
  * Zollner and Fabian [2020] K. Zollner and J. Fabian, Heterostructures of graphene and topological insulators Bi2Se3, Bi2Te3, and Sb2Te3, physica status solidi (b) 258 (2020).
  * Zollner and Fabian [2019] K. Zollner and J. Fabian, Single and bilayer graphene on the topological insulator Bi2Se3: Electronic and spin-orbit properties from first principles, Phys. Rev. B 100, 165141 (2019).
  * Dang  _et al._ [2010] W. Dang, H. Peng, H. Li, P. Wang, and Z. Liu, Epitaxial heterostructures of ultrathin topological insulator nanoplate and graphene, Nano letters 10, 2870 (2010).
  * Song  _et al._ [2010] C.-L. Song, Y.-L. Wang, Y.-P. Jiang, Y. Zhang, C.-Z. Chang, L. Wang, K. He, X. Chen, J.-F. Jia, Y. Wang, _et al._ , Topological insulator Bi2Se3 thin films grown on double-layer graphene by molecular beam epitaxy, Applied Physics Letters 97, 143118 (2010).
  * Lee  _et al._ [2015] P. Lee, K.-H. Jin, S. J. Sung, J. G. Kim, M.-T. Ryu, H.-M. Park, S.-H. Jhi, N. Kim, Y. Kim, S. U. Yu, _et al._ , Proximity effect induced electronic properties of graphene on Bi2Te2Se, Acs Nano 9, 10861 (2015).
  * Steinberg  _et al._ [2015] H. Steinberg, L. A. Orona, V. Fatemi, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunneling in graphene–topological insulator hybrid devices, Physical Review B 92, 241409(R) (2015).
  * Zhang  _et al._ [2016] L. Zhang, Y. Yan, H.-C. Wu, D. Yu, and Z.-M. Liao, Gate-tunable tunneling resistance in graphene/topological insulator vertical junctions, ACS nano 10, 3816 (2016).
  * Bian  _et al._ [2016] G. Bian, T.-F. Chung, C. Chen, C. Liu, T.-R. Chang, T. Wu, I. Belopolski, H. Zheng, S.-Y. Xu, D. S. Sanchez, N. Alidoust, J. Pierce, B. Quiliams, P. P. Barletta, S. Lorcy, J. Avila, G. Chang, H. Lin, H.-T. Jeng, M.-C. Asensio, Y. P. Chen, and M. Z. Hasan, Experimental observation of two massless Dirac-fermion gases in graphene-topological insulator heterostructure, 2D Materials 3, 021009 (2016).
  * Vaklinova  _et al._ [2016] K. Vaklinova, A. Hoyer, M. Burghard, and K. Kern, Current-induced spin polarization in topological insulator–graphene heterostructures, Nano letters 16, 2595 (2016).
  * Rajput  _et al._ [2016] S. Rajput, Y.-Y. Li, M. Weinert, and L. Li, Indirect interlayer bonding in graphene–topological insulator van der Waals heterostructure: Giant spin–orbit splitting of the graphene Dirac states, ACS nano 10, 8450 (2016).
  * Zhang  _et al._ [2017] L. Zhang, B.-C. Lin, Y.-F. Wu, H.-C. Wu, T.-W. Huang, C.-R. Chang, X. Ke, M. Kurttepeli, G. V. Tendeloo, J. Xu, _et al._ , Electronic coupling between graphene and topological insulator induced anomalous magnetotransport properties, ACS nano 11, 6277 (2017).
  * Zalic  _et al._ [2017] A. Zalic, T. Dvir, and H. Steinberg, High-density carriers at a strongly coupled interface between graphene and a three-dimensional topological insulator, Physical Review B 96, 10.1103/physrevb.96.075104 (2017).
  * Khokhriakov  _et al._ [2018] D. Khokhriakov, A. W. Cummings, K. Song, M. Vila, B. Karpiak, A. Dankert, S. Roche, and S. P. Dash, Tailoring emergent spin phenomena in Dirac material heterostructures, Science advances 4, eaat9349 (2018).
  * Chae  _et al._ [2019] J. Chae, S.-H. Kang, S. H. Park, H. Park, K. Jeong, T. H. Kim, S.-B. Hong, K. S. Kim, Y.-K. Kwon, J. W. Kim, _et al._ , Closing the surface bandgap in thin Bi2Se3/graphene heterostructures, ACS nano 13, 3931 (2019).
  * Kariyado and Vishwanath [2019] T. Kariyado and A. Vishwanath, Flat band in twisted bilayer bravais lattices, Phys. Rev. Research 1, 033076 (2019).
  * Kennes  _et al._ [2020] D. M. Kennes, L. Xian, M. Claassen, and A. Rubio, One-dimensional flat bands in twisted bilayer germanium selenide, Nature communications 11, 1 (2020).
  * Luo  _et al._ [2021] Z.-X. Luo, C. Xu, and C.-M. Jian, Magic continuum in a twisted bilayer square lattice with staggered flux, Phys. Rev. B 104, 035136 (2021).
  * Lopes dos Santos  _et al._ [2012b] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Physical Review B 86, 10.1103/physrevb.86.155449 (2012b).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. II. stable symmetry anomaly, Physical Review B 103, 10.1103/physrevb.103.205412 (2021).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface, Nature physics 5, 438 (2009).
  * Xia  _et al._ [2009] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, _et al._ , Observation of a large-gap topological-insulator class with a single Dirac cone on the surface, Nature physics 5, 398 (2009).
  * Tateishi and Hirayama [2021] I. Tateishi and M. Hirayama, Quantum spin hall effect from multi-scale band inversion in twisted bilayer bi2(te1-xsex)3 (2021).
  * Takagi [1924] T. Takagi, On an algebraic problem reluted to an analytic theorem of Carathéodory and Fejér and on an allied theorem of Landau, in  _Japanese journal of mathematics: transactions and abstracts_ , Vol. 1 (The Mathematical Society of Japan, 1924) pp. 83–93.
  * Autonne [1915] L. Autonne,  _Sur les matrices hypohermitiennes et sur les matrices unitaires_ , Vol. 38 (A. Rey, 1915).

##  Appendix A General Hamiltonians for Dirac Materials

In this appendix, we discuss the Dirac Hamiltonian in a single layer. In Sec.
A.1, we discuss our basis choice for the Dirac Hamiltonian. In Sec. A.2 we
consider the representations of rotational and time-reversal symmetry on this
Hamiltonian, and show that, after a further change of basis, time-reversal
symmetry can be brought to a canonical form of $i\sigma_{y}K$. Finally, in
Sec. A.3 we elaborate on the details and symmetry representation of the
graphene Hamiltonian, along with explicit matrices that perform the
aforementioned changes of basis.

###  A.1 Dirac Hamiltonians

Let us begin by considering a two-state Hamiltonian in two dimensions which
identically vanishes in energy at $k=0$. The most general expansion to linear
order in $k$ can be expressed as:

$H(k)=(a\cdot k)\mathbb{I}+k\cdot b\cdot\sigma+(c\cdot k)\sigma_{z}$ (14)

If the spectrum is isotropic, then the vectors $a$ and $c$ must vanish, while
the tensor $b$ is constrained to the trace and antisymmetric parts (which
commute with rotations). Then the Hamiltonian can be written as a sum of two
terms:

$H(k)=b_{1s}(\bar{k}\cdot\sigma)+b_{1a}(\bar{k}\times\sigma)$ (15)

We can then perform a unitary transformation described by the matrix:

$\exp(i\arctan(b_{1a}/b_{1s})\sigma_{z}/2)$ (16)

This will transform our Hamiltonian into:

$\sqrt{b_{1s}^{2}+b_{1a}^{2}}(\bar{k}\cdot\sigma)$ (17)

Thus, an isotropic Dirac cone at $k=0$ can generically be written as
$v\bar{k}\cdot\sigma$ in the appropriate basis.

A Dirac cone not at the origin will instead take the form
$v(k-k_{0})\cdot\sigma$. For this reason, we now define $\bar{k}$ as the
difference in momentum space between $k$ and the center of the Dirac cone
$k_{0}$:

$\bar{k}=k-k_{0}$ (18)

This leads to the more general assertion that any isotropic, gapless Dirac
cone at charge-neutrality can have an effective Hamiltonian written in the
form $H_{D}(\bar{k})$, with $H_{D}$ as defined in Eq. (1).

More generally, for a Hamiltonian with multiple Dirac cones indexed by the
Dirac points $k_{0,i}$, we can do the same operations for each cone in terms
of $\bar{k}_{i}:=k-k_{0,i}$.

At the level of the perturbative small-angle calculation performed here, we
can drop the $i$ sub-index, and regard $\bar{k}$ as a formal parameter rather
than a point in momentum space per se. This is computationally convenient
because it allows for valley-mixing basis changes (which we will use in the
following section); accordingly, we will use $\bar{k}$ without sub-index in
the rest of this paper’s exposition. (To extract the physical meaning of the
terms we compute, the basis change must be reversed to disentangle the
valleys.)

Therefore, under the assumptions we make in this paper, if we index our Dirac
cones by $i$ and $j$, we can always write a single-layer Hamiltonian as:

$H(\bar{k})=H_{D}(\bar{k})\delta_{ij}.$ (19)

###  A.2 Time-reversal and rotational symmetry

We now also assume that our materials have both time-reversal symmetry
$\mathcal{T}$ and in-plane $n$-fold rotational symmetry $C_{n}$. Note both of
these operations leave each layer invariant, so we can discuss them before
coupling two layers; however, in materials with multiple Dirac cones, the
symmetries do not have to leave each Dirac cone invariant - the Dirac cones
may be mixed under the symmetry transformations.

Time-reversal symmetry can be expressed as the condition that:

$\mathcal{T}H_{\bar{k}}\mathcal{T}^{-1}=H_{-\bar{k}}$ (20)

Similarly, if we denote by $\mathfrak{R}$ the usual action of $C_{n}$ on the
plane, then the rotation symmetry requires:

$C_{n}H_{\bar{k}}C_{n}^{\dagger}=H_{\mathfrak{R}\bar{k}}$ (21)

We now derive the matrix forms of $\mathcal{T}$ and $C_{n}$, which will be
implemented by anti-unitary and unitary matrices, respectively, when acting on
our Dirac Hamiltonian in Eq. (19). Given that $\mathcal{T}$ is antiunitary, it
can generally be written as:

$\mathcal{T}=\tilde{T}_{ij}(i\sigma_{y}\mathcal{K}),$ (22)

where $\tilde{T}$ is a matrix in $ij$ space (trivial in $\sigma$ space) whose
structure will depend on the specific problem at hand. Since $\mathcal{T}$ is
anti-unitary and satisfies $\mathcal{T}^{2}=-1$, it must be that
$\tilde{T}\tilde{T}^{\dagger}=\tilde{T}\tilde{T}^{*}=1$, and therefore that
$\tilde{T}$ is complex-symmetric.

The rotation matrix can be written as:

$C_{n}=\tilde{U}_{ij}\exp(i\pi\sigma_{z}/n),$ (23)

Similar to time-reversal symmetry, there is a matrix $\tilde{U}$ which acts
only in $ij$ space and will depend on the specific problem. $\tilde{U}$ must
be unitary and have the property that $\tilde{U}^{n}=1$.

To simplify our calculations of twisted multi-Dirac band structures (Appendix
C), we would like to simplify the form of $\tilde{T}$ and $\tilde{U}$ as much
as possible. As these unitary transformations should also preserve the form of
our Hamiltonian, we constrain ourselves to only operating on the $i$ and $j$
indices (rather than on those of $H_{D}(\bar{k})$).

Since $\tilde{T}$ is complex-symmetric, we can perform an Autonne-Takagi
factorization [53, 54] to diagonalize it into $\delta_{ij}$. (Note the
condition that $\mathcal{T}^{2}=-1$ plays a vital role in this decomposition
by ensuring the symmetry of $\tilde{T}$; therefore, this decomposition does
not directly apply to spinless models.)

This change of basis yields new rotation matrices $U_{ij}$ from our previous
$\tilde{U}_{ij}$, such that the rotation in the new basis is
$U_{ij}\exp(i\pi\sigma_{z}/n)$. These matrices $U$ cannot generally be
simplified by a further basis transformation while preserving the properties
established thus far. (We still have the ability to act on the $ij$ indices
with a real orthogonal matrix without disturbing either our Hamiltonian or the
form of $\mathcal{T}$, but this is of comparatively limited use.)

In some cases, it is beneficial to diagonalize in the basis of
$\mathcal{T}^{\prime}:=C_{2}\mathcal{T}$, instead of $\mathcal{T}$. This will
most often the case if the model has a subset of Dirac cones (decoupled from
the rest) invariant under $\mathcal{T}^{\prime}$ but not $\mathcal{T}$, which
occurs, for instance, in twisted bilayer graphene, where the $K$ and
$K^{\prime}$ Dirac cones are decoupled from each other.

In this circumstance, $\mathcal{T}^{\prime}$ takes the matrix form:

$\mathcal{T}^{\prime}=\tilde{T}_{ij}^{\prime}(\sigma_{x}\mathcal{K})$ (24)

The condition that $\mathcal{T}^{\prime 2}=+1$ results in $T_{ij}^{\prime}$
also being complex-symmetric, and therefore the same Autonne-Takagi
factorization applies to reduce this to $\delta_{ij}$. However, this
diagonalization cannot (generally) be performed simultaneously with the
previous one, so one basis or the other should be chosen for the computation
on a case-by-case basis.

###  A.3 Effective graphene Hamiltonian and symmetry representations

We now present the application of the above methodology to graphene as a
concrete example. In the process, we develop a model of “half-graphene” which
will serve an effective simpler stepping stone in our calculations.

####  A.3.1 Graphene fundamentals

We can write an effective graphene Hamiltonian as four Dirac cones coupled to
sublattice, indexed by spin and valley:

$H_{G}=\begin{bmatrix}vk\cdot\tau&0&0&0\\\ 0&vk\cdot\tau&0&0\\\
0&0&-vk\cdot\tau^{*}&0\\\ 0&0&0&-vk\cdot\tau^{*}\end{bmatrix}$ (25)

We will here use $\tau$ to refer to the sublattice degree of freedom, $\sigma$
to refer to the spin degree of freedom, and $\mu$ to refer to the valley
degree of freedom. Therefore, the above Hamiltonian can be equivalently
written as:

$H_{G}=v[k_{x}\tau_{x}\mu_{z}+k_{y}\tau_{y}\mu_{0}]\sigma_{0}$ (26)

In this basis, time-reversal symmetry will be written as:

$\mathcal{T}=i\sigma_{y}\mu_{x}\mathcal{K}$ (27)

For rotations, it will be computationally convenient to separately consider
$C_{2}$ and $C_{3}$ (rather than $C_{6}$). These take the form:

$C_{2}=i\sigma_{z}\tau_{x}\mu_{x}$ (28)
$C_{3}=-\exp(i\pi\sigma_{z}/3)\exp(i\pi\tau_{z}\mu_{z}/3)$ (29)

####  A.3.2 Half-graphene model

A computationally simpler model would have only one Dirac cone each at $K$ and
$K^{\prime}$. Considering a spinless model of graphene will not work if we
want to preserve time-reversal symmetry with $\mathcal{T}^{2}=-1$. However,
what we can do is write a Hamiltonian for a material like graphene with spin-
valley locking. In other words, we take a material with Dirac cones coupled to
sublattice at $K$ and $K^{\prime}$, but take a single Dirac cone at $K$ with
one spin and one at $K^{\prime}$ with the opposite spin. This setup breaks
$C_{2}$ symmetry, but preserves both $C_{3}$ and time-reversal.

The effective Hamiltonian of such a model, before any basis change, could be
written by picking out the first and fourth cones of our previous graphene
model:

$H_{HG}=\begin{bmatrix}vk\cdot\tau&0\\\ 0&-vk\cdot\tau^{*}\end{bmatrix}$ (30)

We can therefore inherit our symmetry representations from that model as well.
Let $\sigma$ denote the combined spin-valley degree of freedom. Then, time
reversal takes the form:

$\mathcal{T}=i\sigma_{y}\mathcal{K}$ (31)

Similarly, threefold rotations take the form:

$C_{3}=-\exp(i\pi\sigma_{z}(\tau_{z}+\tau_{0})/3)$ (32)

####  A.3.3 Basis change in half-graphene

We now take our simplified model of half-graphene and put it into the standard
form with a Hamiltonian of $(vk\cdot\tau)\sigma_{0}$ with
$\mathcal{T}=i\sigma_{y}\mathcal{K}$.

First, we transform the Hamiltonian to our standard form with a $\sigma_{y}$
transformation on the second Dirac cone only. Time-reversal then takes the
form $(i\tau_{y}\mathcal{K})\sigma_{x}$. We can easily transform $\sigma_{x}$
into $\sigma_{z}$ with a rotation by $e^{i\pi\sigma_{y}/4}$, after which we
transform the second Dirac cone by a factor of $i$ to eliminate the minus sign
(exploiting antiunitarity to eliminate the cone-dependent phase in
$\mathcal{T}$).

Combining these transformations together yields the matrix:

$\left(\frac{1+i}{2}\right)\begin{bmatrix}i\tau_{0}&-\tau_{y}\\\
-\tau_{0}&i\tau_{y}\end{bmatrix}$ (33)

Due to this complicated transformation, the $\tau$ and $\sigma$ matrices in
this basis refer to some combination of the original degrees of freedom in a
non-obvious way. (I.e.: we use these symbols to refer to the indices in the
matrix structure rather than their original degrees of freedom.)

With our Hamiltonian and time-reversal symmetry both in standard form, only
the $C_{3}$ rotation operator remains nontrivial. That the operator is a
symmetry of the Hamiltonian determines how it acts on $\tau$ space, so what
remains is to determine how it acts on $\sigma$ space, which is precisely
specifying the matrix $U$ in Eq. (23)). Computing $U$ in this case yields a
standard rotation matrix:

$U_{3,HG}=\exp(i\pi\sigma_{y}/3)$ (34)

####  A.3.4 Basis change in graphene

We now perform an analogous series of steps in our full graphene model. There
is some additional freedom in our basis change here that we exploit to make
our rotational symmetries as simple as possible. In the end, our basis
transformation is:

$\left(\frac{1+i}{2}\right)\begin{bmatrix}i\tau_{0}&0&0&-\tau_{y}\\\
-\tau_{0}&0&0&i\tau_{y}\\\ 0&i\tau_{0}&\tau_{y}&0\\\
0&\tau_{0}&i\tau_{y}&0\end{bmatrix}$ (35)

With our basis transformations performed, we can now write the rotations in
the new basis. Again, we take Eq. (23) and specify the $U$ matrices defined
therein, this time both for $C_{3}$ and $C_{2}$ symmetry:

$U_{3}=\exp(i\pi\sigma_{y}/3)\mu_{0}$ (36) $U_{2}=\sigma_{0}\mu_{x}$ (37)

That these expressions factor is precisely the motivation for working with
$C_{3}$ and $C_{2}$ rather than $C_{6}$: $C_{3}$ now acts trivially on the
$\mu$ indices and $C_{2}$ acts trivially on the $\sigma$ indices, whereas
$C_{6}$ acts nontrivially on both indices simultaneously.

Moreover, we can now see the exact sense in which our half-graphene model was,
indeed, half of graphene: our full model of graphene can be written as two
$C_{2}$-related copies of our half-graphene model. This will allow us to work
out the details of the $C_{3}$ symmetry in our half-graphene model in Sec.
C.3, then add the details of $C_{2}$ symmetry in the full graphene model in
Sec. C.4.

##  Appendix B Interlayer Couplings

In this appendix, we derive the coupling between two layers stacked with a
relative twist angle. We begin with a general derivation of interlayer
couplings in a tight-binding model given the lattice vectors of each layer and
then specialize to the case of small twist angle about a commensurate crystal
structure. We then specialize to the case of coupled Dirac cones and
incorporate the constraints of time-reversal and rotational symmetry.

###  B.1 General two-layer tight-binding couplings

We begin by considering the interlayer coupling between two layers with
different lattices, with no constraints on the Hamiltonian of each layer.

Letting $I,J$ index all non-momentum degrees of freedom (sublattice, orbital,
spin, etc., but not valley, since it is included in the integral over $k$) of
layer 1 and $I^{\prime},J^{\prime}$ of layer 2, the two-layer Hamiltonian can
generally be written:

$H=\sum_{I,J}\int_{\text{BZ}\
1}d^{2}k\psi_{1,I,k}^{\dagger}H_{1}^{IJ}(k)\psi_{1,J,k}+\sum_{K,L}\int_{\text{BZ}\
2}d^{2}k^{\prime}\psi_{2,I^{\prime},k^{\prime}}^{\dagger}H_{2}^{I^{\prime}J^{\prime}}(k^{\prime})\psi_{2,J^{\prime},k^{\prime}}\\\
+\left[\sum_{I,J^{\prime}}\int_{\text{BZ}\ 1}d^{2}k\int_{\text{BZ}\
2}d^{2}k^{\prime}\left(\psi_{1,I,k}^{\dagger}T^{IJ^{\prime}}(k,k^{\prime})\psi_{2,J^{\prime},k^{\prime}}\right)+h.c.\right]$
(38)

We now simplify $T(k,k^{\prime})$ by imposing some minimal assumptions. We
start by Fourier transforming $T^{IJ^{\prime}}(k,k^{\prime})$:

$T^{IJ^{\prime}}(k,k^{\prime})=\sum_{R,R^{\prime}}e^{-ik\cdot(R+r_{I})}e^{ik^{\prime}\cdot(R^{\prime}+r_{J^{\prime}})}T^{IJ^{\prime}}(R,R^{\prime})$
(39)

In the above equation, $R$ and $R^{\prime}$ denote the centers of the unit
cells of the two layers, whereupon $r_{I}$ and $r_{J}^{\prime}$ denote the
location of sublattices $I$ and $J^{\prime}$ relative to the centers of their
respective unit cells.

We then make a standard moiré tight-binding assumption that the crystal
locally has the translation invariance of the original lattice, i.e.:

$T^{IJ^{\prime}}(R,R^{\prime})=t^{IJ^{\prime}}(R+r_{I}-R^{\prime}-r_{J^{\prime}})$
(40)

Combining this equation with the previous one and Fourier transforming again
allows us to simplify significantly:

$\begin{split}T^{IJ^{\prime}}(k,k^{\prime})=&\sum_{R,R^{\prime}}e^{-ik\cdot(R+r_{I})}e^{ik^{\prime}\cdot(R^{\prime}+r_{J^{\prime}})}t^{IJ^{\prime}}(R+r_{I}-R^{\prime}-r_{J^{\prime}})\\\
=&\sum_{R,R^{\prime}}e^{-ik\cdot(R+r_{I})}e^{ik^{\prime}\cdot(R^{\prime}+r_{J^{\prime}})}\sum_{k^{\prime\prime}}e^{ik^{\prime\prime}\cdot(R+r_{I}-R^{\prime}-r_{J^{\prime}})}t^{IJ^{\prime}}_{k^{\prime\prime}}\\\
=&\sum_{R,R^{\prime}}\sum_{k^{\prime\prime}}e^{-i(k-k^{\prime\prime})\cdot
R}e^{i(k^{\prime}-k^{\prime\prime})\cdot
R^{\prime}}e^{i(k^{\prime\prime}-k)\cdot
r_{I}}e^{-i(k^{\prime\prime}-k^{\prime})\cdot
r_{J^{\prime}}}t^{IJ^{\prime}}_{k^{\prime\prime}}\\\
=&\sum_{G,G^{\prime}}\sum_{k^{\prime\prime}}\delta_{k+G,k^{\prime\prime}}\delta_{k^{\prime}+G^{\prime},k^{\prime\prime}}e^{i(k^{\prime\prime}-k)\cdot
r_{I}}e^{-i(k^{\prime\prime}-k^{\prime})\cdot
r_{J^{\prime}}}t^{IJ^{\prime}}_{k^{\prime\prime}}\\\
=&\sum_{G,G^{\prime}}\delta_{k+G,k^{\prime}+G^{\prime}}e^{i(G\cdot
r_{I}-G^{\prime}\cdot r_{J^{\prime}})}t^{IJ^{\prime}}_{k+G}\end{split}$ (41)

In the above expressions, $k^{\prime\prime}$ is allowed to run over all of
momentum space, whereas $G$ and $G^{\prime}$ are the reciprocal lattice
vectors of the two layers.

This final expression relates the coupling between two materials of
arbitrarily mismatched lattices to the Fourier transform of the tight-binding
interlayer hopping $t^{IJ}(r)$.

###  B.2 Small-angle twisting

We now look specifically at the case where the two layers are twisted with a
small angle relative to a commensurate cell and, further, that the low-energy
physics is well described by a $k\cdot p$ expansion about a particular point
$k_{0}$ in the BZ. (The derivation proceeds similarly for a small lattice
mismatch instead of a small twist angle.)

Making the substitutions $k\rightarrow k_{0}+\bar{k}$,
$k^{\prime}=k_{0}^{\prime}+\bar{k}^{\prime}$, Eq. (41) becomes:

$T^{IJ^{\prime}}(k,k^{\prime})=\sum_{G,G^{\prime}}\delta_{\bar{k}+k_{0}+G,\bar{k}^{\prime}+k_{0}^{\prime}+G^{\prime}}e^{i(G\cdot
r_{I}-G^{\prime}\cdot r_{J^{\prime}})}t^{IJ^{\prime}}_{\bar{k}+k_{0}+G}$ (42)

From here, we explicitly separate out twist angle dependence: we split $k_{0}$
as $k_{0}=\hat{k}_{0}+\delta k_{0}$, where $\hat{k}_{0}$ is twist-angle-
independent and $|\delta k_{0}|\sim\theta|k_{0}|$. Therefore, for small twist
angle, $\delta k_{0}$ will be small. (As an example, in graphene,
$\hat{k}_{0}$ would be the position of the $K$ point before twisting, and
$\delta k_{0}$ would be its deviation after twisting.) We do a similar
decomposition with $G$, $k_{0}^{\prime}$, and $G^{\prime}$.

We now suppose that the interlayer spacing is much greater than the
interatomic spacing, which motivates an assumption that $t$ decays slowly in
position space and therefore quickly in momentum space [5]. Then, the smallest
values for $\bar{k}+k_{0}+G$ (i.e., within the first few BZ) dominate.
Combined with our assumption for small $\bar{k}$, this justifies the following
decomposition:

$\delta_{\bar{k}+k_{0}+G,\bar{k}^{\prime}+k_{0}^{\prime}+G^{\prime}}\simeq\delta_{\hat{k}_{0}+\hat{G},\hat{k}_{0}^{\prime}+\hat{G}^{\prime}}\delta_{\bar{k}+\delta
k_{0}+\delta G,\bar{k}^{\prime}+\delta k_{0}^{\prime}+\delta G^{\prime}}$ (43)

The first term is $\theta$-independent: it implies that cones will only couple
after twisting if they fold onto the same point in the BZ before twisting,
allowing us a simple way to determine what valleys do or do not couple. The
second term then specifies the twist angle dependence, as it contains the
$\theta$-dependent terms.

To simplify, we explicitly input the rotation by declaring $\delta
G^{\prime}=M\hat{G}^{\prime}$, $\delta k_{0}^{\prime}=M\hat{k}_{0}^{\prime}$,
and $\delta G=\delta k_{0}=0$, where $M$ is the difference between the
rotation matrix and the identity. Then, with an approximation in $t$ to lowest
order in $\bar{k}$, we end up with the following expression for
$T^{IJ^{\prime}}(k,k^{\prime})$:

$\sum_{\hat{G},\hat{G}^{\prime}}\underbrace{\delta_{\hat{k}_{0}+\hat{G},\hat{k}_{0}^{\prime}+\hat{G}^{\prime}}}_{\text{Valley
matching}}\underbrace{\delta_{\bar{k},\bar{k}^{\prime}+M(\hat{k}_{0}^{\prime}+\hat{G}^{\prime})}}_{\theta/Q\text{-dependence}}\underbrace{e^{i(\hat{G}\cdot\hat{r}_{I}-\hat{G}^{\prime}\cdot\hat{r}_{J^{\prime}})}t^{IJ^{\prime}}_{\hat{k}_{0}+\hat{G}}}_{\text{Coupling
strength}}$ (44)

where $\hat{r}_{I,J^{\prime}}$ denote the positions of the sublattice degrees
of freedom before being twisted. Finally, we define
$Q:=\bar{k}^{\prime}-\bar{k}$ and express this interlayer hopping as $T_{Q}$,
giving the notation from the main text in Eq. (3).

###  B.3 Symmetry constraints on $T_{Q}$

We now derive the effect of symmetry on the interlayer coupling terms. In so
doing, we explicitly assume the Dirac Hamiltonian is written in the basis
discussed in Appendix A, including the specified form of the symmetry
operators.

For simplicity of expression, we explicitly break the $T_{Q}$ into a set of
2-by-2 matrices $T_{Q,ij^{\prime}}$. Note that $ij^{\prime}$ used here are not
the same as $IJ^{\prime}$ in the previous section: $IJ^{\prime}$ run over all
non-valley degrees of freedom, whereas $ij^{\prime}$ run over time-reversal-
invariant Dirac cones (which may include valley degrees of freedom).

Time-reversal symmetry relates $Q$ to $-Q$, and gives the relation:

$\sigma^{y}T^{*}_{ik,-Q}\sigma^{y}=T_{ik,Q}$ (45)

For in-plane rotations, we return to the $U$ matrices defined following Eq.
(23), which may be different for each layer. Call the matrix for one layer $U$
and for the other $U^{\prime}$. Implementing the rotation in momentum space by
$Q\mapsto\mathfrak{R}Q$, the resulting constraint is:

$e^{i\theta\sigma_{z}/2}[U_{ij}T_{jk^{\prime},\mathfrak{R}Q}U^{\prime\dagger}_{k^{\prime}l^{\prime}}]e^{-i\theta\sigma_{z}/2}=T_{il^{\prime},Q}$
(46)

Notice that the constraints imposed by time-reversal and rotation symmetry in
Eqs. (45) and (46) relate the interlayer coupling at different $Q$. Thus, once
$T_{Q}$ is known for one choice of $Q$, it is determined for all other $Q$
related by time-reversal or rotation symmetry. The exception to this is the
combined symmetry operation $C_{2}\mathcal{T}$, which leaves $Q$ invariant and
(due to antilinearity) acts as a reality constraint on each individual
$T_{Q}$.

In the case of identical materials, there may also be layer-interchanging
rotations, like $C_{2x}$ or $C_{2y}$. These follow a similar formula to the
in-plane rotations, except instead of relating $T$ to $T$, they relate $T$ to
$T^{\dagger}$. Letting $\phi$ denote the angle between the axis of rotation
and the $x$-axis, define:

$\sigma_{\phi}=\cos(\phi)\sigma_{x}+\sin(\phi)\sigma_{y}$ (47)

Then, letting $\mathfrak{L}Q$ denote the action of the rotation on $Q$ and $V$
a material-dependent transformation matrix representing the action of the
symmetry on the $ij$ indices:

$\sigma_{\phi}[V_{ij}T_{jk^{\prime},-\mathfrak{L}Q}^{\dagger}V^{\dagger}_{k^{\prime}l^{\prime}}]\sigma_{\phi}=T_{il^{\prime},Q}$
(48)

There is a subtlety here in the formula involving $T_{-\mathfrak{L}Q}$ instead
of $T_{\mathfrak{L}Q}$. Specifically, $Q$ is defined as the difference in
momentum from layer 1 to layer 2 (see below Eq. (44)). When the two layers are
exchanged, the meaning of $Q$ has to change accordingly, which swaps
$Q\rightarrow-Q$.

In practice, it is only necessary to consider one out-of-plane rotation
symmetry, since the others are generated by products with the in-plane
rotations. When considering a particular $Q$, it is prudent to consider an
out-of-plane rotation that preserves that $Q$ (as done in Sec. III.1), so that
the out-of-plane rotation imposes a constraint on $T_{Q}$, rather than
relating two distinct terms $T_{Q}$ and $T_{-\mathfrak{L}Q}$. This then (being
an antilinear constraint on a single $Q$, akin to $C_{2}\mathcal{T}$) imposes
a reality constraint on individual $T_{Q}$. A material with both
$C_{2}\mathcal{T}$ and such an out-of-plane symmetry may have conflicting
reality constraints that can cause certain interlayer hopping terms to vanish,
as again can be seen in Sec. III.1.

These symmetries are the only ones that we expect to arise in generic models,
since reflection symmetries are (generically) broken in twisted bilayers and
internal symmetries are model-specific. Without further constraints, we expect
all interlayer hopping terms that satisfy Eqs. (45), (46) and (48) to arise,
although symmetry does not constrain their amplitudes.

##  Appendix C Calculations of Self-Energy

In the following appendix, we compute the self-energy of a Dirac cone in one
layer of material that results from hopping into the other layer. We show that
the self-energy takes the general form given by Eq. (4) and derive the
coefficients defined in Eq. (8). We also derive the more general results that
arise when there is more than one Dirac cone in each layer.

We begin by discussing the general case in the special basis described in
Appendix A (in particular parts A.1 and A.2). We then discuss the extent to
which this calculation requires this change of basis and the extent to which
it can be ignored.

With the general case worked out, we proceed to compute the quadratic
corrections to graphene (and our simplified half-graphene model) described in
A.3.

###  C.1 General calculation

We first calculate the self-energy of a Dirac cone in layer 1 due to its
coupling to a Dirac cone in layer 2; the self-energy of the Dirac cone in
layer 2 from hopping into layer 1 is identical up to the definition of $Q$
(which is the momentum shift in going from layer 1 to layer 2, not the
converse) and an interchange of the roles of $T$ and $T^{\dagger}$. The
calculation of self-energy, to lowest-order in perturbation theory, proceeds
as follows:

$\displaystyle\Sigma_{ij}(\bar{k},\omega)=$
$\displaystyle\sum_{Q,l^{\prime}}T_{Q,il^{\prime}}[\omega-v(\bar{k}+Q)\cdot\sigma]^{-1}T_{Q,jl^{\prime}}^{\dagger}$
(49) $\displaystyle\simeq$
$\displaystyle-\sum_{Q,l^{\prime}}\frac{1}{v^{2}Q^{4}}T_{Q,il^{\prime}}[Q^{2}\omega+v(Q^{2}(\bar{k}\cdot\sigma)-2(Q\cdot\bar{k})(Q\cdot\sigma)+Q^{2}(Q\cdot\sigma))]T_{Q,jl^{\prime}}^{\dagger}$
(50) $\displaystyle=$
$\displaystyle-\sum_{Q,l^{\prime}}\frac{1}{2v^{2}Q^{4}}\left\\{\left(T_{Q,il^{\prime}}T^{\dagger}_{Q,jl^{\prime}}+T_{-Q,il^{\prime}}T^{\dagger}_{-Q,jl^{\prime}}\right)Q^{2}\omega+vQ^{2}\left(T_{Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{Q,jl^{\prime}}-T_{-Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{-Q,jl^{\prime}}\right)\right.$
$\displaystyle\qquad\left.+v\left[T_{Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{Q,jl^{\prime}}+T_{-Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{-Q,jl^{\prime}}\right]\right\\}$ (51)
$\displaystyle=$
$\displaystyle-\sum_{Q,l^{\prime}}\frac{1}{2v^{2}Q^{4}}\left\\{\left(T_{Q,il^{\prime}}T^{\dagger}_{Q,jl^{\prime}}+\sigma_{y}(T_{Q,il^{\prime}}T^{\dagger}_{Q,jl^{\prime}})^{*}\sigma_{y}\right)Q^{2}\omega+vQ^{2}\left(T_{Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{Q,jl^{\prime}}+\sigma_{y}(T_{Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{Q,jl^{\prime}})^{*}\sigma_{y}\right)\right.$
$\displaystyle\qquad\left.+v\left(T_{Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{Q,jl^{\prime}}-\sigma_{y}\left(T_{Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{Q,jl^{\prime}}\right)^{*}\sigma_{y}\right)\right\\}$
(52)

Let us now explain each equality: Eq. (50) is derived by expanding Eq. (49) to
linear order in the dimensionless parameters $\omega/v|Q|$ and
$|\bar{k}|/|Q|$. Eq. (51) is derived by averaging the terms $T_{Q}$ and
$T_{-Q}$, and Eq. (52) relates $T_{-Q,il^{\prime}}$ to $T_{Q,il^{\prime}}$ by
the constraint imposed by time reversal symmetry in Eq. (45).

The above equation for the self-energy expanded to linear order in $\bar{k}$
and $\omega$ can be decomposed as:

$\Sigma_{ij}(\bar{k},\omega)=\mathcal{A}_{ij}\omega+\mathcal{B}_{ij}+\bar{k}\cdot\mathcal{V}_{ij},$
(53)

where the coefficients can be further decomposed as:

$\displaystyle\mathcal{A}_{ij}=$ $\displaystyle
A_{ij}\sigma_{0}+iA^{\prime}_{ij}\sigma_{z}+iA^{\prime\prime}_{ij}\cdot\sigma$
(54a) $\displaystyle\mathcal{B}_{ij}=$ $\displaystyle
B_{ij}\sigma_{0}+iB^{\prime}_{ij}\sigma_{z}+iB^{\prime\prime}_{ij}\cdot\sigma$
(54b) $\displaystyle\mathcal{V}_{ij}=$ $\displaystyle
V_{ij}\sigma_{z}+iV^{\prime}_{ij}\sigma_{0}+M_{ij}\cdot\sigma$ (54c)

In the above equations, aside from the $ij$ indices, some of the coefficients
have implicit indices implied by the dot products, i.e., $A$, $A^{\prime}$,
$B$, and $B^{\prime}$ are scalars; $A^{\prime\prime}$, $B^{\prime\prime}$,
$V$, and $V^{\prime}$ are two-component vectors; and $M$ is a rank-2 tensor.

The factors of $i$ have been chosen so that all are real due to time-reversal
symmetry. Since the self-energy is Hermitian, all unprimed coefficients are
symmetric in the $ij$ indices, whereas all primed or double-primed indices are
antisymmetric. (In particular, the antisymmetric terms vanish on the $ij$
diagonal, and therefore do not appear in the case of a single-Dirac cone, for
which $i=j=1$).

We rewrite the above equations in terms of traces over the $T$ matrices using
the following identities for 2x2 Hermitian matrices $H$:

$\displaystyle H+\sigma_{y}H^{*}\sigma_{y}=$ $\displaystyle\Tr[H]\sigma_{0}$
(55a) $\displaystyle H-\sigma_{y}H^{*}\sigma_{y}=$
$\displaystyle\Tr[H\sigma^{\mu}]\sigma_{\mu}$ (55b)

By applying these identities to Eq. (52) and isolating the coefficients
defined in Eqs (53)-(54), we find the following explicit expressions for the
individual coefficients defined in Eq. (54):

$\displaystyle A_{ij}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2v^{2}Q^{2}}\real\Tr[T_{Q,ii^{\prime}}T_{Q,ji^{\prime}}^{\dagger}]$
(56a) $\displaystyle A_{ij}^{\prime}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2v^{2}Q^{2}}\imaginary\Tr[T_{Q,ii^{\prime}}T_{Q,ji^{\prime}}^{\dagger}\sigma^{z}]$
(56b) $\displaystyle A_{ij}^{\prime\prime\mu}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2v^{2}Q^{2}}\imaginary\Tr[T_{Q,ii^{\prime}}T_{Q,ji^{\prime}}^{\dagger}\sigma^{\mu}]$
(56c)

$\displaystyle B_{ij}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\real\Tr[T_{Q,ii^{\prime}}(Q\cdot\sigma)T_{Q,ji^{\prime}}^{\dagger}]$
(57a) $\displaystyle B_{ij}^{\prime}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\imaginary\Tr[T_{Q,ii^{\prime}}(Q\cdot\sigma)T_{Q,ji^{\prime}}^{\dagger}\sigma^{z}]$
(57b) $\displaystyle B_{ij}^{\prime\prime\mu}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\imaginary\Tr[T_{Q,ii^{\prime}}(Q\cdot\sigma)T_{Q,ji^{\prime}}^{\dagger}\sigma^{\mu}]$
(57c)

$\displaystyle V_{ij}^{\mu}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\left(\delta^{\mu}_{\
\lambda}-2\frac{Q^{\mu}Q_{\lambda}}{Q^{2}}\right)\real\Tr[T_{Q,ii^{\prime}}\sigma^{\lambda}T_{Q,ji^{\prime}}^{\dagger}\sigma^{z}]$
(58a) $\displaystyle V_{ij}^{\prime\mu}$
$\displaystyle=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\left(\delta^{\mu}_{\
\lambda}-2\frac{Q^{\mu}Q_{\lambda}}{Q^{2}}\right)\imaginary\Tr[T_{Q,ii^{\prime}}\sigma^{\lambda}T_{Q,ji^{\prime}}^{\dagger}]$
(58b)

$M_{ij}^{\mu\nu}=\sum_{Q,i^{\prime}}\frac{1}{2vQ^{2}}\left(\delta^{\mu}_{\
\lambda}-2\frac{Q^{\mu}Q_{\lambda}}{Q^{2}}\right)\real\Tr[T_{Q,ii^{\prime}}\sigma^{\lambda}T_{Q,ji^{\prime}}^{\dagger}\sigma^{\nu}]$
(59)

So far, we have only imposed time-reversal symmetry (in the form noted in Sec.
A.2, using the constraints derived in Sec. B.3). To further constrain the form
of the self-energy, we should impose rotational symmetry. While we cannot
further simplify without knowing the specific form of the rotation matrices
$U$ (defined in Appendix A.2), which are material-dependent, we now outline
the procedure for simplifying further.

The core aspect of the simplification is relating the terms in the sums over
$Q$ in Eqs (56)-(59) using Eq. (23). Writing, for instance, Eq. (56a) as
$A_{ij}=\sum_{Q}A_{ij}(Q)$ and noting
$A_{il}(\mathfrak{R}Q)=U_{ij}A_{jk}(Q)U_{kl}^{\dagger}$, it follows that
$A_{ij}$ can be expressed in terms of only a single representative $Q_{0}$ (of
each symmetry-related-set of $Q$), i.e., we need only compute $A_{ij}(Q_{0})$.
Summing over different $Q$ is then equivalent to picking out the rotation-
invariant part of $A_{ij}(Q_{0})$ (under the representation of rotations given
by the $U$ matrices) and multiplying by an appropriate symmetry factor.

The same procedure can be applied to the rest of Eqs. (56)-(59) very
similarly. The only modification required is incorporating the $\mu$ and $\nu$
indices, where present: $A^{\prime\prime}(Q)$, etc., pick up a (vector)
transformation on those indices in addition to the $U$ matrices that transform
the $ij$ indices. This makes the decomposition into irreps different, but the
philosophy is the same: find the rotation-invariant part of those matrices and
multiply by an appropriate symmetry factor.

It is often useful to decompose $M$ into the irreps of the $\mu\nu$ indices
(before considering the $ij$ indices): the trace $C$, the antisymmetric part
$D$, and the symmetric trace-free part $N$, i.e.,

$M^{\mu\nu}_{ij}=C_{ij}\delta^{\mu\nu}+D_{ij}\epsilon^{\mu\nu}+N_{ij}^{\mu\nu}$
(60)

We can be more concrete in the specific case where there is only one Dirac
cone. In this case, all primed coefficients vanish as discussed below Eq.
(54). In addition, $U$ is simply a phase factor. Thus, we only need consider
the rotational transformations of $A$, $B$, $V$, and $M$, all of which are
specified by $U$ being a phase: $A(Q)$ and $B(Q)$ are identical for all $Q$;
$V(Q)$ transforms like a vector; and $M(Q)$ transforms like a two-index
tensor.

In the sum over $Q$, therefore, the $V(Q)$ will cancel for any rotational
symmetry (including $C_{2}$): denoting the rotation operation on a vector by
$\mathfrak{R}$, our above calculation reveals that
$V(\mathfrak{R}^{n}Q)=\mathfrak{R^{n}}V(Q)$, and so the sum of rotation-
related $Q$ terms vanishes ($\sum_{n}V(\mathfrak{R}^{n}Q)=0$).

A similar but more complicated representation theory argument reveals that $M$
is constrained to only its trace and antisymmetric parts $C$ and $D$ (i.e.,
$N$ vanishes). This yields the single-cone result given in Eq. (4), with
coefficients as given in Eq. (8).

###  C.2 $U(N)$ symmetry

The above calculation, in principle, only holds in the special basis in which
both materials have a $k\cdot\sigma$ Hamiltonian and time reversal takes the
form $i\sigma_{y}K$. This requires a change of basis that is inconvenient for
developing intuition in materials where the cones do not naturally come in
such a form, such as graphene (where time-reversal does not naturally flip the
internal degree of freedom of a Dirac cone).

However, using our final results (56)–(59), we find that the expression for
self-energy of the first layer is invariant under a change of basis in the
second layer that acts only on the $i^{\prime}j^{\prime}$ indices, since they
are summed over. Thus, our result is invariant under the change of basis
indicated in A.2 (for the second layer), although not under the change of
basis described in Appendix A.1.

In this paper, this invariance allows us to avoid changing the graphene basis
in Sec. IV; otherwise, we would be forced to make all the changes of basis
highlighted in Appendix A.3. Instead, we only need transform the Hamiltonian
to a $vk\cdot\sigma$ form (as described in Appendix A.1), which is trivial for
the purposes of our calculation.

However, note the full basis change on the first layer (the one for which we
are calculating the self-energy) is still necessary: a unitary change of basis
on the $ij$ indices will not necessarily leave $\Sigma_{ij}$ invariant. Hence,
for the computation of the self-energy of the graphene cones Appendices C.3
and C.4, the complex changes of basis are still necessary.

###  C.3 Half-graphene corrections

We now apply the above computational method to the case of the corrections to
graphene from tunneling into the Dirac cone on the surface of a TI, assuming
for simplicity that only the smallest set of $Q$s contribute nonnegligibly to
the sums in Eqs. (56)-(59). We begin with the simpler half-graphene model
presented in Sec. A.3.

For each of the terms presented in Eqs. (56)-(59), we can divide the
contributions into rotational irreps (where the action of the rotations are
defined by the action of $U$ on the $ij$ indices combined with the usual
action on any vector indices $\mu\nu$). For the trivial irreps, the summands
for different (symmetry-related) $Q$s will all be identical, and hence we can
simply take a representative $Q$ and multiply by a symmetry factor of $6$. For
nontrivial irreps, the summands will cancel.

Consider first the corrections from $A$, $A^{\prime}$, $B$, $B^{\prime}$, $C$,
and $D$ (as defined in Eqs. (54) and (60)). For these terms, since there are
no vector indices, any irreps with $ij$ indices that transform nontrivially
under $C_{3}$ will vanish in the sum over $Q$, which means (using the
representation of rotational symmetry given in Eq. (34)) these terms can only
contribute $\sigma_{0}$ or $\sigma_{y}$ terms; the pair
$(\sigma_{x},\sigma_{z})$ transform like a vector instead. By symmetry in
$ij$, the unprimed terms contribute $\sigma_{0}$ and the primed $\sigma_{y}$.
Therefore, the contributions from these terms to the self-energy in Eq. (53)
can be written as:

$(A_{0}\omega+B_{0})\tau_{0}\sigma_{0}+(A_{0}^{\prime}\omega+B_{0}^{\prime})\tau_{z}\sigma_{y}+C_{0}(\bar{k}\cdot\tau)\sigma_{0}+D_{0}(\bar{k}\times\tau)\sigma_{0}$

Next, consider the terms with vector indices: $A^{\prime\prime\mu}_{ij}$,
$B^{\prime\prime\mu}_{ij}$, $V^{\mu}_{ij}$, and $V^{\prime\mu}_{ij}$, (as
defined in eq. (54)). Again, we find the irreps that transform trivially. As
$\mu$ is a vector index, another vector index from the $ij$ is required, and
such an index can be found in the vector combination
$(\sigma^{x},\sigma^{z})$. However, since these matrices are symmetric in
$ij$, the only term they can contribute to is $V$, and so
$A^{\prime\prime}=B^{\prime\prime}=V^{\prime}=0$.

The surviving term, $V$, we have now decomposed as
$V^{\mu}_{ij}=V^{\mu}_{a}\sigma^{a}_{ij}$ where $a=x,z$ and $\mu=x,y$. The
rotationally-invariant terms are those that contract the $\mu$ and $a$ indices
via dot and cross products, as follows: for any matrix $S^{a}_{\mu}$, the
rotation-invariant parts are the contractions of $V^{\mu}_{a}$ with
$S_{\mu}^{a}$ and $\epsilon_{\mu}^{\ \nu}S_{\nu}^{a}$.

Therefore, the contribution to the self-energy from $V$ decomposes as:

$(\bar{k}\cdot
V_{ij})=V_{1}(\bar{k}_{x}\sigma^{x}_{ij}+\bar{k}_{y}\sigma^{z}_{ij})+V_{2}(\bar{k}_{x}\sigma^{z}_{ij}-\bar{k}_{y}\sigma^{x}_{ij})$

This choice of $V_{1}$ and $V_{2}$ is, of course, not unique, but the
invariant subspace will be the same for all choices of decomposition.

We can similarly decompose our last coefficient, $N^{\mu\nu}_{ij}$: the
$\mu\nu$ indices are a combination of $\alpha^{x}$ and $\alpha^{z}$, which
transform as a vector, and therefore we must take a corresponding combination
of $\sigma^{x}$ and $\sigma^{z}$, as before. Therefore, we can express our
solution as:

$N^{\mu\nu}_{ij}=N_{1}(\alpha_{x}^{\mu\nu}\sigma^{x}_{ij}+\alpha_{z}^{\mu\nu}\sigma^{z}_{ij})+N_{2}(\alpha_{x}^{\mu\nu}\sigma^{z}_{ij}-\alpha_{z}^{\mu\nu}\sigma^{x}_{ij})$

Concatenating these all together, we find the corrections to the self-energy
for our half-graphene model (in the special basis of Appendix A.3.3) can be
written in terms of ten coefficients:

$\begin{split}\Sigma_{HG}(\bar{k},\omega)=&(A_{0}\omega+B_{0})\tau_{0}\sigma_{0}+(A_{0}^{\prime}\omega+B_{0}^{\prime})\tau_{z}\sigma_{y}\\\
+&C_{0}\sigma_{0}(\bar{k}\cdot\tau)+D_{0}\sigma_{0}(\bar{k}\times\tau)\\\
+&\left[V_{1}(\bar{k}_{x}\sigma^{x}+\bar{k}_{y}\sigma^{z})+V_{2}(\bar{k}_{x}\sigma^{z}-\bar{k}_{y}\sigma^{x})\right]\tau_{z}\\\
+&(N_{1}\sigma^{x}+N_{2}\sigma^{z})(\bar{k}_{x}\tau_{y}+\bar{k}_{y}\tau_{x})\\\
+&(N_{1}\sigma^{z}-N_{2}\sigma^{x})(\bar{k}_{x}\tau_{x}-\bar{k}_{y}\tau_{y})\end{split}$
(61)

In the original basis, the self-energy takes the form:

$\begin{split}\Sigma_{HG}(\bar{k},\omega)=&(A_{0}\omega+B_{0})\tau_{0}\sigma_{0}+(A_{0}^{\prime}\omega+B_{0}^{\prime})\tau_{z}\sigma_{0}\\\
+&C_{0}\begin{bmatrix}\bar{k}\cdot\tau&0\\\
0&-\bar{k}\cdot\tau^{*}\end{bmatrix}+D_{0}\begin{bmatrix}\bar{k}\times\tau&0\\\
0&-\bar{k}\times\tau^{*}\end{bmatrix}\\\
+&\left[V_{1}(\bar{k}_{x}\sigma^{y}+\bar{k}_{y}\sigma^{x})+V_{2}(\bar{k}_{x}\sigma^{x}-\bar{k}_{y}\sigma^{y})\right]\tau_{x}\\\
+&(\tau_{0}-\tau_{z})(N_{1}(k\cdot\sigma)-N_{2}(k\times\sigma))\end{split}$
(62)

Some of these terms may vanish to lowest order in perturbation theory; Eq.
(61) is the most general form constrained by symmetry.

For spin-preserving sublattice-independent hopping with amplitude $t$, we can
explicitly compute these terms by by multiplying the term from $Q=|Q|\hat{x}$
(for twist) or $Q=|Q|\hat{x}$ (for lattice mismatch) by a symmetry factor of
six. We find that (in either case) the only nonvanishing coefficient in the
self-energy is:

$A_{0}=\frac{3t^{2}}{v^{2}Q^{2}}$ (63)

Hence, in this half-graphene model with spin-preserving sublattice-independent
hopping to the TI layer, the effect of the lowest-order self-energy
corrections is an overall energy rescaling by a factor of $1+A_{0}$.

####  C.3.1 Corrections to TI cone from spin-preserving hopping

For sake of completeness, we here include the self-energy corrections to the
TI Dirac cone with spin-preserving hopping in this model. Similar to graphene
(as discussed in Sec. IV.1), we can split $T_{Q}$ into $T_{Q,K}$ and
$T_{Q,K^{\prime}}$ (with no extra spin index in this case). Taking a $Q$ for
which $T_{Q,K}$ is nonvanishing, we find that it takes one of the forms:

$\displaystyle T_{Q,K,\uparrow}=\begin{bmatrix}a&b\\\ 0&0\end{bmatrix}$ (64a)
or $\displaystyle T_{Q,K,\downarrow}=\begin{bmatrix}0&0\\\ a&b\end{bmatrix}$
(64b)

All other $T_{Q}$ are symmetry-related, as in the graphene case. Note unlike
graphene, the two coefficients $a$ and $b$ can be different due to the lack of
$C_{2}$ symmetry.

Like graphene, the $C$ and $D$ coefficients vanish, but the $A$ and $B$
coefficients do not. As a consequence of these coefficients, the energy shift
is:

$\Delta E=-\frac{6abv_{HG}|Q|}{v_{HG}^{2}Q^{2}+3(a^{2}+b^{2})}$ (65)

###  C.4 Full graphene corrections

We now work with our full model of graphene, using the half-graphene case as a
starting point. We work in the basis described in Appendix A.3.4.

We first decompose the indices $ij$ into sub-indices for the $\mu$ and
$\sigma$ degrees of freedom. For example, $A_{ij}=A_{pa,qb}$, where the $\mu$
indices are $\mu_{pq}$ and the $\sigma$ indices are $\sigma_{ab}$. As
discussed below Eqs. (36) and (37), the unitary rotation matrix $U_{3}$ only
acts nontrivially on the $ab$ indices, while $U_{2}$ only acts nontrivially on
the $pq$ indices, which allows us to straightforwardly separate the two
symmetries (rather than working with a more complicated
$U_{6}=U_{2}U_{3}^{-1}$ matrix).

As we did for the half-graphene case, we begin by considering the scalar
coefficients $A$, $B$, $C$, $D$, $A^{\prime}$, and $B^{\prime}$. We now filter
down to the subrepresentations trivial under both $C_{3}$ and $C_{2}$. $C_{3}$
implies each coefficient must commute with $\sigma_{y}$ and $C_{2}$ implies it
must commute with $\tau_{x}$, by similar logic to that used in Appendix C.3.
Considering these in the context of the symmetry properties of the different
matrices yields the following contribution to the self-energy:

$[(A_{1}\mu_{0}+A_{2}\mu_{x})\sigma_{0}\tau_{0}+(A_{1}^{\prime}\mu_{0}+A_{2}^{\prime}\mu_{x})\sigma_{y}\tau_{z}]\omega\\\
+(B_{1}\mu_{0}+B_{2}\mu_{x})\sigma_{0}\tau_{0}+(B_{1}^{\prime}\mu_{0}+B_{2}^{\prime}\mu_{x})\sigma_{0}\tau_{z}\\\
+(C_{1}\mu_{0}+C_{2}\mu_{x})\sigma_{0}(k\cdot\tau)+(D_{1}\mu_{0}+D_{2}\mu_{x})\sigma_{0}(k\times\tau)$

We next consider the vector parts $A^{\prime\prime}$, $B^{\prime\prime}$, $V$,
and $V^{\prime}$. The vector index must pair with a combination of $ab$ and
$pq$ indices which transforms like a vector under both $C_{3}$ and $C_{2}$. We
use $(\sigma_{x},\sigma_{z})$ as a vectorlike pair as before for our $ab$
indices, and either $\mu_{y}$ or $\mu_{z}$ in our $pq$ indices (to anticommute
with $\mu_{x}$ for $C_{2}$). Combining this with the knowledge that $V$ is
symmetric and the others are antisymmetric in $ij$, the resulting allowed
terms are:

$[A_{1}^{\prime\prime}(\sigma_{x}\tau_{x}+\sigma_{z}\tau_{y})+A_{2}^{\prime\prime}(\sigma_{x}\tau_{y}-\sigma_{z}\tau_{x})]\mu_{y}\omega\\\
+[B_{1}^{\prime\prime}(\sigma_{x}\tau_{x}+\sigma_{z}\tau_{y})+B_{2}^{\prime\prime}(\sigma_{x}\tau_{y}-\sigma_{z}\tau_{x})]\mu_{y}\\\
+[V_{1}^{\prime}(\bar{k}_{x}\sigma_{x}+\bar{k}_{y}\sigma_{z})+V_{2}^{\prime}(\bar{k}_{x}\sigma_{z}-\bar{k}_{y}\sigma_{x})]\tau_{0}\mu_{y}\\\
+[V_{1}(\bar{k}_{x}\sigma_{x}+\bar{k}_{y}\sigma_{z})+V_{2}(\bar{k}_{x}\sigma_{z}-\bar{k}_{y}\sigma_{x})]\tau_{z}\mu_{z}$

Finally, we consider the two-index $N$ matrix, which again is decomposable
into $\alpha_{x}$ and $\alpha_{z}$. These transform like a vector under
$C_{3}$ but a scalar under $C_{2}$, resulting in the allowed terms:

$N=(N_{1}\mu_{0}+N_{2}\mu_{x})(\alpha_{x}\sigma_{x}+\alpha_{z}\sigma_{z})\\\
+(N_{3}\mu_{0}+N_{4}\mu_{x})(\alpha_{x}\sigma_{z}-\alpha_{z}\sigma_{x})$

Therefore, the full self-energy of this system, in the special basis presented
in Appendix A.3.4, is:

$\begin{split}\Sigma_{G}&(\bar{k},\omega)=[(A_{1}\mu_{0}+A_{2}\mu_{x})\sigma_{0}\tau_{0}+(A_{1}^{\prime}\mu_{0}+A_{2}^{\prime}\mu_{x})\sigma_{y}\tau_{z}]\omega\\\
+&(B_{1}\mu_{0}+B_{2}\mu_{x})\sigma_{0}\tau_{0}+(B_{1}^{\prime}\mu_{0}+B_{2}^{\prime}\mu_{x})\sigma_{0}\tau_{z}\\\
+&(C_{1}\mu_{0}+C_{2}\mu_{x})\sigma_{0}(k\cdot\tau)+(D_{1}\mu_{0}+D_{2}\mu_{x})\sigma_{0}(k\times\tau)\\\
+&[A_{1}^{\prime\prime}(\sigma_{x}\tau_{x}+\sigma_{z}\tau_{y})+A_{2}^{\prime\prime}(\sigma_{x}\tau_{y}-\sigma_{z}\tau_{x})]\mu_{y}\omega\\\
+&[B_{1}^{\prime\prime}(\sigma_{x}\tau_{x}+\sigma_{z}\tau_{y})+B_{2}^{\prime\prime}(\sigma_{x}\tau_{y}-\sigma_{z}\tau_{x})]\mu_{y}\\\
+&[V_{1}^{\prime}(\bar{k}_{x}\sigma_{x}+\bar{k}_{y}\sigma_{z})+V_{2}^{\prime}(\bar{k}_{x}\sigma_{z}-\bar{k}_{y}\sigma_{x})]\tau_{0}\mu_{y}\\\
+&[V_{1}(\bar{k}_{x}\sigma_{x}+\bar{k}_{y}\sigma_{z})+V_{2}(\bar{k}_{x}\sigma_{z}-\bar{k}_{y}\sigma_{x})]\tau_{z}\mu_{z}\\\
+&(N_{1}\mu_{0}+N_{2}\mu_{x})[(k_{x}\tau_{y}+k_{y}\tau_{x})\sigma_{x}+(k_{x}\tau_{x}-k_{y}\tau_{y})\sigma_{z}]\\\
+&(N_{3}\mu_{0}+N_{4}\mu_{x})[(k_{x}\tau_{y}+k_{y}\tau_{x})\sigma_{z}-(k_{x}\tau_{x}-k_{y}\tau_{y})\sigma_{x}]\end{split}$
(66)

This 24-coefficient equation is the fully general self-energy for all possible
varieties of coupling. Some terms may vanish to leading order in perturbation
theory.

For the special case of spin-preserving sublattice-independent coupling, which
is the most physically intuitive tunneling term, the only nonvanishing term to
first order (both for lattice mismatch and for twist) is
$A_{0}=\frac{3t_{0}^{2}}{v^{2}Q^{2}}$, and hence (like the half-graphene case)
the practical effect of the lowest-order corrections is a rescaling of the
self-energy.

##  Appendix D Extensions

We now discuss straightforward but nontrivial extensions of the self-energy
computed in Appendix C.1.

###  D.1 $C_{2}\mathcal{T}$ symmetry

In a material that lacks $\mathcal{T}$ symmetry but possesses
$C_{2}\mathcal{T}$ symmetry, it is more convenient to diagonalize in the basis
of $C_{2}\mathcal{T}$ than $\mathcal{T}$. This is also the case in twisted
bilayer graphene, where each valley is invariant under $C_{2}\mathcal{T}$.
Fortunately, the calculation of the self-energy with $C_{2}\mathcal{T}$
symmetry is methodologically similar to the calculation with $\mathcal{T}$
symmetry shown in Sec. C.1, although, in general, more terms appear.

Specifically, the calculation is identical up to Eq. (50), but then instead of
averaging over $T_{Q}$ and $T_{-Q}$ and using $\mathcal{T}$ symmetry to relate
them, we use the fact that
$T_{Q}=\frac{1}{2}(T_{Q}+(C_{2}\mathcal{T})T_{Q}(C_{2}\mathcal{T})^{-1})$, at
which point the calculation proceeds similarly.

If we use Autonne-Takagi factorization on $C_{2}\mathcal{T}$ as discussed in
Sec. A.2, simplifying the operator to $\sigma_{x}\mathcal{K}$, then we instead
find a self-energy of:

$\displaystyle\Sigma_{ij}(\bar{k},\omega)=$
$\displaystyle\sum_{Q,l^{\prime}}T_{Q,il^{\prime}}[\omega-v(\bar{k}+Q)\cdot\sigma]^{-1}T_{Q,jl^{\prime}}^{\dagger}$
(67) $\displaystyle\simeq$
$\displaystyle-\sum_{Q,l^{\prime}}\frac{1}{v^{2}Q^{4}}T_{Q,il^{\prime}}[Q^{2}\omega+v(Q^{2}(\bar{k}\cdot\sigma)-2(Q\cdot\bar{k})(Q\cdot\sigma)+Q^{2}(Q\cdot\sigma))]T_{Q,jl^{\prime}}^{\dagger}$
(68) $\displaystyle=$
$\displaystyle-\sum_{Q,il^{\prime}}\frac{1}{2v^{2}Q^{4}}\left\\{\left(T_{Q,il^{\prime}}T^{\dagger}_{Q,il^{\prime}}+\sigma_{x}(T_{Q,il^{\prime}}T^{\dagger}_{Q,il^{\prime}})^{*}\sigma_{x}\right)Q^{2}\omega+vQ^{2}\left(T_{Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{Q,il^{\prime}}+\sigma_{x}(T_{Q,il^{\prime}}(Q\cdot\sigma)T^{\dagger}_{Q,il^{\prime}})^{*}\sigma_{x}\right)\right.$
$\displaystyle\qquad\left.+v\left(T_{Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{Q,il^{\prime}}+\sigma_{x}\left(T_{Q,il^{\prime}}\left(\bar{k}\cdot(Q^{2}\mathbb{I}-2Q\otimes
Q)\cdot\sigma\right)T^{\dagger}_{Q,il^{\prime}}\right)^{*}\sigma_{x}\right)\right\\}$
(69)

Then the self-energy takes the same form as Eq. (53), except the symmetry
decomposition of the matrices in Eq. (54) are replaced with:

$\displaystyle\mathcal{A}_{ij}=$ $\displaystyle
A_{ij}\sigma_{0}+iA^{\prime}_{ij}\sigma_{z}+A^{\prime\prime}_{ij}\cdot\sigma$
(70a) $\displaystyle\mathcal{B}_{ij}=$ $\displaystyle
B_{ij}\sigma_{0}+iB^{\prime}_{ij}\sigma_{z}+B^{\prime\prime}_{ij}\cdot\sigma$
(70b) $\displaystyle\mathcal{V}_{ij}=$ $\displaystyle
V_{ij}\sigma_{0}+iV^{\prime}_{ij}\sigma_{z}+M_{ij}\cdot\sigma$ (70c)

In this case (contrary to the terms in Eq. (54)), unprimed and double-primed
terms are real symmetric matrices in $ij$, whereas (only) single-primed terms
are antisymmetric in those indices. (Also note the corresponding interchange
of $V$ and $V^{\prime}$, since $C_{2}\mathcal{T}$ imposes different
constraints.)

The equations for these coefficients are generally the same as Eqs. (56)-(59),
except that real and imaginary parts of the trace are changed.

For a single-cone material, the self-energy simplifies to:

$\Sigma^{C_{2}\mathcal{T}}(\bar{k},\omega)\simeq(A_{0}\mathbb{I}+A_{1}\cdot\sigma)\omega\\\
+(B_{0}\mathbb{I}+B_{1}\cdot\sigma)+\bar{k}\cdot
M\cdot\sigma+(\bar{k}\cdot\hat{V})\mathbb{I}$ (71)

We have not yet applied rotational symmetry. In a generic model, we naively
expect more terms allowed by $C_{2}\mathcal{T}$ (without $\mathcal{T}$) than
by $\mathcal{T}$ (without $C_{2}\mathcal{T}$): the $ij$-symmetric (i.e.,
unprimed) parts of (54) allow two scalars, one vector, and a tensor; here we
have three vectors instead of one (and none of the three vectors in the
$C_{2}\mathcal{T}$ care are the same as the vector in the $\mathcal{T}$ case).

However, when rotation symmetry is present, Eq. (71) again simplifies to Eq.
(4). In particular, any rotational symmetry will cause the $A_{1}$, $B_{1}$,
and $V$ terms to vanish (and a rotation of order greater than $2$ will cause
the non-scalar portion of $M$ to vanish), in exactly the same fashion as the
$\mathcal{T}$-based calculation. Then, the remaining coefficients ($A_{0}$,
$B_{0}$, and the remaining parts of $M$) are precisely analogous to the
coefficients in (4). Thus, we find a similar condition on magic angles.

As an example of where this would be useful, consider (a spinless model of)
TBLG. While TBLG can in principle be accounted for by a $\mathcal{T}$-based
transformation, it is much simpler to use this $C_{2}\mathcal{T}$-based result
because $C_{2}\mathcal{T}$ leaves each valley invariant. The ultimate
consequence of this is the the analogy between AB hopping in TBLG and spin-
flipping hopping in our TI-TI model discussed in Sec. III.

###  D.2 More general Dirac Hamiltonians

We now offer a few comments on extending to a broader class of Dirac
Hamiltonians, including Dirac cones not at charge neutrality, gapped Dirac
cones, and anisotropic Dirac cones. In all cases, the calculation is
straightforward to adapt, but it is worth highlighting the extent to which
these alternative Hamiltonians may yield qualitative changes.

For Dirac cones not at charge neutrality, we find that the coefficients $A$
and $B$ mix, i.e., if $A_{0}$ and $B_{0}$ are the coefficients at charge
neutrality, then away from charge neutrality $A$ and $B$ are both linear
combinations of $A_{0}$ and $B_{0}$.

For Dirac cones with mass terms, several new contributions to $V$ and $M$
appear (potentially yielding interesting new physics). To compute the self-
energy, it is most convenient to choose a different basis than the one
outlined in Sec. A: one would rather have the Hamiltonian (with mass term) in
the form $v(k\cdot\sigma)\tau_{0}+m\sigma_{z}\tau_{z}$ with time-reversal
symmetry in the form $\mathcal{T}=i\tau_{x}\sigma_{y}\mathcal{K}$.

This can always be done: performing the sequence of steps through (A.2) puts
the mass term as $im\sigma_{z}\hat{A}_{ij}$, for $\hat{A}$ a real
antisymmetric matrix. This $A_{ij}$ can always be remapped to a block-diagonal
$\tau_{y}$ while preserving the form of time-reversal symmetry (by a Youla
decomposition), at which point we have our Hamiltonian in the form
$v(k\cdot\sigma)\tau_{0}+m\sigma_{z}\tau_{y}$, with
$\mathcal{T}=i\sigma_{y}\tau_{0}\mathcal{K}$. From here, a quarter-rotation
about $\tau_{x}$ will put both matrices in the above-specified form.

Finally, anisotropies also produce additional contributions to $M$ and $V$,
which may lead to different kinds of behavior. For example, anisotropy may
lead to “partial magic angles,” wherein one direction of a Dirac cone has a
vanishing Fermi velocity and the other does not (as discussed in [45, 46]).