# Strain Disorder and Gapless Intervalley Coherent Phase in Twisted Bilayer
Graphene

Gal Shavit  Department of Condensed Matter Physics, Weizmann Institute of
Science, Rehovot, Israel 7610001  Kryštof Kolář  Dahlem Center for Complex
Quantum Systems and Fachbereich Physik, Freie Universität Berlin, 14195
Berlin, Germany Christophe Mora  Université Paris Cité, CNRS, Laboratoire
Matériaux et Phénomènes Quantiques, 75013 Paris, France  Felix von Oppen
Dahlem Center for Complex Quantum Systems and Fachbereich Physik, Freie
Universität Berlin, 14195 Berlin, Germany Yuval Oreg  Department of Condensed
Matter Physics, Weizmann Institute of Science, Rehovot, Israel 7610001

###### Abstract

Correlated insulators are frequently observed in magic angle twisted bilayer
graphene at even fillings of electrons or holes per moiré unit-cell. Whereas
theory predicts these insulators to be intervalley coherent excitonic phases,
the measured gaps are routinely much smaller than theoretical estimates. We
explore the effects of random strain variations on the intervalley coherent
phase, which have a pair-breaking effect analogous to magnetic disorder in
superconductors. We find that the spectral gap may be strongly suppressed by
strain disorder, or vanish altogether, even as intervalley coherence is
maintained. We discuss predicted features of the tunneling density of states,
show that the activation gap measured in transport experiments corresponds to
the diminished gap, and thus offer a solution for the apparent discrepancy
between the theoretical and experimental gaps.

Introduction.—In recent years, magic-angle twisted bilayer graphene (MATBG)
[1] has emerged as an exciting and versatile platform for strong-correlation
physics. Its rich phase diagram prominently features correlated insulators,
Chern insulators, superconductivity, and strange metallicity [2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12].

Correlated insulators consistently found at fillings $\nu=\pm 2$ electrons per
moiré unit cell relative to charge neutrality (CN), and more rarely at CN,
have been suggested to originate from Kramers intervalley-coherent order
(K-IVC) [13, 14, 15, 16]. The K-IVC phase is in some sense analogous to a
superconductor, with the condensate made up of intervalley electron-holes
pairs. Interestingly, theoretical predictions of the K-IVC gap energy, by both
Hartree-Fock and numerically exact methods, overestimate the measured
activation gap of $\sim{\cal O}\left(1\,{\rm meV}\right)$ [2, 4, 5, 6, 9] by
more than an order of magnitude.

It has been experimentally well-established that the moiré lattice formed in
realistic MATBG devices is not pristine, presumably due to substantial
relaxation effects of the underlying graphene lattice [17]. Disorder in the
local twist angle between the graphene layers has been observed [18, 7],
leading to domains with slightly different effective moiré unit-cell sizes.
Local measurements have also shown significant strain effects, consistent with
a moiré lattice distortion of $0.1\%-0.7\%$ [19, 20, 21].

As the two graphene layers may be subjected to different strain fields, the
strain tensor applied to the bilayer is comprised of a layer-symmetric part
(homostrain) and a layer-antisymmetric contribution (heterostrain). Uniform
heterostrain was suggested to have an important role in weakening the even-
filling correlated insulator states in MATBG [13, 22]. This is mainly due to a
large increase of the flat-bands bandwidth [23, 24], leading to diminished
effective interactions.

Whereas heterostrain shifts the two Dirac points within each valley with
respect to each other, homostrain subjects both layers to identical
distortions and mainly acts as a pseudo-gauge field. This field acts
oppositely in the two valleys (as illustrated in Fig. 1a), thus maintaining
time-reversal symmetry (TRS). In this manuscript, we explore the effects of
spatially random homostrain (Fig. 1b) on the properties of the K-IVC phase. We
show that this perturbation induces “pair-breaking” effects in striking
resemblance to magnetic impurities in spin-singlet superconductors. Effects of
other types of strain on the K-IVC phase are discussed in the supplementary
materials (SM) [25], and various disorder perturbations are classified by
their impact on this phase in Ref. [26].

We find that modest strain fluctuations reduces both the K-IVC order parameter
and the spectral gap, but enable the two to be dramatically different from one
another. We propose that this effect may be responsible for the surprisingly
small activation gap observed in transport experiments. Moreover, we show that
intervalley coherence can persist even when the system becomes gapless to
single-particle excitations, a phase that we dub “gapless K-IVC”. This phase
could explain the haphazard appearance of a correlated insulator at CN.

Figure 1:  Schematic description of our model. (a) The band structure in each
valley ($K,K^{\prime}$) is approximated by two Dirac cones, one in each “mini-
valley”. The momentum separating the valleys, $\mathbf{Q}$, is modified by
random homostrain fluctuations acting on the graphene layers, as represented
by colored dotted arrows. (b) Schematic of a possible random strain
configuration. Color indicates the strain strength $\left|{\mathbfcal
A}\right|$, whereas arrows indicate the directions of the local distortions of
the graphene lattice.

Model.—We consider the following simplified spinless model of MATBG,

$H_{0}=\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}\left[uk_{x}\sigma_{x}\tau_{z}+uk_{y}\sigma_{y}\right]\Psi^{\phantom{\dagger}}_{\mathbf{k}},$
(1)

where $\Psi^{\dagger}_{\mathbf{k}}$ is an 8-spinor of fermionic annihilation
operators at momentum $\mathbf{k}$ relative to their respective origin in
momentum space, and with sublattice, valley, and ”mini-valley” degrees of
freedom, described by Pauli matrices $\sigma_{i}$, $\tau_{i}$, and $\rho_{i}$,
respectively. This model features four Dirac cones with linear dispersion
$\epsilon_{\mathbf{k}}=\pm u\left|{\mathbf{k}}\right|$. $H_{0}$ preserves the
TRS ${\cal T}=\tau_{x}\rho_{x}{\cal K}$ (${\cal K}$ implements complex
conjugation), with ${\cal T}^{2}=1$.

The spinless model effectively describes the fillings $\nu=\pm 2$, where two
electronic flavors of opposite valleys are entirely filled or empty, while the
remaining pair is “active” and forms intervalley coherence. Conversely, two
copies of our spinless model may be used describe the K-IVC phase at CN [13,
14, 15, 16].

We introduce local density-density repulsive interactions,

$H_{\rm
int}=\frac{U}{2\Omega}\sum_{{\mathbf{k}},{\mathbf{k}^{\prime}},{\mathbf{q}}}\Psi^{\dagger}_{\mathbf{k+q}}\Psi^{\phantom{\dagger}}_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k^{\prime}-q}}\Psi^{\phantom{\dagger}}_{\mathbf{k^{\prime}}},$
(2)

which induce spontaneous symmetry-breaking in our model on a mean-field level.
($\Omega$ is the system area.) This simplified form of $H_{\rm int}$ suffices
to illustrate the phenomenon we are interested in. Specifically, we examine
the K-IVC phase, which has been argued to be the likely ground state of MATBG
at even fillings. The K-IVC state is characterized by intervalley coherence,
i.e., formation of an exciton condensate with intervalley electron-hole pairs,
a finite gap to charge excitations, and TRS breaking. The corresponding mean-
field Hamiltonian at a given $\mathbf{k}$ has the form

$h_{\rm
MF}\left({\mathbf{k}}\right)=u\left(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\right)+\Delta_{\rm
ivc}\sigma_{x}\tau_{x}\rho_{z}.$ (3)

$h_{\rm MF}$ preserves a Kramers-like TRS, ${\cal
T}^{\prime}=\tau_{y}\rho_{x}{\cal K}$, which concatenates $\mathcal{T}$ with a
valley rotation. It also preserves chiral and particle-hole symmetries,
represented by ${\cal S}=\sigma_{z}$ and ${\cal C}={\cal S}{\cal T}^{\prime}$,
respectively 111We note that a more realistic model of the band structure of
MATBG contains terms which break particle-hole symmetry, as well as terms
proportional to $\tau_{z}$. These break $\cal S$ and $\cal C$, reducing the
symmetry of $h_{\rm MF}$ from symmetry class CII to AII. .

We now introduce random homostrain variations, which enter the model as a
random gauge field acting in opposite directions in opposite valleys, see Fig.
1b. Concretely, the strain Hamiltonian may be written as

$H_{\rm str}=u\int
d{\mathbf{r}}\Psi^{\dagger}\left({\mathbf{r}}\right)\left[{\cal
A}_{x}\left({\mathbf{r}}\right)\sigma_{x}+{\cal
A}_{y}\left({\mathbf{r}}\right)\sigma_{y}\tau_{z}\right]\Psi\left({\mathbf{r}}\right),$
(4)

