# In-plane orbital magnetization as a probe for symmetry breaking in strained
twisted bilayer graphene

Ohad Antebi  Department of Condensed Matter Physics, Weizmann Institute of
Science, Rehovot 76100, Israel  Ady Stern  Department of Condensed Matter
Physics, Weizmann Institute of Science, Rehovot 76100, Israel  Erez Berg
Department of Condensed Matter Physics, Weizmann Institute of Science, Rehovot
76100, Israel

###### Abstract

Three symmetries prevent a twisted bilayer of graphene from developing an in-
plane spontaneous magnetization in the absence of a magnetic field - time
reversal symmetry, $C_{2}$ symmetry to $\pi$ rotation and $C_{3}$ symmetry to
$2\pi/3$ rotation. In contrast, there are experimental and theoretical
indications that, at certain electron densities, time reversal and $C_{2}$
symmetries are broken spontaneously, while the $C_{3}$ symmetry is often
broken due to strain. We show that in-plane orbital magnetization is a very
sensitive probe to the simultaneous breaking of these three symmetries,
exhibiting surprisingly large values (of the order of one Bohr magneton per
moiré unit cell) for valley polarized states at rather small values of
heterostrain. We attribute these large values to the large magnitude of the
characteristic magnetization of individual Bloch states, which we find to
reflect the fast Dirac velocity of single layer graphene, rather than the slow
velocity of the twisted bilayer. These large values are forced to mutually
cancel in valley-symmetric states, but the cancellation does not occur in
valley-polarized states. Our analysis is carried out both for non-interacting
electrons and within a simplified Hartree-Fock framework.

##  I Introduction

Stacking two graphene layers with a relative twist near the “magic
angle”Bistritzer and MacDonald (2011) has produced a flurry of unexpected
physical phenomenaCao  _et al._ (2018a, b); Yankowitz  _et al._ (2019); Arora
_et al._ (2020); Choi  _et al._ (2019); Wong  _et al._ (2020); Saito  _et al._
(2020); Stepanov  _et al._ (2020); Uri  _et al._ (2020); Xie  _et al._ (2019);
Rozen  _et al._ (2021); Sharpe  _et al._ (2019); Park  _et al._ (2021); Jiang
_et al._ (2019); Kerelsky  _et al._ (2019); Lu  _et al._ (2019); Zondiner  _et
al._ (2020); Serlin  _et al._ (2020); Cao  _et al._ (2021); Tschirhart  _et
al._ (2021); Saito  _et al._ (2021); Rozen  _et al._ (2021), including a
variety of symmetry broken phasesJiang  _et al._ (2019); Kerelsky  _et al._
(2019); Lu  _et al._ (2019); Zondiner  _et al._ (2020); Serlin  _et al._
(2020); Cao  _et al._ (2021); Tschirhart  _et al._ (2021); Saito  _et al._
(2021); Rozen  _et al._ (2021), such as spin and valley polarized statesSharpe
_et al._ (2019); Zondiner  _et al._ (2020); Park  _et al._ (2021). An
important question is how are these patterns of broken symmetry manifested in
measurable physical quantities. At first glance, it seems that although an in-
plane magnetic field is a good mean to probe spin polarization, it is an
unlikely candidate for probing _orbital_ response in twisted bilayer graphene
(TBG). This is so because the physical separation between the two layers is of
atomic scale, and at attainable magnetic fields, the flux through the cross
sectional area of a unit cell is a small fraction of a flux quantum.

In this work, we find that, in fact, the single-electron Bloch states of TBG
generically possess a large in-plane orbital magnetization, of the order of
several Bohr magnetons. We trace this observation to the nature of these Bloch
states, which consist of superpositions of components localized at the two
layers. These components move at the Dirac speed of a single graphene layer,
but their directions of motion are nearly opposite, resulting in the slow
group velocity of the narrow TBG bands. This state of affairs is reminiscent
of the familiar “zitterbewegung” phenomenon of relativistic massive Dirac
electrons Breit (1928), which jitter rapidly close to the speed of light even
when their average velocity is non-relativistic. In TBG, the rapid inter-layer
relative motion produces a large in-plane orbital magnetization.

Despite the large magnetic moment of the individual Bloch states, the net in-
plane magnetization vanishes if the system possesses either time reversal or
$C_{3}$ rotational symmetry. However, both of these symmetries are frequently
broken in TBG, either externally or spontaneously. The $C_{3}$ symmetry can be
broken when the system is subjected to strain. Time reversal symmetry can be
broken spontaneously due to electron-electron interactions. Both forms of
symmetry breaking were observed in experiments Choi  _et al._ (2019); Kerelsky
_et al._ (2019); Xie  _et al._ (2019); Kazmierczak and Van Winkle (2021);
Jiang  _et al._ (2019); Sharpe  _et al._ (2019); Serlin  _et al._ (2020);
Tschirhart  _et al._ (2021). Here, we show that in the presence of inter-layer
strain as small as $\sim 10^{-4}$, a spontaneously valley polarized state can
have an in-plane magnetic moment as large as a few $\mu_{B}$ per moiré unit
cell. The in-plane magnetization can be extracted from measurements of the
electronic compressibility in an in-plane magnetic field, by using
thermodynamic relations Zondiner  _et al._ (2020).

Our analysis consists of several parts. In the first part, we compute the in-
plane orbital magnetization of individual Bloch states within the Bistritzer-
MacDonald (BM) model. We then show that in the presence of an inter-layer
strain, the total in-plane moment of each valley is non-zero, and calculate
the moment as a function of the magnitude and direction of the strain.
Finally, we consider the effect of electron-electron interactions within a
simplified Hartree-Fock analysis, and show that if the interactions lead to
the breaking of the valley symmetry, the resulting many-body state has a large
in-plane magnetic moment.

This paper is organized as follows: Section II describes the non-interacting
Hamiltonian in the presence of strain and an in-plane magnetic field. Section
III presents analytical results derived from symmetry considerations and
illustrates them in a toy model. Section IV presents numerical results for the
non-interacting model, and presents a simplified interacting model and
corresponding results. Section V concludes the paper.

##  II Setup

In this section we present a non-interacting model for strained twisted
bilayer graphene with in-plane magnetic fields, which is a generalization of
the BM Hamiltonian. For a pedagogical construction of the effects of strain on
the BM Hamiltonian we refer readers to Refs. [Bi  _et al._ , 2019; Balents,
2019], and give the result here for completeness and in order to set the
notation.

We consider the following continuum Hamiltonian for valley $K$ in the layer
basis:

$H_{\text{BM},K}=\begin{pmatrix}h^{(1)}(\mathbf{k})&T(\mathbf{r})\\\
T^{\dagger}(\mathbf{r})&h^{(2)}(\mathbf{k})\end{pmatrix},$ (1)

where $h^{(1)},h^{(2)}$ are the single layer Hamiltonians of the deformed
graphene sheets near valley $K$, and $T(\mathbf{r})$ is the interlayer
coupling term (see below). The opposite valley’s Hamiltonian
$H_{\text{BM},K^{\prime}}$ is obtained by time-reversal, and is assumed to be
decoupled from $H_{\text{BM},K}$. The momentum space representation of the
continuum intralayer Hamiltonian in layer $l$ is given to linear order in the
(unsymmetrized) strain tensor as follows:

$h^{(l)}(\mathbf{k})=\hbar
v_{F}\left[(\mathbb{1}+\mathcal{U}^{(l)})^{T}(\mathbf{k}-\delta\mathbf{K}^{(l)})\right]\cdot\begin{pmatrix}\sigma_{x},\sigma_{y}\end{pmatrix}.$
(2)

Here $v_{F}=0.9\times 10^{6}\frac{m}{s}$ is the bare Dirac velocity,
$\mathbf{k}$ is measured with respect to the $K$ point of the undeformed
graphene layer. $\delta\mathbf{K}^{(l)}$ is the shift in position of the Dirac
node of layer $l$ from its undeformed position and will be elaborated upon
shortly. The matrix $\mathcal{U}^{(l)}$ is the unsymmetrized strain tensor
whose components are $u_{ij}^{(l)}=\partial_{i}u_{j}^{(l)}$, with
$u_{j}^{(l)}(\mathbf{r})$ being the displacement field along the $j$ axis. The
antisymmetric and symmetric components of this strain tensor will be hereby
denoted as:

$\mathcal{U}^{(l)}=\begin{pmatrix}s_{xx}^{(l)}&s_{xy}^{(l)}-\theta^{(l)}\\\
s_{xy}^{(l)}+\theta^{(l)}&s_{yy}^{(l)}\end{pmatrix}=\mathcal{R}(\theta^{(l)})+\mathcal{S}^{(l)},$
(3)

where the antisymmetric part $\mathcal{R}(\theta^{(l)})$ describes a rotation
of the layer, and the symmetric part $\mathcal{S}^{(l)}$ includes the
traceless strain that breaks $C_{3}$ symmetry. We twist the layer
symmetrically by choosing $\theta^{(1)}=-\theta^{(2)}=\theta/2$. The
discussion of the strain tensor is restricted to systems with
$\mathcal{S}^{(1)}=\chi\mathcal{S}^{(2)}$ where $\chi=-1$ ($\chi=+1$)
describes heterostrain (homostrain). It is assumed throughout that the strain
is homogeneous and weak compared to the twist angle, i.e. $s_{ij}\ll\theta\ll
1$.

The position of the Dirac node in a deformed graphene monolayer is determined
by a combination of geometric and intra-layer energetics considerations.
Geometric considerations cause the original Dirac node to be rotated by the
twist and then shifted by the strain, summarized as
$\mathbf{K}\rightarrow(1-\mathcal{U}^{(l)})^{T}\mathbf{K}$. Here
$\mathbf{K}=\left(\frac{4\pi}{3\sqrt{3}a},0\right)$ is the momentum of the
original Dirac node, and $a=0.142$nm is the distance between carbon atoms. By
accounting for the effect of strain on the intralayer hopping, the Dirac node
is shifted further, an effect that can be described using an artificial gauge
field Suzuura and Ando (2002); Castro Neto  _et al._ (2009):

$\mathbf{A}_{s}^{(l)}=\frac{\beta}{2a}\begin{pmatrix}s_{xx}^{(l)}-s_{yy}^{(l)},-2s_{xy}^{(l)}\end{pmatrix},$
(4)

with $\beta\approx 3.14$ being a dimensionless parameter Nam and Koshino
(2017). Note that the gauge field couples oppositely to the two valleys, and
hence does not break time-reversal symmetry. The combined shift in position of
the Dirac node due to twist and strain, $\delta\mathbf{K}^{(l)}$, is therefore
given by:

$\displaystyle\begin{split}\delta\mathbf{K}^{(l)}&=-(\mathcal{U}^{(l)})^{T}\mathbf{K}-\mathbf{A}_{s}^{(l)}\\\
&=(\mathcal{R}(\theta^{(l)})-\mathcal{S}^{(l)})\mathbf{K}-\mathbf{A}_{s}^{(l)}.\end{split}$
(5)

Note that, in the absence of strain, this shift is simply
$\delta\mathbf{K}^{(l)}=-\theta^{(l)}\hat{z}\times\mathbf{K}$.

The tunneling amplitude between the layers will have the periodicity of the
strained moiré pattern. We allow for different tunneling amplitudes between AA
and AB sites, denoted by $t_{AA}$ and $t_{AB}$ respectively. Accounting for
lattice relaxation Koshino  _et al._ (2018); Carr  _et al._ (2019), we set
$t_{AA}=0.7t_{AB}$ with $t_{AB}=0.11\,\rm{eV}$. The tunneling in the
sublattice basis can be written as:

$T(\mathbf{r})=\sum_{j=0,1,2}T_{j}e^{-i\mathbf{g}_{j}\cdot\mathbf{r}},$ (6)

where
$T_{j}=t_{AA}+t_{AB}\left[\textrm{Re}{(\zeta^{j})}\sigma_{x}-\textrm{Im}{(\zeta^{j})}\sigma_{y}\right]$,
and $\zeta=e^{i2\pi/3}$. The moiré reciprocal lattice vectors are given by
$\mathbf{g}_{j}=\mathcal{U}^{T}\mathbf{G}_{j}$ where
$\mathcal{U}=\mathcal{U}^{(1)}-\mathcal{U}^{(2)}$ is the relative twist and
strain, and
$\mathbf{G}_{0}=0,\mathbf{G}_{1}=\sqrt{3}K\left(\frac{\sqrt{3}}{2},\frac{1}{2}\right),\mathbf{G}_{2}=\sqrt{3}K\left(\frac{\sqrt{3}}{2},-\frac{1}{2}\right)$
are the graphene reciprocal lattice vectors. While the Dirac node moves for
either homostrain or heterostrain, to leading order in the strain and twist
angle only heterostrain modifies the moiré pattern and hence the tunneling
between the layers.

Next, we consider the application of an in-plane magnetic field. A uniform in-
plane field $\mathbf{B}_{\parallel}=|B|(\cos{\phi_{B}},\sin{\phi_{B}})$ can be
described by a $z$-dependent, in-plane vector potential
$\mathbf{A}(z)=z|B|(\sin{\phi_{B}},-\cos{\phi_{B}})$. We choose $z=0$ to be
halfway between the layers, which are separated by distance
$d=0.335\,\rm{nm}$, leading to constant but opposite gauge fields for
different layers. Defining

$\mathbf{A}_{\text{em}}\equiv\mathbf{A}\left(z=\frac{d}{2}\right)=\frac{1}{2}|B|d(\sin{\phi_{B}},-\cos{\phi_{B}}),$
(7)

we can minimally couple the in-plane magnetic field to the TBG as follows:

$H_{\text{BM},K}=\begin{pmatrix}h^{(1)}\left(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)&T(\mathbf{r})\\\
T^{\dagger}(\mathbf{r})&h^{(2)}\left(\mathbf{k}+\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)\end{pmatrix}-\bm{\mu}_{s}\cdot\mathbf{B}_{\parallel},$
(8)

where $\bm{\mu}_{s}$ is the electron’s spin magnetic moment. Note that the
local tunneling amplitude is unaffected by the in-plane field. This completes
our non-interacting description of the strained TBG with in-plane magnetic
field.

##  III General considerations

###  III.1 Symmetry analysis

The zero-temperature, zero-field in-plane magnetization per moiré unit cell is
defined as:

$\mathbf{M}_{\parallel}=-\left.\frac{1}{N_{u.c.}}\frac{\partial
E_{G.S.}}{\partial\mathbf{B}_{\parallel}}\right|_{\mathbf{B}_{\parallel}=0}$
(9)

where $E_{G.S.}=\bra{\Psi_{G.S.}}H\ket{\Psi_{G.S.}}$ is the energy of the
ground state $\ket{\Psi_{G.S.}}$, and $N_{u.c.}$ is the number of moiré unit
cells. We are interested in the restrictions symmetry imposes on this
magnetization. We will show that both valley polarization and traceless
heterostrain are necessary for a non-vanishing orbital magnetization.
Furthermore, to linear order in strain, the magnitude of the magnetization is
independent of the strain direction.

We start by listing the symmetries of the unstrained TBG, in the absence of
any magnetic field. The symmetries of a single unstrained graphene layer are
given by two and three-fold rotations about the $z$ axis, $C_{2}$, $C_{3}$,
two mirror symmetries $M_{y,z}$ about the $x-z$ and $x-y$ planes,
respectively, time-reversal $\mathcal{T}$, and $SU(2)$ spin rotation symmetry.
The twisted bilayer as a whole breaks $M_{y}$ and $M_{z}$, but has two
additional symmetries given by $R_{\pi}^{x}$ and $R_{\pi}^{y}$, which are
$\pi$-rotations about the $x$ and $y$ axes, respectively, in the plane halfway
between the layers (note that these symmetries interchange the two layers).

In order to avoid an overly complicated description, we restrict our attention
to a subset of strains and flavor symmetry breaking configurations of
interest. We restrict the discussion of strain to that of either heterostrain
or homostrain, i.e. $\mathcal{S}^{(1)}=\chi\mathcal{S}^{(2)}\equiv\mathcal{S}$
with $\chi=-1$ or $\chi=+1$, respectively. For the flavor symmetry, we
anticipate any $C_{2}$ or $\mathcal{T}$ symmetric state to yield
$\mathbf{M}_{\parallel}=0$. To highlight orbital magnetization we consider
variable fillings of the different valleys, $\\{\nu_{K},\nu_{K^{\prime}}\\}$,
with
$\nu_{K(K^{\prime})}=\nu_{K(K^{\prime})\uparrow}+\nu_{K(K^{\prime})\downarrow}$
(where $\uparrow$, $\downarrow$ indicate the two spin directions). Equipped
with these parameters, we write the ground state energy as
$E\left(\mathcal{S},\mathbf{B}_{\parallel},\nu_{K},\nu_{K^{\prime}}\right)$,
and explore how the different symmetries constrain this energy.

The $C_{2}$ symmetry constrains the energy to the form:

$E\left(\mathcal{S},\mathbf{B}_{\parallel},\nu_{K},\nu_{K^{\prime}}\right)=E\left(\mathcal{S},-\mathbf{B}_{\parallel},\nu_{K^{\prime}},\nu_{K}\right).$
(10)

We note that this symmetry flips the valleys and the magnetic field, but does
not affect the strain tensor. Time reversal symmetry imposes the same
constraint. By examining (9) and (10) it is clear that finite in-plane
magnetization requires flavor symmetry breaking in the form of valley
polarization, i.e. $\nu_{K}\neq\nu_{K^{\prime}}$.

Next, we consider the effect of $C_{3}$ rotation. This symmetry preserves the
valley fillings, but rotates both magnetic field and strain:

$E\left(\mathcal{S},\mathbf{B}_{\parallel},\nu_{K},\nu_{K^{\prime}}\right)=E\left(C_{3}\mathcal{S}C^{T}_{3},C_{3}\mathbf{B}_{\parallel},\nu_{K},\nu_{K^{\prime}}\right).$
(11)

with $C_{3}=\exp{\left(i2\pi\eta_{y}/3\right)}$ and $\eta_{y}$ is the $y$
Pauli matrix. Expanding (11) in components of $\mathcal{S}$ and
$\mathbf{B}_{\parallel}$ introduces different constraints on different powers
of strain and magnetic field. The zero-field magnetization is given by the
coefficient of the first order in magnetic field. By considering the
$\mathcal{S}=0$ case, we immediately see that zero field in-plane
magnetization is prohibited for an unstrained TBG, regardless of flavor
symmetry breaking, due to the $C_{3}$ symmetry.

We wish to construct a magnetization vector that is linear in strain elements
and rotates under $C_{3}$ as dictated by Eq. (11). When considering rotations
of a strain tensor it is most natural to think of rotating its principal axis.
We therefore continue with the diagonalized form of the strain tensor:

$\mathcal{S}=\begin{pmatrix}\cos{\phi_{S}}&-\sin{\phi_{S}}\\\
\sin{\phi_{S}}&\cos{\phi_{S}}\end{pmatrix}\begin{pmatrix}s^{\prime}&0\\\
0&s^{\prime\prime}\end{pmatrix}\begin{pmatrix}\cos{\phi_{S}}&\sin{\phi_{S}}\\\
-\sin{\phi_{S}}&\cos{\phi_{S}}\end{pmatrix}.$ (12)

where $\phi_{S}$ is the angle between the principal axis and the $x$ axis, and
$s^{\prime},s^{\prime\prime}$ denote the principal strains. The $C_{3}$
rotation of the strain tensor only transforms
$\phi_{S}\rightarrow\phi_{S}+2\pi/3$, leaving the principal strains invariant.

Any rotationally invariant component of the strain tensor would not contribute
to the magnetization vector at linear order. The strain tensor in (12) can be
split into the dilatation and the traceless part as

$\mathcal{S}=\frac{s^{\prime}+s^{\prime\prime}}{2}\mathbb{1}+\frac{s^{\prime}-s^{\prime\prime}}{2}(\eta_{x}\sin
2\phi_{S}+\eta_{z}\cos 2\phi_{S}),$ (13)

where we denote the identity matrix by $\mathbb{1}$, and $\eta_{x,z}$ are
Pauli matrices. Clearly, the dilatation is rotationally invariant, and
therefore does not contribute to the magnetization at linear order, leaving
only the traceless part at this order. We conclude that the magnetization is
proportional to $s^{\prime}-s^{\prime\prime}$, the difference between
principal strains. Since the angle $2\phi_{S}$ transforms under $C_{3}$
rotation as $2\phi_{S}\rightarrow 2\phi_{S}+4\pi/3$, while the magnetization
rotates by $2\pi/3$, we find that the following form of the energy satisfies
Eq. (11):

$\displaystyle\begin{split}&E\left(\mathcal{S},\mathbf{B},\nu_{K},\nu_{K^{\prime}}\right)\cong-\mathbf{M}_{\parallel}\cdot\mathbf{B}_{\parallel}\\\
&\cong\lambda(\nu_{K},\nu_{K^{\prime}})|B|(s^{\prime}-s^{\prime\prime})\cos{(2\phi_{S}+\phi_{B}+\phi_{0})},\end{split}$
(14)

where $\lambda(\nu_{K},\nu_{K^{\prime}})=-\lambda(\nu_{K^{\prime}},\nu_{K})$
is an odd function under exchange of its arguments, $\phi_{0}\in[0,\pi)$ is a
constant angle soon to be determined, and
$\mathbf{B}_{\parallel}=|B|\left(\cos{\phi_{B}},\sin{\phi_{B}}\right)$. Note
that the magnitude of the magnetization is independent of the strain’s
principal axis at linear order. This observation does not hold to higher
orders in strain, as shown in detail in Appendix B.

Eq. (14) clearly shows how the magnetization respects the $C_{3}$ symmetry. We
still, however, have a degree of freedom in choosing $\phi_{0}$. This phase
will turn out to depend on whether the system is homostrained or
heterostrained, and it is determined by symmetries that relate the two
graphene layers.

Using the $R_{\pi}^{x}$ symmetry which preserves the valley fillings, we find
the energy must be symmetric under the transformation
$\mathcal{S}\rightarrow\chi\mathcal{S}$ (due to layer flip),
$\phi_{B}\rightarrow-\phi_{B}$ and $\phi_{S}\rightarrow-\phi_{S}$. For
heterostrain ($\chi=-1$), this constraint enforces $\phi_{0}=\pi/2$, yielding
the following form for the magnetization energy:

$\mathbf{M}_{\parallel}\cdot\mathbf{B}_{\parallel}\cong\lambda(\nu_{K},\nu_{K^{\prime}})|B|(s^{\prime}-s^{\prime\prime})\sin{(2\phi_{S}+\phi_{B})}.$
(15)

For the homostrained case we obtain $\phi_{0}=0$ and one can replace the sine
above by a cosine. It can be verified that the symmetry $R_{\pi}^{y}$ does not
add any other constraint (see Appendix B). The magnetization immediately
follows, and is given here for completeness:

$\mathbf{M}_{\parallel,x}+i\mathbf{M}_{\parallel,y}\cong\lambda(\nu_{K},\nu_{K^{\prime}})(s^{\prime}-s^{\prime\prime})e^{-i\left(2\phi_{S}+\phi_{0}\right)},$
(16)

with $\phi_{0}=\pi/2$ ($\phi_{0}=0$) for heterostrain (homostrain),
respectively. Motivated by the result above, we will proceed by considering
only traceless strain tensors $\mathcal{S}$ with
$s^{\prime}=-s^{\prime\prime}$ such that the leading order contribution does
not vanish.

###  III.2 Qualitative discussion of orbital magnetization

The in-plane orbital magnetic moment of a Bloch state can be pictorially
described as resulting from a current loop consisting of components in
opposite layers, moving in opposite directions close to the single-layer Dirac
velocity $v_{F}$. Close to the magic angle, the net group velocity of the
Bloch state is much smaller than $v_{F}$. We will show, however, that the
flatness of the resulting bands is still important for the total magnetization
in the presence of strain. As shown in Refs. Huder  _et al._ , 2018; Qiao  _et
al._ , 2018; Bi  _et al._ , 2019; Liu  _et al._ , 2021; Parker  _et al._ ,
2021, the flat bands are extremely susceptible to strain, thus allowing for a
significant rotation symmetry breaking near the magic angle.

With this picture in mind, we estimate the magnetic moment of a Bloch state by
constructing a current loop from the physical parameters of the bilayer. The
current in the loop is estimated to be $I\sim ev_{F}/L$, where $L$ is a length
scale that will soon cancel out. The loop area is given by the same length
scale $L$ times the distance between the layers $d$. We therefore find that
orbital magnetic moment associated with such a loop is given by
$|M_{\parallel}|\sim ev_{F}d\sim 5\mu_{B}$, where $\mu_{B}$ is the Bohr
magneton. We show below that this overestimates the magnetization by a factor
of $2-3$.

The simplest model that allows for a quantitative analytical calculation of
the magnetization is the “Tripod model”Bistritzer and MacDonald (2011);
Bernevig  _et al._ (2021). Within this model, the effective Hamiltonian with
applied heterostrain and in-plane magnetic field is given by:

$\displaystyle\begin{split}H(\mathbf{k})&=h^{(1)}\left(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)\\\
&-\sum_{j=0,1,2}T_{j}\left[h^{(2)}\left(\mathbf{k}+\frac{e}{\hbar}\mathbf{A}_{\text{em}}+\mathbf{g}_{j}\right)-E(\mathbf{k})\right]^{-1}T_{j}^{\dagger},\end{split}$
(17)

where the energy is to be solved self-consistently using
$H(\mathbf{k})\ket{\Psi(\mathbf{k})}=E(\mathbf{k})\ket{\Psi(\mathbf{k})}$.
This model is valid in the vicinity of layer $1$’s Dirac point, i.e. where
$\mathbf{p}=\mathbf{k}-\delta\mathbf{K}^{(1)}$ is small compared to the moiré
BZ. We expand this model to first order in $p/(K\theta)$ and in strain in
Appendix A. Omitting a $\mathbf{p}$-independent constant, the result for
valley $\tau=K$ is given by:

$H(\delta\mathbf{K}^{(1)}+\mathbf{p})\cong\hbar
v_{F}\left(\frac{v_{F}^{*}}{v_{F}}(\mathbf{p}-\mathbf{p}_{0})-2\frac{v_{F,A}}{v_{F}}\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)\cdot\mathbf{\sigma},$
(18)

with the renormalized velocities given by:

$\displaystyle\begin{split}\frac{v_{F}^{*}}{v_{F}}&=\frac{1-3w^{2}}{1+3(u^{2}+w^{2})},\\\
\frac{v_{F,A}}{v_{F}}&=\frac{1}{1+3(u^{2}+w^{2})},\end{split}$ (19)

where $u=t_{AA}/(\hbar v_{F}K\theta)$ and $w=t_{AB}/(\hbar v_{F}K\theta)$ are
the dimensionless tunneling parameters.

The tripod model predicts a shift of the Dirac points in the presence of
strain due to the breaking of $C_{3}$ symmetry, given by:

$\mathbf{p}_{0}=\frac{3w^{2}}{1-3w^{2}}\left(1+\frac{u^{2}}{w^{2}}\frac{Ka}{\beta}\right)\mathbf{A}_{s},$
(20)

