We consider the chiral model of twisted bilayer graphene introduced by Tarnopolsky, Kruchkov, and Vishwanath (TKV). TKV proved that for inverse twist angles *α* such that the effective Fermi velocity at the moiré *K* point vanishes, the chiral model has a perfectly flat band at zero energy over the whole Brillouin zone. By a formal expansion, TKV found that the Fermi velocity vanishes at *α* ≈ 0.586. In this work, we give a proof that the Fermi velocity vanishes for at least one *α* between 0.57 and 0.61 by rigorously justifying TKV’s formal expansion of the Fermi velocity over a sufficiently large interval of *α* values. The idea of the proof is to project the TKV Hamiltonian onto a finite-dimensional subspace and then expand the Fermi velocity in terms of explicitly computable linear combinations of modes in the subspace while controlling the error. The proof relies on two propositions whose proofs are computer-assisted, i.e., numerical computation together with worst-case estimates on the accumulation of round-off error, which show that round-off error cannot possibly change the conclusion of computation. The propositions give a bound below on the spectral gap of the projected Hamiltonian, an Hermitian 80 × 80 matrix whose spectrum is symmetric about 0, and verify that two real 18-th order polynomials, which approximate the numerator of the Fermi velocity, take values with a definite sign when evaluated at specific values of *α*. Together with TKV’s work, our result proves the existence of at least one perfectly flat band of the chiral model.

## I. INTRODUCTION

### A. Outline

Twisted bilayer graphene (TBG) is formed by stacking one layer of graphene on top of another in such a way that the Bravais lattices of the layers are twisted relative to each other. For generic twist angles, the lattices will be incommensurate, so the resulting structure will not be periodic. Bistritzer and MacDonald (BM)^{1} introduced an approximate model (BM model) for the electronic states of TBG, which is periodic over the scale of the bilayer moiré pattern, where the twist angle enters as a parameter. Using this model, BM showed that the Fermi velocity, the velocity of electrons at the Fermi level, vanishes at particular twist angles known as “magic angles.” The largest of these angles, known as the first magic angle, is at *θ* ≈ 1.1°. Numerical computations on the BM model show the stronger result that at magic angles, the Bloch band of the BM model at zero energy is approximately flat *over the whole Brillouin zone*.^{1,2} The flatness of the zero energy Bloch band is thought to be a critical ingredient for recently observed superconductivity of TBG,^{3} although the precise mechanism for superconductivity in TBG is not yet settled.

Aiming at a simplified model that explains the nearly flat band of TBG, Tarnopolsky, Kruchkov, and Vishwanath (TKV)^{4} introduced a simplification of the BM model, which has an additional “chiral” symmetry, known as the chiral model. TKV showed analytically that at magic angles (of the chiral model, still defined by vanishing of the Fermi velocity), the chiral model has *exactly* flat bands over the whole Brillouin zone. Using a formal perturbation theory (for the chiral model, the natural parameter is the reciprocal of the twist angle up to a constant), TKV derived approximate values for the magic angles of the chiral model. It is worth noting that the first magic angles of the chiral model and the BM model are nearby, but the higher magic angles are not very close.

Becker *et al.* introduced a spectral characterization of magic angles of the TKV model where the role of a non-normal operator is emphasized [the operator *D*^{α} appearing in (1)]. Using this characterization, they numerically computed precise values for the magic angles of the TKV model [see the discussion below (14)].^{5} In the same work, they also proved that the lowest band of the TKV model becomes exponentially close to flat even away from magic angles as the natural small parameter tends to zero. The same authors also investigated flat bands of the TKV model with more general interlayer coupling potentials and the spectrum of other special cases of the BM model.^{6}

In this work, we study the chiral model introduced by TKV and consider the problem of rigorously proving the existence of the first magic angle. We do this by justifying the formal perturbation theory of TKV to make a rigorous expansion of the Fermi velocity to high enough order, and over a large enough parameter range, so that we can prove the existence of a zero. By numerically verifying that the resulting expansion attains a negative value and proving that the result continues to hold when the effect of round-off error is included (Proposition II.2), we obtain existence of the magic angle (Theorem II.2).

The proof of validity of the expansion is challenging because the reciprocal of the twist angle at the zero of the Fermi velocity is large relative to the spectral gap of the unperturbed Hamiltonian, which means that the magic angle falls outside of the interval of twist angles where the perturbation series for the Fermi velocity is obviously convergent. To overcome this difficulty, we start by representing the chiral model Hamiltonian in a basis that takes full advantage of model symmetries. Then, using a rigorous bound on the high frequency components of the error, we reduce the error analysis to analysis of the eigenvalues of the chiral model projected onto finitely many low frequencies. The final stage of the error analysis (Theorem II.1) is to prove a proposition about the eigenvalues of the projected chiral model by a numerical computation that we prove continues to hold when the accumulation of round-off error is considered (Proposition IV.5). We discuss the limitations of our methods and, in particular, whether our methods might be generalized to the more general settings considered by Becker *et al.*^{5,6} in Remarks II.1–II.3.

### B. Code availability

We have made code for the numerical computations used in our proofs available at github.com/abwats/magic_angle. We give references to specific scripts in the text.

## II. STATEMENT OF RESULTS

### A. Tarnopolsky–Kruchkov–Vishwanath’s chiral model

The chiral model, such as the Bistritzer–MacDonald model (B–M model) from which it is derived, is a formal continuum approximation to the atomistic tight-binding model of twisted bilayer graphene. The BM and chiral models aim to capture physics over the length-scale of the bilayer moiré pattern, which is, for small twist angles, much longer than the length-scale of the individual graphene layer lattices. Crucially, even when the graphene layers are incommensurate so that the bilayer is aperiodic on the atomistic scale, the chiral model and BM model are periodic (up to phases) with respect to the moiré lattice so that they can be analyzed via Bloch theory.

We define the moiré lattice to be the Bravais lattice

generated by the moiré lattice vectors

and we denote a fundamental cell of the moiré lattice by Ω. The moiré reciprocal lattice is the Bravais lattice

generated by the moiré reciprocal lattice vectors defined by *a*_{i} · *b*_{j} = 2*πδ*_{ij}, given explicitly by

We define $q1=0,\u22121$, which is the (re-scaled) difference of the *K* points (Dirac points) of each layer, and

We write Ω^{*} for a fundamental cell of the moiré reciprocal lattice and refer to such a cell as the Brillouin zone.

Let $\varphi \u22542\pi 3$. Tarnopolsky–Kruchkov–Vishwanath’s chiral Hamiltonian is defined as

where $\u2202\u0304=12(\u2202x+i\u2202y)$, $U(r)=e\u2212iq1\u22c5r+ei\varphi e\u2212iq2\u22c5r+e\u2212i\varphi e\u2212iq3\u22c5r$, † denotes the adjoint (Hermitian transpose), and *α* is a real parameter, which we will take to be positive *α* ≥ 0 throughout [see (3)]. The chiral Hamiltonian *H*^{α} is an unbounded operator on $H=L2(R2;C4)$ with domain $H1(R2;C4)$. We will write functions in $H$ as

where $|\psi \tau \sigma (r)|2$ represents the electron density near to the *K* point (in momentum space) on sublattice *σ* and on layer *τ*. The diagonal terms of *D*^{α} arise from Taylor expanding the single layer graphene dispersion relation about the *K* point of each layer, while the off-diagonal terms of *D*^{α} couple the *A* and *B* sublattices of layers 1 and 2. The chiral model is identical to the BM model except that inter-layer coupling between sublattices of the same type is turned off in the chiral model. The precise form of the interlayer coupling potential *U* can be derived under quite general assumptions on the real space interlayer hopping.^{1,7} The parameter *α* is, up to unimportant constants, the ratio

Although the limit *α* → 0 can be thought of as the limit of vanishing interlayer hopping strength at fixed twist, it is physically more interesting to view the limit as modeling increasing twist angle at a fixed interlayer hopping strength.

### B. Rigorous justification of TKV’s formal expansion of the Fermi velocity and proof of existence of first magic angle

Bistritzer and MacDonald studied the effective Fermi velocity of electrons in twisted bilayer graphene modeled by the BM model and computed values of the twist angle such that the Fermi velocity vanishes, which they called “magic angles.” One can similarly define an effective Fermi velocity for the chiral model and refer to values of *α* such that the Fermi velocity vanishes as “magic angles” [although technically *α* is related to the reciprocal of the twist angle (3)].

TKV proved the remarkable result that, at magic angles, the chiral model has a perfectly flat Bloch band at zero energy. Let $LK2$ denote the *L*^{2} space on a single moiré cell Ω with moiré *K* point Bloch boundary conditions. The starting point of TKV’s proof is an expression for the Fermi velocity as a function of *α*, *v*(*α*), as a functional of one of the Bloch eigenfunctions, $\psi \alpha \u2208LK2$, of *H*^{α},

where $..$ denotes the $LK2$ inner product. We give precise definitions of $LK2$, *ψ*^{α}, and *v*(*α*) in Definition III.2, Proposition III.6, and Definition III.3, respectively. We give a systematic formal derivation of why (4) is the effective Fermi velocity at the moiré *K* point in Appendix A. To complete the proof, TKV showed that zeros of *v*(*α*) imply zeros of *ψ*^{α} at special “stacking points” of Ω and that such zeros of *ψ*^{α} allow for Bloch eigenfunctions with zero energy to be constructed for all ** k** in the moiré Brillouin zone.

To derive approximate values for magic angles, TKV computed a formal perturbation series approximation of *ψ*^{α},

and then substituted this expression into the functional for *v*(*α*) to obtain an expansion of *v*(*α*) in powers of *α*,

By setting *v*(*α*) = 0, one obtains an approximation for the smallest magic angle: *α* ≈ 0.586.

Although TKV proved that flat bands occur at magic angles, they did not prove the existence of magic angles, and hence, they did not prove the existence of flat bands. The contribution of the present work is to prove rigorous estimates on the error in approximation (5), which are sufficiently high order and precise that, once substituted into (4), they suffice to rigorously prove the existence of a zero of *v*(*α*) and, hence, via TKV’s proof, the existence of at least one perfectly flat band.

