# Supplemental Material: Critical Filaments and Superconductivity in
Quasiperiodic Twisted Bilayer Graphene

Xinghai Zhang  Department of Physics and Astronomy, Rice University, Houston,
Texas 77005, USA  Justin H. Wilson  Department of Physics and Astronomy,
Louisiana State University, Baton Rouge, LA 70803, USA  Center for Computation
and Technology, Louisiana State University, Baton Rouge, LA 70803, USA
Matthew S. Foster  Department of Physics and Astronomy, Rice University,
Houston, Texas 77005, USA  Rice Center for Quantum Materials, Rice University,
Houston, Texas 77005, USA

###### Contents

  1. I Numerical method for DOS, wave functions, and superconductivity
    1. I.1 Basics on Chebyshev polynomials and the kernel polynomial method
    2. I.2 Density of states
    3. I.3 Chebyshev filter methods for wave functions
      1. I.3.1 Iterative filter method
      2. I.3.2 Filter and shake method
    4. I.4 Self-consistent BCS theory based on KPM
      1. I.4.1 Self-consistent BdG theory for a lattice model
      2. I.4.2 KPM Implementation of the self-consistent BdG theory
      3. I.4.3 Continuum model with local interaction
    5. I.5 Calculation of superfluid stiffness with KPM
      1. I.5.1 Superfluid stiffness of lattice model
      2. I.5.2 KPM implementation
      3. I.5.3 Superfluid stiffness for Dirac superconductor
  2. II More Details on DOS and critical states
    1. II.1 DOS and multifractal properties in the lakes between filaments
    2. II.2 The third multifractal dimension $\tau_{3}$
  3. III Silver Ratio QBM Model
  4. IV More details on the superconductivity
    1. IV.1 Numerical results
    2. IV.2 Pairing of exact eigenstates approximation

##  I Numerical method for DOS, wave functions, and superconductivity

The two kinds of Chebyshev orthogonal polynomials are widely used to expand
functions defined in the interval $[-1,1]$. Physical quantities of a condensed
matter system are usually expressed as functions of the bounded Hamiltonian,
and are in principle well-suited to be expanded in terms of Chebyshev
polynomials. However, the functions that determine many interesting physical
properties such as the spectral and Green’s functions are not smooth or have
poles. The truncated series expansion suffers from the so-called Gibbs
phenomenon near discontinuities, where the expansion shows oscillation
behavior. The Gibbs oscillation of a finite-series expansion cannot be
resolved by simply increasing the degree of the expansion. As a resolution,
various kernels are introduced to suppress the oscillation while keeping the
series expansion positive and normalized. This is the idea of the so-called
kernel polynomial method (KPM) [1]. Here we give a brief review of the KPM
method, and related methods for wave functions, self-consistent mean-field
theory, and superfluid stiffness.

###  I.1 Basics on Chebyshev polynomials and the kernel polynomial method

The Chebyshev polynomial of the first kind is most utilized in KPM, and its
explicit expression is given by

$T_{n}\left(x\right)=\cos\left(n\arccos x\right)\,.$ (S1)

The polynomials can be obtained recursively with

$\displaystyle T_{0}\left(x\right)$ $\displaystyle=1\,,$ (S2a) $\displaystyle
T_{1}\left(x\right)$ $\displaystyle=x\,,$ (S2b) $\displaystyle
T_{n+1}\left(x\right)$
$\displaystyle=2xT_{n}\left(x\right)-T_{n-1}\left(x\right)\,.$ (S2c)

A function $f\left(x\right)$ defined on the interval $x\in\left[-1,1\right]$
can be expanded as

$\displaystyle f\left(x\right)$
$\displaystyle=\frac{1}{\pi\sqrt{1-x^{2}}}\left[\mu_{0}+2\sum_{n=1}^{\infty}\mu_{n}T_{n}(x)\right]\,,$
(S3) $\displaystyle\mu_{n}$
$\displaystyle=\int_{-1}^{1}f\left(x\right)T_{n}\left(x\right)dx\,.$ (S4)

The Gibbs oscillation of the Chebyshev expansion can be suppressed by
introducing kernels

$\displaystyle f\left(x\right)$
$\displaystyle=\frac{1}{\pi\sqrt{1-x^{2}}}\left[\tilde{\mu}_{0}+2\sum_{n=1}^{N_{\mathcal{C}}-1}\tilde{\mu}_{n}T_{n}\left(x\right)\right]\,,$
(S5a) $\displaystyle\tilde{\mu}_{n}$ $\displaystyle=\mu_{n}g_{n}\,.$ (S5b)

Here $N_{\mathcal{C}}$ is the maximum degree of the Chebyshev polynomials,
$\mu_{n}$’s are the Chebyshev moments, and $g_{n}$’s are the kernel
coefficients. The Jackson kernel is commonly employed and suitable for most
applications,

$g_{n}=\frac{1}{N_{\mathcal{C}}+1}\left[\left(N_{\mathcal{C}}-n+1\right)\cos\left(\frac{\pi
n}{N_{\mathcal{C}}+1}\right)+\sin\left(\frac{\pi
n}{N_{\mathcal{C}}+1}\right)\cot\left(\frac{\pi}{N_{\mathcal{C}}+1}\right)\right]\,.$
(S6)

In practice, we can obtain $f\left(x\right)$ at discrete points using the
discrete Fourier transform. For example, we can choose

$x_{k}=\cos\left(\frac{k+\frac{1}{2}}{N_{p}}\pi\right)\,,\qquad
k=0,1,\cdots,N_{p}-1\,.$ (S7)

Here $N_{p}$ is the number of discrete points for $x$ and we can choose
$N_{p}\geq N_{\mathcal{C}}$ to utilize all the available Chebyshev moments
$\mu_{n}$. Then we have

$\gamma_{k}\equiv\pi\sqrt{1-x_{k}^{2}}\,f\left(x_{k}\right)=\tilde{\mu}_{0}+2\sum_{n=1}^{N_{\mathcal{C}}-1}\tilde{\mu}_{n}\cos\left[\frac{\pi
n\left(k+\frac{1}{2}\right)}{N_{p}}\right]\,,$ (S8)

which can be obtained simultaneously with the fast-Fourier-transform
algorithm. Another choice of $x$ discretization is given by

$\displaystyle x_{k}$
$\displaystyle=\cos\left(\frac{k\pi}{N_{p}-1}\right)\,,\qquad
k=0,1,\cdots,N_{p}-1\,,$ (S9) $\displaystyle\gamma_{k}$
$\displaystyle=\pi\sqrt{1-x_{k}^{2}}\,f\left(x_{k}\right)=\tilde{\mu}_{0}+2\sum_{n=1}^{N_{\mathcal{C}}-1}\tilde{\mu}_{n}\cos\left(\frac{\pi
nk}{N_{p}-1}\right)\,.$ (S10)

The advantage of the second one is that it accesses $x_{k}=0$ exactly with odd
$N_{p}$, which is desired when we are interested in the scaling of the zero-
energy DOS. The two choices of $x_{k}$ corresponds to the two kinds of
Chebyshev-Gauss quadrature in numerical integration of functions involving
$f\left(x\right)$. With the first choice, we have

$\int_{-1}^{1}f\left(x\right)g\left(x\right)dx\approx\frac{1}{N_{p}}\sum_{k=0}^{N_{p}-1}\gamma_{k}\,g\left(x_{k}\right)\,.$
(S11)

###  I.2 Density of states

The density of states of a system can be obtained with KPM efficiently. The
global DOS of a system is given by

$\rho\left(E\right)=\frac{1}{V}\sum_{k}\delta\left(E-E_{k}\right)\,,$ (S12)

with $V$ the volume of the system. Here and below, we let $\ket{i}$ and
$\ket{k}$ denote the basis vectors in the position and energy bases,
respectively. In general, the spectrum of a system is bounded by
$E_{\mathsf{min}}$ and $E_{\mathsf{max}}$, which can be obtained via the
Lanczos algorithm. Then we can rescale the Hamiltonian,

$\tilde{H}=aH+b\,,$ (S13)

with

$\displaystyle a$
$\displaystyle=\frac{2\left(1-\epsilon\right)}{E_{\mathsf{max}}-E_{\mathsf{min}}}\,,$
(S14) $\displaystyle b$
$\displaystyle=-\frac{\left(1-\epsilon\right)\left(E_{\mathsf{max}}+E_{\mathsf{min}}\right)}{E_{\mathsf{max}}-E_{\mathsf{min}}}\,,$
(S15)

with $0<\epsilon\ll 1$ to guarantee the numerical stability. Thus we have

$\rho\left(E\right)=\frac{a}{V}\sum_{k}\delta\left(\tilde{E}-\tilde{E}_{k}\right)\,,$
(S16)

with $E=\left(\tilde{E}-b\right)/a$. The Chebyshev moments are then given by

$\mu_{n}=\int_{-1}^{1}d\tilde{E}\,\rho\left(E\right)T_{n}\left(\tilde{E}\right)=\frac{a}{V}\sum_{k}T_{n}\left(\tilde{E}_{k}\right)=\frac{a}{V}\mathsf{Tr}\left[T_{n}\left(\tilde{H}\right)\right]\,.$
(S17)

The trace can be efficiently evaluated by averaging over a few independent
random vectors,

$\mu_{n}\approx\frac{a}{VR}\sum_{r=1}^{R}\bra{\psi_{r}}T_{n}\left(\tilde{H}\right)\ket{\psi_{r}}\,.$
(S18)

Here $R$ is the number of random vectors $\left\\{\psi_{ri}\right\\}$
satisfying

$\displaystyle\left\langle\psi_{ri}\right\rangle$ $\displaystyle=0\,,$ (S19a)
$\displaystyle\left\langle\psi_{ri}\,\psi_{r^{\prime}j}\right\rangle$
$\displaystyle=0\,,$ (S19b)
$\displaystyle\left\langle\psi_{ri}^{*}\,\psi_{r^{\prime}j}\right\rangle$
$\displaystyle=\delta_{rr^{\prime}}\delta_{ij}\,.$ (S19c)

The angle brackets $\left\langle\cdots\right\rangle$ denote the average over
the ensemble of random vectors. $\psi_{ri}$ is the $i$-th component of the
$r$-th random vector, subject to some particular distribution. In our
numerics, we use the uniform distribution for $\psi_{ri}$ or the uniform
distribution for $\theta$ with $\psi_{ri}=e^{i\theta}$. The vector
$\ket{\psi^{n}}\equiv T_{n}\left(\tilde{H}\right)\ket{\psi}$ can be obtained
recursively using the relation in Eq. (S2),

$\displaystyle\ket{\psi^{0}}$
$\displaystyle=T_{0}\left(\tilde{H}\right)\ket{\psi}=\ket{\psi}\,,$ (S20a)
$\displaystyle\ket{\psi^{1}}$
$\displaystyle=T_{1}\left(\tilde{H}\right)\ket{\psi}=\tilde{H}\ket{\psi}\,,$
(S20b) $\displaystyle\vdots$ $\displaystyle\phantom{=}\,\,\vdots$
$\displaystyle\ket{\psi^{n+1}}$
$\displaystyle=T_{n+1}\left(\tilde{H}\right)\ket{\psi}=2\tilde{H}\ket{\psi^{n}}-\ket{\psi^{n-1}}\,.$
(S20c)

