# Mathematical models of topologically protected transport in twisted bilayer
graphene

Guillaume Bal Departments of Statistics and Mathematics and CCAM, University
of Chicago, Chicago, IL 60637<EMAIL_ADDRESS>Paul Cazeaux
Department of Mathematics, Virginia Tech, Blacksburg, VA 24060<EMAIL_ADDRESS>Daniel Massatt  Department of Mathematics, Louisiana State University, Baton
Rouge, LA 70803<EMAIL_ADDRESS>Solomon Quinn Department of Statistics and
CCAM, University of Chicago, Chicago, IL 60637<EMAIL_ADDRESS>
###### Abstract

Twisted bilayer graphene gives rise to large moiré patterns that form a
triangular network upon mechanical relaxation. If gating is included, each
triangular region has gapped electronic Dirac points that behave as bulk
topological insulators with topological indices depending on valley index and
the type of stacking. Since each triangle has two oppositely charged valleys,
they remain topologically trivial.

In this work, we address several questions related to the edge currents of
this system by analysis and computation of continuum PDE models. Firstly, we
derive the bulk invariants corresponding to a single valley, and then apply a
bulk-interface correspondence to quantify asymmetric transport along the
interface. Secondly, we introduce a valley-coupled continuum model to show how
valleys are approximately decoupled in the presence of small perturbations
using a multiscale expansion, and how valleys couple for larger defects.
Thirdly, we present a method to prove for a large class of continuum
(pseudo-)differential models that a quantized asymmetric current is preserved
through a junction such as a triangular network vertex.

We support all of these arguments with numerical simulations using spectral
methods to compute relevant currents and wavepacket propagation.

##  1 Introduction

Twisted bilayer graphene (tBLG) is widely studied for its unique mechanical
and electronic properties including the magic angle superconductivity [13,
10]. Upon gating, it acts as host of a network of topological interface
channels, which can be seen experimentally and theoretically [32, 31, 14, 1].
tBLG is constructed by taking two periodic 2D sheets of graphene and stacking
them with a relative twist, typically small. The atoms relax to minimize
energy, forming large triangular regions of the energetically favorable AB and
BA Bernal stacking [33, 14, 15]. It is relevant to note this asymmetric
transport under gating is a separate phenomena from magic angle
superconductivity. Indeed, the asymmetric transport phenomena only requires
sufficiently small twist angles and vertical gating, i.e. by inducing a
potential difference between the two layers, while superconductivity can only
occur precisely at the magic twist angles. In this section, we make use of
Figure 1 for discussion of the geometry.

The interior of one of these triangles can be considered as approximately an
infinite periodic material as the side of one of these triangles scales
inversely proportional to the small twist angle [14]. We note that in this
work we will consider only continuum models and not lattice models, but for
understanding of the purpose and origins of the continuum models in this work,
it is useful to see the lattice model origins. If $\mathcal{L}$ is the lattice
matrix corresponding to the periodicity of the material, we define the
reciprocal lattice unit cell as $\Gamma^{*}:=2\pi\mathcal{L}^{-T}[0,1)^{2}$.
Consider periodic Hamiltonian $H$ describing the infinitely extended bulk
corresponding to one of the triangular regions. Its spectra, which in turn
describes electronic properties, can be described by the Bloch states
satisfying

$H(\xi)\psi^{(j)}(\xi)=E^{(j)}(\xi)\psi^{(j)}(\xi),\qquad\xi\in\Gamma^{*}$
(1.1)

where $H(\xi)=e^{-i\xi\cdot r}He^{i\xi\cdot r}$. This expression is understood
as operator composition. We assume $E^{(1)}(\xi)\leq E^{(2)}(\xi)\leq\cdots$
are the ordered eigenvalues (or energies). Gating opens a band gap at $E=0$,
meaning all $E^{(j)}(\xi)$ are bounded away from $0$. A topological integer
index can be associated with the operator $H$ by the formula

$I[H]:=\sum_{j\geq
n}\frac{i}{2\pi}\int_{\Gamma^{*}}d(\psi^{(j)},d\psi^{(j)}).$ (1.2)

Here $n$ is the smallest index such that $E^{n}(\xi)>0$ for all
$\xi\in\Gamma^{*}$. Many of these more complex Hamiltonians can be well
approximated by simpler continuum models that capture the relevant dispersion
relations [9, 20, 36]. In the case of a bilayer graphene, $I[H]$ can be
computed by expanding only around the Dirac points $K,K^{\prime}\in\Gamma^{*}$
as they possess the topological information [9, 17, 20]. Such effective moiré-
scale continuum Hamiltonians for low-energy physics of bilayer graphene, in
particular the Bistritzer-MacDonald model [10], have been recently derived in
the small twist angle regime from atomic-scale models: in [12], using a formal
variational approximation of the tBLG Kohn-Sham Hamiltonian, and in [35], by
identifying an appropriate asymptotic parameter regime of the discrete tight-
binding model. A local in momentum Hamiltonian $H(\xi)$ can be constructed for
each of these valleys through a $4\times 4$ model as follows:

$H(\xi):=\begin{pmatrix}\Omega
I+\xi_{1}\sigma_{1}+\eta\xi_{2}\sigma_{2}&\lambda U^{*}\\\ \lambda U&-\Omega
I+\xi_{1}\sigma_{1}+\eta\xi_{2}\sigma_{2}\end{pmatrix}.$ (1.3)

Here $\eta=1$ for valley $K$ and $\eta=-1$ for valley $K^{\prime}$,
${{{2}}}\Omega\in\mathbb{R}$ measures the electrostatic potential difference
between the two layers due to vertical gating, and $U\in\\{A,A^{*}\\}$ models
an interlayer coupling term with strength $0\not=\lambda\in\mathbb{R}$, where
$A=\begin{pmatrix}0&1\\\ 0&0\end{pmatrix}$ [26]: $A^{*}$ and $A$ correspond to
whether the triangle has $AB$ or $BA$ stacking geometry respectively. This
Hamiltonian has been rescaled to ensure Fermi velocity is $1$ and unitless,
i.e. there is no prefactor in front of the Dirac term
$\xi_{1}\sigma_{1}+\xi_{2}\sigma_{2}$. A standard result in topological
insulators is the bulk-interface correspondence [9, 20], which states that if
two topological materials are glued together at an interface, then there is
quantized edge transport given by the difference of the topological index on
either side of the interface; see [5, 6, 17, 29] for the derivation of the
correspondence for differential models. In bilayer graphene the $K$ and
$K^{\prime}$ valleys contribute opposite signed topological charge, so gated
bilayer graphene is technically trivial. However, it is understood that the
two valleys in many settings decouple, with each valley Hamiltonian having
possibly nontrivial topology, which is the basis of valleytronics [37, 38].

In this work, we address gated tBLG via numerics and analysis of continuum
models built to address three different questions corresponding to the two
regions highlighted in Figure 1. For the first question addressed in Section
2, we use the $4\times 4$ single-valley bulk model, compute its topological
invariant, and then using the corresponding interface model apply the bulk-
interface correspondence at an edge between two triangles (region $1$). In
particular, we make use of the bulk-difference invariant [5] between AB and BA
stacked graphene for a single valley to find a quantized current of $\pm 2$.
In Section 3, we focus on analyzing the separation of the $K$ and $K^{\prime}$
valley at the interface (region $1$). Since the $4\times 4$ bilayer graphene
model is only accurate near the valleys, we study valley coupling with a toy
$2\times 2$ model that exhibits two conical Dirac points. To model defects, we
consider a two-scale model to introduce fast small fluctuations, and show the
valley coupling is weak in this regime. We use numerics and scattering theory
to further show that as the defects are tuned to couple the valleys, then
valley separation breaks and quantized edge current is lost. In Section 4, we
define a conductivity that quantifies the asymmetric transport through a
junction (e.g. region $2$). Following pseudodifferential calculus arguments
from [29], we derive an analytic formula for this junction conductivity and
show that the latter is immune to perturbations. This allows for explicit
evaluations of the conductivity for both the $4\times 4$ bilayer graphene and
$2\times 2$ Dirac models. Finally in Section 5, we include numerical
verification of results from the preceding sections using spectral methods to
properly capture the Dirac operator. Direct computation of traces is used to
compute currents either isolated to the edge (Region 1) and currents through a
junction (Region 2), and wavepacket propagation is used to illustrate the
connection between the two regions by evolving a state localized on an edge
through a network junction.

Figure 1: Triangular domains form in mechanically relaxed tBLG. We highlight
two regions of interest for our analysis. Region $1$ corresponds to an edge,
where a standard bulk-interface correspondence will be studied as well as a
continuum valley coupling model. Region 2 is at a junction, where we will
analyze current conservation.

##  2 Gated tBLG model

In this section, we analyze region 1 from Figure 1 and use the bulk-interface
correspondence to show quantized asymmetric transport for a single valley. To
begin, we consider the bulk Hamiltonian in (1.3) described by the continuum
model

$H=\begin{pmatrix}\Omega I+D\cdot\sigma^{(\eta)}&\lambda U^{*}\\\ \lambda
U&-\Omega I+D\cdot\sigma^{(\eta)}\end{pmatrix}.$

Here $\sigma^{(\eta)}:=(\sigma_{1},\eta\sigma_{2})$, and
$\xi\cdot\sigma^{(\eta)}:=\xi_{1}\sigma_{1}+\eta\xi_{2}\sigma_{2}$ where we
recall $\eta\in\pm 1$ is the valley index. We will sometimes denote
$\sigma^{(1)}$ by $\sigma$. This operator in two-dimensional physical
variables is an unbounded operator on $L^{2}(\mathbb{R}^{2};\mathbb{C}^{4})$.
The operator in the Fourier variables is given as in (1.3) for each
$\xi\in\mathbb{R}^{2}$ by the $4\times 4$ matrix:

$H(\xi)=\begin{pmatrix}\Omega I+\xi\cdot\sigma^{(\eta)}&\lambda U^{*}\\\
\lambda U&-\Omega I+\xi\cdot\sigma^{(\eta)}\end{pmatrix}.$

These models are defined on a plane $\mathbb{R}^{2}$, and hence are not
defined over a compact manifold. To construct topological invariants, we
introduce as in [5] the bulk-difference invariant describing the difference
between AB and BA stacked graphene, the two bulks enclosing an interface. The
two bulk Hamiltonians are given by

$\displaystyle H_{+}(\xi)=\begin{pmatrix}\Omega
I+\xi\cdot\sigma^{(\eta)}&\lambda A^{*}\\\ \lambda A&-\Omega
I+\xi\cdot\sigma^{(\eta)}\end{pmatrix},$ $\displaystyle
H_{-}(\xi)=\begin{pmatrix}\Omega I+\xi\cdot\sigma^{(\eta)}&\lambda A\\\
\lambda A^{*}&-\Omega I+\xi\cdot\sigma^{(\eta)}\end{pmatrix}.$

These $4\times 4$ Hermitian matrices may be diagonalized with eigenvectors
denoted by $\psi_{\pm}^{(j)}(\xi)$ for $1\leq j\leq 4$. The eigenprojectors
$(\pm,\xi)\to\Pi^{(j)}_{\pm}(\xi):=|\psi_{\pm}^{(j)}(\xi)\rangle\langle\psi_{\pm}^{(j)}(\xi)|$
are uniquely defined as the eigenvalues are all simple
$E_{1}(\xi)<E_{2}(\xi)<0<E_{3}(\xi)<E_{4}(\xi)$ and independent of the phase
$\pm$ as we will show in the appendix. Since $\xi\in\mathbb{R}^{2}$ is simply
connected, the normalized eigenvectors $\psi_{\pm}^{(j)}(\xi)$, a priori
defined up to a multiplicative phase, may be chosen smoothly in $\xi$; an
explicit expression is presented in Appendix A. Associated to each projector
$\Pi$ is a curvature defined by ${\rm tr}\,\Pi\,d\Pi\wedge d\Pi$. We introduce
the following integrals

$W_{\pm}^{j}:=\frac{i}{2\pi}\int_{\mathbb{R}^{2}}{\rm
tr}\,\Pi^{(j)}_{\pm}d\Pi^{(j)}_{\pm}\wedge
d\Pi^{(j)}_{\pm}=\frac{i}{2\pi}\int_{\mathbb{R}^{2}}d(\psi_{\pm}^{(j)},d\psi_{\pm}^{(j)})=\frac{i}{2\pi}\lim_{R\to\infty}\displaystyle\int_{{\mathbb{S}}_{R}^{1}}(\psi_{\pm}^{(j)},d\psi_{\pm}^{(j)}).$
(2.1)

Here, ${\mathbb{S}}_{R}^{1}$ is the circle of radius $R$ with counterclockwise
orientation. The second equality can be obtained by explicit calculations. The
third equality is a consequence of the Stokes theorem and the smoothness of
the vector field $\psi_{\pm}^{(j)}(\xi)$.

We will show that $H_{\pm}(\xi)$ are gaped at energy $E=0$. As a consequence,
we expect $W_{\pm}:=W_{\pm}^{3}+W_{\pm}^{4}$ corresponding to the integral of
the curvature of the rank-two projectors $\Pi_{\pm}=\chi(H_{\pm}>0)$ to have
topological significance. A main difference with (1.2), however, is that
unlike $\Gamma^{*}$, $\mathbb{R}^{2}$ is not compact. When $H$ is a Dirac
operator for instance, we observe that
$W_{\pm}\in\\{-\frac{1}{2},\frac{1}{2}\\}$ is not an integer [5]. This
motivates the definition of the bulk-difference invariant

$W:=W_{+}-W_{-}=W_{+}^{3}+W_{+}^{4}-(W_{-}^{3}+W_{-}^{4}).$ (2.2)

The invariant is modeling a transition between the bulk Hamiltonians $H_{+}$
and $H_{-}$. We will show that $W$ is indeed an integral-valued invariant
unlike its components $W_{\pm}^{j}$. Heuristically, this shows that it is
easier to define a phase transition between insulators rather than absolute
phases, which are often ill-defined for partial differential models such as
Dirac equations.

We will compare this bulk invariant to an edge current, which describes
asymmetric transport along an interface separating the bulk phases and is
defined next. Firstly, we introduce an edge model by making the bulk above
into a continuum model with an interface:

$H_{e}:=\begin{pmatrix}\Omega I+D\cdot\sigma^{(\eta)}&\lambda U^{*}(y)\\\
\lambda U(y)&-\Omega I+D\cdot\sigma^{(\eta)}\end{pmatrix}$ (2.3)

with $U(y)=\frac{1}{2}((1+m(y))A+(1-m(y))A^{*})$ and
$D=(D_{x},D_{y})=\frac{1}{i}(\partial_{x},\partial_{y})$. We choose $m$ such
that for large $y$ we obtain BA stacking, while for large negative $y$ we
obtain AB stacking. In order to do this, we define the switch function space
$\bar{\mathfrak{S}}(c_{1},c_{2};y_{1},y_{2})$ as follows. Given
$f_{1},f_{2}\in\mathbb{C}$ and $c_{1},c_{2}\in\mathbb{R}$ with $c_{1}\leq
c_{2}$, we say that a function $f:\mathbb{R}\rightarrow\mathbb{R}$ belongs to
$\bar{\mathfrak{S}}(f_{1},f_{2};c_{1},c_{2})$ if and only if

$\displaystyle f(y)=\begin{cases}f_{1},&y<c_{1}\\\ f_{2},&y>c_{2}\end{cases}.$
(2.4)

We then define the space of _smooth switch functions_ by
$\mathfrak{S}(f_{1},f_{2};c_{1},c_{2})=\bar{\mathfrak{S}}(f_{1},f_{2};c_{1},c_{2})\cap\mathcal{C}^{\infty}(\mathbb{R})$.

We assume $m\in\mathfrak{S}(-1,1;-y_{0},y_{0})$ for some $y_{0}>0$. We define
an edge state capturing energy function
$\varphi\in\mathfrak{S}(0,1;-E_{0},E_{0})$ where $E_{0}>0$ is chosen to be
smaller than the bulk band gap. We have not yet shown there is a bulk band
gap, but prove this below in Theorem 2.1. We note that $\varphi^{\prime}$ is
supported in the gap, and thus $\varphi^{\prime}(H_{e})$ captures edge state
spectra only. Finally we define $P\in\mathfrak{S}(0,1;-x_{0},x_{0})$ for some
$x_{0}>0$, which in the limit $x_{0}\to 0$ approximates the Heaviside
function. The rate of propagation of current through the edge of $P$ is given
by [3, 5, 19, 27]

$\sigma_{I}:=\text{Tr}\,i[H_{e},P]\varphi^{\prime}(H_{e}).$ (2.5)

The main objective of this section is to show that
$2\pi\sigma_{I}=W=-2\eta\operatorname{\text{sign}}{\Omega}$. The first
equality is a bulk-interface correspondence, which has been derived in a
variety of contexts [3, 5, 19, 27, 29], and in particular in [5, 29] for the
operators of interest here. The advantage of such a correspondence is that the
explicit computation of $W$ is often simpler than that of $\sigma_{I}$. We
will obtain the former by adapting results derived in [8].

The Fourier transform of $H_{e}$ can be taken along the interface, obtaining
the Hamiltonian

$H_{e}(\xi_{1})=\begin{pmatrix}\Omega I+\xi_{1}\sigma_{1}+\eta
D_{y}\sigma_{2}&\lambda U^{*}(y)\\\ \lambda U(y)&-\Omega
I+\xi_{1}\sigma_{1}+\eta D_{y}\sigma_{2}\end{pmatrix}.$ (2.6)

In Kronecker product notation, this yields

$H_{e}(\xi_{1})=\Omega\sigma_{3}\otimes
I_{2}+I_{2}\otimes(\xi_{1}\sigma_{1}+\eta
D_{y}\sigma_{2})+\frac{1}{2}(\sigma_{1}\otimes\sigma_{1}+m(y)\sigma_{2}\otimes\sigma_{2}).$
(2.7)

In Figure 3, we plot the bulk band structure and superimpose the eigenvalues
within the bulk band gap of a periodized version of $H_{e}(\xi_{1})$, where we
only include eigenvalues weighted towards one interface. This is necessary as
the periodization creates a second domain wall with oppositely directed edge
states. We then see that there are two right-ward propagating states, as
expected (two curves crossing the bulk band gap with positive group velocity
$\partial E/\partial\xi$). We note that the symmetry in the edge states
$E(\xi_{1})=-E(-\xi_{1})$ comes from the following symmetry operation. Let
$K=\sigma_{2}\otimes\sigma_{1}$. Then

$K^{-1}H_{e}(\xi_{1})K=-H_{e}(-\xi_{1}).$

To prove that $W$ is an invariant, the first step will be verifying there is
indeed a bulk gap, which turns out to have a ring shape. We then use this to
calculate the bulk-difference invariant $W$. See Figure 2 to see concentration
of curvature on the ring and the corresponding band structure. Combining the
above yields our main result:

Figure 2: Here we plot the eigenvalues of $H_{+}(\xi)$, i.e. the band
structure on the left. On the right, we plot the curvature, i.e., the
integrand of $W_{+}^{3}+W_{+}^{4}$ in (2.1). Here we set $\lambda=0.2$ and
$\Omega=1$. As expected, curvature concentrates on the line $|\xi|=1$ for
$|\lambda|$ small.  Figure 3: Bulk band structure (filled) is plotted along
with calculated edge state eigenvalues (points) of a periodized
$H_{e}(\xi_{1})$. In particular, the plotted interface eigenvalues are the
nearest two eigenvalues to $E=0$.

######  Theorem 2.1.

The bulk system $H$ has a band gap for all $\lambda\neq 0$. Further,

${{{-}}}W_{\pm}=\pm\eta\,\text{sign}(\Omega).$ (2.8)

The bulk-interface correspondence then gives us

$2\pi\sigma_{I}=W={{{-}}}2\eta\,\text{sign}(\Omega).$ (2.9)

Figure 4:  $K$ and $K^{\prime}$ valley current directions are shown for each
edge. Note that the currents per valley are always moving in opposite
directions.

The proof of this result is postponed to Appendix A. The proof provides a more
detailed expression for the components describing $W$. In particular, assume
$\Omega>0$ and $\lambda>0$ for concreteness and define
$\beta=2\Omega+\sqrt{4\Omega^{2}+\lambda^{2}}$. We will show that