It turns out to be relatively straightforward to prove that series (5) and (6) are uniformly convergent and to derive precise bounds on the error in truncating the series for $|\alpha |<13$; see Proposition IV.3. The basic challenge, then, is to derive similar error bounds for *α* over an interval, which includes the expected location of the first magic angle, at $\u224813$. The first main theorem we will prove, roughly stated, is the following. See Theorem IV.1 for the more precise statement. The theorem relies on the existence of a spectral gap for an 80 × 80 Hermitian matrix, which requires numerical computation for its proof; see Proposition IV.5.

*The*

*K*

*point Bloch function*

*ψ*

^{α}

*satisfies*

*where*$\eta \alpha \u22a5\u2211n=08\alpha n\Psi n$

*with respect to the*$LK2$

*inner product and*

The functions Ψ^{n} for 0 ≤ *n* ≤ 8 are derived recursively: see Appendix C. We stop at eighth order in the expansion because this is the minimal order such that we can guarantee the existence of a zero of *v*(*α*), but the functions Ψ^{n} are well-defined by a recursive procedure for arbitrary positive integers *n*; see Proposition IV.1.

Substituting (7) into the functional for Fermi velocity (4) and using $\eta \alpha \u22a5\u2211n=18\alpha n\Psi n$, we find that

where

and

where $..$ denotes the $LK2$ inner product and *η*^{α} satisfies (8). The following is a straightforward calculation.

*The following identities hold:*

We prove Proposition II.1 in Appendix F. Naïvely, expansions (10) and (11) approximate the formal infinite series expansions of $\u2211n=0\u221e\alpha n\Psi n*(\u2212r)\u2211n=0\u221e\alpha n\Psi n(r)$ and $\u2211n=0\u221e\alpha n\Psi n\u2211n=0\u221e\alpha n\Psi n$ up to terms of order *α*^{9}. We prove in Proposition F.2 that because of some simplifications, expansions (10) and (11) agree with the infinite series up to terms of order *α*^{10}.

We are now in a position to state and prove our second result. This result also relies on a proposition, which requires numerical computation for its proof: that one real 18-th order polynomial in *α* attains a negative value and another attains a positive value, when evaluated at specific values of *α*; see Proposition II.2.

*There exist positive numbers* *α*_{min} *and* *α*_{max} *with* 0.57 < *α*_{min} < *α*_{max} < 0.61 *such that the Fermi velocity* *v*(*α*) *defined by* *(4)* *has a zero* *α*^{*} *satisfying* *α*_{min} ≤ *α*^{*} ≤ *α*_{max}*.*

*v*(

*α*) is bounded below for all $0\u2264\alpha \u2264710$ by the polynomial

Proposition II.2 obviously implies by continuity (12), (13), and $vN(\alpha )$, each has at least one zero in the interval 0.57 < *α* < 0.61. We denote the largest zero of (12) in the interval by *α*_{max} and the smallest zero of (13) in the interval by *α*_{min}. Since the zeroes of $vN(\alpha )$ must lie between those of (12) and (13), we are done.□

(*computer-assisted*). We will first prove that (12) attains a negative value at 0.61 and then explain the modifications necessary to prove that (13) is positive at 0.57. Evaluating using the double-precision floating-point arithmetic, we find that at *α* = 0.61, (12) attains the negative value −0.020 263 (five significant figures; this value was computed by running the script compute_expansion_symbolically.py in the Github repo). It is straightforward to bound the numerical error, which accumulates when evaluating an 18-th order polynomial using the floating-point arithmetic. Even the simplest exact bound, which does not account for error cancellation (see, e.g., Eq. (8) of Oliver),^{8} yields an upper bound on the possible accumulated round-off error in the evaluation of an *n*th order polynomial $\u2211j=0npj\alpha j$, for *α* ∈ [−1, 1], as $(n+1)e(2n+1)\u03f5\u22121sup0\u2264j\u2264n|pj|$, where *ϵ* is “machine epsilon”: roughly speaking, the maximum possible round-off error generated in a single arithmetic operation. Bounding the maximum coefficient in (12) by 1000, taking *n* = 18, and bounding *ϵ* by 3 × 10^{−16} (which was easily attained working in Python on our machine), the maximum possible numerical error in the evaluation of (12) is ≈10^{−11}, which is much smaller than 0.020 263. We conclude that the first claim of Proposition II.2 must hold. Regarding the second, evaluating at *α* = 0.57, we find that (13) equals 0.029 138 (5sf). The same argument as before now shows that accumulated round-off error in the evaluation cannot possibly change the sign of (13) at *α* = 0.57.□

We do not attempt to rigorously estimate *α*_{min} and *α*_{max} precisely in this work, but numerically computing roots of polynomials (12) and (13) suggests *α*_{min} ≈ 0.576 83 (5sf) and *α*_{max} ≈ 0.601 77 (5sf), respectively, where 5sf is an abbreviation for five significant figure. Numerical computation of the first zero of $\u2211n=08\alpha n\Psi n*(\u2212r)\u2211n=08\alpha n\Psi n(r)$ gives 0.585 97 (5sf); see Fig. 1 (the zero values were computed by running the script compute_expansion_symbolically.py in the Github repo).

Using Proposition C.1 and the package Sympy^{9} for symbolic computation, we can compute the formal expansion of *v*(*α*) up to arbitrarily high order in *α*. In particular, we find the higher-order terms in expansion (6) to be

Truncating the numerator after order *α*^{40} and setting the numerator equal to zero yield *α* = 0.585 663 558 389 56 (14sf) for the first zero of the Fermi velocity (to compute this value, run compute_expansion_symbolically.py in the Github repo with *N* = 40). This is consistent with the numerical computation of Becker *et al.*,^{5} who found *α* = 0.585 663 558 389 55 truncated (not rounded) to 14 digits, by diagonalizing a non-normal but compact operator whose reciprocal eigenvalues correspond to magic angles. Note that we do not attempt to rigorously justify series (14) to such large values of *α* and to such high order in this work; see Remark II.4.

(Higher magic angles). *The chiral model has been conjectured to have infinitely many magic angles,*^{5} *but it is not straightforward to extend our methods to prove the existence of such higher magic angles. The problem is that calculating the perturbation series centered at* *α* = 0 *requires diagonalizing the unperturbed operator* *H*^{0}*. In principle, it might be possible to calculate the perturbation series to higher order in order to get an accurate approximation of the Fermi velocity near to the higher magic angles. However, this would require significantly more calculation compared with the present work, and we have no guarantee that the error can be made small enough to prove the existence of another zero in that case.*

(More general interlayer hopping potentials). *The chiral model* *(1)* *is an approximation to the full Hamiltonian of the twisted bilayer, even in the chiral limit where coupling between sublattices of the same type is turned off because the interlayer hopping potential* *U* *only allows for hopping between nearest neighbors in the momentum lattice (see* Fig. 6 *). More general interlayer hopping potentials have been studied by Becker et al.*^{6} *In principle, such models should be amenable to the analysis of this work, but longer-range hopping would lead to much more involved calculations, and the construction of the finite-dimensional subspace* Ξ *of Proposition IV.4 would require more care: the fact that we can choose* Ξ *so that* $\Vert P\Xi H1P\Xi \u22a5\Vert =1$ *depends on* *H*^{1} *only coupling nearest neighbors in the momentum lattice. Locality of hopping in the momentum space lattice has been exploited for efficient computation of density of states*^{10} *of twisted bilayers.*

(Generalization to BM model). *Parts of our analysis should also apply to the full Bistritzer–MacDonald model. Specifically, one could study perturbation series for Bloch functions near to zero energy in powers of the inter-layer hopping strength, derive an equivalent expression for the Fermi velocity in terms of that series, and then study the zeroes of that series. However, there are various complications because of the lack of “chiral” symmetry. First, there is no reason for the continuation of the zero eigenvalue of the unperturbed operator to remain at zero. Second, the expression for the Fermi velocity in terms of the associated eigenfunction could be more complicated. Since zeros of the BM model Fermi velocity do not imply the existence of flat bands for that model, we do not consider these complications in this* *work.*

*(Expanding to higher order). Our methods could, in principle, be continued to justify the expansion of the Fermi velocity to arbitrarily high order and potentially over larger intervals of* *α* *values. However, these extensions are not immediate: pushing the expansion to higher order or to a larger interval of* *α* *values would require a larger set* Ξ *in Lemma IV.1, and Proposition IV.5 would have to be re-proved for the new set* Ξ*. Note that the essential difficulty is justifying the perturbation series for large* *α**: the series are easily justified to all orders for* $|\alpha |<13$; *see Proposition IV.3.*

### C. Structure of paper

We review the symmetries, Bloch theory, and symmetry-protected zero modes of TKV’s chiral model in Sec. III. We prove Theorem II.1 in Sec. IV, postponing most details of the proofs to Appendixes A– F. In Appendix A, we show why (4) corresponds to the effective Fermi velocity at the moiré *K* point. In Appendix B, we construct an orthonormal basis, which we refer to as the chiral basis, which allows for efficient computation and analysis of TKV’s formal expansion. We re-derive TKV’s formal expansions in Appendix C. We give details of the Proof of Theorem II.1 in Appendixes D and E. We prove Proposition II.1 in Appendix F. In the supplementary material, we list the basis functions of the subspace onto which we project the TKV Hamiltonian, give the explicit forms of the higher-order corrections in expansion (7), and present a derivation of the TKV Hamiltonian from the Bistritzer–MacDonald model.

## III. SYMMETRIES, BLOCH THEORY, AND ZERO MODES OF TKV’s CHIRAL MODEL

### A. Symmetries of the TKV model

In this section, we review the symmetries of the TKV model for the reader’s convenience and to fix notation. Becker *et al.*^{5} gave a group theoretical account of these symmetries, and further reviews can be found in the physics literature.^{11–13} Recall that $\varphi =2\pi 3$, and let *R*_{ϕ} denote the matrix, which rotates vectors counter-clockwise by *ϕ*, i.e.,

We define the following.

*For any*

**∈ Λ,**

*v**we define a phase-shifted translation operator acting on functions*$f\u2208H$

*by*

*We define a phase-shifted version of the operator, which rotates functions*$f\u2208H$

*clockwise by*

*ϕ*

*by*

*For any*$f\u2208H$,

*we finally define the “chiral” symmetry operator*

We then have the following.

**∈ Λ,**

