Dirac points for twisted bilayer graphene with in-plane magnetic field

We study Dirac points of the chiral model of twisted bilayer graphene (TBG) with constant in-plane magnetic field. For a fixed small magnetic field, we show that as the angle of twisting varies between magic angles, the Dirac points move between $ K, K' $ points and the $ \Gamma $ point. The Dirac points for zero magnetic field and non magic angles lie at $ K $ and $ K'$, while in the presence of a non-zero magnetic field and near magic angles, they lie near the $ \Gamma $ point. For special directions of the magnetic field, we show that the Dirac points move, as the twisting angle varies, along straight lines and bifurcate orthogonally at distinguished points. At the bifurcation points, the linear dispersion relation of the merging Dirac points disappears and exhibit a quadratic band crossing point (QBCP). The results are illustrated by links to animations suggesting interesting additional structure.


Introduction
We consider the chiral model for twisted bilayer graphene in the form considered by Tarnopolsky-Kruchkov-Vishwanath [TKV19] and then studied mathematically by Becker et al [Be*22,BHZ22a,BHZ22b].Following Kwan et al [KPS20] and Qin-MacDonald [QiMa21], see also [RY13] for bilayer graphene, we introduce an additional term B = be 2πiθ corresponding to an in-plane magnetic field of strength b and direction 2πθ -see (2.1).
We concentrate on the case of simple magic α's (α is a dimensionless parameter roughly corresponding to the reciprocal of the angle of twisting of the two graphene sheets; see §3.7 for the discussion of simplicity).For the Bistritzer-MacDonald potential U BM (z) (see the caption to Figure 1) the real magic angles are expected to be simple (see Remark 1 after Theorem 2).
We have the following combination of mathematical and numerical observations: • We show (Theorem 2) that a small in-plane magnetic field destroys flat bands corresponding to simple magic α's (under an additional non-degeneracy assumption); • For small magnetic fields, the motion of Dirac points appears quasi-periodic for α ∈ [α j , α j+1 ] where α j are the magic angles for the Bistritzer-MacDonald ].In the case on the left α passes two simple magic α's; on the right, it passes two double magic α's.The Γ point corresponds to 0 and K, K ′ points to ±i.The boundary of the Brillouin zone, a fundamental domain of Λ * , is outlined in black.See https://math.berkeley.edu/~zworski/B01.mp4and https://math.berkeley.edu/~zworski/B01_double.mp4 for the corresponding animations.potential [TKV19].That is most striking for θ = 0, 2 3 for which the motion is linear -see the Remark after Theorem 1 and also Figure 6.
• Theorem 1 shows that most of the action takes place near the magic angles: the Dirac points get close to Γ point (Theorem 2; they meet there for θ = 0, Proposition 5.1 and θ = 2 3 , Proposition 5.3) at simple magic angles -see https://math.berkeley.edu/~zworski/magic_billiard.mp4 for an animation.When the Dirac cones meet, they exhibit a quadratic band crossing point (QBCP), see Figure 3 and Proposition 5.2 (its formulation requires introduction of Bloch-Floquet spectra in §3.1) -for the discussion of such phenomena in the physics literature see [dGGM12,KPS20,MLFP18].
• Figure 2 (right) shows that for fixed α's and varying directions of the magnetic field, we have "fixed points" at Γ and K, K ′ with "normal crossings" and the vertices and middle of points of edges of the boundary of the Brillouin .The magnetic field given by B = B 0 e 2πiθ with B 0 = 0.1 On the left different colours correspond to different values of θ shown in the colour bar and α varies between 0.1 and 0.9 (this is a colour map version of the left panel of Figure 1).On the right, the colours correspond to different values of α shown in the colour bar and θ varies.The predominance of green (corresponding to the range between 0.5 and 0.6) means that most of the motion happens near the (first) magic alpha -see https://math.berkeley.edu/~zworski/ first_band.mp4 for E 1 (α, k)/ max k E 1 (α, k) for fixed B as α varies.
zone.These points are precisely the intersection of the rectangles (other than Γ, K, K ′ ).
• The situation is more complicated for double (protected) magic angles: see the right panel in Figure 1: at magic α's, Dirac points are now close to K and K ′ .
This note is organized as follows: we present the Hamiltonian and the definition of Dirac points in §2.We also establish basic symmetry properties of Dirac points and a perturbation result valid away from magic α's.The next section reviews the theory of magic angles following [Be*22, BHZ22b] but in a more invariant and general way.In §4 we set up Grushin problems needed for the understanding the small in-plane magnetic fields as a perturbation.We then specialize, in §5, to directions of the magnetic field for which the Dirac points move linearly as α changes.In particular, they meet at special points and we describe the resulting quadratic band crossing.We conclude in §6 with the proofs of the main theorems.