where $\mathbfcal{A}$ is the strain-induced perturbation, and
$\Psi\left({\mathbf{r}}\right)=\frac{1}{\sqrt{\Omega}}\sum_{\mathbf{k}}e^{i{\mathbf{k}}\cdot{\mathbf{r}}}\Psi_{\mathbf{k}}$.
In Eq. (4) we have implicitly assumed that the strain potential is smooth on
the moiré length scale. Otherwise one should replace the renormalized velocity
$u$ by the much larger bare graphene Fermi velocity $v_{F}$. Notice that
uniform strain constitutes a shift of $\mathbf{k}$ in $H_{0}$ [28] and
redefines the momentum connecting the two valleys $\mathbf{Q=K-K^{\prime}}$.
This only changes the momentum carried by the condensed electron-hole K-IVC
pairs, so that we can assume ${\cal A}_{i}$ has zero spatial mean.

The perturbation in Eq. (4) does not break any of the symmetries of $h_{\rm
MF}$. However, as we will show below (see also Ref. [26]), the fact that it
commutes with the K-IVC operator $\sigma_{x}\tau_{x}\rho_{z}$, enables a
drastic reduction of the gap which opens in the K-IVC spectrum due to the
random strain disorder. For this reason, we have also neglected inter-
minivalley scattering, which have additional $\rho_{x},\rho_{y}$ factors, and
thus anticommute with the order parameter.

Considering a simple point-like perturbation,
$\mathbfcal{A}\left({\mathbf{r}}\right)=\mathbfcal{A}_{0}\delta\left({\mathbf{r}}\right)$,
one may employ a $T$-matrix formalism to find the bound-state spectrum inside
the mean-field gap [25], in analogy to Yu-Shiba-Rusinov states induced by
magnetic impurities in a superconductor [29, 30, 31]. We find that as the
perturbation strength increases, the in-gap-state energy is reduced and the
two bound-state energies cross at zero when the impurity strength becomes an
appreciable fraction of the bandwidth $W$.

Moreover, we find that when approximating the density-of-states (DOS) as a
constant around the Fermi level (rather than linear as appropriate in our
case), one recovers – apart from additional degeneracies – precisely the
bound-state spectrum of a magnetic impurity inside a singlet superconductor.
This outcome may be traced to the analogous algebraic structure of the two
problems, i.e., the impurity operator commuting with the order parameter. This
analogy enables the treatment of random strain fluctuations in MATBG by tools
similar to those employed in superconductors with magnetic impurities.

Abrikosov-Gor’kov approach.—We therefore treat the random strain fluctuations
within the self-consistent Born approximation (SCBA), inspired by the
Abrikosov-Gor’kov theory of superconductivity in magnetically disordered
alloys [32]. A similar method was also used to study exciton condensates in
the presence of potential impurities [33, 34].

The main object of interest is the Green’s function,

$G\left({\mathbf{k}},\omega\right)=\left(i\omega-h_{\rm MF}-\hat{\Sigma}_{\rm
SCBA}\right)^{-1}.$ (5)

Within the SCBA, the self-energy matrix $\hat{\Sigma}_{\rm SCBA}$ can be
written as

$\hat{\Sigma}_{\rm
SCBA}\left({\mathbf{k}},\omega\right)=\left\langle\sum_{\mathbf{p}}{\cal
U}_{{\mathbf{k}}-{\mathbf{p}}}G\left({\mathbf{p}},\omega\right){\cal
U}_{{\mathbf{p}}-{\mathbf{k}}}\right\rangle_{\rm dis.},$ (6)

where the matrix $\cal U$ represents the random strain perturbation in
momentum space, $H_{\rm
str}=\sum_{\mathbf{kq}}\Psi^{\dagger}_{\mathbf{k+q}}{\cal
U}_{\mathbf{q}}\Psi_{\mathbf{k}}$, and $\left\langle...\right\rangle_{\rm
dis.}$ stands for disorder averaging.

Upon standard manipulation of the Green’s function, it can be written as

$G=-\frac{i\tilde{\omega}+u\left(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\right)+\tilde{\Delta}\sigma_{x}\tau_{x}\rho_{z}}{\epsilon_{\mathbf{k}}^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}}.$
(7)

The parameters $\tilde{\omega},\tilde{\Delta}$ are related to
${\omega},{\Delta_{\rm ivc}}$ by the self-consistency equations

$\begin{pmatrix}\tilde{\omega}\\\
\tilde{\Delta}\end{pmatrix}=\begin{pmatrix}\omega\\\ \Delta_{{\rm
ivc}}\end{pmatrix}+\Gamma{\cal
F}\left(\tilde{\omega},\tilde{\Delta}\right)\begin{pmatrix}\tilde{\omega}\\\
-\tilde{\Delta}\end{pmatrix},$ (8)

where ${\cal F}\left(\tilde{\omega},\tilde{\Delta}\right)=\frac{W}{2}\int
d\epsilon\frac{{\cal
N}\left(\epsilon\right)}{\epsilon^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}}$,
${\cal N}\left(\epsilon\right)$ is the DOS per unit cell of area
$\Omega_{u.c.}$, and $W$ is the bandwidth. The disorder energy scale $\Gamma$
is

$\Gamma=\frac{2\Omega}{W\Omega_{u.c.}}\left\langle\int d\theta
u^{2}\mathbf{A}_{\theta}^{\dagger}\cdot\mathbf{A}_{\theta}^{\phantom{\dagger}}\right\rangle_{\rm
dis.},$ (9)

where $\mathbf{A}_{\mathbf{q}}$ is the Fourier transform of
$\mathbfcal{A}\left(\mathbf{r}\right)$, and we have used the standard
approximation that the scattering mostly depends on the angle between incoming
and outgoing momenta $\theta$ [33, 34]. As for magnetic impurities in
superconductors and potential scatterers in excitonic condensates, the form of
Eq. (8) is due to the K-IVC order parameter commuting with the random
perturbation. Thus, the equations we find are identical to the Abrikosov-
Gor’kov equations for superconductors with magnetic impurities, with one
important difference. In MATBG, we cannot assume a constant DOS near the Fermi
energy, but should account for the fact that the DOS vanishes linearly at the
Dirac point. This leads to important qualitative differences in the results.

By relating the local strain to the effective gauge-field $\mathbfcal{A}$, one
may obtain an order-of-magnitude estimate of $\Gamma$. For root-mean-square
strains of ${\cal E}\sim 0.1\%$ and disorder correlation lengths of few unit
cells, we find $\Gamma/W$ values of $0.1-0.3$ [25]. As will be shown, such
values are sufficient to dramatically reduce the spectral gap or even close it
completely.

The strength of K-IVC order in the presence of disorder is obtained by
combining Eq. (8) with the gap equation

$\Delta_{\rm ivc}=-2\frac{U}{\beta\Omega}\sum_{\omega{\mathbf{k}}}{\rm
Tr}\left\\{\sigma_{x}\tau_{x}\rho_{z}G\right\\},$ (10)

which we solve numerically. (Here, $\beta$ is inverse temperature.) Figure 2a
shows results for the order parameter $\Delta_{\rm ivc}$ as a function of
temperature and disorder. We find that both the order parameter and the
critical temperature deteriorate with increasing disorder.

Figure 2:  (a) The order parameter $\Delta_{\rm ivc}$ as a function of
temperature and disorder energy scale $\Gamma$. For a given temperature, there
is a critical $\Gamma$ at which the order parameter vanishes. (b) The ratio
between the gap in the TDOS, $\omega_{g}$, and $\Delta_{\rm ivc}$. As $\Gamma$
increases, $\omega_{g}$ deviates from $\Delta_{\rm ivc}$ appreciably. There
exists a gapless K-IVC region bounded by the vanishing of the single-particle
gap (black dashed line) and of the order parameter (blue dashed line). We used
$U=0.7W$ and a linear DOS ${\cal N}_{\rm lin.}$. (c) Linecuts for different
values of $\Gamma$ showing $\Delta_{\rm ivc}$ (solid) and $\omega_{g}$ (dash-
dotted).

When assuming a DOS which is linear in energy, ${\cal N}_{\rm
lin.}=4\left|\epsilon\right|/W^{2}$ with cutoff energy
$\left|\epsilon\right|<W/2$, we can also make analytical progress. In
particular, we may calculate the critical disorder scale $\Gamma_{c}$, at
which $\Delta_{\rm ivc}$ vanishes for $T=0$. In the regime $\Gamma_{c}\ll
W/2$, we find [25]

$\Gamma_{c}=\frac{U_{c}}{\log 8}\left(1-\frac{U_{c}}{U}\right),$ (11)

where $U_{c}=W/4$ is the critical interaction $U$ below which the K-IVC order
vanishes at $\Gamma,T=0$. The finite $U_{c}$ as well as the form of the
critical “pair-breaking” parameter $\Gamma_{c}$ originate in the DOS vanishing
linearly at zero. This suppresses the analog of the Cooper-instability for
arbitrarily weak interactions, which requires a finite DOS at the Fermi level.

Having found the self-consistent Green’s function, we can calculate the
tunneling density of states (TDOS) $\tilde{\cal
N}\left(\epsilon\right)=\frac{1}{\pi}{\rm
Im}\frac{1}{\Omega}\sum_{\mathbf{k}}{\rm Tr}G\left({\mathbf{k}},\omega\to
i\epsilon\right)$. In clean systems, the gap $\omega_{g}$ in the TDOS is equal
to the order parameter $\Delta_{\rm ivc}$. In the presence of finite disorder
$\Gamma$, $\omega_{g}$ is in general smaller than $\Delta_{\rm ivc}$.