*v*The “chiral” symmetry (19) implies that the spectrum of *H*^{α} is symmetric about zero because

The same calculation implies that zero modes of *H*^{α} can always be chosen without loss of generality to be eigenfunctions of $S$.

### B. Bloch theory for the TKV Hamiltonian

We now want to reduce the eigenvalue problem for *H*^{α} using the symmetries just introduced. The symmetry (18) means that eigenfunctions of *H*^{α} can be chosen without loss of generality to be simultaneous eigenfunctions of *τ*_{v} for all ** v** ∈ Λ. It therefore suffices to seek solutions of

for ** r** in a fundamental cell $\Omega \u2254R2/\Lambda $ of the moiré lattice in the symmetry-restricted spaces

where ** k** is known as the quasimomentum. Since $Lk+w2=Lk2$ for any

**∈ Λ**

*w*^{*}, it suffices to restrict attention to

**in a fundamental cell of Λ**

*k*^{*}, which we denote $\Omega *\u2254R2/\Lambda *$ and refer to as the Brillouin zone. We also define symmetry-restricted Sobolev spaces $Hks$ for each

**∈ Ω**

*k*^{*}and positive integer

*s*by

We claim the following.

*For each fixed* ** k** ∈ Ω

^{*}

*and*

*α*≥ 0

*,*

*H*

^{α}

*, defined on the domain*$Hk1$

*, extends to an unbounded self-adjoint elliptic operator*$Lk2\u2192Lk2$

*with compact resolvent. In a complex neighborhood of every*

*α*≥ 0

*, the family*

*H*

^{α}

*is a holomorphic family of type (A) in the sense of Kato.*

^{14}Ellipticity is immediate since the principal symbol of *H*^{α} is invertible. Self-adjointness is clear using the Fourier transform when *α* = 0 and for *α* ≠ 0 because *αH*^{1} is a bounded symmetric perturbation of *H*^{0} (see, e.g., Theorem 1.4 of Cycon *et al.*^{15}). Elliptic regularity implies that the resolvent maps $Lk2\u2192Hk1$, and compactness of the resolvent then follows by Rellich’s theorem (see, e.g., Proposition 3.4 of Taylor^{16}). The family *H*^{α} is holomorphic of type (A) since the domain of *H*^{α} is independent of *α*, and *H*^{α}*f* is holomorphic for every $f\u2208Hk1$ (see Chap. 7 of Kato^{14}).□

We now claim the following.

*Let* $f\u2208Lk2$*. Then,* $Rf\u2208LR\varphi *k2$*.*

**∈ Λ,**

*v*

*b*_{2}·

**= 0 mod 2**

*v**π*for all

**∈ Λ.□**

*v*In particular, whenever $R\varphi *k=k$ mod Λ^{*}, we have $RLk2=Lk2$. Regarding such ** k**, the following is a simple calculation.

*The moiré* *K* *and* *K*′ *points* ** k** = 0

*and*

**= −**

*k*

*q*_{1}

*and the moiré*Γ

*point*

**=**

*k*

*q*_{1}+

*b*_{1}

*satisfy*$R\varphi *k=k$ mod Λ

^{*}

*.*

The moiré *K*, *K*′, and Γ points are shown in Fig. 2. Note that the moiré *K*, *K*′, and Γ points should not be confused with the single layer *K*, *K*′, and Γ points. The moiré *K* point corresponds to the *K* point of layer 1, while the moiré *K*′ point corresponds to the *K* point of layer 2. Interactions with the *K*′ points of layers 1 and 2 are formally small for small twist angles and are hence ignored.

In this work, we will be particularly interested in Bloch functions at the moiré *K* and *K*′ points. We therefore define the following.

Let *ω* = *e*^{iϕ}. Since the spaces $LK2$ and $LK\u20322$ are invariant under $R$, they can be divided up into invariant subspaces corresponding to the eigenvalues of $R$,

where

and $LK\u2032,\sigma 2,\sigma =1,\omega ,\omega *$, are defined similarly.

The following, which is trivial to prove, will be important for studying the zero modes of *H*^{α}.

*The operator* $S$ *commutes with* *τ*_{v} *and* $R$ *and hence maps the* $LK,\sigma 2$ *and* $LK\u2032,\sigma 2$ *spaces to themselves for* *σ* = 1, *ω*, *ω*^{*}*.*

Since $S$ has eigenvalues ±1, we can define the spaces

and spaces $LK\u2032,\sigma ,\xb112,\sigma =1,\omega ,\omega *$ similarly.

*Note that* +1 *eigenspaces of* $S$ *correspond to wave-functions, which vanish in their third and fourth entries, which correspond, through* *(2)**, to wave-functions supported only on* *A* *sites of the layers. Similarly,* −1 *eigenspaces of* $S$ *correspond to wave-functions, which vanish in their first and second entries, which are supported only on* *B* *sites of the layers.*

### C. Zero modes of the chiral model

We now want to investigate zero modes of *H*^{α} in detail. When *α* = 0, there are exactly four zero modes given by *e*_{j}, *j* = 1, 2, 3, 4, where *e*_{j} equals 1 in its *j*th entry and 0 in its other entries. It is easy to check that

and hence, 0 is a simple eigenvalue of *H*^{α} when restricted to each of these subspaces. Recall that zero modes can always be chosen as eigenfunctions of $S$, and indeed, we have

We now claim that these zero modes persist for all *α*. This was already established by TKV,^{4} and the following proposition is also similar to Proposition 3.1 of Becker *et al.*^{5} We re-state it using our notation and give the proof for completeness.

*There exist smooth functions* *ψ*^{α} *with* ‖*ψ*^{α}‖ = 1 *in each of the spaces* $LK,1,12$*,* $LK\u2032,1,12$*,* $LK,\omega *,\u221212$*, and* $LK\u2032,\omega *,\u221212$ *such that* *ψ*^{0} *is as in* *(21)**,* *α* ↦ *ψ*^{α} *is real-analytic, and* *H*^{α}*ψ*^{α} = 0 *for all* *α**. The dimension of* ker *H*^{α} *restricted to each of the spaces* $LK,12$*,* $LK\u2032,12$*,* $LK,\omega *2$*, and* $LK,\omega *2$ *is always odd-dimensional.*

Since $S$ preserves each of the spaces $LK,12$, $LK\u2032,12$, $LK,\omega *2$, and $LK,\omega *2$ and anti-commutes with *H*^{α}, the spectrum of *H*^{α} restricted to each space must be symmetric about 0 for all *α*. Since *H*^{α} restricted to each space has compact resolvent and *H*^{α} is a holomorphic family of type (A), the spectrum of *H*^{α} consists of finitely degenerate isolated eigenvalues depending real-analytically on *α*, with associated eigenfunctions also real-analytically depending on *α* (although the real-analytic choice of eigenfunction at an eigenvalue crossing may not respect ordering); see Theorem 3.9 of Chap. 7 of Kato.^{14} The null space of *H*^{α} in each of the spaces is one-dimensional at *α* = 0 by explicit calculation, with the zero modes given by (up to non-zero constants) (21). For small *α* > 0, real-analyticity and the chiral symmetry force the null space to remain simple and it is clear how to define *ψ*^{α}. For large *α* > 0, the non-zero eigenvalues of *H*^{α} may cross 0 at isolated values of *α*, and in this case, we define *ψ*^{α} to be the real-analytic continuation of the zero mode through the crossings. Note that real-analyticity prevents non-zero eigenvalues from equaling zero except at isolated points so that the real-analytic continuation of the zero mode through the crossing must indeed be a zero mode. At crossings, the null space must be odd-dimensional in order to preserve symmetry of the spectrum of *H*^{α} about 0. It remains to check that if *ψ*^{0} is in, say, $LK,1,12$, then *ψ*^{α} must remain in $LK,1,12$ for all *α* > 0. However, this must hold because the $S$-eigenvalue of *ψ*^{α} cannot change abruptly while preserving real-analyticity. Smoothness of *ψ*^{α} follows from elliptic regularity.□

In this work, we will restrict attention to the moiré *K* point and, especially, the family $\psi \alpha \u2208LK,1,12$. We expect that our analysis would go through with only minor modifications if we considered instead the moiré *K*′ point. The zero modes in $LK,1,12$ and $LK,\omega *,\u221212$ are related by the following symmetry.

*Let* $\psi 1\alpha $ *and* $\psi \u22121\alpha $ *denote the zero modes of* *H*^{α} *in the spaces* $LK,1,12$ *and* $LK,\omega *,\u221212$, *respectively. Then,* $\psi 1\alpha =(\Phi \alpha ,0)\u22a4$, *where* $\Phi \alpha \u2208L2(\Omega ;C2)$*,* $\Phi \alpha (r+v)=diag(1,eiq1\u22c5v)\Phi \alpha (r)$ *for all* ** v** ∈ Λ

*,*Φ

^{α}(

*R*

_{ϕ}

**) = Φ**

*r*^{α}(

**)**

*r**. Up to gauge transformations*$\psi \u22121\alpha \u21a6ei\varphi (\alpha )\psi \u22121\alpha $,

*which preserve real-analyticity of*$\psi \u22121\alpha $

*, we have*$\psi \u22121\alpha (r)=(0,\Phi \alpha *(\u2212r))\u22a4$

*.*

Since $S\psi 1\alpha =\psi 1\alpha $, the last two entries of $\psi 1\alpha $ must vanish, so we can write $\psi 1\alpha =(\Phi \alpha ,0)\u22a4$. That Φ^{α} satisfies the stated symmetries follows immediately from $\psi 1\alpha \u2208LK,12$. It is straightforward to check using the definitions of $R$ and *τ*_{v} that $(0,\Phi \alpha *(\u2212r))\u22a4\u2208LK,\omega *,\u221212$. To see that $(0,\Phi \alpha *(\u2212r))\u22a4$ is a zero mode, note that Φ^{α} satisfies *D*^{α}Φ^{α} = 0, which implies that $D\alpha \u2020\Phi \alpha *(\u2212r)=0$ by a simple manipulation. To see that $\psi \u22121\alpha (r)=(0,\Phi \alpha *(\u2212r))\u22a4$ (up to real-analytic gauge transformations) for all *α*, note first that this clearly holds for *α* = 0 [the zero modes are explicit (22)]. For *α* > 0, the identity must continue to hold by uniqueness (up to real-analytic gauge transformations) of the real-analytic continuation of $\psi \u22121\alpha $ starting from *α* = 0 and continuing first along the non-zero interval where $\psi 1\alpha $ is non-degenerate in $LK,1,12$ and then through eigenvalue crossings as in the Proof of Proposition III.6.□