$-W^{3}_{+}=+W^{3}_{-}=\frac{\lambda^{2}+3\beta^{2}}{2(\lambda^{2}+\beta^{2})},\quad-W^{4}_{+}=+W^{4}_{-}=\frac{\lambda^{2}-\beta^{2}}{2(\lambda^{2}+\beta^{2})}.$

This is a surprising result. If $W^{j}_{\pm}$ involved an integration over a
compact domain such as $\Gamma^{*}$, then each of these numbers would be
guaranteed to be an integer (a Chern number). We already mentioned that for
the Dirac operator $W_{\pm}=\pm\frac{1}{2}$. Here, however, the terms
$W^{j}_{\pm}$ take a continuous range of values as $(\Omega,\lambda)$ vary.
Moreover, we show in Appendix A using the gluing procedure proposed in [5]
that $W^{3}_{+}-W^{4}_{-}$ and $W^{4}_{+}-W^{3}_{-}$ are both topological
invariants (Chern numbers) with integral values (here equal to $1$). We thus
obtain the peculiar result that the gluing of the projector of the
fourth/third energy band of the AB Hamiltonian with (an appropriate rotation
of) the projector of the third/fourth energy band of the BA Hamiltonian
defines a topological invariant. The details of this computation are provided
in Appendix A.

The results of Appendix A also provide a spectral gap about $0$ of size
$2E_{0}=\frac{2|\Omega||\lambda|}{\sqrt{4\Omega^{2}+\lambda^{2}}}$. So the
above theorem applies for a spectral density $\varphi^{\prime}(h)$ in (2.5)
supported in the interval $(-E_{0},E_{0})$.

We note that we assumed an ordering to the AB and BA stacking. If we switch
the order, we would swap the sign of $W$, and hence $\sigma_{I}$. Hence each
edge has an edge current dependent on the order of the AB-BA regions and also
the valley index of interest (see Figure 4).

We also note that the structure of the Hamiltonians $H$ and $H_{e}$ are
identical to the so called 2-replica model static Hamiltonian used to
approximate the dynamics of the circularly polarized laser-driven graphene
system, which is a Floquet Topological Insulator [8, 34]. The reason these two
very different physical systems share a common Hamiltonian is that in both
situations the bulk systems are described as a coupling between Dirac cones
shifted relative to each other in energy. The overlapping cones’ intersection
forms a ring, and the coupling between these Dirac cones hybridizes the bands
to form ring-shaped band gaps. We remark that we expect bulk multilayer
stacked graphene with all layers either AB stacked or BA stacked form the
larger replica systems corresponding to a collection of uniformly shifted
Dirac cones. Here the number of replicas or shifted Dirac cones corresponds to
the number of sheets, while in the Floquet setting the number corresponded to
the number of included time Fourier frequencies. In keeping with the Floquet
terminology, we also refer to our Hamiltonians $H$ and $H_{e}$ as a 2-replica
model as it corresponds to two shifted Dirac cones.

##  3 Valley Coupling

As we saw from the previous section, both valley Hamiltonians are
topologically non-trivial. However, if valley coupling is allowed, then the
total asymmetric current summed over the two valleys vanishes making the total
system topologically trivial.

In section 3.1, we introduce a simplified differential operator that captures
the two-valley physics, and construct one-dimensional and two-dimensional
effective valley-coupling models using a multiscale expansion. In section 3.2,
we then define a valley conductivity object. In section 3.3, we use numerics
to calculate currents demonstrating what regime of defects yield good valley
separation, and when does valley decoupling fail. In section 3.4 we use
scattering theory to explain the numeric results from section 3.3.

###  3.1 Valley coupling models

To model a two-valley system and inter-valley coupling, consider first the
bulk Hamiltonian without mass

$H_{b}=\Phi(D_{x})\sigma_{1}+D_{y}\sigma_{2}.$

Here $\Phi\in C(\mathbb{R})$ is a function with exactly two zeros
corresponding to two Dirac points. We will assume $\Phi$ is twice
differentiable at these zeros. In frequency space, we denote the wavenumber
pair $(\xi,\zeta)\in\mathbb{R}^{2}$. We have changed the notation from section
2 as now we wish to emphasize the direction along the edge $\xi$. We have

$H_{b}(\xi,\zeta)=\Phi(\xi)\sigma_{1}+\zeta\sigma_{2},$

which admits band structure. We will assume for simplicity of presentation
throughout this section that the Dirac points are located at frequencies
$(\xi,\zeta)=(\pm 1,0)$.

Figure 5: Band structure of two-valley linear Dirac Hamiltonian,
$\Phi(\xi)=1-|\xi|$.

For example, see Figure 5 for $\Phi(\xi)=|\xi|-1$. Another natural choice for
$\Phi$ leading to a differential operator $H_{b}$ is
$\Phi(\xi)=\frac{1}{2}(\xi^{2}-1)$. In both cases, the valley at Dirac points
$\xi=1$ and $\xi=-1$ have topological indices given by $1$ and $-1$
respectively. We will primarily focus on small highly oscillatory defects as a
model for distribution of small defects.

Consider a semiclassical theory with $0<\varepsilon\ll 1$ scaling near Dirac
points so that the two valleys are separated by wavenumbers of order
$2\varepsilon^{-1}$. We focus on a small energy window near $E=0$ and thus
analyze an equation of the form

$\displaystyle\biggl{(}\Phi(\varepsilon D_{x})\sigma_{1}+\varepsilon
D_{y}\sigma_{2}+\lambda\frac{y}{\varepsilon}\sigma_{3}+\varepsilon
V(x,\frac{x}{\varepsilon},\frac{y}{\varepsilon})-\varepsilon
E\biggr{)}\psi=0.$ (3.1)

In the absence of a domain wall ($\lambda=0$) and perturbing potential
($V=0$), assuming $|\Phi^{\prime}(\pm 1)|=1$, this system has a circular
dispersion relation level set at small energies as is typical for Dirac
points.

We consider a formal derivation in the limit $\varepsilon\ll 1$ of an inter-
valley coupling model whose qualitative predictions will be confirmed by
numerical simulations. While many steps can undoubtedly be justified more
rigorously using methodology such as in [18], we will not do so here.

The domain wall $\lambda\frac{y}{\varepsilon}\sigma_{3}$ is assumed to be
steep with $\lambda\gg 1$ as this prevents coupling between various confined
modes in $y$. Introduce the microscopic variable $Y=y/\varepsilon$ and look
for solutions of the form $\psi(Y=\frac{y}{\varepsilon},x)$. Applied to such
functions, (3.1) becomes

$\displaystyle\big{(}D_{Y}\sigma_{2}+\lambda Y\sigma_{3}+\Phi(\varepsilon
D_{x})\sigma_{1}+\varepsilon V(x,\frac{x}{\varepsilon},Y)-\varepsilon
E\big{)}\psi=0.$ (3.2)

When $\lambda$ is large, then the two first terms dominate in the sense that

$(D_{Y}\sigma_{2}+\lambda
Y\sigma_{3})\phi_{0}(Y)=0,\quad\phi_{0}(Y)=(\lambda\pi)^{-1/4}e^{-\frac{1}{2}\lambda
Y^{2}}\cdot\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\ -1\end{pmatrix},$

while $D_{Y}\sigma_{2}+\lambda Y\sigma_{3}$ applied to any other Hermite
function provides a term of order $|\lambda|\gg 1$ formally much larger than
the term $\Phi(\varepsilon D_{x})\sigma_{1}$.

Decomposing $\psi(Y,x)=\phi_{0}(Y)\psi(x)$ for a scalar function $\psi(x)$
thus neglecting higher-order Hermite functions in $Y$, we obtain from (3.2)
(multiplying by $\phi_{0}^{T}$ and integrating in $Y$) a reduced Hamiltonian
for $\psi(x)$ given by

$\displaystyle\Big{(}-\Phi(\varepsilon D_{x})+\varepsilon
V(x,\frac{x}{\varepsilon})-\varepsilon E\Big{)}\psi=0.$ (3.3)

Here, $V(x,X)=(\phi_{0}(Y),V(x,X,Y)\phi_{0}(Y))_{Y}$. We also used that the
current

$(\phi_{0},\sigma_{1}\phi_{0})_{Y}=-1.$

In both cases, $(\cdot,\cdot)_{Y}$ is the usual inner product on
$L^{2}(\mathbb{R}_{Y};\mathbb{C}^{2})$.

Define now multiscale functions $f(x,X=\frac{x}{\varepsilon})$ in the variable
$x$ so that $D_{x}$ becomes $D_{x}+\varepsilon^{-1}D_{X}$. Then the multiscale
operator is

$\displaystyle-\Phi(D_{X})+(-\Phi^{\prime}(D_{X})D_{x}+V(x,X)-E)\varepsilon+O(\varepsilon^{2}).$
(3.4)

Going to the Fourier domain $X\to K$ we obtain

$\displaystyle-\Phi(K)+(-\Phi^{\prime}(K)D_{x}+\hat{V}(x,K)*-E)\varepsilon+O(\varepsilon^{2}).$
(3.5)

Here $*$ is convolution in the Fourier variable $K$. To leading order we get

$\displaystyle\Phi(K)\psi_{0}=0,\qquad\mbox{ i.e.,
}\quad\psi_{0}=\psi_{+}(x)\delta(K-1)+\psi_{-}(x)\delta(K+1).$ (3.6)

The next order is

$\displaystyle-\Phi(K)\psi_{1}-\Phi^{\prime}(K)D_{x}\psi_{0}+\int\hat{V}(x,K-K^{\prime})\psi_{0}(x,K^{\prime})dK^{\prime}-E\psi_{0}=0.$
(3.7)

We thus find the coupled equations

$\displaystyle-\Phi^{\prime}(1)D_{x}\psi_{+}+\hat{V}(x,0)\psi_{+}+\hat{V}(x,2)\psi_{-}$
$\displaystyle=E\psi_{+}$ (3.8)
$\displaystyle-\Phi^{\prime}(-1)D_{x}\psi_{-}+\hat{V}(x,0)\psi_{-}+\hat{V}(x,-2)\psi_{+}$
$\displaystyle=E\psi_{-}.$ (3.9)

This generates inter-valley coupling if by valleys we mean the wavenumbers
close to $K=\pm 1$. Letting

$\displaystyle
H_{2}=\begin{pmatrix}-\Phi^{\prime}(1)D_{x}+\hat{V}(x,0)&\hat{V}(x,2)\\\
\hat{V}(x,-2)&-\Phi^{\prime}(-1)D_{x}+\hat{V}(x,0)\end{pmatrix},$ (3.10)

we obtain the eigenproblem $H_{2}\psi=E\psi$, $\psi=(\psi_{+},\psi_{-})^{T}$,
as a one-dimensional model incorporating valley coupling.

We now build a two-dimensional Hamiltonian that incorporates a domain wall and
valley coupling. Generalizing the above Hamiltonian to include a $y$-dependent
potential,

$\displaystyle
H_{2}^{y}=\begin{pmatrix}-\Phi^{\prime}(1)D_{x}+\hat{V}(x,0,y)&\hat{V}(x,2,y)\\\
\hat{V}(x,-2,y)&-\Phi^{\prime}(-1)D_{x}+\hat{V}(x,0,y)\end{pmatrix},$ (3.11)

we then define the $4\times 4$ system

$H_{4}=D_{y}I_{2}\otimes\sigma_{2}+\lambda
yI_{2}\otimes\sigma_{3}+H_{2}^{y}\otimes\sigma_{1}.$ (3.12)

We choose this Pauli matrix representation such that $H_{4}$ is consistent
with $H_{2}$ when $V$ is $y$-independent. In particular, for $e_{\pm}$ being
the standard basis of $\mathbb{R}^{2}$ corresponding to valley index,
$e_{\pm}\otimes\phi_{0}(y)$ are still in the kernel of the domain wall terms
(first two terms), and $\phi_{0}(y)$ is an eigenfunction of $\sigma_{1}$,
meaning when $V$ is $y$-independent the edge state eigenfunctions of $H_{2}$
naturally generalize to eigenfunctions of $H_{4}$.

The operators in (3.10) and (3.12) are used for the numerical simulations of
valley coupling effects.

######  Remark 1.

Instead of $\varepsilon V(x,\frac{x}{\varepsilon})$, we could consider a
larger-amplitude lower-frequency perturbation $V(x)$ in (3.1), another form of
defect studied numerically in [28, 38]. Let us make a few brief remarks on
this setting. In the semiclassical limit of large valley separation of order
$2\varepsilon^{-1}$, we derive as above the one-dimensional semiclassical
problem

$\displaystyle(-\Phi(\varepsilon D_{x})+V(x)-E)\psi=0.$ (3.13)

Consider the well-studied case $\Phi(\xi)=1-\xi^{2}$ for concreteness and see,
e.g., [16] for an analysis of such semiclassical problems.

We may find to arbitrary algebraic accuracy in powers of $\varepsilon$
solutions of the form

$\displaystyle\psi(x)=\sum_{\pm}e^{iS_{\pm}(x)/\varepsilon}a_{\pm}(x;\varepsilon),\quad-\Phi(\partial_{x}S_{\pm})+V(x)-E=0.$
(3.14)

The latter eikonal equation for $S_{\pm}$ admits two solutions

$\displaystyle S_{\pm}(x)=\pm\int_{0}^{x}(\Phi(0)+E-V(z))^{\frac{1}{2}}dz.$
(3.15)

We assume $\Phi(0)+E-V(0)>0$ so that propagating modes exist at $x_{0}$ with
local wavenumber given by $\xi=\sqrt{\Phi(0)+E-V(0)}$ at the energy of
interest $E$. Assuming $\Phi(0)+E-V(x)$ remains positive for all values of
$x$, then $\psi$ decomposes into two components
$e^{iS_{\pm}(x)/\varepsilon}a_{\pm}(x;\varepsilon)$ that do not interact,
resulting in extremely limited valley coupling.

Now, if $\Phi(0)+E-V(x)<0$ is possible, then such values of $x$ are forbidden
energetically. Points $x_{0}$ such that $\Phi(0)+E-V(x_{0})=0$, assuming they
are isolated (i.e., $V^{\prime}(x_{0})\not=0$), are (caustic) turning points.
Any signal arriving at $x_{0}$ from the area where $\Phi(0)+E-V(x)>0$ will
fully turn around so that the amplitudes $a_{\pm}$ are now fully coupled. See,
e.g., [21] for a detailed analysis of coupling (scattering) coefficients in
various semiclassical settings.

We thus obtain a situation of extremely limited coupling if the eikonals
$S_{\pm}$ remain real valued, and full valley-coupling at turning points if
the eikonals $S_{\pm}$ may become complex-valued. We do not consider such a
case further here.

###  3.2 Valley Conductivity

Now that we have constructed valley Hamiltonians $H_{2}$ and $H_{4}$, we
define valley contributions to conductivity so that we can study the
separation of valleys under varying forms of perturbation. Recall we
originally considered the conductivity object
$\sigma_{I}=\operatorname{Tr}{i[H,P]\varphi^{\prime}(H)}$ with $P(x)$ a
regularized Heaviside function. Here we consider $H$ either $H_{2}$ or
$H_{4}$. However, for the two-valley system $\sigma_{I}$ will always be $0$ as
when the two valleys are combined, we have a topologically trivial system.

We thus introduce a conductivity corresponding to a single valley so we can
see the separation in valley contributions to conductivity. $H_{2}$ and
$H_{4}$ are naturally structured to separate the valleys, making the
definition of single valley conductivities reasonably straightforward. To
motivate our definition of valley conductivity, consider a wavefunction
$(\psi_{+},\psi_{-})^{t}$ and define valley switch function

$P_{+}=\begin{pmatrix}P&0\\\ 0&0\end{pmatrix}.$ (3.16)

We consider the current on one side of the interface generated by $P$ via the
expression
$(\psi_{+},\dot{P}_{+}\psi_{+})=\operatorname{Tr}{\dot{P}_{+}\psi_{+}\otimes\psi_{+}^{*}}$.
If $U(t)=e^{-iHt}$ is the evolution of the Schrödinger equation, then

$\dot{P}_{+}:=\partial_{t}(U(t)P_{+}U^{*}(t))=U(t)[H,P_{+}]U^{*}(t).$ (3.17)

We wish to rewrite including all wavefunctions corresponding to the edge. We
pick $\varphi\in\mathfrak{S}(0,1;-E_{0},E_{0})$ with $[-E_{0},E_{0}]$ strictly
in the bulk gap to capture spectral interface information. In particular,
$\varphi^{\prime}(H)$ captures all relevant modes. Decomposing the Hamiltonian
$H$ as $(H_{ij})$ for $1\leq i,j\leq 2$ valley indices, we obtain

$[H,P_{+}]=\begin{pmatrix}[H_{11},P]&-PH_{12}\\\ H_{21}P&0\end{pmatrix}$

so that, with the same decomposition of $\varphi^{\prime}(H)$,

$[H,P_{+}]\varphi^{\prime}(H)=\begin{pmatrix}[H_{11},P]\varphi^{\prime}_{11}(H)-PH_{12}\varphi^{\prime}_{21}(H)&*\\\
*&H_{21}P\varphi^{\prime}_{12}(H)\end{pmatrix}.$

The matrix trace therefore gives

$\sigma_{I,+}(H)=\operatorname{Tr}i[H,P_{+}]\varphi^{\prime}(H)=\operatorname{Tr}i[H_{11},P]\varphi^{\prime}_{11}(H)+\operatorname{Tr}i[H_{21},P]\varphi^{\prime}_{12}(H).$
(3.18)

The conductivity for the opposite valley is

$\sigma_{I,-}(H):=\operatorname{Tr}i[H,P_{-}]\varphi^{\prime}(H)=\operatorname{Tr}i[H_{22},P]\varphi^{\prime}_{22}(H)+\operatorname{Tr}i[H_{12},P]\varphi^{\prime}_{21}(H).$

For the above equation, $H$ is first-order with an off-diagonal term, which is
a multiplication operator in the case of $H_{2}$ and $H_{4}$. This yields
$[H_{21},P]=0$. Hence

$\displaystyle\sigma_{I,+}(H)=\operatorname{Tr}i[H_{11},P]\varphi_{11}^{\prime}(H),\qquad\sigma_{I,-}(H)=\operatorname{Tr}i[H_{22},P]\varphi_{22}^{\prime}(H).$
(3.19)

###  3.3 Numerical example

First, we illustrate the deviation of the valley conductivity (3.18) from the
quantized value $\pm 1$ in the presence of a compactly supported potential $V$
which couples the two valleys, see (3.1). Consider the leading ansatz (3.6)
with

$\displaystyle\psi_{0}(x)=e^{ix/\varepsilon}\psi_{+}(x)+e^{-ix/\varepsilon}\psi_{-}(x)$
(3.20)

with the two Dirac cones at frequencies $\pm\varepsilon^{-1}$. Plugging the
ansatz into (3.3), dividing by $\varepsilon$ and multiplying by $e^{\pm
ix/\varepsilon}$, we obtain, consistently with (3.8) after neglecting highly
oscillatory terms $e^{\pm 2ix/\varepsilon}D_{x}\psi_{\pm}$ and $e^{\pm
2ix/\varepsilon}E$ that formally average to $0$ according to (3.8), the
following one-dimensional coupled Hamiltonian:

$H^{\varepsilon}_{2}=\begin{pmatrix}-D_{x}&0\\\
0&D_{x}\end{pmatrix}+V(x,x/\varepsilon)\begin{pmatrix}1&e^{-2ix/\varepsilon}\\\
e^{2ix/\varepsilon}&1\end{pmatrix},$ (3.21)

so that to leading order $H^{\varepsilon}_{2}\psi=E\psi$ for
$\psi=(\psi_{+}(x),\psi_{-}(x))^{T}$. We choose a two-scale potential of the
form:

$V(x,X)=V_{0}\chi(x)\cos(\omega X),$ (3.22)

where $\chi(x)$ is a smooth, compactly supported function. Model (3.21) allows
us to inspect the effects of the frequency of the fast component of the
potential at a finite value of $\varepsilon$, such that some coupling is
achieved at frequencies different from $\pm 2$, and homogenizes to (3.10) in
the limit $\varepsilon\to 0$, when $\Phi^{\prime}(\pm 1)=\pm 1$.