Using the product relation of Chebyshev polynomials,

$2T_{m}\left(x\right)T_{n}\left(x\right)=T_{m+n}\left(x\right)+T_{m-n}\left(x\right)\,.$
(S21)

Then we have

$\displaystyle\mu_{2n}$
$\displaystyle=\frac{2a}{V}\braket{\psi^{n}}{\psi^{n}}-\mu_{0}\,,$ (S22)
$\displaystyle\mu_{2n+1}$
$\displaystyle=\frac{2a}{V}\braket{\psi^{n+1}}{\psi^{n}}-\mu_{1}\,.$ (S23)

The global DOS can be obtained with discrete Fourier transform of
$\tilde{\mu}_{n}=\mu_{n}g_{n}$. Typically, averaging a few ($\sim 10$) random
vectors is enough to evaluate the global DOS accurately for large enough
systems. As a result, the computational effort required for calculating the
global DOS is of the order of $NN_{\mathcal{C}}$ for a sparse Hamiltonian,
with $N$ denoting the dimension of the Hilbert space and $N_{\mathcal{C}}$ the
degree of the Chebyshev polynomials.

###  I.3 Chebyshev filter methods for wave functions

####  I.3.1 Iterative filter method

Various numerical methods based on the KPM have been proposed to obtain the
eigenvalues and eigenvectors of a large sparse matrix. One of them is the
iterative Chebyshev filter diagonalization method proposed in Ref. [2]. The
method aims to obtain the eigenvalue and eigenvectors in a small interval
$\left[E_{1},E_{2}\right]$. It starts with an ensemble of random vectors,
whose dimension is greater than the number of states in the interval estimated
from the DOS. At each step of the iteration, the vectors are filtered with the
kernel polynomial expansion of a filter function. Then the set of filtered
vectors are orthogonalized to obtain an orthonormal basis in the interval. The
eigenvalues, eigenvectors and residual vectors are obtained with Rayleigh-Ritz
method. The steps can be repeated until the residual of all the eigenpairs in
the interval satisfy the preset criteria. Various filter functions can be used
for the method. One of them is the rectangular function,

$W\left(x\right)=\Theta\left(x-E_{1}\right)\Theta\left(E_{2}-x\right)\,.$
(S24)

A vector $\ket{\psi}$ can be filtered to amplify its components in the
interval via

$\ket{\psi^{\mathsf{F}}}=\sum_{n=0}^{N_{\mathcal{C}}}c_{n}\,g_{n}\,T\left(\tilde{H}\right)\ket{\psi}\,.$
(S25)

Here $g_{n}$ is the kernel and $c_{n}$ is an expansion coefficient. The latter
can be obtained exactly with

$c_{n}=\frac{2-\delta_{n0}}{\pi}\int_{-1}^{1}\frac{W\left(\tilde{x}\right)T_{n}\left(\tilde{x}\right)}{\sqrt{1-x^{2}}}dx\,.$
(S26)

####  I.3.2 Filter and shake method

Instead of the iterative Chebyshev filter method, we use the filter and shake
method proposed in Ref. [3] to evaluate the eigenvalues and eigenvectors in an
interval. The method uses the so-called exponential semicircle filter to
obtain a vector $\ket{\psi_{0}}$ with numerically vanishing overlap with
eigenvectors outside of the interval. Then the vector is “shaken” by Chebyshev
iteration to obtain a set of vectors, which are superpositions of eigenvectors
in the desired interval. The eigenvalues, eigenvectors and residual vectors
are then obtained via the Rayleigh-Ritz method using the orthonormal basis
vectors.

The exponential semicircle filter exploits the fact that the Chebyshev
polynomials $T_{n}(x)$ diverge exponentially with $n$ when $x$ is outside of
the defining interval $\left[-1,1\right]$. When $\left|x\right|>1$, we have

$\arccos x=\begin{cases}i\cosh^{-1}x&x>1\,,\\\
\pi-i\cosh^{-1}x&x<-1\,.\end{cases}$ (S27)

Thus, we have

$T_{n}\left(x\right)=\cos\left(n\arccos
x\right)=\begin{cases}\cosh\left(n\cosh^{-1}x\right)\,,&x>1\,,\\\
\left(-1\right)^{n}\cosh\left(n\cosh^{-1}x\right)\,,&x<-1\,.\end{cases}$ (S28)

Using $\cosh^{-1}\left(1+x\right)\approx\sqrt{2x}$ with $x>0$, we have

$\left|T_{n}\left(x\right)\right|\sim\frac{1}{2}e^{n\sqrt{2\left(\left|x\right|-1\right)}}\,,$

with $\left|x\right|>1$ and $n\sqrt{2\left(\left|x\right|-1\right)}\gg 1$. We
can see that $\left|T_{n}\left(x\right)\right|$ grows exponentially with the
degree of Chebyshev iteration when $\left|x\right|>1$, while
$\left|T_{n}\left(x\right)\right|\ll 1$ for $\left|x\right|<1$. Therefore, the
Chebyshev iteration itself provides an efficient filter by projecting the
desired energy window out of the defining range of the Chebyshev polynomials.

Consider a system with Hamiltonian $H$ and spectrum in the range
$\left[E_{\mathsf{min}},E_{\mathsf{max}}\right]$. We seek the eigenvalues and
eigenstates in an energy window $\left[E_{0}-\delta,E_{0}+\delta\right]$. We
can use a parabolic function $F\left(H\right)$ [See Fig. S1] to project the
energy spectrum of $H$ so that $F\left(E\right)>1$ when
$E\in\left[E_{0}-\delta,E_{0}+\delta\right]$ and $F\left(E\right)<1$ when
$E\notin\left[E_{0}-\delta,E_{0}+\delta\right]$. The parabolic function
$F\left(x\right)$ is given by

$\displaystyle F\left(H\right)$
$\displaystyle=a^{2}\left(H-E_{0}\right)^{2}-c^{2}=\left[a\left(H-E_{0}\right)+c\right]\left[a\left(H-E_{0}\right)-c\right]\,,$
(S29a) $\displaystyle a$
$\displaystyle=\sqrt{\frac{2-\epsilon}{\mathcal{E}^{2}-\delta^{2}}}\,,\qquad
c=\sqrt{\frac{\left(2-\epsilon\right)\delta^{2}}{\mathcal{E}^{2}-\delta^{2}}+1}\,,$
(S29b) $\displaystyle\mathcal{E}$
$\displaystyle=\max\left\\{\left|E_{\mathsf{min}}-E_{0}\right|,\left|E_{\mathsf{max}}-E_{0}\right|\right\\}\,.$
(S29c)

Figure S1: Illustration of the exponential semicircle filter.

Performing Chebyshev iteration on $F\left(H\right)$ on a random state
$\ket{\psi}$, we have

$\displaystyle\ket{\psi^{0}}$ $\displaystyle=\ket{\psi}\,,$ (S30a)
$\displaystyle\ket{\psi^{1}}$ $\displaystyle=F\left(H\right)\ket{\psi^{0}}\,,$
(S30b) $\displaystyle\cdots$ $\displaystyle\ket{\psi^{n}}$
$\displaystyle=2F\left(H\right)\ket{\psi^{n}}-\ket{\psi^{n-1}}\,.$ (S30c)

The initial random vector can be written as

$\ket{\psi}=\sum_{m=1}^{N}c_{m}\ket{m}\,,$ (S31)

with $\ket{m}$ denoting the eigenstate of $H$ with energy $E_{m}$ and $c_{m}$
a random number with $\left|c_{m}\right|\sim\frac{1}{\sqrt{N}}$. At order $n$
of the Chebyshev iteration, we have

$\displaystyle\ket{\psi^{n}}$ $\displaystyle=\sum_{m}\tilde{c}_{m}\ket{m}\,,$
(S32a) $\displaystyle\tilde{c}_{m}$
$\displaystyle=c_{m}T_{n}\left[F\left(E_{m}\right)\right]\,.$ (S32b)

When $E_{m}$ falls out of the interval
$\left[E_{0}-\delta,E_{0}+\delta\right]$, $T_{n}=\cos\left[n\arccos
F\left(E_{m}\right)\right]\in\left[-1,1\right]$ and $\tilde{c}_{m}$ does not
increase with the degree of Chebyshev expansion $n$. However, when
$E_{m}\in\left[E_{0}-\delta,E_{0}+\delta\right]$, we have

$\displaystyle T_{n}\left[F\left(E_{m}\right)\right]$
$\displaystyle=\left(-1\right)^{n}\cosh\left[n\cosh^{-1}F\left(E_{m}\right)\right]$
$\displaystyle\approx\left(-1\right)^{n}\frac{1}{2}\,e^{n\sqrt{2\left[-a^{2}\left(E_{m}-E_{0}\right)^{2}+c^{2}-1\right]}}$
$\displaystyle=\frac{\left(-1\right)^{n}}{2}e^{\frac{2n}{\sqrt{\mathcal{E}^{2}-\delta^{2}}}\sqrt{\delta^{2}-\left(E_{m}-E_{0}\right)^{2}}}\,.$
(S33)

Thus $\tilde{c}_{m}$ grows as the exponential of a semicircle when
$E_{m}\in\left[E_{0}-\delta,E_{0}+\delta\right]$ and $n$ is large. Focusing on
the interval within the half-width of the semi-circle, we have

$\left|\frac{\tilde{c}_{m}}{c_{m}}\right|>\frac{1}{2}e^{\frac{n\delta}{\mathcal{E}}}\,.$
(S34)

It turns out that the exponential semicircle filter is much more efficient
than the rectangular one and it can be used to accelerate the iterative method
in Sec. I.3.1. Moreover, Ref. [3] proposed that a complete basis for the
states in the desired interval can be obtained by shaking the filtered state
with Chebyshev iteration. In numerics, a double-precision float number has
about $16$ significant figures. When
$n>\frac{\mathcal{E}}{\delta}\log\left(10^{16}\right)$, the $\ket{\psi^{n}}$
contains numerically vanishing components of eigenstates outside of the
desired interval and

$\ket{\psi^{n}}\approx\sum_{E_{0}-\delta<E_{m}<E_{0}+\delta}\tilde{c}_{m}\ket{m}\,.$
(S35)

Now we take Chebyshev iteration with scaled Hamiltonian $\tilde{H}$ and
$\ket{\psi^{n}}$ as initial state. We have

$\displaystyle\ket{\phi^{0}}$ $\displaystyle=\ket{\psi^{n}}\,,$ (S36a)
$\displaystyle\ket{\phi^{1}}$ $\displaystyle=\tilde{H}\ket{\phi^{0}}\,,$
(S36b) $\displaystyle\cdots$ $\displaystyle\ket{\phi^{l+1}}$
$\displaystyle=2\tilde{H}\ket{\phi^{l}}-\ket{\phi^{l-1}}\,.$ (S36c)

Then

$\ket{\phi^{l}}\approx\sum_{E_{0}-\delta<E_{m}<E_{0}+\delta}\tilde{c}_{m}\cos\left(l\arccos\tilde{E}_{m}\right)\ket{m}\,.$
(S37)

Thus the amplitude of $\ket{m}$ in the $\ket{\phi^{l}}$ oscillates with its
energy and the phase is given by

