The First Magic Angle Beyond the Chiral Limit in Twisted Bilayer Graphene

Introduction

Twisted van der Waals (vdW) materials provide a versatile platform for controlling electronic structure through geometry. A small relative rotation between adjacent layers produces a moiré superlattice with a period much larger than the atomic lattice constant, strongly modifying the low-energy spectrum and enhancing the role of interactions, topology, and quantum interference effects [1, 2, 3, 4, 5, 6, 7]. Among these systems, twisted bilayer graphene (TBG) remains the paradigmatic example. Near the first magic angle, approximately 1.1°, the electronic bandwidth close to charge neutrality is strongly suppressed, producing narrow bands in which Coulomb interactions become comparable to, or larger than, the kinetic energy scale [4]. This regime has enabled correlated insulating states, superconductivity, orbital magnetism, and topological phases [2, 8, 9, 10, 11, 12].

The broader relevance of moiré materials follows from the same geometric mechanism. By enlarging the unit cell, the moiré pattern rescales the characteristic kinetic, interaction, and magnetic-field energy scales, allowing regimes that are difficult to access in ordinary two-dimensional crystals [1, 2, 4, 5]. This tunability has motivated extensive studies not only of TBG, but also of twisted transition-metal dichalcogenides, multilayer graphene systems, heterobilayers, hybrid moiré structures, and one-dimensional moiré systems such as double-walled carbon nanotubes [13, 14, 15, 16, 17, 18, 19, 20]. These platforms exhibit a wide range of correlated and topological phenomena, including generalized Wigner crystals, Chern insulators, fractional Chern insulators, heavy-fermion-like behavior, magnetic textures, and superconducting phases [7, 10, 11, 14, 15, 21, 22]. The same control over twist angle, pressure, electrostatic gating, strain, and external fields has also motivated proposals for quantum simulation and quantum-information applications based on moiré flat bands and interaction-driven topological states [2, 7, 23, 24].

A central theoretical description of TBG is the continuum model of Bistritzer and MacDonald, in which two rotated Dirac Hamiltonians are coupled by spatially modulated interlayer tunneling terms [4]. This model captures the emergence of narrow bands at a sequence of magic angles. A particularly important idealization is the chiral limit, obtained by suppressing the same-sublattice tunneling amplitude while keeping the opposite-sublattice tunneling finite. In this limit, the Hamiltonian has an additional chiral symmetry and supports exactly flat bands at a discrete set of magic angles [25]. The chiral model has revealed a highly structured mathematical problem, including zero-mode wave functions, a hierarchy of magic angles, a characteristic scaling rule, and close analogies with Landau-level physics [25, 26]. These developments connect TBG to broader ideas in quantum geometry, generalized Landau-level descriptions of Chern bands, and effective gauge-field formulations of moiré flat bands [27, 28, 29, 30].

Realistic TBG, however, is not exactly chiral. Atomic relaxation changes the relative area of local stacking regions, reducing the effective tunneling through AA regions while enlarging the energetically favorable AB and BA domains [31, 32, 33]. This effect is commonly incorporated in the continuum model through a finite ratio between same-sublattice and opposite-sublattice interlayer tunneling amplitudes. In addition, realistic systems may contain strain, particle-hole asymmetry, interlayer-coupling anisotropy, and other symmetry-breaking perturbations [31, 32, 33]. These corrections modify the ideal chiral spectrum. In particular, numerical and analytical studies indicate that the higher-order chiral magic angles are strongly destabilized, whereas the first magic angle remains comparatively robust near experimentally relevant parameters [25, 33, 34, 35]. Understanding why this first angle survives more effectively than the higher-order members of the chiral hierarchy remains an important theoretical question.

One useful route to this problem is provided by the squared-Hamiltonian formulation of the chiral model. In this representation, the original Dirac-like continuum Hamiltonian is mapped to an effective Schrödinger-like operator containing confinement terms, pseudo-magnetic gauge fields, and interlayer matrix elements. Previous work has related this structure to effective magnetic flux scaling, current-like interlayer contributions, and quantum-Hall-type descriptions in momentum or moiré space [36, 37, 38, 39, 40, 41, 42]. These ideas suggest that the loss of exact flatness in the non-chiral model should not be viewed only as a numerical shift of the band energies, but rather as a modification of the effective confinement and interlayer coupling structure inherited from the chiral problem.

In this work, we study this evolution from the chiral to the non-chiral Bistritzer–MacDonald model within an effective single-particle framework. We focus on how a finite same-sublattice tunneling ratio, motivated by lattice relaxation, modifies the squared Hamiltonian and the associated flat-band hierarchy. We show that finite same-sublattice tunneling renormalizes the diagonal confinement potential, introduces additional off-diagonal interlayer matrix elements, enhances real-space localization around AA regions, and increases hybridization between the central bands and remote bands. A key outcome is that the commonly used relaxation-renormalized value $\kappa\simeq0.7$ lies close to $\kappa= 1/\sqrt{2}$, where the oscillatory part of the symmetric confinement in the squared Hamiltonian cancels. These effects suppress the higher-order chiral magic-angle structure while leaving the first magic angle comparatively robust. Rather than invoking many-body physics, the analysis identifies the single-particle mechanisms by which non-chiral tunneling modifies spectral squeezing, wave-function localization, and current-like interlayer channels in magic-angle TBG.

Chiral and non-chiral TBG Hamiltonians

The Bistritzer–MacDonald (BM) Hamiltonian provides the standard continuum description of the low-energy electronic structure of twisted bilayer graphene (TBG) [43]. It is constructed from the Dirac Hamiltonians of two rotated graphene layers coupled by spatially modulated interlayer tunneling terms. The dominant tunneling processes are parametrized by the amplitudes $w_{AA}$ and $w_{AB}$, which describe same-sublattice tunneling through local AA regions and opposite-sublattice tunneling through AB/BA regions, respectively.

Lattice relaxation plays an essential role in realistic TBG. Atomic reconstruction reduces the spatial extent of AA regions and enlarges the energetically favorable AB/BA domains, thereby suppressing the effective AA tunneling amplitude [31, 32, 33]. The ideal chiral limit is obtained by setting $w_{AA}=0$ while keeping $w_{AB}$ finite. This limit, introduced by Tarnopolsky, Kruchkov, and Vishwanath, gives the Hamiltonian an exact chiral symmetry and supports exactly flat bands at a discrete set of magic angles [44]. It therefore provides the natural reference point for analyzing how finite same-sublattice tunneling modifies the magic-angle hierarchy.

We work in a single-valley continuum theory. In the sublattice basis, the four-component bi-spinor is written as

where the indices 1, 2 label the graphene layers. The amplitudes $\psi_j(\mathbf{r})$ and $\chi_j(\mathbf{r})$ correspond to the continuum envelope functions associated with the two graphene sublattices. It is useful to introduce the layer spinors