Here, we perform computations on the interval $[-50,50]$ with periodic
boundary conditions comparing all three models: (3.21) with $\varepsilon=1$
and potential (3.22), and the models (3.10) (one-dimensional) and (3.12) (two-
dimensional) obtained in the limit $\varepsilon\to 0$ with $\omega=2$ with
respective effective potentials described by

$\widehat{V}(x,0)=0,\quad\widehat{V}(x,\pm
2)=\frac{V_{0}}{2}\chi(x)\quad\text{and}\quad\widehat{V}(x,0,y)=0,\quad\widehat{V}(x,\pm
2,y)=\frac{V_{0}}{2}\chi(x)\xi(y).$ (3.23)

We choose functions $\chi(x)$ and $\xi(y)$ supported in the interval
$[-20,20]$ and $[-10,10]$ and taking constant value $1$ in the interval
$[-10,10]$ and $[-5,5]$, respectively. Details of the pseudo-spectral
discretization approach are postponed to section 5. The resulting conductivity
$\sigma^{+}$ is plotted on Figure 6 as a function of $x_{0}$, the center of
the spatial switch function (3.16) and $\omega$, the spatial frequency of the
defect (3.22). In particular, we observe that at low frequencies
($\omega\approx 0$) the conductivity stays very close to the quantized value
$\sigma^{+}=1$, while at $\omega=2$, a maximum coupling between the valleys is
achieved, locally reducing the conductivity measure by a large factor. Note
that as $\varepsilon\to 0$, the limit model (3.10) shows coupling only at
frequency $\omega=2$. On the other hand, away from the defect the perturbation
of the conductivity decays extremely fast regardless of the frequency
$\omega$.

As a second experiment, we compare the valley conductivity as a function of
perturbation amplitude $V_{0}$ and energy $E$ for Hamiltonians (3.10) and
(3.12), that is by numerically computing the conductivity (3.18) when shifting
the Hamiltonian as $H\to H-E$. We observe that the effect of the intervalley
coupling term is to open a band gap at the Dirac cones $E=0$, reducing the
conductivity in a small range of energies around zero for both the one-
dimensional, two-valley Hamiltonian (3.10) and the two-dimensional, two-valley
Hamiltonian (3.12), showing excellent qualitative agreement between the two
models.

Figure 6: Valley conductivity $\sigma_{I,+}$ using the 1D models (3.21),
(3.10) (left) and 2D model (3.12) (right), as a function of: spatial switch
function center $x_{0}$ and frequency $\omega$ of the perturbation with
amplitude $V_{0}=0.1$, using model (3.21) (top left); spatial switch function
center $x_{0}$, with a perturbation with amplitude $V_{0}=0.2$, using model
(3.10) (top right); as a function of amplitude $V_{0}$ and energy $E$ with
frequency $\omega=2$ and at $x_{0}=0$ for the 1D model (3.10) (bottom left)
and 2D model (3.12) (bottom right).

###  3.4 One dimensional scattering theory

The conductivity $\sigma_{I,+}$ in Fig.6 (left) is close to $1$ for $\omega$
away from $2$ as expected from minimal valley coupling and drops significantly
for $\omega=2$, which is consistent with strong valley-coupling. For
intermediate values of $\omega$, we also observe a conductivity larger than
$1$. We now briefly consider a simplified one-dimensional model that exhibits
such a behavior.

We revisit the one dimensional scattering problem (3.3), which we rewrite for
$\Phi(\xi)=\frac{1}{2}(\xi^{2}-1)$ at $\varepsilon=1$ and for a slightly
different expression for the potential as the Helmholtz equation

$\displaystyle(D_{x}^{2}+V(x)-E)\psi=0.$ (3.24)

Assume $E=k^{2}>0$ so that when $V=0$, we obtain $e^{\pm ikx}$ as standing
(propagating) modes associated with each valley. We recall that
$-D_{x}^{2}=\partial^{2}_{x}$. In the presence of $V\not=0$, these modes are
coupled as described by (3.10) in the semiclassical regime. We wish to derive
an explicit expression for valley conductivities using a scattering matrix
formalism in the setting where $P^{\prime}=\delta_{x_{0}}$ for
$x_{0}\in\mathbb{R}$ and where $\varphi^{\prime}=\delta_{E}$ so that we
estimate conductivities at fixed $x_{0}$ and energy $E$.

The two valley conductivities estimate the signals moving from left to right
proportional to $e^{ikx}$ and from right to left proportional to $e^{-ikx}$,
respectively. Let us derive explicit expressions for those in terms of
coefficients of scattering matrices. We wish to compute
$\sigma_{I}^{\pm}[x_{0},E]$ for $\pm$ corresponding to left moving and right
moving signals, respectively.

Consider (3.24) with $E=k^{2}$ and $V$ compactly supported on $[x_{L},x_{R}]$.
We may then construct a standard scattering matrix such that
$\psi(x<x_{L})=e^{ikx}+R_{+}e^{-ikx}$ and $\psi(x>x_{R})=T_{+}e^{ikx}$ for an
incoming plane wave from the left of the support of $V$ and such that
$\psi(x<x_{L})=T_{-}e^{-ikx}$ and $\psi(x>x_{R})=e^{-ikx}+R_{-}e^{ikx}$ for an
incoming plane wave from the right of the support of $V$.

Let us now decompose $V$ into two subcomponents $V_{L,R}$ supported on
$[x_{L},x_{0}]$ and $[x_{0},x_{R}]$, respectively when $x_{L}<x_{0}<x_{R}$.
For each subcomponent, define the scattering matrices

$S^{L}=\begin{pmatrix}R_{+}^{L}&T_{-}^{L}\\\
T_{+}^{L}&R_{-}^{L}\end{pmatrix},\qquad\mbox{ and }\qquad
S^{R}=\begin{pmatrix}R_{+}^{R}&T_{-}^{R}\\\ T_{+}^{R}&R_{-}^{R}\end{pmatrix},$

where the coefficients $R$ and $T$ are interpreted as above.

These two matrices could be appropriately combined to construct a scattering
matrix on the interval $[x_{L},x_{R}]$ called $S$; we do not need the details
of this object.

Both $S^{L}$ and $S^{R}$ as well as $S$ for the whole interval are unitary
scattering matrix so that $SS^{*}=S^{*}S=I$. In particular, we observe that
$|T_{+}|=|T_{-}|$ and $|R_{+}|=|R_{-}|$.

Consider the (unique) solution $\psi_{L}(x)$ of (3.24) at energy $E$ equal to
$e^{ikx}+Re^{-ikx}$ at $x_{L}$ and let $(u^{L},v^{L})$ be defined such that
$\psi_{L}(x_{0})=u^{L}e^{ikx_{0}}+v^{L}e^{-ikx_{0}}$. For this to hold, we
formally assume that $V=0$ at $x_{0}$. Using the above matrices, we find

$u^{L}=T^{L}_{+}+R^{L}_{-}v^{L},\quad v^{L}=R^{R}_{+}u^{L},\qquad
u^{L}=\frac{T_{+}^{L}}{1-R^{L}_{-}R^{R}_{+}},\quad
v^{L}=\frac{R^{R}_{+}T_{+}^{L}}{1-R^{L}_{-}R^{R}_{+}}.$

Considering the solution equal to $e^{-ikx}+Re^{ikx}$ at $x_{R}$ instead, we
find for $(u^{R},v^{R})$ the solution at $x_{0}$

$v^{R}=T_{-}^{R}+R_{+}^{R}u^{R},\quad u^{R}=R^{L}_{-}v^{R},\qquad
v^{R}=\frac{T_{-}^{R}}{1-R_{-}^{L}R_{+}^{R}},\quad
u^{R}=\frac{R_{-}^{L}T_{-}^{R}}{1-R_{-}^{L}R_{+}^{R}}.$

The above corresponding solutions $\psi_{L}(x)$ and $\psi_{R}(x)$ of (3.24)
are orthogonal in the following sense. Assume an inner product
$(u,v)_{M}=(2M)^{-1}\int_{-M}^{M}\bar{u}vdx$ on the interval $(-M,M)$. The
contributions of that integral on the support of $V$ vanishes in the limit
$M\to\infty$. The two solutions constructed above are of the form

$e^{ikx}+R_{+}e^{-ikx}\ \mbox{ and }\ T_{+}e^{ikx}$

for the first solution on the left and right of the support of $V$ and

$T_{-}e^{-ikx}\ \mbox{ and }\ e^{-ikx}+R_{-}e^{ikx}$

for the second solution in the same sense. The inner product of these two
functions for $k\not=0$ is asymptotically given by
$R_{+}^{*}T_{-}+T_{+}^{*}R_{-}=0$. In other words, these two solutions
correspond to two simple branches of absolutely continuous spectrum
(parametrized by $E$) of the operator that live in orthogonal subsets of the
Hilbert space.

The conductivity associated to the signal moving from left to right, which we
associate with one of the valleys is given by the signal at $x_{0}$ moving to
the right for the mixture of states having energy $E$. Thanks to the above
orthogonality, the trace (sum over densities of orthogonal modes evaluated at
$x_{0}$) is given by

$\sigma_{I}^{+}[x_{0},E]:=|u^{L}|^{2}+|u^{R}|^{2}=\frac{|T_{+}^{L}|^{2}+|R^{L}_{-}|^{2}|T_{-}^{R}|^{2}}{|1-R^{L}_{-}R^{R}_{+}|^{2}}.$

The conductivity associated to the signal moving from the right to the left is
given by

$\sigma_{I}^{-}[x_{0},E]:=-(|v^{L}|^{2}+|v^{R}|^{2})=-\frac{|T_{-}^{R}|^{2}+|R^{R}_{+}|^{2}|T_{+}^{L}|^{2}}{|1-R^{L}_{-}R^{R}_{+}|^{2}}.$

The sum of the nominators is given by

$|T_{+}^{L}|^{2}+(1-|T^{L}_{-}|^{2})|T_{-}^{R}|^{2}-|T_{-}^{R}|^{2}-(1-|T^{R}_{+}|^{2})|T_{+}^{L}|^{2}=0.$

So, we do verify that $\sigma_{I}^{+}+\sigma_{I}^{-}=0$ as it should be for a
globally topologically trivial material.

When $x_{0}$ is outside the support of $V$, say $x_{0}\geq x_{R}$, then
$T^{R}=1$ and $R^{R}=0$ so that

$\sigma_{I}^{+}=|T_{+}|^{2}+|R_{-}|^{2}=1.$

Similarly, $\sigma_{I}^{-}$ would be equal to $-1$ outside of the support of
$V$. It is therefore only on the support of the fluctuation $V$ that we expect
to observe spatial variations (in $x_{0}$) of the valley conductivities.

Let us consider a simplified setting. Assume that
$x_{0}=\frac{1}{2}(x_{L}+x_{R})$ and that $V$ is even about $x_{0}$. Then by
symmetry, $R_{\pm}^{L}=R^{R}_{\mp}$ and $T_{\pm}^{L}=T_{\mp}^{R}$ so that

$\sigma_{I}^{+}=\dfrac{|T^{L}_{+}|^{2}(1+|R_{-}^{L}|^{2})}{|1-(R_{-}^{L})^{2}|^{2}}=\dfrac{(1-|R_{-}^{L}|^{2})(1+|R_{-}^{L}|^{2})}{|1-(R_{-}^{L})^{2}|^{2}}=\dfrac{1-|R_{-}^{L}|^{4}}{|1-(R_{-}^{L})^{2}|^{2}}.$

Denoting $R^{L}_{-}=e^{i\theta}r$, we find

$\sigma_{I}^{+}=\dfrac{1-r^{4}}{|1-e^{2i\theta}r^{2}|^{2}}\approx 1+2\cos
2\theta\ r^{2},$

for $r$ small (corresponding to $V$ small). The term $\cos 2\theta$ is
positive or negative depending on the value of $\theta$ we obtain from the
scattering theory. When $\theta=\frac{\pi}{2}$ so that $R_{-}^{L}$ is purely
imaginary, then we observe the largest drop in conductivity. When $\theta=0$
however, we observe an increase in the local conductivity for small $V$ before
it decreases of large values of $V$. These observations are consistent with
what we obtained in numerical simulations in Fig.6 and show that in the
presence of intervalley-coupling, the conductivities associated to left-going
and right-going modes are spatially dependent and hence not protected
topologically.

As a concrete example, assume $V=V_{0}\chi_{|[-a,a]}$ as a barrier potential
for the Left domain above leading to the computation of $R_{-}^{L}$. With
$q=\sqrt{k^{2}-V_{0}}$ assumed positive, standard calculations in the
scattering theory of the one-dimensional Schrödinger equation then yield

$R^{L}_{-}=\frac{V_{0}\sin(2qa)e^{-2ika}}{(2k^{2}-V_{0})\sin(2qa)+2ikq\cos(2qa)}.$

When the width $a$ is small, then to leading order in $a$, the phase of
$R^{L}_{-}$ is close to $e^{-i\frac{\pi}{2}}$. It may however, take arbitrary
values for larger width $a$. Note that $R^{L}_{-}=0$ when $2qa=n\pi$ for
$n\in\mathbb{N}$.

##  4 Junction topology

This section addresses asymmetric transport through a junction in the
$xy-$plane. We introduce a junction conductivity and compute it by means of a
Fedosov-Hörmander formula, which extends what exists for the flat interface
case [5, 29]. Pseudo-differential calculus will be used throughout; see
Appendix B for notation and background.

Although our main result will apply to a more general class of models and
geometries, we use Region 2 from Figure 1 as a motivating example. As shown,
the region selects one vertex from the underlying periodic hexagonal
structure. We push all other vertices to infinity (recall that the side length
of each triangle in Figure 1 is inversely proportional to the twist angle,
with the latter assumed to be small), so that our model consists of six slices
alternating between AB and BA stacking.

The described geometry supports surface waves that propagate along the
interfaces between the AB and BA regions. These waves have a preferred
direction, e.g. moving from left to right. (For numerical examples, see
Section 5.) The role of the junction conductivity is to quantify the net
current associated with these waves at a particular energy. One can define the
conductivity so that it measures the contribution from any combination of
interfaces (more details below Corollary 4.5; see also Figure 8).

Our main result (Theorem 4.1 below) is a bulk-interface correspondence, which
equates the conductivity to the difference of two integrals depending only on
the associated insulators. In particular, the conductivity does not depend on
local properties of the material, such as the way the insulators are joined
together at each interface. Theorem 4.1 allows for straight-forward
evaluations of the conductivity using previously developed theory for
interface models, e.g. [8, 29].

The proof of Theorem 4.1 requires familiar stability properties of the
junction conductivity, resembling those of [6, 5, 29, 30]. We tie the
conductivity to a Fredholm index (Theorem 4.4), with the latter known to be
quantized and invariant with respect to a large class of perturbations. This
implies a conservation law (Corollary 4.5) showing (among other things) that
the currents entering and leaving the junction are the same. We also prove
invariance of the junction conductivity with respect to the corresponding
density of states (Theorem 4.6), semiclassical rescaling (Theorem 4.7), and
relatively bounded (Theorem 4.8) and compact (Theorem 4.9) perturbations of
the Hamiltonian. We conclude the section with an application to a $2\times 2$
Dirac model, which is then analyzed numerically in Section 5.

###  4.1 Junction conductivity for tBLG

Let us now introduce more mathematical detail. The junction Hamiltonian
corresponding to region 2 from Figure 1 is

$\displaystyle\operatorname{Op}(\sigma)=H=\begin{pmatrix}\Omega
I+D\cdot\sigma^{(\eta)}&\lambda U^{*}(x,y)\\\ \lambda U(x,y)&-\Omega
I+D\cdot\sigma^{(\eta)}\end{pmatrix},$ (4.1)

where the notation on the left-hand side is defined in Appendix B, and
$U(x,y):=\frac{1}{2}((1+\tilde{m}(x,y))A+(1-\tilde{m}(x,y))A^{*})$. Here,
$\tilde{m}\in\mathcal{C}^{\infty}$ is positive in regions of AB stacking,
negative in regions of BA stacking, and zero on the boundaries between
regions. We assume that $\tilde{m}\in\\{-1,1\\}$ away from these boundaries.
Recall (2.3), where the analogous Hamiltonian for an interface (rather than
junction) model is defined.

To construct such a function $\tilde{m}$ explicitly, define
$f_{\Theta}:\mathbb{T}\rightarrow\mathbb{T}$ by
$f_{\Theta}(\theta)=\sin(3\theta)$ and let
$f\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ such that
$f(x,y)=rf_{\Theta}(\theta)$ for all $r\geq 1$, where $(r,\theta)$ are the
polar coordinates corresponding to $(x,y)$. For concreteness, take
$f(x,y):=\chi(r)rf_{\Theta}(\theta)$, where
$\chi\in\mathfrak{S}(0,1;\varepsilon,1)$ and $0<\varepsilon<1$. (Recall the
definition of $\mathfrak{S}$ below (2.4).) An example of such a function $f$
is plotted in Figure 7 (left panel). Since $\chi$ is smooth and vanishes near
the origin, it is immediate that $f$ is smooth. Let $m\in\mathfrak{S}(-1,1)$
be monotonically increasing with $m(0)=0$, and define
$\tilde{m}(x,y):=m(f(x,y))$. Here, $\mathfrak{S}(a_{1},a_{2})$ is the union of
$\mathfrak{S}(a_{1},a_{2};c_{1},c_{2})$ over all $c_{1}\leq c_{2}$. This
definition of $\tilde{m}$ is consistent with the above constraints, as
$f_{\Theta}>0$ if and only if $2k\pi/3<\theta<(2k+1)\pi/3$ for some
$k\in\\{0,1,2\\}$. Since $m$ and all of its derivatives are bounded, it
follows that $\sigma\in S^{1}_{1,0}$ with the right-hand side defined in
Appendix B. Moreover, the structure of $H$ (with $D\cdot\sigma^{(\eta)}$ on
the diagonal and no other derivatives) implies the existence of a constant
$c>0$ such that $|\sigma_{\min}(x,y,\xi,\zeta)|\geq
c\langle\xi,\zeta\rangle-1$, where $\sigma_{\min}$ is the smallest magnitude
eigenvalue of $\sigma$, and for any vector $u$, we define $\langle
u\rangle=\sqrt{1+|u|^{2}}$. We see that when $f(x,y)$ is sufficiently large
(resp. small), $\sigma=\sigma_{+}$ (resp. $\sigma=\sigma_{-}$), with
$\operatorname{Op}(\sigma_{\pm}):=H_{\pm}$ and $H_{\pm}$ defined in section 2.
By Appendix A, there exists $E>0$ such that $\sigma_{\pm}$ (equivalently
$H_{\pm}$) has a spectral gap in the interval $(E_{1},E_{2}):=(-E,E)$.

Analogous to (2.5) and [3, 6, 5, 29], we define the _junction conductivity_ by

$\displaystyle\sigma_{I}(H,P):=\operatorname{Tr}i[H,P]\varphi^{\prime}(H),$
(4.2)

where $\varphi\in\mathfrak{S}(0,1;-E,E)$. The operator (of point-wise
multiplication by) $P$ is still a regularized indicator function, but it can
no longer depend on only one variable. Indeed, for the operator on the right-
hand side of (4.2) to be trace-class, it is necessary that $\nabla P$ decay to
$0$ along the level curves of $f$ (otherwise, the symbol of
$[H,P]\varphi^{\prime}(H)$ would not decay in $\langle x,y\rangle$). This
condition is satisfied if we set $P(x,y)=\chi_{p}(g(x,y))$, where
$\chi_{p}\in\mathfrak{S}(0,1)$ and $g\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$
such that $C_{1}\langle x,y\rangle\leq\langle f,g\rangle\leq C_{2}\langle
x,y\rangle$ for some positive constants $C_{1}$ and $C_{2}$. To construct such
a function $g$ explicitly, we can for example define
$g_{\Theta}:\mathbb{T}\rightarrow\mathbb{T}$ by
$g_{\Theta}(\theta)=\cos(\theta)+\cos(\pi/6)=\cos(\theta)+\sqrt{3}/2$ (note
that $g_{\Theta}<0$ if and only if $5\pi/6<\theta<7\pi/6$) and take
$g(x,y)=\chi(r)rg_{\Theta}(\theta)$. See Figure 7 (center panel) for a plot of
$g$. Since the zeros of $f_{\Theta}$ and $g_{\Theta}$ are disjoint, it is
clear that the growth condition for $\langle f,g\rangle$ is satisfied. As we
will show in detail below (see Proposition B.1), the symbol of
$\varphi^{\prime}(H)$ (resp. $[H,P]$) decays rapidly in $\langle
f(x,y),\xi,\zeta\rangle$ (resp. $\langle g(x,y)\rangle$), meaning that
$\sigma_{I}$ in (4.2) is well defined.