$l\arccos\tilde{E}_{m}\approx
l\arccos{\tilde{E}_{0}}-\frac{l\left(\tilde{E_{m}}-\tilde{E_{0}}\right)}{\sqrt{1-\tilde{E}_{0}^{2}}}\,,$
(S38)

with $\left|\tilde{E_{m}}-\tilde{E_{0}}\right|\ll 1$. The phases of the states
in the desired interval now distribute across the interval
$l\arccos{\tilde{E}_{0}}-\frac{l}{\sqrt{1-\tilde{E}_{0}^{2}}}\left[-\tilde{\delta},\tilde{\delta}\right]$
with $\tilde{\delta}\ll 1$. When
$l\tilde{\delta}/{\sqrt{1-\tilde{E}_{0}^{2}}}\approx\frac{\pi}{2}\mod\pi$, the
phases of the all the states in the desired interval runs across an interval
of length $\pi$. If the energies of the states uniformly distributed in the
desired interval, the phases are also uniformed distributed in the length
$\pi$ interval around $l\pi/2$ if $\tilde{E_{0}}\approx 0$. Thus for each $l$,
we have a different superposition of eigenstates in the desired interval. Thus
we can keep one state $\ket{\phi^{l}}$ for the basis for every
$\lfloor\frac{\pi\sqrt{1-\tilde{E}_{0}^{2}}}{2\tilde{\delta}}\rfloor$
Chebyshev iteration. Moreover, observing that
$\arccos\tilde{E}_{m}\approx\frac{\pi}{2}+\tilde{E}$ when
$\tilde{E}_{m}\approx 0$ and we have

$\cos\left(l\arccos\tilde{E}_{m}\right)\approx\cos\left(\frac{\pi
l}{2}+l\tilde{E}\right)=\begin{cases}\left(-1\right)^{\frac{l}{2}}\cos\left(l\tilde{E}\right)&l\in\text{even}\,,\\\
\left(-1\right)^{\frac{l+1}{2}}\sin\left(l\tilde{E}\right)&l\in\text{odd}\,.\end{cases}$
(S39)

Thus we can keep both $\ket{\phi^{l}}$ and $\ket{\phi^{l+1}}$ in the basis
when $l=\lfloor\frac{j\pi\mathcal{E}}{2\delta}\rfloor$. Shaking
$\ket{\psi^{n}}$ with Chebyshev iteration, we can generate a basis for the
states in the desired interval. The eigenvalue and eigenstates can then be
generated with this basis via the Rayleigh-Ritz method.

Now we summarize the filter and shake method as follows:

  1. 1.

Calculate $E_{\mathsf{min}}$ and $E_{\mathsf{max}}$ of $H$ with the Lanczos
method.

  2. 2.

Calculate the DOS of the system using KPM and estimate the number of states
$N_{T}$ in the desired interval $\left[E_{0}-\delta,E_{0}+\delta\right]$.

  3. 3.

Generate the filtered state $\ket{\psi^{n}}$ from a random state using the
exponential semicircle filter and
$n>\frac{\mathcal{E}}{\delta}\log\left(10^{16}\right)$ with
$\mathcal{E}=\max\left\\{\left|E_{\mathsf{min}}-E_{0}\right|,\left|E_{\mathsf{max}}-E_{0}\right|\right\\}$.

  4. 4.

Shake $\ket{\psi^{n}}$ with Chebyshev iteration and
$\ket{\phi^{0}}=\ket{\psi^{n}}$. Keep $\ket{\phi^{l}}$ and $\ket{\phi^{l+1}}$
to the basis with $\frac{l\delta}{\mathcal{E}}\approx\frac{\pi}{2}$ (assuming
$E_{0}\approx 0$) and construct a basis $\left\\{\Phi_{j}\right\\}$ containing
$N_{S}>N_{T}$ states. $N_{S}=1.5N_{T}$ or $N_{S}=2N_{T}$ is usually
sufficient.

  5. 5.

Orthogonalize the set of basis vector $\left\\{\Phi_{j}\right\\}$ and obtain
an orthonormal basis containing $N_{B}\leq N_{S}$ states. This can be done
with the SVQB algorithm.

  6. 6.

Calculate the eigenvalue, eigenvectors and residual vectors using the
orthonormal basis via Rayleigh-Ritz method.

  7. 7.

Keep the states with $E\in\left[E_{0}-\delta,E_{0}+\delta\right]$ and residual
smaller than the preset criteria $\eta$.

The method works best when the eigenvalues distribute uniformly across the
desired interval. When there are degeneracies, we should generate a number of
filtered states $\ket{\psi^{n}}$ with independent random initial vectors and
generate a subset of the basis for each filtered state. When the number of
filtered states is greater than the degeneracy of states in the interval, the
filter and shake method works well. Moreover, the method is more efficient by
using an appropriate number of filtered states, balancing the computational
time for the filter step and the shake step.

###  I.4 Self-consistent BCS theory based on KPM

For an inhomogeneous system with interactions, mean-field order parameters
also become inhomogeneous. The spatially varying order parameter generally
prohibits analytical evaluation and exact diagonalization (ED) is often
required to perform the calculation. However, it takes a lot of memory and
computational time to diagonalize a large Hamiltonian matrix and the system
size accessible is limited at the order of $10^{4}$ lattice sites. The KPM
based self-consistent mean field theory allows us to access larger systems
with much less cost of memory and computational time.

####  I.4.1 Self-consistent BdG theory for a lattice model

Now we consider a lattice model with attractive Hubbard interaction,

$H=-t\sum_{\left\langle
ij\right\rangle\sigma}\left(c_{i\sigma}^{\dagger}c_{j\sigma}+H.c.\right)-U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\mu\sum_{i}n_{i}\,.$
(S40)

It can be solved by the mean-field ansatz,

$\Delta_{i}=-U\left\langle
c_{i\downarrow}c_{i\uparrow}\right\rangle\,,\qquad\left\langle
n_{i\sigma}\right\rangle=\left\langle
c_{i\sigma}^{\dagger}c_{i\sigma}\right\rangle=\frac{1}{2}\left\langle
n\right\rangle\,.$ (S41)

The BdG mean-field Hamiltonian can be written as

$H_{\mathsf{eff}}=-t\sum_{\left\langle
i,j\right\rangle\sigma}\left(c_{i\sigma}^{\dagger}c_{j\sigma}+H.c.\right)+\sum_{i}\left(V_{i}-\tilde{\mu}_{i}\right)n_{i\sigma}+\sum_{i}\left(\Delta_{i}c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}+\Delta_{i}^{*}c_{i\uparrow}c_{i\downarrow}\right)\,.$
(S42)

Here $\tilde{\mu}=\mu+U\left<n\right>/2$ incorporates the Hartree shift.
$H_{\mathsf{eff}}$ can be diagonalized by the Bogoliubov transformation,