In this representation, $\Psi$ and $\chi$ define the two sublattice sectors, while the entries inside each spinor refer to the layer degree of freedom. The numerical calculations are performed in a Bloch plane-wave basis. For the $\Psi$ sector, we use

where $a_{m,n}$ and $b_{m,n}$ are Fourier coefficients, $\mathbf{K}_{m,n}=m\mathbf{b}_1+n\mathbf{b}_2$ is a moiré reciprocal vector, and $\mathbf{k}$ is the crystal momentum in the moiré Brillouin zone. The same type of Bloch expansion is used for the $\chi$ sector. The additional phase in the second layer fixes the gauge convention used below for the interlayer coupling functions. In the sublattice-block basis $(\Psi,\chi)$, the non-chiral BM Hamiltonian can be written as [44, 45, 46, 47]

The operators $D$ and $D^\dagger$ contain the Dirac kinetic energy and the opposite-sublattice interlayer tunneling,

and

Here $\bar{\partial}=\partial_x+i\partial_y$ and $\partial=\partial_x-i\partial_y$. The functions $U_m(\mathbf{r})$ describe the moiré modulation of the interlayer coupling,

where $\phi=2\pi/3$ and $m=0,\pm1$. With the conventions used here, these functions satisfy

which relates the tunneling functions to the derivative structure that will appear in the squared Hamiltonian. The diagonal block $M$ contains the same-sublattice interlayer tunneling,

where $\kappa=w_{AA}/w_{AB}$. The parameter $\kappa$ measures the deviation from the chiral limit. When $\kappa=0$, the same-sublattice tunneling vanishes and the Hamiltonian anticommutes with the chiral operator. For finite $\kappa$, the operator $M$ is diagonal in the sublattice-block structure of (4) but off diagonal in layer space. It therefore represents tunneling between identical sublattices in opposite layers, namely $A_1\leftrightarrow A_2$ and $B_1\leftrightarrow B_2$. In realistic relaxed structures, AA tunneling is reduced but not eliminated, and values close to $\kappa\approx0.7$ are commonly used to describe relaxation-renormalized TBG [48].

The moiré coupling vectors entering the interlayer potentials are

where $k_\theta=2k_D\sin(\theta/2)$ is the characteristic moiré momentum scale and $k_D=4\pi/(3a_0)$ is the magnitude of the graphene Dirac momentum, with $a_0$ the monolayer lattice constant. Since $k_\theta$ decreases with decreasing twist angle, the moiré length scale grows rapidly at small angles.

The continuum model is controlled by the dimensionless coupling parameter

where $w_1\equiv w_{AB}$ and $v_0$ is the Fermi velocity of monolayer graphene. Throughout this work we use units such that $v_0=1$ and $k_\theta=1$, so that lengths are measured in units of $k_\theta^{-1}$ and energies in units of $v_0k_\theta$. In these units, the twist-angle dependence enters through $\alpha$: decreasing the twist angle increases $\alpha$ and strengthens the effective interlayer hybridization relative to the intralayer kinetic energy. In reciprocal space, the moiré Brillouin zone is generated by the reciprocal lattice vectors

The high-symmetry points used below are

following the conventions of Ref. [37]. The corresponding moiré real-space lattice vectors can be written compactly as

with the sign chosen consistently with $\mathbf{b}_{1,2}$. In this gauge, the tunneling functions acquire phase factors under primitive moiré translations, while translations by three primitive periods leave them invariant:

This enlarged periodicity is closely related to the magnetic-translation structure encountered in the quantum-Hall interpretation of the chiral model [38]. It also provides a useful way to organize the Fourier harmonics entering the numerical diagonalization.

The BM Hamiltonian possesses crystalline and antiunitary symmetries that strongly constrain its spectrum and eigenstates. As discussed in Refs. [44, 49, 50, 51, 52], the single-valley continuum model is constrained by $C_{3z}$, $C_{2x}$, $C_{2y}T$, and $C_{2z}T$, where $T$ denotes time reversal. These symmetries restrict the band degeneracies, the valley structure, and the topology of the low-energy bands. The interlayer potentials satisfy

while

As a consequence,

Within the symmetry convention adopted here, the eigenstates can therefore be chosen so that the two sublattice sectors are related by inversion and complex conjugation,

This relation will be used later when comparing expectation values computed in the $\Psi$ and $\chi$ sectors.

In the exact chiral limit, the energy spectrum is symmetric with respect to $E = 0$. If $\lvert\Phi_{+}\rangle= (\psi_{1},\psi_{2},\chi_{1},\chi_{2})^{T}$ is an eigenstate with energy $E$, then $\lvert\Phi_{-}\rangle= (\psi_{1},\psi_{2},-\chi_{1},-\chi_{2})^{T}$ is an eigenstate with energy $-E$. The zero-energy states are therefore special. They satisfy

and, because they are degenerate, they can be chosen to be polarized in one sublattice sector:

or alternatively,

At the chiral magic angles, these zero modes form the completely flat bands of the TKV model. Their existence is protected by the chiral structure of the Hamiltonian and by the topology of the zero-mode problem. They have also been connected to topological solitons protected by crystalline symmetries and to compact localized states in flat-band lattice systems [53, 54, 55]. The finite-$\kappa$ problem studied in the following sections can therefore be viewed as a deformation of this chiral zero-mode structure. Same-sublattice tunneling breaks exact flatness, but part of the symmetry and topology inherited from the chiral limit remains visible in the low-energy bands.

Spectral squeezing and wave-function localization

We now analyze the spectrum and wave functions obtained from the continuum Hamiltonian in (4). Figure 1 compares the chiral limit, $\kappa= 0$, with the relaxation-renormalized non-chiral model, $\kappa= 0.7$, for representative values of the coupling parameter $\alpha$ corresponding to the first and higher-order magic-angle regimes. This comparison isolates the effect of finite same-sublattice tunneling on the low-energy bands. In the perturbative regime $\kappa\alpha\ll1$, the numerical results can also be compared with the analytical expansion discussed in the Appendix.

A first diagnostic is provided by the behavior of the band-edge energy at the $\Gamma$ point. In Figure 2 we plot $E^{2}$ as a function of $\alpha$ for the low-energy state that controls the band edge in the numerical spectrum. Together with the band structures in Figure 1, this quantity shows the qualitative difference between the chiral and non-chiral regimes. In the exact chiral limit, the bandwidths $W_n$ associated with the low-energy bands follow the characteristic squeezing behavior

This progressive narrowing of the spectrum is one of the central signatures of the chiral model. The origin of this squeezing mechanism can be understood from the real-space rescaling introduced in the chiral theory [36],

This transformation exposes the quantum-Hall-like structure of the squared chiral Hamiltonian in the large-$\alpha$ regime [36, 41]. In this limit, corrections associated with the effective Zeeman contribution become asymptotically small [36, 37, 41]. Under the same scaling, the dominant operators in the non-chiral Hamiltonian behave as