In Appendix A, we use Proposition III.7 to derive the effective Dirac operator with the *α*-dependent Fermi velocity, which controls the Bloch band structure in a neighborhood of the moiré *K* point. The Fermi velocity of the effective Dirac operator is given by the following. Note that we drop the subscript +1 when referring to the zero mode of *H*^{α} in $LK,1,12$ since the zero mode of *H*^{α} in $LK,\omega *,\u221212$ plays no further role.

*Let*$\psi \alpha \u2208LK,1,12$

*be as in Proposition III.6. Then, we define*

*where*$..$

*denotes the*$LK2$

*inner product.*

## IV. RIGOROUS JUSTIFICATION OF TKV’s EXPANSION OF THE FERMI VELOCITY

### A. Alternative formulation of TKV’s expansion

We now turn to approximating the zero mode $\psi \alpha \u2208LK,1,12$ by a series expansion in powers of *α*. We write *H*^{α} = *H*^{0} + *αH*^{1} and formally expand *ψ*^{α} as a series

where *H*^{0}Ψ^{0} = 0 and

for all *n* ≥ 1. To solve *H*^{0}Ψ^{0} = 0, we take Ψ^{0} = *e*_{1}. We prove the following in Appendix C.

*Let*

*P*

*denote the projection operator in*$LK,12$

*onto*

*e*

_{1}

*and*

*P*

^{⊥}=

*I*−

*P*

*. The sequence of*

*Eq. (25)*

*has a unique solution such that*$\Psi n\u2208LK,1,12$

*for all*

*n*≥ 0

*and*

*P*Ψ

^{n}= 0

*for all*

*n*≥ 1,

*given by*Ψ

^{0}=

*e*

_{1}

*and*

*for each*

*n*≥ 1

*.*

Expansion (24) appears different from the series studied by TKV since we work only with the self-adjoint operators *H*^{0}, *H*^{1}, and *H*^{α} rather than the non-self-adjoint operator *D*^{α} [defined in (1)]. Since functions in $LK,1,12$ vanish in their last two components, there is no practical difference. However, working with only self-adjoint operators allows us to use the spectral theorem, which greatly simplifies the error analysis. We compute the first eight terms in expansion (24) in Proposition C.2 after developing some necessary machinery in Appendix B.

### B. Rigorous error estimates for the expansion of the moiré *K* point Bloch function

In this section, we explain the essential challenge in proving error estimates for series (24) and explain how we overcome this challenge. Our goal is to prove the following.

*Let*$\psi \alpha \u2208LK,1,12$

*be as in Proposition III.6. Then,*

*where*$\eta \alpha \u22a5\u2211n=18\alpha n\Psi n$

*with respect to the*$LK2$

*inner product and*

Proposition IV.1 guarantees that series (24) is well-defined up to arbitrarily many terms. A straightforward bound on the growth of terms in the series comes from the following proposition.

*The spectrum of*

*H*

^{0}

*in*$LK,12$

*is*

*and hence,*

*We also have*

This proposition is a combination of Propositions B.2, B.4, and B.7, proved in Appendix B.□

Proposition IV.2 implies that $\Vert P\u22a5(H0)\u22121P\u22a5H1\Vert LK,12\u2192LK,12\u22643$, which implies the following.

*The formal series* *(24)* *converges to* *ψ*^{α} *in* $LK2$*, with an explicit error rate, for all* $|\alpha |<13$*. The formal series for the Fermi velocity* *v*(*α*) *obtained by substituting the series expansion of* *ψ*^{α} *into* *(23)* *converges for the same range of* *α**, also with an explicit error rate.*

For any non-negative integer *N*, let $\psi N,\alpha \u2254\u2211n=0N\alpha n\Psi n$, where Ψ^{n} are as in (26). Since Ψ^{0} ⊥ Ψ^{n} for all *n* ≥ 1 and ‖Ψ^{0}‖ = 1, we have that ‖*ψ*^{N,α}‖ ≥ 1 for all *N*. Let $\varphi N,\alpha \u2254\psi N,\alpha \Vert \psi N,\alpha \Vert $; then, we can decompose *ϕ*^{N,α} = *cψ*^{α} + *η*^{α} for some constant *c* and where *η*^{α} ⊥ *ψ*^{α}. Applying *H*^{α} to both sides, we have that $H\alpha \varphi N,\alpha =\alpha N+1H1\Psi N\Vert \psi N,\alpha \Vert =H\alpha \eta \alpha $. Now, fix *α* ≥ 0 such that $|\alpha |<13$. Then, *α*‖*H*^{1}‖ < 1, and hence, the first non-zero eigenvalue of *H*^{α} is bounded away from 0 by 1 − 3*α* (recall that the first non-zero eigenvalues of *H*^{0} are ±1). Since *η*^{α} ⊥ *ψ*^{α}, where *ψ*^{α} spans the eigenspace of the zero eigenvalue of *H*^{α}, we have that $\Vert \eta \alpha \Vert \u2264|\alpha N+1|\Vert H1\Psi N\Vert |1\u22123\alpha |\Vert \psi N,\alpha \Vert $. Using the bound ‖Ψ^{N}‖ ≤ (3*α*)^{N} and the bound below on ‖*ψ*^{N,α}‖, we have that $\Vert \eta \alpha \Vert \u2264(3\alpha )N+1|1\u22123\alpha |$, which clearly → 0 as *N* → ∞, so that lim_{N→∞}*ϕ*^{N,α} = *ψ*^{α} (up to a non-zero constant). Now, consider *v*(*α*) defined by (23). Assuming WLOG that ‖*ψ*^{α}‖ = 1 and substituting *ψ*^{α} = *ϕ*^{N,α} + *η*^{α}, we find immediately, using Cauchy–Schwarz, that $v(\alpha )\u2212\varphi N,*(\u2212r)\varphi N(r)\u22642\Vert \eta \alpha \Vert +\Vert \eta \alpha \Vert 2$. In terms of *ψ*^{N}, we have $v(\alpha )\u2212\psi N,*(\u2212r)\psi N(r)\psi N\psi N\u22642\Vert \eta \alpha \Vert +\Vert \eta \alpha \Vert 2$.□

Proposition IV.3 shows that for $|\alpha |<13$, series (24) converges to *ψ*^{α} and can be used to compute the Fermi velocity. However, this restriction is too strong to prove that the Fermi velocity has a zero, which occurs at the larger value $\alpha \u224813$. Of course, Proposition IV.2 establishes only the most pessimistic possible bound on the expansion functions Ψ^{n}, and this bound appears to be far from sharp from explicit calculation of each Ψ^{n}; see Proposition C.3. We briefly discuss a possible route to a tighter bound in Remark C.2 but do not otherwise pursue this approach in this work.

We now explain how to obtain error estimates over a large enough range of *α* values to prove *v*(*α*) has a zero. We seek a solution of *H*^{α}*ψ*^{α} = 0 in $LK,1,12$ with the form

For arbitrary *α*, let *Q*^{α} denote the projection in $LK,12$ onto *ψ*^{N,α}, and *Q*^{α,⊥} ≔ *I* − *Q*^{α} (note that *Q*^{0} = *P*). Note that *Q*^{α} depends on *N*, but we suppress this to avoid clutter. We assume WLOG that *Q*^{α}*η*^{α}(** r**) = 0. It follows that

*η*

^{α}satisfies

To obtain a bound on *η*^{α} in *L*^{2}(Ω), we require a lower bound on the operator $Q\alpha ,\u22a5H\alpha Q\alpha ,\u22a5:Q\alpha ,\u22a5LK,12\u2192Q\alpha ,\u22a5LK,12$. The following lemma gives a lower bound on this operator in terms of a lower bound on the projection of this operator onto the finite-dimensional subspace of $LK,12$ corresponding to a finite subset of the eigenfunctions of *H*^{0}. The importance of this result is that since *H*^{1} only couples finitely many modes of *H*^{0}, for *fixed**N*, by taking the subset sufficiently large, we can always arrange that *ψ*^{N,α} lies in this subspace.

*Let*

*P*

_{Ξ}

*denote the projection onto a subset*Ξ

*of the eigenfunctions of*

*H*

^{0}

*in*$LK,12$

*, and let*

*μ*≥ 0

*be maximal such that*

*(with this notation, the operator*

*P*

*introduced in Proposition IV.1 corresponds to*

*P*

_{Ξ}

*with*Ξ

*being the set*{

*e*

_{1}}

*and*

*μ*= 1

*). Suppose that*

*Q*

^{α}

*P*

_{Ξ}=

*P*

_{Ξ}

*Q*

^{α}=

*Q*

^{α}

*,*i.e.

*, that*

*ψ*

^{N,α}

*lies in*ran

*P*

_{Ξ}

*. Define*

*g*

^{α}

*by*

*We note that*

*P*

_{Ξ}

*Q*

^{α,⊥}

*is the projection onto the subspace of*$P\Xi LK,12$

*orthogonal to*

*ψ*

^{N,α}

*. As long as*

*then*

Note that *g*^{α} would be identically zero if not for the restriction that the matrix acts on $Q\alpha ,\u22a5P\Xi LK,12$ since otherwise *ψ*^{N,α} would be an eigenfunction with eigenvalue zero for all *α*. As it is, *g*^{0} = 1 and *α* ↦ *g*^{α} is real-analytic so that *g*^{α} must be positive for a non-zero interval of positive *α* values.

*Q*

^{α}

*P*

_{Ξ}=

*P*

_{Ξ}

*Q*

^{α}, we have $P\Xi \u22a5Q\alpha ,\u22a5=Q\alpha ,\u22a5P\Xi \u22a5=P\Xi \u22a5$, and hence,

*P*

_{Ξ}

*Q*

^{α,⊥}is a projection, and $\Vert Q\alpha ,\u22a5P\Xi H1P\Xi \u22a5\Vert =\Vert P\Xi \u22a5H1P\Xi Q\alpha ,\u22a5\Vert $. Hence, we can bound

*μ*≥ 3