Figure 7: Plots of $f$ (left) and $g$ (center) in the $xy-$plane, with several
level curves included. The right panel shows $f$ evaluated along the level
curves $\\{g=a\\}$ for three values of $a$. The curves are parametrized by the
variable $t$, with the convention that $y_{a}(t)\rightarrow\pm\infty$ as
$t\rightarrow\pm\infty$.

Due to the factor of $\chi$ in the definition of $g$, the level set
$g^{-1}(0)$ is not a curve in the $xy-$plane. Still, there exists a smooth
function $(x_{0},y_{0}):\mathbb{R}\rightarrow\mathbb{R}^{2}$ whose range
separates the $xy-$plane into two regions, $R_{+}$ and $R_{-}$, such that
$g\geq 0$ in $R_{+}$ and $g\leq 0$ in $R_{-}$. More specifically, we can take

$\displaystyle(x_{0}(t),y_{0}(t))=\begin{cases}(\frac{\sqrt{3}}{2}t,\frac{1}{2}t),&t\leq-1\\\
(-\frac{\sqrt{3}}{2}t,\frac{1}{2}t),&t\geq 1\end{cases}$

and smoothly connect $(x_{0},y_{0})$ for $|t|<1$ so that
$g(x_{0}(t),y_{0}(t))=0$ for all $t\in\mathbb{R}$. For
$(t,\xi,\zeta)\in\mathbb{R}^{3}$, define

$\displaystyle\tau(t,\xi,\zeta):=\sigma(x_{0}(t),y_{0}(t),\xi,\zeta),\qquad\tau_{z}:=z-\tau.$
(4.3)

By assumption, there exist $t_{0}>0$ and symbols

$\tau_{\pm}(\xi,\zeta)\in\Bigg{\\{}\begin{pmatrix}\Omega
I+\xi\sigma_{1}+\eta\zeta\sigma_{2}&\lambda A^{*}\\\ \lambda A&-\Omega
I+\xi\sigma_{1}+\eta\zeta\sigma_{2}\end{pmatrix},\begin{pmatrix}\Omega
I+\xi\sigma_{1}+\eta\zeta\sigma_{2}&\lambda A\\\ \lambda A^{*}&-\Omega
I+\xi\sigma_{1}+\eta\zeta\sigma_{2}\end{pmatrix}\Bigg{\\}}$

such that $\tau(t,\xi,\zeta)=\tau_{\pm}(\xi,\zeta)$ whenever $\pm t\geq
t_{0}$. The main result of this section is the following bulk-interface
correspondence.

######  Theorem 4.1.

Let $\alpha\in(E_{1},E_{2})$ and define $\tau_{\pm,z}:=z-\tau_{\pm}$ with
$z=\alpha+i\omega$. Then

$\displaystyle 2\pi\sigma_{I}=\frac{i}{8\pi^{2}}(I_{+}-I_{-}),\qquad
I_{\pm}=\int_{\mathbb{R}^{3}}\operatorname{tr}[\tau_{\pm,z}^{-1}\partial_{\xi}\tau_{\pm,z},\tau_{\pm,z}^{-1}\partial_{\zeta}\tau_{\pm,z}]\tau_{\pm,z}^{-1}d\omega
d\xi d\zeta.$ (4.4)

We postpone the proof to Appendix B. Although motivated and introduced by the
tBLG model from Section 2, the above result applies to any Hamiltonian $H$ and
corresponding energy interval $(E_{1},E_{2})$ satisfying (H1) below; see
Section 4.2. Moreover, the function $g$ (and hence the regularized indicator
function $P$) does not need to be defined as in the above example.

An immediate consequence of Theorem 4.1 is

######  Corollary 4.2.

Let $\tilde{H}=\operatorname{Op}(\tau)$ with $\tau=\tau(y,\xi,\zeta)$ and
$\tilde{P}(x)=\tilde{P}\in\mathfrak{S}(0,1)$. Then
$\sigma_{I}(H,P)=\sigma_{I}(\tilde{H},\tilde{P})$.

###### Proof.

Recall the definition (4.3) of $\tau$, where $\sigma$ is the symbol of $H$
(that is, $\operatorname{Op}(\sigma)=H$). It follows that the bulk symbol
$\tau_{+}$ is the same for $H$ and $\tilde{H}$, and thus so is the quantity
$I_{+}$ appearing in (4.4). The same holds for $\tau_{-}$ and $I_{-}$. We
apply Theorem 4.1 and the proof is complete. ∎

Observe that Corollary 4.2 reduces the junction conductivity to the simpler
setting of a flat interface separating two insulators. The insulating
materials are described by $\tau_{\pm}=\tau(\pm\infty,\xi,\zeta)$. For the
tBLG Hamiltonian (4.1) and spatial filter $P$ in the above example, we can now
apply Theorem 2.1 and use the fact that

$m(f(x_{0}(y),y_{0}(y)))=:m(y)=m\in\mathfrak{S}(-1,1)$

to obtain

$\displaystyle 2\pi\sigma_{I}(H,P)=-2\eta\,\text{sign}(\Omega).$

If $P_{j}=\chi_{p}(g_{j}(x,y))$ for $j\in\\{0,1,2\\}$ with the $g_{j}$
appropriately chosen to satisfy the growth condition for $\langle f,g\rangle$,
then $\sigma_{I}(H,\sum_{j}P_{j})=\sum_{j}\sigma_{I}(H,P_{j}).$ Thus if
$g_{j}(x,y)=\chi(r)rg_{j_{\Theta}}(\theta)$ with
$g_{j_{\Theta}}=\cos(\theta-2\pi j/3)+\cos(\pi/6)$, then the
$2\pi/3-$rotational symmetry of the tBLG Hamiltonian implies that
$\sigma_{I}(H,\sum_{j}P_{j})=3\sigma_{I}(H,P_{0}).$ That is, each term
$\sigma_{I}(H,P_{j})$ in the sum over $j$ sees the same transition from
$\sigma_{-}$ to $\sigma_{+}$. In words, the superposition of three
appropriately chosen regularized indicator functions increases the
conductivity by a factor of $3$. This makes perfect sense, as the conductivity
is now measured through three distinct regions ($\operatorname{supp}\nabla
P_{j}$ for $j=0,1,2$), each of which contributes the same value of
${{{-}}}\eta\,\text{sign}(\Omega)/\pi$.

###  4.2 General setting

In this subsection, we provide a full class of operators $H$ and $P$ to which
Theorem 4.1 applies. Let us introduce some notation. As above, the
Hamiltonians $H$ act on
$\mathcal{H}:=L^{2}(\mathbb{R}^{2})\otimes\mathbb{C}^{n}$. We label the
spatial coordinates $(x,y)\in\mathbb{R}^{2}$ and the corresponding dual
variables $(\xi,\zeta)\in\mathbb{R}^{2}$. Given functions
$u\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ and
$A\in\mathcal{C}^{\infty}(\mathbb{R}^{2};\mathbb{C}^{n\times n})$, fixed
matrices $A_{1}$ and $A_{2}$ in $\mathbb{C}^{n\times n}$, and fixed constants
$c_{1}\leq c_{2}$, we write $A\in\mathfrak{S}(A_{1},A_{2};c_{1},c_{2};u)$ to
mean that

$\displaystyle A=\begin{cases}A_{1},&u<c_{1}\\\ A_{2},&u>c_{2}\end{cases}.$
(4.5)

We let $\mathfrak{S}(A_{1},A_{2};u)$ denote the union of
$\mathfrak{S}(A_{1},A_{2};c_{1},c_{2};u)$ over all $-\infty<c_{1}\leq
c_{2}<\infty$. If $A=A(\alpha)$ is a function of just one variable, then we
define
$\mathfrak{S}(A_{1},A_{2};c_{1},c_{2}):=\mathfrak{S}(A_{1},A_{2};c_{1},c_{2};\alpha)$
and $\mathfrak{S}(A_{1},A_{2}):=\mathfrak{S}(A_{1},A_{2};\alpha)$. Note that
if $A(x,y)=\chi_{A}(u(x,y))$ for some $\chi_{A}\in\mathfrak{S}(A_{1},A_{2})$,
then $A\in\mathfrak{S}(A_{1},A_{2};u)$.

Let $(r,\theta)$ denote the polar coordinates for $(x,y)$. Fix
$k\in\mathbb{N}$ and let $\Theta_{k}\subset[0,2\pi)$ be any set containing
$2k$ elements. Let $f_{\Theta}\in\mathcal{C}^{\infty}(\mathbb{T})$ such that
$f_{\Theta}(\theta)=0$ if and only if $\theta\in\Theta_{k}$, and
$f^{\prime}_{\Theta}(\theta)\neq 0$ whenever $\theta\in\Theta_{k}$. This means
$f_{\Theta}$ changes sign across every point in $\Theta_{k}$ and is bounded
away from $0$ for all $\theta$ away from $\Theta_{k}$. Let
$f\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ such that
$f(x,y)=f(r,\theta):=rf_{\Theta}(\theta)$ whenever $r\geq 1$. We assume the
following:

(H1) Let $H=\operatorname{Op}(\sigma)$ be a symmetric pseudo-differential
operator, with $\sigma$ independent of $h$. Suppose $\sigma\in S^{m}_{1,0}$
for some $m>0$, and $|\sigma_{\min}(x,y,\xi,\zeta)|\geq
c\langle\xi,\zeta\rangle^{m}-1$ for some $c>0$. Suppose there exist values
$-\infty<E_{1}<E_{2}<\infty$ and symbols $\sigma_{\pm}\in S^{m}_{1,0}$
independent of $(x,y)$ with no spectrum in $(E_{1},E_{2})$, such that
$\sigma(\cdot,\cdot,\xi,\zeta)\in\mathfrak{S}(\sigma_{-}(\xi,\zeta),\sigma_{+}(\xi,\zeta);f)$
for all $(\xi,\zeta)\in\mathbb{R}^{2}$.

+-+-$g=0$

-+-+-+$g=0$

Figure 8: Illustrations of $\sigma$ in the $xy$-plane for $k=2$ (left) and
$k=3$ (right). The labels $\pm$ indicate the sign of $f$. Two example level
curves $\\{g=0\\}$ are presented in each case (solid and dashed curves), for
some arbitrary choices of $g$ satisfying the appropriate assumptions. The
scale of this picture is much smaller than $1$, so that we are not in the
region where $g(r,\theta)=rg_{\Theta}(\theta)$.

Above, $\sigma_{\min}$ denotes the smallest magnitude eigenvalue of $\sigma$.
By [11, 23], we know (H1) implies that
$H:\mathcal{D}(H)\rightarrow\mathcal{H}$ is self-adjoint with
$\mathcal{D}(H)=\mathcal{H}^{m}$ (and the right-hand side the standard Hilbert
space of functions with square-integrable derivatives up to order $m$).
Moreover, for any $z\in\mathbb{C}$ with $\Im z\neq 0$, we have
$(z-H)^{-1}\in\operatorname{Op}(S^{-m}_{1,0})$.

Let $\\{\theta_{1},\theta_{2}\\}\subset\mathbb{T}$ such that
$\theta_{1},\theta_{2}\notin\Theta_{k}$. Let
$g_{\Theta}\in\mathcal{C}^{\infty}(\mathbb{T})$ such that
$g_{\Theta}(\theta)=0$ if and only if $\theta\in\\{\theta_{1},\theta_{2}\\}$,
and $g^{\prime}_{\Theta}(\theta_{j})\neq 0$ for $j\in\\{1,2\\}$. Let
$g\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ such that
$g(x,y)=g(r,\theta)=rg_{\Theta}(\theta)$ for all $r\geq 1$. Our assumptions on
$f$ and $g$ imply the existence of positive constants $C_{1}<C_{2}$ such that

$\displaystyle C_{1}\langle x,y\rangle\leq\langle f(x,y),g(x,y)\rangle\leq
C_{2}\langle x,y\rangle.$ (4.6)

Moreover, there exists a smooth function
$(x_{0},y_{0}):\mathbb{R}\rightarrow\mathbb{R}^{2}$ with
$\lim_{|t|\rightarrow\infty}|(x_{0}(t),y_{0}(t))|=\infty$ and whose range
separates the $xy-$plane into two regions, $R_{+}$ and $R$, such that $g\geq
0$ in $R_{+}$. One can for example take a curve contained in $\\{g=0\\}$, in
which case $(x_{0}(t),y_{0}(t))=(\beta_{1}t,\beta_{2}t)$ for all $|t|$
sufficiently large, with the constants $\beta_{1}$ and $\beta_{2}$ depending
only on $\operatorname{\text{sign}}t$. Our convention is that $R_{+}$ be on
the _right_ side of the curve defined by $(x_{0},y_{0})$ as $t$ increases.
That is, if you were to move along the curve with $t$ increasing, then $R_{+}$
would be on your right. Thus $\tau$ in Theorem 4.1 can always be defined by
(4.3) in this more general setting.

We conclude this subsection with the following result, which verifies that our
theory applies to the tBLG Hamiltonian.

######  Proposition 4.3.

Define $H=\operatorname{Op}(\sigma)$ by (4.1). Then $H$ satisfies (H1).

We refer to Appendix B for the proof.

###  4.3 Stability of the junction conductivity

In this subsection, we show that the junction conductivity $\sigma_{I}(H,P)$
is invariant with respect to a large class of perturbations of the Hamiltonian
$H$, spatial filter $P$ and density of states $\varphi^{\prime}$. While
interesting in their own right, these results are also used in the proof of
Theorem 4.1 (see Appendix B). Throughout, we will assume that
$H=\operatorname{Op}(\sigma)$ satisfies (H1), and let
$\varphi\in\mathfrak{S}(0,1;E_{1},E_{2})$ and $P\in\mathfrak{S}(0,1;g(x,y))$
be smooth switch functions. Let $\bar{P}\in\bar{\mathfrak{S}}(0,1;g(x,y))$ be
a projector ($\bar{P}^{2}=\bar{P}$). Here, $\bar{\mathfrak{S}}(A_{1},A_{2};u)$
is the space of functions $A$ satisfying (4.5) for some $c_{1}\leq c_{2}$.
(Note that $A\in\bar{\mathfrak{S}}(A_{1},A_{2};u)$ need not be continuous.)
Let $\Phi\in\mathcal{C}^{\infty}_{c}(E_{1},E_{2})$ such that
$\varphi^{\prime}\in\mathcal{C}^{\infty}_{c}(\\{\Phi=1\\}^{\circ})$, where
$U^{\circ}$ denotes the interior of $U$. We will tie the junction conductivity
to a Fredholm index, which will be used to prove that the former is stable.

######  Theorem 4.4.

Let $\tilde{H}$ be a self-adjoint pseudo-differential operator with
$\tilde{H}\in\operatorname{Op}(S(\langle\xi,\zeta\rangle^{m}))$ and
$\Phi(\tilde{H})\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))$. Then
$[\tilde{H},P]\varphi^{\prime}(\tilde{H})$ is trace-class so that
$\sigma_{I}(\tilde{H},P)$ is well-defined. Moreover, defining
$U(\tilde{H}):=e^{i2\pi\varphi(\tilde{H})}$, we have that
$\bar{P}U(\tilde{H})\bar{P}$ is a Fredholm operator on the range of $\bar{P}$,
with

$\displaystyle 2\pi\sigma_{I}(\tilde{H},P)={\rm
Index}(\bar{P}U(\tilde{H})\bar{P}).$

Note that (H1) is a strict subset of the class of operators $\tilde{H}$ to
which the above result applies. (That any $H$ satisfying (H1) must also
satisfy the assumptions of Theorem 4.4 will be proved rigorously in Appendix
B; see Proposition B.1.) For a proof of Theorem 4.4, see Appendix B.

Theorem 4.4 is a powerful result with many implications. For example, we
immediately have

######  Corollary 4.5.

If $\tilde{H}$ satisfies the assumptions of Theorem 4.4 and
$P_{1},P_{2}\in\mathfrak{S}(0,1;g(x,y))$, then
$\sigma_{I}(\tilde{H},P_{1})=\sigma_{I}(\tilde{H},P_{2})$.

###### Proof.

Applying Theorem 4.4, we have $2\pi\sigma_{I}(\tilde{H},P_{1})={\rm
Index}(\bar{P}U(\tilde{H})\bar{P})=2\pi\sigma_{I}(\tilde{H},P_{2})$. ∎

For an intuitive explanation of the above corollary, recall that our model is
of two types of materials ($+$ and $-$) that are smoothly glued together at a
junction. See Figure 8 for two examples, where we assume for concreteness that
$P$ transitions from $0$ to $1$ in the vicinity of $\\{g=0\\}$. The solid
curve (in both panels) can see only one transition ($\sigma_{+}$ to
$\sigma_{-}$ and $\sigma_{-}$ to $\sigma_{+}$ in the left and right panels,
respectively), while the dashed level curves contain multiple transitions
($+\rightarrow-\rightarrow+\rightarrow-$ for the left panel and
$-\rightarrow+\rightarrow-\rightarrow+\rightarrow-\rightarrow+$ for the right
panel). We refer also to Figure 7 (right panel) for an illustration of $f$
evaluated along various level curves of $g$ when $k=3$. As expected, the
conductivity only cares about the “starting” and “ending” topology along the
level curve and is unaffected by oscillations in between. Observe that
Corollary 4.5 describes a conservation law, as the solid and dashed curves
respectively measure the conductivity entering and leaving the junction.

We will use Theorem 4.4 to prove invariance of $\sigma_{I}(\tilde{H},P)$ with
respect to $\tilde{H}$. Since a Fredholm index must be integer-valued, we know
that any continuous perturbation of $\tilde{H}$ would leave the junction
conductivity unchanged; see Theorems 4.7-4.9 below. First, we show that the
conductivity is independent of the density of states we choose.

######  Theorem 4.6.

Take $\tilde{H}$ as in Theorem 4.4 and let $\varphi_{1}\in\mathfrak{S}(0,1)$
such that
$\varphi^{\prime}_{1}\in\mathcal{C}^{\infty}_{c}(\\{\Phi=1\\}^{\circ})$. Then

$\displaystyle\sigma_{I}(\tilde{H},P)=\operatorname{Tr}i[\tilde{H},P]\varphi^{\prime}_{1}(\tilde{H}).$

We now show that the conductivity is independent of semiclassical rescaling.
To this end, define $H_{h}:=\operatorname{Op}_{h}(\sigma)$ for $0<h\leq 1$,
where $H=\operatorname{Op}(\sigma)$ is assumed to satisfy (H1). This allows
for an expansion of $\sigma_{I}(H_{h},P)$ in the semiclassical parameter $h$,
which will be used to prove Theorem 4.2 (see Appendix B). We refer to [6, 5,
29], where a similar semiclassical analysis is employed.

######  Theorem 4.7.

For any $h\in(0,1]$ we have $\sigma_{I}(H_{h},P)=\sigma_{I}(H,P)$.

Let $W$ be a symmetric pseudo-differential operator and define
$H^{(\mu)}:=H+\mu W$ for $\mu\in[0,1]$.

######  Theorem 4.8.

If $W\in\operatorname{Op}(S^{m}_{1,0})$, then
$\sigma_{I}(H^{(\mu)},P)=\sigma_{I}(H,P)$ for all $\mu>0$ sufficiently small.

Theorem 4.8 asserts the stability of the junction conductivity with respect to
relatively bounded perturbations that are small enough. To avoid this
smallness condition, we must make the stronger assumption that the
perturbation be relatively compact; that is, we require the symbol of $W(i\pm
H)^{-1}$ to decay at infinity in all variables.

######  Theorem 4.9.

If $W\in\operatorname{Op}(S^{m}_{1,0}\cap
S(\langle\xi,\zeta\rangle^{m-\delta}\langle x,y\rangle^{-\delta}))$ for some
$\delta>0$, then $\sigma_{I}(H^{(\mu)},P)=\sigma_{I}(H,P)$ for all
$\mu\in[0,1]$.