Therefore, finite $\kappa$ introduces terms that scale with the same power of $\alpha$ as the chiral interlayer coupling. These terms compete directly with the chiral scaling structure and obstruct the exact renormalization mechanism responsible for ideal spectral squeezing. This obstruction is visible already in the low-energy bands. In the chiral model, exact flatness is recovered at the magic angles of the TKV theory, while away from those values the bands display the expected squeezing with increasing $\alpha$. Once finite $\kappa$ is included, the central bands acquire residual dispersion and the higher-order magic-angle structure is strongly modified. The magic-angle values are also shifted relative to the chiral prediction. Thus, the realistic model does not preserve the exact flatness of the ideal chiral continuum theory, in agreement with symmetry-based analyses showing that atomic-scale models are not expected to retain the full idealized flat-band structure [56]. The effect is especially pronounced near the $\Gamma$ point, where the band edge is located in the spectra considered here. For the first magic angle, the non-chiral bands remain narrow and close to the chiral result. At higher-order magic angles, however, finite $\kappa$ produces a much stronger deviation: the squeezing of the low-energy spectrum is partially blocked, and the separation between the central bands and the first remote bands is strongly reduced. This indicates that remote-band effects become increasingly important at small twist angles and that isolated-flat-band descriptions become less reliable in this regime. The low-energy states near the high-symmetry points $K$ and $K'$ remain close to charge neutrality even when $\kappa$ is finite. This behavior reflects the fact that the same-sublattice tunneling term preserves the crystalline symmetries of the continuum model while breaking the exact chiral anticommutation symmetry. Hence, part of the symmetry and topology inherited from the chiral limit remains visible in the non-chiral spectrum, although the exact spectral symmetry around zero energy is lost.

A second consequence of finite $\kappa$ is the enhanced hybridization between the central bands and the first remote bands. In the chiral model, the central bands are protected by the zero-mode structure and remain well separated from the remote bands at the magic angles. In the non-chiral model, AA tunneling opens additional same-sublattice tunneling channels and mixes states that are effectively decoupled in the chiral limit. As shown in Figure 1, this mixing becomes particularly strong at larger $\alpha$, where the remote-band gap becomes very small and the central bands are no longer cleanly isolated from the first excited bands. The corresponding real-space effect is shown in Figure 3. At the first magic angle, the chiral model displays a broad and symmetric density distribution around the AA regions. When finite $\kappa$ is introduced, the density profile becomes narrower, indicating stronger localization. The underlying $C_3$ symmetry is preserved, but the electronic weight contracts toward the center of the AA stacking region. This localization becomes much more pronounced at higher-order magic angles. For $\alpha_9 = 12.860$, the non-chiral density is concentrated in a substantially smaller spatial region than in the chiral case. Finite same-sublattice tunneling therefore acts as an effective localization mechanism for the low-energy wave functions. This enhanced localization is expected to increase local interaction matrix elements and may modify the balance between kinetic-energy suppression and Coulomb interactions in realistic moiré systems.

The real-space localization is consistent with the redistribution of Fourier weight shown in Figure 4. We plot the coefficients $\lvert a_{mn}\rvert^2$ for the chiral and non-chiral models at $\alpha_1=0.586$ and $\alpha_9=12.860$. At the first magic angle, the two distributions remain qualitatively similar: most

of the weight is concentrated around a small number of symmetry-related reciprocal-space points. Finite $\kappa$ nevertheless changes the relative intensity and angular distribution of the dominant components. At the higher-order magic angle, the difference becomes much stronger. The dominant Fourier components move to larger $|\mathbf{K}_{mn}|$ in the non-chiral model, indicating that the real-space wave function is more tightly localized. The approximately sixfold pattern observed in the chiral limit evolves into a distribution with stronger trigonal anisotropy, while the underlying $C_3$ symmetry of the continuum Hamiltonian remains preserved. Thus, finite $\kappa$ does not simply broaden the bands; it also reorganizes the wave functions in reciprocal space and concentrates the electronic density in real space.

Overall, finite same-sublattice tunneling qualitatively modifies the chiral flat-band hierarchy. It suppresses the ideal spectral squeezing, enhances hybridization between the central and remote bands, and localizes the low-energy wave functions more strongly around AA regions. These effects show that the realistic BM model remains continuously connected to the chiral theory, but with important finite-$\kappa$ corrections that become increasingly relevant at higher-order magic angles. In the next section, we analyze these corrections directly from the squared Hamiltonian.

Squared TBG Hamiltonian with finite-$\kappa$ corrections

The previous section showed that finite same-sublattice tunneling suppresses the chiral spectral squeezing, enhances hybridization with remote bands, and localizes the low-energy wave functions more strongly around the AA regions. We now analyze the origin of these effects from the squared Hamiltonian. In the chiral limit, the squared formulation reduces the problem to an effective $2 \times2$ operator, where the relation between flat bands, confinement, and emergent pseudo-magnetic fields becomes transparent [36, 37, 38, 39, 41]. Starting from (4), the square of the non-chiral Hamiltonian is

where $\{M,D\}$ denotes the anticommutator between $M$ and $D$. In the exact chiral limit, $M = 0$, the anticommutator terms vanish and the two sublattice sectors decouple. The squared Hamiltonian then consists of the partner operators $D^\dagger D$ and $DD^\dagger$, whose zero modes generate the flat bands of the chiral model.

For finite $\kappa$, this decoupling is no longer exact. Equation (26) contains two distinct non-chiral corrections. The term $M^2$ modifies the diagonal confinement potential, while the anticommutators $\{M,D\}$ and $\{M,D^\dagger\}$ couple the two sectors of the squared problem. These two contributions provide the microscopic mechanism behind the numerical trends observed in Figure 1: the loss of ideal squeezing, the enhanced central–remote-band hybridization, and the stronger real-space localization.

The Schrödinger equations associated with Eq. (26) are

and

Equation (28) can be formally solved for $\chi$ and substituted into Eq. (27). This gives an effective equation for the $\Psi$ sector,

This expression shows that the non-chiral model cannot be reduced to the chiral squared Hamiltonian by a simple parameter renormalization. The diagonal term $M^2$ changes the effective confinement, while the second term in Eq. (29) produces an energy-dependent coupling to the opposite sublattice sector. A perturbative expansion is possible in the regime $\kappa\alpha\ll1$, as discussed in the Appendix, and provides a useful check of the numerical results at small $\alpha$. For the realistic value $\kappa\simeq0.7$, however, the full squared-Hamiltonian structure must be retained.

It is useful to separate the diagonal and off-diagonal parts of Eq. (26) as

where

The two diagonal blocks are

with