where $\mathbf{A}_{s}=\mathbf{A}_{s}^{(1)}-\mathbf{A}_{s}^{(2)}$, as given by
(4), and is proportional to the strain. Note that for homostrain
$\mathbf{A}_{s}$ vanishes, and for heterostrain it does not. The reason is
that, to first order in strain, only the strain difference between layers
modifies the moiré pattern.

We define the magic angle condition to be $v_{F}^{*}=0$, or equivalently
$w^{2}=1/3$. By examining Eq. (20) we find, for a heterostrained TBG, an acute
sensitivity of the Dirac point position to the twist angle in the vicinity of
the magic angle condition. A similar strong sensitivity of the Dirac points to
twist angle in the presence of an in-plane magnetic field was discussed in
Ref. Kwan  _et al._ (2020). The acute sensitivity is manifest also in the full
numerical calculation which extends beyond the tripod model. The large shift
of the Dirac point indicates a strong breaking of $C_{3}$ symmetry, and
therefore allows for a substantial net in-plane orbital magnetization per
moiré unit cell, as section IV will show.

Consider the magnetization of a particular momentum state:

$\mathbf{M}_{\parallel}(\mathbf{p})=-\left.\frac{\partial
E_{-}(\mathbf{p})}{\partial\mathbf{B}_{\parallel}}\right|_{\mathbf{B}_{\parallel}=0}$
(21)

where $E_{-}(\mathbf{p})$ is the lower energy band. The in-plane magnetization
is readily calculated from Eq. (18) to be

$\mathbf{M}_{\parallel}(\mathbf{p})=ev_{F,A}d\frac{\hat{z}\times(\mathbf{p}-\mathbf{p}_{0})}{\left|\mathbf{p}-\mathbf{p}_{0}\right|}.$
(22)

This result is in agreement with our initial discussion of the expected
magnetization magnitude, up to a renormalization of the velocity due to
tunneling. The beauty of this result is that, indeed, the renormalized
velocity $v_{F,A}$ that enters the magnetization does not vanish for any twist
angle. At the magic angle, our new estimate for the magnetization magnitude is
given by $|M_{\parallel}|=ev_{F,A}d\approx 2\mu_{B}$. We will see that this
result is a fair estimate for the total magnetization per moiré unit cell when
the $C_{3}$ rotation and time-reversal $\mathcal{T}$ symmetries are strongly
broken.

###  III.3 Effects of interactions

The vanishing Dirac velocity at the magic angle is a property of the non-
interacting band structure. By including electron-electron interactions, one
expects the Dirac velocity and the bandwidth to modified. Therefore, a natural
question to ask is - How do interactions affect the strain sensitivity of the
$C_{3}$ symmetry breaking near the magic angle? We will demonstrate in section
IV that, within our simplified Hartree-Fock treatment, the necessary strain
for a significant total magnetization remains within reasonable experimental
values. We will also show that the magnetization magnitudes of individual
Bloch states are practically unchanged by the interactions, due to the
insensitivity to interactions of the parameters that determine it, $v_{F,A}$
and $d$.

##  IV Numerical analysis

###  IV.1 Valley resolved magnetization in a non-interacting model

The non-interacting ground state, even in the presence of strain, has a zero
overall orbital magnetization due to time-reversal symmetry. However, it is
still instructive to calculate the contribution of a single valley to the net
orbital magnetization. This calculation will become useful when we discuss
states which are valley polarized due to interactions.

In the non-interacting model, the ground state energy is the sum of energies
of single particle Bloch states,
$E_{G.S.}=\sum_{\alpha,\mathbf{k}}^{\prime}E_{\alpha}(\mathbf{k})$, where
$\alpha=(\tau,s,m)$ is a combined label for valley $\tau$, spin $s$ and band
index $m$, and the prime restricts the summation to occupied states.
Similarly, the total magnetization in Eq. (9) is given by an average over the
magnetizations of individual Bloch states:
$\mathbf{M}_{\parallel}=\sum_{\alpha,\mathbf{k}}^{\prime}\mathbf{m}_{\alpha}(\mathbf{k})/N_{u.c.}$.

The contribution of a single band and flavor combination $\alpha$ to the
magnetization can be visualized by plotting the orbital component (ignoring
the trivial contribution of spin moments) of $\mathbf{m}_{\alpha}(\mathbf{k})$
in the entire moiré Brillouin zone. Focusing on the valence band (the lower
flat band) of valley $K$, we plot in Fig. 1 the magnetization for strain-free
and heterostrained MATBG. We find that the orbital magnetization of each Bloch
state is of order $\sim 1\mu_{B}$ in both cases, and the contribution of the
entire band to the magnetization is determined by the strength of the $C_{3}$
symmetry breaking due to the strain.

Figure 1: Orbital in-plane magnetization vectors of individual Bloch states in
the valence band of valley $\tau=K$, plotted across the moiré BZ for the
noninteracting TBG with twist angle of $\theta=1.1^{\circ}$. Background
contours represent the energy band. a) No strain, b) Heterostrain with
principal strains $s^{\prime}=-s^{\prime\prime}=2\times 10^{-4}$ and angle
$\phi_{S}=\pi/4$. Notice how the magnetization of each state is $\sim
1\mu_{B}$ for both cases, but in the unstrained case the total magnetization
vanishes due to $C_{3}$ symmetry. The hexagon outlines the moiré BZ, which is
effectively unchanged by the strain for $s^{\prime}/\theta\sim 10^{-2}$.

Due to the qualitative considerations of Sec. III.2, we expect the strength of
the symmetry breaking to be enhanced in the vicinity of a magic angle. We
confirm this enhancement by computing the valley resolved magnetization for
different twist angles as a function of filling. We will restrict our sum over
bands to the two flat bands of a single spin direction and valley $\tau=K$.
The filling factor of this flavor with two bands will be denoted by
$\nu_{K}\in[-1,1]$. The magnetization vs. filling is plotted for a typical
magnitude of strain in Fig. 2 for five different twist angles,
$\theta\in\\{0.95^{\circ},1.07^{\circ},1.1^{\circ},1.3^{\circ},1.5^{\circ}\\}$.
With our choice of model parameters, the first magic angle is at
$\theta_{MA}\approx 1.07^{\circ}$, and the flat bands are separated by a gap
from the remote bands for $\theta\gtrsim 0.9^{\circ}$. Fig. 2 shows three key
points: $(i)$ It is possible to obtain a significant in-plane orbital
magnetization in the non-interacting model using strain, provided that flavor
symmetry is broken, $(ii)$ the magnetization is highly sensitive to the twist
angle, peaking in the vicinity of the magic angle, and $(iii)$ the conduction
band exactly cancels the valence band contribution. This final point implies
that filled flavors do not contribute to an overall magnetization, and we find
that a significant response occurs near half filling of a single flavor.

Figure 2: The valley-$K$ resolved magnitude of the magnetization vs filling
for the noninteracting TBG with various twist angles. The filling $\nu_{K}$
denotes the number of electrons per moiré unit cell with respect to charge
neutrality. The strain parameters are specified in Fig. 1b. Note the strong
magnetization near the magic angle $\theta_{MA}\approx 1.07^{\circ}$, and that
the magnetization of the valence band is exactly cancelled by that of the
conduction band.

Near the magic angle, the magnetization peaks when completely filling the
valence band ($\nu_{K}=0$). Therefore, we plot the magnetization vs. strain at
this filling in Fig. 3 for several twist angles. The key result from this
figure is that, near the magic angle, the strain magnitude necessary to reach
$M_{\parallel}\sim 1\mu_{B}$ is extremely small, even in comparison with the
twist angle, which sets a reference scale of $\theta\sim 0.02$. For small
strain one observes (inset) the anticipated linear dependence on strain
magnitude. The linear response is associated with a magnetization magnitude
that is independent of the strain direction (Appendix B). At higher strain
values, the magnitude of the magnetization becomes generally dependent on
strain direction. Extreme strains approaching the twist angle are beyond the
scope of this paper. In particular, for $s^{\prime}=-s^{\prime\prime}=s_{0}$
and $\phi_{S}=45^{\circ}$, the moiré pattern is a square lattice for
$s_{0}/\theta=2-\sqrt{3}\approx 0.268$. This motivates us to consider
$s_{0}/\theta\lesssim 0.1$ as a bound for the strain under which the system is
dominantly “twisted”.

Figure 3: The valley-$K$ resolved magnitude of the magnetization vs strain at
$\nu_{K}=0$ for the noninteracting TBG. The horizontal axis $s_{0}$ denotes
the principal strains ($s^{\prime}=-s^{\prime\prime}=s_{0}$) of the
heterostrain with $\phi_{S}=45^{\circ}$, and is bounded by
$s_{0}/\theta\lesssim 0.1$ such that the system is dominantly “twisted”. Note
that the strain necessary to obtain a significant magnetization is extremely
small near the magic angle. Inset presents the data on a log-log plot, and
shows the linear dependence at small strain.

###  IV.2 Magnetization of valley polarized states in the presence of
interactions

We now turn our attention to the effect of electron-electron interactions. As
previously discussed, $C_{3}$ symmetry breaking is insufficient by itself for
producing a net orbital magnetization - one must also break $C_{2}$ and
$\mathcal{T}$ symmetries. In MATBG, these two symmetries may be spontaneously
broken due to electron-electron interactions, which possibly lead to ground
states that are valley polarized, spin polarized or both. In this section, we
use a simplified Hartree-Fock procedure to calculate the orbital magnetization
of valley polarized states, including effects of band renormalizations due to
interactions.

| Initial: $\overrightarrow{\nu_{i}}$ | Final: $\overrightarrow{\nu_{f}}$
---|---|---
Configuration | $\nu_{K\uparrow}^{i}$ | $\nu_{K^{\prime}\downarrow}^{i}$ | $\nu_{K\downarrow}^{i}$ | $\nu_{K^{\prime}\uparrow}^{i}$ | $\nu_{K\uparrow}^{f}$ | $\nu_{K^{\prime}\downarrow}^{f}$ | $\nu_{K\downarrow}^{f}$ | $\nu_{K^{\prime}\uparrow}^{f}$
Unpolarized | $-1$ | $-1$ | $-1$ | $-1$ | $+1$ | $+1$ | $+1$ | $+1$
VP $\nu=+3$ | $-1$ | $+1$ | $+1$ | $+1$ | $+1$ | $+1$ | $+1$ | $+1$
VP $\nu=+1$ | $-1$ | $+1$ | $-1$ | $-1$ | $+1$ | $+1$ | $+1$ | $+1$
VP $\nu=-1$ | $-1$ | $-1$ | $-1$ | $-1$ | $+1$ | $-1$ | $+1$ | $+1$
VP $\nu=-3$ | $-1$ | $-1$ | $-1$ | $-1$ | $+1$ | $-1$ | $-1$ | $-1$
Table 1: Flavor interpolation scheme for unpolarized and valley-polarized (VP)
states. A completely empty (filled) flavor $\\{\tau,s\\}$ is represented by a
filling factor of $\nu_{\tau,s}^{i/f}=-1$ ($\nu_{\tau,s}^{i/f}=+1$). The
density matrix is constructed for the two limits labeled here by “Initial” and
“Final” using the non-interacting states, and then interpolated for
intermediate fillings. The VP states with $\nu\in\\{-3,-1,+1,+3\\}$ are
obtained by the halfway point between the initial and final points defined by
$\overrightarrow{\nu_{i}}$ and $\overrightarrow{\nu_{f}}$.

####  IV.2.1 Calculation method for the ground state energy and magnetization

The orbital magnetization in Eq. (9) requires us to evaluate the gradient of
the ground state energy with respect to magnetic field. A standard
approximation procedure for obtaining the ground state of an interacting
Hamiltonian $\hat{H}$ involves a self-consistent Hartree-Fock (SCHF)
procedure. In this procedure one minimizes the expectation energy
$E=\langle\hat{H}\rangle$ over a family of Slater determinant wavefunctions by
self-consistently solving for the ground state of the single-particle Hartree-
Fock Hamiltonian:

$\hat{H}_{HF}=\hat{H}_{0}+\hat{\Sigma}^{HF}[\rho],$ (23)

where $\hat{H}_{0}$ is a non-interacting Hamiltonian,
$\hat{\Sigma}^{HF}[\rho]$ is the self energy operator in the Hartree-Fock
approximation, and $\rho$ is the single-particle density matrix corresponding
to the ground state. The explicit form of these operators, including the
scheme to subtract the interaction effects of a pair of decoupled layersXie
and Macdonald (2020, 2021), is given in Appendix C.

The SCHF procedure performed at filling $\nu$ and magnetic field
$\mathbf{B}_{\parallel}$ provides the chemical potential
$\mu(\nu,\mathbf{B}_{\parallel})$ as the lowest quasiparticle excitation
energy. The ground state energy can be obtained by integrating the chemical
potential over the filling, and the orbital magnetization in Eq. (9) can be
calculated by:

$\mathbf{M}_{\parallel}(\nu)=-\int_{\nu_{0}}^{\nu}\left.\frac{\partial\mu(\nu^{\prime},\mathbf{B}_{\parallel})}{\partial\mathbf{B}_{\parallel}}\right|_{\mathbf{B}_{\parallel}=0}d\nu^{\prime}+\mathbf{M}_{\parallel}(\nu_{0}).$
(24)

The self-consistent approach for determining $\rho(\nu)$ in Eq. (23) is
computationally expensive. Instead, we employ an interpolation scheme similar
to the work by Xie and MacDonald (Ref. Xie and Macdonald, 2021) for
approximating the density matrix for arbitrary filling without performing a
self-consistent calculation. In this method, we take two single-particle
density matrices at two filling configurations, denoted by $\rho(\nu_{i})$ and
$\rho(\nu_{f})$, and assume that the density matrix $\rho(\nu)$ at any
intermediate filling is a linear interpolation between the two limits.
Neglecting any interaction-driven hybridization between remote and flat bands,
the many-electron ground state at $\nu=\pm 4$ becomes the non-interacting
single-particle ground state, thereby providing exact density matrices for
$\nu_{i}=-4$ and $\nu_{f}=+4$.