In Fig. 2b, we plot the ratio $\omega_{g}/\Delta_{\rm ivc}$ for the same
parameter range as in Fig. 2a. As $\Gamma$ increases, the ratio gradually
deviates from unity, and eventually reaches zero, before $\Delta_{\rm ivc}$
vanishes. We dub the regime with finite $\Delta_{\mathrm{ivc}}$ and vanishing
$\omega_{g}$ as the gapless K-IVC phase. Similar to gapless superconductivity,
we interpret this regime as one where intervalley coherence exists throughout
a large fraction of the system, yet strain-induced in-gap states form a low-
energy compressible continuum of states. For linear DOS, we find that the
disorder strength at which the gap closes is related to $\Delta_{\rm ivc}$
[25] through $\Gamma_{g}=W/\log\left(1+\frac{W^{2}}{\Delta_{{\rm
ivc}}^{2}}\right)^{2}$.

In Fig. 2c, we plot linecuts of Fig. 2a for several values of $\Gamma$. We
find that even for modest (and realistic) disorder strengths, the spectral gap
is significantly suppressed compared to the naive gap as given by the order
parameter. Thus, devices which host appreciable strain fluctuations may
exhibit an effective gap as seen in a global transport measurement, which is
about an order of magnitude smaller than the expected condensation order
parameter.

Figure 3:  Evolution of the TDOS $\tilde{\cal N}(\epsilon)$ with $\Gamma$. The
gap $\omega_{g}$ gradually closes as the disorder strength increases. The blue
line indicates the evolution of $\Delta_{\rm ivc}/W$ with increasing $\Gamma$.
Notice that its evolution differs from that of the single-particle gap as it
deteriorates more slowly. We used $U=0.7W$, $T=0$ in this figure.

In Fig. 3 we track the evolution of the TDOS with disorder strength. At very
low $\Gamma$ we find two narrow K-IVC bands separated by roughly $2\Delta_{\rm
ivc}$, as expected from previous theoretical investigations of the pristine
K-IVC phase. As the disorder strength $\Gamma$ increases, these bands spread
out in energy, and their separation diminishes. As shown in the figure, this
is very different from the behavior of $\Delta_{\rm ivc}$. While $\omega_{g}$
diminishes already for weak disorder, $\Delta_{\rm ivc}$ remains mostly
unaffected up to intermediate values of $\Gamma$.

The TDOS depicted in Fig. 3 can be measured in planar tunneling junctions with
a large tunneling area, similar to the experimental verification of gapless
superconductivity [35]. The large tunneling area is required for effective
averaging over disorder configurations in a particular device. Such
measurements are expected to show two spread-out TDOS lobes, with centers
separated by $\sim 2\Delta_{\rm ivc}$ and a spectral gap of
$2\omega_{g}<2\Delta_{\rm ivc}$. In contrast, local scanning-tunneling-
microscopy (STM) measurements do not reveal disorder-averaged quantities, yet
we expect the tunneling gap to vary as a function of position, in a manner
correlated with the homostrain variations. Such phenomenology would be a clear
indicator of the proposed disordered K-IVC phase.

Figure 4:  Arrhenius plots of the longitudinal DC conductivity $\sigma_{xx}$.
Dots: $\sigma_{xx}$ as a function of temperature for different strengths of
the strain fluctuations, $\Gamma/W\in\left[0.02,0.1\right]$, increasing from
bottom (red) to top (blue) in steps of $0.02$. We used the same parameters as
in Fig. 2. Solid lines: guides to the eye
$\propto\exp\left[-\omega_{g}\left(T\to 0\right)/T\right]$ for each $\Gamma$
value. Whereas $\omega_{g}$ differs by a factor of $\sim 5$ between the first
and last plots, $\Delta_{\rm ivc}$ changes only modestly by $<25\%$.

To make contact with transport experiments, we now turn to calculating the DC
conductivity within our model. We use the Kubo-Bastin formula [36],
$\sigma_{xx}\propto\int d\epsilon\left(-\frac{df}{d\epsilon}\right){\cal
S}\left(\epsilon\right)$, where $f$ is the Fermi function, and

${\cal S}\left(\epsilon\right)\approx\frac{1}{\Omega}\sum_{\mathbf{k}}{\rm
Tr}\left\\{j_{x}{\rm Im}G\left({\mathbf{k}},\omega\to
i\epsilon\right)j_{x}{\rm Im}G\left({\mathbf{k}},\omega\to
i\epsilon\right)\right\\}$ (12)

is the conductivity kernel. The current operator is
$j_{x}=u\sigma_{x}\tau_{z}$. Equation (12) neglects vertex corrections to
$j_{x}$, which may have quantitative significance, but are beyond the scope of
this work.

Fig. 4 presents Arrhenius plots of the conductivity for different $\Gamma$
values. From these plots, we can extract the activation energy $E_{\rm act}$
by fitting the longitudinal conductivity to
$\sigma_{xx}\propto\exp\left(-E_{\rm act}/T\right)$. Remarkably, we observe
that the low-$T$ behavior is indeed temperature activated with $E_{\rm
act}\approx\omega_{g}$, and not the potentially much larger $\Delta_{\rm
ivc}$. (In the topmost plot of Fig. 4 $\Delta_{\rm ivc}\approx
5.4\omega_{g}$.) This behavior can be traced to the analytic structure of
${\cal S}\left(\epsilon\right)$, which has a gap $\approx 2\omega_{g}$ around
$\epsilon=0$ [25], similar to the TDOS.

Conclusions.—We have explored the consequences of random homostrain
fluctuations on the K-IVC state, believed to describe the insulating phases of
MATBG at even fillings. Using a simplified model for MATBG, we have studied
this problem using the SCBA in conjunction with a mean-field treatment of the
K-IVC order parameter. Homostrain disorder has a pair-breaking effect on the
intervalley coherent condensate, since it locally acts on the two valleys in
opposite ways.

In contrast to similar pair-breaking disorder problems, random homostrain does
not break any symmetries of the K-IVC state. However, it does lead to in-gap
states, gap closing, and order parameter deterioration due to its operator
structure – it commutes with the order parameter. Moreover, the DOS dependence
on energy had to be taken into account, since it vanishes at the Dirac point.
This led to the unique form of the solutions of the Abrikosov-Gor’kov
equations which we derive, and of the critical pair-breaking parameter
$\Gamma_{c}$.

One of our key results is the significant separation between the energy scales
of the K-IVC order ($\Delta_{\mathrm{ivc}}$) and the spectral gap for single-
particle excitations ($\omega_{g}$), even for modest values of disorder.
Borrowing insights from superconductors, the gap reduction stems from in-gap
bound states, which become stronger and more abundant with increasing
$\Gamma$, yet impact the surrounding intervalley-coherent condensate only
weakly.

We suggest that the order-of-magnitude discrepancy between the theorized K-IVC
gap and the activation gap measured in transport experiments can be resolved
within our model. We have demonstrated that the relevant activation energy as
measured via the DC conductivity is the spectral gap, which may be
considerably smaller than the order parameter due to disorder. (Both scales
coincide in the pristine case). The rare appearance of insulators at CN can
also be understood by considering two copies of our model with different spin
labels. Variations of the magnitude of strain disorder between devices may tip
the state at CN from a weakly-insulating K-IVC state to the gapless K-IVC
regime. The relative weakness of the insulating state at CN compared to
fillings $\nu=\pm 2$ has been attributed to bandwidth renormalizations [37,
38, 39], rendering the effective interactions stronger away from CN.

The interplay of strain fluctuations with other sources of K-IVC suppression,
such as twist-angle disorder and uniform heterostrain, remains to be explored.
Additionally, the fact that the considered disorder couples only to
intervalley-ordered states may also be important. This may have ramifications
for the competition between the K-IVC state and other correlated insulating
phases, such as the valley-Hall phase, for which the order parameter
$\propto\sigma_{z}$ anticommutes with the strain fluctuations. Incorporating
such complications, as well as including more intricate aspects of the
(particle-hole asymmetric) band structure will shed much-needed light on the
nature of the insulating phases in MATBG and their variation between different
devices.

###### Acknowledgements.

We gratefully acknowledge funding by Deutsche Forschungsgemeinschaft through
CRC 183 (project C02; YO and FvO), a joint ANR-DFG project (TWISTGRAPH; CM and
FvO), the European Union’s Horizon 2020 research and innovation program (Grant
Agreement LEGOTOP No. 788715; YO), the ISF Quantum Science and Technology
program (2074/19; YO), and a BSF and NSF grant (2018643; YO).