Thus $H^2$ and $\widetilde{H}^2$ are the two squared chiral partner Hamiltonians, while $M^2$ gives the diagonal finite-$\kappa$ correction. The squaring procedure is closely related to supersymmetric mappings of bipartite Hamiltonians [58, 59, 60, 61]; in the present context, its main advantage is that it exposes the confinement and pseudo-magnetic structure of the chiral problem [36, 37].

To keep the matrix expressions compact, we define the diagonal operators

where $X(r)$ denotes the symmetric part of the diagonal confinement. We also define

The notation $\mathcal{J}_{-}$ denotes the lower-left operator appearing in the squared Hamiltonian; it should not be interpreted as the complex conjugate of $\mathcal{J}_{+}$, since the factor $-2i$ is kept fixed by the operator structure. Hermiticity is understood with the derivative operators acting on periodic functions in the moiré unit cell; $\mathcal{J}_{-}$ is the lower-left differential operator in this convention, not the naive complex conjugate of $\mathcal{J}_{+}$.

With this notation, the squared chiral Hamiltonian takes the form [62]

The complex vector quantity entering (35) is

It acts as an emergent pseudo-magnetic vector potential. It satisfies the Coulomb gauge condition $\nabla\cdot\mathbf{A}=0$ and generates an effective field perpendicular to the graphene plane,

We denote the corresponding out-of-plane scalar component by $B$, which in the present convention is

The complex structure of this field follows from the valley structure and the phase winding generated by the moiré tunneling functions [38]. The diagonal potentials of the chiral squared Hamiltonian can be decomposed into a symmetric confinement term and an antisymmetric pseudo-magnetic contribution,

Thus $A^2(r)$ defines the symmetric confinement landscape, while $\Delta(r)$ encodes the antisymmetric part of the effective pseudo-magnetic structure. We now include the diagonal non-chiral correction $M^2$. Since

finite $\kappa$ renormalizes the symmetric confinement term according to

with

This expression explains why finite $\kappa$ changes the localization properties found in the previous section: the same-sublattice tunneling does not only break chiral symmetry, but also reshapes the effective confinement potential experienced by the low-energy states. Using (43), the finite-$\kappa$ block becomes

An especially important reference point occurs for

which is strikingly close to the commonly used relaxation-renormalized value $\kappa\simeq0.7$. At this value, the oscillatory cosine contribution in (43) vanishes exactly. Therefore, finite same-sublattice tunneling does not merely perturb the chiral confinement landscape; near the realistic value of $\kappa$, it cancels the spatial modulation of the symmetric confinement term and leaves a uniform contribution,

where

Here $\sigma$ is the localization scale obtained from the asymptotic scaling analysis [36]. This cancellation provides a direct analytical interpretation of the relaxation-renormalized BM parameter: the realistic non-chiral model lies close to a special confinement point of the squared Hamiltonian. The fact that realistic TBG samples appear naturally close to this special value suggests the intriguing possibility that lattice relaxation favors a nearly optimized confinement regime, in which the oscillatory part of the symmetric potential is strongly suppressed. The remaining spatial structure of the diagonal block is then controlled by the antisymmetric pseudo-magnetic term $\Delta(\mathbf{r})$, while the off-diagonal anticommutator terms discussed below provide the additional coupling between the two squared-Hamiltonian sectors. Under the approximation $\kappa\approx1/\sqrt{2}$, the diagonal block simplifies to

We finally analyze the off-diagonal blocks in Eq. (30). These terms are absent in the chiral limit and are responsible for coupling the two squared-Hamiltonian sectors. Direct evaluation gives

where $p=p_x-ip_y$, with $p_j=-i\partial_j$ in the dimensionless units used here. Similarly,

The scalar potential

is a complex sublattice-sector coupling generated by finite same-sublattice tunneling. Explicitly,

The first term in Eqs. (49) and (50) scales as $\kappa\alpha^2$, while the momentum-dependent term scales as $\kappa\alpha$. These off-diagonal contributions therefore become increasingly relevant at larger $\alpha$, precisely where the numerical spectra show stronger remote-band hybridization and a breakdown of the ideal chiral squeezing.

Figure 5 summarizes the real-space structure of the effective fields entering the squared Hamiltonian. The symmetric potential controls the confinement landscape, the antisymmetric term carries the pseudo-magnetic structure, and the complex potential $V_S(\mathbf{r})$ represents the additional sublattice-sector coupling generated by finite $\kappa$. Thus the two main effects found numerically have a direct origin in the squared Hamiltonian: $M^2$ renormalizes the confinement landscape and changes the localization scale, while the anticommutator terms couple the chiral sectors and enhance hybridization with remote bands.

Energetic contributions in the non-chiral model

In this section, we analyze the energetic contributions associated with the different terms entering the non-chiral squared Hamiltonian. The expectation values provide direct information about the physical mechanisms governing the localization of the electronic states, the suppression of the bandwidth, and the redistribution of interlayer processes generated by finite same-sublattice tunneling. Starting from Eq. (26), the expectation value of the squared Hamiltonian is

The two diagonal sectors are related by the symmetry relation in Eq. (19) for the states considered here. Therefore, the expectation values can be evaluated using one sector and then compared with the corresponding symmetry-related contribution from the other sector. In the numerical results below, all plotted quantities are reported using the same normalization convention, so that the relative weights of the kinetic, confinement, pseudo-magnetic, and interlayer terms can be directly compared.

Using the explicit operators introduced in the previous section, the energetic decomposition can be written as

Here $\Delta(\mathbf{r})=i[A_x,A_y]$, and the signs are inherited from the squared-Hamiltonian convention used in Eq. (44) and in the anticommutator terms of Eqs. (49) and (50). The first line of Eq. (54) contains the kinetic and diagonal confinement energies, while the second line contains the chiral orbital interlayer contribution, the pseudo-Zeeman contribution, and the non-chiral inter-sublattice terms activated by finite $\kappa$.

The kinetic contribution is

with the symmetry-related contribution from the $\chi$ sector included through the normalization convention described above. The diagonal symmetric confinement term is

Using Eq. (43), this becomes

where

is the norm of the sector used in the expectation value. If the sector density is normalized to the moiré unit-cell area, then $\mathcal{N}_{\Psi}=A_M$, with $A_M=8\pi^2/(3\sqrt{3})$. Equation (57) is the energetic counterpart of the confinement result derived in the previous section. The uniform part grows as $3(1+\kappa^2)$, while the oscillatory part is controlled by the coefficient $-1+2\kappa^2$. Therefore, at

the spatially oscillating contribution vanishes exactly. Since the relaxation-renormalized value $\kappa\simeq0.7$ lies extremely close to this point, realistic TBG is naturally placed near a regime in which the symmetric confinement is dominated by a uniform contribution and the oscillatory modulation is strongly suppressed. In this regime,