$\displaystyle c_{i\uparrow}$
$\displaystyle=\sum_{n}^{\prime}\left[u_{n}^{i}\alpha_{n}-\left(v_{n}^{i}\right)^{*}\beta_{n}^{\dagger}\right]\,,\qquad
c_{i\downarrow}=\sum_{n}^{\prime}\left[\left(v_{n}^{i}\right)^{*}\alpha_{n}^{\dagger}+u_{n}^{i}\beta_{n}\right]\,,$
(S43) $\displaystyle\left[\begin{array}[]{c}c_{\uparrow}\\\
\left(c_{\downarrow}^{\dagger}\right)^{\mathsf{T}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{cc}u&-v^{*}\\\
v&u^{*}\end{array}\right]\left[\begin{array}[]{c}\alpha\\\
\left(\beta^{\dagger}\right)^{\mathsf{T}}\end{array}\right]\,.$ (S50)

with $\alpha$ and $\beta{}^{\dagger}$ quasiparticle operators, $\alpha_{i}$
and $\beta_{i}$ correspond to quasiparticles with the same energy and
$\sum^{\prime}$ sums over all the positive energy states. $u,v$ denote
$N\times N$ matrices and
$\left(u_{1n},\cdots,u_{Nn},v_{1n},\cdots,v_{Nn}\right)^{\mathsf{T}}$ is the
$n$-the positive energy $E_{n}$ eigenstate of the $2N\times 2N$ BdG
Hamiltonian,

$\left[\begin{array}[]{cc}\hat{h}_{0}&\Delta\\\
\Delta^{\dagger}&-\hat{h}_{0}^{\mathsf{T}}\end{array}\right]\left[\begin{array}[]{c}u_{n}\\\
v_{n}\end{array}\right]=E_{n}\left[\begin{array}[]{c}u_{n}\\\
v_{n}\end{array}\right]\,,$ (S51)

with $u_{n}=\left(u_{1n},\cdots,u_{Nn}\right)^{\mathsf{T}}$ and
$v_{n}=\left(v_{1n},v_{2n},\cdots,v_{Nn}\right)^{\mathsf{T}}$,
$(\hat{h}_{0})_{ij}=-t\left(\delta_{i,j+1}+\delta_{i+1,j}\right)-t\left(\delta_{i1}\delta_{jN}+\delta_{iN}\delta_{j1}\right)+\left(V_{i}-\tilde{\mu}_{i}\right)\delta_{ij}$
and $\Delta_{ij}=\Delta_{i}\delta_{ij}$. Thus we have

$H_{\mathsf{eff}}=\sum_{n}^{\prime}E_{n}\left(\alpha_{n}^{\dagger}\alpha_{n}-\beta_{n}\beta_{n}^{\dagger}\right)\,.$
(S52)

The self-consistent conditions at $T=0$ are then given by

$\Delta_{i}=U\sum_{n}^{\prime}u_{n}^{i}\left(v_{n}^{i}\right)^{*}\,,\qquad\left\langle
n_{i}\right\rangle=2\sum_{n}^{\prime}\left|v_{n}^{i}\right|^{2}\,.$ (S53)

At finite temperature, we have

$\displaystyle\Delta_{i}$ $\displaystyle=-U\left\langle
c_{i\downarrow}c_{i\uparrow}\right\rangle=-U\sum_{mn}^{\prime}\left\langle\left[\left(v_{m}^{i}\right)^{*}\alpha_{m}^{\dagger}+u_{n}^{i}\beta_{m}\right]\left[u_{n}^{i}\alpha_{n}-\left(v_{n}^{i}\right)^{*}\beta_{n}^{\dagger}\right]\right\rangle$
$\displaystyle=U\sum_{n}^{\prime}u_{n}^{i}\left(v_{n}^{i}\right)^{*}\left[1-2f\left(E_{n}\right)\right]\,$
(S54)

and

$\displaystyle\left\langle n_{i}\right\rangle$ $\displaystyle=\left\langle
c_{i\uparrow}^{\dagger}c_{i\uparrow}\right\rangle+\left\langle
c_{i\downarrow}^{\dagger}c_{i\downarrow}\right\rangle$
$\displaystyle=2\sum_{n}^{\prime}\left[\left|u_{n}\right|^{2}f\left(E_{n}\right)+\left|v_{n}^{i}\right|^{2}\left(1-f\left(E_{n}\right)\right)\right]\,.$
(S55)

Here $f\left(E_{n}\right)=\frac{1}{e^{\beta E_{n}}+1}$ is the Fermi
distribution function for the quasiparticles with energy $E_{n}$.

It turns out to be convenient to write the eigenstates in the following form,

$\left(\begin{array}[]{cc}u&-v^{*}\\\
v&u^{*}\end{array}\right)=\left(\begin{array}[]{cccccc}u_{1}^{i}&\cdots&u_{N}^{i}&-\left(v_{1}^{i}\right)^{*}&\cdots&-\left(v_{N}^{i}\right)^{*}\\\
v_{n}^{i}&\cdots&v_{N}^{i}&\left(u_{1}^{i}\right)^{*}&\cdots&\left(u_{N}^{i}\right)^{*}\end{array}\right)=\left(\begin{array}[]{cccccc}u_{1}^{i}&\cdots&u_{N}^{i}&u_{-1}^{i}&\cdots&u_{-N}^{i}\\\
v_{n}^{i}&\cdots&v_{N}^{i}&v_{-1}^{i}&\cdots&v_{-1}^{i}\end{array}\right)\,,\,,$
(S56)

with the first $N$ columns for the $N$ positive energy states and the last $N$
columns for the negative energy states. Denote the eigenstate matrix as

$\left(\begin{array}[]{cccccc}u_{1}^{i}&\cdots&u_{N}^{i}&u_{-1}^{i}&\cdots&u_{-N}^{i}\\\
v_{n}^{i}&\cdots&v_{N}^{i}&v_{-1}^{i}&\cdots&v_{-1}^{i}\end{array}\right)\,,$
(S57)

with $u_{-n}^{i}=-\left(v_{n}^{i}\right)^{*}$ and
$v_{-n}^{i}=\left(u_{n}^{i}\right)^{*}$. Here

$u_{n}^{i}=\braket{i\uparrow}{n}\,,\qquad
v_{n}=\braket{i\downarrow}{n}\,,\qquad 1\leq i\leq N,\qquad-N\leq n\leq
N\,,n\neq 0\,.$ (S58)

Here $\ket{i\uparrow}=c_{i\uparrow}^{\dagger}\ket{0}$ and
$\ket{i\downarrow}=c_{i\downarrow}\ket{0}$ with $\ket{0}$ being the BCS ground
state without quasiparticle excitations. The corresponding Bogoliubov
transformation is given by

$\displaystyle c_{i\uparrow}$
$\displaystyle=\sum_{n}^{\prime}\left(u_{n}^{i}\alpha_{n}+u_{-n}^{i}\alpha_{-n}^{\dagger}\right)\,,\qquad
c_{i\downarrow}=\sum_{n}^{\prime}\left[\left(v_{n}^{i}\right)^{*}\alpha_{n}^{\dagger}+\left(v_{-n}^{i}\right)^{*}\alpha_{-n}\right]\,,$
(S59) $\displaystyle\left[\begin{array}[]{c}c_{\uparrow}\\\
\left(c_{\downarrow}^{\dagger}\right)^{\mathsf{T}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{cc}u_{+}&u_{-}\\\
v_{+}&v_{-}\end{array}\right]\left[\begin{array}[]{c}\alpha\\\
\left(\alpha^{\dagger}\right)^{\mathsf{T}}\end{array}\right]\,.$ (S66)

In summary, we have

$\displaystyle\Delta_{i}$
$\displaystyle=-U\sum_{n}u_{n}^{i}\left(v_{n}^{i}\right)^{*}f\left(E_{n}\right)\,,\qquad
E_{-n}=-E_{n}\,,$ (S67a) $\displaystyle n_{i}$
$\displaystyle=2\sum_{n}\left|u_{n}^{i}\right|^{2}f\left(E_{n}\right)\,,$
(S67b) $\displaystyle u_{n}^{i}$ $\displaystyle=\braket{i\uparrow}{n}\,,\qquad
v_{n}=\braket{i\downarrow}{n}\,.$ (S67c)

These expressions for $\Delta_{i}$ and $n_{i}$ are easier to implement at
finite temperature and are suitable to be generalized to the KPM method.

####  I.4.2 KPM Implementation of the self-consistent BdG theory

The BdG Hamiltonian has particle-hole symmetry and its spectrum is symmetric
with respective to $0$. Thus we can rescale the BdG Hamiltonian with

$\tilde{H}=aH\,,\qquad a=\frac{1-\epsilon}{\left|E_{g}\right|}\,,$ (S68)

with $E_{g}$ being the smallest eigenenergy of the BdG Hamiltonian.

We define the pairing and density spectral functions,

$\displaystyle\Delta_{i}\left(E\right)$
$\displaystyle=-U\sum_{n}u_{n}^{i}\left(v_{n}^{i}\right)^{*}\delta\left(E-E_{n}\right)=-U\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\downarrow}\delta\left(E-E_{n}\right)\,,$
(S69) $\displaystyle\rho_{i}\left(E\right)$
$\displaystyle=2\sum_{n}\left|u_{n}^{i}\right|^{2}\delta\left(E-E_{n}\right)=2\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\uparrow}\delta\left(E-E_{n}\right)\,.$
(S70)

Therefore, the local pairing amplitude and local density can be obtained via

$\displaystyle\Delta_{i}$ $\displaystyle=\int
dE\,\Delta_{i}\left(E\right)f\left(E\right)\,,$ (S71) $\displaystyle n_{i}$
$\displaystyle=\int dE\,\rho_{i}\left(E\right)f\left(E\right)\,.$ (S72)

The rescaled spectral functions of the scaled Hamiltonian are given by

$\displaystyle\tilde{\Delta}_{i}\left(\tilde{E}\right)$
$\displaystyle=-U\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\downarrow}\delta\left(\tilde{E}-\tilde{E}_{n}\right)=\frac{1}{a}\Delta_{i}\left(E=a\tilde{E}\right)\,,$
(S73) $\displaystyle\tilde{\rho}_{i}\left(\tilde{E}\right)$
$\displaystyle=2\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\uparrow}\delta\left(\tilde{E}-\tilde{E}_{n}\right)=\frac{1}{a}\rho_{i}\left(E=a\tilde{E}\right)\,.$
(S74)

The Chebyshev moments of the rescaled spectral functions are given by

$\displaystyle\mu_{m}^{\Delta_{i}}$
$\displaystyle=\int_{-1}^{1}\tilde{\Delta}_{i}\left(\tilde{E}\right)T_{n}\left(\tilde{E}\right)d\tilde{E}=-U\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\downarrow}T_{m}\left(\tilde{E}_{n}\right)$
$\displaystyle=-U\left(\bra{i\downarrow}T_{m}\left(\tilde{H}\right)\ket{i\uparrow}\right)^{*}\,,$
(S75) $\displaystyle\mu_{m}^{\rho_{i}}$
$\displaystyle=\int_{-1}^{1}\tilde{\rho}_{i}\left(\tilde{E}\right)T_{n}\left(\tilde{E}\right)d\tilde{E}=2\sum_{n}\braket{i\uparrow}{n}\braket{n}{i\uparrow}T_{m}\left(\tilde{E}_{n}\right)$
$\displaystyle=2\bra{i\uparrow}T_{m}\left(\tilde{H}\right)\ket{i\uparrow}\,.$
(S76)

The Chebyshev iteration can be performed starting with $\ket{i\uparrow}$,

$\displaystyle\left|\psi_{0}\right>$ $\displaystyle=\left|i\uparrow\right>\,,$
(S77a) $\displaystyle\left|\psi_{1}\right>$
$\displaystyle=\tilde{H}\left|\psi_{0}\right>\,,$ (S77b) $\displaystyle\cdots$
$\displaystyle\left|\psi_{m+1}\right>$
$\displaystyle=2\tilde{H}\left|\psi_{m}\right>-\left|\psi_{m-1}\right>\,.$
(S77c)

Thus we have the Chebyshev moments,

$\displaystyle\mu_{m}^{\Delta_{i}}$
$\displaystyle=-U\braket{i\downarrow}{\psi_{m}}\,,$ (S78)
$\displaystyle\mu_{m}^{\rho_{i}}$
$\displaystyle=2\braket{i\uparrow}{\psi_{m}}\,.$ (S79)

The Chebyshev moments can be obtained by reading spin and down components of
$\ket{\psi_{m}}$ without performing the inner product.

With $N_{\mathcal{C}}$ Chebyshev moments of $\Delta_{i}$ and $\mu_{i}$
available, we can evaluate the scaled spectral functions at points,

$\tilde{E}_{k}=\cos\left[\frac{\pi\left(k+\frac{1}{2}\right)}{N_{\mathsf{p}}}\right]\,,\qquad
k=0,1,\cdots,N_{\mathsf{p}}-1\,.$ (S80)

Usually we consider $N_{p}\geq N_{\mathcal{C}}$ (e.g.
$N_{p}=2N_{\mathcal{C}}$) and $\tilde{\Delta}_{i}\left(\tilde{E}_{k}\right)$
can be evaluated through the discrete cosine transform from the 1D array
$\left\\{\tilde{\mu}_{0},\tilde{\mu}_{1},\cdots,\tilde{\mu}_{N_{\mu}-1}\underset{N_{p}-N_{\mathcal{C}}}{,\underbrace{0,\cdots,0}}\right\\}$,

$\displaystyle\gamma_{k}^{\Delta_{i}}$
$\displaystyle=\pi\sqrt{1-\tilde{E}_{k}^{2}}\tilde{\Delta}_{i}\left(\tilde{E}_{k}\right)=\tilde{\mu}_{0}^{\Delta_{i}}+2\sum_{m=1}^{N_{\mathcal{C}}-1}\tilde{\mu}_{n}^{\Delta_{i}}\cos\left[\frac{\pi
n\left(k+1/2\right)}{N_{p}}\right]\,,$ (S81)
$\displaystyle\gamma_{k}^{\rho_{i}}$
$\displaystyle=\pi\sqrt{1-\tilde{E}_{k}^{2}}\tilde{\rho}_{i}\left(\tilde{E}_{k}\right)=\tilde{\mu}_{0}^{\rho_{i}}+2\sum_{m=1}^{N_{\mathcal{C}}-1}\tilde{\mu}_{n}^{\Delta_{i}}\cos\left[\frac{\pi
n\left(k+1/2\right)}{N_{p}}\right]\,.$ (S82)

Here

$\tilde{\mu}_{m}=g_{m}\mu_{m}\,,$ (S83)

with $g_{m}$ being the Jackson kernel in Eq. (S6). The integrals involving
$f\left(x\right)$ can be obtained via Chebyshev-Gauss quadrature,

$\displaystyle\int_{-1}^{1}f\left(x\right)g\left(x\right)dx$
$\displaystyle=\int_{-1}^{1}dx\frac{\sqrt{1-x^{2}}f\left(x\right)g\left(x\right)}{\sqrt{1-x^{2}}}=\frac{1}{N_{p}}\sum_{k=0}^{N_{p}-1}\gamma_{k}g\left(x_{k}\right)\,.$
(S84)

Thus we have