In-plane magnetic field
Adding a constant in-plane magnetic field [KPS20,QiMa21] with magnetic vector potential A = z ⊥ B × êz ⊥ , where z ⊥ is the coordinate perpendicular to the twodimensional plane of TBG and êz ⊥ the unit vector pointing in that direction, to the chiral model of TBG [TKV19] results for layers at positions z ⊥ = ±1, in the Hamiltonian H B (α) in (2.7) build from non-normal operators where we make the following assumptions on U : In this convention the Bistritzer-MacDonald potential used in [TKV19, Be*22] corresponds to ω ℓ e i⟨z,ω ℓ K⟩ , K = 4 3 π. (2.3) For a discussion of a perpendicular constant magnetic field in the chiral model of twisted bilayer graphene we refer to [BKZ22].
Remark.We adapt here a more mathematically straightforward convention of coordinates than that of [Be*22, BHZ22a] where we followed [TKV19] (with some, possibly also misguided, small changes; our motivation comes from a cleaner agreement with theta function conventions).The translation between the two conventions is as follows: the operator considered in [Be*22], and rigorously derived in [CGG22, Wa*22] was (2.4) We then have a (twisted) periodicity with respect to 1 3 Γ and periodicity with respect to This means that to switch to (twisted) periodicity with respect to Λ we need a change of variables: The twisted periodicity condition in (2.4) corresponds to the condition in (2.2) since ωa The self-adjoint Hamiltonian built from (2.4) is given by and the Dirac points are given by the spectrum of K⟩ , e i⟨γ,K⟩ )u(x)}, with a similar definition of H 1 0 (replace L 2 loc with H 1 loc ) -see §3.1 for a systematic discussion and explanations.We recall (see §3.4) that there exists a discrete set A ⊂ C such that (2.8) The elements of A are reciprocals of magic angles and the real ones are of physical interest.As recalled in Proposition 3.3, elements of A are characterized by the condition that α −1 ∈ Spec L 2 0 T k , where C \ {K, −K} → T k is a (holomorphic) family of compact operators given in (3.25) (the spectrum is independent of k and so are its algebraic multiplicities).In this paper we will use the following notion of simplicity (see also §3.7): (2.9) Here simplicity of an eigenvalue is meant in the algebraic sense.
The first result is a consequence of simple perturbation theory and of symmetries of D B (α): and k B (α) = K + O(B).In addition, for α, B ∈ C, (2.10) Proof of Theorem 1. Proposition 3.3 shows that for α ∈ Ω the spectrum of D(α) is given by ±K + Λ * and for small B we have two eigenvalues for D B (α).The structure of D(α) implies that and since JB = −BJ we also have that is the spectrum is invariant under reflection k → −k.
Then there exists δ 0 > 0 such that for 0 < |B| < δ 0 and |α − α| < δ 0 , the spectrum of where the elements of the spectrum are included according to their (algebraic) multiplicity.In addition, for a fixed constant a 0 > 0 and for every ε there exists δ such that for 0 < |B| < δ, |α − α| < a 0 δ|B|, where we recall that elements of Λ * , in particular 0, correspond to the Γ point.
Remarks. 1. Existence of the first real magic angle α ≃ 0.585 was proved by Watson-Luskin [WaLu21] and its simplicity (including the simplicity as an eigenvalue of the operator T k defined in (3.25)) in [BHZ22a], with computer assistance in both cases.Numerically, the simplicity is valid at the computed real elements of A for the Bistritzer-MacDonald potential used in [TKV19].
2. The constant g 0 (α) can be evaluated numerically (and its non-vanishing for the first magic angles could be established via a computer assisted proof) and here are the results for the (numerically) simple magic angles for the potential U BM in Figure 1 1.Values of g 0 (α) and g 1 (α) at first magic angles.
3. The combination of Theorems 1 and 2 shows that for any U ⋐ (C \ A) ∪ {α} (with α satisfying the assumptions of Theorem 2) there exists δ = δ(U ) such that 0 From the symmetries in (2.10) we conclude that for special values of θ = 0, ± 2 3 the spectrum of D B (α) has a particularly nice structure as α varies.We state the result for θ = 0, as we can use the first identity in (2.10) to obtain the other two. (2.17) Moreover, if the assumptions of Theorem 2 are satisfied at α ∈ R then for every ε > 0 there are δ 0 , δ 1 > 0 such that (2.18) Remarks. 1.A more precise statement about the behaviour at R is given in Propositions 5.1 and 5.3 -the implicit formulas for λ = 1/α in terms of k and B describe a bifurcation phenomenon.In particular, when B is real, the bifurcation of the eigenvalues of D B (α) at 0 (at the specific value of α) is given by (5.5).For the bifurcation at the vertices of the boundary of the Brillouin zone, see (5.13).
2. The inclusion (2.17) means that the spectrum lies on a grid of straight lines parallel to the x and y axes -see https://math.berkeley.edu/~zworski/Rectangle_1.mp4.To obtain the sets of other rectangles we use the the first identity in (2.10), that is take B = ωb, b > 0.