or equivalently $\langle A_f^2\rangle\simeq(9/2)A_M$ when the sector normalization is chosen as $\mathcal{N}_{\Psi}=A_M$. This is not a minor correction to the chiral model; it identifies a special confinement point selected by the realistic value of the relaxation-renormalized BM parameter. The remaining diagonal term is the antisymmetric pseudo-magnetic contribution,

This contribution remains bounded because $\Delta(\mathbf{r})$ is a bounded periodic function. Thus, in (54), the diagonal confinement terms enter with an overall prefactor $\alpha^2$. The uniformization of $A_f^2$ near $\kappa\simeq1/\sqrt{2}$ therefore has a direct energetic consequence: it changes the leading $\alpha^2$ confinement channel while leaving the remaining spatial dependence of the diagonal block mainly controlled by $\Delta(\mathbf{r})$.

The chiral off-diagonal contribution inside the squared Hamiltonian is

where

This term contains one spatial derivative and therefore becomes increasingly important as the wave functions localize in real space. In the large-$\alpha$ regime, the characteristic gradients increase with the inverse localization length, so the full contribution $-\alpha\langle A \cdot\nabla\rangle$ competes directly with the diagonal $\alpha^2$ confinement energy.

The pseudo-Zeeman contribution is

Unlike Eq. (62), this term contains no derivative. It is therefore most visible at the first magic angle and becomes relatively less important when the derivative-controlled interlayer contributions grow at larger $\alpha$. This makes the first magic angle qualitatively distinct from the higher-order members of the chiral hierarchy [37].

We now consider the non-chiral off-diagonal contribution generated by the anticommutators in Eqs. (49) and (50). The derivative part is

where $p=p_x-ip_y$. This term is activated only when $\kappa\ne0$ and enters Eq. (54) with the prefactor $\kappa\alpha$. Because it also contains a momentum operator, it is strongly enhanced when the wave functions become more localized. This contribution is therefore one of the main energetic channels through which finite same-sublattice tunneling reorganizes the low-energy spectrum.

The anticommutator also contains a local derivative piece, since

The commutator term is proportional to a derivative of the pseudo-magnetic field and can be written as a local sublattice-mixing contribution. We denote its expectation value by

where $D_B$ is the local operator obtained from the derivative of $B$. In the numerical results for the lowest-energy band considered here, this local sublattice-mixing contribution vanishes within numerical precision,

Thus the dominant non-chiral inter-sublattice contribution comes from the derivative anticommutator term $\langle A_{g,s}\rangle$.

The scalar off-diagonal potential generated by the anticommutator is

and its expectation value is

For the $C_3$-symmetric states analyzed here, this expectation value vanishes by the symmetry of the moiré unit cell,

Equivalently, the three first-star harmonic contributions transform into one another under $C_3$, and their phases cancel after integration over the unit cell. Therefore, although $V_S(r)$ is an important operator in the squared Hamiltonian, its direct contribution to the total energy of these states vanishes. The non-chiral effects seen numerically are instead dominated by the confinement renormalization and by the derivative anticommutator channel $\langle A_{g,s}\rangle$.

Equation (54) therefore separates the physical mechanisms entering the non-chiral model. The term $\langle A^2\rangle$ is the symmetric confinement energy. Near $\kappa\simeq0.7$, this term is close to the special uniform-confinement point $\kappa=1/\sqrt{2}$. The term $\langle\Delta\rangle$ gives the remaining antisymmetric pseudo-magnetic structure of the diagonal block. The term $\langle A\cdot\nabla\rangle$ represents the chiral orbital interlayer contribution, while $\langle B\cdot\tau\rangle$ is the pseudo-Zeeman contribution. Finally, $\langle A_{g,s}\rangle$ is the leading non-chiral inter-sublattice contribution generated by finite same-sublattice tunneling.

In the exact chiral limit, $\kappa=0$, Eq. (54) reduces to the expectation-value equation of the TKV model,

For the symmetry-adapted states considered here, the unit-cell average of the antisymmetric contribution does not contribute to the total energy,

so that the chiral energy balance is controlled by kinetic energy, symmetric confinement, and the orbital and pseudo-Zeeman interlayer terms.

Figure 2 shows the behavior of $E^2(\alpha)$ at the $\Gamma$ point for both the chiral and relaxation-renormalized non-chiral models. In the exact chiral limit, the numerical data display the expected squeezing of the spectrum, with the band edge decreasing rapidly near the chiral magic angles. By contrast, in the non-chiral model with $\kappa\simeq0.7$, the squeezing mechanism is partially blocked. The expectation-value decomposition shows why this happens. Finite $\kappa$ places the diagonal confinement close to the uniform-confinement point, but it also activates the derivative anticommutator channel $\langle A_{g,s}\rangle$, which introduces an additional inter-sublattice energy scale. The realistic model therefore does not simply inherit the chiral cancellation mechanism; it replaces it by a new balance between confinement and competing interlayer channels.

More information is obtained from Figure 6, where the low-energy scale is shown as a function of $\kappa$ for the first three chiral magic angles. The first magic angle behaves differently from the higher-order ones. For $\alpha_1 = 0.586$, the energy varies monotonically with $\kappa$, whereas higher-order magic angles display oscillatory behavior. This shows that the first magic angle is not merely the first member of the chiral sequence; it is the member that remains most compatible with the relaxation-renormalized energy balance of the non-chiral model.

To further understand the role of lattice relaxation, Figure 7 shows the energetic contributions as a function of $\kappa$ for the first magic angle. The calculations are performed for the low-energy state associated with the flat-band sector at the $\Gamma$ point. The physically relevant interval $0.6 \lesssim\kappa\lesssim0.8$ is precisely where the confinement channel approaches the special value $\kappa= 1/\sqrt{2}$, and where the chiral and non-chiral interlayer contributions become comparable in magnitude. This is a central result of the expectation-value analysis: realistic relaxation places TBG close to a balanced regime in which the symmetric confinement is nearly uniform and the two interlayer channels compete on comparable energy scales.

The dependence of the energetic contributions on the twist-angle parameter $\alpha$ is shown in Figure 8 for the realistic value $\kappa= 0.7$. Panels (a) and (b) correspond to the high-symmetry points $\mathbf{K}$ and $\Gamma$, respectively. The vertical black lines indicate the chiral-model magic angles.

Several important trends emerge from these plots. First, both the confinement energy and the kinetic energy increase systematically with $\alpha$, reflecting the stronger real-space localization of the electronic states at smaller twist angles. Second, the pseudo-Zeeman contribution is most relevant for $\alpha< 1$, and therefore plays its most visible role near the first magic angle. Third, and most importantly, the chiral and non-chiral interlayer contributions display a pronounced complementary behavior. Whenever the chiral contribution $\langle A \cdot\nabla\rangle$ develops a local minimum, the non-chiral contribution $\langle A_{g,s}\rangle$ tends to develop a corresponding maximum. This anticorrelation is the energetic signature of an effective compensation between the two interlayer channels. It does not imply an exact conservation law at the operator level, but it shows that the realistic non-chiral spectrum is governed by a strong redistribution of weight between the chiral orbital channel and the relaxation-induced inter-sublattice channel.