*α*. We now estimate

*α*≤

*μ*,

*α*≤

*μ*and $\alpha \Vert Q\alpha ,\u22a5P\Xi H1P\Xi \u22a5\Vert \u2264min(g\alpha ,\mu \u22123\alpha )$ upon substituting (33) and (34) into (32).□

For Lemma IV.1 to be useful, we must check that it is possible to choose Ξ so that the bound (31) is non-trivial, i.e., so that the constant is positive. We will prove the following in Appendix D.

The set Ξ constructed in Proposition IV.4 is the set of $LK,12$-eigenfunctions of *H*^{0} with eigenvalues with magnitude $\u226443$, augmented with two extra basis functions to ensure that $\Vert P\Xi H1P\Xi \u22a5\Vert =1$. Including all $LK,12$-eigenfunctions of *H*^{0} with eigenvalue magnitudes up to and including $43$ ensure that *ψ*^{8,α} lies in ran *P*_{Ξ}.

We now require the following.

*Let* Ξ *be as in Proposition IV.4. Then,* $g\alpha \u226534$ *for all* $0\u2264\alpha \u2264710$*.*

*(computer-assisted).* Consider $H\Xi \alpha \u2254Q\alpha ,\u22a5P\Xi H\alpha P\Xi Q\alpha ,\u22a5$ acting on $P\Xi LK,12$. Assuming that *α* is restricted to an interval such that the zero eigenspace of $H\Xi \alpha $ is simple, then, using orthogonality of eigenvectors corresponding to different eigenvalues and the fact that *Q*^{α} is the spectral projection onto the unique zero mode of $H\Xi \alpha $, $H\Xi \alpha $ has the same non-zero eigenvalues as the matrix *Q*^{α,⊥}*P*_{Ξ}*H*^{α}*P*_{Ξ}*Q*^{α,⊥} acting on $Q\alpha ,\u22a5P\Xi LK,12$. The matrix $H\Xi \alpha $ is an 81 × 81 matrix whose spectrum is symmetric about 0 because of the chiral symmetry. When *α* = 0, the spectrum is explicit: 0 is a simple eigenvalue, and the smallest non-zero eigenvalues are ±1, both also simple. Proposition IV.5 is proved if we can prove that the first positive eigenvalue of $H\Xi \alpha $ is bounded away from zero by $34$ for all $0\u2264\alpha \u2264710$. Note that if this holds, the zero eigenspace of $H\Xi \alpha $ must be simple for all $0\u2264\alpha \u2264710$, and hence, our basic assumption is justified. The strategy of the proof is as follows.

Define a grid $G\u22547n10N:n\u2208{0,1,\u2026,N}$, where

*N*is a positive integer taken sufficiently large that the grid spacing $h\u2254710N<1388831$ (the number 388 831 comes from Proposition IV.7).Numerically compute the eigenvalues of $H\Xi \alpha $ for $\alpha \u2208G$. We find that the numerically computed first positive eigenvalues of these matrices are uniformly bounded below by $810>34$.

Perform a backward error analysis that fully accounts for round-off error in the numerical computation in order to prove that the exact first positive eigenvalues of the matrices $H\Xi \alpha $ must also be bounded below by $810$ at each $\alpha \u2208G$.

Use perturbation theory to bound the exact first positive eigenvalue of $H\Xi \alpha $ below by $34$ over the whole interval of

*α*values between 0 and $710$.

When discussing round-off error due to working in the floating-point arithmetic, we will denote “machine epsilon” by *ϵ*. The significance of this number is that we will assume that all complex numbers *a* can be represented by floating-point numbers $a\u0303$ such that $|a\u2212a\u0303|\u2264\u03f5a$. We will also make the standard assumption about creation of round-off error in the floating-point arithmetic operations: if $a\u0303$ and $b\u0303$ are floating-point complex numbers and if $(a\u0303Ob\u0303)comp$ and $a\u0303Ob\u0303$ represent the numerically computed value and exact value of an arithmetic operation on the numbers $a\u0303$ and $b\u0303$, then $(a\u0303Ob\u0303)comp=a\u0303Ob\u0303+e$, where $|e|\u2264(a\u0303Ob\u0303)\u03f5$. In Python, this is indeed the case, for all reasonably sized (such that stack overflow does not occur) complex numbers, with *ϵ* = 2.220 44 × 10^{−16} (5sf). We now present the main points of parts 2–4 of the strategy, postponing proofs of intermediate lemmas to Appendix E.

For part 2 of the strategy, for each $\alpha \u2208G$, we let $H\u0303\Xi \alpha $ denote $H\Xi \alpha $ (which is known exactly) evaluated as floating-point numbers. We generate numerically computed eigenpairs $\lambda \u0303j,v\u0303j$ for 1 ≤ *j* ≤ 81 for each $H\u0303\Xi \alpha $ using numpy’s Hermitian eigensolver eigh. We find that the smallest first positive eigenvalue of $H\u0303\Xi \alpha $ for $\alpha \u2208G$ is 0.814 719 126 144 543 6 (computed using compute_PHalphaP_enclosures.py in the Github repo). Note that the difference between this number and $810$ is bounded below by 0.01.

The main tool for part 3 of the strategy is the following theorem.

*Let*

*m*

*and*

*n*

*denote positive integers with*

*m*≤

*n*

*. Let*

*A*

*be a Hermitian*

*n*×

*n*

*matrix, and let*${vj}1\u2264j\u2264m$

*be orthonormal*

*n*

*-vectors satisfying*(

*A*−

*λ*

_{j}

*I*)

*v*

_{j}=

*r*

_{j}

*for scalars*

*λ*

_{j}

*and*

*n*

*-vectors*

*r*

_{j}

*for each*1 ≤

*j*≤

*m*

*. Then, there are*

*m*eigenvalues ${\alpha j}1\u2264j\u2264m$

*of*

*A*,

*which can be put into one-to-one correspondence with*

*λ*

_{j}’s

*such that*

See Appendix E 1.□

Naïvely, one would hope to be able to calculate enclosure intervals for every eigenvalue of $H\Xi \alpha $ and, in particular, a lower bound on the first positive eigenvalue of $H\Xi \alpha $ by directly applying Theorem IV.2 with $A=H\Xi \alpha $, *m* = 81, and *λ*_{j} and *v*_{j} given for each 1 ≤ *j* ≤ 81 by the approximate eigenpairs $\lambda \u0303j,v\u0303j$ computed in part 2. However, we cannot directly apply the theorem because ${v\u0303j}1\u2264j\u226481$ are not exactly orthonormal because of round-off error. Hence, we will prove the existence of an exactly orthonormal set ${v\u0302j}1\u2264j\u226481$ close to the set ${v\u0303j}1\u2264j\u226481$ and apply Theorem IV.2 to the set ${v\u0302j}1\u2264j\u226481$ (with the same $\lambda \u0303j$) instead. Note that to carry out this strategy, we must bound the residuals $r\u0302j\u2254(H\Xi \alpha \u2212\lambda \u0303j)v\u0302j$. The result we need to implement this strategy is the following. Note that the result requires numerical computation of inner products and residuals, and we account for round-off error in these computations.

*Let*

*m*

*and*

*n*

*be positive integers with*

*m*≤

*n*

*. Let*

*A*

*be an*

*n*×

*n*

*Hermitian matrix, let*$A\u0303$

*denote*

*A*

*evaluated in floating-point numbers, and let*$v\u0303j,\lambda \u0303j$

*for*1 ≤

*j*≤

*m*

*be a set of*

*n*

*-dimensional vectors and real numbers, respectively. Let*$v\u0303iv\u0303jcomp$

*denote their numerically computed inner products, and let*$r\u0303j,comp\u2254A\u0303\u2212\lambda \u0303jIv\u0303jcomp$

*denote their numerically computed residuals. Let*

*ϵ*

*denote machine epsilon, and assume*

*nϵ*< 0.01

*. Let*

*μ*

*be*

*Then, as long as*$m\mu <12$

*, there is an orthonormal set of*

*n*

*-vectors*${v\u0302j}1\u2264j\u2264m$

*whose residuals*$r\u0302j\u2254(A\u2212\lambda \u0303jI)v\u0302j$

*satisfy the bound*

*where*‖

*A*‖

_{max}

*denotes the largest of the absolute values of the elements of the matrix*

*A*.

See Appendix E 2.□

Numerical computation (using the script compute_PHalphaP_enclosures.py in the Github repo) shows that the maximum of $sup1\u2264i\u2264m|v\u0303iv\u0303icomp\u22121|$ and $supi\u2260j1\u2264i,j\u2264m|v\u0303iv\u0303jcomp|$ over $\alpha \u2208G$ is bounded by 7 × 10^{−15}. Hence, we can apply Theorem IV.3 with $A=H\Xi \alpha $ and $\lambda \u0303j,v\u0303j$ given by the numerically computed eigenpairs of $H\u0303\Xi \alpha $ to obtain orthonormal sets ${v\u0302j}1\u2264j\u226481$ whose residuals with respect to $H\Xi \alpha $ satisfy (35). The following is straightforward.

The first estimate follows from ‖*P*_{Ξ}*H*^{0}*P*_{Ξ}‖ ≤ 7 and ‖*H*^{1}‖ ≤ 3. The second estimate follows immediately from writing the matrix $H\Xi \alpha $ in the chiral basis.□

We can now apply Theorem IV.2 with $A=H\Xi \alpha $ and *λ*_{j}, *v*_{j} given by the numerically computed $\lambda \u0303j$ from part 2 and $v\u0302j$ coming from Theorem IV.3, in order to derive rigorous enclosure intervals for every eigenvalue of $H\Xi \alpha $. We find that (using the script compute_PHalphaP_enclosures.py in the Github repo) the suprema over $\alpha \u2208G$ of $sup1\u2264j\u2264m\Vert v\u0303j\Vert \u221e$, $sup1\u2264j\u2264m\Vert r\u0303j,comp\Vert \u221e,\Vert H\u0303\Xi \alpha \Vert max$, and $sup1\u2264j\u2264m|\lambda \u0303j|$ are bounded by 1, 5 × 10^{−14}, 7, and 8, respectively. It is then easy to see that 2 × 81 times the right-hand side of (35) is much smaller than 0.01 and is hence smaller than the distance between the minimum over $\alpha \u2208G$ of the numerically computed first positive eigenvalues of $H\u0303\Xi \alpha $ and $810$. We can therefore conclude that the first positive eigenvalues of $H\Xi \alpha $ are bounded below by $810$ at every $\alpha \u2208G$.