$\displaystyle\Delta_{i}$
$\displaystyle=\int_{-E_{\mathsf{max}}}^{E_{\mathsf{max}}}\Delta_{i}\left(E\right)f\left(E\right)dE=\int_{-1}^{1}\tilde{\Delta}_{i}\left(\tilde{E}\right)f\left(\tilde{E}\right)d\tilde{E}=\frac{1}{N_{p}}\sum_{k=0}^{N_{p}-1}\gamma_{k}^{\Delta_{i}}f\left(E_{k}\right)\,,$
(S85) $\displaystyle n_{i}$
$\displaystyle=\int_{-E_{\mathsf{max}}}^{E_{\mathsf{max}}}\rho_{i}\left(E\right)f\left(E\right)dE=\int_{-1}^{1}\tilde{\rho}_{i}\left(\tilde{E}\right)f\left(\tilde{E}\right)d\tilde{E}=\frac{1}{N_{p}}\sum_{k=0}^{N_{p}-1}\gamma_{k}^{\rho_{i}}f\left(E_{k}\right)\,,$
(S86)

Here $f\left(E_{k}\right)$ is the Fermi distribution function.

The self-consistent BdG theory can be performed until the difference of the
output order parameters $\left\\{\Delta_{i},n_{i}\right\\}$ are close enough
to the input ones. At each iteration, $\Delta_{i}$ and $n_{i}$ need to be
evaluated site by site. The convergence of the iteration can be further
accelerated by Broyden’s method [4, 5]. The KPM-based method requires only a
small amount of memory and computation effort that scales with
$N^{2}N_{\mathcal{C}}$, with $N$ being the number of lattice sites. On the
other hand, the computational effort of ED scales with $N^{3}$. On our test
platform with $2$ CPUs [Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz], the KPM-
based method surpasses the ED-based one when $N\approx 2000$ and
$N_{\mathcal{C}}=4096$. Therefore, the KPM-based mean-field method is superior
to ED-based one as long as we need to deal with systems with more than several
thousand lattice sites. Moreover, the KPM-based method can be easily
accelerated with a parallel computation scheme like MPI. The communication
cost is very low and the acceleration scales linearly with the number of MPI
processes. Thus, very large system sizes can be accessed with enough
resources.

####  I.4.3 Continuum model with local interaction

For the qBM-BCS model, the Hamiltonian can be written as

$H=\int
d\mathbf{r}\,\chi^{\dagger}\left[h\left(\mathbf{r}\right)\mu^{3}+\Delta\left(\mathbf{r}\right)\mu^{1}\right]\chi\,,$
(S87)

where $h\left(\mathbf{r}\right)$ is the qBM one-body Hamiltonian from Eq. (1)
in the main text. Both terms are local in real space, but only the second term
can be expressed as a sparse matrix in position space. The first term involves
the Dirac operator $-i\bm{\sigma}\cdot\nabla$, which cannot be expressed as a
short-range lattice Hamiltonian without introducing the fermion-doubling
problem. Therefore, we work simultaneously with both real- and momentum-space
lattices, switching back and forth as necessary. Starting with a real-space
state $\ket{\psi}$, we have

$H\ket{\psi}=H_{2}\ket{\psi}+\mathcal{F}^{-1}\left(H_{1}\ket{\mathcal{F}\psi}\right)\,.$
(S88)

Here $\mathcal{F}$ is the Fourier transform from real space to momentum space,
and

$\displaystyle H_{1}$ $\displaystyle=\int
d\mathbf{r}\,\chi^{\dagger}h\left(\mathbf{r}\right)\mu^{3}\chi=\intop_{\mathbf{k}_{1}\mathbf{k}_{2}}\chi_{\mathbf{k}_{1}}^{\dagger}\hat{h}_{\mathbf{k}_{1}\mathbf{k}_{2}}\chi_{\mathbf{k}_{2}}\,,$
(S89) $\displaystyle H_{2}$ $\displaystyle=\intop
d\mathbf{r}\,\chi^{\dagger}\Delta\left(\mathbf{r}\right)\mu^{1}\chi\,.$ (S90)

The Hamiltonian matrix of the first term is sparse in momentum space and the
action of $H_{1}$ on the momentum space wave function $\ket{\mathcal{F}\psi}$
is cheap to calculate. On the other hand, the Hamiltonian matrix of the second
term is sparse in real space, as it acts locally on the real-space wave
function. With this procedure, we can deal with a system with Dirac cones
exactly without involving the fermion-doubling problem. The most computational
challenging part now lies on the Fourier transform of wave functions, which
can be done by fast-Fourier algorithm provided in vendor software packages.

###  I.5 Calculation of superfluid stiffness with KPM

####  I.5.1 Superfluid stiffness of lattice model

Consider the monochromatic electromagnetic field with vector potential
$\mathbf{A}\left(\mathbf{r}_{i},t\right)$,

$\mathbf{A}\left(\mathbf{r}_{i},t\right)=\Re\left[\mathbf{A}\left(\omega,\mathbf{q}\right)e^{-i\omega
t+i\mathbf{q}\cdot\mathbf{r}_{i}}\right]\,,\qquad\mathbf{E}=-\partial_{t}\mathbf{A}\,.$
(S91)

With Peierls substitution, the hopping terms are modified as

$t_{ij}\to
t_{ij}e^{ie\phi_{ij}}\,,\qquad\phi_{ij}=A_{x}\left(\mathbf{r}_{i},t\right)\left(x_{i}-x_{j}\right)\,.$
(S92)

Here we assume the slowly-varying electric field is applied in $x$-direction
and the hopping is short-ranged. Then the attractive Hubbard model under
vector field is given by

$H=H_{0}-e\sum_{i}j_{x}^{P}\left(\mathbf{r}_{i}\right)A_{x}\left(\mathbf{r}_{i}\right)-\frac{e^{2}}{2}\sum_{i}A_{x}^{2}\left(\mathbf{r}_{i}^{2}\right)K_{x}\left(\mathbf{r}_{i}\right)\,.$
(S93)

The paramagnetic current and kinetic-energy densities are

$\displaystyle j_{x}^{P}\left(\mathbf{r}_{i}\right)$
$\displaystyle=i\sum_{j\sigma}t_{ij}\left(x_{i}-x_{j}\right)c_{i\sigma}^{\dagger}c_{j\sigma}\,,$
(S94) $\displaystyle K_{x}\left(\mathbf{r}_{i}\right)$
$\displaystyle=-\sum_{j\sigma}t_{ij}\left(x_{i}-x_{j}\right)^{2}c_{i\sigma}^{\dagger}c_{j\sigma}\,.$
(S95)

The total current is given by

$j_{x}\left(\mathbf{r}_{i}\right)=-\frac{\delta H}{\delta
A_{x}\left(\mathbf{r}_{i}\right)}=eJ_{x}^{P}\left(\mathbf{r}_{i}\right)+e^{2}A_{x}\left(\mathbf{r}_{i}\right)K_{x}\left(\mathbf{r}_{i}\right)\,.$
(S96)

With linear response theory, we have

$\displaystyle J_{x}^{P}\left(\omega,\mathbf{q}\right)$
$\displaystyle=-e\Pi_{xx}^{R}\left(\omega,\mathbf{q}\right)A_{x}\left(\omega,\mathbf{q}\right)\,,$
(S97) $\displaystyle j_{x}\left(\omega,\mathbf{q}\right)$
$\displaystyle=e^{2}\left[-\Pi_{xx}^{R}\left(\omega,\mathbf{q}\right)+K_{x}\right]A_{x}\left(\omega,\mathbf{q}\right)\,.$
(S98)

The superfluid stiffness is then given by [6]

$\frac{D_{s}}{\pi e^{2}}=\Pi_{xx}^{R}\left(\omega=0,q_{x}=0,q_{y}\to
0\right)-\left\langle K_{x}\right\rangle\,.$ (S99)

It can be shown that the superfluid stiffness for the lattice model can be
evaluated via

$\displaystyle\frac{D_{s}}{\pi e^{2}}$
$\displaystyle=-\Pi_{xx}^{R}\left(\omega=0,q_{x}=0,q_{y}\to
0\right)-\left\langle K_{x}\right\rangle\,,$ (S100a)
$\displaystyle\Pi_{xx}^{R}\left(\omega,\mathbf{q}\right)$
$\displaystyle=\frac{2}{V}\sum_{nm}\frac{f\left(E_{n}\right)-f\left(E_{m}\right)}{\omega+E_{n}-E_{m}+i\eta}\bra{n}J_{e}^{x}\left(\mathbf{q}\right)\ket{m}\bra{m}J^{x}\left(-\mathbf{q}\right)\ket{n}\,,$
(S100b) $\displaystyle\hat{J}_{e}^{x}\left(\mathbf{q}\right)$
$\displaystyle=i\sum_{ij}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}t_{ij}\left(x_{i}-x_{j}\right)\ket{i\uparrow}\bra{j\uparrow}\,,$
(S100c) $\displaystyle\hat{J}^{x}\left(\mathbf{q}\right)$
$\displaystyle=i\sum_{ij}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}t_{ij}\left(x_{i}-x_{j}\right)\ket{i\uparrow}\bra{j\uparrow}+i\sum_{ij}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}t_{ij}\left(x_{i}-x_{j}\right)\ket{i\downarrow}\bra{j\downarrow}\,,$
(S100d) $\displaystyle\left\langle K_{x}\right\rangle$
$\displaystyle=-\frac{2}{V}\sum_{ij}t_{ij}\left(x_{i}-x_{j}\right)^{2}\sum_{n}\left(u_{n}^{i}\right)^{*}u_{n}^{j}f\left(E_{n}\right)=\sum_{n}\bra{n}\hat{K}_{x}f\left(E_{n}\right)\ket{n}\,,$
(S100e)
$\displaystyle\hat{K}_{x}=-\frac{2}{V}\sum_{ij}t_{ij}\left(x_{i}-x_{j}\right)^{2}\ket{i\uparrow}\bra{j\uparrow}\,.$
(S100f)

####  I.5.2 KPM implementation

We define the two-point spectral function,

$\pi\left(\omega_{1},\omega_{2},\mathbf{q}\right)=\frac{2}{V}\sum_{nm}\bra{n}\hat{J}_{e}^{x}\left(\mathbf{q}\right)\ket{m}\bra{m}\hat{J}^{x}\left(-\mathbf{q}\right)\ket{n}\delta\left(\omega_{1}-E_{n}\right)\delta\left(\omega_{2}-E_{m}\right)\,.$
(S101)

Then the rescaled spectral function can be expressed using the kernel
polynomial method,

$\tilde{\pi}\left(\omega_{1},\omega_{2}\right)=\sum_{kl=0}^{N_{\mathcal{C}}-1}\frac{\mu_{kl}h_{kl}g_{k}g_{l}T_{k}\left(\omega_{1}\right)T_{l}\left(\omega_{2}\right)}{\pi^{2}\sqrt{\left(1-\omega_{1}^{2}\right)\left(1-\omega_{2}^{2}\right)}}\,.$
(S102)

Here $\tilde{\pi}$ stands for the function with argument rescaled to the
interval $\left[-1,1\right]$ and

$h_{kl}=\frac{2}{1+\delta_{k,0}}\frac{2}{1+\delta_{l,0}}\,,$ (S103)

and $g_{k}$ are the kernel damping factors. We use the Jackson kernel for the
spectral function.

The moments are obtained with