## References

  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proc. Nat. Acad. Sci. 108, 12233 (2011).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Phys. 16, 926 (2020).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Nature 582, 203 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198 (2020).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, Science 371, 1261 (2021).
  * Chen  _et al._ [2020] G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature 579, 56 (2020).
  * 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, Science 367, 900 (2020).
  * Stepanov  _et al._ [2021] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Phys. Rev. Lett. 127, 197701 (2021).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020).
  * Hofmann  _et al._ [2021] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, Fermionic monte carlo study of a realistic model of twisted bilayer graphene (2021), arXiv:2105.12112 .
  * Shavit  _et al._ [2021] G. Shavit, E. Berg, A. Stern, and Y. Oreg, Phys. Rev. Lett. 127, 247703 (2021).
  * Wagner  _et al._ [2021] G. Wagner, Y. H. Kwan, N. Bultinck, S. H. Simon, and S. A. Parameswaran, Global phase diagram of twisted bilayer graphene above $t_{c}$ (2021), arXiv:2109.09749 .
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).
  * 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, Nature 581, 47 (2020).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 572, 101 (2019).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Nature 572, 95 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Nature Phys. 15, 1174 (2019).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Huder  _et al._ [2018] L. Huder, A. Artaud, T. Le Quang, G. T. de Laissardière, A. G. M. Jansen, G. Lapertot, C. Chapelier, and V. T. Renard, Phys. Rev. Lett. 120, 156405 (2018).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Phys. Rev. B 100, 035448 (2019).
  * [25] See Supplemental Material including Refs. [40, 41, 42, 43] for details regarding strain effects, mean-field approach, local impurity calculations, and the Abrikosov-Gor’kov results.
  * Kolář  _et al._ [2022] K. Kolář, G. Shavit, C. Mora, Y. Oreg, and F. von Oppen, Anderson’s theorem for correlated insulating states in twisted bilayer graphene (2022), arXiv:2207.11281 .
  * Note [1] We note that a more realistic model of the band structure of MATBG contains terms which break particle-hole symmetry, as well as terms proportional to $\tau_{z}$. These break $\cal S$ and $\cal C$, reducing the symmetry of $h_{\rm MF}$ from symmetry class CII to AII.
  * Vozmediano  _et al._ [2010] M. Vozmediano, M. Katsnelson, and F. Guinea, Physics Rep. 496, 109 (2010).
  * Yu [1965] L. Yu, Acta Phys. Sin. 21, 75 (1965).
  * Shiba [1968] H. Shiba, Prog. Theor. Phys. 40, 435 (1968).
  * Rusinov [1969] A. I. Rusinov, Zh. Eksp. Teor. Fiz. 9, 146 (1969), [JETP Lett. 9, 85 (1969)].
  * Abrikosov and Gor’kov [1959] A. Abrikosov and L. Gor’kov, Sov. Phys. JETP 8, 1090 (1959).
  * Zittartz [1967] J. Zittartz, Phys. Rev. 164, 575 (1967).
  * Bistritzer and MacDonald [2008] R. Bistritzer and A. H. MacDonald, Phys. Rev. Lett. 101, 256406 (2008).
  * Woolf and Reif [1965] M. A. Woolf and F. Reif, Phys. Rev. 137, A557 (1965).
  * Bastin  _et al._ [1971] A. Bastin, C. Lewiner, O. Betbeder-Matibet, and P. Nozieres, J. Phys. Chem. Sol. 32, 1811 (1971).
  * Calderón and Bascones [2020] M. J. Calderón and E. Bascones, Phys. Rev. B 102, 155149 (2020).
  * Goodwin  _et al._ [2020] Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Hartree theory calculations of quasiparticle properties in twisted bilayer graphene (2020), arXiv:2004.14784 .
  * Xie and MacDonald [2021] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 127, 196401 (2021).
  * Balatsky  _et al._ [2006] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
  * Pientka  _et al._ [2013] F. Pientka, L. I. Glazman, and F. von Oppen, Phys. Rev. B 88, 155420 (2013).
  * Maki [1969] K. Maki, in  _Superconductivity_ , Vol. 2, edited by R. D. Parks (Marcel Dekker, 1969) p. 1035.
  * Hirschfeld  _et al._ [1986] P. Hirschfeld, D. Vollhardt, and P. Wölfle, Solid State Commun. 59, 111 (1986).

## Supplemental material

##  S.1 Components of the strain tensor

Consider the two graphene layers which make up the twisted bilayer. One may
apply a spatially-dependent displacement field to the carbon atoms of each
layer, $u^{\alpha}_{\rm displacement}\left(\ell,{\mathbf{r}}\right)$. Here,
$\mathbf{r}$ is the position of the displaced carbon atom, $\ell={\rm t},{\rm
b}$ corresponds to either the top or bottom layer, and $\alpha=x,y$ is the
axis along which the atom is displaced.

We then define the strain tensor

$\epsilon^{\ell}_{\alpha\beta}=\frac{\partial_{\alpha}u_{\rm
displacement}^{\beta}\left(\ell,{\mathbf{r}}\right)+\partial_{\beta}u_{\rm
displacement}^{\alpha}\left(\ell,{\mathbf{r}}\right)}{2}.$ (S1)

The homostrain and heterostrain tensor are now simply the layer-symmetric and
antisymmetric components of $\epsilon^{\ell}_{\alpha\beta}$,

$\epsilon^{\rm homo}_{\alpha\beta}=\frac{\epsilon^{\rm
t}_{\alpha\beta}+\epsilon^{\rm b}_{\alpha\beta}}{2},$ (S2) $\epsilon^{\rm
hetero}_{\alpha\beta}=\frac{\epsilon^{\rm t}_{\alpha\beta}-\epsilon^{\rm
b}_{\alpha\beta}}{2},$ (S3)

and one therefore expects both components to exist for a generic random strain
field configuration.

###  S.1.1 Uniform strain

Uniform heterostrain.—The effects of uniform heterostrain on both the band
structure of magic angle twisted bilayer graphene (MATBG) and on the
correlated insulator appearing at even fillings of electron/holes per moiré
unit-cell have been considered in Refs. [13, 22, 23, 24]. The main notable
effect of this perturbation is the distortion of the moiré superlattice,
leading to a significant increase of the flat-bands’ bandwidth. This then
leads to weakening of the correlated insulator, due to an effective weaker
interaction (as compared to the bandwidth of the electrons).

Uniform homostrain.—Uniform strain which is applied identically to both
graphene layers will shift the band structure in the two valleys in opposite
directions in momentum space [28],

${\mathbf{K}}\to{\mathbf{K}}+{\mathbf{A}},\,\,\,\,{\mathbf{K}^{\prime}}\to{\mathbf{K}^{\prime}}-{\mathbf{A}},$
(S4)

where ${\mathbf{K}},{\mathbf{K}^{\prime}}$ are the two monolayer graphene
valleys and ${\mathbf{A}}$ is the pseudo-gauge field.

Assuming electron-electron interactions which are weakly dependent on
momentum, consider a spontaneous intervalley-coherent order parameter, which
has an oscillatory spatial dependence $\Delta_{\rm
ivc}\left({\mathbf{r}}\right)\propto\exp\left(i{\mathbf{Q}}\cdot{\mathbf{r}}\right)$
in the absence of strain, with
${\mathbf{Q}}={\mathbf{K}}-{\mathbf{K}^{\prime}}$ the intervalley separation.
Then, applying uniform homostrain would have no effect on the instability
towards the intervalley order, though the “new” order parameter will now have
a slightly different spatial periodicity, $\Delta_{\rm
ivc}\left({\mathbf{r}}\right)\propto\exp\left[i\left({\mathbf{Q}+2{\mathbf{A}}}\right)\cdot{\mathbf{r}}\right]$.
We therefore do not include this type of perturbation in our considerations in
this work.

###  S.1.2 Strain disorder

Let us now discuss the effects of the disordered strain potentials within our
effective model. As a reminder, its mean-field form is

$h_{\rm
MF}\left({\mathbf{k}}\right)=u\left(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\right)+\Delta_{\rm
ivc}\sigma_{x}\tau_{x}\rho_{z},$ (S5)

where sublattice, valley, and ”mini-valley” degrees of freedom are described
by Pauli matrices $\sigma_{i}$, $\tau_{i}$, and $\rho_{i}$, respectively.

In this language, homostrain disorder is proportional to the operators
$\sigma_{x}$, $\sigma_{y}\tau_{z}$, or a linear combination of the two.
Therefore, the homostrain operator commutes with the order parameter
$\propto\sigma_{x}\tau_{x}\rho_{z}$. This fact is crucial for allowing this
disorder term to close the gap induced by the intervalley order parameter, as
we discuss throughout this work.

Heterostrain disorder however, is proportional to the operators
$\sigma_{x}\tau_{z}\rho_{z}$, $\sigma_{y}\rho_{z}$ and their combinations. One
easy way to understand why, is to recognize that the same mini-valley in
opposite valleys originates in the Dirac points corresponding to opposite
layers. Therefore, when it comes to response to a strain potential, the
operator $\tau_{z}\rho_{z}$ effectively acts as an “originating-layer” index.
Heterostrain operators are then naturally given by multiplying the homostrain
operators by $\tau_{z}\rho_{z}$, enforcing opposite pseudo-gauge fields in
opposite layers.