The main tool for part 4 of the strategy is the following.

*Let*

*A*

^{α}

*be an*

*n*×

*n*

*Hermitian matrix real-analytically depending on a real parameter*

*α*

*. Denote the ordered eigenvalues of*

*A*

^{α}

*by*$\lambda j\alpha $

*for*1 ≤

*j*≤

*n*

*. Then, for any*

*α*

*and*

*α*

_{0},

See Appendix E 3.□

We would like to apply Theorem IV.4 to bound the variation of eigenvalues of $H\Xi \alpha $. To this end, we require the following proposition, which bounds the derivative of $H\Xi \alpha $ with respect to *α* over the interval $0\u2264\alpha \u2264710$.

See Appendix E 4.□

Proposition IV.7 combined with Theorem IV.4 explains the choice of distance $h=1388831$ between grid points. Assuming that the first positive eigenvalue of $H\Xi \alpha $ is bounded below by $810$ at grid points between 0 and $710$ separated by *h*, we see that as long as

Proposition IV.7 and Theorem IV.4 guarantee that the first eigenvalue of $H\Xi \alpha $ must be greater than $34$ over the whole interval $0\u2264\alpha \u2264710$.□

*In the Proof of Proposition IV.5, we bound the round-off error, which can occur in our numerical computations in order to draw rigorous conclusions. A common approach to this is the interval arithmetic; see Rump*^{17} *and references therein. Our approach applies directly to the present problem and is just as rigorous.*

The results of a computation of the eigenvalues of $H\Xi \alpha $ are shown in Fig. 3.

Assuming Proposition IV.4 and Proposition IV.5, the bound (31) becomes, for all $0\u2264\alpha \u2264710$,

We now assume the following, proved in Appendix C.

$\Vert H1\Psi 8\Vert \u2264320$*.*

We can now give the Proof of Theorem IV.1.

The proof follows immediately from Lemma IV.1, Proposition IV.4, Assumption IV.5, and Proposition IV.8.□

## SUPPLEMENTARY MATERIAL

In the supplementary material, we list the chiral basis functions, which span the space Ξ, list terms Ψ^{5} − Ψ^{8} in the formal expansion of *ψ*^{α}, and derive the TKV Hamiltonian from the Bistritzer–MacDonald model.

## ACKNOWLEDGMENTS

M.L. and A.B.W. were supported, in part, by ARO MURI Award No. W911NF-14-0247 and NSF DMREF Award No. 1922165.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

### APPENDIX A: DERIVATION OF EXPRESSION FOR FERMI VELOCITY IN TERMS OF $LK,1,12$ ZERO MODE OF *H*^{α}

The Bloch eigenvalue problem for the TKV Hamiltonian at quasi-momentum ** k** is

where *H*^{α} is as in (1) and

By Propositions III.6 and III.7, 0 is a twofold (at least) degenerate eigenvalue at the moiré *K* point ** k** = 0, with associated eigenfunctions $\psi \xb11\alpha $ as in Proposition III.7. In what follows, we assume that 0 is

*exactly*twofold degenerate so that $\psi \xb11\alpha $ form a basis of the degenerate eigenspace. This assumption is clearly true for small

*α*but could, in principle, be violated for

*α*> 0.

Introducing $\chi k\alpha \u2254e\u2212ik\u22c5r\psi k\alpha $, we derive the equivalent Bloch eigenvalue problem with ** k**-independent boundary conditions

where

where *D*_{x,y} ≔ −*i∂*_{x,y} and

Clearly, $\psi \xb11\alpha $ remain a basis of the zero eigenspace for the problem (A1) at ** k** = 0.

Differentiating the operator $Dk\alpha $, we find $\u2202kxDk\alpha =I2$ and $\u2202kyDk\alpha =iI2$, where *I*_{2} denotes the 2 × 2 identity matrix, so that

By degenerate perturbation theory,^{18} for small ** k**, we have that eigenfunctions $\chi k\alpha $ of (A1) are given by

where the coefficients *c*_{σ,k} and associated eigenvalues *E*_{k} ≈ *ϵ*_{k} are found by solving the matrix eigenvalue problem

Using (A2) and the explicit forms of $\psi \xb11\alpha $ given by Proposition III.7, we find that the matrix on the left-hand side of (A3) can be simplified to

It follows that, for small ** k**, we have

*E*

_{k}≈ ±

*v*(

*α*)|

**|, where**

*k**v*(

*α*) = |

*λ*(

*α*)| is as in (23).

### APPENDIX B: THE CHIRAL BASIS OF $LK,12$ AND ACTION OF *H*^{0} AND *H*^{1} WITH RESPECT TO THIS BASIS

#### 1. The spectrum and eigenfunctions of *H*^{0} in $LK2$

The first task is to understand the spectrum and eigenfunctions of *H*^{0} in $LK2$. In the Appendix B 2, we will discuss the spectrum and eigenfunctions of *H*^{0} in $LK,12$. Recall that

where $\u2202\u0304=12(\u2202x+i\u2202y)$. To describe the eigenfunctions of *H*^{0} in $LK2$, we introduce some notation. Let $v=v1,v2$ be a vector in $R2$. Then, we will write

Finally, let *V* denote the area of the moiré cell Ω.

*The zero eigenspace of*

*H*

^{0}

*in*$LK2$

*is spanned by*

*For all*

**≠ 0**

*G**in the reciprocal lattice, then*

*are eigenfunctions with eigenvalues*±|

**|**

*G**. For all*

*G**in the reciprocal lattice,*

*are eigenfunctions with eigenvalues*±|

*q*_{1}+

**|**

*G**. The operator*

*H*

^{0}

*has no other eigenfunctions in*$LK2$

*other than linear combinations of these, and hence, the spectrum of*

*H*

^{0}

*in*$LK2$

*is*

The proof is a straightforward calculation taking into account the $LK2$ boundary conditions given by (20) with ** k** = 0. For example,

*e*

_{2}and

*e*

_{4}are zero eigenfunctions of

*H*

^{0}but in $LK\u20322$, not $LK2$.□

Note that (as it must be because of the chiral symmetry) the spectrum is symmetric about 0 and all of the eigenfunctions with negative eigenvalues are given by applying $S$ to the eigenfunctions with positive eigenvalues.

The union of the lattices Λ^{*} and Λ^{*} + *q*_{1} has the form of a honeycomb lattice in momentum space, where the lattice Λ^{*} corresponds to “A” sites and Λ^{*} + *q*_{1} corresponds to “B” sites (or vice versa); see Fig. 4.

#### 2. The spectrum and eigenfunctions of *H*^{0} in $LK,12$

We now discuss the spectrum of *H*^{0} in $LK,12$.

*The zero eigenspace of*

*H*

^{0}

*in*$LK,12$

*is spanned by*

*For all*

**≠ 0**

*G**in the reciprocal lattice*Λ

^{*}

*,*

*are eigenfunctions of*

*H*

^{0}

*in*$LK,12$

*with associated eigenvalues*±|

**|**

*G**. For all*

*G**in the reciprocal lattice*Λ

^{*}

*,*

*are eigenfunctions of*

*H*

^{0}

*in*$LK,12$

*with associated eigenvalues*±|

*q*_{1}+

**|**

*G**. The operator*

*H*

^{0}

*has no other eigenfunctions in*$LK,12$

*other than linear combinations of these, and hence, the spectrum of*

*H*

^{0}

*in*$LK,12$

*is*

The proof is another straightforward calculation starting from Proposition B.1.□

For an illustration of the support of the $LK,12$-eigenfunctions of *H*^{0} on the momentum space lattice, see Fig. 5. It is important to note that the notation introduced in Proposition B.2 is not one-to-one because, for example,

for any ** G** ≠ 0 in Λ

^{*}.

#### 3. The chiral basis of $LK,12$

Recall that zero modes of *H*^{α} can be assumed to be eigenfunctions of the chiral symmetry operator $S$. It follows that the most convenient basis for our purposes is not be the spectral basis just introduced but the basis of $LK,12$ consisting of eigenfunctions of $S$. We call this basis the chiral basis.

*The chiral basis of*$LK,12$

*is defined as the union of the functions*

*and*

The following is straightforward.

*The chiral basis is an orthonormal basis of* $LK,12$*. The modes* $\chi 0\u0303$*,* $\chi G\u0303,1$*, and* $\chi q1+G\u0303,1$ *are* +1 *eigenfunctions of* $S$*, while the modes* $\chi G\u0303,\u22121$ *and* $\chi q1+G\u0303,\u22121$ *are* −1 *eigenfunctions of* $S$*.*

Written out, chiral basis functions have a very simple form. We have

and for all ** G** ∈ Λ

^{*}\{

**0**},

and for all ** G** ∈ Λ

^{*},

We use the chiral basis to divide up $LK,12$ as follows.

*We define spaces* $LK,1,\xb112$ *to be the spans of the* ±1 *eigenfunctions of* $S$ *in* $LK,12$*, respectively.*

Clearly, we have

We can divide up the chiral basis more finely as follows.

*We define*

*Note that notations* *A* *and* *B* *in Definition B.3 refer to* *A* *and* *B* *sites of the momentum space lattice, not to the* *A* *and* *B* *sites of the real space lattice. Recalling Remark III.1 and comparing* *(B2)* *and* *(B3)* *with* *(2)**, we see that* $LK,1,1,A2$ *corresponds to wave-functions supported on* *A* *sites of layer 1,* $LK,1,1,B2$ *corresponds to wave-functions supported on* *A* *sites of layer 2,* $LK,1,\u22121,A2$ *corresponds to wave-functions supported on* *B* *sites of layer 1, and* $LK,1,\u22121,B2$ *corresponds to wave-functions supported on* *B* *sites of layer 2.*

Clearly, we have

The following propositions are straightforward to prove. For the first claim, note that ${S,H0}=0$.

*Operator*

*H*

^{0}

*maps*$LK,1,\xb11,\sigma 2\u2192LK,1,\u22131,\sigma 2$

*for*

*σ*=

*A*,

*B*

*. The action of*

*H*

^{0}