We refer to Appendix B for the proofs of Theorems 4.6-4.9. Recall that $m=1$
for the tBLG Hamiltonian (4.1).

###  4.4 Dirac model

We conclude this section with an application of Theorem 4.1 to the $2\times 2$
Dirac Hamiltonian given by

$H=D_{x}\sigma_{1}+D_{y}\sigma_{2}+\tilde{m}(x,y)\sigma_{3},$ (4.7)

where $\tilde{m}(x,y)=m(f(x,y))$ for some $m\in\mathfrak{S}(-1,1)$, and

$\displaystyle
f(r,\theta)=\chi(r)rf_{\Theta}(\theta),\qquad\chi\in\mathfrak{S}(0,1;\varepsilon,1),\qquad
f_{\Theta}(\theta)=\sin(k\theta)$

for some $0<\varepsilon<1$ and $k\in\mathbb{N}_{+}$. Note that when $k=1$, we
recover the setting of a flat interface analyzed in detail in [3, 4, 6, 5,
29], while $k=3$ yields the same hexagonal structure analyzed at the beginning
of this section. As with Proposition 4.3, it is straightforward to verify that
$H$ defined by (4.7) satisfies (H1), and thus the above theory applies. The
bulk spectral gap in this case is $(E_{1},E_{2})=(-1,1)$.

Take $P(x,y)=\chi_{p}(g(x,y))$ for $\chi_{p}\in\mathfrak{S}(0,1)$, where
$g(x,y)=g(r,\theta)=\chi(r)rg_{\Theta}(\theta)$ and
$g_{\Theta}(\theta)=\cos(\theta-\theta_{1})-\cos(\theta_{0})$, for some
$0<\theta_{0}<\pi$ and $-\frac{\pi}{k}<\theta_{1}<\frac{\pi}{k}$ satisfying

$\theta_{+}\notin\Theta_{k},\quad
2\pi-\theta_{+}\notin\Theta_{k},\quad\Theta_{k}=\Big{\\{}\frac{j}{k\pi}:j\in\mathbb{Z}\Big{\\}},\quad\theta_{+}:=\theta_{0}+\theta_{1}.$

This way, the zeros of $g_{\Theta}$ and $f_{\Theta}$ are disjoint. Indeed,
$g_{\Theta}(\theta)=0$ if and only if
$\cos(\theta-\theta_{1})=\cos(\theta_{0})$, which occurs exactly when
$\\{\theta+2\pi
j:j\in\mathbb{Z}\\}\cap\\{\theta_{+},2\pi-\theta_{+}\\}\neq\emptyset$. It
follows that
$\tau(t,\xi,\zeta)=\xi\sigma_{1}+\zeta\sigma_{2}+\mu(t)\sigma_{3}$ for some
$\mu\in\cup_{\varepsilon_{1},\varepsilon_{2}\in\\{-1,1\\}}\mathfrak{S}(\varepsilon_{1},\varepsilon_{2})$,
where the $\varepsilon_{j}$ are determined by $\theta_{0}$ and $\theta_{1}$.
We have reduced the problem to computing the interface conductivity for the
translation-invariant $2\times 2$ Dirac system [3, 5, 29], and hence
$2\pi\sigma_{I}=\frac{1}{2}(\varepsilon_{1}-\varepsilon_{2})$. For the case
$(k,\theta_{1})=(3,0)$ which is analyzed numerically in the following section,
we have

$\displaystyle 2\pi\sigma_{I}=\begin{cases}-1,&0<\theta_{0}<\pi/3\\\
1,&\pi/3<\theta_{0}<2\pi/3\\\ -1,&2\pi/3<\theta_{0}<\pi\end{cases}.$

If $g$ from Figure 8 were to have the above form, then $\theta_{1}=0$ and
$\theta_{0}$ would be the angle between [the line making up the top of the
level curve $\\{g=0\\}$] and [the positive $x-$axis]. So
$\pi/2<\theta_{0}<\pi$ and $2\pi/3<\theta_{0}<\pi$ for the left and right
panels, respectively.

##  5 Numerical Simulations

We now present results of numerical simulations illustrating the stability of
the interface and junction conductivities analyzed in the preceding sections.
We augment these spectral results with simulations of wavepackets propagating
across the junction considered in section 4 (case $k=3$ in Fig. 8).

###  5.1 Methods and Discretization

All computations are conducted on a finite interval of size $L$ for one-
dimensional problems such as (3.10) (see section 3.3) or a rectangular domain
of size $L_{x}\times L_{y}$ for two-dimensional problems such as
conductivities for the Hamiltonian (3.12) or (4.7), equipped with periodic
boundary conditions. We use a pseudo-spectral discretization where all
functions are approximated using truncated Fourier series:

$u(x)=\sum_{k=-K}^{K}\widehat{u}_{k}e^{2i\pi kx/L}\quad\text{or}\quad
u(x,y)=\sum_{k=-K_{x}}^{K_{x}}\sum_{l=-K_{y}}^{K_{y}}\widehat{u}_{k,l}e^{2i\pi
kx/L_{x}+2i\pi ly/L_{y}},$

such that derivation operators are represented as diagonal matrices. Pointwise
multiplication operators such as $Vu(x)=v(x)u(x)$ are computed using Fourier
interpolation on a uniform real-space grid of size $3(K+1)$ in the one-
dimensional case and $3(K_{x}+1)\times 3(K_{y}+1)$ in the two-dimensional one
employing the discrete Fourier transform $\mathcal{F}$, avoiding the aliasing
of products for small values of $K$ and sharp or highly oscillatory mass or
potential profiles $m(x,y),v(x,y)$.

Computations of conductivities using finite domains are accomplished using the
approximation

$2\pi\widetilde{\sigma}(H)=2\pi\operatorname{Tr}iQ[H,P]\phi^{\prime}(H),$
(5.1)

where $P=p(x,y)$ and $Q=q(x,y)$ are pointwise multiplication operators where,
as before, $p(x,y)$ is a spatial smooth switch function across a certain
contour, and $q(x,y)$ is a spatial filter designed to mask the boundaries
which may host artificial, unwanted domain walls (necessary to have space-
dependent coefficients such as $m(x,y)$ remain smoothly connected across the
periodic boundaries of the torus). More precisely, we set here:

$\displaystyle p_{x_{0},\delta}(x,y)$
$\displaystyle=\begin{cases}1,&\max\left\\{x-x_{0},\quad-\cos(\theta)(y-L_{y}/2)\pm\sin(\theta)(x-x_{0}-L_{x}/4)\right\\}\geq\delta,\\\
0,&\max\left\\{x-x_{0},\quad-\cos(\theta)(y-L_{y}/2)\pm\sin(\theta)(x-x_{0}-L_{x}/4)\right\\}\leq-\delta,\end{cases}$
(5.2) and $\displaystyle q_{\delta}(x,y)$
$\displaystyle=\begin{cases}1,&\max\left\\{|x-L_{x}/2|-L_{x}/4,\quad|y-L_{y}/2|-3L_{y}/8\right\\}\leq-\delta,\\\
0,&\max\left\\{|x-L_{x}/2|-L_{x}/4,\quad|y-L_{y}/2|-3L_{y}/8\right\\}\geq\delta,\end{cases}$

such that $\delta$ is a small parameter controlling the width of transition
regions, $\theta=\pi-\frac{\pi}{12}$ is as in Section 4 and $x_{0}$ is the
abscissa at which the support of $\nabla p_{x_{0},\delta}$ crosses the
horizontal axis: see Figure 9 for such a profile. Valley-projected numerical
conductivities $\widetilde{\sigma}_{\pm}$ are defined similarly (3.18).

Figure 9: Top left: mass profile $m(x,y)$ yielding a three-fold symmetric
topological junction. The contour $m(x,y)=0$ is highlighted in dashed black
line. Top right: spatial switch function $p_{x_{0},\delta}(x,y)$ for
$x_{0}=L_{x}/2$ used in the computation of the conductivity (5.1). Bottom row:
magnitude profiles of the gradient of the spatial switch function $\nabla
p_{x_{0},\delta}(x,y)$ multiplied by the filter $q_{0}(x,y)$ (5.2) for the
three values (from left to right) $x_{0}=3L_{x}/8,\ L_{x}/2,\ 5L_{x}/8$
presented in Table 1.

Practically, the evaluation of (5.1) is achieved by direct diagonalization of
the discretized Hamiltonian $H$, resulting in the sum-over-states formula:

$2\pi\widetilde{\sigma}(H)=2\pi\sum_{j}\phi^{\prime}(\lambda_{j})\left\langle
u_{j}\left|Q_{\delta}\begin{pmatrix}0&\partial_{x}p_{x_{0},\delta}-i\partial_{y}p_{x_{0},\delta}\\\
\partial_{x}p_{x_{0},\delta}+i\partial_{y}p_{x_{0},\delta}&0\end{pmatrix}\right|u_{j}\right\rangle,$

where $\lambda_{j}$, $u_{j}$ are the eigenvalues and eigenvectors of $H$ and
the scalar product is evaluated by quadrature on the uniform real-spaced grid
discussed earlier.

The spectral accuracy of the numerical conductivity (5.1) using our spectral
discretization scheme can be established rigorously by observing that the
relevant edge states vanish exponentially fast in the bulk. We do not pursue
this analysis here and refer to [29] for the case of a single straight domain
wall for elliptic partial differential operators, including the massive Dirac
operator.

###  5.2 Invariant Calculations for Junctions

As a first numerical experiment illustrating the results of Section 4, we
compute the conductivity in a hexagonal junction setup as in region 2 from
Figure 1, i.e. a vertex with three-fold symmetry with the topology as in
Figure 8. Specifically, our choice of spatial projector (5.2) with
$\theta=\pi-\frac{\pi}{12}$ corresponds to a function $g$ whose level set
$g^{-1}(0)$ fully encloses only the left horizontal branch of the boundary set
$m^{-1}(0)$, i.e. the negative real axis.

| $N=8$ | $N=16$ | $N=32$ | $N=64$
---|---|---|---|---
$x_{0}=3L_{x}/8$ | $0.36144$ | $0.92310$ | $0.99923$ | $0.99983$
$x_{0}=L_{x}/2$ | $0.53574$ | $0.72033$ | $0.99612$ | $0.99997$
$x_{0}=5L_{x}/8$ | $0.05648$ | $0.88647$ | $0.99901$ | $0.99993$
Table 1: Conductivity values (5.1) computed for different values of the
spatial switch function center $x_{0}$ along the $x$-axis, and discretization
parameter $N=K_{x}=K_{y}$.

Results from Section 4 predict that the quantity $2\pi\sigma_{I}$ takes the
value $1$ independently of the exact choice of such a spatial projector, in
particular whether $x_{0}<L_{x}/2$ and the level set $g^{-1}(0)$ intersects
the negative real axis only, or $x_{0}>L_{x}/2$ and it crosses the other five
boundaries. Our numerical computations, presented in Table 1, are consistent
with the theory, with the converged results ($N=64$) agreeing up to four
digits of accuracy with the predicted value.

Table 1 also illustrates the stability and efficiency of our pseudo-spectral
discretization approach. For all three choices of parameter $x_{0}$, we note
the exponential convergence of the computed conductivity with the frequency
cutoff parameter $N$, consistent with the robustness of the conductivity
demonstrated in earlier sections. Numerical experiments (not presented) also
indicate fast exponential convergence in the size of the domain $L_{x}\times
L_{y}$, consistent with the exponential localization of eigenmodes in the
energy window of interest around boundaries of interest as well as artificial
domain walls induced by periodization.

###  5.3 Wavepacket Propagation through Junctions

Figure 10: Scattering of a wavepacket incoming from the left along a
horizontal branch towards a hexagonal cross, with a rotationally symmetric
mass profile.

Next, we study numerically the time-domain propagation of wavepackets through
a hexagonally symmetric junction given by the same mass profile as plotted in
Figure 9. Let us note that the stable topological invariants presented in
Section 4, which can be interpreted as counting the net number of modes or
currents traveling in certain directions along the interfaces [27, 4, 3], only
partly predict the dynamics of localized wavepackets traveling across a
junction. In particular, the sign of the conductivity identifies certain
branches of the interface $m^{-1}(0)$ as incoming channels on which net
propagation of current is towards the junction, while on others current is
outgoing, directed outwards from the junction (with the sum of current on all
channels always being zero). On the other hand, the steady-state current
generated by the state $\phi^{\prime}(H)$ is equally distributed between all
channels, and does not allow to distinguish how the current ’originating’ from
one of the incoming channels distributes itself.

Semiclassical analysis, in the case of non-intersecting, straight or curved
interfaces, shows that for Dirac operators such as the one considered here
(4.7), these currents have a local interpretation via the construction of
long-lived localized wavepackets travelling unidirectionally [7], with some
dispersion due to the curvature of the interface. Hence, wavepackets
travelling along an interface towards the junction are expected to scatter
into any number wavepackets travelling on the outgoing channels.

While such rigorous analysis has not been carried out in the case of
junctions, our numerical simulations of such scattering across a hexagonal
junction show that a Gaussian wavepacket, initially travelling along the
incoming channel supported by the horizontal edge to the left of the junction,
splits three-fold at the intersection along the three outgoing branches of the
junction: see snapshots of the evolution presented in Figure 10. More
precisely, an initial Gaussian packet is prepared at $t=0$ on the left
horizontal branch which forms an incoming channel towards the intersection,
and we compute its time evolution until $t=50$ where the packet has completed
its travel through the intersection. In contrast to the numerical results
presented in [7], we observe remarkably little dispersion of the packet
despite the sharp turns in the direction of propagation.

As mentioned above, the topology measured by the quantized junction
conductivity, while heuristically linked to the robustness and stability of
the propagating edge modes, does not distinguish between channels and thus
does not constrain the entries of the scattering matrix between modes
supported by incoming and outgoing channels. We observe on Figure 10 that the
wavepacket does not split equally between all three branches, with a smaller
amount travelling straight across the intersection.

Figure 11: Scattering of a wavepacket through a hexagonal cross with broken
symmetry profiles, allowing to steer a wavepacket incoming from the left
towards any outgoing branch.

To illustrate the absence of stability of the scattering with respect to
perturbations, we present in Figure 11 the final state of simulations with the
same initial wavepacket as in Figure 10, but where a small, local perturbation
of the mass term breaks the rotational symmetry: namely we set

$m(x,y)=m_{0}(x,y)-\frac{1}{4}\exp\left(-\frac{(x-L_{x}/2)^{2}+(y-L_{y}/2)^{2}}{2\sigma^{2}}\right)\begin{pmatrix}\sin(\theta_{m})\\\
\cos(\theta_{m})\end{pmatrix}\cdot\begin{pmatrix}x-L_{x}/2\\\
y-L_{y}/2\end{pmatrix},$

where $m_{0}(x,y)$ is the symmetric mass profile as in Figure 9, $\sigma=5$ is
a width parameter and $\theta_{m}=2\pi/3,\ 0,\ -2\pi/3$ respectively from left
to right. These simulations indicate that the scattering can be completely
manipulated by such local mass perturbations, and is not a protected property
of the system: indeed, the wavepacket can be steered towards any of the three
outgoing branches.

##  6 Conclusion

We consider continuum partial (and pseudo) differential models of twisted
bilayer graphene (tBLG) that accurately describe macroscopic transport
properties in tBLG. The twisted layers generate a macroscopic triangular moiré
pattern separating regions of AB and BA stacking as shown in Fig. 1.

The primary focus of this work is to analyze the asymmetric transport observed
at the interface between insulating AB and BA stackings (region 1 in Fig. 1)
and at the junctions of the triangular pattern (region 2 in Fig. 1).

Asymmetry is quantified by observables of the system referred to as
conductivities. A familiar line conductivity analyzed for such pseudo-
differential equations in [3, 6, 5, 29] describes the asymmetry in region 1.
We adapted the derivations in [8] in the context of Floquet topological
insulators to compute the line conductivity and show that it took integral
values by means of a bulk-interface correspondence.

The aforementioned asymmetry takes opposite values in each valley of the
problem so that asymmetric transport is observed only when inter-valley
coupling is negligible. We presented a class of simplified Hamiltonians
modeling such valley coupling. We defined partial line conductivities for the
two valleys and showed formally and numerically that these conductivities were
no longer quantized and hence no longer stable with respect to perturbations.

Finally, to constrain transport across junctions in region 2 of Fig. 1, we
introduced a notion of junction conductivity defined for a large class of
domain walls modeling the insulating regions near the junction and for a large
class of (pseudo-)differential operators including the Dirac operators
modeling tBLG. Adapting techniques developed in [29], we demonstrated that the
junction conductivity was indeed quantized and obtained by a Fedosov-Hörmander
formula generalizing that for line conductivities. Numerical simulations of
junction conductivities and other transport phenomena for a simple model of
Dirac equation confirmed the theoretical predictions.

## Acknowledgments

This research was partially supported by the National Science Foundation,
Grants DMS-1908736, DMS-1819220 and EFMA-1641100.