Importantly, the heterostrain operators anticommute with the intervalley-
coherent order parameter, so that the gap cannot be closed by weak disorder of
that kind. Their effect is thus analogous to that of non-magnetic charge
impurities in conventional spin-singlet superconductors. Therefore, we do not
consider this type of disorder in this work.

##  S.2 Mean-field results in the pristine case

Let us briefly discuss the model

$H=H_{0}+H_{\rm int},$ (S6)
$H_{0}=\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}\left[uk_{x}\sigma_{x}\tau_{z}+uk_{y}\sigma_{y}\right]\Psi_{\mathbf{k}},$
(S7) $H_{\rm
int}=\frac{U}{2\Omega}\sum_{{\mathbf{k}},{\mathbf{k}^{\prime}},{\mathbf{q}}}\Psi^{\dagger}_{\mathbf{k+q}}\Psi_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k^{\prime}-q}}\Psi_{\mathbf{k^{\prime}}}.$
(S8)

$\Psi^{\dagger}_{\mathbf{k}}$ is an 8-spinor of fermionic annihilation
operators at momentum $\mathbf{k}$ relative to their respective origin in
momentum space, and with sublattice, valley, and ”mini-valley” degrees of
freedom, described by Pauli matrices $\sigma_{i}$, $\tau_{i}$, and $\rho_{i}$,
respectively.

We note that the basis used in this model is quite different from the commonly
used “sublattice basis” (see, e.g., Refs. [13, 15]), since we effectively
separate the electronic species in each valley into two species near the two
Dirac cones in each valley. In this case, the momentum $\mathbf{k}$ is
measured relative to the appropriate Dirac point labeled by both the quantum
numbers $\tau_{z}$ and $\rho_{z}$. As we show momentarily, this also
translates to a different form of the intervalley-coherent order parameter.

Our model is meant to provide an effective description of the physics near
$\nu=\pm 2$ filling in MATBG. It is widely believed that in these regimes, the
flavor (spin and valley) symmetry is spontaneously broken due to strong
electron-electron interactions [7, 8]. Consequently, two flavors are either
completely full or completely empty, whereas the remaining two flavors have
their chemical potential at the Dirac point. Eq. (S7) describes precisely the
two remaining active flavors in this scenario.

We are interested in the order parameter

$\Delta_{{\rm
ivc}}=-\frac{U}{2\Omega}\sum_{\mathbf{k}}\left\langle\Psi^{\dagger}\left(\mathbf{k}\right)\sigma_{x}\tau_{x}\rho_{z}\Psi\left(\mathbf{k}\right)\right\rangle$
(S9)

for Kramers intervalley-coherent (K-IVC) order. We note that $\Delta_{{\rm
ivc}}$ spontaneously breaks the valley-$U\left(1\right)$ symmetry in our
model, manifested in $\tau_{x}\to\tau_{x}e^{i\alpha_{\rm valley}\tau_{z}}$.

Performing a mean-field decomposition on the Hamiltonian, we find

$H_{{\rm
MF}}=\sum_{\mathbf{k}}\Psi^{\dagger}\left(\mathbf{k}\right)\left[uk_{x}\sigma_{x}\tau_{z}+uk_{y}\sigma_{y}+\Delta_{{\rm
ivc}}\sigma_{x}\tau_{x}\rho_{z}\right]\Psi\left(\mathbf{k}\right)+2\frac{\Omega}{U}\Delta_{{\rm
ivc}}^{2}.$ (S10)

Importantly, the order parameter $\Delta_{\rm ivc}$ is the only possible K-IVC
order parameter (up to a rotation in the $\tau_{x}$–$\tau_{y}$ plane), i.e.,
it fulfills the three conditions: (i) it is valley off-diagonal, (ii) it opens
a gap in the spectrum (anti-commutes with $H_{0}$), and (iii) it breaks the
time-reversal symmetry ${\cal T}=\tau_{x}\rho_{x}{\cal K}$, while preserving
the Kramers one ${\cal T}^{\prime}=\tau_{y}\rho_{x}{\cal K}$ ($\cal K$ is the
complex conjugation operator).

Defining $\epsilon_{\mathbf{k}}=u\left|{\mathbf{k}}\right|$,
$E_{\mathbf{k}}=\sqrt{\epsilon_{\mathbf{k}}^{2}+\Delta_{\rm ivc}^{2}}$, we
find the mean-field thermodynamic potential $\Phi$ at temperature $T$,

$\Phi=2\frac{\Omega}{U}\Delta_{{\rm
ivc}}^{2}-8T\sum_{\mathbf{k}}\log\left[2\cosh\left(\frac{E_{\mathbf{k}}}{2T}\right)\right].$
(S11)

The self-consistent gap equation follows from the condition
$\frac{\partial\Phi}{\partial\Delta_{{\rm ivc}}}=0$, which leads to

$\frac{1}{U}=\int d\epsilon{\cal
N}\left(\epsilon\right)\frac{\tanh\left(\frac{\sqrt{\epsilon^{2}+\Delta_{{\rm
ivc}}^{2}}}{2T}\right)}{\sqrt{\epsilon^{2}+\Delta_{{\rm ivc}}^{2}}},$ (S12)

where ${\cal N}\left(\epsilon\right)$ is the density of states per unit cell
and flavor (spin, valley, minivalley) of $H_{0}$ at energy $\epsilon$.

We may extract some analytical results by assuming a strictly linear density
of states,

${\cal
N}\left(\epsilon\right)=\begin{cases}\frac{4\left|\epsilon\right|}{W^{2}}&\left|\epsilon\right|<W/2,\\\
0&{\rm otherwise},\end{cases}$ (S13)

with $W$ the total bandwidth. Setting $T=0$, we find

$\Delta_{00}=\Delta_{\rm
ivc}\left(T=0\right)=U\left[1-\left(\frac{U_{c}}{U}\right)^{2}\right].$ (S14)

Here, we have defined the critical interaction strength $U_{c}=W/4$, below
which there is no gap opening. The existence of a critical coupling contrasts
with the familiar behavior of the standard Bardeen-Cooper-Schrieffer (BCS) gap
equation, which has the same form as Eq. (S12). The difference originates from
the density of states ${\cal N}\left(\epsilon\right)$. In the present case,
the density of states vanishes at the Dirac point ($\epsilon=0$), which
eliminates the Cooper instability for infinitesimally weak interaction
strengths.

We may also extract the critical temperature $T_{c0}$ by setting $\Delta_{\rm
ivc}=0$ in the gap equation. This leads to the transcendental equation

$\frac{U_{c}}{U}=\frac{T_{c0}}{U_{c}}\log\left[\cosh\left(\frac{U_{c}}{T_{c0}}\right)\right],$
(S15)

which does not have a $T_{c0}>0$ solution for $U<U_{c}$, as expected. Close to
the transition, when $T_{c0}\ll U_{c}$, we may approximate
$T_{c0}\approx\frac{1}{\log 2}\frac{U_{c}}{U}\left(U-U_{c}\right)$.

##  S.3 Point-like perturbation

In this section, we analyze the effect of a localized strain “impurity” on the
mean-field K-IVC state. We assume that the strain does not suppress the
formation of intervalley coherent order for the sake of this discussion. Thus,
we investigate the Hamiltonian

$\displaystyle H_{\rm imp}$
$\displaystyle=\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}\left[uk_{x}\sigma_{x}\tau_{z}+uk_{y}\sigma_{y}+\Delta_{\rm
ivc}\sigma_{x}\tau_{x}\rho_{z}\right]\Psi_{\mathbf{k}}$
$\displaystyle+\frac{1}{\Omega}\sum_{{\mathbf{k}},{\mathbf{k}^{\prime}}}\Psi^{\dagger}_{\mathbf{k}}\left[v_{F}{\cal
A}_{x}\sigma_{x}+v_{F}{\cal
A}_{y}\sigma_{y}\tau_{z}\right]\Psi_{\mathbf{k}^{\prime}},$ (S16)

where ${\cal A}_{x,y}$ are constants, reflecting the delta-potential like
spatial dependence of ${\mathbfcal A}\left({\mathbf{r}}\right)$, to be
determined by the strain tensor, see Eq. (S40). Notice that the velocity in
this impurity strain potential is the bare Fermi velocity of graphene $v_{F}$,
and not the renormalized one $u$. This is because the deformed potential is
not smooth on the moiré scale, as we have assumed it is pointlike by
definition. We note that much like in the case of the strain perturbation we
consider in the main text, $H_{\rm imp}$ preserves both $\cal T$ and ${\cal
T}^{\prime}$ time-reversal symmetries.

The Green’s function for the ${\cal A}_{x,y}=0$ system can be written as

$G^{0}\left({\mathbf{k}},\omega\right)=-\frac{\omega+uk_{x}\sigma_{x}\tau_{z}+uk_{y}\sigma_{y}+\Delta_{{\rm
ivc}}\sigma_{x}\tau_{x}\rho_{z}}{\epsilon_{\mathbf{k}}^{2}+\Delta_{{\rm
ivc}}^{2}-\omega^{2}},$ (S17)