We are interested in valley-polarized (VP) states, which cannot be adequately
described by filling all flavors equally. Therefore we consider individual
flavor filling components
$\overrightarrow{\nu}=(\nu_{K\uparrow},\nu_{K^{\prime}\downarrow},\nu_{K\downarrow},\nu_{K^{\prime}\uparrow})$,
such that the total filling $\nu$ is the sum of the four components of
$\overrightarrow{\nu}$. We perform a linear interpolation of the density
matrix in this four dimensional parameter space between two configuration
points $\overrightarrow{\nu_{i}}$ and $\overrightarrow{\nu_{f}}$:

$\rho(\nu)=\frac{\nu_{f}-\nu}{\nu_{f}-\nu_{i}}\rho(\overrightarrow{\nu_{i}})+\frac{\nu-\nu_{i}}{\nu_{f}-\nu_{i}}\rho(\overrightarrow{\nu_{f}}),\quad\nu\in[\nu_{i},\nu_{f}],$
(25)

where $\rho(\overrightarrow{\nu_{i}})$ and $\rho(\overrightarrow{\nu_{f}})$
are diagonal in flavor, with an overall filling factor $\nu_{i}$ and
$\nu_{f}$, respectively. If the individual flavors in
$\overrightarrow{\nu_{i}}$ (or $\overrightarrow{\nu_{f}}$) are either empty or
full, each flavor block in the density matrix $\rho(\overrightarrow{\nu_{i}})$
(or $\rho(\overrightarrow{\nu_{f}})$) is given by the corresponding density
matrix of the non-interacting BM model.

The choice of initial and final density matrices for interpolation determines
the nature of the ground state under consideration. Consider the VP state at
$\nu=+1$ which is characterized by $\overrightarrow{\nu}=(0,+1,0,0)$, and
consists of one filled flavor and three half filled flavors. Within our
scheme, the density matrix of such a state can be obtained by interpolating
between a $\nu_{i}=-2$ initial state in which three flavors are empty and one
is full, $\overrightarrow{\nu_{i}}=(-1,+1,-1,-1)$, and a $\nu_{f}=+4$ final
state in which all flavors are full, $\overrightarrow{\nu_{f}}=(+1,+1,+1,+1)$.
The interpolation scheme used to describe other valley-polarized states is
summarized in Table. 1.

The calculation method for the magnetization is given by the following steps.
First, one chooses a VP state (e.g. $\nu=+1$). Next, the initial and final
configuration points for the interpolarion are selected from Table. 1, and the
density matrix is calculated for every filling $\nu\in[\nu_{i},\nu_{f}]$ using
Eq. (25). Given the interpolated density matrix $\rho(\nu)$, we compute the
Hartree-Fock Hamiltonian in Eq. (23) without self-consistency. The
quasiparticle spectrum of this Hartree-Fock Hamiltonian provides our
approximate chemical potential $\mu(\nu)$, which is then integrated to find
the ground state energy. Repeating the procedure for different magnetic fields
allows one to compute the orbital magnetization using Eq. (24). Finally, the
integration constant is set to nullify the magnetization when all interpolated
flavors are completely empty or filled, which in Table 1 is always the case
for either $\nu=\nu_{i}$ or $\nu=\nu_{f}$. Numerical parameters and further
technical details are given in Appendix E.

In every interpolation configuration in Table 1 we set one pair of flavors
related by $\mathcal{T}$ to have equal filling at both initial and final
points, setting $\nu_{K\downarrow}(\nu)=\nu_{K^{\prime}\uparrow}(\nu)$ at
every intermediate filling $\nu$. By construction, these flavor’s
contributions to the magnetization cancel each other. A third flavor filling
$\nu_{K^{\prime}\downarrow}$ is kept fixed throughout the interpolation, and
therefore can only contribute to the magnetization’s integration constant.
Therefore, for our purpose of calculating the orbital magnetization, this
scheme allows us to focus on a single flavor block in the Hartree-Fock
Hamiltonian, with points of subtlety discussed in Appendix D.

####  IV.2.2 Magnetization results

Figure 4: Orbital in-plane magnetization vectors of individual Hartree-Fock
eigenstates at filling $\nu=+1$ in the valence band of valley $\tau=K$,
plotted across the moiré BZ for the VP $\nu=+1$ state of TBG with twist angle
of $\theta=1.1^{\circ}$. In this state, $\nu_{K^{\prime}\downarrow}=+1$ and
all other flavors are at zero filling. Background contours represent the
quasiparticle energy band. a) No strain, b) Heterostrain as specified in Fig.
1b. Individual magnetization vectors are similar in magnitude to those
obtained in the absence of interactions (Fig. 1), but their orientations are
somewhat less sensitive to heterostrain once interactions are introduced. The
hexagon outlines the moiré BZ.  Figure 5: The in-plane orbital magnetization
magnitude per moiré unit cell vs strain, for $\nu=+1$ VP states at differents
twist angle, compared with their non-interacting counterparts. The horizontal
axis $s_{0}$ denotes the principal strains as described in Fig. 3. Note that
the plots of the non-interacting model are valley-$K$ resolved, while the VP
states display a net magnetization. We find that in the vicinity of the magic
angle the interactions reduce sensitivity to strain, yet significant
magnetization is still obtained at small strains.

Fig. 4 compares between the magnetizations of individual Hartree-Fock
eigenstates of the unstrained and strained TBGs at filling $\nu=+1$, across
the moiré BZ. The calculation is done for the VP $\nu=+1$ configuration and
presented for the valence band of valley $K$. We find that, similarly to the
non-interacting case, individual states contribute $\sim 1\mu_{B}$ with and
without strain, but the 2D orientations seem more resilient to heterostrain in
the presence of interactions. In a sense, it seems that interactions stabilize
the $C_{3}$ symmetry.

In the presence of interactions, the total magnetization of the TBG is not the
sum of the magnetizations of the individual states plotted in Fig. 4. The
reason is that the individual Hartree-Fock eigenvalues depend on filling,
which is varied with each additional particle. Thus, the correct way to
calculate the magnetization is by means of Eq. (24).

In Fig. 5 we plot a comparison between the magnetizations of $\nu=+1$ VP
states at different twist angles and their non-interacting counterparts, as a
function of strain. We find that, in the vicinity of the magic angle,
interactions make the magnetization less sensitive to strain. Furthermore, the
dependencies of the magnetizations on strain for the three twist angles we
examined become closer to each other. However, note that even when
interactions are taken into account, we obtain a significant orbital
magnetization (larger than $1\mu_{B}$ per electron), even for strains as small
as $s_{0}\approx 4\cdot 10^{-4}$. We comment that the $\nu=+3$ VP states
behave similarly to the $\nu=+1$ VP states. Due to our approximate particle-
hole symmetry, the negative fillings states $\nu=-1$ and $\nu=-3$ show the
same magnetization magnitude as their symmetric counterparts.

We find an empirical correspondence between the bandwidth of the individual
(quasi-)particle spectrum and the strength of the magnetization response,
regardless of whether the spectrum originated from the interacting or non-
interacting model (Appendix C). This could be interpreted as a consequence of
the change in the Dirac velocity caused by interactions, which mimics a change
in the twist angle.

##  V Discussion and conclusions

In this paper we have shown that valley-polarized states in heterostrained
MATBG have a significant in-plane orbital magnetization, which may approach
the order of $\sim 1\mu_{B}$ per moiré unit cell. This orbital magnetization
can be thought of as arising from individual magnetic moments of Bloch states
partially aligning due to $C_{3}$ symmetry breaking, with their time-reversed
partners unoccupied due to the state being valley-polarized. The flatness of
the bands in the vicinity of the magic angle plays a dual role in this
picture: it increases the relative strength of electron-electron interaction,
which allows for the $C_{2}$ and $\mathcal{T}$-broken valley polarized states,
and it increases the sensitivity of the band structure to $C_{3}$ symmetry
breaking due to heterostrain. Band renormalization due to interactions is
accounted for using a simplified model, and near the magic angle its main
effect is an increase in bandwidth and consequently a slight decrease in
sensitivity to strain. Even in the presence of this band renormalization,
typical strain values of $s_{0}\sim 10^{-4}-10^{-3}$ Kerelsky  _et al._
(2019); Choi  _et al._ (2019); Xie  _et al._ (2019); Kazmierczak and Van
Winkle (2021) are enough to generate an orbital magnetization of $\sim
1\mu_{B}$ per moiré unit cell for twist angles near the magic angle, a
surprising result considering that $s_{0}/\theta\sim 0.01-0.1$ in this
parameter regime.

Our conclusions are of direct relevance to experiments which exhibit evidence
of flavor symmetry breaking in the form of spin and valley polarized states at
$\nu=+1,+3$Zondiner  _et al._ (2020); Park  _et al._ (2021); Sharpe  _et al._
(2019); Serlin  _et al._ (2020). In these states, the total magnetic moment
comprises of orbital and spin components, with the latter being at most
$1\mu_{B}$ per moiré unit cell. Any magnetization in excess of this value is a
strong indication of an orbital component. There is a further signature that
our model predicts - the orbital component is determined by the direction of
the strain, as opposed to the isotropic spin magnetization.

In fact, it is possible that this orbital component has already been observed
in Ref. Zondiner  _et al._ , 2020, where the in-plane magnetization was
extracted from electronic compressibility measurements in the presence of a
field. Close to the magic angle and near $\nu=1$, the in-plane magnetization
exceeds $1\mu_{B}$ per unit cell. As expected, this effect disappears away
from the magic angle. However, the direction of the magnetization in the plane
was not studied in the experiment. Our analysis predicts the relation of this
direction to the strain.

Our result suggests that in the presence of strain, an in-plane magnetic field
induces a valley-Zeeman effect, which favours valley-polarized states. In a
sample that contains domains with different strains, we expect the in-plane
magnetization to be different in each domain, resulting in an out-of-plane
field at the domain boundaries. We emphasize that in our model we assumed that
our system has $C_{2}\mathcal{T}$ symmetry (there is no sublattice potential),
and hence there is no out-of-plane magnetization in the bulk, in contrast with
the system described in Refs. Zhang  _et al._ , 2020; He  _et al._ , 2020;
Bultinck  _et al._ , 2020.

For future works, we suggest exploring the in-plane magnetization of nematic
states where the $C_{3}$ symmetry is broken spontaneously Liu  _et al._
(2021); Cao  _et al._ (2021); Choi  _et al._ (2019); Jiang  _et al._ (2019).
In such a phase, the direction of the magnetization is presumably related to
the orientation of the nematic order parameter. Therefore, one could imagine
training the nematic order using an in-plane magnetic field.

###### Acknowledgements.

We thank Keshav Pareek and Achim Rosch for useful discussions. We acknowledge
support from the Israeli Science Foundation Quantum Science and Technology
grant no. 2074/19, and from the CRC 183 of the Deutsche Forschungsgemeinschaft
(Project C02). This work has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme [grant agreements No. 788715 and 817799, Projects LEGOTOP (AS) and
HQMAT (EB)].