*on chiral basis functions is as follows:*

*for all*

**∈ Λ**

*G*^{*}

*with*

**≠ 0,**

*G**and for all*

**∈ Λ**

*G*^{*},

*Let*

*P*

*denote the projection operator onto*$\chi 0\u0303$

*in*$LK,12$

*, and*

*P*

^{⊥}= 1 −

*P*

*. Then, the operator*$P\u22a5(H0)\u22121P\u22a5$

*maps*$LK,1,\xb11,\sigma 2\u2192LK,1,\u22131,\sigma 2$

*for*

*σ*=

*A*,

*B*

*, and*

*for all*

**∈ Λ**

*G*^{*}

*with*

**≠ 0**

*G**, and*

*for all*

**∈ Λ**

*G*^{*}

*.*

In the Appendixes B 4 and 5, we will study the action of the operator *H*^{1} on $LK,12$ with respect to the chiral basis.

#### 4. The spectrum of *H*^{1} in $LK2$ and $LK,12$

Recall that

where $U(r)=e\u2212iq1\u22c5r+ei\varphi e\u2212iq2\u22c5r+e\u2212i\varphi e\u2212iq3\u22c5r$. We claim the following.

*For each*

*r*_{0}∈ Ω

*,*±|

*U*(

*r*_{0})|

*and*±|

*U*(−

*r*_{0})|

*are eigenvalues of*$H1:LK2\u2192LK2$

*. For*

*r*_{0}

*such that*

*U*(

*r*_{0}) ≠ 0

*, the*±|

*U*(

*r*_{0})|

*eigenvectors are*

*For*

*r*_{0}

*such that*

*U*(−

*r*_{0}) ≠ 0

*, the*±|

*U*(−

*r*_{0})|

*eigenvectors are*

*When*

*U*(

*r*_{0}) = 0

*, zero is a degenerate eigenvalue with associated eigenfunctions*

*e*

_{2}

*δ*(

**−**

*r*

*r*_{0})

*and*

*e*

_{3}

*δ*(

**−**

*r*

*r*_{0})

*. When*

*U*(−

*r*_{0}) = 0

*, zero is a degenerate eigenvalue with associated eigenfunctions*

*e*

_{1}

*δ*(

**−**

*r*

*r*_{0})

*and*

*e*

_{4}

*δ*(

**−**

*r*

*r*_{0})

*. Finally,*

*H*

^{1}must be contained in the interval [−3, 3]. To see that the spectrum actually equals [−3, 3], note that if $r0=4\pi 33,0$; then,

*U*(

*r*_{0}) = 3. On the other hand, when

*r*_{0}= 0, we have

*U*(

*r*_{0}) = 0 so that the spectrum of

*H*

^{1}in $LK2$ equals [−3, 3].□

By taking linear combinations of rotated copies of the *H*^{1} eigenfunctions, just as we did with the *H*^{0} eigenfunctions, it is straightforward to prove an analogous result to Proposition B.6 in $LK,12$. We record only the following.

#### 5. The action of *H*^{1} on $LK,12$ with respect to the chiral basis

We now want to study the action of *H*^{1} on $LK,12$ with respect to the chiral basis. We will prove two propositions, which parallel Proposition B.4.

*The operator*

*H*

^{1}

*maps*$LK,1,1,A2\u2192LK,1,\u22121,B2$

*and*$LK,1,1,B2\u2192LK,1,\u22121,A2$

*. The action of*

*H*

^{1}

*on chiral basis functions is as*

*follows:*

*and*

*For all*

**∈ Λ**

*G*^{*}\{

**0**}

*,*

*For all*

**∈ Λ**

*G*^{*}\{

**0**}

*,*

Note that *H*^{1} exchanges chirality ($S$ eigenvalue) and the *A* and *B* momentum space sublattices, while *H*^{0} only exchanges chirality. Proposition B.8 has a simple interpretation in terms of nearest neighbor hopping in the momentum space lattice; see Figs. 6 and 7.

*At first glance,*

*(B5) and (B6*)

*appear different from*

*(B7)*

*and*

*(B8)*

*because they appear to violate*$2\pi 3$

*rotation symmetry. However, this is not the case since every chiral basis function individually respects this symmetry. For example, using*$\chi q1\u0303,\u22121=\chi q2\u0303,\u22121=\chi q3\u0303,\u22121$

*and*$z\u0302q1\u0304=ei\varphi z\u0302q2\u0304=e\u2212i\varphi z\u0302q3\u0304$

*, we can re-write*

*(B5)*

*in a way that manifestly respects the*$2\pi 3$

*rotation symmetry as*

*Equation (B6)*

*can also be written in a manifestly rotationally invariant way, but the expression is long, and hence, we omit it. Note that*

*(B6)*

*cannot have a term proportional to*$\chi 0\u0303$

*since*$\chi 0\u0303\u2208LK,1,12$

*and*

*H*

^{1}

*maps*$LK,1,12\u2192LK,1,\u221212$

*.*

*The operator*

*H*

^{1}

*maps*$LK,1,\u22121,A2\u2192LK,1,1,B2$

*and*$LK,1,\u22121,B2\u2192LK,1,1,A2$

*. The action of*

*H*

^{1}

*on chiral basis functions is as follows:*

*For all*

**∈ Λ**

*G*^{*}\{

**0**}

*,*

*For all*

**∈ Λ**

*G*^{*}

*,*

The proof is similar to that of Proposition B.8 and is hence omitted.□

### APPENDIX C: FORMAL EXPANSION OF THE ZERO MODE

We now bring to bear the developments of the Appendix B on the asymptotic expansion of the zero mode $\psi \alpha \u2208LK,1,12$ starting from $\Psi 0=e1=\chi 0\u0303$. We first give the Proof of Proposition IV.1.

*H*

^{0}. The general solution of

*H*

^{0}Ψ

^{1}= −

*H*

^{1}Ψ

^{0}is

*C*is an arbitrary constant, which is in $LK,1,12$ by Proposition B.4. To ensure that Ψ

^{1}is orthogonal to Ψ

^{0}, we take

*C*= 0. It is clear that this procedure can be repeated to derive an expansion to all orders satisfying the conditions of Proposition IV.1.□

Our goal is to calculate $\Psi n\u2208LK,1,12$ satisfying the conditions of Proposition B.4 up to *n* = 8. This amounts to calculating, for *n* = 1 to *n* = 8,

We do this algorithmically by repeated application of the following proposition, which combines Proposition B.8 and Proposition B.5.

*The operator*$\u2212P\u22a5(H0)\u22121P\u22a5H1$

*maps*$LK,1,1,A2\u2192LK,1,1,B2$

*and*$LK,1,1,B2\u2192LK,1,1,A2$

*. Its action on chiral basis functions is as follows:*

*and*

*For all*

**∈ Λ**

*G*^{*}\{

**0**}

*,*

*For all*

**∈ Λ**

*G*^{*}\{

**0**}

*,*

We now claim the following.

*Let*Ψ

^{n}

*be the sequence defined by Proposition IV.1. Then,*

*q*_{2}=

*q*_{1}+

*b*_{1}and

*q*_{3}=

*q*_{1}+

*b*_{2}. The derivation of Eq. (C6) is more involved, so we give details. Using linearity and applying (C3) twice, we find that

*R*

_{ϕ}(

*q*_{1}+

*b*_{1}−

*b*_{2}) =

*q*_{1}+

*b*_{2}−

*b*_{1}, we have $\chi q1+b1\u2212b2\u0303,1=\chi q1+b2\u2212b1\u0303,1$. These terms also cancel, leaving (C6). The derivation of (C7) (and the higher corrections) is involved but does not depend on any new ideas and is therefore omitted.□

We give the explicit forms of Ψ^{5}-Ψ^{8} in the supplementary material.

*Written out,*

*(C4)*

*and*

*(C5)*

*become*

*and*

*which agree with*

*Eq. (24)*

*of Tarnopolsky et al.*

^{4}*up to an overall factor of*$V$

*(this factor cancels in the Fermi velocity, so there is no discrepancy).*

Using orthonormality of the chiral basis functions, it is straightforward to calculate the norms of each of the Ψ^{n}. We have the following.

*Note that the sequence of norms of the expansion functions grows much slower than the pessimistic bound* ‖Ψ^{N+1}‖ ≤ 3‖Ψ^{N}‖, *N* = 0, 1, 2, … *guaranteed by Proposition IV.2. The reason is that bounds* *(27)* *and* *(28)* *are never attained. As* *N* *becomes larger, the bound* *(27)* *is very pessimistic because* Ψ^{N} *is mostly made up of eigenfunctions of* *H*^{0} *with eigenvalues strictly larger than 1. The bound* *(28)* *is also very pessimistic because it is attained only at delta functions, which can only be approximated with a superposition of a large number of eigenfunctions of* *H*^{0}*. It seems possible that a sharper bound could be proved starting from these observations, but we do not pursue this in this work.*

We finally give the Proof of Proposition IV.8.

### APPENDIX D: PROOF OF PROPOSITION IV.4

We choose Ξ as

Part 1 of Proposition IV.4 follows immediately from observing that $\chi q1\u22122b1\u22122b2\u0303,\xb11$ is not in Ξ, but |*q*_{1} − 2*b*_{1} − 2*b*_{2}| = 7. That *μ* = 7 is optimal can be seen from Fig. 8.

Part 2 follows from the fact that *ψ*^{8,α} depends only on eigenfunctions of *H*^{0} with eigenvalues with magnitude less than or equal to $43$. The largest eigenvalue is $43$, coming from dependence of Ψ^{8} on $\chi \u22124b2\u0303,1$, since $|\u22124b2|=43$.

Part 3 can be seen from Fig. 8.

### APPENDIX E: PROOF OF PROPOSITION IV.5

#### 1. Proof of Theorem IV.2

We will prove Theorem IV.2 starting from Theorem 11.5.1 of Parlett,^{19} where the proof can be found.

*Let*

*Q*

*be an*

*n*×

*m*

*matrix, which satisfies*

*Q*

^{†}

*Q*=

*I*

_{m}

*. Define*

*H*=

*Q*

^{†}

*AQ*

*and*

*R*=

*AQ*−

*QH*

*. Let*${\theta j}1\u2264j\u2264m$

*denote the eigenvalues of*

*H*