## References

  * [1] M. Anđelković, L. Covaci, and F. M. Peeters.  DC conductivity of twisted bilayer graphene: Angle-dependent transport properties and effects of disorder.  Phys. Rev. Materials, 2:034004, Mar 2018.
  * [2] J. E. Avron, R. Seiler, and B. Simon.  Charge deficiency, charge transport and comparison of dimensions.  Communications in mathematical physics, 159:399–422, 1994.
  * [3] G. Bal.  Continuous bulk and interface description of topological insulators.  J. Math. Phys., 60(8):081506, 2019.
  * [4] G. Bal.  Topological protection of perturbed edge states.  Comm. Mathematical Sciences, 17(1):193–225, 2019.
  * [5] G. Bal.  Topological invariants for interface modes.  Comm. Partial Differential Equations, 47:1636–1679, 2022.
  * [6] G. Bal.  Topological charge conservation for continuous insulators.  To appear in J. Math. Phys. and arXiv:2106.08480, 2023.
  * [7] G. Bal, S. Becker, A. Drouot, C. F. Kammerer, J. Lu, and A. Watson.  Edge state dynamics along curved interfaces.  To appear in SIAM J. Math. Anal. and arXiv:2106.00729, 2021.
  * [8] G. Bal and D. Massatt.  Multiscale invariants of Floquet topological insulators.  Multiscale Modeling & Simulation, 20(1):493–523, 2022.
  * [9] B. A. Bernevig and T. L. Hughes.  Topological Insulators and Topological Superconductors.  Princeton University Press, 2013.
  * [10] R. Bistritzer and A. H. MacDonald.  Moiré bands in twisted double-layer graphene.  Proceedings of the National Academy of Sciences - PNAS, 108(30):12233–12237, 2011.
  * [11] J.-M. Bony.  On the Characterization of Pseudodifferential Operators (Old and New), pages 21–34.  Springer New York, New York, NY, 2013.
  * [12] E. Cancès, L. Garrigue, and D. Gontier.  Simple derivation of moiré-scale continuous models for twisted bilayer graphene.  Phys. Rev. B, 107:155403, Apr 2023.
  * [13] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero.  Unconventional superconductivity in magic-angle graphene superlattices.  Nature, 556(7699):43–50, 2018.
  * [14] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras.  Relaxation and domain formation in incommensurate two-dimensional heterostructures.  Phys. Rev. B, 98:224102, Dec 2018.
  * [15] S. Dai, Y. Xiang, and D. J. Srolovitz.  Twisted bilayer graphene: Moiré with a twist.  Nano Letters, 16(9):5923–5927, 09 2016.
  * [16] M. Dimassi and J. Sjostrand.  Spectral Asymptotics in the Semi-Classical Limit.  London Mathematical Society Lecture Note Series. Cambridge University Press, 1999.
  * [17] A. Drouot.  The bulk-edge correspondence for continuous honeycomb lattices.  Communication in Partial Differential Equations, 44(12):1406–1430, 2019.
  * [18] V. Duchêne and M. I. Weinstein.  Scattering, homogenization, and interface effects for oscillatory potentials with strong singularities.  Multiscale Modeling & Simulation, 9(3):1017–1063, 2011.
  * [19] P. Elbau and G. Graf.  Equality of bulk and edge hall conductance revisited.  Communications in mathematical physics, 229(3):415–432, 2002.
  * [20] M. Fruchart and D. Carpentier.  An introduction to topological insulators.  Comptes Rendus Physique, 14(9):779 – 815, 2013.  Topological insulators / Isolants topologiques.
  * [21] C. Gérard and A. Grigis.  Precise estimates of tunneling and eigenvalues near a potential barrier.  Journal of differential equations, 72(1):149–177, 1988.
  * [22] L. Hörmander.  The weyl calculus of pseudo-differential operators.  Communications on Pure and Applied Mathematics, 32(3):359–443, 1979\.
  * [23] L. Hörmander.  The analysis of linear partial differential operators III: Pseudo-differential operators.  Springer Science & Business Media, 2007.
  * [24] N. Kalton.  Trace-class operators and commutators.  Journal of Functional Analysis, 86(1):41–74, 1989.
  * [25] N. Lerner.  Metrics on the phase space and non-selfadjoint pseudo-differential operators, volume 3.  Springer Science & Business Media, 2011.
  * [26] E. McCann and M. Koshino.  The electronic properties of bilayer graphene.  Reports on Progress in Physics, 76(5):056503, apr 2013.
  * [27] E. Prodan and H. Schulz-Baldes.  Bulk and boundary invariants for complex topological insulators.  Springer, Berlin, 2016.
  * [28] Z. Qiao, J. Jung, Q. Niu, and A. H. MacDonald.  Electronic highways in bilayer graphene.  Nano Letters, 11(8):3453–3459, 08 2011.
  * [29] S. Quinn and G. Bal.  Approximations of interface topological invariants.  arXiv preprint arXiv:2112.02686, 2021.
  * [30] S. Quinn and G. Bal.  Asymmetric transport for magnetic dirac equations.  arXiv preprint arXiv:2211.00726, 2022.
  * [31] P. Rickhaus, J. Wallbank, S. Slizovskiy, R. Pisoni, H. Overweg, Y. Lee, M. Eich, M.-H. Liu, K. Watanabe, T. Taniguchi, T. Ihn, and K. Ensslin.  Transport through a network of topological channels in twisted bilayer graphene.  Nano Letters, 18(11):6725–6730, 11 2018.
  * [32] P. San-Jose and E. Prada.  Helical networks in twisted bilayer graphene under interlayer bias.  Phys. Rev. B, 88:121408, Sep 2013.
  * [33] H. Tang, S. Carr, and E. Kaxiras.  Geometric origins of topological insulation in twisted layered semiconductors.  Physical Review B, 104(15), Oct 2021.
  * [34] G. Usaj, P. M. Perez-Piskunow, L. E. F. Foa Torres, and C. A. Balseiro.  Irradiated graphene as a tunable floquet topological insulator.  Phys. Rev. B, 90:115423, Sep 2014.
  * [35] A. B. Watson, T. Kong, A. H. MacDonald, and M. Luskin.  Bistritzer–MacDonald dynamics in twisted bilayer graphene.  Journal of Mathematical Physics, 64(3):031502, 2023.
  * [36] E. Witten.  Three lectures on topological phases of matter.  La Rivista del Nuovo Cimento, 39(7):313–370, 2016.
  * [37] F. Zhang, A. H. MacDonald, and E. J. Mele.  Valley chern numbers and boundary modes in gapped bilayer graphene.  Proceedings of the National Academy of Sciences, 110(26):10546–10551, 2013.
  * [38] J. Zhou, S. Cheng, W.-L. You, and H. Jiang.  Effects of intervalley scattering on the transport properties in one-dimensional valleytronic devices.  Scientific Reports, 6(1):23211, 2016.
  * [39] M. Zworski.  Semiclassical analysis, volume 138.  American Mathematical Soc., 2012.

##  Appendix A Proof of Theorem 2.1

For the bilayer problem with coupling
$A=\varepsilon\frac{1}{2}(\sigma_{1}+i\sigma_{2})$, we recall that the tBLG
Hamiltonian is

$H_{+}(\xi)=\begin{pmatrix}\Omega&\bar{\xi}&0&0\\\ \xi&\Omega&\varepsilon&0\\\
0&\varepsilon&-\Omega&\bar{\xi}\\\ 0&0&\xi&-\Omega\end{pmatrix}.$

Consider first the case $\Omega>0$ and $\varepsilon>0$. Note that for $\xi=0$,
$\pm\Omega$ are eigenvalues associated to eigenvectors $e_{1}$ and $e_{4}$
while $\pm\sqrt{\Omega^{2}+\varepsilon^{2}}$ are the last eigenvalues with
eigenspaces supported on the span of $(e_{2},e_{3})$.

The constraint ${\rm det}(H-E)=0$ on the eigenvalues is given by

$(|\xi|^{2}-(E-\Omega)^{2})(|\xi|^{2}-(E+\Omega)^{2})=\varepsilon^{2}(E^{2}-\Omega^{2}).$

The four eigenvalues (called $E_{j}$ for $1\leq j\leq 4$ in the main text)
come in pairs with opposite signs, are all different, and are given explicitly
via

$E_{\pm}^{2}=(\Omega^{2}+\frac{1}{2}\varepsilon^{2}+|\xi|^{2})\pm\sqrt{(4\Omega^{2}+\varepsilon^{2})|\xi|^{2}+\frac{1}{4}\varepsilon^{4}}.$

To reflect this symmetry, we denote the eigenvalues as $-E_{+}=E_{1}$,
$-E_{-}=E_{2}$, $E_{-}=E_{3}$ and $E_{+}=E_{4}$. We observe that $E^{2}_{-}$
is large for $|\xi|$ large while $E=0$ implies
$(|\xi|^{2}-\Omega^{2})^{2}+\varepsilon^{2}\Omega^{2}=0$ which is impossible,
and hence the presence of a spectral gap. More precisely, we obtain that
$|\xi|^{2}\to E_{-}^{2}$ is a convex function with a minimum equal to

$E^{2}_{\rm
min}=\frac{\Omega^{2}\varepsilon^{2}}{4\Omega^{2}+\varepsilon^{2}}\quad\mbox{and
attained at
}\quad|\xi|^{2}=\frac{2\Omega^{2}(2\Omega^{2}+\varepsilon^{2})}{4\Omega^{2}+\varepsilon^{2}}.$

The equations for the eigenvectors are

$(\Omega-E)u_{1}+\bar{\xi}u_{2}=0\quad\xi u_{1}+(\Omega-E)u_{2}+\varepsilon
u_{3}=0$ $\varepsilon u_{2}-(\Omega+E)u_{3}+\bar{\xi}u_{4}=0,\quad\xi
u_{3}-(\Omega+E)u_{4}=0.$

We wish to construct eigenvectors that are smooth in the variable $\xi$. For
$E=E_{+}$ the largest eigenvalue, we can choose the smooth eigenvector

$u(\xi)=c_{u}\Big{(}\frac{\bar{\xi}}{E-\Omega}(\Omega+E-\frac{|\xi|^{2}}{\Omega+E}),(\Omega+E-\frac{|\xi|^{2}}{\Omega+E}),\varepsilon,\frac{\varepsilon\xi}{\Omega+E}\Big{)}$

with $c_{u}=c_{u}(\xi)$ chosen such that $|u|=1$. Since $E_{+}\pm\Omega>0$ is
bounded below by a positive constant, $\mathbb{C}\ni\xi\to c_{u}(\xi)$ is a
smooth function bounded above and below by positive constants.

This vector converges as $|\xi|\to\infty$ to an expression that depends on
$(\varepsilon,\Omega)$ in a non-trivial manner. Namely, in the limit
$|\xi|\to\infty$:

$0<E_{\pm}^{2}=|\xi|^{2}\pm\sqrt{4\Omega^{2}+\varepsilon^{2}}|\xi|+O(1),$

and $u(\xi)$ converges in the same sense to

$u_{+}^{A}(\hat{\xi})=c_{\infty}(\hat{\bar{\xi}}\beta,\beta,\varepsilon,\varepsilon\hat{\xi}),\quad\beta=2\Omega+\sqrt{4\Omega^{2}+\varepsilon^{2}}\approx\Omega+E-\frac{|\xi|^{2}}{\Omega+E},\quad
c_{\infty}=(2\varepsilon^{2}+2\beta^{2})^{-\frac{1}{2}}.$

We obtain similarly that the differential 1-form $du$ converges to
$du_{+}^{A}$ as $|\xi|\to\infty$. As a result, we have

$W_{+}^{4}=\frac{i}{2\pi}\lim_{R\to\infty}\displaystyle\int_{\mathbb{S}^{1}_{R}}(u,du)=\frac{i}{2\pi}c_{\infty}^{2}\displaystyle\int_{\mathbb{S}^{1}}\big{(}(\beta
e^{-i\theta},\beta,\varepsilon,\varepsilon e^{i\theta}),(-i\beta
e^{-i\theta},0,0,i\varepsilon e^{i\theta})\big{)}d\theta.$

This implies that

$W_{+}^{4}=\frac{i}{2\pi}c_{\infty}^{2}(-i\beta^{2}+i\varepsilon^{2})2\pi=\dfrac{\beta^{2}-\varepsilon^{2}}{2(\varepsilon^{2}+\beta^{2})}.$

Let us now look at the second positive eigenvalue $E=E_{-}>0$. The (smooth)
function $\xi\to E_{-}(\xi)-\Omega$ vanishes at $\xi=0$ (and in fact at two
other values of $\xi$). Thus we may no longer divide by $(\Omega-E)$ in the
construction of a smooth vector $u(\xi)$. However, it is easy to verify that
$\xi^{-1}(E(\xi)-\Omega)$ is a smooth function. Eliminating $u_{2}$ and
$u_{4}$ from the above system for $u$, we obtain

$u=c_{u}\Big{(}\varepsilon,\frac{E-\Omega}{\bar{\xi}}\varepsilon,-(\xi-\frac{(\Omega-E)^{2}}{\bar{\xi}}),\frac{-\xi}{\Omega+E}(\xi-\frac{(\Omega-E)^{2}}{\bar{\xi}})\Big{)}.$

This converges as $|\xi|\to\infty$ to

$u_{-}^{A}=c_{\infty}(\varepsilon,\varepsilon\hat{\xi},-\beta\hat{\xi},-\beta\hat{\xi}^{2}).$

Following the same procedure as above, we deduce that

$W_{+}^{3}=\frac{-\varepsilon^{2}-3\beta^{2}}{2(\varepsilon^{2}+\beta^{2})},\qquad\qquad
W_{+}=W_{+}^{3}+W_{+}^{4}=-1.$

We now consider the second Hamiltonian and replace $A$ by $A^{*}$:

$H_{-}(\xi)=\begin{pmatrix}\Omega&\bar{\xi}&0&\varepsilon\\\ \xi&\Omega&0&0\\\
0&0&-\Omega&\bar{\xi}\\\ \varepsilon&0&\xi&-\Omega\end{pmatrix}.$

We observe that the equation for the eigenvalues $E$ is the same as in the
previous case. For the largest positive eigenvalue, one finds

$u=c_{u}\Big{(}\Omega+E-\frac{|\xi|^{2}}{\Omega+E},\frac{\xi}{E-\Omega}(\Omega+E-\frac{|\xi|^{2}}{\Omega+E}),\frac{\hat{\xi}}{\Omega+E}\varepsilon,\varepsilon\Big{)}.$

In the limit $|\xi|\to\infty$, we have

$u_{+}^{A^{*}}=c_{\infty}(\beta,\beta\hat{\xi},\hat{\bar{\xi}}\varepsilon,\varepsilon),\quad\mbox{so
that}\quad
W_{-}^{4}=\dfrac{\varepsilon^{2}-\beta^{2}}{2(\varepsilon^{2}+\beta^{2})}.$

For the eigenvalue $0<E=E_{-}$, we find

$u=c_{u}\Big{(}\frac{E-\Omega}{|\xi|}\varepsilon,\varepsilon\frac{\bar{\xi}}{\Omega+E}\big{(}\frac{(\Omega-E)^{2}}{\xi}-\bar{\xi}\big{)},\frac{(\Omega-E)^{2}}{\xi}-\bar{\xi}\Big{)}.$

In the limit $|\xi|\to\infty$,

$u_{-}^{A^{*}}=c_{\infty}\big{(}\hat{\bar{\xi}}\varepsilon,\varepsilon,-\beta\hat{\bar{\xi}}^{2},-\beta\hat{\bar{\xi}}\big{)}\quad\mbox{so
that}\quad
W_{-}^{3}=\dfrac{\varepsilon^{2}+3\beta^{2}}{2(\varepsilon^{2}+\beta^{2})}\quad\mbox{
and }\quad W_{-}=W_{-}^{3}+W_{-}^{4}=1.$

Figure 12: A family of two projectors
$(\pm,\xi)\to\Pi_{\pm}(\xi):=|\psi_{\pm}(\xi)\rangle\langle\psi_{\pm}(\xi)|$
is mapped to the sphere by mapping each plane $\pm$ onto a hemisphere while
preserving continuity along the ‘circle at infinity’, i.e., the limits
$\lim_{|\xi|\to\infty}\Pi_{+}(|\xi|\hat{\xi})=\lim_{|\xi|\to\infty}\Pi_{-}(|\xi|\hat{\xi})$
exist and agree for all $\hat{\xi}\in\mathbb{S}^{1}$.

Upon inspection of the limiting vectors $u_{\pm}^{A}$ and $u_{\pm}^{A^{*}}$,
we observe that $u_{+}^{A}$ and $u_{+}^{A^{*}}$ are vectors that are not
defined up to a multiplicative phase. The same holds true for $u_{-}^{A}$ and
$u_{-}^{A^{*}}$. However, we observe that

$v_{-}^{A^{*}}=\sigma_{2}\otimes
I_{2}u_{-}^{A^{*}}=i\big{(}\beta\hat{\bar{\xi}}^{2},\beta\hat{\bar{\xi}},\hat{\bar{\xi}}\varepsilon,\varepsilon\big{)}=i\hat{\bar{\xi}}\big{(}\beta\hat{\bar{\xi}},\beta,\varepsilon,\varepsilon\hat{\xi}\big{)}=i\hat{\bar{\xi}}u_{+}^{A}.$

So the following gluing procedure is possible. Let
$\Pi_{+}(\xi)=u_{+}^{4}\otimes u_{+}^{4}(\xi)$ be the projector where
$u_{+}^{4}(\xi)$ converges to $u_{+}^{A}(\hat{\xi})$ as $|\xi|\to\infty$. Let
similarly $\Pi_{-}(\xi)=v_{-}^{3}\otimes v_{-}^{3}(\xi)$ with
$v_{-}^{3}(\xi)=\sigma_{2}\otimes I_{2}u_{-}^{3}$ converging to
$v_{-}^{A^{*}}$ as $|\xi|\to\infty$. Since $v_{-}^{A^{*}}$ is proportional to
$u_{+}^{A}$, then $\Pi_{+}(|\xi|\hat{\xi})$ and $\Pi_{-}(|\xi|\hat{\xi})$ have
the same (well defined) limit as $|\xi|\to\infty$ for each
$\hat{\xi}\in\mathbb{S}_{1}$. We can then apply the gluing procedure of [5] to
obtain that

$-1=W_{+}^{4}-W_{{{-}}}^{3}=\frac{i}{2\pi}\displaystyle\int_{\mathbb{R}^{2}}{{{\operatorname{tr}}}}\Pi_{+}d\Pi_{+}\wedge
d\Pi_{+}-{{{\operatorname{tr}}}}\Pi_{-}d\Pi_{-}\wedge
d\Pi_{{{-}}}\in\mathbb{Z}$

is an integer-valued Chern number integrated over a sphere (the
compactification of two planes $\mathbb{R}^{2}$ over which the two projectors
$\Pi_{\pm}$ are defined); see Figure 12.

We can similarly glue the $E_{+}$ branch of $A^{*}$ with $E_{-}$ branch of $A$
and get another invariant $W_{+}^{3}-W_{-}^{4}=1$ using that

$v_{+}^{A^{*}}=\sigma_{2}\otimes
I_{2}u_{+}^{A^{*}}=-i(\hat{\bar{\xi}}\varepsilon,\varepsilon,-\beta,-\beta\hat{\xi})=-i\hat{\bar{\xi}}u_{-}^{A}.$

In some sense, the second Hamiltonian with coupling term $A^{*}$ should be
replaced by the equivalent operator $\sigma_{2}\otimes
I_{2}\hat{H}_{A^{*}}\sigma_{2}\otimes I_{2}$. No matter the basis, we observe
that in order to get bulk-difference invariants, we need to glue projectors in
the AB and BA sectors corresponding to different eigenvalues.

This concludes the proof of Theorem 2.1 when $\Omega>0$,
$\varepsilon\equiv\lambda>0$ and $\eta=1$. It is reasonably straightforward to
observe that the sign of $\varepsilon$ has no influence on the derivation. All
we need is $\varepsilon\not=0$ so that a gap opens. The sign of $\Omega$,
however, matters since for $\Omega<0$ and $E>0$, it is $(\Omega+E)^{-1}$ that
is no longer bounded.

To address how $W$ changes depending on system parameters $\Omega,\varepsilon$
and $\eta$, define

$H_{\pm}(\xi;\Omega,\varepsilon,\eta)=\Omega\sigma_{3}\otimes
I_{2}+I_{2}\otimes(\xi_{1}\sigma_{1}+\eta\xi_{2}\sigma_{2})+\frac{\varepsilon}{2}(\sigma_{1}\otimes\sigma_{1}\pm\sigma_{2}\otimes\sigma_{2}).$
(A.1)

We then define the following symmetry operations:

$\displaystyle S_{1}=\sigma_{1}\otimes I_{2},\qquad
S_{2}\psi(\xi)=\psi(-\xi),\qquad S_{3}=\sigma_{1}\otimes\sigma_{1}.$

We let $W(\Omega,\varepsilon,\eta)$ denote the bulk-difference invariant as
described above. For $\Omega,\varepsilon,\eta>0$, by the above calculation we
have

$W(\Omega,\varepsilon,\eta)=-2.$

We have the symmetry relations

$\displaystyle
S_{1}H_{\pm}(\xi;\Omega,\varepsilon,\eta)S_{1}=H_{\mp}(\xi;-\Omega,\varepsilon,\eta),$
$\displaystyle
S_{2}H_{\pm}(\xi;\Omega,\varepsilon,\eta)S_{2}=-H_{\pm}(\xi;-\Omega,-\varepsilon,\eta),$
$\displaystyle
S_{3}H_{\pm}(\xi;\Omega,\varepsilon,\eta)S_{3}=H_{\pm}(\xi;-\Omega,\varepsilon,-\eta).$

Since $(S_{j}\psi,d(S_{j}\psi))=(S_{j}\psi,S_{j}d\psi)=(\psi,d\psi)$ for
$j\in\\{1,3\\}$, and likewise $\int(\psi,d\psi)=\int(S\psi,dS_{2}\psi)$, the
connection is invariant to the symmetry operations. However, when the sign of
the Hamiltonian changes, i.e. if $S_{j}HS_{j}=-\tilde{H}$, the invariant of
$H$ is minus the invariant of $\tilde{H}$ as the order of the bands is
reversed, and the sum of the four connnections (including positive and
negative energies) is zero. We thus conclude

$W(\Omega,\varepsilon,\eta)=-W(-\Omega,\varepsilon,\eta)=W(\Omega,-\varepsilon,\eta)=W(-\Omega,\varepsilon,-\eta)=-W(\Omega,\varepsilon,-\eta).$

This concludes the proof of (2.8) in the theorem.

For (2.9), we can apply the bulk-edge correspondence in [29] since $H_{e}$
satisfies the required ellipticity conditions. That is,
$H_{e}=\operatorname{Op}(\sigma)$, where $\sigma\in S^{1}_{1,0}$ with all of
its singular values bounded below by $\langle\xi,\zeta\rangle$ whenever
$\langle y,\xi,\zeta\rangle$ is sufficiently large. Moreover, the bulk
Hamiltonians $H_{\pm}$ have a spectral gap at $0$, as demonstrated by the
above calculations. Thus $H_{e}$ satisfies (H1) from [29], meaning that (2.9)
holds.

##  Appendix B Proofs from Section 4 and relevant notation

###  B.1 Notation and functional setting

We first briefly define the notation used in Section 4 regarding pseudo-
differential operators. For a more detailed exposition, see [29, section 2]
and references therein.