The expectation-value analysis therefore identifies the energetic origin of the robustness of the first magic angle beyond the chiral limit. Finite $\kappa$ brings the symmetric confinement close to a special nearly uniform regime, while simultaneously activating a non-chiral derivative channel that competes with the chiral orbital contribution. At the first magic angle these effects remain balanced, whereas at higher-order magic angles the balance becomes oscillatory and more strongly hybridized with remote bands. This motivates the current interpretation developed in the next section, where the two interlayer channels are analyzed directly.

Currents between layers and sublattices

The expectation-value analysis of the previous section shows that finite $\kappa$ activates an additional interlayer channel through the derivative anticommutator term $\langle A_{g,s}\rangle$. We now interpret this contribution in terms of current-like Fourier components of the squared Hamiltonian. As detailed in [37], the off-diagonal terms of the squared chiral Hamiltonian can be identified with the sum of the Fourier components of an interlayer current at the three first-star momenta $\boldsymbol{q}_1$, $\boldsymbol{q}_2$, and $\boldsymbol{q}_3$. In the present notation, the chiral orbital contribution is

This quantity measures the current-like Fourier weight associated with the chiral interlayer channel. It is already present at $\kappa= 0$, and it is controlled by the emergent gauge field $A(r)$ that appears in the squared chiral Hamiltonian.

Finite same-sublattice tunneling introduces a second current-like contribution through the non-chiral anticommutator term. We define

This term is absent in the chiral limit and is activated only when $\kappa\ne0$. It couples the two sublattice sectors $\Psi$ and $\chi$ through the derivative anticommutator $\{B,p\}$. For this reason, $J_{AB}$ represents a relaxation-induced inter-sublattice channel of the squared Hamiltonian.

In what follows, we refer to $J_{AA}$ and $J_{AB}$ as currents for brevity. More precisely, they are not macroscopic transport currents, but sums of current-like Fourier components at $\boldsymbol{q}_1$, $\boldsymbol{q}_2$, and $\boldsymbol{q}_3$. This distinction is important because the quantities are extracted from the squared Hamiltonian and measure how the off-diagonal operators redistribute spectral weight between layer and sublattice sectors.

The physical interpretation is direct. The chiral channel $J_{AA}$ is the orbital interlayer contribution inherited from the TKV limit. It preserves the chiral sector structure and participates in the cancellation mechanism responsible for spectral squeezing. The non-chiral channel $J_{AB}$, by contrast, is generated by finite same-sublattice tunneling in the original BM Hamiltonian and appears in the squared Hamiltonian as an inter-sublattice derivative channel. It therefore provides a leakage path between the two sublattice sectors. This leakage breaks the strict chiral cancellation pattern and is one of the mechanisms by which higher-order chiral magic angles are destabilized in the relaxation-renormalized model.

To make this behavior explicit, Figure 9 shows the two current-like contributions for the realistic non-chiral model at the Dirac point $K$. Each contribution oscillates strongly as a function of $\alpha$, especially near the chiral magic angles. The important result is that these oscillations are not independent. The two channels display a pronounced complementary behavior: when the chiral contribution is enhanced, the non-chiral contribution is reduced, and vice versa. As a consequence, their sum remains much smoother than either component separately.

The total current-like contribution,

is found numerically to follow the approximate scaling

This empirical relation is one of the clearest signatures of the energy redistribution induced by finite $\kappa$. Although $J_{AA}$ and $J_{AB}$ separately show strong oscillatory structure, their sum tracks a smooth negative contribution with approximately quadratic scaling in $\alpha$.

This behavior can be understood from the energy balance in (53). Near the magic-angle regime, the total squared energy $E^2$ remains small compared with the individual kinetic, confinement, and interlayer contributions. The diagonal terms of the squared Hamiltonian provide the positive kinetic and confinement channels, including the dominant $\alpha^2$ confinement scale. The off-diagonal current-like terms must therefore supply a compensating negative contribution. In the chiral model this compensation is carried only by the orbital channel $J_{AA}$. In the relaxation-renormalized model, the compensation is redistributed between $J_{AA}$ and the non-chiral inter-sublattice channel $J_{AB}$. This redistribution prevents the chiral cancellation pattern from remaining intact at higher-order magic angles and explains why the realistic spectrum deviates from the ideal TKV hierarchy.

The current-balance relation is therefore not an exact conservation law imposed at the operator level. Rather, it is an emergent energetic compensation visible in the expectation values of the low-energy states. Its importance is that it identifies how relaxation modifies the flat-band mechanism without completely destroying the chiral structure. At the first magic angle, the two current-like channels remain sufficiently balanced to preserve narrow central bands. At higher-order magic angles, the oscillatory redistribution becomes stronger, remote-band hybridization increases, and the exact chiral squeezing mechanism is lost.

Breakdown of the 3/2 flux quantization rule

In the chiral squared-Hamiltonian formulation, the higher-order magic angles can be organized in terms of an effective flux-quantization condition associated with the pseudo-magnetic confinement and current balance of the TKV limit. In this description, the chiral hierarchy follows the characteristic 3/2 flux-scaling rule derived for the ideal model [36, 37, 38]. Finite same-sublattice tunneling modifies both the confinement landscape and the off-diagonal current-like channels, so the breakdown of the higher-order magic-angle sequence can be viewed as the breakdown of this chiral 3/2 quantization structure.

The previous sections showed that the non-chiral model remains continuously connected to the chiral TKV limit, but with important corrections generated by finite same-sublattice tunneling. An important consequence is that the first magic angle retains many of the characteristic features of the chiral model, whereas higher-order magic angles are much more strongly modified. This robustness is particularly significant because the strongest correlated and superconducting phenomena in twisted bilayer graphene are observed precisely near the first magic angle. A central question, therefore, is why the first magic angle is so special.

Part of the answer follows directly from the perturbative analysis. As discussed in the Appendix, the leading non-chiral corrections at small $\alpha$ scale as $(\kappa\alpha)^2$. These corrections are therefore least important in the vicinity of the first magic angle, $\alpha_1 \approx0.586$, and become progressively more relevant for higher-order magic angles. This already explains why the first member of the chiral hierarchy is the most robust against non-chiral effects. However, the current-channel analysis reveals a more specific mechanism.

Figure 10 provides additional information. Panels (a) and (b) show the two current-like contributions defined in the previous section, $J_{AA}$ and $J_{AB}$, as functions of momentum for the first and second magic angles, respectively. At the first magic angle, $\alpha_1 \approx0.586$, the two channels become comparable near the $\Gamma$ point.