$\displaystyle\mu_{kl}$
$\displaystyle=\int_{-1}^{1}d\omega_{1}\int_{-1}^{1}d\omega_{2}\tilde{\pi}_{xx}\left(\omega_{1},\omega_{2}\right)T_{k}\left(\omega_{1}\right)T_{l}\left(\omega_{2}\right)$
$\displaystyle=\frac{2}{N}\sum_{nm}\bra{n}\hat{J}_{e}^{x}\left(\mathbf{q}\right)\ket{m}\bra{m}\hat{J}^{x}\left(-\mathbf{q}\right)\ket{n}T_{k}\left(\tilde{E}_{n}\right)T_{l}\left(\tilde{E}_{m}\right)$
$\displaystyle=\frac{2}{V}\sum_{nm}\bra{n}T_{k}\left(\tilde{H}\right)\hat{J}_{e}^{x}\left(\mathbf{q}\right)\ket{m}\bra{m}T_{l}\left(\tilde{H}\right)\hat{J}^{x}\left(-\mathbf{q}\right)\ket{n}$
$\displaystyle=\frac{2}{V}\mathsf{Tr}\left[T_{k}\left(\tilde{H}\right)\hat{J}_{e}^{x}\left(\mathbf{q}\right)T_{l}\left(\tilde{H}\right)\hat{J}^{x}\left(-\mathbf{q}\right)\right]\,.$
(S104)

The trace can be efficiently and accurately averaged over a small number of
random vectors.

In practice, the evaluation of $\tilde{\Pi}_{xx}$ in Eq. (S102) can be
efficiently implemented with the help of the discrete FFT, similar to the one-
point case. We would evaluate
$\tilde{\Pi}_{xx}\left(\omega_{1},\omega_{2}\right)$ at discrete points,

$\omega_{1}=\omega_{m}=\cos\left[\frac{\pi\left(m+1/2\right)}{N_{p}}\right]\,,\qquad\omega_{2}=\omega_{n}=\cos\left[\frac{\pi\left(n+1/2\right)}{N_{p}}\right]\,.$
(S105)

Here $N_{p}$ is the number of points evaluated for the integral and a
reasonable choice is $N_{p}=2N_{\mathcal{C}}$. Then we have

$\gamma_{mn}=\pi^{2}\sqrt{\left(1-\omega_{m}^{2}\right)\left(1-\omega_{n}^{2}\right)}\,\tilde{\pi}\left(\omega_{m},\omega_{n}\right)=\sum_{kl=0}^{N_{\mathcal{C}}-1}\mu_{kl}h_{kl}g_{k}g_{l}\cos\left[\frac{\pi
k\left(m+\frac{1}{2}\right)}{N_{p}}\right]\cos\left[\frac{\pi
l\left(n+\frac{1}{2}\right)}{N_{p}}\right]\,.$ (S106)

The summation is actually the discrete Fourier transform of the 2D array
$\tilde{\mu}_{kl}=\mu_{kl}g_{k}g_{l}$, and it can be efficiently evaluated via
the fast-Fourier transform. Thus we have

$\displaystyle\Pi_{xx}\left(\omega,\mathbf{q}\right)$
$\displaystyle=\int_{-1}^{1}d\omega_{1}\int_{-1}^{1}d\omega_{2}\frac{f\left(\omega_{1}\right)-f\left(\omega_{2}\right)}{\omega+\omega_{1}-\omega_{2}+i\eta}\tilde{\pi}\left(\omega_{1},\omega_{2}\right)$
$\displaystyle=\frac{1}{N_{p}^{2}}\sum_{mn=0}^{N_{p}-1}\gamma_{mn}\frac{f\left(\omega_{n}\right)-f\left(\omega_{m}\right)}{\omega+\omega_{n}-\omega_{m}+i\eta}\,.$
(S107)

Here the integral is evaluated with the Gaussian quadrature.

For the kinetic energy term, we can defined the kinetic spectral function,

$k_{x}\left(\omega\right)=\sum_{n}\bra{n}\hat{K}_{x}\ket{n}\delta\left(\omega-
E_{n}\right)\,.$ (S108)

This term can also be evaluated with the KPM in the same fashion. At the end
of the day, we would have

$\left\langle
K_{x}\right\rangle=\frac{1}{N_{p}}\sum_{m=0}^{N_{p}-1}\gamma_{m}f\left(\omega_{m}\right)\,.$
(S109)

In Fig. S2, we show the superfluid stiffness of a Dirac superconductor on the
honeycomb lattice with onsite disorder potential evaluated with the KPM. The
superfluid stiffness $D_{s}/\pi$ follows the analytical prediction
$D_{s}/\pi=2\Delta/\pi$ when the pairing amplitude is weak in the clean
system. Interestingly, the superfluid stiffness is robust against a weak
disorder potential and is only weakened by strong disorder.

Figure S2: Superfluid stiffness of a spin-singlet $s$-wave superconductor on
the honeycomb lattice with uniform distributed on-site disorder potential
$W_{i}\in\left[-\frac{W}{2},\frac{W}{2}\right]$. The numerical data are
evaluated on a honeycomb lattice with $100\times 100$ unit cells and
$N_{\mathcal{C}}=1024$. The dashed line shows the analytical prediction
$\frac{D_{s}}{\pi}=\frac{2}{\pi}\Delta$ for a clean 2D Dirac superconductor.

####  I.5.3 Superfluid stiffness for Dirac superconductor

Now we consider Dirac superconductor with spin-singlet, layer-singlet,
pseudospin-singlet pairing,

$H=\sum_{\mathbf{r}}c_{\mathbf{r}}^{\dagger}\left(-iv_{F}\bm{\sigma}\cdot\nabla\right)c_{\mathbf{r}}+\Delta\sum_{\mathbf{r}}\left[ic_{\mathbf{r}\uparrow}^{\dagger}\sigma^{2}\kappa^{2}\left(c_{\mathbf{r}\downarrow}^{\dagger}\right)^{\mathsf{T}}+\textrm{H.c.}\right]\,.$
(S110)

Here $c_{\mathbf{r}}$ carries pseudospin ($\sigma$), layer ($\tau$), and
physical spin-1/2 ($s$) indices. Introducing the Nambu spinor

$\chi=\left[\begin{array}[]{c}\psi_{\uparrow}\\\
i\sigma^{2}\kappa^{2}\left(\psi_{\downarrow}^{\dagger}\right)^{\mathsf{T}}\end{array}\right]_{\mu}\,,$
(S111)

we have

$H=\sum_{\mathbf{r}}\chi^{\dagger}\,\hat{h}_{\mathsf{BdG}}\,\chi\,,\qquad\hat{h}_{\mathsf{BdG}}=\left(-iv_{F}\bm{\sigma}\cdot\nabla\right)\mu^{3}+\Delta\mu^{1}\,.$
(S112)

Here $\mu^{i}$ is the Pauli matrix in the particle-hole sector. The electric
current operator is then given by

$\mathbf{J}=\chi^{\dagger}\bm{\sigma}\chi\,.$ (S113)

The current-current correlation function is then given by

$\Pi_{xx}^{R}\left(\omega,\mathbf{q}\right)=\frac{1}{V}\sum_{nm}\frac{f\left(E_{n}\right)-f\left(E_{m}\right)}{\omega+E_{n}-E_{m}+i\eta}\bra{n}J_{x}\left(\mathbf{q}\right)\ket{m}\bra{m}J_{x}\left(-\mathbf{q}\right)\ket{n}\,,$
(S114)

which can be evaluated with KPM in the same fashion as the lattice model.
However, unlike the lattice model, the diamagnetic term $\left\langle
K_{x}\right\rangle$ vanishes in the expression for the superfluid stiffness
and the paramagnetic current-current correlation functions diverges if no
momentum cutoff is enforced. The superfluid stiffness should vanish when
$\Delta=0$, i.e. there is no Meissner effect for normal state. Therefore, we
can regularize the current-current correlation function to obtain the physical
superfluid stiffness via

$\frac{D_{s}}{\pi}=-\Pi_{xx}^{R}\left(\omega=0,q_{x}=0,q_{y}\to 0,\Delta\neq
0\right)+\Pi_{xx}^{R}\left(\omega=0,q_{x}=0,q_{y}\to 0,\Delta=0\right)\,.$
(S115)

##  II More Details on DOS and critical states

###  II.1 DOS and multifractal properties in the lakes between filaments

Figure S3: The zero energy DOS $\nu(E=0)$, $N_{\mathcal{C}}$ scaling of
$\nu(E=0)$ and multifractal dimension of zero energy wave functions in the
lake of DOS map. In (b) and (d), we have $M=233$ and $N_{\Lambda}=15$ and the
$N_{\mathcal{C}}$ dependence are fitted with $\sim N_{\mathcal{C}}^{-0.352}$
and $\sim N_{\mathcal{C}}^{-0.209}$, respectively.

The density of states of a system is given by

$\rho\left(E\right)=\frac{1}{L}\sum_{i}\delta\left(E-\varepsilon_{i}\right)\,,$
(S116)

with $\\{\varepsilon_{i}\\}$ denoting the set of eigenenergies. In the KPM
[1], the Dirac delta functions are approximated by smooth functions. Using the
Jackson kernel, these are Gaussians

$\delta_{\text{KPM}}^{J}\left(x\right)\approx\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left[-\frac{\left(x-a\right)^{2}}{2\sigma^{2}}\right]\,.$
(S117)

The variance

$\sigma^{2}=\left(\frac{\pi}{N_{\mathcal{C}}}\right)^{2}\left[1-a^{2}+\frac{3a^{2}-2}{N_{\mathcal{C}}}\right]\,.$
(S118)

Here $-1\leq a\leq 1$ and $\sigma=\pi/N_{\mathcal{C}}$ at $a=0$,
$\sigma=\pi/N_{\mathcal{C}}^{3/2}$ at $a=\pm 1$, with $N_{\mathcal{C}}$ being
the order of the Chebyshev expansion.

Suppose we have a system whose density of states has power-law dependence on
energy near $E=0$,

$\rho\left(E\right)\propto E^{\alpha}\,.$ (S119)

Then in the KPM with Jackson kernel, we have

$\displaystyle\rho_{KPM}\left(E=0\right)$
$\displaystyle\approx\frac{1}{\sqrt{2\pi\sigma^{2}}}\int
dE\,\rho\left(E\right)e^{-\frac{E^{2}}{2\sigma^{2}}}=\frac{1}{\sqrt{2\pi\sigma^{2}}}\int
dE\,E^{\alpha}e^{-\frac{E^{2}}{2\sigma^{2}}}=\frac{\sigma^{\alpha+1}}{\sqrt{2\pi\sigma^{2}}}\int
dx\,x^{\alpha}e^{-\frac{x^{2}}{2}}$ $\displaystyle\propto\sigma^{\alpha}\sim
N_{\mathcal{C}}^{-\alpha}\,.$ (S120)

The scaling behavior holds when $N_{\mathcal{C}}$ is not too small or too
large. When $N_{\mathcal{C}}$ is too small, the higher-order corrections are
significant and spoil the simple power-law scaling. When $N_{\mathcal{C}}$ is
large enough that $\sigma$ is at the order of level spacing of the system, the
DOS resolves into isolated Gaussian peaks.

Fig. S3 shows the $N_{\mathcal{C}}$ scaling of the zero-energy DOS and
multifractal dimensions of the zero-energy wave functions at two points
$A_{1}=0.7$, $A_{2}=0.75$ and $A_{1}=1$, $A_{2}=0.6$ in the lakes between the
filaments. When $N_{\mathcal{C}}$ is not very large, the $\nu(E=0)$ has power-
law dependence on $N_{\mathcal{C}}$. It implies that the DOS near zero energy
has power-law dependence on energy with $\nu(E)\sim\left|E\right|^{0.352}$ for
$A_{1}=0.7$, $A_{2}=0.75$, and $\nu(E)\sim\left|E\right|^{0.209}$ for
$A_{1}=1$, $A_{2}=0.6$, respectively.