Review of magic angle theory
We start with a general discussion of operators arising in chiral TBG models.
3.1.Bloch-Floquet theory.We recall that (The dual basis of {1, ω} is given by We then consider a generalization of (2.1): Let ρ : Λ → U (n) be a unitary representation and assume that We note that without loss of generality (amounting to a basis change on C n ) we can assume that which is a periodic operator.In view of this standard Bloch-Floquet theory applies, which can be presented using modified translations: We have Thus, we can define a generalized Bloch transform We check that and that is the inverse of B. We now define We have a unitary operator identifying where used the notation of (3.4).
In view of this, Spec L 2 0 (H k (α)) (with the domain given by H 1 0 ) is discrete and Finally, we use (3.4) and Spec L 2 (C/Λ;C) (2D z ) = Λ * (with simple eigenvalues) to see that (for ρ given by (3.2)) (3.9) 3.2.Rotational symmetries.We now introduce and in addition to (3.1) assume that (3.10) (We do not have many options here as ΩD z = ωD z Ω).Then and We have the following commutation relation A natural case to consider is given by (3.12) In the notation of (3.2), condition (3.11) means that We see that K/Λ * is the subgroup of fixed points of multiplication ω : C/Λ * → C/Λ * and it is isomorphic to Z 3 .
Since (3.11) implies that we follow [Be*22, §2.1] combine the two actions into a group of unitary action which commute with H(α): By taking a quotient by 3Λ we obtain a finite group which acts unitarily on L 2 (C/3Λ), and that action commutes with H(α): By restriction to the first two components, G and G 3 act on C n -valued function and use the same notation for those actions.
The key fact (hence the name chiral model ) is that (3.16) 3.3.Protected states.We now make the assumption (3.11) and consider the question of protected states.We are looking for the set This condition is equivalent to , where we used the notation of (3.4).Putting α = 0 we see that K 0 ⊂ K.
The following simple lemma is used a lot.To formulate it we introduce the following spaces: (with the corresponding definition of L 2 k ).
Hence if (D(α) + k)u = 0 and L γ u = u then w := e i⟨z,k⟩ u ∈ H 1 (C/3Λ; C 2n ), D(α)w = 0, and We are interested in the case of n = 2 and obtain the following reinterpretation of earlier protected states statements -see [TKV19].
Proposition 3.2.If n = 2 (in the notation of (3.2) and (3.17)) and Proof.We use (3.19) and decompose ker H 1 (C/3Λ;C 4 ) H(α) into representations of G 3 given by (3.14).From (3.16) we see that the spectrum of H(α) restricted to a representation of G 3 is symmetric with respect to the origin.If (see [Be*22, §2.2] for a review of representations of G 3 ) , (with the corresponding definition of L 2 k,p ) then the constant functions (given by the standard basis vectors in C 4 ) satisfy and since k 1 ̸ ≡ k 2 mod Λ * , all these spaces are different.The spectrum of is even (see (3.16)) and ker Existence of protected shows that we have a natural labelling for the eigenvalues of where the eigenvalues are included according to their multiplicities (and Z * := Z\{0}).
3.4.Magic angles.We recall the main result of [Be*22], the spectral characterization of magic angles.See also proof of [BHZ22b, Proposition 2.2].
Proposition 3.3.Suppose that n = 2 and that the condition (3.11) holds.Then, in the notation of Proposition 3.2 there exists a discrete set A such that where T k is a compact operator given by 3.5.Antilinear symmetry.We will make the following assumption: (3.26) A calculation based on the definition of L γ gives In particular if (as we assume)

and consequently
Since (we put α = 1 to streamline notation; that amounts to absorbing α into V ) for (3.26) to hold we need V 11 = −V 22 =: W 1 .From (3.3) we see that W 1 is Λ-periodic and there exists Λ-periodic W 0 such that (From (3.10) we see that W 1 (ωz) = ωW 1 (z) and hence the integral of W 1 over C/Λ is equal to 0; this shows that we can find W 0 , which is unique up to an additive constant.) We conclude that if we insist on (3.26) then we can, without loss of generality assume that (3.29) To verify the latter we check that, with w = (w 1 , w 2 ), Remarks. 1.The antilinear symmetry is closely related to the C 2z T symmetry in the physics literature.
We now define Then, using (3.31) and differentiating in the sense of distributions, The following Lemma is now immediate.It reinterprets the theta function argument in [TKV19].
Lemma 3.4.Suppose that p ∈ K and u ∈ ker where c(k) is given in (3.33).In particular, if u(z(k ′′ )) = 0 then (3.35) 3.7.Multiplicity one.The definition of the set of magic α's based on Proposition 3.3 does not involve the notion of multiplicity.Here we will discuss the case of multiplicity one † .One natural definition of multiplicity of magic angles is given in terms of eigenvalues of H k (α) in (3.22).We first note that We then say, that the magic angle α ∈ A is simple/has multiplicity one, if and only if As stated in (2.9) we use a stronger definition in this paper.
The operators Then for α ∈ A we have In particular, α ∈ C is a simple magic angle (in the sense of (3.37)) if and only if (3.39) † A more general discussion is presented in [BHZ23] -generic simplicity presented there is modified in view of protected multiplicity two magic angles -see the proof of Proposition 3.6.
We recall that the proof is based on Proposition 3.3 and theta function arguments reviewed in §3.6.
A symmetric choice of ρ in (3.2) is given by: This corresponds to Γ = 0 in the physics notation.In [Be*22] we followed [TKV19] and used a non-symmetric (equivalent) choice.This corresponds to the assumptions in (2.2) with k 1 = K.
Proposition 3.6.Suppose that (3.39) holds.Then, in the notation of Lemma 3.1, that is, in the notation of (3.20), u 0 ∈ L 2 k 0 ,2 .In addition, Remark.The key insight in [TKV19] was to use vanishing of u ∈ ker In [Be*22, Theorems 1] this was shown equivalent to the spectral definition based on Proposition 3.3.Here we take a direct approach: only at magic α's we have ker H 1 k 0 D(α) ̸ = {0} and (3.41) shows that its elements have to vanish at 0. Lemma 3.34 then implies vanishing of other eigenfunctions.
Proof of Proposition 3.6.From Lemma 3.1 and (3.39) we conclude that ker k 0 ,j we can decompose the kernel using these subspaces.Since D(0) is invertible with the inverse given by R(0).It then follows that (see (3.25)) (We do use ellipticity of D(α) here: the element of the kernel on L 2 must automatically be smooth.)Hence if ker This means that dim ker H 1 k 0 D(α) > 1, contradicting the simplicity assumption.The simplicity and uniqueness of the zero of u 0 (3.42) follows from [BHZ22b,Theorem 3].□ For an α ∈ A, we assume that (3.39) holds.In that case Proposition 3.6 and Lemma 3.4 show that ker Using (3.26) we see that (since A 2 = −I) Remark.From [BHZ22b, (6.6)] we see that (note the difference of notation: u(k) there is not normalized) for the basis of Λ * satisfying z(e 1 ) = 1, z(e 2 ) = ω, we have where the unitary operator τ (p) was defined in (3.7).

Grushin problems
In this section we construct Grushin problems (see [SZ07] and [TaZw23,§6]) which allow us to treat small in-plane magnetic fields as perturbations.In §5 we combine that with the spectral characterization of magic angles (Proposition 3.3) to analyze the behaviour at the Γ point and at the vertices of the boundary of the Brillouin zone.4.1.Grushin problem for D B (α). Suppose α ∈ A is simple, in the sense that (3.39) holds.We then put, in the notation of (3.43) and (3.44), We have , We now apply [TaZw23, Proposition 2.12] to obtain and, if u 0 = (ψ, φ) t , and where c(k), c * (k) > 0 come from L 2 -normalizations of u and u * . Hence, In fact G(k) is a multiple of θ(z(k)) 2 which follows from a theta function identity (see [KhZa15,(4.7a)] or [Mu83, §I.5, (A.3)]): this integral vanishes, and (4.4) gives Numerical evidence, see Table 1, suggests that for the Bistritzer-MacDonald potential and the first magic angle, (The number g 0 is determined up to phase which we can choose arbitrarily by modifying u 0 → e iθ u 0 .)Table 1 shows approximate values of |g 0 | for higher magic angles for the same potential.

4.2.
Grushin problem for the self-adjoint Hamiltonian.We now turn to the corresponding Grushin problem for H B k (α) given in (3.5) (note the irrelevant change of sign of k) where R ± (k) are the same as in (4.1).The operator H B k (α, z) is invertible for all k, |B| ≪ 1, |α − α| ≪ 1 and |z| ≪ 1.We denote the components of the inverse by E B • (k, α, z) and we have Using [TaZw23, Proposition 2.12] again we see that (in the notation of (4.6)) (Here we used the fact that E 0 − (k)E 0 − (k) * ≡ 1 and E 0 + (k) * E 0 + (k) ≡ 1 which follows from (4.2) and normalization of u(k) and u * (k).) where (under the assumption that g 0 (α) ̸ = 0) γ 0 ̸ = 0, γ 1 ̸ = 0. (The exact symmetry of signs follows from the extension of the chiral symmetry (3.16) to the Grushin problem (4.8) which shows that det

Bifurcation
This section is devoted to showing (2.18) and giving a stronger version of Theorem 3. In view of (3.9), for k / For k near ±K (that is, near K and K ′ ), and 0 < B ≪ 1, we can use a modified operator for which the equivalence in (5.1) still holds (near K and K ′ ): (5.2) In the notation of (3.25) T k (0) = T k and its spectrum is given by reciprocals of the magic α's -see Proposition 3.3.
3. Numerical evidence suggests (see Figure 7) that λ 0 (0) < 0 for the Bistritzer-MacDonald potential.If B = B 0 e 2πiθ that means the Γ point (corresponding (5.7) From (5.1) and the periodicity of the spectrum of D B (α) with respect to Λ * (see Lemma 3.1) we also note that provided that k, k + p ∈ U .This allows an extension of λ(k, B) to Ω ε in the statement of the proposition, provided that |B| < δ for some sufficiently small δ.The properties of the expansion (5.4) come from the fact that individual terms in the Taylor expansion satisfy the symmetries (5.7): ak p B q = a(−1) p k p B q = ak p B q = aω p+q k p B q =⇒ a ∈ R, p ∈ 2N, −p ≡ q mod 3. (5.9) This proves (5.3) and (5.4) except for the non-degeneracy of ∂ B ∂ 2 k λ(0, 0).To prove it, we compare the Taylor expansion of λ(k, B) with the effective Hamiltonian (4.6).Hence, with µ := 1/α, k ∈ Spec L 2 0 D B (α) then (with some modification of notation) (5.10) where F P (0, 0, λ) ̸ = 0 and k P = 0, k p > 0, p < P .That has to be so by noting that 0 / ∈ Spec L 2 (D B (α)) for 0 < |α − α| ≪ 1.Hence, (note that λ(k, B) is even in k) Combined with (5.4) this gives which implies that P = 1 and c 1 = c 0 , and concludes the proof.
We also note that comparing this conclusion with (4.6) shows that g 1 (α) ̸ = 0. □ At the bifurcation point, the Bloch eigenvalues exhibit a quadratic well, see Figure 3.
Proposition 5.2.Under the assumptions, and in the notation, of Proposition 5.1 and (5.5), let Proof.This follows from (4.9) and (5.5). □ The next proposition deals with the vertices of the boundary of the Brillouin zone.In view of (5.8) it is enough to consider one of the vertices, say, (5.11) We then have Proposition 5.3.Suppose that λ is a simple eigenvalue of T k = T k (0) and that assumptions of Theorem 2 hold for α = 1/λ.Then, for k near k 1 given in (5.11), where q ≥ 1, c 2 ∈ R \ {0} and λ 2 (z) = λ 2 (z).
2. If we know (which can be checked numerically for the Bistritzer-MacDonald potential) that g 1 (k 1 , α) ̸ = 0 and q = 1 then This follows from a comparison with the effective Hamiltonian in (4.6).
Proof of Proposition 5.3.From (5.7) and periodicity, λ(k + Λ * , B) = λ(k, B), we conclude that (note that 2k (5.14) We also note that for k / Suppose that c 2 = 0. {Proposition 5.1 already shows that λ(k, B) cannot be independent of k.This and and c 3 = 0 imply that for some q ≥ 1, p > 1, Then for a fixed 0 < |B| ≪ 1 we would have 2p solutions k near k 1 (recall that k → λ(k, B) is holomorphic near k 1 ).However, Theorem 2 shows that there are at most two solutions for a fixed λ(k, B).This gives a desired contradiction and shows that c 2 ̸ = 0. □

Proofs of Theorems 2 and 3
Combining the results of previous of sections we can now prove the main results of this paper.
Proof of Theorem 2. In the notation of (4.6) we see the effective Hamiltonian for D B (α) for B small Since θ(z(k)) ̸ = 0 for k / ∈ Λ * (see §3.6) there exists a constant a 1 such that if |α − α| < a 1 |B| then E B −+ (k) is not identically 0 (provided that B is small enough).This shows invertibility at some k and hence discreteness of the spectrum (by the analytic Fredholm theory applied to k → On the other hand, we can put k = 0 and recall from the proof of Proposition 5.1 (see (4.6)) that Again, that implies discreteness of the spectrum.We now note that there exists δ 0 > 0 such that (D(0, δ and this proves discreteness of the spectrum of D B (α) for 0 < |B| < δ 0 and |α−α| < δ 0 .We also see that (6.1) implies (2.16): for U ⋐ C, for any epsilon there exits ρ > 0 such that |θ(z(k)) 2 | > ρ for z ∈ U \ (Λ * + D(0, ε)).But then provided B is small enough (depending on ε, note that α = α in the calculation; the answer has to be an integer).
We now need to account for the possibility that D B (α) has an eigenvalue on ∂F .Periodicity of the spectrum shows that if k 1 ∈ Spec D B (α) ∩ ∂F then k 1 + γ ∈ ∂F for a finite number of γ ∈ Λ * (from the definition of a fundamental domain).Only one of these points can be in the fundamental domain F and a small deformation includes it in the interior of (the new) F , while excludes all others from ∂F .The previous argument shows that the number of eigenvalues remains 2. □ Since Λ * = Λ * this means that Spec L 2 0 D B (α) ⊂ (R + Λ * ) ∪ (iR + Λ * ) which is the same as (2.17).
To prove (2.18) we recall that C × (C \ K 0 ) ∋ (B, k) → T k (B), is a holomorphic family of compact operators with simple eigenvalue µ = 1/α ∈ Spec(T k (0)).We define K := R \ k ′ ∈K 0 B(k ′ , ε), then by periodicity of the spectrum of D B (α) it suffices to restrict us to a fundamental domain: Since K /Λ * is a compact set, the spectrum of K ∋ k → T k (B) is uniformly continuous on B in compact sets.Thus, for 0 < |B| < δ 0 small enough, the operator T k (B) has precisely one eigenvalue in a δ 1 neighbourhood of µ for every k.This implies that for every k ∈ K /Λ * there is precisely one µ k such that µ k ∈ Spec(T k (B)) and |µ k − µ| < δ 1 .From Propositions 5.1 and 5.3 we conclude that µ k ∈ R and the result follows.□ Remark.While our proof does not show that for B ∈ R \ {0} the points K, K ′ are also in the spectrum of D B (α) for some real α between successive magic angle, the bottom figure in Figure 8 shows that this is indeed the case.For general B / ∈ R this is however false, as the top figure in Figure 8 shows.Both figures exhibit an interesting universal pattern for |α| large.

Figure 2 .
Figure 2. The dynamics of Dirac points for the Bistritzer-MacDonald potential U (z) = U BM = 2 k=0 ω k e 1 2 (z ωk −zω k ).The magnetic field given by B = B 0 e 2πiθ with B 0 = 0.1 On the left different colours correspond to different values of θ shown in the colour bar and α varies between 0.1 and 0.9 (this is a colour map version of the left panel of Figure1).On the right, the colours correspond to different values of α shown in the colour bar and θ varies.The predominance of green (corresponding to the range between 0.5 and 0.6) means that most of the motion happens near the (first) magic alpha -see https://math.berkeley.edu/~zworski/ first_band.mp4 for E 1 (α, k)/ max k E 1 (α, k) for fixed B as α varies.

Figure 3 .
Figure 3.When B is real (in the convention of (2.1)), two Dirac cones approach Γ point as α → α * = α + O(B 3 ) (α a simple real magic parameter) on the line Im k = 0 (left).For α = α * , the quasi-momentum k at which the bifurcation happens are the boundary of the Brillouin zone and the Γ-point which is shown in the figure (right).The animation https://math.berkeley.edu/~zworski/Rectangle_1.mp4 shows the motion of Dirac points in this case.

Figure 5 .
Figure 5. Dirac point trajectory for B = 0.1 (left) and B = 0.1ω (right).The bifurcation happens at Γ and one additional point (modulo Λ * ) in each figure, respectively.The colors indicate the position of the Dirac cones for given values of α.The exclusion of K and K ′ points in the statement of Theorem 3 seems to be a technical issue, as shown in https://math.berkeley.edu/~zworski/Rectangle_2.mp4 (for the case of the figure on the right).
.21) which in view of Lemma 3.1 concludes the proof.□ Remark.Under the assumptions of Proposition 3.2 the corresponding −k 1 , −k 2 ∈ C/Λ * are called the K and K ′ points in the physics literature.The remaining element of K/Λ * is called the Γ point.

Figure 6 .
Figure 6.Bifurcation for B = 0.1.(Top): The color-coding indicates the position of the Dirac points for given values of α ∈ R. The right figure illustrates the bifurcation at Γ and the left figure at a non-equivalent (modulo Λ * ) bifurcation point that is a vertex of the boundary of the Brillouin zone, see Figure 1.

Figure 7 .
Figure 7. Plot of B → λ 0 (B 3 ) of Proposition 5.1 for the first (left) and second to fourth magic angle (right).

Figure 8 .
Figure 8. Top figure showing α∈ C such that 1/α ∈ Spec L 2 0 ( T K (B)) or K ∈ Spec L 2 0 (D B (α)).We see that indeed for B ∈ R\{0} the trajectory of Dirac points passes through K, K ′ .Bottom figure showing α ∈ C such that 1/α ∈ Spec L 2 0 ( T K (B)) or K ∈ Spec L 2 0 (D B (α)).For general B / ∈ R the trajectory of Dirac points for varying α ∈ R does not pass through K between successive real magic angles.