with $\epsilon_{\mathbf{k}}=u\left|{\mathbf{k}}\right|$. We also calculate

$g^{0}\left(\omega\right)\equiv\frac{1}{\Omega}\sum_{\mathbf{k}}G^{0}\left({\mathbf{k}},\omega\right)=-\left(\omega+\Delta_{{\rm
ivc}}\sigma_{x}\tau_{x}\rho_{z}\right)C,$ (S18) $C=\int d\epsilon\frac{{\cal
N}\left(\epsilon\right)}{\epsilon^{2}+\Delta_{{\rm ivc}}^{2}-\omega^{2}},$
(S19)

which we will use momentarily.

###  S.3.1 $T$-matrix formalism

Our goal is to find the in-gap bound state spectrum due to the strain
perturbation. In order to do so, we shall employ the $T$-matrix formalism [43,
40]. Within this formalism, in the case of a point-like perturbation at the
origin, the real-space Green’s function may be written as

$G\left({\mathbf{r}},{\mathbf{r}^{\prime}};\omega\right)=G^{0}\left({\mathbf{r}},{\mathbf{r}^{\prime}};\omega\right)+G^{0}\left({\mathbf{r}},{\mathbf{0}};\omega\right)T\left(\omega\right)G^{0}\left({\mathbf{0}},{\mathbf{r}^{\prime}};\omega\right),$
(S20)

with $G^{0}\left({\mathbf{r}},{\mathbf{r}^{\prime}};\omega\right)$ the
unperturbed Green’s function and the $T$-matrix

$T\left(\omega\right)=\frac{U_{0}}{1-U_{0}g^{0}\left(\omega\right)}.$ (S21)

In our case $U_{0}$ is given by the matrix $U_{0}=v_{F}{\cal
A}_{x}\sigma_{x}+v_{F}{\cal A}_{y}\sigma_{y}\tau_{z}$. Notice
$\left[U_{0},\Delta_{\rm ivc}\sigma_{x}\tau_{x}\rho_{z}\right]=0$.

The quasiparticle spectrum may be recovered by finding the poles of $G$. Thus,
to calculate the perturbation-induced in-gap bound state spectrum it is
sufficient to find the poles of the $T$-matrix. Thus, the in-gap states have
energy $\omega$ such that the matrix

$Q\equiv 1-U_{0}g^{0}\left(\omega\right),$ (S22)

has a zero eigenvalue. The eigenvalues of $Q$ are

$\lambda_{Q}=1\pm v_{F}{\cal A}C\left(\Delta_{\rm ivc}\pm\omega\right),$ (S23)

with ${\cal A}=\sqrt{{\cal A}_{x}^{2}+{\cal A}_{y}^{2}}$. Each of these four
eigenvalues is doubly degenerate. Notice the particle-hole symmetry. If
$\omega^{*}$ makes one eigenvalue go to zero, then $-\omega^{*}$ makes another
eigenvalue go to zero as well.

In order to make further progress, we must address $C$, which is not a
constant, but rather a function of $\omega$ and $\Delta_{\rm ivc}$.

###  S.3.2 Linear-in-$\epsilon$ density of states

Using the linear density of states we employed in Eq. (S13), we find

$C=\frac{1}{\left(W/2\right)^{2}}\log\left[1+\frac{\left(W/2\right)^{2}}{\Delta_{{\rm
ivc}}^{2}-\omega^{2}}\right].$ (S24)

It is instructive to find the impurity strength at which the gap closes ,
i.e., $\omega=0$ leads to a zero $\lambda_{Q}$. This is given by ${\cal
A}_{0}$,

$\frac{v_{F}{\cal A}_{0}}{W/2}=\frac{W/2}{\Delta_{{\rm
ivc}}}\frac{1}{\log\left[1+\left(\frac{W/2}{\Delta_{{\rm
ivc}}}\right)^{2}\right]}.$ (S25)

Additionally, it is interesting to find the behavior of the in-gap state at
very small strain, $v_{F}{\cal A}\ll W$. In this limit, we find

$\omega\approx\pm\Delta_{{\rm
ivc}}\left[1-\left(\frac{W/2}{2\Delta}\right)^{2}\exp\left(-\frac{W/2}{2\Delta}\frac{W/2}{v_{F}{\cal
A}}\right)\right],$ (S26)

i.e., the bound state is exponentially close to the gap edge. This behavior
can indeed be seen in Fig. S1 at small $\cal A$.

Figure S1:  Bound state spectrum of the strain impurity problem. The two
branches of the spectrum appear in blue and red, showing the particle-hole
symmetry. Each of the branches is doubly degenerate, reflecting the Kramers-
like time-reversal symmetry ${\cal T}^{\prime}$ of the Hamiltonian, even in
the presence of the strain perturbation. For this plot we used $\Delta_{\rm
ivc}=W/5$.

###  S.3.3 Constant density of states

In order to shed some light on the results of the last section, as well as to
further facilitate the comparison to a magnetic impurity embedded in a singlet
superconductor, let us repeat the calculation with a constant density of
states, ${\cal N}\left(\epsilon\right)={\cal N}_{0}$. In this case,

$C=\frac{\pi{\cal N}_{0}}{\sqrt{\Delta_{{\rm ivc}}^{2}-\omega^{2}}}.$ (S27)

The bound state spectrum is thus given by

$\omega=\pm\Delta_{\rm ivc}\frac{1-s^{2}}{1+s^{2}},$ (S28)

with the dimensionless impurity strength $s=\pi{\cal N}_{0}v_{F}{\cal A}$. It
is worth mentioning that the result of Eq. (S28) is identical to the bound
state energies of in-gap states induced in a singlet superconductor in the
presence of a magnetic impurities, see, e.g., Ref. [41]. Thus, although the
perturbation breaks no symmetries of the mean-field Hamiltonian, due to its
matrix form which commutes with the order parameter – an in-gap state (which
may eventually reach all the way to zero energy) becomes possible.

##  S.4 Details of the Abrikosov-Gorkov calculations

Calculations in this section closely resemble the techniques used to treat
magnetic impurities in singlet superconductors [32, 42] and excitonic
condensation in semiconductors [33] and graphene [34]. We combine self-
consistent mean-field calculations for the order parameter $\Delta_{\rm ivc}$,
together with a self-consistent Born approximation, which allows us to handle
the random strain fluctuations.

The key player in this treatment of the problem is the Green’s function, which
satisfies (within the mean-field approximation)

$G\left({\mathbf{k}},\omega\right)=\left(i\omega-h_{\rm MF}-\Sigma_{\rm
SCBA}\right)^{-1},$ (S29)

with the mean-field Hamiltonian

$h_{\rm MF}=u\left(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\right)+\Delta_{\rm
ivc}\sigma_{x}\tau_{x}\rho_{z}.$ (S30)

Within the self-consistent Born approximation, the self-energy is given by

$\Sigma_{\rm
SCBA}\left({\mathbf{k}},\omega\right)=\left\langle\sum_{|\mathbf{p}|<\Lambda_{c}}{\cal
U}_{{\mathbf{k}}-{\mathbf{p}}}G\left({\mathbf{p}},\omega\right){\cal
U}_{{\mathbf{p}}-{\mathbf{k}}}\right\rangle_{\rm disorder}$ (S31)

where the matrix $\cal U$ represents the random strain perturbation in
momentum space, $H_{\rm
str}=\sum_{\mathbf{kq}}\Psi^{\dagger}_{\mathbf{k+q}}{\cal
U}_{\mathbf{q}}\Psi_{\mathbf{k}}$, and $\left\langle...\right\rangle_{\rm
disorder}$ stands for averaging over disorder configurations. Notice that Eq.
(S31) contains within it the full Green’s function, such that Eq. (S29) and
Eq. (S31) must be solved together self-consistently.

The strain operator has the form

${\cal
U}_{\mathbf{q}}=u\left(A_{x,{\mathbf{q}}}\sigma_{x}+A_{y,{\mathbf{q}}}\sigma_{y}\tau_{z}\right).$
(S32)

Assuming the disorder averaging eliminates any anisotropy introduced by a
specific disorder configuration, we may assume $\left\langle
A_{x,{\mathbf{q}}}^{*}A_{y,{\mathbf{q}}}\right\rangle_{\rm
disorder}=\left\langle
A_{y,{\mathbf{q}}}^{*}A_{x,{\mathbf{q}}}\right\rangle_{\rm disorder}=0$, and
$\left\langle
A_{x,{\mathbf{q}}}^{*}A_{x,{\mathbf{q}}}-A_{y,{\mathbf{q}}}^{*}A_{y,{\mathbf{q}}}\right\rangle_{\rm
disorder}=0$. Using the fact that the real-space strain perturbation is real,
we have $A_{i,-{\mathbf{q}}}=A^{*}_{i,{\mathbf{q}}}$. Employing the ansatz,

$G_{\rm
SCBA}\left({\mathbf{k}},\omega\right)=\left(i\tilde{\omega}-uk_{x}\sigma_{x}\tau_{z}-uk_{y}\sigma_{y}-\tilde{\Delta}\sigma_{x}\tau_{x}\rho_{z}\right)^{-1},$
(S33)