Given a parameter $h\in(0,1]$ and a symbol
$a(x,\xi;h)\in\mathcal{S}^{\prime}(\mathbb{R}^{d}\times\mathbb{R}^{d})\otimes\mathbb{M}_{n}$,
we define the Weyl quantization of $a$ as the operator

$\displaystyle\operatorname{Op}_{h}(a)\psi(x):=\frac{1}{(2\pi
h)^{d}}\int_{\mathbb{R}^{2d}}e^{i(x-y)\cdot\xi/h}a(\frac{x+y}{2},\xi;h)\psi(y)dyd\xi,\qquad\psi\in\mathcal{S}(\mathbb{R}^{d})\otimes\mathbb{C}^{n}.$
(B.1)

Here, $\mathcal{S}$ denotes the Schwartz space (with $\mathcal{S}^{\prime}$
its dual) and $\mathbb{M}_{n}$ the space of Hermitian $n\times n$ matrices. We
define $\operatorname{Op}(a):=\operatorname{Op}_{1}(a)$.

A function $u:\mathbb{R}^{2d}\rightarrow[0,\infty)$ is called an order
function if there exist constants $C_{0}>0$, $N_{0}>0$ such that $u(X)\leq
C_{0}\langle X-Y\rangle^{N_{0}}u(Y)$ for all $X,Y\in\mathbb{R}^{2d}$. Here we
use the notation $\langle X\rangle:=\sqrt{1+|X|^{2}}$. Note that if $u_{1}$
and $u_{2}$ are order functions, then so is $u_{1}u_{2}$.

We say that $a\in S(u)$ if for every $\alpha\in\mathbb{N}^{2d}$, there exists
$C_{\alpha}>0$ such that $|\partial^{\alpha}a(X;h)|\leq C_{\alpha}u(X)$ for
all $X\in\mathbb{R}^{2d}$ and $h\in(0,1]$. We write $S(u^{-\infty})$ to denote
the intersection over $s\in\mathbb{N}$ of $S(u^{-s})$. For $\delta\in[0,1]$
and $k\in\mathbb{R}$, we say that $a(X;h)\in S^{k}_{\delta}(u)$ if for every
$\alpha\in\mathbb{N}^{2d}$, there exists $C_{\alpha}>0$ such that

$\displaystyle|\partial^{\alpha}a(X;h)|\leq
C_{\alpha}u(X)h^{-\delta|\alpha|-k},$ (B.2)

uniformly in $X\in\mathbb{R}^{2d}$ and $h\in(0,1]$. If either $k$ or $\delta$
are omitted, they are assumed to be zero.

By [16, Chapter 7], we know that if $a\in S(u_{1})$ and $b\in S(u_{2})$, then
$\operatorname{Op}_{h}(c):=\operatorname{Op}_{h}(a)\operatorname{Op}_{h}(b)$
is a pseudo-differential operator, with

$\displaystyle
c(x,\xi)=(a\sharp_{h}b)(x,\xi):=\Big{(}e^{i\frac{h}{2}(\partial_{x}\cdot\partial_{\zeta}-\partial_{y}\cdot\partial_{\xi})}a(x,\xi)b(y,\zeta)\Big{)}\Big{|}_{y=x,\zeta=\xi}$

and $c\in S(u_{1}u_{2})$. We write $A\in\operatorname{Op}_{h}(S(u))$ to mean
that $A=\operatorname{Op}_{h}(a)$ for some $a\in S(u)$.

Following [11, 22, 25], we define the Hörmander class $S^{m}_{1,0}$ to be the
space of symbols $a(x,\xi)$ that satisfy

$\displaystyle|(\partial^{\alpha}_{\xi}\partial^{\beta}_{x}a)(x,\xi)|\leq
C_{\alpha,\beta}\langle\xi\rangle^{m-|\alpha|};\qquad\alpha,\beta\in\mathbb{N}^{d}.$
(B.3)

By [11, 23], we know that if $H$ satisfies (H1), then
$H:\mathcal{D}(H)\rightarrow\mathcal{H}$ is self-adjoint with
$\mathcal{D}(H)=\mathcal{H}^{m}$.

By [11], we know that if $a\in S^{m}_{1,0}$ is Hermitian-valued and satisfies
that the smallest singular value $|a_{\min}(x,\xi)|\geq
c\langle\xi,\zeta\rangle^{m}-1$ for some $c>0$, then for all $\Im(z)\neq 0$,
$R_{z}:=(z-\operatorname{Op}(a))^{-1}$ is a bijection of
$L^{2}(\mathbb{R}^{d})\otimes\mathbb{C}^{n}$ onto $\mathcal{H}^{m}$, and
$R_{z}=\operatorname{Op}(r_{z})$ with $r_{z}\in S^{-m}_{1,0}$. Defining
$z=:\lambda+i\omega$, we apply [16, Proposition 8.4] to obtain

$\displaystyle|\partial^{\alpha}_{x}\partial^{\beta}_{\xi}r_{z}(x,\xi;h)|\leq
C_{Z_{0},\alpha,\beta}\Big{(}1+\frac{\sqrt{h}}{|\omega|}\Big{)}^{2d+1}|\omega|^{-(|\alpha|+|\beta|)-1},\qquad
z\in Z_{0}$

for any compact set $Z_{0}\subset\mathbb{C}$. To derive the above result, it
helps to write $(z-H)^{-1}=(i-H)^{-1}(I+(z-i)(z-H)^{-1}),$ where
$(i-H)^{-1}\in S^{-m}_{1,0}$ is independent of $z$, and $(z-H)^{-1}$ is
bounded with $\left\lVert(z-H)^{-1}\right\rVert\leq\frac{C}{|\omega|}$. The
result from [11] does not address how symbolic bounds of $r_{z}$ depend on the
semi-classical parameter $h$. That is, given
$\operatorname{Op}(r_{z,h})=R_{z,h}:=(z-\operatorname{Op}_{h}(a))^{-1}$, we
would like bounds on $r_{z,h}$ that are uniform in $h$. This is where we use
[16, Chapter 8], which tells us there exists $h_{0}>0$ such that for all
$h\in(0,h_{0}]$, we have
$\operatorname{Op}(r_{z,h})=\operatorname{Op}_{h}(\tilde{r}_{z,h})$, with
$\tilde{r}_{z,h}\in S(\langle\xi\rangle^{-m})$. Since the symbolic bounds from
[11] must be continuous in $h$, this means $\tilde{r}_{z,h}\in
S(\langle\xi\rangle^{-m})$ uniformly in $h\in(0,1]$, as desired.

Given $\phi\in\mathcal{C}^{\infty}_{0}(\mathbb{R})$, there exists an almost
analytic extension $\tilde{\phi}\in\mathcal{C}^{\infty}_{0}(\mathbb{C})$ that
satisfies

$\displaystyle|\bar{\partial}\tilde{\phi}|\leq C_{N}|\Im z|^{N},$
$\displaystyle\quad
N\in\\{0,1,2,\dots\\};\qquad\tilde{\phi}(\lambda)=\phi(\lambda),\quad\lambda\in\mathbb{R}.$
(B.4)

We now recall [16, Theorem 8.1]. If $H$ is a self-adjoint operator on a
Hilbert space, then

$\displaystyle\phi(H)=-\frac{1}{\pi}\int\bar{\partial}\tilde{\phi}(z)(z-H)^{-1}d^{2}z,$
(B.5)

where $\bar{\partial}:=\frac{1}{2}\partial_{\Re z}+\frac{i}{2}\partial_{\Im
z}$ and $d^{2}z$ is the Lebesgue measure on $\mathbb{C}$. (B.5) is known as
the Helffer-Sjöstrand formula.

Suppose $u\in L^{1}(\mathbb{R}^{2d})$, and $|\partial^{\alpha}a(x,\xi;h)|\leq
C_{\alpha}u(x,\xi)$ for all $\alpha\in\mathbb{N}^{2d}$ and $h\in(0,1]$
(meaning that $a\in S(u)$). Then by [16, Theorem 9.4],
$\operatorname{Op}_{h}(a)$ is trace-class with
$\left\lVert\operatorname{Op}_{h}(a)\right\rVert_{1}\leq C\max_{|\alpha|\leq
2d+1}C_{\alpha}\left\lVert u\right\rVert_{L^{1}}$ and

$\displaystyle\operatorname{Tr}\operatorname{Op}_{h}(a)=\frac{1}{(2\pi
h)^{d}}\int_{\mathbb{R}^{2d}}\operatorname{tr}a(x,\xi;h)dxd\xi,$ (B.6)

where $C$ depends only on $d$ and $\operatorname{tr}$ is the standard matrix
trace.

###  B.2 Proofs

###### Proof of Proposition 4.3.

We have

$\displaystyle\sigma=\begin{pmatrix}\Omega
I+(\xi,\zeta)\cdot\sigma^{(\eta)}&\lambda U^{*}(x,y)\\\ \lambda U(x,y)&-\Omega
I+(\xi,\zeta)\cdot\sigma^{(\eta)}\end{pmatrix},$

with $U:\mathbb{R}^{2}\rightarrow\mathbb{R}$ smooth and bounded. Thus it is
clear that $\sigma$ is Hermitian-valued, with $\sigma\in S^{1}_{1,0}$ and
$|\sigma_{\min}(x,y,\xi,\zeta)|\geq c\langle\xi,\zeta\rangle-1$ for some
$c>0$. By assumption, $U(x,y)=A$ (resp. $U(x,y)=A^{*}$) whenever $f(x,y)$
(resp. $-f(x,y)$) is sufficiently large. Thus $\sigma_{\pm}$ from (H1) are
well defined, with the calculations in Appendix A proving that $\sigma_{\pm}$
both have a spectral gap at $0$. This completes the proof. ∎

###### Proof of Theorem 4.4.

We will prove that

$\displaystyle
2\pi\sigma_{I}(\tilde{H},P)=\operatorname{Tr}[U,P]U^{*}=\operatorname{Tr}[U,\bar{P}]U^{*}={\rm
Index}(\bar{P}U\bar{P}),$ (B.7)

where we use the shorthand $U:=U(\tilde{H})$.

We begin by proving the first equality in (B.7). Let $V:=U-I$, and observe
that $V\in\mathcal{C}^{\infty}_{c}(\\{\Phi=1\\}^{\circ})$. We have

$\displaystyle\operatorname{Tr}[U,P]U^{*}=\operatorname{Tr}[V,P]U^{*}=\operatorname{Tr}[V,P]V^{*}+\operatorname{Tr}[V,P].$

Using that
$[(z-\tilde{H})^{-1},P]=(z-\tilde{H})^{-1}[\tilde{H},P](z-\tilde{H})^{-1}$,
the Helffer-Sjöstrand formula (B.5) and cyclicity of the trace [24] imply that

$\displaystyle\operatorname{Tr}[V,P]V^{*}=\operatorname{Tr}\Big{(}-\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}\tilde{V}(z)[\tilde{H},P](z-\tilde{H})^{-2}d^{2}zV^{*}(\tilde{H})\Big{)}.$

After integrating by parts in $\partial$ and using that
$\partial\tilde{V}=\tilde{V^{\prime}}$ for some almost analytic extension
$\tilde{V^{\prime}}$ of $V^{\prime}$, we see that

$\displaystyle\operatorname{Tr}[V,P]V^{*}=\operatorname{Tr}[\tilde{H},P]V^{\prime}V^{*}=2\pi\sigma_{I}(\tilde{H},P)-\operatorname{Tr}2\pi
i[\tilde{H},P]\varphi^{\prime}U=2\pi\sigma_{I}(\tilde{H},P)-\operatorname{Tr}[\tilde{H},P]V^{\prime}.$
(B.8)

We have thus shown that

$\displaystyle\operatorname{Tr}[U,P]U^{*}=2\pi\sigma_{I}(\tilde{H},P)-\operatorname{Tr}[\tilde{H},P]V^{\prime}+\operatorname{Tr}[V,P].$

By the same logic used in (B.8), we obtain that
$\operatorname{Tr}[\tilde{H},P]V^{\prime}=\operatorname{Tr}[\tilde{H},P]V^{\prime}\Phi=\operatorname{Tr}[V,P]\Phi$,
hence

$\displaystyle\operatorname{Tr}[U,P]U^{*}=2\pi\sigma_{I}(\tilde{H},P)+\operatorname{Tr}[V,P](1-\Phi).$

Again using cyclicity of the trace, this means

$\displaystyle\operatorname{Tr}[U,P]U^{*}=2\pi\sigma_{I}(\tilde{H},P)+\operatorname{Tr}\Psi_{1}[V,P]\Psi_{2},$

for some $\Psi_{j}=\psi_{j}(\tilde{H})$, where
$\psi_{j}\in\mathcal{C}^{\infty}(\mathbb{R})$ vanishes on
$\operatorname{supp}(V)$. The first equality in (B.7) follows.

We now prove the second equality of (B.7). Let $\chi(x,y)=\chi_{1}(g(x,y))$
where $\chi_{1}\in\mathcal{C}_{c}^{\infty}$ such that
$\chi(\bar{P}-P)=\bar{P}-P$. Then

$\displaystyle[U,\bar{P}-P]U^{*}=[V,\bar{P}-P]U^{*}=V\chi(\bar{P}-P)U^{*}-(\bar{P}-P)\chi
VU^{*}$

is trace class, as $V\chi$ and $\chi V$ are both trace-class by the $\Psi$DO
calculus with $P-P_{1}$ and $U^{*}$ bounded. Since we know that $[U,P]U^{*}$
is trace-class, this proves that

$[U,\bar{P}]U^{*}=[U,\bar{P}-P]U^{*}+[U,P]U^{*}$

is trace-class. Now,

$\displaystyle\operatorname{Tr}[U,\bar{P}-P]U^{*}=\operatorname{Tr}V\chi(\bar{P}-P)U^{*}-\operatorname{Tr}(\bar{P}-P)\chi
VU^{*},$

with

$\displaystyle\operatorname{Tr}V\chi(\bar{P}-P)U^{*}=\operatorname{Tr}(\bar{P}-P)U^{*}V\chi$

and

$\displaystyle\operatorname{Tr}(\bar{P}-P)\chi
VU^{*}=\operatorname{Tr}\chi(\bar{P}-P)\chi
VU^{*}=\operatorname{Tr}(\bar{P}-P)\chi
VU^{*}\chi=\operatorname{Tr}(\bar{P}-P)VU^{*}\chi.$

Since $[V,U^{*}]=0$, we have proven that
$\operatorname{Tr}[U,\bar{P}-P]U^{*}=0$, which verifies the second equality of
(B.7).

Finally, the last equality of (B.7) follows immediately from [2, Proposition
2.4] and the fact that $[U,\bar{P}]U^{*}$ is trace-class. ∎

###### Proof of Theorem 4.6.

For $\mu\in[0,1]$, define $\varphi_{\mu}:=\varphi+\mu(\varphi_{1}-\varphi)$.
Theorem 4.4 implies that
$\operatorname{Tr}i[\tilde{H},P]\varphi^{\prime}_{\mu}(\tilde{H})={\rm
Index}(\bar{P}U_{\mu}\bar{P})$ for all $\mu\in[0,1]$, where
$U_{\mu}:=e^{i2\pi\varphi_{\mu}(\tilde{H})}$. With $V_{\mu}:=U_{\mu}-I$, the
Helffer-Sjöstrand formula (B.5) implies

$\displaystyle
U_{\mu_{2}}-U_{\mu_{1}}=V_{\mu_{2}}-V_{\mu_{1}}=-\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}(\tilde{V}_{\mu_{2}}(z)-\tilde{V}_{\mu_{1}}(z))(z-\tilde{H})^{-1}d^{2}z.$

Since $|\bar{\partial}(\tilde{V}_{\mu_{2}}(z)-\tilde{V}_{\mu_{1}}(z))|\leq
C|\mu_{2}-\mu_{1}|$ uniformly in $\mu_{1},\mu_{2}\in[0,1]$ and
$z\in\mathbb{C}$, it follows that $\left\lVert
U_{\mu_{2}}-U_{\mu_{1}}\right\rVert\rightarrow 0$ as
$\mu_{2}-\mu_{1}\rightarrow 0$. Therefore, by [22, Theorem 19.1.5], ${\rm
Index}(\bar{P}U_{\mu}\bar{P})$ is independent of $\mu\in[0,1]$, and the result
is complete. ∎

We now state the following regularity result.

######  Proposition B.1.

Suppose $H=\operatorname{Op}(\sigma)$ satisfies (H1), define

$H_{h}:=\operatorname{Op}_{h}(\sigma(x,y,\xi,\zeta))\mbox{ for }\ 0<h\leq 1,$
(B.9)

and let $\Phi\in\mathcal{C}^{\infty}_{c}(E_{1},E_{2})$. Then for all
$P\in\mathfrak{S}(0,1;g(x,y))$, we have

$\displaystyle\Phi(H_{h})\in\operatorname{Op}_{h}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))\qquad\text{and}\qquad[H_{h},P]\Phi(H_{h})\in\operatorname{Op}_{h}(S(\langle
x,y,\xi,\zeta\rangle^{-\infty})).$

By the paragraph above (B.6), Proposition B.1 immediately implies that
$\sigma_{I}(H_{h},P)$ is well defined. By Theorem 4.4, Proposition B.1 also
implies that $2\pi\sigma_{I}(H)={\rm Index}(\bar{P}U(H)\bar{P})$ for any $H$
satisfying (H1).

###### Proof.

We emulate the proofs of [30, Lemmas 3.3 and 3.4].

For any $p>0$, we can write $\Phi(H_{h})=(i-H_{h})^{-p}\Phi_{p}(H_{h})$ with
$\Phi_{p}\in\mathcal{C}^{\infty}_{c}(E_{1},E_{2})$. Since $H$ satisfies (H1),
we know that $(i-H_{h})^{-1}\in S^{-m}_{1,0}$ (see Appendix B.1 for more
details). It follows from the composition calculus that
$\Phi(H_{h})\in\operatorname{Op}_{h}(S(\langle\xi,\zeta\rangle^{-\infty}))$.
By assumption on $\sigma_{\pm}$, we know that
$H_{\pm,h}:=\operatorname{Op}_{h}(\sigma_{\pm})$ has a spectral gap in
$(E_{1},E_{2})$ for all $h\in(0,1]$, hence $\Phi(H_{\pm,h})=0$. We can thus
write
$\Phi(H_{h})=\phi(x,y)(\Phi(H_{h})-\Phi(H_{+,h}))+(1-\phi(x,y))(\Phi(H_{h})-\Phi(H_{-,h}))$,
for some $\phi\in\mathfrak{S}(0,1;f)$. The Helffer-Sjöstrand formula (B.5)
implies that

$\displaystyle\Phi(H_{h})-\Phi(H_{+,h})=\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}\tilde{\Phi}(z)(z-H_{h})^{-1}(H_{h}-H_{+,h})(z-H_{+,h})^{-1}d^{2}z.$

Since $\sigma-\sigma_{+}$ vanishes whenever $f(x,y)$ is sufficiently large, it
follows that $\Phi(H_{h})-\Phi(H_{+,h})\in\operatorname{Op}_{h}(S(\langle
f_{+}(x,y)\rangle^{-\infty}))$, where $f_{+}:=\max\\{f,0\\}$. Since $\phi$
vanishes whenever $-f(x,y)$ is sufficiently large, we conclude that
$\phi(x,y)(\Phi(H_{h})-\Phi(H_{+,h}))\in\operatorname{Op}_{h}(S(\langle
f(x,y)\rangle^{-\infty}))$. The same reasoning shows that
$(1-\phi(x,y))(\Phi(H_{h})-\Phi(H_{-,h}))\in\operatorname{Op}_{h}(S(\langle
f(x,y)\rangle^{-\infty}))$. We have thus shown that

$\Phi(H_{h})\in\operatorname{Op}_{h}(S(\langle\xi,\zeta\rangle^{-\infty})\cap
S(\langle f(x,y)\rangle^{-\infty}))=\operatorname{Op}_{h}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty})),$

which proves the first result.

For the second claim, observe that
$[H_{h},P]=(1-P)H_{h}P-PH_{h}(1-P)\in\operatorname{Op}_{h}(S(\langle
g(x,y)\rangle^{-\infty}\langle\xi,\zeta\rangle^{m}))$. The composition
calculus and (4.6) then imply that
$[H_{h},P]\Phi(H_{h})\in\operatorname{Op}_{h}(S(\langle
x,y,\xi,\zeta\rangle^{-\infty})),$ and the proof is complete. ∎