## References

  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” Proceedings of the National Academy of Sciences of the United States of America 108, 12233–12237 (2011), arXiv:1009.4203 .
  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80–84 (2018a).
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018b).
  * Yankowitz  _et al._ (2019) Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn, Yuxuan Zhang, K Watanabe, T Taniguchi, David Graf, Andrea F Young, and Cory R Dean, “Tuning superconductivity in twisted bilayer graphene,” Science (2019).
  * Arora  _et al._ (2020) Harpreet Singh Arora, Robert Polski, Yiran Zhang, Alex Thomson, Youngjoon Choi, Hyunjin Kim, Zhong Lin, Ilham Zaky Wilson, Xiaodong Xu, Jiun-Haw Chu, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, Stevan Nadj-Perge, and Yiran Zhang, “Superconductivity in metallic twisted bilayer graphene stabilized by WSe 2,” Nature 583, 379 (2020).
  * Choi  _et al._ (2019) Youngjoon Choi, Jeannette Kemmer, Yang Peng, Alex Thomson, Harpreet Arora, Robert Polski, Yiran Zhang, Hechen Ren, Jason Alicea, Gil Refael, Felix von Oppen, Kenji Watanabe, Takashi Taniguchi, and Stevan Nadj-Perge, “Electronic correlations in twisted bilayer graphene near the magic angle,” Nature Physics 15, 1174–1180 (2019), arXiv:1901.02997 .
  * Wong  _et al._ (2020) Dillon Wong, Kevin P Nuckolls, Myungchul Oh, Biao Lian, Yonglong Xie, Sangjun Jeon, Kenji Watanabe, Takashi Taniguchi, B Andrei Bernevig, and Ali Yazdani, “Cascade of electronic transitions in magic-angle twisted bilayer graphene Sequence of spectroscopic transitions,” Nature 582 (2020), 10.1038/s41586-020-2339-0.
  * Saito  _et al._ (2020) Yu Saito, Jingyuan Ge, Kenji Watanabe, Takashi Taniguchi, and Andrea F Young, “Independent superconductors and correlated insulators in twisted bilayer graphene,” Nature Physics (2020), 10.1038/s41567-020-0928-3.
  * Stepanov  _et al._ (2020) Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe, Takashi Taniguchi, Frank H L Koppens, Johannes Lischner, Leonid Levitov, and Dmitri K Efetov, “Untying the insulating and superconducting orders in magic-angle graphene,” Nature 583, 375 (2020).
  * Uri  _et al._ (2020) A Uri, S Grover, Y Cao, J A Crosse, K Bagani, D Rodan-Legrain, Y Myasoedov, K Watanabe, T Taniguchi, P Moon, M Koshino, P Jarillo-Herrero, and E Zeldov, “Mapping the twist-angle disorder and Landau levels in magic-angle graphene,” Nature 581 (2020), 10.1038/s41586-020-2255-3.
  * Xie  _et al._ (2019) Yonglong Xie, Biao Lian, Berthold Jäck, Xiaomeng Liu, Cheng Li Chiu, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Ali Yazdani, “Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene,” Nature 572, 101–105 (2019), arXiv:1906.09274 .
  * Rozen  _et al._ (2021) Asaf Rozen, Jeong Min Park, Uri Zondiner, Yuan Cao, Daniel Rodan-Legrain, Takashi Taniguchi, Kenji Watanabe, Yuval Oreg, Ady Stern, Erez Berg, Pablo Jarillo-Herrero, and Shahal Ilani, “Entropic evidence for a Pomeranchuk effect in magic-angle graphene,” Nature 592 (2021), 10.1038/s41586-021-03319-3.
  * Sharpe  _et al._ (2019) Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, Joe Finney, Kenji Watanabe, Takashi Taniguchi, M. A. Kastner, and David Goldhaber-Gordon, “Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene,” Science 608, 605–608 (2019).
  * Park  _et al._ (2021) Jeong Min Park, Yuan Cao, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, “Flavour Hund’s coupling, Chern gaps and charge diffusivity in moiré graphene,” Nature 592, 43 (2021).
  * Jiang  _et al._ (2019) Yuhang Jiang, Xinyuan Lai, Kenji Watanabe, Takashi Taniguchi, Kristjan Haule, Jinhai Mao, and Eva Y. Andrei, “Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene,” Nature 573, 91–95 (2019).
  * Kerelsky  _et al._ (2019) Alexander Kerelsky, Leo J. McGilly, Dante M. Kennes, Lede Xian, Matthew Yankowitz, Shaowen Chen, K. Watanabe, T. Taniguchi, James Hone, Cory Dean, Angel Rubio, and Abhay N. Pasupathy, “Maximized electron interactions at the magic angle in twisted bilayer graphene,” Nature 572, 95–100 (2019).
  * Lu  _et al._ (2019) Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, Adrian Bachtold, Allan H. MacDonald, and Dmitri K. Efetov, “Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene,” Nature 574, 653–657 (2019).
  * Zondiner  _et al._ (2020) U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, “Cascade of phase transitions and Dirac revivals in magic-angle graphene,” Nature 582, 203–208 (2020), arXiv:1912.06150 .
  * Serlin  _et al._ (2020) M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, “Intrinsic quantized anomalous Hall effect in a moiré heterostructure,” Science 367, 900–903 (2020), arXiv:1907.00261 .
  * Cao  _et al._ (2021) Yuan Cao, Daniel Rodan-Legrain, Jeong Min Park, Noah F.Q. Yuan, Kenji Watanabe, Takashi Taniguchi, Rafael M. Fernandes, Liang Fu, and Pablo Jarillo-Herrero, “Nematicity and competing orders in superconducting magic-angle graphene,” Science 372, 264–271 (2021).
  * Tschirhart  _et al._ (2021) C. L. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. E. Huber, and A. F. Young, “Imaging orbital ferromagnetism in a moiré Chern insulator,” Science 372, 1323–1327 (2021), arXiv:2006.08053 .
  * Saito  _et al._ (2021) Yu Saito, Fangyuan Yang, Jingyuan Ge, Xiaoxue Liu, Takashi Taniguchi, Kenji Watanabe, J I A Li, Erez Berg, and Andrea F Young, “Isospin Pomeranchuk effect in twisted bilayer graphene,” Nature 592 (2021), 10.1038/s41586-021-03409-2.
  * Breit (1928) Gregory Breit, “An Interpretation of Dirac’s Theory of the Electron,” Proceedings of the National Academy of Sciences 14, 553–559 (1928).
  * Kazmierczak and Van Winkle (2021) Nathanael P Kazmierczak and Madeline Van Winkle, “Strain fields in twisted bilayer graphene,” Nature Materials (2021), 10.1038/s41563-021-00973-w.
  * Bi  _et al._ (2019) Zhen Bi, Noah F Q Yuan, and Liang Fu, “Designing flat bands by strain,” Physical Review B 100, 35448 (2019).
  * Balents (2019) Leon Balents, “General continuum model for twisted bilayer graphene and arbitrary smooth deformations,” SciPost Phys 7, 48 (2019).
  * Suzuura and Ando (2002) Hidekatsu Suzuura and Tsuneya Ando, “Phonons and electron-phonon scattering in carbon nanotubes,” Physical Review B - Condensed Matter and Materials Physics 65, 1–15 (2002).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M.R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Reviews of Modern Physics 81, 109–162 (2009), arXiv:0709.1163 .
  * Nam and Koshino (2017) Nguyen N T Nam and Mikito Koshino, “Lattice relaxation and energy band modulation in twisted bilayer graphene,” Physical Review B 96, 75311 (2017).
  * Koshino  _et al._ (2018) Mikito Koshino, Noah F Q Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally Localized Wannier Orbitals and the Extended Hubbard Model for Twisted Bilayer Graphene,” Physical Review X (2018), 10.1103/PhysRevX.8.031087.
  * Carr  _et al._ (2019) Stephen Carr, Shiang Fang, Ziyan Zhu, and Efthimios Kaxiras, “Exact continuum model for low-energy electronic states of twisted bilayer graphene,” Phys. Rev. Research 1, 013001 (2019).
  * Huder  _et al._ (2018) Loïc Huder, Alexandre Artaud, Toai Le Quang, Guy Trambly De Laissardière, Aloysius G.M. Jansen, Gérard Lapertot, Claude Chapelier, and Vincent T. Renard, “Electronic Spectrum of Twisted Graphene Layers under Heterostrain,” Physical Review Letters 120, 1–5 (2018), arXiv:1803.03505 .
  * Qiao  _et al._ (2018) Jia Bin Qiao, Long Jing Yin, and Lin He, “Twisted graphene bilayer around the first magic angle engineered by heterostrain,” Physical Review B 98, 1–8 (2018).
  * Liu  _et al._ (2021) Shang Liu, Eslam Khalaf, Jong Yeon Lee, and Ashvin Vishwanath, “Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality,” Physical Review Research 3, 13033 (2021).
  * Parker  _et al._ (2021) Daniel E Parker, Tomohiro Soejima, Johannes Hauschild, Michael P Zaletel, and Nick Bultinck, “Strain-Induced Quantum Phase Transitions in Magic-Angle Graphene,” Physical Review Letters 127 (2021), 10.1103/PhysRevLett.127.027601.
  * Bernevig  _et al._ (2021) B Andrei Bernevig, Zhi-Da Song, Nicolas Regnault, and Biao Lian, “Twisted bilayer graphene. I. Matrix elements, approximations, perturbation theory, and a $k\cdot p$ two-band model,” Physical Review B 103, 205411 (2021).
  * Kwan  _et al._ (2020) Yves H. Kwan, S. A. Parameswaran, and S. L. Sondhi, “Twisted bilayer graphene in a parallel magnetic field,” Phys. Rev. B 101, 205116 (2020).
  * Xie and Macdonald (2020) Ming Xie and A. H. Macdonald, “Nature of the Correlated Insulator States in Twisted Bilayer Graphene,” Physical Review Letters 124, 97601 (2020).
  * Xie and Macdonald (2021) Ming Xie and A H Macdonald, “Weak-Field Hall Resistivity and Spin-Valley Flavor Symmetry Breaking in Magic-Angle Twisted Bilayer Graphene,” Physical Review Letters 127 (2021), 10.1103/PhysRevLett.127.196401.
  * Zhang  _et al._ (2020) Yu Zhang, Zhe Hou, Ya-Xin Xin Zhao, Zi-Han Han Guo, Yi-Wen Wen Liu, Si-Yu Yu Li, Ya-Ning Ning Ren, Qing-Feng Feng Sun, and Lin He, “Correlation-induced valley splitting and orbital magnetism in a strain-induced zero-energy flatband in twisted bilayer graphene near the magic angle,” Physical Review B 102 (2020), 10.1103/PhysRevB.102.081403.
  * He  _et al._ (2020) Wen Yu He, David Goldhaber-Gordon, and K. T. Law, “Giant orbital magnetoelectric effect and current-induced magnetization switching in twisted bilayer graphene,” Nature Communications 11, 1–8 (2020), arXiv:1908.11718 .
  * Bultinck  _et al._ (2020) Nick Bultinck, Shubhayu Chatterjee, and Michael P Zaletel, “Mechanism for Anomalous Hall Ferromagnetism in Twisted Bilayer Graphene,” Physical Review Letters 124 (2020), 10.1103/PhysRevLett.124.166601.

##  Appendix A Tripod model with strain

In this section we derive Eqs. (18-20) within the “tripod model”. Our starting
point is the following Hamiltonian:

$\displaystyle\begin{split}H(\mathbf{k})&=h^{(1)}\left(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)\\\
&-\sum_{j=0,1,2}T_{j}\left[h^{(2)}\left(\mathbf{k}+\frac{e}{\hbar}\mathbf{A}_{\text{em}}+\mathbf{g}_{j}\right)-E(\mathbf{k})\right]^{-1}T_{j}^{\dagger},\end{split}$
(26)

which is valid at small $\mathbf{p}=\mathbf{k}-\delta\mathbf{K}^{(1)}$ and
small parallel magnetic field. To simplify our derivation, we consider a
shifted momentum
$\mathbf{p}^{\prime}=\mathbf{p}+\frac{e}{\hbar}\mathbf{A}_{\text{em}}$ such
that all the magnetic field terms appear on the first layer:

$\displaystyle
h^{(1)}\left(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)=\hbar
v_{F}\left(\mathbf{p}^{\prime}-2\frac{e}{\hbar}\mathbf{A}_{\text{em}}\right)\cdot\bm{\sigma},$
(27a) $\displaystyle
h^{(2)}\left(\mathbf{k}+\frac{e}{\hbar}\mathbf{A}_{\text{em}}+\mathbf{g}_{j}\right)=\hbar
v_{F}(\mathbf{p}^{\prime}+\mathbf{g}_{j}+\delta\mathbf{K})\cdot\bm{\sigma},$
(27b)

where $\delta\mathbf{K}=\delta\mathbf{K}^{(1)}-\delta\mathbf{K}^{(2)}$, and we
neglected the rotation of the Pauli matrices as we focus on leading order
effects. The reference scale for momentum and energy is set by introducing the
following dimensionless parameters:

$\tilde{\mathbf{p}}=\frac{\mathbf{p}^{\prime}}{K\theta},\quad\tilde{E}=\frac{E}{\hbar
v_{F}K\theta},\quad\mathbf{d}_{j}=\frac{1}{K\theta}\left(\mathbf{g}_{j}+\delta\mathbf{K}\right).$
(28)

The corresponding dimensionless tunneling is given:

$\tilde{T}_{j}=u+w\bm{\zeta}_{j}\cdot\bm{\sigma},\quad\bm{\zeta}_{j}=\left(\text{Re}\left[\zeta^{j}\right],-\text{Im}\left[\zeta^{j}\right]\right),\quad\zeta=e^{i\frac{2\pi}{3}}.$
(29)

In terms of these definitions, our tripod model can be rewritten as:

$\displaystyle\begin{split}\tilde{H}(\tilde{\mathbf{p}})&=\left(\tilde{\mathbf{p}}-2\frac{e}{\hbar
K\theta}\mathbf{A}_{\text{em}}\right)\cdot\bm{\sigma}\\\
&-\sum_{j=0,1,2}\tilde{T}_{j}\left[\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}\right)\cdot\bm{\sigma}+\tilde{E}(\tilde{\mathbf{p}})\right]^{-1}\tilde{T}_{j}.\end{split}$
(30)

We are interested in expanding Eq. (30) to first order in
$\tilde{\mathbf{p}}$, $\tilde{E}$ and in components of $\mathcal{S}/\theta$.
Specifically, the task at hand is expanding the inverse of the second layer’s
Hamiltonian. We start by explicitly writing the form of $\mathbf{d}_{j}$:

$\displaystyle\begin{split}\mathbf{d}_{j}&=\frac{1}{K\theta}\left[\mathcal{U}^{T}\left(\mathbf{G}_{j}-\mathbf{K}\right)-\left(\mathbf{A}_{s}^{(1)}-\mathbf{A}_{s}^{(2)}\right)\right]\\\
&=\frac{1}{\theta}\mathcal{R}(\theta)\bm{\zeta}_{j}-\frac{1}{\theta}\left[\mathcal{S}\bm{\zeta}_{j}+\frac{1}{K}\left(\mathbf{A}_{s}^{(1)}-\mathbf{A}_{s}^{(2)}\right)\right]\\\
&=\mathbf{d}_{j}^{(0)}+\mathbf{d}_{j}^{(1)},\end{split}$ (31)

where we used $\mathbf{K}-\mathbf{G}_{j}=K\bm{\zeta}_{j}$, and separated the
$\mathcal{O}(1)$ contribution,
$\mathbf{d}_{j}^{(0)}=\theta^{-1}\mathcal{R}(\theta)\bm{\zeta}_{j}$, from the
$\mathcal{O}(\mathcal{S}/\theta)$ corrections,
$\mathbf{d}_{j}^{(1)}=\mathbf{d}_{j}-\mathbf{d}_{j}^{(0)}$. It will be useful
later to introduce a dimensionless artificial gauge field
$\tilde{\mathbf{A}}_{s}$ such that
$\mathbf{d}_{j}^{(1)}=-\theta^{-1}\mathcal{S}\zeta_{j}-\tilde{\mathbf{A}}_{s}$.
Next, we invert the operator between the tunneling matrices, and keep only
first order correction in strain and momentum. Using
$\left|\mathbf{d}_{j}^{(0)}\right|^{2}=1$, we find:

$\displaystyle\begin{split}&\left[\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}\right)\cdot\bm{\sigma}+\tilde{E}\right]^{-1}\\\
&\approx\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}^{(1)}+\left[1-2\mathbf{d}_{j}^{(0)}\cdot\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}^{(1)}\right)\right]\mathbf{d}_{j}^{(0)}\right)\cdot\bm{\sigma}-\tilde{E}\\\
&=\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E},\end{split}$ (32)

where the last line defines $\mathbf{f}_{j}$.

Using the form of the dimensionless tunneling in Eq. (29), we have the
following identity:

$\displaystyle\begin{split}\tilde{T}_{j}(\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E})\tilde{T}_{j}&=\left[(u^{2}-w^{2})\mathbf{f}_{j}+2w^{2}(\mathbf{f}_{j}\cdot\bm{\zeta}_{j})\bm{\zeta}_{j}\right]\cdot\bm{\sigma}\\\
&-2uw\tilde{E}\bm{\zeta}_{j}\cdot\bm{\sigma}+2uw\mathbf{f}_{j}\cdot\bm{\zeta}_{j}-(u^{2}+w^{2})\tilde{E}.\end{split}$
(33)

Noting that the following cancel - $\sum_{j}\bm{\zeta}_{j}=0$,
$\sum_{j}\mathbf{d}_{j}^{(0)}=0$ and
$\mathbf{d}_{j}^{(0)}\cdot\bm{\zeta}_{j}=0$, we insert Eq. (32) into Eq. (33)
and obtain:

$\displaystyle\begin{split}&\sum_{j}\tilde{T}_{j}(\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E})\tilde{T}_{j}\\\
&=(u^{2}-w^{2})\sum_{j}\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}^{(1)}-2\left[\mathbf{d}_{j}^{(0)}\cdot\left(\tilde{\mathbf{p}}+\mathbf{d}_{j}^{(1)}\right)\right]\mathbf{d}_{j}^{(0)}\right)\cdot\bm{\sigma}\\\
&+2w^{2}\sum_{j}((\tilde{\mathbf{p}}+\mathbf{d}_{j}^{(1)})\cdot\bm{\zeta}_{j})\bm{\zeta}_{j}\cdot\bm{\sigma}\\\
&-3(u^{2}+w^{2})\tilde{E}+2uw\sum_{j}\mathbf{d}_{j}^{(1)}\cdot\bm{\zeta}_{j}.\end{split}$
(34)

Next, we perform the summation over $j$. For brevity, we divide it into two
steps. First, we reinstate the explicit form of $\mathbf{d}_{j}^{(1)}$ given
by Eq. (31). Making use of
$\sum_{j}(\tilde{\mathbf{p}}\cdot\bm{\zeta}_{j})\cdot\bm{\zeta}_{j}=\sum_{j}(\tilde{\mathbf{p}}\cdot\mathbf{d}_{j}^{(0)})\cdot\mathbf{d}_{j}^{(0)}=\frac{3}{2}\tilde{\mathbf{p}}$,
we rewrite Eq. (34) as:

$\displaystyle\begin{split}&\sum_{j}\tilde{T}_{j}(\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E})\tilde{T}_{j}\\\
&=2(u^{2}-w^{2})\theta^{-1}\sum_{j}\left[\mathbf{d}_{j}^{(0)}\cdot\left(\mathcal{S}\bm{\zeta}_{j}\right)\right]\mathbf{d}_{j}^{(0)}\cdot\bm{\sigma}\\\
&+3w^{2}(\tilde{\mathbf{p}}-\tilde{\mathbf{A}}_{s})\cdot\bm{\sigma}-2w^{2}\theta^{-1}\sum_{j}\left[\left(\mathcal{S}\bm{\zeta}_{j}\right)\cdot\bm{\zeta}_{j}\right]\bm{\zeta}_{j}\cdot\bm{\sigma}\\\
&-3(u^{2}+w^{2})\tilde{E}-2uw\theta^{-1}\sum_{j}\left(\mathcal{S}\bm{\zeta}_{j}\right)\cdot\bm{\zeta}_{j}.\end{split}$
(35)

Second, we perform the summations involving the vectors
$\mathbf{d}_{j}^{(0)}$, $\bm{\zeta}_{j}$ and the tensor $\mathcal{S}$. These
require a bit more algebra, and result in the following expression:

$\displaystyle\begin{split}&\sum_{j}\tilde{T}_{j}(\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E})\tilde{T}_{j}\\\
&=3w^{2}(\tilde{\mathbf{p}}-\tilde{\mathbf{A}}_{s})\cdot\bm{\sigma}-3(u^{2}+w^{2})\tilde{E}-3uw(s^{\prime}+s^{\prime\prime})\theta^{-1}\\\
&-\frac{3}{2}u^{2}(s^{\prime}-s^{\prime\prime})\theta^{-1}\hat{\mathbf{A}}_{s}\cdot\bm{\sigma},\end{split}$
(36)

where the last term contains a unit vector $\hat{\mathbf{A}}_{s}$ (with the
direction of $\tilde{\mathbf{A}}_{s}$) that originated from manipulations over
$\mathcal{S}$. Using $\tilde{\mathbf{A}}_{s}=\frac{\beta}{2\theta
Ka}(s^{\prime}-s^{\prime\prime})\hat{\mathbf{A}}_{s}$, we can rewrite the
above as:

$\displaystyle\begin{split}&\sum_{j}\tilde{T}_{j}(\mathbf{f}_{j}\cdot\bm{\sigma}-\tilde{E})\tilde{T}_{j}\\\
&=3w^{2}\left[\tilde{\mathbf{p}}-\left(1+\frac{u^{2}}{w^{2}}\frac{Ka}{\beta}\right)\tilde{\mathbf{A}}_{s}\right]\cdot\bm{\sigma}\\\
&-3(u^{2}+w^{2})\tilde{E}-3uw(s^{\prime}+s^{\prime\prime})\theta^{-1}.\end{split}$
(37)

The outcome of this calculation is given by inserting Eq. (37) back into the
Eq. (30) for the effective Hamiltonian:

$\displaystyle\begin{split}&H_{\text{eff}}(\tilde{\mathbf{p}})=\\\
&\frac{\left((1-3w^{2})\tilde{\mathbf{p}}-2\frac{e}{\hbar
K\theta}\mathbf{A}_{\text{em}}+3w^{2}\left(1+\frac{u^{2}}{w^{2}}\frac{Ka}{\beta}\right)\tilde{\mathbf{A}}_{s}\right)\cdot\bm{\sigma}}{1+3(u^{2}+w^{2})}\\\
&+\frac{3uw}{1+3(u^{2}+w^{2})}\frac{s^{\prime}+s^{\prime\prime}}{\theta}.\end{split}$
(38)

Restoring the dimensions of the parameters and rearranging to emphasize the
shift of the Dirac point, Eq. (38) is rewritten as:

$\displaystyle\begin{split}&H_{\text{eff}}(\mathbf{p})=\\\ &\frac{\hbar
v_{F}(1-3w^{2})\left(\mathbf{p}-\frac{3w^{2}}{1-3w^{2}}\left(1+\frac{u^{2}}{w^{2}}\frac{Ka}{\beta}\right)\mathbf{A}_{s}\right)\cdot\bm{\sigma}}{1+3(u^{2}+w^{2})}\\\
&+\frac{2}{1+3(u^{2}+w^{2})}ev_{F}\mathbf{A}_{\text{em}}\cdot\bm{\sigma}\\\
&+\hbar
v_{F}K\frac{3uw}{1+3(u^{2}+w^{2})}(s^{\prime}+s^{\prime\prime}).\end{split}$
(39)

One can clearly see from the result in Eq. (39) the sensitivity of the
crossing point on the magic angle condition $w^{2}=1/3$ in the denominator,
and that the response to the magnetic field is modified due to tunneling by a
prefactor that does not vanish at the magic angle.

##  Appendix B Symmetry constraints on magnetization for higher orders in
strain

###  B.1 Effect of $C_{3}$ symmetry on the energy

For the purpose of deriving the constraints imposed on the energy by the
$C_{3}$ symmetry we can omit the valley fillings as they are conserved.
Focusing on the strain and magnetic field, we rewrite (11) here as:

$E\left(\mathcal{S},\mathbf{B}_{\parallel}\right)=E\left(C_{3}\mathcal{S}C^{T}_{3},C_{3}\mathbf{B}_{\parallel}\right).$
(40)

For the zero-field magnetization, we are interested only in terms of Eq. (40)
that are linear in $\mathbf{B}_{\parallel}$.

The constraint in Eq. (40) allows only for combinations of strain tensor
components $s_{\alpha\beta}$ and magnetic field vector components $B_{\gamma}$
that remain invariant under $C_{3}$. The transformation laws for these
components are given by:

$\displaystyle\begin{split}s_{\alpha\beta}&\rightarrow\left[C_{3}\right]_{\alpha\alpha^{\prime}}\left[C_{3}\right]_{\beta\beta^{\prime}}s_{\alpha^{\prime}\beta^{\prime}}\\\
B_{\gamma}&\rightarrow\left[C_{3}\right]_{\gamma\gamma^{\prime}}B_{\gamma^{\prime}}\end{split}$
(41)

where the $C_{3}$ rotation tensor components are
$\left[C_{3}\right]_{\alpha\alpha^{\prime}}=\left[e^{i\frac{2\pi}{3}\eta_{y}}\right]_{\alpha\alpha^{\prime}}$.

The $C_{3}$ rotation tensor has the following eigenvectors and corresponding
eigenvalues:

$\displaystyle\begin{split}\left[C_{3}\right]_{\alpha\alpha^{\prime}}(\delta_{\alpha^{\prime},x}+i\delta_{\alpha^{\prime},y})&=e^{i2\pi/3}(\delta_{\alpha,x}+i\delta_{\alpha,y}),\\\
\left[C_{3}\right]_{\alpha\alpha^{\prime}}(\delta_{\alpha^{\prime},x}-i\delta_{\alpha^{\prime},y})&=e^{-i2\pi/3}(\delta_{\alpha,x}-i\delta_{\alpha,y}).\end{split}$
(42)

Eq. (42) forces the magnetic field to enter our invariant combination through
either $B_{x}+iB_{y}$ or $B_{x}-iB_{y}$. Focusing on the former, we have an
eigenvalue of $e^{i2\pi/3}$ corresponding to the magnetic field. For the
energy to be invariant under $C_{3}$, the strain field must present an
eigenvalue of $e^{-i2\pi/3}$. This implies that a zeroth order term in strain
is prohibited from coupling to a linear term in magnetic field, and thus there
is no in-plane magnetization without strain. We will limit the discussion from
here on to first and second order in strain, but the generalization is
trivial.

To linear order in strain, the required eigenvalue of $e^{-i2\pi/3}$ can be
obtained only by the following combination:

$\displaystyle\begin{split}(\delta_{\alpha,x}+i\delta_{\alpha,y})(\delta_{\beta,x}+i\delta_{\beta,y})s_{\alpha\beta}&=s_{xx}-s_{yy}+i2s_{xy}\\\
&=(s^{\prime}-s^{\prime\prime})e^{i2\phi_{S}},\end{split}$ (43)

where the second line is written in terms of principal strains
$s^{\prime},s^{\prime\prime}$ and angle $\phi_{S}$ as defined in Eq. (12), and
elucidates the rotational nature of this term. A similar argument couples
$B_{x}-iB_{y}$ to $(s^{\prime}-s^{\prime\prime})e^{-i2\phi_{S}}$. Rewriting
the magnetic field as
$\mathbf{B}=|B|\left(\cos{\phi_{B}},\sin{\phi_{B}}\right)$, and demanding the
energy to be real, we find that to linear order in strain and magnetic field:

$E(\mathcal{S},\mathbf{B}_{\parallel})\approx\lambda_{1}(s^{\prime}-s^{\prime\prime})|B|\cos{\left(2\phi_{S}+\phi_{B}+\phi_{0}\right)},$
(44)

for some $\phi_{0}\in\left[0,\pi\right)$ which will be determined by other
symmetries [see below Eq. (16)]. The prefactor $\lambda_{1}$ is the coupling
coefficient at this order.

The second order in strain must also present the same eigenvalue of
$e^{-i2\pi/3}$. In contrast with the linear order, we have two ways of
obtaining this eigenvalue:

$\displaystyle\begin{split}&(\delta_{\alpha,x}+i\delta_{\alpha,y})(\delta_{\beta,x}+i\delta_{\beta,y})(\delta_{\mu,x}+i\delta_{\mu,y})(\delta_{\nu,x}-i\delta_{\nu,y})s_{\alpha\beta}s_{\mu\nu}\\\
&=\left((s^{\prime})^{2}-(s^{\prime\prime})^{2}\right)e^{i2\phi_{S}},\end{split}$
(45)

and

$\displaystyle\begin{split}&(\delta_{\alpha,x}-i\delta_{\alpha,y})(\delta_{\beta,x}-i\delta_{\beta,y})(\delta_{\mu,x}-i\delta_{\mu,y})(\delta_{\nu,x}-i\delta_{\nu,y})s_{\alpha\beta}s_{\mu\nu}\\\
&=\left(s^{\prime}-s^{\prime\prime}\right)^{2}e^{-i4\phi_{S}}.\end{split}$
(46)

Following similar arguments as with the linear order in strain, we find the
energy to have the form:

$\displaystyle\begin{split}E(\mathcal{S},\mathbf{B}_{\parallel})&\approx\lambda_{1}(s^{\prime}-s^{\prime\prime})|B|\cos{\left(2\phi_{S}+\phi_{B}+\phi_{0}\right)}\\\
&+\lambda_{2}\left((s^{\prime})^{2}-(s^{\prime\prime})^{2}\right)|B|\cos{\left(2\phi_{S}+\phi_{B}+\phi_{1}\right)}\\\
&+\lambda_{3}\left(s^{\prime}-s^{\prime\prime}\right)^{2}|B|\cos{\left(-4\phi_{S}+\phi_{B}+\phi_{2}\right)},\end{split}$
(47)

for three general angles $\phi_{0},\phi_{1},\phi_{2}\in\left[0,\pi\right)$,
and we find two new coupling coefficients $\lambda_{2}$ and $\lambda_{3}$. The
angles $\phi_{0,1,2}$ are fixed by $R_{\pi}^{x,y}$ symmetries, as explained in
the next subsection.

###  B.2 Effects of $R_{\pi}^{x}$ and $R_{\pi}^{y}$ symmetries on the energy

We now turn to explore further constraints imposed by the $R_{\pi}^{x}$ and
$R_{\pi}^{y}$ symmetries. To do so, we must specify how the strain tensor
changes due to layer flip, and reintroduce the valley filling dependence in
the coupling coefficients $\lambda_{i}(\nu_{K},\nu_{K^{\prime}})$. In order to
account for both heterostrain and homostrain, we take $\chi=\pm 1$ to be the
eigenvalue of the layer flip applied to the strain tensor, i.e. $\chi=-1$
($\chi=+1$) for the heterostrain (homostrain). The $R_{\pi}^{x}$ symmetry
preserves the valleys and takes $\mathcal{S}\rightarrow\chi\mathcal{S}$,
$\phi_{S}\rightarrow-\phi_{S}$ and $\phi_{B}\rightarrow-\phi_{B}$. For an
energy in the form of Eq. (47), we find that $\phi_{1}=\phi_{2}=0$, and that
for heterostrain (homostrain) we have $\phi_{0}=\frac{\pi}{2}$ ($\phi_{0}=0$).
We can combine both the case heterostrain and of homostrain into a single
equation by writing:

$\displaystyle\begin{split}E(\mathcal{S},\mathbf{B}_{\parallel})&\approx\\\
&\left[\lambda_{1}(\nu_{K},\nu_{K^{\prime}})(s^{\prime}-s^{\prime\prime})\cos{\left(2\phi_{S}+\phi_{B}+\frac{(1-\chi)\pi}{4}\right)}\right.\\\
&+\lambda_{2}(\nu_{K},\nu_{K^{\prime}})\left((s^{\prime})^{2}-(s^{\prime\prime})^{2}\right)\cos{\left(2\phi_{S}+\phi_{B}\right)}\\\
&+\left.\lambda_{3}(\nu_{K},\nu_{K^{\prime}})\left(s^{\prime}-s^{\prime\prime}\right)^{2}\cos{\left(-4\phi_{S}+\phi_{B}\right)}\right]|B|,\end{split}$
(48)

and selecting the desired value of $\chi$.

Next, we consider the $R_{\pi}^{y}$ symmetry. This symmetry flips the valleys
and takes $\mathcal{S}\rightarrow\chi\mathcal{S}$,
$\phi_{S}\rightarrow\pi-\phi_{S}$ and $\phi_{B}\rightarrow\pi-\phi_{B}$.
Changing the valleys is associated with a sign flip of the magnetization as
seen by Eq. (10). This valley flip simply adds the rule
$\lambda_{i}\rightarrow-\lambda_{i}$ for $i=1,2,3$. Using these transformation
rules, one finds that the $R_{\pi}^{y}$ symmetry does not impose any
additional constraints, for both the heterostrained and homostrained TBG.
Therefore, our final result is given by Eq. (48). The effects of $C_{2}$ and
$\mathcal{T}$ are described in the main text. The magnetization immediately
follows from this energy, as it is given to linear order in magnetic field:

$\displaystyle\begin{split}M_{\parallel,x}+iM_{\parallel,y}&\approx-\lambda_{1}(\nu_{K},\nu_{K^{\prime}})(s^{\prime}-s^{\prime\prime})e^{i(\chi-1)\pi/4}e^{-i2\phi_{S}}\\\
&-\lambda_{2}(\nu_{K},\nu_{K^{\prime}})\left((s^{\prime})^{2}-(s^{\prime\prime})^{2}\right)e^{-i2\phi_{S}}\\\
&-\lambda_{3}(\nu_{K},\nu_{K^{\prime}})\left(s^{\prime}-s^{\prime\prime}\right)^{2}e^{i4\phi_{S}}.\end{split}$
(49)

The result in Eq. (49) is important, as it tells us that the magnetization
magnitude $\left|\mathbf{M}_{\parallel}\right|$ is no longer $\phi_{S}$
independent at second order in strain. There is no way for symmetry arguments
to predict the amount of variance expected, nor when higher orders become
significant. We can, however, plot the magnetization magnitude and orientation
as a function of angle $\phi_{S}$ for a strain tensor with
$s^{\prime}=-s^{\prime\prime}=s_{0}$, and compare different values of $s_{0}$.
For this choice of strain, the term with $\lambda_{2}$ drops, and we have two
competing couplings - $\lambda_{0}$ and $\lambda_{1}$. The results are
presented in Fig. 6, which displays how the magnetization magnitude varies
with strain angle. For weak strain, the variation is small. Increasing the
strain displays strong oscillations with frequency $6\phi_{S}$, in agreement
with Eq. (49). For even larger strain, higher order oscillations are observed,
and the magnitude saturates. The orientation of the magnetization $\phi_{M}$
is plotted against the strain angle in Fig. 7. A linear phase of $2\phi_{S}$
is added to compensate for the rotation of the orientation due to the first
order contribution. The oscillatory behaviour of the second order is clear
even for weak strains, and the amplitude increases with strain until higher
order effects become significant.

Figure 6: The valley-$K$ resolved magnitude of the magnetization vs strain
angle for the non-interacting TBG with $\theta=1.1^{\circ}$. For very weak
strain the magnitude does not depend on strain. Increasing strain introduces
angular dependence, and further increasing introduces higher order effects.
Figure 7: The valley-$K$ resolved magnetization orientation $\phi_{M}$ shifted
by $2\phi_{S}$ as a function of strain angle $\phi_{S}$, for the non-
interacting TBG with $\theta=1.1^{\circ}$. The oscillatory behaviour is a
result of higher orders in strain.

##  Appendix C Hartree-Fock Approximation

###  C.1 The Hartree-Fock Hamiltonian

We start by considering the complete interaction model. The interacting
Hamiltonian $\hat{H}$ is given by a sum of a single-particle Hamiltonian
$\hat{H}_{0}$ and an interaction term $\hat{V}$:

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{0}+\hat{V},$ (50)
$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\left(H_{\text{BM}}(\mathbf{k})-\Sigma^{HF}[\rho_{0}](\mathbf{k})\right)c_{\mathbf{k}},$
(51) $\displaystyle\hat{V}$
$\displaystyle=\frac{1}{2A}\sum_{\mathbf{q}}V(\mathbf{q}):\hat{n}_{\mathbf{q}}\hat{n}_{-\mathbf{q}}:,$
(52)

where $c_{\mathbf{k}}^{\dagger}$ is a row vector of creation operators
$c_{\mathbf{k},\alpha}^{\dagger}$ for Bloch eigenstates
$\ket{u_{\mathbf{k},\alpha}}$ of $H_{BM}(\mathbf{k})$, which are labeled by
the combined index $\alpha=(\tau,s,m)$ as defined in the main text. The total
area of the system is denoted by $A$.

The single-particle Hamiltonian contains two parts: the BM Hamiltonian,
$H_{BM}(\mathbf{k})$, and the Hartree-Fock self-energy of a reference state
$\rho_{0}$, $\Sigma^{HF}[\rho_{0}](\mathbf{k})$, which is subtracted to avoid
double-counting, as discussed below. The BM Hamiltonian is given for valley
$\tau=K$, by Eq. (8), and for valley $\tau=K^{\prime}$ by its time-reversal.
The parameters in the BM Hamiltonian implicitly include interaction effects,
which necessitates a subtraction scheme in order to avoid double counting. To
this end, we deduct a Hartree-Fock self energy operator
$\Sigma^{HF}[\rho_{0}](\mathbf{k})$ evaluated at a reference state with a
single-particle density matrix $\rho_{0}$ that corresponds to two decoupled
graphene layers Xie and Macdonald (2020, 2021). The resulting $\hat{H}_{0}$ is
designed such that the Hartree-Fock Hamiltonian that corresponds to
$\hat{H}_{0}+\hat{V}$ in the absence of inter-layer tunneling reproduces
$H_{BM}$ with $t_{AA}=t_{AB}=0$. The density matrix $\rho_{0}$ is also
responsible for deducting the background charge density, as it corresponds to
a uniform system at charge neutrality.

The interaction term $\hat{V}$ is chosen to be a double-gated Coulomb
potential, $V(\mathbf{q})=e^{2}/(2\epsilon\epsilon_{0}q)\tanh(qd_{g})$ with
$\epsilon_{0}$ being the vacuum permittivity, $d_{g}=40\rm{nm}$ being the
distance to the gates. The dielectric constant is chosen to be $\epsilon=30$,
for reasons to be explained in section C.2. The operator
$\hat{n}_{\mathbf{q}}$ is the electron density operator with wavevector
$\mathbf{q}$, which can be written in the basis of the BM eigenstates as

$\hat{n}_{\mathbf{q}}=\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\Lambda_{\mathbf{q}}(\mathbf{k})c_{\mathbf{k}+\mathbf{q}},\quad\left[\Lambda_{\mathbf{q}}(\mathbf{k})\right]_{\alpha,\beta}=\braket{u_{\mathbf{k},\alpha}}{u_{\mathbf{k}+\mathbf{q},\beta}}.$
(53)

The Hartree-Fock self energy operator $\hat{\Sigma}^{HF}$ for a given Slater
determinant with a single-particle density matrix $\rho$ is given by:

$\displaystyle\begin{split}\hat{\Sigma}^{HF}[\rho](\mathbf{k})&=\frac{1}{A}\sum_{\mathbf{g}}V(\mathbf{g})c_{\mathbf{k}}^{\dagger}\Lambda_{\mathbf{g}}(\mathbf{k})c_{\mathbf{k}}\sum_{\mathbf{k}^{\prime}}\text{tr}\left[{\rho(\mathbf{k}^{\prime})\Lambda_{\mathbf{g}}^{*}(\mathbf{k}^{\prime})}\right]\\\
&-\frac{1}{A}\sum_{\mathbf{q}=\mathbf{g}+\mathbf{k}^{\prime}}V(\mathbf{q})c_{\mathbf{k}}^{\dagger}\Lambda_{\mathbf{q}}(\mathbf{k})\rho^{T}(\mathbf{k}+\mathbf{q})\Lambda_{\mathbf{q}}^{\dagger}(\mathbf{k})c_{\mathbf{k}},\end{split}$
(54)

where we adopted a periodic gauge, i.e.
$c_{\mathbf{k}+\mathbf{g}}^{\dagger}=c_{\mathbf{k}}^{\dagger}$ for any moiré
reciprocal lattice vector $\mathbf{g}$. The wavevector $\mathbf{k}^{\prime}$
is limited to the first BZ, while $\mathbf{q}$ is any wavevector compatible
with the boundary conditions. The first and second terms in Eq. (54)
correspond to the Hartree and Fock self energies, respectively. In this work
we neglect any flavor coherence, and Bloch wavefunctions are assumed
orthogonal between different flavors. This sets all form factors, interpolated
density matrices and their corresponding self energies to be diagonal in
flavor.

The Hartree-Fock procedure is an approximation method for obtaining the ground
state of (50). The Hartree-Fock Hamiltonian presented in Eq. (23) can be
rewritten in terms of $\hat{H}_{BM}$ and $\hat{\Sigma}^{HF}$ as:

$\hat{H}^{HF}=\hat{H}_{BM}+\hat{\Sigma}^{HF}[\rho-\rho_{0}].$ (55)

The purpose of our interpolation scheme is to provide an approximated form for
$\hat{\Sigma}^{HF}$ without self-consistently determining $\rho$. Each
interpolation configuration in Table 1 fixes the density matrix $\rho(\nu)$ in
terms of the filling, as defined in Eq. (25).

###  C.2 Quasi-particle excitation spectrum

In Fig. 8 we plot the flat bands of the Hartree-Fock spectrum of Eq. (23) for
different interpolation configurations. We compare between the non-interacting
spectrum of the unpolarized configuration at $\nu=0$, and the valley polarized
configurations $\nu=+1$ and $\nu=+3$. We note that the increase in Dirac
velocity is a consequence of the Fock self energy, and the relative shift
between the Dirac points and the $\Gamma$ point is due to the Hartree self
energy, as described in Ref. [Xie and Macdonald, 2021].

Figure 8: Hartree-Fock flat band spectrum of flavor $(K,\uparrow)$ and
$\theta=1.1^{\circ}$ for different interpolation configurations. The non-
interacting spectrum is slightly particle-hole asymmetric, leading to a
deviation of the Dirac points from zero energy. The rest of the plots are
shifted in energy such that at the $\Gamma$ point the two bands sit at exactly
opposite energies.

The dielectric constant $\epsilon$ is treated as a phenomenological parameter
that includes both the screening due to the substrate and the remote bands. We
fix it by requiring that the obtained magic angle roughly matches the
experimental one. We expect a minimal bandwidth for $\theta\sim
1.07^{\circ}-1.1^{\circ}$, and a wider one for $\theta\sim 1.3^{\circ}$. This
is accomplished by our choice of $\epsilon=30$, for which the bandwidths are
presented in the quasi-particle excitation spectrums in Fig. 9. Note that for
this choice the band is actually narrower for $\theta=1.1^{\circ}$ than for
$\theta=1.07^{\circ}$, and that the interaction decreases the bandwidth for
$\theta=1.3^{\circ}$. This figure therefore complements Fig. 5 in the main
text by demonstrating how the bandwidth (or Dirac velocity) is related to the
magnetization response.