the self-energy takes the form

$\Sigma_{\rm
SCBA}\left({\mathbf{k}},\omega\right)-=\left\langle\sum_{\mathbf{p}}\left|u\vec{A}_{\mathbf{p-q}}\right|^{2}\frac{i\tilde{\omega}+\tilde{\Delta}\sigma_{x}\tau_{x}\rho_{z}}{\epsilon_{\mathbf{p}}^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}}\right\rangle_{\rm
disorder}.$ (S34)

For simplicity, we now make the standard assumption that the disorder
potential $\vec{A}_{\mathbf{p-q}}$ depends mostly on the relative angle
between the initial and final momenta, i.e.,
$\vec{A}_{\mathbf{p-q}}\approx\vec{A}_{\theta}$, where the angle is defined by
$\tan\theta=\frac{p_{y}-q_{y}}{p_{x}-q_{x}}$. This assumption is due to the
singularity of the Green’s function at the Fermi level, such that the main
contribution in the sum appearing in Eq. (S34) is by momenta with
$\left|{\mathbf{p}}\right|\approx p_{F}$, with $p_{F}$ the Fermi momentum.
With this approximation, we find the self-consistent equations for
$\tilde{\omega}$ and $\tilde{\Delta}$,

$\tilde{\omega}=\omega+\tilde{\omega}\Gamma\frac{W}{2}\int
d\epsilon\frac{{\cal
N}\left(\epsilon\right)}{\epsilon^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}},$
(S35) $\tilde{\Delta}=\Delta-\tilde{\Delta}\Gamma\frac{W}{2}\int
d\epsilon\frac{{\cal
N}\left(\epsilon\right)}{\epsilon^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}},$
(S36)

where ${\cal N}\left(\epsilon\right)$ is the density of states per moiré unit-
cell $\Omega_{u.c.}$, $W$ is the bandwidth, and the scattering rate $\Gamma$
is

$\Gamma=\frac{2\Omega}{W\Omega_{u.c.}}\left\langle\int d\theta
u^{2}\vec{A}_{\theta}^{\dagger}\cdot\vec{A}_{\theta}\right\rangle_{\rm
disorder}.$ (S37)

Finally, to complete the picture, we must include the self-consistent
calculation of the order parameter itself, with the gap equation

$\displaystyle\Delta_{{\rm ivc}}$
$\displaystyle=-2\frac{U}{\beta\Omega}\sum_{\omega\mathbf{k}}{\rm
Tr}\left\\{\sigma_{x}\tau_{x}\rho_{z}G_{\rm
SCBA}\left(\mathbf{k},\omega\right)\right\\}$
$\displaystyle=2U\frac{1}{\beta}\sum_{\omega}\int d\epsilon\frac{{\cal
N}\left(\epsilon\right)}{\epsilon^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}}\tilde{\Delta}.$
(S38)

Equations (S35),(S36), and (S38) form a set of closed equations that can be
solved together to recover $\Delta_{\rm ivc}$ and $G_{\rm SCBA}$ in the
presence of strain fluctuations.

Lastly, we may also calculate the tunneling density of states, given by

$\displaystyle\tilde{\cal N}\left(\epsilon\right)$
$\displaystyle=\frac{1}{\pi}{\rm Im}\frac{1}{\Omega}\sum_{\mathbf{k}}{\rm
Tr}G_{\rm SCBA}\left({\mathbf{k}},\omega\to i\epsilon\right)$
$\displaystyle=\frac{1}{\pi}{\rm Im}\int d\epsilon{\cal
N}\left(\epsilon\right)\frac{i\tilde{\omega}}{\epsilon^{2}+\tilde{\Delta}^{2}+\tilde{\omega}^{2}}|_{\omega\to
i\epsilon}.$ (S39)

When $\Delta_{\rm ivc}$ is finite, one expects to find a gap in $\tilde{\cal
N}\left(\epsilon\right)$ near zero energy. We denote the value of this gap as
$\omega_{g}$.

###  S.4.1 Estimating the disorder scale $\Gamma$

The disorder parameter $\Gamma$ determines completely the properties of the
disorder-averaged Green’s function within our SCBA treatment. Given the
homostrain tensor $\epsilon^{\rm homo}_{\alpha\beta}$ introduced in Eq. (S2),
the effective pseudo-gauge field can be expressed as [17]

${\mathbfcal A}=\frac{\sqrt{3}\beta}{a}\left(\frac{\epsilon^{\rm
homo}_{xx}-\epsilon^{\rm homo}_{yy}}{2},-\epsilon^{\rm homo}_{xy}\right),$
(S40)

where $a=0.246\,{\rm nm}$ is the monolayer graphene lattice constant, and
$\beta\approx 3.14$ is a dimensionless parameter of graphene, which quantifies
the sensitivity of nearest-neighbor hopping strength to a change in their
relative distance.

As an example, let us consider a scenario where root-mean-square strain is
$\epsilon_{\rm RMS}$, approximate $u\approx 0.1v_{F}$ [2], and use $\hbar
v_{F}/a\approx 2.68eV$, such that

$u\left|{\cal A}_{\rm RMS}\right|\approx 14.5\,{\rm meV}\times\epsilon_{\rm
RMS}\left[\%\right].$ (S41)

In order to continue, we examine the specific example of Gaussian correlated
disorder, i.e., assume the following correlations,

$\left\langle{\cal A}_{i}\left(\mathbf{r}\right){\cal
A}_{j}\left(\mathbf{r^{\prime}}\right)\right\rangle_{\rm
disorder}=\delta_{ij}{\cal A}_{\rm
RMS}^{2}\exp\left({-\frac{\left|\mathbf{r}-\mathbf{r^{\prime}}\right|^{2}}{2\xi_{\rm
dis.}^{2}}}\right),$ (S42)

where $\xi_{\rm dis.}$ is the disorder correlation length. Following a
straightforward Fourier transform, we find

$\left\langle{A}_{i}\left(\mathbf{q}\right){A}_{j}\left(\mathbf{-q}\right)\right\rangle_{\rm
disorder}=\delta_{ij}{\cal A}_{\rm RMS}^{2}\frac{1}{\Omega}\xi_{\rm
dis.}^{2}\exp\left(-\frac{\left|\mathbf{q}\right|^{2}\xi_{\rm
dis.}^{2}}{2}\right).$ (S43)

Since we began with the expression in Eq. (S34), our main focus is on momenta
near the Fermi level. We thus replace $\mathbf{q}\approx k_{F}=0$ in our
estimates below. Estimating the MATBG bandwidth as $W\approx 20\,{\rm meV}$,
and using Eq. (S37) we find

$\Gamma/W\approx\frac{4\xi_{{\rm dis.}}^{2}}{A_{{\rm u.c.}}}\left(\frac{u{\cal
A}_{{\rm RMS}}}{W}\right)^{2}\approx 2.1\times\left(\epsilon_{\rm
RMS}\left[\%\right]\times\xi_{\rm dis.}\left[a_{\rm m}\right]\right)^{2},$
(S44)

where $a_{\rm m}\approx\sqrt{A_{\rm u.c.}}$ is the moiré lattice constant.

For example, for $\epsilon_{\rm RMS}=0.1\%$ and $\xi_{\rm dis.}=3a_{\rm m}$ we
find $\Gamma/W\approx 0.19$. We can thus conclude that for appreciable yet
reasonable homostrain disorder, $\Gamma$ may indeed be large enough to
diminish both $\Delta_{\rm ivc}$ and $\omega_{g}$ significantly, and possibly
completely close the gap in the TDOS.

###  S.4.2 Critical disorder strength

If the energy scale associated with the random strain fluctuations $\Gamma$ is
sufficiently large, the K-IVC order may be completely washed away, even at
zero temperature. Using the linear density of states approximation, we can
recover this critical $\Gamma_{c}$, by setting $\Delta_{\rm
ivc},\tilde{\Delta}\to 0$ and $T=0$. We note that it is important not to set
$\Delta_{\rm ivc}=\tilde{\Delta}$ exactly, although both quantities are taken
to zero. This is because in the gap equation we must divide both sides by
$\Delta_{\rm ivc}$ Under these assumptions,

$\tilde{\omega}=\omega+\frac{\Gamma_{c}}{W/2}\log\left(1+\frac{\left(W/2\right)^{2}}{\tilde{\omega}^{2}}\right)\tilde{\omega},$
(S45)
$\tilde{\Delta}=\frac{\Delta}{1+\frac{\Gamma_{c}}{W/2}\log\left(1+\frac{\left(W/2\right)^{2}}{\tilde{\omega}^{2}}\right)},$
(S46)

from which we also calculate

$d\omega=d\tilde{\omega}\left[1-\frac{\Gamma_{c}}{W/2}\log\left(1+\frac{\left(W/2\right)^{2}}{\tilde{\omega}^{2}}\right)+2\Gamma_{c}\frac{W/2}{\tilde{\omega}^{2}+\left(W/2\right)^{2}}\right].$
(S47)