This overlap is absent at the second magic angle, $\alpha_2 \approx2.221$, where the two channels are clearly separated. Thus, the first magic angle corresponds to a special balance point between the chiral orbital current channel and the relaxation-induced inter-sublattice current channel. This balance is lost once one moves to higher-order magic angles.

This result suggests a useful analogy with heavy-fermion systems, where an itinerant channel hybridizes with a more localized one. In the present case, $J_{AA}$ plays the role of the itinerant current channel inherited from the chiral model, whereas $J_{AB}$ behaves as a more localized, relaxation-induced channel that is activated only in the non-chiral model. The overlap between both channels near the $\Gamma$ point at the first magic angle suggests a natural mechanism for the formation of strongly correlated narrow bands. In this sense, the relaxation-induced channel can be viewed as an effective localized component embedded within the itinerant topological flat-band structure, and the overlap between both channels near high-symmetry points is reminiscent of a Kondo-like hybridization mechanism [63, 64]. We stress that this is an analogy at the level of channel hybridization, but it provides a suggestive framework for interpreting why the first magic angle is particularly favorable for correlated phenomena.

The role of higher-energy bands is also crucial. Panels (c) and (d) of Figure 10 show the contributions of the three lowest bands to $J_{AA}$ and $J_{AB}$ at the first magic angle. For the chiral current $J_{AA}$, the flat band dominates near the $\Gamma$ point, while the higher-energy bands give only relatively small corrections. By contrast, the non-chiral channel $J_{AB}$ receives substantially stronger contributions from the remote bands. This shows that remote-band hybridization selectively enhances the non-chiral inter-sublattice channel.

This distinction is important for understanding the breakdown of the $3/2$ flux quantization rule. In the ideal chiral theory, the flat-band condition is protected by a delicate balance encoded in the chiral current channel. Once finite $\kappa$ is introduced, the additional non-chiral channel $J_{AB}$ grows through hybridization with remote bands and progressively disrupts that balance. The first magic angle remains special because the overlap between $J_{AA}$ and $J_{AB}$ near $\Gamma$ still allows a balanced low-energy regime. At higher-order magic angles, however, this overlap is lost and the non-chiral channel becomes increasingly dominated by remote-band processes. The resulting imbalance provides a direct microscopic mechanism for the breakdown of the chiral $3/2$ flux quantization rule in realistic non-chiral twisted bilayer graphene.

Conclusions

In this work, we analyzed why the first magic angle of twisted bilayer graphene remains robust beyond the ideal chiral limit, while higher-order magic angles are strongly destabilized by realistic same-sublattice tunneling. We worked within the non-chiral Bistritzer–MacDonald continuum Hamiltonian, where lattice relaxation is incorporated through the finite tunneling ratio $\kappa= w_{AA}/w_{AB}$. By extending the squared-Hamiltonian formulation previously developed for the chiral model, we obtained an effective description in which the role of finite $\kappa$ can be separated into diagonal confinement renormalization and off-diagonal inter-sublattice coupling terms.

A central result is that finite $\kappa$ does not simply perturb the chiral flat-band problem. Instead, it reorganizes the effective confinement landscape of the squared Hamiltonian. The diagonal confinement potential is renormalized according to

with an oscillatory component proportional to $-1 + 2\kappa^2$. Therefore, the special value

cancels the spatial modulation of the symmetric confinement term. This value is strikingly close to the relaxation-renormalized value $\kappa\simeq0.7$ commonly used for realistic TBG. Thus, realistic lattice relaxation places the continuum model close to a special confinement point of the squared Hamiltonian, where the symmetric confinement becomes nearly uniform and the remaining diagonal spatial structure is controlled mainly by the antisymmetric pseudo-magnetic term. This provides a direct analytical interpretation of the relaxation-renormalized BM parameter and suggests that lattice relaxation favors a nearly optimized confinement regime.

This confinement mechanism has clear spectral and wave-function consequences. In the chiral limit, the magic-angle hierarchy is governed by an exact scaling structure that produces ideal spectral squeezing and exactly flat bands at the TKV magic angles. Once finite $\kappa$ is included, this scaling structure is partially broken. The first magic angle remains narrow and continuously connected to the chiral flat-band regime, but the higher-order magic angles develop substantial residual dispersion and stronger hybridization with remote bands. The loss of ideal squeezing is therefore not a numerical accident; it is the spectral consequence of the finite-$\kappa$ terms that compete with the chiral cancellation mechanism.

The real-space and Fourier-space analyses support the same conclusion. Finite $\kappa$ squeezes the low-energy wave functions around the AA regions. In reciprocal space, this appears as a shift of the dominant Fourier components toward larger moiré momenta, especially at higher-order magic angles. Therefore, lattice relaxation simultaneously enhances real-space localization and increases coupling to remote bands. This dual effect explains why isolated-flat-band descriptions remain most reliable near the first magic angle but become progressively less controlled for higher-order chiral magic angles.

The expectation-value analysis identified the energetic channels responsible for this behavior. The kinetic energy, the renormalized confinement energy, the antisymmetric pseudo-magnetic contribution, and the off-diagonal interlayer terms all contribute to the balance of the squared energy. Near $\kappa\simeq0.7$, the system lies close to the special uniform-confinement point, while the chiral and non-chiral interlayer contributions become comparable in magnitude. This is a key result: realistic TBG is not merely a weakly broken version of the chiral model, but a system in which the relaxation-renormalized confinement and the interlayer current-like channels enter a balanced regime.

In the squared-Hamiltonian language, the two relevant interlayer channels are the chiral orbital channel $J_{AA}$ and the relaxation-induced inter-sublattice channel $J_{AB}$. The first is inherited from the chiral TKV model, while the second is activated only by finite $\kappa$. These two contributions exhibit a pronounced complementary behavior: their individual values oscillate strongly as a function of $\alpha$, but their sum remains comparatively smooth and follows an approximate quadratic scaling. This current-balance relation is not an exact operator conservation law, but it is an emergent energetic compensation in the low-energy states. It reveals how the non-chiral model redistributes the current-like spectral weight that, in the chiral limit, is carried only by the orbital channel.

This redistribution also clarifies the breakdown of the $3/2$ flux-quantization rule at higher-order magic angles. At the first magic angle, the two current-like channels overlap near the $\Gamma$ point, producing a balanced low-energy regime. At the second and higher magic angles, this overlap is lost. Moreover, remote bands contribute more strongly to the non-chiral channel $J_{AB}$ than to the chiral channel $J_{AA}$. As a result, hybridization with remote bands selectively enhances the relaxation-induced inter-sublattice channel and destabilizes the delicate chiral balance required for the higher-order magic-angle hierarchy.