###### Proof of Theorem 4.7.

For $h\in(0,1]$, define $U_{h}:=e^{i2\pi\varphi(H_{h})}$. Proposition B.1 and
Theorem 4.4 imply that $2\pi\sigma_{I}(H_{h},P)={\rm
Index}(\bar{P}U_{h}\bar{P})$ for all $h\in(0,1]$. It thus suffices to show
that for any fixed $h\in(0,1]$, ${\rm Index}(\bar{P}U_{h^{\prime}}\bar{P})$ is
constant over $h^{\prime}$ in an open neighborhood of $h$. Using the Helffer-
Sjöstrand formula (B.5), we write

$\displaystyle
U_{h^{\prime}}-U_{h}=\frac{1}{\pi}\int_{\mathbb{R}}\bar{\partial}\tilde{V}(z)(z-H_{h^{\prime}})^{-1}(H_{h^{\prime}}-H_{h})(z-H_{h})^{-1}d^{2}z.$

Since $\left\lVert(z-H_{h^{\prime}})^{-1}\right\rVert\leq|\Im z|^{-1}$ and
$H_{h^{\prime}}-H_{h}\in\operatorname{Op}(|h^{\prime}-h|S(\langle\xi,\zeta\rangle^{m}))$
with $(z-H_{h})^{-1}\in\operatorname{Op}(S(\langle\xi,\zeta\rangle^{-m}))$ for
all $\Im z\neq 0$ (with bounds growing at most algebraically in $|\Im
z|^{-1}$), the rapid decay of $\bar{\partial}\tilde{V}(z)$ near the real axis
implies that $\left\lVert U_{h^{\prime}}-U_{h}\right\rVert\leq
C|h^{\prime}-h|$. The result then follows from [22, Theorem 19.1.5]. ∎

###### Proof of Theorem 4.8.

We will apply Theorem 4.4. First, observe that
$H^{(\mu)}=\operatorname{Op}(\sigma^{(\mu)})$ with
$\sigma^{(\mu)}\in\operatorname{Op}(S^{m}_{1,0})$. Moreover, (H1) and the
assumption $W\in\operatorname{Op}(S^{m}_{1,0})$ implies that
$|\sigma^{(\mu)}_{\min}(x,y,\xi,\zeta)|\geq c\langle\xi,\zeta\rangle^{m}-1$
whenever $\mu>0$ is sufficiently small, where $\sigma^{(\mu)}_{\min}$ denotes
the smallest magnitude eigenvalue of $\sigma^{(\mu)}$. Since $\sigma^{(\mu)}$
is Hermitian-valued, it follows that whenever $\mu$ is sufficiently small,
$H^{(\mu)}$ is self-adjoint with the same domain of definition
$\mathcal{D}(H^{(\mu)})=\mathcal{D}(H)=\mathcal{H}^{m}$, and
$(z-H^{(\mu)})^{-1}\in\operatorname{Op}(S^{-m}_{1,0})$ for all $\Im z\neq 0$.
We refer to [11] for more details, particularly Corollary 2 and the paragraph
following Theorem 3.

Next, we verify that $\Phi(H^{(\mu)})\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))$ for all $\mu>0$ sufficiently small. Let
$\Phi_{0}\in\mathcal{C}^{\infty}_{c}(E_{1},E_{2})$ such that
$\Phi\in\mathcal{C}^{\infty}_{c}(\\{\Phi_{0}=1\\}^{\circ})$. With
$\Theta:=\Phi(H^{(\mu)})-\Phi(H)$ and
$\Theta_{0}:=\Phi_{0}(H^{(\mu)})-\Phi_{0}(H)$, it follows that

$\displaystyle\Phi(H^{(\mu)})=\Phi(H)+\Theta,\qquad\Theta=\Theta\Theta_{0}+\Phi(H)\Theta_{0}+\Theta\Phi_{0}(H).$

By assumption on $W$, $1-\Theta_{0}$ is invertible whenever $\mu>0$ is
sufficiently small, with $(1-\Theta_{0})^{-1}\in\operatorname{Op}(S(1))$. By
Proposition B.1,
$\Theta=(\Phi(H)\Theta_{0}+\Theta\Phi_{0}(H))(1-\Theta_{0})^{-1}\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))$ and thus
$\Phi(H^{(\mu)})\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))$ for all $\mu>0$ sufficiently small, as
desired.

We have shown that $H^{(\mu)}$ satisfies the assumptions of Theorem 4.4, hence
$2\pi\sigma_{I}(H^{(\mu)},P)={\rm Index}(\bar{P}U^{(\mu)}\bar{P})$ whenever
$\mu$ is sufficiently small, with $U^{(\mu)}:=e^{i2\pi\varphi(H^{(\mu)})}$. As
in the proof of Theorem 4.7 above, we use the Helffer-Sjöstrand formula and
the property

$\displaystyle(z-H^{(\mu)})^{-1}\in\operatorname{Op}(S(\langle\xi,\zeta\rangle^{-m})),\qquad\Im
z\neq 0$

to verify that $U^{(\mu)}$ is continuous in $\mu$ (with respect to operator
norm). The result follows. ∎

###### Proof of Theorem 4.9.

We again must show that $H^{(\mu)}$ satisfies the assumptions of Theorem 4.4,
this time for all $\mu\in[0,1]$. The hypothesis (H1) and our growth constraint
on $W$ imply that $H^{(\mu)}=\operatorname{Op}(\sigma^{(\mu)})$ with
$\sigma^{(\mu)}$ Hermitian-valued and satisfying $\sigma^{(\mu)}\in
S^{m}_{1,0}$ and $|\sigma^{(\mu)}_{\min}(x,y,\xi,\zeta)|\geq
c\langle\xi,\zeta\rangle^{m}-1$. As in the proof of Theorem 4.8, this means
$H^{(\mu)}$ is self-adjoint with domain of definition
$\mathcal{D}(H^{(\mu)})=\mathcal{H}^{m}$, and
$(z-H^{(\mu)})^{-1}\in\operatorname{Op}(S^{-m}_{1,0})$ whenever $\Im z\neq 0$
[11].

To obtain decay of the symbol of $\Phi(H^{(\mu)})$, we again use the identity

$\displaystyle\Theta(1-\Theta_{0})=\Phi(H)\Theta_{0}+\Theta\Phi_{0}(H),$
(B.10)

with $\Theta$ and $\Theta_{0}$ defined in the proof of Theorem 4.8.
Proposition B.1 and (B.10) imply that
$\Theta(1-\Theta_{0})\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty}))$. By assumption on $W$, an application of
the Helffer-Sjöstrand formula reveals that
$\Theta_{0}\in\operatorname{Op}(S(\langle x,y,\xi,\zeta\rangle^{-\delta}))$.
Thus there exist $\Theta_{00}\in\mathcal{C}^{\infty}_{c}(\mathbb{R}^{4})$ and
$\Theta_{01}\in S(1)$ as small as necessary such that
$\Theta_{0}=\Theta_{00}+\Theta_{01}$. It follows that

$\displaystyle\Theta(1-\Theta_{01})=\Theta(1-\Theta_{0})+\Theta\Theta_{00}\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty})),$

hence

$\displaystyle\Theta=(\Theta(1-\Theta_{0})+\Theta\Theta_{00})(1-\Theta_{01})^{-1}\in\operatorname{Op}(S(\langle
f(x,y),\xi,\zeta\rangle^{-\infty})).$

We have shown that $\Theta$, and therefore $\Phi(H^{(\mu)})$, has symbol with
decay required by Theorem 4.4.

As in the proof of Theorem 4.8, we can now conclude that
$2\pi\sigma_{I}(H^{(\mu)},P)={\rm Index}(\bar{P}U^{(\mu)}\bar{P})$ for all
$\mu\in[0,1]$, and the result follows. ∎

We now prove Theorem 4.1, which states a bulk-interface correspondence. The
strategy is to use an asymptotic expansion in the semiclassical parameter $h$
and apply Theorem 4.7 to eliminate terms that are not $O(1)$. A similar
technique is used in [6, 5, 29].

###### Proof of Theorem 4.1.

By Theorem 4.6, $\sigma_{I}(H,P)$ is independent of
$\varphi\in\mathfrak{S}(0,1;E_{1},E_{2})$, thus we can take
$\varphi^{\prime}\in\mathcal{C}^{\infty}_{c}(\alpha_{0},\alpha)$ for some
$\alpha_{0}>E_{1}$. We will use the shorthand $\sigma_{I}:=\sigma_{I}(H,P)$.

Let $\operatorname{Op}_{h}\nu_{h}:=\varphi^{\prime}(H_{h})$. Using [16,
equation (7.19) and the preceding paragraph] and [39, Theorems 3.14 & 4.16 and
their proofs] as in [29], we verify that

$\displaystyle\nu_{h}+\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\tilde{q}_{z,h}d^{2}z\in
S^{-2}(\langle f(x,y),\xi,\zeta\rangle^{-\infty}),$

where

$\displaystyle\tilde{q}_{z,h}=\sigma_{z}^{-1}+\frac{ih}{2}\\{\sigma_{z}^{-1},\sigma_{z}\\}\sigma_{z}^{-1},\qquad\\{a,b\\}:=\partial_{\xi}a\partial_{x}b+\partial_{\zeta}a\partial_{y}b-\partial_{x}a\partial_{\xi}b-\partial_{y}a\partial_{\zeta}b$

and $\sigma_{z}:=z-\sigma$. With
$\operatorname{Op}_{h}(\kappa_{h}):=[H_{h},P]$, we have that

$\displaystyle\kappa_{h}+ihk_{1}-\frac{h^{2}}{4}k_{2}\in S^{-3}(\langle
g(x,y)\rangle^{-\infty}\langle\xi,\zeta\rangle^{m}),$

where

$\displaystyle
k_{1}:=\partial_{\xi}\sigma\partial_{x}P+\partial_{\zeta}\sigma\partial_{y}P,\qquad
k_{2}:=\partial_{\xi\xi}\sigma\partial_{xx}P+2\partial_{\xi\zeta}\sigma\partial_{xy}P+\partial_{\zeta\zeta}\sigma\partial_{yy}P.$
(B.11)

Since $\nu_{h}\in S(\langle f(x,y),\xi,\zeta\rangle^{-\infty})$ and
$\kappa_{h}\in S^{-1}(\langle
g(x,y)\rangle^{-\infty}\langle\xi,\zeta\rangle^{m})$, the composition calculus
implies that

$\displaystyle\kappa_{h}\sharp_{h}\nu_{h}-\kappa_{h}\nu_{h}+\frac{ih}{2}\\{\kappa_{h},\nu_{h}\\}\in
S^{-3}(\langle x,y,\xi,\zeta\rangle^{-\infty}),$

with $S^{-3}$ (rather than $S^{-2}$) above because $\kappa_{h}$ is $O(h)$ in
$S(\langle g(x,y)\rangle^{-\infty}\langle\xi,\zeta\rangle^{m})$. Therefore,

$\displaystyle\sigma_{I}$ $\displaystyle=\frac{i}{(2\pi
h)^{2}}\operatorname{tr}\int_{\mathbb{R}^{4}}\kappa_{h}\sharp_{h}\nu_{h}dR_{4}=\frac{i}{(2\pi
h)^{2}}\operatorname{tr}\int_{\mathbb{R}^{4}}\Big{(}\kappa_{h}\nu_{h}-\frac{ih}{2}\\{\kappa_{h},\nu_{h}\\}\Big{)}dR_{4}+o(1)$

as $h\rightarrow 0$, with $dR_{4}:=dxdyd\xi d\zeta$. Since

$\displaystyle\kappa_{h}\nu_{h}=(ihk_{1}-\frac{h^{2}}{4}k_{2})\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\Big{(}\sigma_{z}^{-1}+\frac{ih}{2}\\{\sigma_{z}^{-1},\sigma_{z}\\}\sigma_{z}^{-1}\Big{)}d^{2}z+h^{3}a_{h},\qquad
a_{h}\in S(\langle x,y,\xi,\zeta\rangle^{-\infty})$

and

$\displaystyle\\{\kappa_{h},\nu_{h}\\}=\Big{\\{}ihk_{1},\frac{1}{\pi}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\sigma_{z}^{-1}d^{2}z\Big{\\}}+h^{2}b_{h},\qquad
b_{h}\in S(\langle x,y,\xi,\zeta\rangle^{-\infty}),$

it follows that

$\displaystyle\sigma_{I}$ $\displaystyle=\frac{i}{(2\pi
h)^{2}}\frac{1}{\pi}\operatorname{tr}\int_{\mathbb{R}^{4}}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\Big{(}ihk_{1}\sigma_{z}^{-1}-\frac{h^{2}}{2}k_{1}\\{\sigma_{z}^{-1},\sigma_{z}\\}\sigma_{z}^{-1}$
$\displaystyle\hskip
85.35826pt-\frac{h^{2}}{4}k_{2}\sigma_{z}^{-1}+\frac{h^{2}}{2}\\{k_{1},\sigma_{z}^{-1}\\}\Big{)}d^{2}zdR_{4}+o(1)$

as $h\rightarrow 0$. Since $\sigma_{I}$ is independent of $h$, it follows that
the $O(h^{-1})$ term above vanishes, and thus

$\displaystyle\sigma_{I}=\frac{i}{(2\pi)^{3}}\operatorname{tr}\int_{\mathbb{R}^{4}}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\Big{(}-k_{1}\sigma_{z}^{-1}\\{\sigma_{z},\sigma_{z}^{-1}\\}-\frac{1}{2}k_{2}\sigma_{z}^{-1}+\\{k_{1},\sigma_{z}^{-1}\\}\Big{)}d^{2}zdR_{4}.$
(B.12)

Observe that whenever $\langle f(x,y),\xi,\zeta\rangle$ is sufficiently large,
$z\mapsto\sigma_{z}^{-1}$ is holomorphic and thus the above integral over $z$
vanishes (this is verified via an integration by parts in $\bar{\partial}$).
Using that $k_{1}$ and $k_{2}$ vanish whenever $\langle g(x,y)\rangle$ is
sufficiently large, we can replace the above integration limit
$\mathbb{R}^{4}$ by a sufficiently large rectangle $B:=B_{xy}\times
B_{\xi\zeta}\subset\mathbb{R}^{2}\times\mathbb{R}^{2}$. We require that $B$
contain all points $(x_{0}(t),y_{0}(t),\xi,\zeta)$ for which
$\sigma(x_{0}(t),y_{0}(t),\xi,\zeta)$ has an eigenvalue of $\alpha$. Moreover,
assume that $B_{xy}$ contains $(x_{0}(t_{0}),y_{0}(t_{0}))$ and
$(x_{0}(-t_{0}),y_{0}(-t_{0}))$, where we recall the definitions of $t_{0}$
above Theorem 4.1 and $(x_{0},y_{0})$ below (4.6).

At this point, $\\{k_{1},\sigma_{z}^{-1}\\}$ can be written in divergence form
and the corresponding term converted (via integration by parts in
$(x,y,\xi,\zeta)$) to an integral over the surface $\partial B$. An
integration by parts in $\bar{\partial}$ then reveals that the contribution of
this term vanishes. Thus we are left with

$\displaystyle\sigma_{I}=-\frac{i}{(2\pi)^{3}}\operatorname{tr}\int_{B}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)\Big{(}k_{1}\sigma_{z}^{-1}\\{\sigma_{z},\sigma_{z}^{-1}\\}+\frac{1}{2}k_{2}\sigma_{z}^{-1}\Big{)}d^{2}zdR_{4}.$
(B.13)

We now simplify the first term above. Integrating by parts in $\bar{\partial}$
with $z=:\lambda+i\omega$, we see that

$\displaystyle\sigma_{I,1}:=\frac{-i}{(2\pi)^{3}}\operatorname{tr}\int_{B}\int_{\mathbb{C}}\bar{\partial}\tilde{\varphi^{\prime}}(z)k_{1}\sigma_{z}^{-1}\\{\sigma_{z},\sigma_{z}^{-1}\\}d^{2}zdR_{4}=\frac{-1}{2(2\pi)^{3}}\operatorname{tr}\int_{B}\int_{\alpha_{0}}^{\alpha}\varphi^{\prime}(\lambda)k_{1}\sigma_{z}^{-1}\\{\sigma_{z},\sigma_{z}^{-1}\\}\Big{|}^{\omega=0^{+}}_{\omega=0^{-}}d\lambda
dR_{4}.$

Let $\mathcal{L}:=\\{x_{0}(t),y_{0}(t):t\in\mathbb{R}\\}$ be the range of
$(x_{0},y_{0})$. For simplicity, assume that $\mathcal{L}\subset g^{-1}(0)$.
By cyclicity of the trace and recalling the definition (B.11) of $k_{1}$, it
follows that

$\displaystyle\operatorname{tr}k_{1}\\{\sigma_{z}^{-1},\sigma_{z}\\}\sigma_{z}^{-1}=-\operatorname{tr}(\partial_{x}P\partial_{\xi}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\zeta,y}+\partial_{y}P\partial_{\zeta}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\xi,x})\sigma_{z}^{-1},$

where we have defined
$\\{a,b\\}_{\alpha,\beta}:=\partial_{\alpha}a\partial_{\beta}b-\partial_{\beta}a\partial_{\alpha}b$
and used the fact that
$(\partial_{\xi}\sigma,\partial_{\zeta}\sigma)=-(\partial_{\xi}\sigma_{z},\partial_{\zeta}\sigma_{z})$.
Integrating by parts in $x$ (first term) and $y$ (second term), we obtain

$\displaystyle\sigma_{I,1}$
$\displaystyle=-\frac{1}{2(2\pi)^{3}}\operatorname{tr}\int_{B}\int_{\alpha_{0}}^{\alpha}\varphi^{\prime}(\lambda)P\Big{(}\partial_{x}(\partial_{\xi}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\zeta,y}\sigma_{z}^{-1})+\partial_{y}(\partial_{\zeta}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\xi,x}\sigma_{z}^{-1})\Big{)}\Big{|}^{\omega=0^{+}}_{\omega=0^{-}}d\lambda
dR_{4}$
$\displaystyle\qquad+\frac{1}{2(2\pi)^{3}}\operatorname{tr}\int_{B_{\xi\zeta}}\int_{B_{xy}\cap\mathcal{L}}\int_{\alpha_{0}}^{\alpha}\varphi^{\prime}(\lambda)(\partial_{\xi}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\zeta,y}\nu_{x}+\partial_{\zeta}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\xi,x}\nu_{y})\sigma_{z}^{-1}\Big{|}^{\omega=0^{+}}_{\omega=0^{-}}d\lambda
d\ell dR_{2}$ $\displaystyle=:\sigma_{I,10}+\sigma_{I,11},$

where $dR_{2}:=d\xi d\zeta$ and $d\ell$ is the integration measure on
$\mathcal{L}$. Here, $\nu$ is the unit vector (outwardly) normal to the
surface $\partial(\\{g(x,y)\leq 0\\}\cap B_{xy})$. Note that the other surface
terms do not contribute; over $(\partial(\\{g(x,y)\leq 0\\}\cap
B_{xy}))\setminus(B_{xy}\cap\mathcal{L})$, either $P=0$ or the map
$z\mapsto\sigma_{z}^{-1}$ is holomorphic and thus the difference between
$\omega=0^{+}$ and $\omega=0^{-}$ vanishes.

We now verify that the volume term $\sigma_{I,10}$ vanishes. First, observe
that

$\displaystyle\int_{B}\int_{\alpha_{0}}^{\alpha}\varphi^{\prime}(\lambda)P\Big{(}\partial_{\xi}(\partial_{x}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\zeta,y}\sigma_{z}^{-1})+\partial_{\zeta}(\partial_{y}\sigma_{z}\\{\sigma_{z}^{-1},\sigma_{z}\\}_{\xi,x}\sigma_{z}^{-1})\Big{)}\Big{|}^{\omega=0^{+}}_{\omega=0^{-}}d\lambda
dR_{4}=0,$

as we can use the fact that $P$ is independent of $(\xi,\zeta)$ to integrate
the above left-hand side by parts in $\xi$ (first term) and $\zeta$ (second
term) to obtain integrals over $\partial B_{\xi\zeta}$, over which
$z\mapsto\sigma_{z}^{-1}$ is holomorphic. Therefore,