*(the Ritz values). Then,*

*m*

*of*

*A*

*’s eigenvalues*${\alpha j}1\u2264j\u2264m$

*can be put into one-to-one correspondence with the*${\theta j}1\u2264j\u2264m$

*in such a way*

*that*

*Q*be the matrix whose columns are

*v*

_{1}, …,

*v*

_{m}. Using orthonormality of the

*v*

_{j},

*Q*

^{†}

*Q*=

*I*

_{m}and

*H*, denoted by

*θ*

_{j}, are close to

*λ*

_{j}’s. By the Gershgorin circle theorem, we have

*v*

_{j}‖

_{2}= 1,

*λ*

_{j}s and exact eigenvalues

*α*

_{j},

*R*≔

*AQ*−

*QH*= (

*I*−

^{†})

*AQ*. Since

^{†}projects onto

*v*

_{j},

*R*simplifies to

^{†}is a projection, so is

*I*−

^{†}, and hence, ‖

*R*‖

_{2}≤ ‖

*R*′‖

_{2}. To bound ‖

*R*′‖

_{2}, note that for any

*v*with ‖

*v*‖

_{2}= 1, we have

*e*

_{j}denote the standard orthonormal basis vectors. The result now follows.□

#### 2. Proof of Theorem IV.3

We start with the following lemma, which guarantees that numerically computed approximately orthonormal sets can be approximated by exactly orthonormal sets.

*Let*$v\u03031,\u2026,v\u0303m$

*be*

*n*

*-dimensional vectors, let*$v\u0303iv\u0303jcomp$

*for*1 ≤

*i*,

*j*≤

*m*

*denote their numerically computed inner products, let*

*ϵ*

*denote machine epsilon, and assume that*

*nϵ*< 0.01

*. Define*

*Then, as long as*$m\mu <12$,

*there is a set of*

*n*-

*dimensional orthonormal vectors*$v\u03021,\u2026,v\u0302m$,

*which satisfy*

^{20}) and assuming that

*nϵ*< 0.01, we have that for each 1 ≤

*i*,

*j*≤

*m*, $v\u0303iv\u0303j=v\u0303iv\u0303jcomp+eij$, where $|eij|\u2264(1.01)n\u03f5|v\u0303i|\u22a4|v\u0303j|\u2264(1.01)n\u03f5\Vert v\u0303i\Vert 2\Vert v\u0303j\Vert 2$. Letting $Q\u0303$ denote the matrix whose columns are $v\u0303i$’s, then $Q\u0303\u2020Q\u0303\u2212Im=E$, where, for all

*i*≠

*j*, $|Eij|\u2264|v\u0303iv\u0303jcomp|+(1.01)n\u03f5\Vert v\u0303i\Vert 2\Vert v\u0303j\Vert 2$, and for all

*i*, $|Eii|\u2264|1\u2212v\u0303iv\u0303icomp|+(1.01)n\u03f5\Vert v\u0303i\Vert 22$. Paying the price of factors of $n$ to replace ‖·‖

_{2}norms by ‖·‖

_{∞}norms, we can obtain a trivial bound on the maximal element of

*E*, denoted as ‖

*E*‖

_{max}, by

*μ*in the statement of the theorem. Using the Gershgorin circle theorem, we then have that the eigenvalues

*λ*of $Q\u0303\u2020Q\u0303$ satisfy |

*λ*− 1| ≤

*m*‖

*E*‖

_{max}. We claim that there are exact orthonormal vectors $v\u0302j$ near (in the ‖·‖

_{2}-norm) to the $v\u0303j$. To see this, note that $Q\u0302\u2254Q\u0303(Q\u0303\u2020Q\u0303)\u22121/2$ satisfies $Q\u0302\u2020Q\u0302=Im$ and

*λ*

_{max}and

*λ*

_{min}denote the maximum and minimum eigenvalues of $Q\u0303\u2020Q\u0303$, respectively. Then, $(Q\u0303\u2020Q\u0303)1/2\u2212Im2\u2264max|\lambda min1/2\u22121|,|\lambda max1/2\u22121|$. Since

*λ*

_{min}is bounded below by 1 −

*m*‖

*E*‖

_{max}and

*λ*

_{max}is bounded above by 1 +

*m*‖

*E*‖

_{max}, we have

*x*)

^{1/2}− 1| ≤ 2

^{−1/2}|

*x*| and |(1 −

*x*)

^{1/2}− 1| ≤ 2

^{−1/2}|

*x*|. Since by assumption $m\Vert E\Vert max<12$, we conclude

*j*≤

*m*, the result is proved.□

Using Lemma E.2, we have that there exists an exactly orthonormal set ${v\u0302j}1\u2264j\u2264m$ nearby to the set ${v\u0303j}1\u2264j\u2264m$. We now want to bound the residuals of the $v\u0302j$ in terms of numerically computable quantities. We start with the following easy lemma whose proof is a straightforward manipulation.

*Let*

*A*

*be an*

*n*×

*n*

*Hermitian matrix, and suppose that*$r\u0302\u2254(A\u2212\lambda \u0303I)v\u0302$

*and*$r\u0303\u2254(A\u2212\lambda \u0303I)v\u0303$

*. Then,*

The following lemma quantifies the error in approximating exact residuals by numerically computed values.

*Let*

*A*

*be a Hermitian*

*n*×

*n*

*matrix, and let*$A\u0303$

*denote the matrix whose entries are those of*

*A*

*evaluated as floating-point numbers. Let*$A\u0303\u2212\lambda \u0303Iv\u0303comp$

*denote the numerically computed value of*$A\u0303\u2212\lambda \u0303Iv\u0303$

*in the floating-point arithmetic. Then,*$r\u0303\u2254(A\u2212\lambda \u0303I)v\u0303$

*satisfies*

*A*and

*B*with entries

*A*

_{ij}and

*B*

_{ij}, we will write |

*A*| to denote the matrix with entries |

*A*

_{ij}| for all

*i*,

*j*, and |

*A*| ≤ |

*B*| if |

*A*

_{ij}| ≤ |

*B*

_{ij}| for all

*i*,

*j*. It is straightforward to see that (see Chap. 2.7 of Golub and Van Loan

^{20}) $A\u0303=A+F$, where |

*F*| <

*ϵ*|

*A*|. In addition, $(A\u0303\u2212\lambda \u0303I)v\u0303=(A\u0303\u2212\lambda \u0303I)v\u0303comp+g$, where $|g|\u2264(1.01)n\u03f5|(A\u0303\u2212\lambda \u0303I)||v\u0303|$. Now, note that $(A\u2212\lambda \u0303I)v\u0303=(A\u0303\u2212\lambda \u0303I)v\u0303\u2212Fv\u0303$ so that

We now prove estimate (35). Applying Lemma E.2 to the set ${v\u0303j}1\u2264j\u2264m$ yields an orthonormal set ${v\u0302j}1\u2264j\u2264m$ such that $\Vert v\u0302j\u2212v\u0303j\Vert 2\u22642\u22121/2m\mu $, where *μ* is as in (E1). By Lemma E.3, we have that

The estimate now follows easily upon applying Lemma E.4 and taking the sup over *j*.□

#### 3. Proof of Theorem IV.4

*U*denotes a subspace of $Cn$),

*v*≠ 0, we have by Taylor’s theorem

#### 4. Proof of Proposition IV.7

*α*< 1, we have ‖

*P*

_{Ξ}

*H*

^{α}

*P*

_{Ξ}‖

_{2}≤ 10, and ‖

*H*

^{1}‖

_{2}≤ 3. It remains only to estimate ‖

*∂*

_{α}

*Q*

^{α}‖

_{2}. Using the Dirac notation to represent $LK,12$-projections, we have

*α*≤ 1,

### APPENDIX F: PROOF OF PROPOSITION II.1

We can now prove Proposition II.1. We start by proving (11).

#### 1. Proof of (11)

We now prove (11). It is straightforward to derive

We now make two observations which simplify the computation. First, recall that the operator $\u2212P\u22a5(H0)\u22121P\u22a5H1$ maps $LK,1,1,A2\u2192LK,1,1,B2$ and $LK,1,1,B2\u2192LK,1,1,A2$. It follows that $\Psi 0\u2208LK,1,1,A2$, $\Psi 1\u2208LK,1,1,B2$, $\Psi 2\u2208LK,1,1,A2$, and so on, and hence,

It follows that all terms in (F1) with odd powers of *α* vanish. Second, note that since Ψ^{0} ∈ ran*P*, while Ψ^{n} ∈ ran*P*^{⊥} for all *n* ≥ 1, we have that

Deriving (11) is then just a matter of computation using the properties of the chiral basis. For the leading term, we have

For the *α*^{2} term, the only non-zero term is

using (C4). For the *α*^{4} term, the possible non-zero terms are

but Ψ^{3} and Ψ^{1} depend on orthogonal chiral basis vectors [see (C4) and (C6)], so we are left with

using (C5) and orthgonality of $\chi \u2212b1\u0303,1$ and $\chi \u2212b2\u0303,1$. We omit the derivation of the higher terms since the derivations do not require any new ideas.

#### 2. Proof of (10)

It is straightforward to derive

We now note the following.

*Let* *χ* *be a chiral basis function in* $LK,1,12$*. Then,* *χ*^{*}(−** r**) =

*χ*(

**)**

*r**.*

Using Proposition F.1 and the same two observations as in the Appendix F 1, we have that the only non-zero terms in (F2) are those with even powers of *α* and that other than the leading term, terms involving Ψ^{0} do not contribute. The calculation is then similar to the previous case. For the leading order term, we have

The only non-zero *α*^{2} term is

The only non-zero *α*^{4} term is

We omit the derivation of the higher terms since the derivations do not require any new ideas.

Proposition IV.1 implies that the series expansion of *ψ*^{α} exists up to any order. We can therefore define formal infinite series by

We then have the following.

The series agree exactly without any simplifications up to terms of *α*^{9}. However, because the even and odd terms in the expansion of *ψ*^{α} are orthogonal (since they lie in $LK,1,1,A2$ and $LK,1,1,B2$, respectively), all terms with odd powers of *α* vanish in expansions (F3) and (F4). The series may disagree at order *α*^{10} because the infinite series includes terms arising from inner products of Ψ^{1} and Ψ^{9}.□