Plugging this relation into the zero temperature gap equation and changing
variables, we find

$\frac{U_{c}}{U}=\int_{-\infty}^{\infty}\frac{dx}{2\pi}\frac{1-g\log\left(1+\frac{1}{x^{2}}\right)+2g\frac{1}{x^{2}+1}}{1+g\log\left(1+\frac{1}{x^{2}}\right)}\log\left(1+\frac{1}{x^{2}}\right),$
(S48)

where $g=2\Gamma_{c}/W$. The self-consistent Eq. (S48) relates the interaction
strength $U$ to the critical disorder strength via a dimensionless integral.
For relatively weak disorder, $g\ll 1$, we may expand Eq. (S48) to linear
order in $g$, to find the relation

$\Gamma=\frac{U_{c}}{\log 8}\left(1-\frac{U_{c}}{U}\right),$ (S49)

which is valid close for $U\approx U_{c}$, such that $g\ll 1$.

###  S.4.3 Gapless K-IVC

Unlike in the translation-symmetry invariant case, the gap in the tunneling
density of states, $2\omega_{g}$, does not necessarily equal the self-
consistent order parameter $\Delta_{\rm ivc}$. More concretely, there exists a
regime in which $\Delta_{\rm ivc}$ is finite, yet $\omega_{g}=0$ and the
tunneling density of states is gapless.

Let us find the disorder strength $\Gamma_{g}$ for which this gap closes. To
do so, we calculate the tunneling density of states at zero energy
$\tilde{\cal N}\left(\epsilon=0\right)$. Plugging $\omega\to 0$ in the
expressions for $\tilde{\omega}$ and $\tilde{\Delta}$, we find

$\tilde{\omega}^{2}=\frac{\left(W/2\right)^{2}}{\exp{\frac{W/2}{\Gamma}}-1}-\left(\Delta_{{\rm
ivc}}/2\right)^{2}.$ (S50)

The tunneling density of states is then given by

${\cal\tilde{N}}\left(\epsilon=0\right)=\frac{1}{2\pi}\frac{1}{\Gamma
W}\times{\rm Re}\left\\{\tilde{\omega}_{\omega\to 0}\right\\}.$ (S51)

Thus, the gap closes, i.e., there is a finite value of $\tilde{\cal N}$ at
zero energy, for

$\Gamma\geq\Gamma_{g}=\frac{W/2}{\log\left(1+\frac{W^{2}}{\Delta_{{\rm
ivc}}^{2}}\right)}.$ (S52)

Notice Eq. (S52) is not a closed-form expression for $\Gamma_{g}$, since in
the right hand side $\Delta_{\rm ivc}$ depends on $\Gamma$ itself.

Figure S2:  (a) The order parameter $\Delta_{\rm ivc}$ as a function of
interaction strength and the disorder energy scale $\Gamma$. At a given
interaction strength $U$, a critical $\Gamma$ exists at which the order
parameter vanishes. (b) The ratio of the gap in the TDOS, $\omega_{g}$, to
$\Delta_{\rm ivc}$. As $\Gamma$ increases, $\omega_{g}$ deviates from
$\Delta_{\rm ivc}$ appreciably. Moreover, there exists a gapless K-IVC region
bounded by the vanishing of the gap (black dashed line) and of the order
parameter (blue dashed line). Calculations were done at zero temperature,
using the linear DOS ${\cal N}_{\rm lin.}$.

In the main text, we have shown that the gapless K-IVC regime exists over a
range of temperatures and disorder strengths. For completeness, we also show
how the gapless K-IVC at zero temperature evolves as a function of the
repulsive interaction strength $U$. This is plotted in Fig. S2.

###  S.4.4 Conductivity calculations

The primary method by which the K-IVC gap is extracted in experiments is by
derivation of the activation energy from transport measurements. Namely, the
longitudinal conductivity $\sigma_{xx}$ is measured as a function of
temperature at the appropriate filling corresponding to the insulating phase,
and fitted to $\propto\exp\left(-E_{\rm act.}/T\right)$, where the activation
energy $E_{\rm act.}$ presumably coincides with the K-IVC order parameter
$\Delta_{\rm ivc}$.

In order to put our results thus far in a similar context, we employ the Kubo-
Bastin formula for the dc conductivity [36],

$\sigma_{xx}\propto\int d\epsilon\left(-\frac{df}{d\epsilon}\right){\cal
S}\left(\epsilon\right),$ (S53)

where $f$ is the Fermi function, and the kernel of the integral is

${\cal S}\left(\epsilon\right)=\frac{1}{\Omega}\sum_{\mathbf{k}}{\rm
Tr}\left\\{\left\langle j_{x}{\rm Im}G_{\rm SCBA}\left({\mathbf{k}},\omega\to
i\epsilon\right)j_{x}{\rm Im}G_{\rm SCBA}\left({\mathbf{k}},\omega\to
i\epsilon\right)\right\rangle_{\rm disorder}\right\\},$ (S54)

where the appropriate current operator in our case is
$j_{x}=\partial_{k_{x}}h_{\rm MF}=u\sigma_{x}\tau_{z}$.

The full calculation of Eq. (S54), including the so-called vertex corrections
to the current operator, are beyond the scope of this work. However,
qualitatively most of the disorder effects stemming from the Abrikosov-Gor’kov
treatment we perform may be captured by approximating

${\cal S}\left(\epsilon\right)\approx\frac{1}{\Omega}\sum_{\mathbf{k}}{\rm
Tr}\left\\{j_{x}{\rm Im}G_{\rm SCBA}\left({\mathbf{k}},\omega\to
i\epsilon\right)j_{x}{\rm Im}G_{\rm SCBA}\left({\mathbf{k}},\omega\to
i\epsilon\right)\right\\},$ (S55)

where we emphasize that the Green’s function $G_{\rm SCBA}$ is calculated
self-consistently accounting for the strain disorder.

Plugging in the Green’s function [Eq. (S33)], we find by a straightforward
calculation that the conductivity kernel ${\cal S}\left(\epsilon\right)$ has
the following form,

$\displaystyle{\cal S}\left(\epsilon\right)$
$\displaystyle\approx\frac{1}{\Omega}\sum_{\mathbf{k}}\frac{{\rm
Re}\left(\tilde{\omega}\right)^{2}-{\rm
Im}\left(\tilde{\Delta}\right)^{2}}{\left(\epsilon_{\mathbf{k}}^{2}+\left|\tilde{\omega}\right|^{2}+\left|\tilde{\Delta}\right|^{2}\right)^{2}-4\left[{\rm
Im}\left(\tilde{\omega}\tilde{\Delta}\right)\right]^{2}-4\left[{\rm
Im}\left(\tilde{\omega}\right)\right]^{2}\epsilon_{\mathbf{k}}^{2}}\times$
$\displaystyle\times\left\\{1+\frac{4\left[{\rm
Im}\left(\tilde{\omega}\right)\right]^{2}\epsilon_{\mathbf{k}}^{2}}{\left(\epsilon_{\mathbf{k}}^{2}+\left|\tilde{\omega}\right|^{2}+\left|\tilde{\Delta}\right|^{2}\right)^{2}-4\left[{\rm
Im}\left(\tilde{\omega}\tilde{\Delta}\right)\right]^{2}-4\left[{\rm
Im}\left(\tilde{\omega}\right)\right]^{2}\epsilon_{\mathbf{k}}^{2}}\right\\}.$
(S56)

Notice that $\tilde{\omega}$ and $\tilde{\Delta}$ are generally both complex
numbers and functions of energy $\epsilon$, through the analytic continuation
$\omega\to i\epsilon$, preceding the iterative solution of Eqs. (S35)–(S36).

Importantly we find that the functions ${\rm Re}\tilde{\omega}$ and ${\rm
Im}\tilde{\Delta}$ both vanish for $\left|\epsilon\right|<\omega_{g}$, i.e.,
the conductivity kernel ${\cal S}\left(\epsilon\right)$ has a gap of the size
$2\omega_{g}$. This fact, in conjunction with the relation in Eq. (S53),
explains our observation given in the main text, that the activation energy
actually coincides with $\omega_{g}$. To see why, let us approximate ${\cal
S}\left(\epsilon\right)\approx{\cal
S}_{0}\Theta\left(\left|\epsilon\right|-\omega_{g}\right)$ ($\Theta$ is the
Heaviside function), such that

$\sigma_{xx}\propto{\cal
S}_{0}\int_{\omega_{g}}^{\infty}d\epsilon\left(-\frac{df}{d\epsilon}\right)=\int_{\omega_{g}}^{\infty}d\epsilon\frac{{\cal
S}_{0}}{4T\cosh^{2}\frac{\epsilon}{2T}}=\frac{{\cal
S}_{0}}{2}\left(1-\tanh\frac{\omega_{g}}{2T}\right)=\frac{{\cal
S}_{0}}{1+e^{\frac{\omega_{g}}{T}}},$ (S57)

which for sufficiently low temperatures $T\ll\omega_{g}$, leads to the
$\propto\exp\left(-\omega_{g}/T\right)$ activated behavior.