The main physical picture that emerges is therefore the following. Lattice relaxation drives realistic TBG close to a special confinement regime of the squared Hamiltonian, while finite same-sublattice tunneling activates an additional inter-sublattice channel that competes with the chiral orbital current. At the first magic angle these effects remain balanced, allowing narrow central bands and preserving important remnants of the chiral theory. At higher-order magic angles the balance becomes unstable: remote-band hybridization grows, the current-channel compensation becomes oscillatory, and the exact chiral squeezing mechanism is lost.

The squared-Hamiltonian framework developed here provides a compact way to understand these effects from a single-particle perspective. It connects lattice relaxation, pseudo-magnetic confinement, real-space localization, remote-band hybridization, and interlayer current redistribution within one effective description. More broadly, this approach suggests that realistic moiré systems beyond the chiral approximation can be understood in terms of emergent gauge-field structures and symmetry-constrained current-like channels, rather than only through numerical band-structure deformation.

Acknowledgements

This work was supported by CONAHCyT project 1564464 and UNAM DGAPA project IN101924. L.A.N.L. was supported by a CONAHCyT PhD scholarship No. 1564464. P.A.P. acknowledges support from Grant No. JSF-24-05-0002 of the Julian Schwinger Foundation for Physics Research.

Appendix

Appendix A: Perturbative treatment of the non-chiral term

In this appendix we derive the leading analytical correction produced by the finite same-sublattice tunneling term in the small-$\alpha$ regime. The purpose is twofold.

First, the calculation verifies the numerical results at small twist-angle coupling. Second, it explains why the first magic angle remains comparatively robust beyond the chiral limit: the projected first-order correction vanishes, and the leading non-chiral shift appears only at order $\kappa^2\alpha^2$.

We decompose the non-chiral BM Hamiltonian as

with

Here $H_D$ is the chiral BM Hamiltonian, while $H_M$ contains the finite same-sublattice tunneling generated by lattice relaxation. The block $M$ is the same operator defined in (9). In the notation of the main text, it is diagonal in the $(\Psi,\chi)$ sublattice-block structure and off diagonal in layer space,

Therefore $H_M$ scales as $\kappa\alpha$.

At the $\Gamma$ point, $\mathbf{k}=\mathbf{q}_1$, the chiral zero-mode wave function used in the perturbative calculation can be written as

where

The label $\mu_\alpha=\pm1$ distinguishes the two symmetry-related zero-mode branches, and

is the normalization factor in the first-star approximation. Figure 11 and Figure 12 show that this perturbative description is reliable only at small $\alpha$. Near higher-order magic angles, the expectation values of $H_D$ and $H_M$ become strongly mixed, consistent with the breakdown of the simple chiral scaling picture discussed in the main text.

Vanishing of the first-order correction

The first-order correction is the matrix element of $H_M$ projected onto the unperturbed chiral zero-mode subspace,

Using Eqs. (81) and ((82)), this becomes

Substituting the explicit zero-mode components gives, to leading order in the first-star expansion,

The unit-cell integral of a product of moiré harmonics is nonzero only if the total momentum carried by the product is a reciprocal lattice vector of the moiré superlattice. Equivalently,

only when

For the $\Gamma$-point chiral zero mode, the products appearing in Eq. (86) do not contain a $C_3$-invariant zero-momentum component. The three first-star contributions are cyclically related by $C_3$, and their phases cancel in the unit-cell average. Thus the projection of $H_M$ onto the chiral zero-mode subspace vanishes,

This is the key perturbative result. Finite same-sublattice tunneling is present already at order $\kappa\alpha$ in the Hamiltonian, but its first-order energy correction to the $\Gamma$-point chiral zero mode vanishes by the momentum-selection and $C_3$-symmetry structure of the unperturbed state.

Leading second-order correction

Because the first-order correction vanishes, the leading non-chiral shift is second order in $H_M$. Standard Schrödinger perturbation theory gives

where $\left|\Phi_n^{(0)}\right\rangle$ are eigenstates of the chiral Hamiltonian $H_D$. Since $H_M \propto\kappa\alpha$, (90) immediately implies the leading scaling

This scaling is the perturbative origin of the robustness of the first magic angle: for $\alpha\lesssim\alpha_1$, the finite-$\kappa$ correction remains controlled, whereas for higher-order magic angles the same correction is amplified by the larger value of $\alpha$ and by stronger remote-band hybridization.

To obtain an explicit analytical approximation, we introduce the projector $Q$ onto the remote-band subspace,

and replace the weighted remote-band denominators by an effective gap $\Delta E_{\mathrm{eff}}$. This gives the closure approximation

Using (89), this reduces to

In the dimensionless BM units used throughout the paper, the leading remote-band denominators at small $\alpha$ are of order unity. Therefore the dominant dependence of (94) is controlled by the matrix element of $H_M^2$, and the leading scaling remains $\kappa^2\alpha^2$.

The square of the relaxation block is

Using $U_0(-\mathbf{r}) = U_0^\dagger(\mathbf{r})$, the diagonal product defines the scalar moiré field

Thus

The density of the unperturbed state contains the harmonic expansion

with an analogous expression for $\lvert\psi_2^{(0)}\rvert^2$. After integration over the moiré unit cell, only the zero-momentum harmonic components survive. The resulting overlap coefficients determine the analytical numerator of the perturbative correction. In the gauge and first-star convention used here, the explicit harmonic evaluation gives the compact structure

where the $\alpha^2$ term is the leading non-chiral correction and the $\alpha^3$ term is the first higher-harmonic correction within the same expansion. The factor $\mu_\alpha= \pm1$ encodes the branch dependence of the $\Gamma$-point chiral zero mode.

To compare directly with the exact numerical data, we keep the leading normalization and overlap corrections generated by the same harmonic expansion. This gives the effective normalization factor

The numerical coefficients in Eq. (100) are fixed by the unit-cell harmonic overlaps of the $U_m$ functions in the convention of Eq. (82); they are not fitting parameters. Combining Eqs. (99) and (100) gives the analytical approximation

This expression makes explicit that the leading non-chiral correction is controlled by $\kappa^2\alpha^2$, with higher-order corrections entering at larger $\alpha$.

Figure 13 compares (101) with the numerical calculation. The agreement is excellent for $\alpha\lesssim0.3$, confirming the validity of the perturbative expansion in the small-$\alpha$ regime.

Deviations at larger $\alpha$ are expected and physically meaningful: higher harmonics, remote-band hybridization, and the non-chiral current channel become increasingly important as the system approaches the first magic angle and beyond.

The important conclusion is that the finite-$\kappa$ correction is not linear in $\kappa\alpha$ at the $\Gamma$-point chiral zero mode. The first-order term vanishes by symmetry, and the leading correction enters as

This result explains why the first magic angle preserves several features of the chiral model despite the realistic value $\kappa\simeq0.7$. At the same time, it also explains why higher-order magic angles are much less robust: increasing $\alpha$ pushes the system out of the perturbative regime, where confinement reshaping, current-channel redistribution, and remote-band hybridization dominate the spectrum.