Figure 9: Hartree-Fock flat band spectrum for different twist angles. There is
a correspondence between bandwidth (or Dirac velocity) and the sensitivity to
strain, implying the effect of interactions on magnetization is mostly
determined by the Fock term bandwidth renormalization.

##  Appendix D Reduction to a single flavor

As described in the main text, the in-plane magnetization in Eq. (24) for the
valley-polarized states can be computed by focusing on a single flavor block
in the Hartree-Fock interpolation. However, as seen in Eq. (54), the Hartree
term consists of a sum over charge contributions from all the flavors.
Therefore, in the general case, all the flavor blocks in the density matrices
$\rho,\rho_{0}$ are required to obtain the self energy for a specific flavor.
However, using our interpolation scheme, we can simplify the computation and
calculate the Hartree term with a single flavor block.

Consider the following term that appears in the Hartree self energy:

$\braket{n_{-\mathbf{g}}}_{\rho-\rho_{0}}=\sum_{\mathbf{k}^{\prime}}\text{tr}\left[{\left(\rho(\nu)-\rho_{0}(\mathbf{k}^{\prime})\right)\Lambda_{\mathbf{g}}^{*}(\mathbf{k}^{\prime})}\right],$
(56)

where the trace is taken over the spin, valley, and band indices. Since both
the density matrices and the form factors are diagonal in flavor, we can split
the trace into a sum over flavor sectors denoted by a subscript $(\tau,s)$ for
valley and spin:

$\displaystyle\begin{split}\braket{n_{-\mathbf{g}}}_{\rho-\rho_{0}}&=\sum_{\mathbf{k}^{\prime}}\sum_{\tau,s}\text{tr}\left[{\left(\rho_{\tau,s}(\nu_{\tau,s})-\rho_{0;\tau,s}(\mathbf{k}^{\prime})\right)\Lambda_{\mathbf{g}}^{*}(\mathbf{k}^{\prime})}\right]\\\
&=\sum_{\tau,s}\left(\braket{n_{-\mathbf{g}}}_{\rho;(\tau,s)}-\braket{n_{-\mathbf{g}}}_{\rho_{0};(\tau,s)}\right),\end{split}$
(57)

where here the trace is over the band indices.

For the density matrix $\rho_{0;\tau,s}(\mathbf{k})$ the expression above is
simple: there is no charge at any nonzero wavevector, and all flavors are
symmetric,
$\braket{n_{-\mathbf{g}}}_{\rho_{0};(\tau,s)}=n_{0}\delta_{\mathbf{g}=\mathbf{0}}$.
For the interpolated density matrix $\rho_{\tau,s}(\nu_{\tau,s})$, each flavor
is a combination of density matrices of filled and empty flavors. Suppose that
all flavors are filled to $\nu_{\tau,s}=\nu_{0}$, such that the flavor blocks
in the density matrix are related by symmetry. The spin $SU(2)$ symmetry
implies that
$\braket{n_{-\mathbf{g}}}_{\rho;(\tau,\uparrow)}=\braket{n_{-\mathbf{g}}}_{\rho;(\tau,\downarrow)}$.
The $C_{2}$ symmetry (which holds in the presence of strain) implies that
$\braket{n_{-\mathbf{g}}}_{\rho;(K,s)}=\braket{n_{+\mathbf{g}}}_{\rho;(K^{\prime},s)}$.
Together with $\hat{n}_{-\mathbf{g}}=\hat{n}_{+\mathbf{g}}^{\dagger}$ we find
that
$\braket{n_{-\mathbf{g}}}_{\rho;(K,s)}=\braket{n_{-\mathbf{g}}}_{\rho;(K^{\prime},s)}^{*}$.
The last symmetry we invoke is the anti-unitary
$C_{2}\mathcal{T}=\sigma_{x}\mathcal{K}$, where $\sigma_{x}$ flips the
sublattice and $\mathcal{K}$ is the complex conjugation operator. This
symmetry implies that
$\ket{u_{\mathbf{k},\alpha}}=e^{i\phi_{\mathbf{k},\alpha}}\sigma_{x}(\ket{u_{\mathbf{k},\alpha}})^{*}$,
with a gauge dependent phase $\phi_{\mathbf{k},\alpha}$. In the periodic gauge
it is easy to show that we have
$\phi_{\mathbf{k}+\mathbf{g},\alpha}=\phi_{\mathbf{k},\alpha}$ for any lattice
vector $\mathbf{g}$. This symmetry immediately implies that the diagonal
elements of the form factors,
$\left[\Lambda_{\mathbf{g}}(\mathbf{k}^{\prime})\right]_{\alpha,\alpha}$, are
strictly real. Since $\rho_{\tau,s}(\nu_{\tau,s})$ is diagonal as well, we
find that $\braket{n_{-\mathbf{g}}}_{\rho;(\tau,s)}$ is strictly real, hence
$\braket{n_{-\mathbf{g}}}_{\rho;(K,s)}=\braket{n_{-\mathbf{g}}}_{\rho;(K^{\prime},s)}$.

Combining the arguments above, we find that $\rho_{K,\uparrow}$ and
$\rho_{0;K,\uparrow}$ can be used to calculate any flavor’s contribution to
the overall charge density by using them with each flavor’s respective
filling:

$\displaystyle\begin{split}\braket{n_{-\mathbf{g}}}_{\rho-\rho_{0}}&=\sum_{\mathbf{k}^{\prime}}\sum_{\tau,s}\text{tr}\left[{\left(\rho_{K,\uparrow}(\nu_{\tau,s})-\rho_{0;K,\uparrow}(\mathbf{k}^{\prime})\right)\Lambda_{\mathbf{g}}^{*}(\mathbf{k}^{\prime})}\right].\end{split}$
(58)

Eq. (58) implies that the Hartree-Fock self energy for the flavor
$(\tau,s)=(K,\uparrow)$ in Eq. (55) can be written solely with
$\rho_{K,\uparrow}$ and $\rho_{0;K,\uparrow}$ by using the following effective
density matrix for the Hartree and Fock terms separately:

$\displaystyle\begin{split}\Sigma^{HF}[\rho-\rho_{0}]_{K,\uparrow}&=\Sigma^{H}\left[\left(\Sigma_{\tau,s}\rho_{K,\uparrow}(\nu_{\tau,s})\right)-4\rho_{0;K,\uparrow}(\mathbf{k}^{\prime})\right]\\\
&+\Sigma^{F}[\rho_{K,\uparrow}(\nu_{K,\uparrow})-\rho_{0;K,\uparrow}(\mathbf{k}^{\prime})].\end{split}$
(59)

We therefore restrict all density matrices in Appendix E to be in the
$K\uparrow$ flavor block, and solve the Hartree-Fock Hamiltonian for this
flavor.

##  Appendix E Numerical procedure for calculating the chemical potential

The purpose of this appendix is to detail technical aspects involving
calculating the chemical potential as a function of filling and its gradient
with respect to magnetic field.

We start by detailing the numerical parameters used in the computation. The
non-interacting BM Hamiltonian of flavor $(\tau,s)=(K,\uparrow)$ is
constructed symmetrically around the $\Gamma$ point in the moiré BZ, with the
momentum space cutoff taken to include all hopping sites within a hexagon of
edge size of $3\sqrt{3}[K\theta]$, i.e. an edge size of $R_{BZ}=3$ in units of
a reciprocal lattice wavevector magnitude. This cutoff includes 27 hopping
sites in each layer. Multiplying by the sublattice degree of freedom, we
obtain $N_{\text{bands}}=108$ bands for each quasimomentum point in the moiré
BZ. Note that, by construction, this Hamiltonian is not exactly $C_{3}$
symmetric about any point beside the $\Gamma$ point. The quasimomentum
sampling grid within the first BZ is centered around the $\Gamma$ point with
the number of unit cells taken to be $N_{u.c.}=133$. The grid is sampled first
along the unstrained moiré lattice vectors, and then transformed under the
desired strain, in order to preserve $N_{u.c.}$ across different strains.

The chemical potential for a particular filling factor $\mu(\nu)$ is
calculated from the spectrum of the HF Hamiltonian in Eq. (23) using the
procedure outlined in the main text. We will now describe this procedure with
further technical details. The first step is to use the reduction scheme in
Eq. (59) to restrict our attention of the HF Hamiltonian to a single flavor
block. Next, we simplify the diagonalization of this Hamiltonian by projecting
the self energy $\hat{\Sigma}^{HF}[\rho](\mathbf{k})$ to the two flat bands in
the non-interacting Hamiltonian, as in Ref. [Xie and Macdonald, 2021]. We
emphasize that all the remote bands are accounted for by considering the
density matrix $\rho(\mathbf{k})$ and the form factors
$\Lambda_{\mathbf{g}}(\mathbf{k})$ to be the entire $N_{\text{bands}}\times
N_{\text{bands}}$ matrices. The flat-band projection truncates the Hartree
Fock operator in Eq. (23) to be a $2\times 2$ matrix for each quasimomentum.

After projection, the two Hartree-Fock bands consist of $2N_{u.c.}$ energy
eigenvalues from all quasimomenta in the moiré BZ combined. The filling factor
range $\nu\in(\nu_{e},\nu_{f})$ is discretized accordingly to $\nu(m)$ where
$m=1,2,\dots,2N_{u.c.}$. Starting with the lowest filling $\nu(1)=\nu_{e}$, we
compute the sorted energy eigenvalues $\epsilon_{i}(\nu(1))$,
$i=1,2,\dots,2N_{u.c.}$, such that $\epsilon_{i}<\epsilon_{i+1}$. The chemical
potential $\mu(\nu=\nu(1))$ is given by choosing the lowest energy eigenvalue
$\epsilon_{1}$. We iterate this procedure over $m$ by computing the density
matrix at the next filling factor $\nu(m)$, computing its corresponding
projected HF Hamiltonian, sorting its eigenvalues and selecting the new $m$-th
energy $\epsilon_{m}$. In summary, the chemical potential $\mu(\nu)$ is a list
of $2N_{u.c.}$ values, whose $m$-th entry corresponds to the $m$-th eigenvalue
of a HF Hamiltonian interpolated to the $m$-th filling. We note that the
information about the specific value of $\mathbf{k}$ which is used for
calculating $\mu(\nu)$ is lost.

For the purpose of calculating the magnetization we computed the chemical
potential at four magnetic fields,
$\mathbf{B}_{\parallel}\in\\{-B\hat{x},-B\hat{y},B\hat{x},B\hat{y}\\}$, and
extracted the gradient of the chemical potential
$\mu(\nu,\mathbf{B}_{\parallel})$ at $\mathbf{B}_{\parallel}=0$. Next, one
might expect that integrating over filling would yield the desired
magnetization, as Eq. (24) requires. Unfortunately, there is one critical
problem with immediately performing this integration. The problem is most
clearly explained for TBG without strain. In the absence of strain, we expect
the magnetization to be exactly zero as a function of filling. Since quasi-
particle states are filled one at a time, the $C_{3}$ symmetry is explicitly
broken. We therefore find that the magnetization per moiré unit cell varies as
a function of filling, reaching magnitudes of the order of $1\mu_{B}/N_{u.c.}$
that are cancelled only after filling all three $C_{3}$-symmetric states. The
problem persists at weak strains, but becomes irrelevant once all single-
particle magnetization contributions are aligned.

To remedy this we introduce a finite temperature to smooth out the chemical
potential. Given a filling $\nu$ and temperature $T$, the chemical potential
is set to fix the Fermi-Dirac distribution $f_{FD}$ to the number of particles
$N_{\text{p}}(\nu)$:

$N_{\text{p}}(\nu)=\sum_{\mathbf{k}}f_{FD}\left(\epsilon_{\mathbf{k}}(\nu,\mathbf{B}_{\parallel})-\mu(\nu,\mathbf{B}_{\parallel})\right),$
(60)

where $\epsilon_{\mathbf{k}}(\nu,\mathbf{B}_{\parallel})$ are the HF
eigenvalues for the density matrix corresponding to filling $\nu$, evaluated
in the presence of a magnetic field $\mathbf{B}_{\parallel}$. Taking the
gradient with respect to magnetic field, we find the self consistent equation
for $\mu(\nu,\mathbf{B}_{\parallel})$:

$\frac{\partial\mu}{\partial\mathbf{B}_{\parallel}}=\left(\sum_{\mathbf{k}}\frac{\partial
f_{FD}}{\partial\epsilon_{\mathbf{k}}}\right)^{-1}\sum_{\mathbf{k}}\frac{\partial
f_{FD}}{\partial\epsilon_{\mathbf{k}}}\frac{\partial\epsilon_{\mathbf{k}}}{\partial\mathbf{B}_{\parallel}}.$
(61)

We use the following equation as a definition for the “smoothed” chemical
potential gradient at $\mathbf{B}_{\parallel}=0$, and evaluate it using the
chemical potential $\mu(\nu,\mathbf{B}_{\parallel}=0)$. This chemical
potential gradient is the one used for all the numerical results. The specific
temperature $T$ varies for each plot but is of the order of average level
spacing, i.e. $T\sim E_{BW}/N_{u.c.}\sim 0.1\rm{meV}$, where $E_{BW}$ is the
combined bandwidth of the HF flat bands.

We have verified that our results for magnetization magnitudes
$|M_{\parallel}|\gtrsim 0.1\mu_{B}$ are sufficiently converged for our choice
of $R_{BZ}$, $N_{u.c.}$, and stable in a range of temperatures $T$. For lower
magnetization values we have recovered the linear dependence of the
magnetization vs. strain, but the precise slope has a small dependence on the
numerical parameters.