Figure S4: The third multifractal dimension $\tau_{3}$ averaged for states
near zero energy (a)–(d) or finite energies (e), (f), along with the zero-
energy DOS $\nu(E=0)$ for the golden-ratio qBM model. (a) and (b): along the
vertical cuts at $A_{1}=0.4$ and $1$, respectively; (c) and (d) along the
horizontal cuts at $A_{2}=0.2$ and $0.5$, respectively; (e) and (f) $\tau_{3}$
versus energy for $A_{1}=0.4$ and $A_{2}=0.36$ (at a filament), $A_{2}=0.6$
(in the lake between filaments). The navy blue dashed and orange dot-dashed
lines indicate $\tau_{3}=5/2$ for WZW states and $13/4$ for SQHT states,
respectively [7].

###  II.2 The third multifractal dimension $\tau_{3}$

In Figs. 3 and 4 of the main text, we show the second multifractal dimension
$\tau_{2}$ for the qBM model. Fig. S4 shows the third multifractal dimension
$\tau_{3}$, which is similarly defined via Eqs. (5) and (6) of the main text.

##  III Silver Ratio QBM Model

Figure S5: Zero-energy DOS and the second multifractal dimension for the
silver-ratio qBM model with $\mathbf{q}_{i}=\beta\,\mathbf{q}_{i}^{\prime}$,
with $\beta$ a commensurate approximant to the inverse silver ratio
$1/\left(\sqrt{2}+1\right)$.

The critical filaments and multifractal properties of wave functions are
universal for the quasiperiodic Bistritzer-MacDonald model. In the main text,
we focus on the case $\mathbf{q}_{i}=\beta\,\mathbf{q}^{\prime}_{i}$ with
$\beta$ the inverse golden ratio. In this section, we consider instead
$\beta=1/(\sqrt{2}+1)$, which is the inverse silver ratio. The silver ratio
can be approximately by $P_{n+1}/P_{n}$ with $P_{n}=1,2,5,12,29,70,\dots$ the
Pell numbers. Fig. S5 shows the zero-energy DOS map, and the second
multifractal dimension $\tau_{2}$ along two horizontal cuts $A_{2}=0.5$ and
$A_{2}=1$ for the silver-ratio qBM model. There are filaments with large-zero
energy DOS connecting the magic angle points of $A_{1}$ and $A_{2}$. Between
the filaments, the nonzero zero-energy DOS shows that the system is no longer
semimetallic. At the filaments and in the lakes between the filaments, the
wave functions near the zero energy become multifractal. The $\tau_{2}$ is
close to that of WZW state ($\tau_{2}=3/2$) for zero-energy wave functions in
the lake along the cut $A_{1}=0.5$. In the lake along the $A_{2}=1$ cut,
$\tau_{2}$ is between that of WZW and SQHT states ($\tau_{2}=7/4$). At the
filaments, $\tau_{2}$ is smaller than $3/2$, indicating that these states are
strongly multifractal. The results are very similar to the case of golden-
ratio qBM model, indicating that our conclusions in the main text are general.

##  IV More details on the superconductivity

###  IV.1 Numerical results

We consider superconductivity in the qBM model with a local pairing
interaction,

$H=\int
d\mathbf{r}\,\psi^{\dagger}(\mathbf{r})\,h(\mathbf{r})\,\psi(\mathbf{r})-U\int
d\mathbf{r}\,\mathcal{B}^{\dagger}\,\mathcal{B}\,.$ (S121)

Here $h(\mathbf{r})$ is the qBM Hamiltonian given by Eq. (1) of the main text,
$\psi\rightarrow\psi_{\sigma,\tau,s}$ carries indices in sublattice
($\sigma$), layer ($\tau$), and physical spin-1/2 ($s$) spaces, $U$ gives the
strength of the local attractive interaction, and $\mathcal{B}$ is the
bilinear for pseudospin-singlet, layer-singlet, and spin-singlet $s$-wave
pairing

$\mathcal{B}=\psi^{\mathsf{T}}\sigma^{2}\tau^{2}s^{2}\psi.$ (S122)

The local interaction is well-known to be RG irrelevant for 2+1-D Dirac
fermions and superconductivity only arises for an interaction strength
comparable to the ultraviolet energy cutoff, see Eq. (S138) below. In the
presence of the usual Bistritzer-MacDonald (monochromatic) moiré potential,
the low-energy spectrum remains linear and strong interactions are required
for superconductivity, except in the close vicinity of the magic angle where a
nearly flat band results in a much smaller critical interaction strength.

The attractive interaction can be decoupled with the mean-field ansatz,

$\Delta(\mathbf{r})=-U\left<\mathcal{B}(\mathbf{r})\right>.$ (S123)

Because we work with the quasiperiodic variant of the single-valley BM model,
the $s$-wave pairing considered here would correspond to an FFLO state
relative to the microscopic honeycomb lattices. We expect that the physics of
pairing enhancement due to Chalker scaling of quantum-critical wave functions
[8, 9] and (approximate) topology-protected criticality will hold for
intervalley pairing as well, but we leave this as a study for future work.

The mean-field equations [Eq. (8) in the main text and Eq. (S123)] are then
readily solved self-consistently with the KPM based method. The KPM-based
method allows us to deal with system sizes inaccessible with exact
diagonalization. In Fig. 5 of the main text, the data are obtained with $M=8$
and $N_{\Lambda}=11$, and the dimension of the Hamiltonian matrix is $61952$,
which is very hard to diagonalize exactly even with the help of modern GPUs.

Fig. S6 shows the dependence of spatially averaged pairing amplitude
$\left<\Delta\right>$ on the system size and the degree of the kernel
polynomials. We show the data with three system sizes $M=5$, $8$, and $M=13$
with $N_{\Lambda}=11$, which is solved on a triangular lattice with $55\times
55$, $88\times 88$ and $143\times 143$ sites, respectively. The
$\left<\Delta\right>$ only slightly changes with $M$. It indicates that the
system size employed is large enough and the results apply for systems in
thermodynamic limit. Fig. S6(b) show the $N_{\mathcal{C}}$ dependence of
$\left<\Delta\right>$ and the results converges when $N_{\mathcal{C}}$ is
larger than $10^{3}$. In Fig. 5 of the main text, we use $M=8$,
$N_{\Lambda}=11$ and $N_{\mathcal{C}}=2048$.

Figure S6: The average pairing amplitude for qBM-BCS model with different
system size $M$ and degree of kernel polynomials $N_{\mathcal{C}}$.

###  IV.2 Pairing of exact eigenstates approximation

Consider Dirac fermions with local attractive interactions and a spatially
inhomogeneous potential,

$H=\sum_{\mathbf{r}}c_{\mathbf{r}}^{\dagger}\,h(\bm{\mathrm{r}})\,c_{\mathbf{r}}-U\sum_{\mathbf{r}}\left[c_{\mathbf{r}\uparrow}^{\dagger}\left(-i\sigma^{2}\tau^{2}\right)\left(c_{\mathbf{r}\downarrow}^{\dagger}\right)^{\mathsf{T}}\right]\left[\left(c_{\mathbf{r}\downarrow}\right)^{\mathsf{T}}i\sigma^{2}\tau^{2}c_{\mathbf{r}\uparrow}\right]\,.$
(S124)

Here $c_{\mathbf{r}}$ carries pseudospin ($\sigma$), layer ($\tau$) and spin
($s$) indices, and the one-body Dirac Hamiltonian $h(\bm{\mathrm{r}})$
incorporates random or structured matrix potentials in some symmetry class,
e.g. the class-CI chiral qBM model in Eq. (1) of the main text. As in that
case, we assume spin SU(2) and time-reversal symmetries in the normal state.
Suppose the non-interacting part is solved by

$h\ket{\alpha}=\varepsilon_{\alpha}\ket{\alpha}\,.$ (S125)

Then we have

$c_{\mathbf{r}\sigma\tau
s}=\sum_{\alpha}\braket{\mathbf{r\sigma}\tau}{\alpha}a_{\alpha
s}=\sum_{\alpha}\psi_{\alpha}\left(\mathbf{r}\sigma\tau\right)a_{\alpha s}\,.$
(S126)

Expressing the Eq. (S124) in this basis gives

$H=\sum_{\alpha,s}\varepsilon_{\alpha}a_{\alpha s}^{\dagger}a_{\alpha
s}-U\sum_{\mathbf{r}}\sum_{\alpha\beta\gamma\delta}\left[\psi_{\alpha}^{\dagger}\left(\mathbf{r}\right)\left(-i\sigma^{2}\tau^{2}\right)\left(\psi_{\beta}^{{}^{\dagger}}\left(\mathbf{r}\right)\right)^{\mathsf{T}}\psi_{\gamma}^{\mathsf{T}}\left(\mathbf{r}\right)\left(i\sigma^{2}\tau^{2}\right)\psi_{\delta}\left(\mathbf{r}\right)\right]a_{\alpha\uparrow}^{\dagger}a_{\beta\downarrow}^{\dagger}a_{\gamma\downarrow}a_{\delta\uparrow}\,.$
(S127)

The Cooper pairs are formed by pairing the time-reversal counterparts of spin-
up and -down electrons. The time-reversed conjugate of $\ket{\alpha}$ is given
by

$\ket{\bar{\alpha}}=-i\sigma^{2}\tau^{2}\ket{\alpha^{*}}\,.$ (S128)

We neglect interactions in Eq. (S127) except those responsible for creating
such pairs. Thus we obtain the reduced BCS Hamiltonian

$\displaystyle H_{\mathsf{red}}$
$\displaystyle=\sum_{\alpha,s}\varepsilon_{\alpha}a_{\alpha
s}^{\dagger}a_{\alpha
s}-U\sum_{\alpha\beta}\left[\sum_{\mathbf{r}}\psi_{\alpha}^{*}\left(\mathbf{r}\right)\left(-i\sigma^{2}\tau^{2}\right)\psi_{\bar{\alpha}}^{*}\left(\mathbf{r}\right)\psi_{\bar{\beta}}\left(\mathbf{r}\right)\left(i\sigma^{2}\tau^{2}\right)\psi_{\beta}\left(\mathbf{r}\right)\right]a_{\alpha\uparrow}^{\dagger}a_{\bar{\alpha}\downarrow}^{\dagger}a_{\bar{\beta}\downarrow}a_{\beta\uparrow}$
$\displaystyle=\sum_{\alpha,s}\varepsilon_{\alpha}a_{\alpha
s}^{\dagger}a_{\alpha
s}-\sum_{\alpha\beta}M_{\alpha\beta}\,a_{\alpha\uparrow}^{\dagger}a_{\bar{\alpha}\downarrow}^{\dagger}a_{\bar{\beta}\downarrow}a_{\beta\uparrow}\,.$
(S129)

Here

$M_{\alpha\beta}=U\sum_{\mathbf{r}}\left|\psi_{\alpha}\left(\mathbf{r}\right)\right|^{2}\left|\psi_{\beta}\left(\mathbf{r}\right)\right|^{2}\,,$
(S130)

where we have used

$\psi_{\alpha}=-i\sigma^{2}\tau^{2}\psi_{\alpha}^{*}\,,\qquad\left|\psi_{\alpha}\left(\mathbf{r}\right)\right|^{2}\equiv\sum_{\sigma\tau}\left|\psi_{\alpha}\left(\mathbf{r}\sigma\tau\right)\right|^{2}\,.$
(S131)

Then we have the mean-field gap equation,

$\Delta_{\alpha}=\sum_{\lambda}M_{\alpha\lambda}\tanh\left(\frac{\beta\sqrt{\varepsilon_{\lambda}^{2}+\Delta_{\lambda}^{2}}}{2}\right)\frac{\Delta_{\lambda}}{2\sqrt{\varepsilon_{\lambda}^{2}+\Delta_{\lambda}^{2}}}\,.$
(S132)

In the _clean limit_ , the wave functions are plane waves and the spatial
overlap of wave functions at different energies is the same,

$M_{\alpha\beta}\approx U\left(\frac{a}{L}\right)^{d}\,.$ (S133)

Here $a$ is the lattice constant and $L$ is the system size. The mean field
gap equation is reduced to

$\Delta_{\alpha}=U\sum_{\lambda}\left(\frac{a}{L}\right)^{d}\tanh\left(\frac{\beta\sqrt{\varepsilon_{\lambda}^{2}+\Delta_{\lambda}^{2}}}{2}\right)\frac{\Delta_{\lambda}}{2\sqrt{\varepsilon_{\lambda}^{2}+\Delta_{\lambda}^{2}}}\,.$
(S134)

Assuming that the pairing amplitude is energy independent,
$\Delta\left(\varepsilon_{\alpha}\right)=\Delta$, we have

$1\approx
Ua^{d}\int_{-\Lambda}^{\Lambda}d\varepsilon\,\nu\left(\varepsilon\right)\tanh\left(\frac{\beta\sqrt{\varepsilon^{2}+\Delta^{2}}}{2}\right)\frac{1}{2\sqrt{\varepsilon^{2}+\Delta^{2}}}\,.$
(S135)

For Dirac fermions in 2D, the DOS is given by

$\nu(\varepsilon)=\frac{1}{\pi v_{F}^{2}}\left|\varepsilon\right|\,.$ (S136)

Thus, we have the zero-temperature gap equation

$\frac{2}{U}=a^{2}\int_{-\Lambda}^{\Lambda}d\varepsilon\frac{\left|\varepsilon\right|}{\pi
v_{F}^{2}}\frac{1}{\sqrt{\varepsilon^{2}+\Delta^{2}}}=\frac{a^{2}}{\pi
v_{F}^{2}}\left(\sqrt{\Lambda^{2}+\Delta^{2}}-\Delta\right)\,.$ (S137)

Here $\Lambda$ is the energy cutoff of the Dirac fermions. The equation has a
solution for nonzero $\Delta$ only when the interaction strength is above the
critical value

$U_{c}=\frac{2\pi v_{F}^{2}}{\Lambda a^{2}}=\frac{2\Lambda}{\pi}\,,$ (S138)

where we have assumed that $\Lambda=v_{F}(\pi/a)$. Thus we obtain the well-
known result that superconductivity only arises above the critical interaction
strength for clean 2D Dirac fermions [10, 11, 12, 13, 14].

For _critical states_ , the spatial overlap of the wave functions entering the
matrix element is determined by Chalker scaling [15, 16, 17], and we have

$M_{\alpha\beta}\sim
U\mathfrak{b}\left|\varepsilon_{\alpha}-\varepsilon_{\beta}\right|^{(\tau_{2}-d)/z}\,.$
(S139)

Here $z$ is the dynamical critical exponent of Dirac fermions with random or
structural spatial inhomogeneity. The DOS of critical Dirac fermions near
$\varepsilon=0$ is modified by the disorder or quasiperiodic potential,

$\nu(\varepsilon)\approx\mathfrak{c}\left|\varepsilon\right|^{(d-z)/z}\,.$
(S140)

Thus we have the gap equation at $T=0$,

$1=U\int_{-\Lambda}^{\Lambda}d\varepsilon\,\nu(\varepsilon)\,M_{\alpha\beta}\tanh\left(\frac{\beta\sqrt{\varepsilon^{2}+\Delta^{2}}}{2}\right)\frac{1}{2\sqrt{\varepsilon^{2}+\Delta^{2}}}\approx
U\mathfrak{a}\int_{-\Lambda}^{\Lambda}d\varepsilon\left|\varepsilon\right|^{(\tau_{2}-z)/z}\tanh\left(\frac{\beta\sqrt{\varepsilon^{2}+\Delta^{2}}}{2}\right)\frac{1}{2\sqrt{\varepsilon^{2}+\Delta^{2}}}.$
(S141)

Here $\mathfrak{a}=\mathfrak{b}\mathfrak{c}$. For 2D disordered Dirac fermions
in the chiral/WZNW version of class CI with winding number $2$, we have [15]

$z=\frac{7}{4}$ (S142)

If we further assume that, due to SWQC [7, 17], the appropriate wave-function
fractal dimension is that of the class-C SQHT

$\displaystyle d_{2}\approx 7/4,$ (S143)

then we have

$1\approx
U\mathfrak{a}\int_{-\Lambda}^{\Lambda}d\varepsilon\tanh\left(\frac{\beta\sqrt{\varepsilon^{2}+\Delta^{2}}}{2}\right)\frac{1}{2\sqrt{\varepsilon^{2}+\Delta^{2}}}.$
(S144)

This is exactly the gap equation for a BCS superconductor of electrons with a
Fermi surface if $\mathfrak{a}$ is replaced by the constant DOS at the Fermi
level. At zero temperature, we have

$\Delta=\frac{\Lambda}{\sinh\left(\frac{1}{U\mathfrak{a}}\right)}\approx
2\Lambda\,e^{-1/(U\mathfrak{a})}.$ (S145)

By contrast, if we instead employ the fractal dimension for low-energy WZNW
wave functions $d_{2}=3/2$, we would get the non-BCS result if we neglect the
UV cutoff [8, 9]

$\displaystyle\Delta\propto U^{1/7}.$ (S146)

In both scenarios, the finite-threshold strength in the clean limit $U_{c}$
[Eq. (S138)] is lowered to zero for quantum-critical wave functions, _despite_
the vanishing density-of-states for class CI
$\nu(\varepsilon)\approx|\varepsilon|^{1/7}.$

Figure S7: The dependence of the pairing amplitude $\Delta$ on the interaction
strength $U$ for Dirac fermions in both the clean case and in the chiral qBM
model.

The critical exponents in the qBM model deviate from the case of dirty Dirac
fermions. But as is shown in Fig. S3, the DOS in the lake shows power-law
dependence on energy. In Fig. S7, we show the dependence of the pairing
amplitude on the interaction strength for Dirac fermions in the clean limit
and in the critical lake of the chiral qBM model. In the clean limit, $\Delta$
becomes nonzero when $U$ is above the critical interaction strength. In the
lake with $A_{1}=1$ and $A_{2}=0.6$, the $\Delta$ vs. $1/U$ curve is
approximately linear for weak interactions in the log-linear plot, which is
consistent with the prediction in Eq. (S145).

## References

  * [1] Alexander Weiße, Gerhard Wellein, Andreas Alvermann, and Holger Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).
  * [2] Andreas Piepera, Moritz Kreutzerb, Andreas Alvermanna , Martin Galgonc, Holger Fehskea, Georg Hagerb, Bruno Langc, and Gerhard Welleinb, _High-performance implementation of Chebyshev filter diagonalization for interior eigenvalue computations_ , J. Comput. Phys. 325, 226 (2016).
  * [3] Haoyu Guan, Wenxian Zhang, Dual applications of Chebyshev polynomials method: Efficiently finding thousands of central eigenvalues for many-spin systems, SciPost Phys. 11, 103 (2021).
  * [4] G. P. Srivastava, _Broyden’s method for self-consistent field convergence acceleration_ , J. Phys. A: Math. Gen. 17 L317 (1984).
  * [5] D. D. Johnson, _Modified Broyden’s method for accelerating convergence in self-consistent calculations_ , Phys. Rev. B 38, 12807 (1988).
  * [6] D. J. Scalapino, S. R. White, and S. C. Zhang, _Superfluid Density and the Drude Weight of the Hubbard Model_ , Phys. Rev. Lett. 68, 2830 (1992).
  * [7] S. A. A. Ghorashi, Y. Liao, and M. S. Foster, _Critical Percolation without Fine-Tuning on the Surface of a Topological Superconductor_ , Phys. Rev. Lett. 121, 016802 (2018).
  * [8] M. V. Feigel’man, L. B. Ioffe, V. E. Kravtsov, and E. A. Yuzbashyan, _Eigenfunction Fractality and Pseudogap State near the Superconductor-Insulator Transition_ , Phys. Rev. Lett. 98, 027001 (2007).
  * [9] M. V. Feigel’man, L. B. Ioffe, V. E. Kravtsov, and E. Cuevas, _Fractal superconductivity near the localization threshold_ , Ann. Phys. (NY) 325, 1390 (2010).
  * [10] B. Uchoa, G. G. Cabrera, and A. H. Castro Neto, _Nodal liquid and $s$-wave superconductivity in transition metal dichalcogenides_, Phys. Rev. B 71, 184509 (2005).
  * [11] Bruno Uchoa and A. H. Castro Neto, _Superconducting States of Pure and Doped Graphene_ , Phys. Rev. Lett. 98, 146801 (2007).
  * [12] N. B. Kopnin and E. B. Sonin, _BCS Superconductivity of Dirac Electrons in Graphene Layers_ , Phys. Rev. Lett. 100, 246808 (2008); B. Uchoa and A. H. Castro Neto, _Comment on “BCS Superconductivity of Dirac Electrons in Graphene Layers”_ , _ibid._ 102, 109701 (2009); N. B. Kopnin and E. B. Sonin, _Kopnin and Sonin Reply:_ , _ibid._ 102, 109702 (2009).
  * [13] Rahul Nandkishore, Joseph Maciejko, David A. Huse, and S. L. Sondhi, _Superconductivity of disordered Dirac fermions_ , Phys. Rev. B 87, 174511 (2013).
  * [14] Ionut-Dragos Potirniche, Joseph Maciejko, Rahul Nandkishore, and S. L. Sondhi, _Superconductivity of disordered Dirac fermions in graphene_ , Phys. Rev. B 90, 094516 (2014).
  * [15] Matthew S. Foster, Hong-Yi Xie, and Yang-Zhi Chou, _Topological protection, disorder, and interactions: Survival at the surface of three-dimensional topological superconductors_ , Phys. Rev. B 89, 155140 (2014).
  * [16] Yang-Zhi Chou and Matthew S. Foster, _Chalker scaling, level repulsion, and conformal invariance in critically delocalized quantum matter: Disordered topological superconductors and artificial graphene_ Phys. Rev. B 89, 165136 (2014).
  * [17] J. F. Karcher and M. S. Foster, _How spectrum-wide quantum criticality protects surface states of topological superconductors from Anderson localization: Quantum Hall plateau transitions (almost) all the way down_ , Ann. Phys. 435, 168439 (2021).