# Asymmetric transport computations in Dirac models of topological insulators

Guillaume Bal Departments of Statistics and Mathematics and CCAM, University
of Chicago, Chicago, IL 60637<EMAIL_ADDRESS>Jeremy G Hoskins
Department of Statistics and CCAM, University of Chicago, Chicago, IL 60637;
<EMAIL_ADDRESS>Zhongjian Wang Department of Statistics and CCAM,
University of Chicago, Chicago, IL 60637<EMAIL_ADDRESS>
###### Abstract

This paper presents a fast algorithm for computing transport properties of
two-dimensional Dirac operators with linear domain walls, which model the
macroscopic behavior of the robust and asymmetric transport observed at an
interface separating two two-dimensional topological insulators. Our method is
based on reformulating the partial differential equation as a corresponding
volume integral equation, which we solve via a spectral discretization scheme.

We demonstrate the accuracy of our method by confirming the quantization of an
appropriate interface conductivity modeling transport asymmetry along the
interface, and moreover confirm that this quantity is immune to local
perturbations. We also compute the far-field scattering matrix generated by
such perturbations and verify that while asymmetric transport is topologically
protected the absence of back-scattering is not.

##  1 Introduction

Two-dimensional Dirac operators appear naturally in the analysis of
topological phases of matter [1, 2] and one-particle models of topological
insulators and topological superconductors [3, 4, 5, 6, 7, 8, 9, 10]. This
paper focuses on efficient numerical simulation of the transport properties of
one such operator, which for concreteness we take as

$\displaystyle
H_{V}=D_{x}\sigma_{1}+D_{y}\sigma_{2}+m(y)\sigma_{3}+V=\begin{pmatrix}m(y)+V_{11}&D_{x}-iD_{y}+V_{12}\\\
D_{x}+iD_{y}+V_{12}^{*}&-m(y)+V_{22}\end{pmatrix}.$ (1)

Here, $\sigma_{1,2,3}$ are the standard Pauli matrices (which along with the
identity matrix $\sigma_{0}=I_{2}$ form a basis of $2\times 2$ Hermitian
matrices), $D_{a}=-i\partial_{a}$ for $a=x,y$, $m(y)=y$ is a linear domain
wall, and $V(x,y)$ is an arbitrary local Hermitian-valued perturbation with
compact support. Physically, the scalar component of $V$ corresponds to an
electric potential, the components of $V$ in front of $\sigma_{1,2}$ to a
magnetic potential, and the component in front of $\sigma_{3}$ to a local
perturbation of the domain wall $m(y)=y$. We also use the notation $H_{V}=H+V$
with $H$ the unperturbed Dirac operator.

The domain wall $m(y)$ models a topological phase transition between two bulk
phases [3, 7, 9], one in the domain $y>0$ where $m(y)>0$ and one in the domain
$y<0$ with $m(y)<0$. A striking feature of topologically non-trivial materials
is that transport (say, electronic) along the interface $y\sim 0$ is
asymmetric, in the sense that we observe stronger transport towards negative
values of $x$ than transport in the opposite direction. Moreover, this
transport asymmetry is quantized and related to a difference of bulk
topological invariants following the celebrated bulk-edge correspondence [3,
4] (more precisely in fact a bulk-difference invariant [11] in the context of
Dirac operators). Finally, the topological nature of the asymmetry ensures its
robustness to perturbations. In our context, the transport asymmetry is
independent of any perturbation $V(x,y)$ (no matter how large or how
oscillatory) that vanishes at infinity sufficiently rapidly.

The robustness of this asymmetric transport is one of the main practical
appeals of such materials. It stands in sharp contrast to topologically
trivial materials, for instance modeled by a one-dimensional wave equation
along the interface $y=0$, where transmission is exponentially suppressed by
Anderson localization in the presence of strong perturbations [4, 12, 13] (so
that incoming signals may be fully back-scattered). However, robust asymmetric
transport does not imply the absence of back-scattering; rather, merely that
some quantized transmission is guaranteed. See [4, 14, 15, 16] for references
on topologically-protected asymmetric transport in numerous applications.

In order to understand the transport properties of models such as (1), we
develop an efficient algorithm to compute its generalized eigenvectors. This
allows us to estimate a conductivity quantifying asymmetric transport and
scattering matrices generated by the perturbation $V$. For $\psi_{\rm in}$ an
incoming plane wave solution of $(H-E)\psi_{\rm in}=0$, we look for outgoing
solutions of

$\displaystyle(H+V-E)\psi_{\rm out}=-V\psi_{\rm in}.$ (2)

We will construct explicitly the (unique) outgoing Green’s function, i.e., the
kernel of the operator $(H-E)^{-1}\equiv(H-(E-i0^{+}))^{-1}$. The above
equation can then be reduced to a volume integral equation for a new unknown
$\rho$ which we will refer to as the density corresponding to $\psi_{\rm
out}.$ Specifically,

$\displaystyle\psi_{\rm
out}=(H-E)^{-1}\rho\quad\implies\quad\rho+V(H-E)^{-1}\rho=-V\psi_{\rm in}.$
(3)

A computational advantage of the latter integral equation is that by
construction, $\rho(x,y;E)$ is supported on the same domain as $V$. Numerical
methods to solve such Fredholm integral equations of the second kind are
discussed, e.g., in [17].

The Green’s function $G(x,y;x_{0},y_{0};E)$, the $2\times 2$ matrix-valued
kernel of $(H-(E-i0^{+}))^{-1}$, does not admit a simple, explicit expression.
It may, however, be efficiently written in a basis of Hermite functions in the
$y-$variable. The construction is detailed in section 2.2. The Green’s
function admits two main features. One in the singularity at the source
location where the Green’s function is inversely proportional to the distance
between $(x,y)$ and $(x_{0},y_{0})$ as for any two-dimensional Dirac operator.
Another feature is the behavior of the Green’s function ‘at infinity’ along
the $x$ axis reflecting the wave-guide nature of the model generated by the
domain wall as well as the transport asymmetry along the edge $y\sim 0$.

For $V\in C^{\infty}_{c}(\mathbb{R}^{2};\mathbb{M}_{2})$ a smooth compactly
supported matrix-valued perturbation, solving the volume integral equation for
$\rho$ requires inverting an operator of the form $I+K$ with $K$ a compact
operator (from $H^{k}(\mathbb{R}^{2};\mathbb{C}^{2})$ to itself for each
integer $k\geq 0$). When $-1$ is not an eigenvalue of $K$ (this is independent
of $k$), we thus obtain that the unique solution $\rho$ of (3) is itself
smooth in $C^{\infty}_{c}(\mathbb{R}^{2};\mathbb{C}^{2})$. The numerical
simulation of (3) thus boils down to efficiently approximating the integral
operator $V(x,y)\int_{X}G(x,y;x_{0},y_{0})\rho(x_{0},y_{0})dx_{0}dy_{0}$ for
$X$ a domain including the support of $V$.

Due to the nature of the problem, it is natural to discretize $\rho$ using a
finite Cartesian product basis of Legendre polynomials in the $x$ direction
and Hermite functions in the $y$ direction. The details of the construction
and demonstration of convergence are presented in section 3.1.

We accelerate the computation of the density $\rho$ using a hierarchical
merging approach in the $x$-direction. First, the $x-$axis is subdivided into
adjacent small slabs which we refer to as leaves. The full transmission-
reflection matrix (solution operator) for each small leaf is then computed
explicitly. A standard algebraic merging procedure, described in section 3.2,
is then used to solve for $\rho(x,y)$ globally. This approach is a slight
generalization of the method in [18] for solving two point boundary value
problems. More broadly, it is closely related to fast direct solver methods
for constructing compressed representations of Green’s functions for linear
PDEs ([19, 20, 21, 22, 23, 24, 25, 26, 27, 28]). Additionally, similar
techniques are frequently employed in electromagnetic waveguide problems to
solve the PDE directly, often though not always with the additional assumption
that $V$ is piecewise constant in the direction of propagation. See, for
example, the $S$-matrix propagation algorithm [29], the reflection-
transmission coefficient matrix (RTCM) method [30], and the transfer-matrix
method (TMM) [31], as well as the references therein. In particular, see [32]
for an example of TMM applied to the Dirac equation in a different context.
Finally, we also note that equation (1) shares some similarities with the
gravity Helmholtz equation, which arises in the study of the ionosphere, nano-
scale optical devices, seismology, and acoustics (see [33] and the references
therein).

We next describe applications of the above algorithm to compute accurate
approximations of several relevant transport properties of (1). The interface
conductivity, a physical observable describing asymmetric transport may be
defined as follows. Let $P_{x_{0}}(x,y)=P_{x_{0}}(x)$ be the Heaviside
function equal to $1$ for $x>x_{0}\in\mathbb{R}$ and equal to $0$ for
$x<x_{0}$. For notational convenience, in the following we will suppress the
dependence of $P$ on $x_{0}.$ For $\psi$ a solution of the Schrödinger
equation $i\partial_{t}\psi=H\psi$, we define
$\langle{P}\rangle=\langle{\psi|P|\psi}\rangle,$ which physically corresponds
to the probability that the particle with wavefunction $\psi$ is on the right
of the line $x=x_{0}$. The time derivative

$\displaystyle\dfrac{d{}}{d{t}}\langle{P}\rangle=\langle{\psi|i[H,P]|\psi}\rangle$

then describes the probability current crossing the line $x=x_{0}$.

Let $(E_{1},E_{2})$ be an energy interval and $\chi(h)$ a smooth non-negative
scalar function such that $\chi(h)=0$ for $h<E_{1}$ and $\chi(h)=1$ for
$h>E_{2}$. Thus, $\chi^{\prime}(h)$ may be interpreted as a density of states
confined to the energy range $(E_{1},E_{2})$. The conductivity modeling
asymmetric transport for this choice of density of states is then given by

$\displaystyle 2\pi\sigma_{I}=2\pi{\rm Tr}\ i[H,P]\chi^{\prime}(H).$ (4)

For the Dirac operator, $2\pi\sigma_{I}$ is quantized and equals $-1$ [9, 34],
indicating that the overall current integrated over the states present in the
system always equals $-1$ independent of the choice of $V$ and of the energy
interval $(E_{1},E_{2})$. See [9, 11, 3, 35, 36, 37, 38] for details and
context on the above observable describing asymmetric transport. The
theoretical computation of (4) is often greatly simplified by appealing to the
bulk-edge correspondence, relating it to the difference of topological
invariants of the two insulating phases; for derivations and applications of
such a correspondence see [3, 39, 38].

In section 6.1, we show how $\sigma_{I}$ may be related to the generalized
eigenfunctions $\psi_{\rm out}(E)$ computed in (2); see (38) below. The
explicit expression involves a number of propagating modes that depends on the
energy range $(E_{1},E_{2})$. To simplify the presentation, we will assume
$E_{2}=E_{1}$ there.

Another natural transport property of (1) is the scattering matrix generated
by the perturbation $V$. The scattering matrix describes the far field
behavior of $\psi_{\rm out}$ away from the support of $V$ associated in (2) to
any incoming plane wave $\psi_{\rm in}$. The description of the $E-$dependent
plane waves $\psi_{\rm in}$ is described in section 2.1 while scattering
matrices are defined and computed numerically in section 6.2.

This paper focuses on the Dirac model (1) with a linear (unbounded) domain
wall $m(y)=y$ so that the Green’s function given by the kernel of $(H-E)^{-1}$
is explicitly constructed. We note that a method to compute $\sigma_{I}$ for a
large class of partial differential operators was developed in [40] along with
a rate of convergence estimate. While the latter method does not require
knowledge of an unperturbed Green’s function, it is based on embedding a
domain of interest in a larger box augmented with periodic boundary
conditions. The method proposed in this paper is considerably faster
numerically as the solutions $\rho$ in (3) are computed only on the support of
$V$ and no artificial boundary conditions are necessary; see also [41] and
references there for a computation of invariants for discrete (tight-binding)
Hamiltonians.

Generalizations of the above integral formulation to bounded domain wall
profiles [40], where high-energy bulk propagation is possible, and other
partial differential models [15, 16] will be considered elsewhere. We also
note that the above transport descriptions are purely spectral. For an
analysis of localized wavepackets propagating along curved domain walls (the
domain wall $y$ above is replaced by a more general function $m(x,y)$), we
refer the reader to [42, 43, 44].

An outline for the rest of the paper is as follows. Section 2 presents the
integral equation (3) for the density $\rho$ following a spectral
decomposition of the unperturbed Dirac operator in section 2.1 and a
computation of the associated Green’s function in 2.2. Numerical
approximations of the objects appearing in the integral equation (3) are given
in section 3 and an efficient merging algorithm in section 3.2. The algorithm
to solve (3) is then described in detail in section 4. Numerical results on
accuracy and convergence are displayed in section 5. Applications to transport
calculations, including the computation of the quantized line conductivity and
of scattering matrices of illustrative choices of the perturbation $V$ are
given in section 6. Evidence for the presence of point spectrum (localized
modes) embedded in the continuous spectrum is also presented in section 6.3.
Section 7 offers some concluding remarks.

##  2 Integral formulation of the perturbed system

This section derives the integral equation (3) and presents several relevant
properties of its solutions. We first provide a spectral decomposition of the
unperturbed ($V=0$) operator in section 2.1 and construct the associated
Green’s function in section 2.2. The integral equation for $\rho$ when
$V\not=0$ and its properties are detailed in section 2.3.

To simplify calculations, we perform a unitary transformation replacing $H$ in
(1) by $QHQ^{*}$, which, with some slight abuse of notation we still call $H$,
and where $Q=\frac{1}{\sqrt{2}}(\sigma_{1}+\sigma_{3})$. We observe that
$Q^{*}=Q$ and that
$Q(\sigma_{1},\sigma_{2},\sigma_{3})Q=(\sigma_{3},-\sigma_{2},\sigma_{1}).$ In
this new basis, we thus observe that

$\displaystyle H=D_{x}\sigma_{3}-D_{y}\sigma_{2}+m(y)\sigma_{1}$ (5)

We also still denote by $V$ the transformed perturbation $QVQ^{*}$ so that
$H_{V}=H+V$ in the new basis as well.

###  2.1 Spectral decomposition of the unperturbed operator

The objective of this section is to find solutions in the kernel of $H-E$ for
$E\in\mathbb{R}$ when the perturbation $V=0$. This provides a description of
the propagating and evanescent modes we need later on and of the branches of
absolutely continuous spectrum of the operator.

The unperturbed operator $H$ in (5) is invariant under translations in the
$x-$variable. Denoting by $\xi$ the dual Fourier variable to $x$ and
$\mathcal{F}$ the corresponding transform, we observe that
$H=\mathcal{F}^{-1}\hat{H}(\xi)\mathcal{F}$ satisfies

$\displaystyle\hat{H}(\xi)-E=\xi\sigma_{3}-D_{y}\sigma_{2}+y\sigma_{1}-E.$

Let ${\mathfrak{a}}=\partial_{y}+y$ and define
$\varphi_{n}(y)=a_{n}({\mathfrak{a}}^{*})^{n}\varphi_{0}(y)$, the normalized
(in $L^{2}(\mathbb{R}_{y})$, $a_{n}$ is normalizing factor) Hermite functions
satisfying,

$\displaystyle{\mathfrak{a}}^{*}{\mathfrak{a}}\varphi_{n}=2n\varphi_{n},\quad{\mathfrak{a}}\varphi_{n}=\sqrt{2n}\varphi_{n-1},\quad{\mathfrak{a}}^{*}\varphi_{n}=\sqrt{2n+2}\varphi_{n+1},\quad\varphi_{0}(y)=\pi^{-\frac{1}{4}}e^{-\frac{1}{2}y^{2}}.$

They form an orthonormal basis of $L^{2}(\mathbb{R}_{y})$.

We define the countable set $M$ as the union of the following (pairs of)
indices $m$. For $\mathbb{N}\ni n\geq 1$, we define $m=(n,\epsilon_{m})$ with
$\epsilon_{m}=\pm 1$, while for $n=0$, we define $m=(0,-1)$. In particular,
$(0,1)\not\in M$.

While a description of the continuous spectrum of $H$ involves the branches
$\xi\to E_{m}(\xi)$, we will need a characterization of the propagating and
evanescent modes described by the ‘inverse’ maps $E\to\xi_{m}(E)$.

We thus fix energy $E\in\mathbb{R}$. The solutions to $E^{2}=\xi^{2}+2n$ are
given explicitly by

$\displaystyle\xi_{m}=\epsilon_{m}(E^{2}-2n)^{\frac{1}{2}},\qquad m\in M,$ (6)

with the choice $(E^{2}-2n)^{\frac{1}{2}}=\sqrt{E^{2}-2n}$ when $E^{2}\geq 2n$
and $(E^{2}-2n)^{\frac{1}{2}}=i\sqrt{-E^{2}+2n}$ when $E^{2}\leq 2n$.

We then define $\phi_{m}$ to be the solution of

$\displaystyle(\hat{H}(\xi_{m})-E)\phi_{m}=\begin{pmatrix}\xi_{m}-E&{\mathfrak{a}}\\\
{\mathfrak{a}}^{*}&-(\xi_{m}+E)\end{pmatrix}\phi_{m}=0.$

The normalized solutions are given by $\xi_{0}=-E$ and
$\phi_{0}=(0,\varphi_{0})^{t}$ when $m=(0,-1)$ and for $n\geq 1$ otherwise by

$\displaystyle\phi_{m}=c_{m}\begin{pmatrix}{\mathfrak{a}}\varphi_{n}\\\
(E-\xi_{m})\varphi_{n}\end{pmatrix}=c_{m}\begin{pmatrix}\sqrt{2n}\varphi_{n-1}\\\
(E-\xi_{m})\varphi_{n}\end{pmatrix},\qquad c_{m}^{-2}=2n+|E-\xi_{m}|^{2}.$ (7)

We observe that for $m=(n,\varepsilon_{m})$ and $q=(p,\varepsilon_{q})$ such
that $n\not=p$,

$\displaystyle(\phi_{m},\phi_{q})_{y}=0,$

where $(\cdot,\cdot)_{y}$ is the usual inner product on
$L^{2}(\mathbb{R}_{y};\mathbb{C}^{2})$. When $n=p$ while
$\varepsilon_{m}\not=\varepsilon_{q}$, $\phi_{m}$ and $\phi_{q}$ are linearly
independent vectors in

$\displaystyle\text{span}\Big{\\{}\varphi_{n-1}\begin{pmatrix}1\\\
0\end{pmatrix},\varphi_{n}\begin{pmatrix}0\\\ 1\end{pmatrix}\Big{\\}}.$

The functions $\phi_{m}(y;E)$ may be shown to be complete and thus form a
basis of $L^{2}(\mathbb{R}_{y};\mathbb{C}^{2})$.

In particular, we find that if $\psi_{m}$ is defined via the following formula

$\displaystyle\psi_{m}(x,y;E)=e^{i\xi_{m}x}\phi_{m}(y;E)$ (8)

then $(H-E)\psi_{m}(x,y;E)=0$. Moreover, when $\xi_{m}$ is real-valued (i.e.,
when $E^{2}\geq 2n$), $\psi_{m}(x,y;E)$ is a (bounded) generalized eigenvector
(plane wave) of $H$ associated with energy $E$. These generalized eigenvectors
form the set of incoming plane waves $\psi_{\rm in}$ mentioned in the
introduction, section 1. We refer to them as propagating modes.

When $\xi_{m}$ is purely imaginary, the corresponding $\psi_{m}(x,y;E)$ is an
evanescent mode. Accounting for evanescent modes is necessary in the presence
of perturbations $V$ and in the computation of the Green’s function described
in the next section.

The above calculations described propagating and evanescent modes depending on
the reality of $\xi_{m}$. The same calculations provide a complete spectral
decomposition of $H$ as follows. Define

$\displaystyle
E_{m}(\xi)=\epsilon_{m}\sqrt{\xi^{2}+2n},\qquad\phi_{m}=c_{m}\begin{pmatrix}{\mathfrak{a}}\varphi_{n}\\\
(E_{m}-\xi)\varphi_{n}\end{pmatrix},\quad c_{m}^{-2}=2n+|E_{m}-\xi|^{2}.$ (9)

Each $\xi\to E_{m}(\xi)$ is a branch of absolutely continuous spectrum with
generalized eigenvectors solving $(H-E_{m}(\xi))\psi_{m}=0$ and given by
$\psi_{m}(x,y;\xi)=e^{i\xi x}\phi_{m}(y;\xi)$.

Thus, $\sigma(H)=\sigma_{\rm ac}(H)=\mathbb{R},$ i.e., the whole real line is
absolutely continuous spectrum, with the following energy-dependent structure;
see also Fig.1.

In the range $E^{2}<2$, we obtain a unique branch of absolutely continuous
spectrum $m=(0,-1)$ parameterized by a simple branch of generalized
eigenvectors $\psi_{(0,-1)}(x,y;\xi)$ with group velocity equal to $\partial
E(\xi)/\partial\xi=-1$ (straight, blue curve in Fig.1). In this energy range,
no mode displays a positive group velocity. A perturbation $V$ of $H$ can only
generate a phase shift that does not perturb the particle density. Back-
scattering is therefore entirely suppressed in this energy range, both for
topological and energetic reasons.

The spectrum is symmetric in $E\to-E$ and it therefore suffices to assume
$E\geq 0$. For $2<E^{2}<4$, the absolutely continuous spectrum has a
degeneracy equal to $3$ corresponding to three propagating modes $\psi_{\rm
in}^{j}$ for $1\leq j\leq 3$ and generating corresponding outgoing modes
$\psi_{\rm out}^{j}$ solution of (2) (see level set $E=1.8$ in Fig.1). When
$V=0$, two modes propagate towards negative values of $x$ while one propagates
in the opposite direction. In terms of asymmetric transport,
$-2+1=2\pi\sigma_{I}$ is valid. When $V\not=0$, the three propagating modes
interact via the perturbation as we will demonstrate in section 6.2. Back-
scattering is then clearly present even though the ‘topology’ of the Dirac
operator is still non-trivial. We will present numerical simulations showing
that the structure of the scattering matrix strongly depends on the choice of
perturbation $V$ while $2\pi\sigma_{I}=-1$ holds up to arbitrary (i.e., at
least $14$ digits here) accuracy even in the presence of strong perturbations
$V$.

Figure 1: Three branches of absolutely continuous spectrum of the unperturbed
operator $H$ with dots describing the propagating modes for different fixed
energy levels $E$. While no backscattering is present for $|E|<\sqrt{2}$, this
is no longer true for $|E|>\sqrt{2}$.

More generally, the degeneracy of the spectrum is equal to $2n+1$ for energies
in the range $2n<E^{2}<2(n+1)$. The perturbation $V$ then generates
$(2n+1)\times(2n+1)$ (far field) scattering matrices; see level set $E=2.2$ in
Fig.1 when $n=2$. We illustrate such computations for $n=0,1,2$ in section
6.2.

###  2.2 Green’s Function of the unperturbed system

We now construct the outgoing Green’s function of the operator $(H-E)$ for
$E\in\mathbb{R}$ (i.e., the Green’s function of $H-(E-i0^{+})$) when the
perturbation $V=0$. We will in fact avoid a set of (Lebesgue) measure zero and
assume that $E^{2}\not=2n$ for $n\in\mathbb{N}_{*}$.

Consider the Green’s function $G$ which is the solution of the following
equation

$\displaystyle(H-E)G=\delta(x-x_{0})\delta(y-y_{0})I.$

Without loss of generality, since $H$ is translationally invariant in $x$, we
can assume that $x_{0}=0.$ Then,

$\displaystyle(H+E)(H-E)G$
$\displaystyle=(H^{2}-E^{2})G=(H+E)\delta(x)\delta(y-y_{0})I,$
$\displaystyle\implies G$
$\displaystyle=(H+E)(H^{2}-E^{2})^{-1}\delta(x)\delta(y-y_{0})I,$

where we have used the fact that clearly $(H-E)$ and $(H+E)$ commute.
Additionally, we note that here

$\displaystyle
H^{2}-E^{2}=\begin{pmatrix}D_{x}^{2}-E^{2}+{\mathfrak{a}}{\mathfrak{a}}^{*}&0\\\
0&D_{x}^{2}-E^{2}+{\mathfrak{a}}^{*}{\mathfrak{a}}\end{pmatrix}$

is a diagonal matrix of scalar operators. We thus need to solve two similar
equations

$\displaystyle(-\partial_{x}^{2}-\partial_{y}^{2}+y^{2}\pm
1-E^{2})G_{\pm}=\delta(x)\delta(y-y_{0}),$

and then observe that the complete Green’s function is given by

$\displaystyle G=(H+E){\rm Diag}(G_{+},G_{-}).$

We recall that ${\mathfrak{a}}^{*}{\mathfrak{a}}\varphi_{n}=2n\varphi_{n}$ and
${\mathfrak{a}}{\mathfrak{a}}^{*}\varphi_{n}=(2n+2)\varphi_{n}$ and assume
that $E^{2}\not=2n$ for each $n$.

Expanding $G_{-}$ in the basis $\phi_{n}$, and applying the operator
$D_{x}^{2}-E^{2}+2n$, we see that if

$\displaystyle G_{-}=\displaystyle\sum_{n}G_{-,n}(x)\varphi_{n}(y),\quad{\rm
then}\quad(D_{x}^{2}-E^{2}+2n)G_{-,n}(x)=\delta(x)\varphi_{n}(y_{0}).$

When $2n>E^{2}$, we obtain from the jump condition at $x=0$ the evanescent
modes

$\displaystyle
G_{-,n}(x)=\frac{1}{2\sqrt{2n-E^{2}}}e^{-\sqrt{2n-E^{2}}|x|}\varphi_{n}(y_{0}),\qquad
2n>E^{2}.$

When $2n<E^{2}$, two propagating modes satisfy the above scalar equation for
$G_{-}$. Keeping only the outgoing mode since we are interested in computing
the outgoing Green’s function, (i.e., keeping $e^{ikx}$ for $k$ and $x$
positive while $e^{-ikx}$ would be considered incoming), we thus observe that

$\displaystyle
G_{-,n}(x)=\frac{1}{-2i\sqrt{E^{2}-2n}}e^{i\sqrt{E^{2}-2n}|x|}\varphi_{n}(y_{0}),\qquad
2n<E^{2}.$ (10)

To combine the above calculations into a single expression, we define

$\displaystyle\theta_{n}=i\sqrt{E^{2}-2n}=i\xi_{(n,1)}.$ (11)

The complex numbers $\theta_{n}$ are thus related to $\xi_{m}$ introduced
earlier in the analysis of the Dirac operator.

We thus find

$\displaystyle G_{-}(x,y;y_{0})=\displaystyle\sum_{n\geq
0}\frac{-1}{2\theta_{n}}e^{\theta_{n}|x|}\varphi_{n}(y)\varphi_{n}(y_{0}).$
(12)

For $G_{+}$, we have the same expression except that $2n$ is replaced by
$2n+2$:

$\displaystyle G_{+}(x,y;y_{0})=\displaystyle\sum_{n\geq
0}\frac{-1}{2\theta_{n+1}}e^{\theta_{n+1}|x|}\varphi_{n}(y)\varphi_{n}(y_{0}).$
(13)

It remains to apply $H+E$ to get the final result

$\displaystyle G=\begin{pmatrix}(D_{x}+E)G_{+}&{\mathfrak{a}}G_{-}\\\
{\mathfrak{a}}^{*}G_{+}&(-D_{x}+E)G_{-}\end{pmatrix}.$ (14)

Note that $G$ is a $2\times 2$ matrix-valued function. We also use the
notation $G(x,y;x_{0},y_{0})=G(x-x_{0},y;y_{0})$. See also (22) below for an
explicit expression of each entry of the matrix.

#### Numerical Green’s Function.

As a numerical illustration, consider the Green’s function
$G(x,y,x_{0},y_{0})$ with $E=1.8$. The source is located at
$(x_{0},y_{0})=(0,1)$. In Fig.2, we plot contours of absolute values of the
four entries of $G$. For the first plots, we truncate the series expansion for
$G$ at $n_{y}=10^{5},$ and then show how the $|x-x_{0}|^{-1}$ singularity near
the source [45] is captured as $n_{y}$ varies in Fig.3. We observe in Fig.3
that the Green’s function is asymmetric in the $x$ direction, which is a
manifestation of the transport asymmetry of the unperturbed operator. In the
$y$ direction, which is confined by the domain wall, the Green’s function
decreases exponentially as $|y|$ increases.

(a) $|G_{11}|$

(b) $|G_{12}|$

(c) $|G_{21}|$

(d) $|G_{22}|$

Figure 2:  Absolute value (truncated to $[0,1]$) of the components of the
Green’s function with $x_{0}=0$, $y_{0}=1$, and $n_{y}=100000$.

In Fig.3, we investigate the local behaviour of the Green’s function. It is
known that $G$ has a $\frac{1}{r}$ singularity. In (14), the spectral series
expansions of $G_{+}$, see (13), and $G_{-}$, see (12), converge slowly and
have oscillations when $n$ is small near the source.

(a) in $x$ direction, $y=y_{0}+10^{-2}$

(b) in $y$ direction, $x=x_{0}+10^{-2}$

Figure 3:  Local behaviour of approximations of $|G_{11}|$ near
$(x_{0},y_{0})=(0,1)$ for several values of $n_{y}$.

###### Remark 2.1

The convergence of the spectral decomposition of $G$ is necessarily slow since
$G$ is singular. However, in the algorithm described below to compute
generalized eigenfunctions and densities $\rho$, only the integration of $G$
against smooth functions is involved.

###  2.3 Generalized eigenfunctions and integral formulation

We now use the above unperturbed Green’s function to solve for generalized
eigenfunctions of the perturbed operator $H_{V}=H+V$. Let $\psi_{\rm in}$ be
an incoming plane wave, e.g. $\psi_{m}$ as described in (8) with $E^{2}>2n$ so
that $\xi_{m}$ is real-valued. We then look for solutions of $(H+V-E)\psi=0$
of the form

$\displaystyle\psi=\psi_{\rm in}+\psi_{\rm out},$

so that

$\displaystyle(H+V-E)\psi_{\rm out}=-V\psi_{\rm in},$ (15)

where we assume $\psi_{\rm out}$ to be outgoing. In other words, using the
unperturbed outgoing Green’s function constructed in the preceding section, we
look for

$\displaystyle\psi_{\rm out}(x,y)=\int
G(x,y;x_{0},y_{0})\rho(x_{0},y_{0})dx_{0}dy_{0},$ (16)

where $\rho(x,y)$ is an unknown density. Since $G$ is the kernel of the
integral operator $(H-E)^{-1}$ (with outgoing conditions), we readily obtain
that $\rho$ satisfies the integral equation

$\displaystyle\rho(x,y)+V(x,y)\int
G(x,y;x_{0},y_{0})\rho(x_{0},y_{0})dx_{0}dy_{0}=-V(x,y)\psi_{\rm in}(x,y).$
(17)

In what follows, we assume that $V$ is a smooth Hermitian matrix-valued
compactly supported function. As is apparent from the above formulation, the
density $\rho$ vanishes outside of the support of $V$. The above formulation
is therefore convenient numerically as only the support of $V$ needs to be
discretized in practice.

###### Remark 2.2

For ease of exposition, here we only consider the case in which the right-hand
side of (17) is generated by an incoming propagating mode. It is relatively
straightforward to modify the approach to include more general (compactly-
supported) right-hand sides.

#### Properties of the density $\rho(x,y)$.

The above equation for $\rho$ can be re-written in the form

$\displaystyle\rho+V\mathcal{G}\rho=-V\psi_{\rm in}$ (18)

with $\mathcal{G}=(H-E)^{-1}$ the linear operator with kernel
$G(x,y;x_{0},y_{0})$. Since the Dirac operator is elliptic, the inverse
$\mathcal{G}$ maps $H^{k}(\mathbb{R}^{2};\mathbb{C}^{2})$ to
$H^{k+1}(\mathbb{R}^{2};\mathbb{C}^{2})$, where
$H^{k}(\mathbb{R}^{2};\mathbb{C}^{2})$ is the space of $\mathbb{C}^{2}$ valued
functions with square integrable derivatives of order up to $k$. Since $V$ is
smooth and compactly supported, we obtain that $V\mathcal{G}$ is a compact
operator from $H^{k}(\mathbb{R}^{2};\mathbb{C}^{2})$ to itself for any $k\geq
0$.

The above equation thus admits a unique solution provided $-1$ is not an
eigenvalue of $V\mathcal{G}$ (in any $H^{k}(\mathbb{R}^{2};\mathbb{C}^{2})$).
The complement of the latter condition is then equivalent to the existence of
a normalized (in $L^{2}(\mathbb{R}^{2};\mathbb{C}^{2})$) solution of the
problem

$\displaystyle(H+V-E)\psi=0.$ (19)

While solutions of (19) for Schrödinger and Dirac (with constant mass term)
operators $H$ are known not to exist for $E$ in the continuous spectrum of $H$
when $V$ is compactly supported, see for instance [46] and [47, Theorem
XIII.58], the confinement afforded by the domain may in fact generate such
point spectrum. We revisit the question in section 6.3, where we show the
presence of localized modes for certain choices of $V$ which are compactly
supported in $x$.

When $-1$ is not in the spectrum of the compact operator $V\mathcal{G}$, we
apply the Fredholm alternative to (18) and formally define

$\displaystyle\rho=-(I+V\mathcal{G})^{-1}V\psi_{\rm in}.$ (20)

Since $V\psi_{\rm in}$ is smooth and supported on the support of $V$, we
obtain that $\rho$ is itself smooth and compactly supported and therefore an
element in $C^{\infty}_{c}(\mathbb{R}^{2};\mathbb{C}^{2})$.

We now show that $-1$ may be in the spectrum of $V\mathcal{G}$ only for
discrete values of $E$ when $V$ is compactly supported (exponential decay
being sufficient). Indeed, following the proof in, e.g., [48, Theorem XI.45],
$E\to V\mathcal{G}$ is analytic in the vicinity $\mathcal{V}$ of any interval
$[E_{1},E_{2}]\subset\mathbb{R}$ such that $\\{2n;\
n\in\mathbb{N}\\}\cap[E_{1}^{2},E_{2}^{2}]=\emptyset$. Since $H+V$ is self-
adjoint, $-1$ cannot be an eigenvalue of $V\mathcal{G}$ as soon as
$E\in\mathcal{V}$ is not purely real. The analytic Fredholm theory [49,
Theorem VI.14] then states that $-1$ is an eigenvalue of $V\mathcal{G}$ only
for discrete (real) values of $E$. For the rest of the paper, we assume that
$E$ does not belong to any of these discrete sets or to $\\{\sqrt{2n};\
n\in\mathbb{N}\\}$.

#### Basis expansion of the density $\rho$.

In the following derivation, we denote

$\displaystyle\rho(x,y)=\sum_{i,n}\rho_{i,n}P_{i}(x)\varphi_{n}(y),$ (21)

where $P_{i}$ denotes the i-th Legendre polynomial and $\varphi_{n}$ denotes
$n$-th order Hermite functions, while $\rho_{i,n}\in\mathbb{C}^{2}$ are the
corresponding generalized Fourier coefficients. The generalized eigenfunction
solution of $(H+V-E)\psi=0$ is then given by $\psi=\psi_{\rm
in}+\mathcal{G}\rho$.

From the expression (14) for the Dirac Green’s function, we observe that the
integral operator in (17) involves explicit integrals given by

$\displaystyle\int G(x,y;x_{0},y_{0})\begin{pmatrix}\varphi_{n}(y_{0})\\\
0\end{pmatrix}dy_{0}=\begin{pmatrix}\frac{1}{2}(i\,{\rm
sign}(x-x_{0})-\frac{E}{\theta_{n+1}})e^{\theta_{n+1}|x-x_{0}|}\varphi_{n}(y)\\\
-\frac{\sqrt{2n+2}}{2\theta_{n+1}}e^{\theta_{n+1}|x-x_{0}|}\varphi_{n+1}(y)\end{pmatrix}$
(22) $\displaystyle\int G(x,y;x_{0},y_{0})\begin{pmatrix}0\\\
\varphi_{n}(y_{0})\end{pmatrix}dy_{0}=\begin{pmatrix}-\frac{\sqrt{2n}}{2\theta_{n}}e^{\theta_{n}|x-x_{0}|}\varphi_{n-1}(y)\\\
\frac{1}{2}(-i\,{\rm
sign}(x-x_{0})-\frac{E}{\theta_{n}})e^{\theta_{n}|x-x_{0}|}\varphi_{n}(y)\end{pmatrix}.$
(23)

We may then compute $\psi_{\rm out}$ from $\rho$ using (16). Outside of the
support of $\rho$, we observe that the approximation of $\psi_{\rm out}$
decays exponentially as $|y|$ increases and is a sum of propagating modes that
oscillate in the $x-$variable along the interface $y\equiv 0$ and of
evanescent modes that decay exponentially away from the support of $\rho$.

##  3 Numerical preliminaries

###  3.1 Global discretization

In order to obtain a numerical approximation of the density $\rho,$ we
truncate the expansion of $\rho$ in (21), keeping a finite number of
coefficients. In the following, we denote the solution of this truncated
system by $\hat{\rho}.$ Additionally, we denote the finite product space
consisting of products of linear combinations of the first $n_{x}-1$ Legendre
polynomials in $x$ with linear combinations of the first $n_{y}$ Hermite
functions in $y$ by $\mathcal{V}_{n_{x},n_{y}}.$ In particular,
$\mathcal{V}_{n_{x},n_{y}}=\mathcal{P}_{n_{x}}\otimes(\Phi_{n_{y}}\oplus\Phi_{n_{y}})$,
where $\mathcal{P}_{n_{x}}$ denotes the span of the first $n_{x}$ Legendre
polynomials in $x$, and $\Phi_{n_{y}}$ denotes the span of the first $n_{y}$
Hermite functions in $y.$ In the sequel we use $\mathcal{V}$ for short.
Empirically, we observe that in this basis the solution of the projected
system converges spectrally (super-algebraically) in $n_{x}$ and $n_{y}$ to
the solution of the original system, given sufficient regularity of the
potential $V.$ Interested readers can refer to [50] for applicable theorems on
convergence.

Clearly, there are two linear operators in (17) which must be truncated and
whose action must be computed numerically: integration with $G$ and
multiplication of $V$.

  * •

Kernel integration with $G$. The integration with respect to $y_{0}$ is
analytically available, given by equations (22) and (23). In the $x$
direction, the integrands have a singularity in their derivatives when
$x_{0}=x$. Thus, using the $n_{x}$ point quadrature rule for Legendre
polynomials (i.e. the $n_{x}$ Gauss-Legendre quadrature rule) leads to slow
convergence of the integral over the support of $V$ in $x$, which we denote by
$[x_{L},x_{R}].$ To circumvent this issue, for any fixed $x\in[x_{L},x_{R}]$,
we divide the interval $[x_{l},x_{R}]$ into two subintervals $[x_{L},x]$ and
$[x,x_{R}]$. On both intervals, the integrand is smooth and we can apply
standard smooth quadrature rules to compute the integral of the Green’s
function (after integrating over $y_{0}$ analytically) against the basis
functions (i.e. suitably rescaled, shifted, and dilated Legendre polynomials).
Furthermore, we observe that by choosing the values of $x$ at which to
evaluate the integrals as the nodes of the $n_{x}$-point Gauss-Legendre
quadrature rule, the map from values at these nodes to the coefficients in an
$n_{x}$-th order Legendre expansion is invertible and stable. We let $\hat{G}$
denote the $2n_{x}n_{y}\times 2n_{x}n_{y}$ square matrix obtained by
discretizing the application of the Green’s function in this way.

  * •

Multiplication by $V$. We recall that by assumption, $V$ is smooth and
compactly supported in both the $x$ and $y$ directions. Evidently, pointwise
multiplication by $V$ is a linear operator, and in the orthonormal basis
$\mathcal{V},$ has matrix entries given by
$(\mathcal{V}_{i},V\mathcal{V}_{j})_{x,y}$, where $\mathcal{V}_{i}$ and
$\mathcal{V}_{j}$ are two basis functions from $\mathcal{V}$. To compute it,
we first evaluate $\mathcal{V}_{j}$ on a tensor product of Legendre and
Hermite quadrature nodes. The orders of these quadrature rules should be large
than $n_{x}$ and $n_{y}$, respectively. Indeed, using standard estimates from
approximation theory it is possible to derive explicit upper bounds for the
oversampling required as a function of the smoothness and support of $V.$ In
practice, we observe that for most of the examples considered here, increasing
the integration orders by no more than 20 suffices, though obviously this
would need to be increased for rougher $V$’s. One then computes the product of
$\mathcal{V}_{i},$ $\mathcal{V}_{j},$ and $V,$ at these nodes and sums against
the quadrature weights. This process generates a $2n_{x}n_{y}\times
2n_{x}n_{y}$ square matrix $\hat{V}$ which encodes this transformation.

In addition to the previously defined $\hat{G}$ and $\hat{V}$, we denote
$\Pi_{\mathcal{V}}$ as the quadrature projection to the product basis
$\mathcal{V}$. Using these matrices, the discretization of our integral
formulation is given by

$\displaystyle(I+\hat{V}\hat{G})\hat{\rho}=-\hat{V}\Pi_{\mathcal{V}}\psi_{\rm
in}.$

The last step in computing the generalized eigenfunctions is to recover
$\psi_{\rm out}$ from the integral representation (16). The eigenfunction
$\psi$ is given by

$\displaystyle\hat{\psi}(x,y)=\psi_{\rm in}(x,y)+\int
G(x,y;x_{0},y_{0})\hat{\rho}(x_{0},y_{0})dx_{0}dy_{0}.$ (24)

It is worth emphasizing that after the computation of $\rho$ we are not
restricted to evaluate $\psi_{out}$ only within the support of $V.$ Indeed,
again using the integral representation (16), $\psi_{\rm out}$ can be computed
directly by numerically evaluating $\int
G(x,y;x_{0},y_{0})\rho(x_{0},y_{0})dx_{0}dy_{0}$ at any point $(x,y)$. The
discretization of this integral is similar to the one employed in the
construction of $\hat{G}$, i.e., in the $x$ direction, we can split the domain
of integration, if necessary, at $x_{0}=x,$ while in the $y$ direction, (22),
(23) give analytic expressions for the integrals.

###  3.2 Computing eigenfunctions by merging TR matrices

In this section, we use the transmission reflection (TR) matrix formalism (see
[30] and the references therein, for example) to solve the integral equation
((17)) on large domains.

We begin by defining the TR matrices. To that end, first suppose that $V$ is
compactly supported in the x direction on the interval $[x_{L},x_{R}]$ and
$\psi$ is the eigenfunction of $H-E+V$ given some incoming condition. We
decompose $\psi$ as

$\displaystyle\psi(x,y)=\displaystyle\sum_{m\in M}\alpha_{m}(x)\phi_{m}(y),$
(25)

with $\alpha_{m}(x)=(\phi_{m}(y),\psi(x,y))_{y}$ proportional to
$e^{i\xi_{m}x}$ outside of the support of $V$ and $\phi_{m}$ defined in (7).
In this way, instead of mapping the incoming field to the outgoing field in
the far field, we consider the linear map right at the boundary of
$[x_{L},x_{R}]$. The incoming condition is then the coefficients of the right
traveling modes $\phi_{m}$ ($\epsilon_{m}>0$) at left boundary and the left
traveling modes $\phi_{m}$ ($\epsilon_{m}<0$) at the right boundary, namely
$\alpha_{+}(x_{L})$ and $\alpha_{-}(x_{R})$. The outgoing solution consists of
the coefficients of the left traveling modes $\phi_{m}$,($\epsilon_{m}<0$), at
left boundary and the right traveling modes $\phi_{m}$ ($\epsilon_{m}>0$) at
right boundary, namely $\alpha_{-}(x_{L})$ and $\alpha_{+}(x_{R})$. In
summary, the TR matrix for $V$ restricted to the interval $[x_{L},x_{R}]$,
denoted by $M$, can be defined as,

$\displaystyle\begin{pmatrix}\alpha_{-}(x_{L})\\\
\alpha_{+}(x_{R})\end{pmatrix}=\begin{pmatrix}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{pmatrix}\begin{pmatrix}\alpha_{-}(x_{R})\\\
\alpha_{+}(x_{L})\end{pmatrix}.$ (26)

In numerical computation, if we restrict the density $\rho$ to $\Phi_{n_{y}}$
in the $y$ direction, in which case the TR matrix is a
$(2n_{y}-1)\times(2n_{y}-1)$ matrix. $M_{11}$ ($M_{22}$) denotes the left
(right) transmission coefficients, and $M_{21}$ ($M_{12}$) denotes the left
(right) reflection coefficients.

#### Merging Two TR Matrices

It is of particular interest in physics to consider the case when there is an
elongated perturbation near the interface. The TR matrix in this case will
depend on the length of interface, $l$, in the $x$ direction. Obtaining the TR
matrix requires solving for the total field $\psi$ for each of the $2n_{y}-1$
independent modes as incoming conditions. When solving $\psi$, if the support
of $V$ in the $x$ direction is large, accurate computation of $\rho$ requires
very high order quadrature rules in the $x$ direction in order to resolve the
local structure of $V$ and $\rho$ in a relatively long interval. Clearly, a
naive implementation then leads to a significant computational cost (both in
time and memory), as well as possible numerical instabilities. To avoid this,
we break the problem up into small intervals or leaves in the $x$ direction
and compute the TR matrices of these sub-problems. It is then possible to
merge two TR matrices for adjacent intervals to obtain a TR matrix for their
union. We can continue merging in this way until the TR matrix for the entire
domain is obtained.

To that end, let $L$ and $R$ denote the TR matrices of two adjacent interval
in the $x$ direction, $I_{L}$ and $I_{R}$, respectively. Then the TR matrix
for $I_{L}\cup I_{R}$ is

$\displaystyle\begin{pmatrix}L_{11}(I-R_{12}L_{21})^{-1}R_{11}&L_{11}(I-R_{12}L_{21})^{-1}R_{12}L_{22}+L_{12}\\\
R_{22}(I-L_{21}R_{12})^{-1}L_{21}R_{11}+R_{21}&R_{22}(I-L_{21}R_{12})^{-1}L_{22}\end{pmatrix}$
(27)

Furthermore, given incoming waves at the left and right boundary, we can also
calculate the Fourier coefficient at the intersection of two merging
intervals. Specifically, the intersection coefficient matrix $\mathcal{M}$ is
defined by,

$\displaystyle\begin{pmatrix}\alpha_{-,\mathcal{M}}\\\
\alpha_{+,\mathcal{M}}\end{pmatrix}=\begin{pmatrix}\mathcal{M}_{11}&\mathcal{M}_{12}\\\
\mathcal{M}_{21}&\mathcal{M}_{22}\end{pmatrix}\begin{pmatrix}\alpha_{-,R}\\\
\alpha_{+,L}\end{pmatrix},$ (28)

where $\begin{pmatrix}\alpha_{-,\mathcal{M}}\\\
\alpha_{+,\mathcal{M}}\end{pmatrix}$ denotes the Fourier coefficient projected
to unperturbed eigenfunction basis $\\{\phi_{m}\\}$ at the intersection.
$\mathcal{M}$ can be directly calculate by,

$\displaystyle\mathcal{M}=\begin{pmatrix}(I-R_{12}L_{21})^{-1}R_{11}&(I-R_{12}L_{21})^{-1}R_{12}L_{22}\\\
(I-L_{21}R_{12})^{-1}L_{21}R_{11}&(I-L_{21}R_{12})^{-1}L_{22}\end{pmatrix}.$
(29)

#### Binary merging and recovering the eigenfunctions inside the material

Merging TR matriices is independent of interval position and length, so for a
general partition of a single interval we can compute the TR matrix first for
each leaf by Alg.1 and merge them. Here we apply a binary merging strategy. To
be more precise, we first partition the whole interval into $2^{L}$ leaves
$\\{I_{i}\\}_{i=1:2^{L}}$ of equal length. Here $L$ denotes the total number
of levels. We compute the TR matrices on each interval $I_{i}$ and merge two
adjacent intervals $I_{2k-1}$ and $I_{2k}$ where $k=1\cdots 2^{L-1}$. We
continue this merging $L$ times (i.e. at each level) until we have the TR
matrix for the whole interval. When merging two adjacent intervals, the TR
matrix reflects the outgoing waves as a linear transform of incoming waves.
The intersection coefficient matrix then reveals the Fourier coefficient
($\alpha$) at intersection is a linear function of incoming waves at the
boundaries. In this way, after merging all the intervals, we can recover
$\alpha$ from the last level of merging back to the first level. The full
algorithm is summarized in Alg.2. An illustration of the algorithm is given by
Fig.4, where, for concreteness, we compute the TR matrix of $[0,8]$.

##  4 The algorithm

In this section we summarize our algorithms for computing generalized
eigenfunctions for compact smooth perturbations of Dirac equations with linear
domain walls. The first algorithm (Alg.1) computes the density and
eigenfunction in a single slab without merging. As input, it takes in a
perturbation $V$ (a compactly supported smooth, $2\times 2$ Hermitian matrix-
valued function), an interval $[x_{l},x_{R}]$ containing the support of $V$ in
the $x$-direction, an integer $n_{y}$ giving the order of the Hermite
expansion to be used, an integer $n_{x}$ giving the order of the Legendre
expansion to be used It also requires as input two vectors
$(\alpha_{+},\alpha_{-1})\in\mathbb{C}^{n_{y}-1}\times\mathbb{C}^{n_{y}}$
corresponding to the amplitudes of the incoming waves. As output it returns a
tensor-product (piecewise) Legendre - Hermite approximation to the density
$\rho$ and outgoing amplitudes $\alpha_{-}(x_{L})$ and $\alpha_{+}(x_{R})$.
Algorithm 1 Computing density and eigenfunction in a single slab (leaf)
1:Potential Field $V$; Interval of $V$ that is compacted supported
$I=[x_{L},x_{R}]$; Level of binary merging $L$; Incoming wave condition
$\alpha_{+}(x_{L})$ and $\alpha_{-}(x_{R})$; discretization configuration
$(n_{x},n_{y})$.  2:Density $\hat{\rho}=\hat{\rho}^{n_{x},n_{y}}$ and
eigenfunction $\hat{\psi}$ given incoming condition; outgoing wave
$\alpha_{-}(x_{L})$ and $\alpha_{+}(x_{R})$.  3:Construct $\hat{V}$ and
$\hat{G}$ projected to $n_{x}$ Legendre polynomials and $n_{y}$ Hermite
functions.  4:Compute $\Pi_{\mathcal{V}}\psi_{\rm in}$ with
$\alpha_{-}(x_{L})$ and $\alpha_{+}(x_{R})$.  5:Solve $\hat{\rho}$ by
$\hat{\rho}=-(I+\hat{V}\hat{G})^{-1}\hat{V}\Pi_{\mathcal{V}}\psi_{\rm
in}$.111$(I+\hat{V}\hat{G})^{-1}$ can be re-used when constructing TR matrix.
6: * optional step  7:Recover $\psi$ by $\hat{\psi}(x,y)=\psi_{\rm
in}(x,y)+\int G(x,y;x_{0},y_{0})\hat{\rho}(x_{0},y_{0})dx_{0}dy_{0}$.*
8:Extract $\alpha_{-}(x_{L})$ by
$\alpha_{n,-}(x_{L})=\int\phi_{n,-}(y)G(x_{L},y;x_{0},y_{0})\hat{\rho}(x_{0},y_{0})dx_{0}dy_{0}dy$.*
9:Extract $\alpha_{+}(x_{R})$ by
$\alpha_{n,+}(x_{R})=\int\phi_{n,+}(y)G(x_{R},y;x_{0},y_{0})\hat{\rho}(x_{0},y_{0})dx_{0}dy_{0}dy$.*

The second algorithm (Alg.2) describe the hierarchical merging that compute
generalized eigenfunctions efficiently. In addition to inputs that are
required by Alg.1, we used $L$ to specified total level of merging.

Algorithm 2 Computing eigenfunctions in a slab with potential $V$

1:Potential Field $V$; Interval of $V$ that is compacted supported
$I=[x_{L},x_{R}]$; Level of binary merging $L$; Incoming wave condition
$\alpha_{+}(x_{L})$ and $\alpha_{-}(x_{R})$.

2:TR matrix of interval $I$, $M_{I}$; eigenfunctions given incoming wave,
$\alpha$.

3:Partition $I$ to be $2^{L}$ intervals. Denote the intervals as
$\\{I_{k}\\}$, $k=1,2,\cdots,2^{L}$ and the grid points as $x_{k}$,
$k=0,1,\cdots,2^{L}$.

4:for $k$ in $1,\cdots,2^{L}$ do

5: Compute TR matrix $M_{0,k}$ with potential field $V$ limited on interval
$I_{k}$ by Alg.1.

6:end for

7:for $l$ in $1\to L$ do

8: for $k$ in $1\to 2^{(L-l)}$ do

9: Merge $M_{l,2k-1}$ with $M_{l,2k}$:

10: Calculate the TR matrix $M_{l+1,k}$ for $[x_{(k-1)2^{l}},x_{k2^{l}}]$ from
(27)

11: Calculate intersection coefficient matrix $\mathcal{M}_{l+1,k}$ at
$x_{k2^{l}-2^{l-1}}$ from (29).

12: end for

13:end for$M_{I}=M_{L+1,1}$

14:Assign $\alpha_{+}(x_{0})=\alpha_{+}(x_{L})$,
$\alpha_{-}(x_{0})=0$,$\alpha_{+}(x_{2^{L}})=0$,$\alpha_{-}(x_{2^{L}})=\alpha_{-}(x_{R})$.

15:for $l$ in $L\to 1$ do

16: for $k$ in $1,\cdots,2^{(L-l)}$ do

17: $\begin{pmatrix}\alpha_{-}(x_{(2k-1)2^{l-1}})\\\
\alpha_{+}(x_{(2k-1)2^{l-1}})\end{pmatrix}=\mathcal{M}_{l+1,k}\begin{pmatrix}\alpha_{-}(x_{k2^{l}})\\\
\alpha_{+}(x_{(k-1)2^{l}})\end{pmatrix}$

18: end for

19:end for

20:for $k$ in $1,\cdots,2^{L}$ do

21: Recover $\rho|_{[x_{k-1},x_{k}]}$ by $\alpha_{+}(x_{k-1})$ and
$\alpha_{-}(x_{k})$ via Alg.1.222optional step, we can then compute values of
$\psi_{\rm out}$ at any point by (16).

22:end for

Figure 4: An Illustration of Merging/Recovering Algorithm

###### Remark 4.1

For ease of exposition, here we have restricted our attention to uniform leaf
size. The method described above generalizes naturally to the non-uniform
case, in which the leaf size is chosen to depend on the local smoothness of
$V.$ Moreover, it can also be easily extended to the fully adaptive case in
which, after a coarse solve, leaves are refined adaptively until the solution
$\rho$ is fully resolved. See [51] for a description of this approach in the
context of ODE solvers.

##  5 Numerical results

#### Numerical solution in one interval.

As a first example, we consider a single interval with potential
$V=V_{1}\chi_{I}$ where

$\displaystyle
V_{1}=\exp\left(-y^{2}\right)y\cos\left(\left(\xi_{(0,-1)}-\xi_{(1,1)}\right)x\right),$
(30)

and it is supported on $I=[-1,1]$. For the incoming field we choose
$\psi_{in}=e^{i\xi_{(0,-1)}x}\phi_{(0,-1)}(y)$. It is a left traveling plane
wave in $x$ direction. We compute $\rho$ via Alg.1 with $n_{x}=20$ and
$n_{y}=108$ (we will justify the accuracy in the following paragraph). In
Fig.5, we recover the generalized eigenfunction obtained via the formula (24).

Looking at Fig.5, we note that our integral formulation preserves the plane
wave structure outside the support of $V$, i.e. the interval $I$. Moreover, in
comparing the plot of $\psi_{1}$ with $\psi_{2}$, one can clearly observe the
difference in wavelength. This is expected, since the $(0,-1)$ mode vanishes
in the first entry of $\psi$ and hence the dominant wavelength in $\psi_{2}$
is $\frac{2\pi}{\sqrt{E^{2}-2}}\approx 5.64$ and the dominant wavelength in
$\psi_{1}$ is $\frac{2\pi}{E}\approx 3.49$.

(a) Real part of $\psi_{1}$

(b) Imaginary part of $\psi_{1}$

(c) Real part of $\psi_{2}$

(d) Imaginary part of $\psi_{2}$

(e) Amplitude of Fourier coefficients, $|\alpha|$; Slender lines denotes
evanescent modes.

Figure 5: Solution $\psi=\psi_{in}+\psi_{out}$

#### Verification of spectral convergence.

Next, we discuss the results of a self-convergence study to verify the
spectral convergence in computing $\rho$ by Alg.1. In Fig.6 we show the error
between the reference solution and the solutions with different discretization
orders ( $n_{x}$ or $n_{y}$). The error in our experiment is an approximated
$L_{2}$ relative error. Recalling definition of $\rho_{i,n}$ in (21), we
denote
$\hat{\rho}_{in}^{n_{x},n_{y}}=(P_{i}\phi_{n},\hat{\rho}^{n_{x},n_{y}})_{x,y}$,
where $\hat{\rho}^{n_{x},n_{y}}$ is the numerical solution computed in Alg.1
on the product space of $n_{x}$ Legendre polynomial and $n_{y}$ Hermite basis.
Then the error in the convergence study is defined as,

$\displaystyle
e_{n_{x},n_{y}}=\sqrt{\frac{\displaystyle\sum_{i,n}(\hat{\rho}_{in}^{n_{x},n_{y}}-\rho_{in}^{n_{x}^{\prime},n_{y}^{\prime}})^{2}}{\displaystyle\sum_{i,n}(\rho_{in}^{n_{x}^{\prime},n_{y}^{\prime}})^{2}}},$
(31)

where $\rho^{n_{x}^{\prime},n_{y}^{\prime}}$ is the solution computed from
finer discretization, which is considered as ground truth. In our study, we
take $n_{x}^{\prime}=16$, $n_{y}^{\prime}=300$ on the interval
$I=[0,100/1024]$.

Clearly, $\hat{\rho}$ converges to around machine precision ($10^{-15}$) when
$n_{x}\geq 10$. We remark that the convergence with respect to $n_{x}$ is
significantly faster than the convergence with respect to $n_{y}$, since we
are considering a relatively thin leaf in the $x$ direction, so the decay of
coefficients is faster in the $x$ direction. In Fig.6(b), in addition to the
error, we also plot the norm of the Fourier coefficients of the corresponding
Hermite function basis. Clearly the error decays at a similar rate to the
ground truth coefficients.

(a) Convergence with respect to $n_{x}$, i.e. $e_{n_{x},n_{y}^{\prime}}$
versus $n_{x}$

(b) Convergence with respect to $n_{y}$, i.e. $e_{n_{x}^{\prime},n_{y}}$
versus $n_{y}$

Figure 6: Convergence of $\rho$

#### Convergence of scattering matrices by Alg.2

Now we are going to investigate the error of computing scattering matrix by by
Alg.2. We consider $V=V_{1}\chi_{I}$, see (30), on the interval $I=[-1,1]$ and
$E=1.8$. The the scattering matrix under consideration is a $3\times 3$ square
matrix representing the incoming/outgoing coefficients between $(0,-1)$,
$(1,1)$ and $(1,-1)$ mode. We apply Alg.2 with $n_{x}^{\prime}=30$,
$n_{y}^{\prime}=200$ and $L=0$ to generate ground truth scattering matrix
matrix. In Fig.7 we show the Frobenius norm of error in computing scattering
matrix with various $n_{x}$ and $n_{y}$ configuration and different levels of
merging.

In Fig.7, we compare the performance of our merging algorithm, Alg.2 with
direct solving on the whole interval, Alg.1. We the claim that the leaf-and-
merge algorithm is very efficient and accurate in computing eigenfunctions in
a long slab. This is because, we can divide the slab to arbitrary small
pieces. For each piece, the discretization in $x$ direction only takes few
modes ($6$, see dashed line in Fig.7(a)). Meanwhile the overall accuracy is
high as the merge/recover procedure is purely algebraic and its error comes
from roundoff in the computation procedure (see difference of dashed and solid
line when $n_{x}$ or $n_{y}$ is sufficiently large in Fig.7(a) and Fig.7(b).
In Fig.7(c), compare the outcome of the algorithm with different overall level
configuration while keeping $n_{x}=6$ to be the same in each leaf.

(a) Convergence with respect to $n_{x}$

(b) Convergence with respect to $n_{y}$

(c) Convergence with respect to level of merging $L$

Figure 7: Convergence of scattering matrix with Alg.2, error measured by
Frobenius norm

##  6 Applications to conductivity and scattering matrix computations

This section presents three applications for the above algorithm. Firstly, in
section 6.1, we derive an expression for the conductivity $\sigma_{I}$ in (4)
in terms of generalized eigenfunctions of the perturbed system $(H+V-E)\psi=0$
and show that it is indeed quantized and equal to $-1$ with high accuracy
independent of the energy range described by $\varphi^{\prime}(E)$ and of the
perturbation $V$.

In a second application, described in section 6.2, we study the far-field
scattering matrix associated with the generalized eigenvectors of
$(H+V-E)\psi=0$ in the presence of a compactly supported perturbation $V$. We
derive the relationship between scattering matrix and conductivity in (39). We
also investigate how the scattering matrix is related to length of the slab.

Lastly, in section 6.3, we present a special case that $V\mathcal{G}$ has an
eigenvalue equal to $-1$. In such situation, the corresponding eigenfunction
$\rho_{loc}$ recovers an eigenfunction of $H+V-E$, namely $\psi_{loc}$, that
vanishes outside the support of $V$ in $x$ direction.

###  6.1 Asymmetric transport conductivity

#### Current conservation.

Before presenting and computing the line conductivity modeling asymmetric
transport, we consider the related notion of current conservation.

For $\psi(x,y)$ a generalized eigenfunction of the problem $(H+V-E)\psi=0$
(including both incoming and outgoing components), we define the current as

$\displaystyle\tilde{j}_{\psi}(x)=(\psi(x,\cdot),\sigma_{3}\psi(x,\cdot))_{y},$
(32)

where we recall that $(\cdot,\cdot)_{y}$ is the standard inner product on
$L^{2}(\mathbb{R}_{y};\mathbb{C}^{2})$.

For a (non-local) Hermitian perturbation $V$ that only couples propagating
modes, it was shown in [13] that the current $\tilde{j}_{\psi}(x)$ was
conserved, i.e., independent of $x\in\mathbb{R}$. We now generalize this
current conservation to local (pointwise-multiplication) Hermitian
perturbations $V$.

First we decompose $\psi$ as,

$\displaystyle\psi(x,y)=\displaystyle\sum_{m\in M}\beta_{m}(x)\psi_{m}(x,y),$
(33)

with $\alpha_{m}(x)=\beta_{m}(x)e^{i\xi_{m}x}$ in (25) so that $\beta_{m}(x)$
is constant when $V=0$. For $M\ni m=(n,\epsilon_{m})$ and $M\ni
q=(p,\epsilon_{q})$, we observe that

$\displaystyle(\psi_{m},\sigma_{3}\psi_{q})_{y}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{e}e^{i\xi_{q}x}(\phi_{m},\sigma_{3}\phi_{q})_{y}$

When $n\not=p$, then $(\psi_{m},\sigma_{3}\psi_{q})_{y}=0$ by orthogonality of
the Hermite functions $\varphi_{n}$. When $n=p$, then

$\displaystyle(\psi_{m},\sigma_{3}\psi_{q})_{y}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{e}e^{i\xi_{q}x}c_{m}c_{q}(2n-(E-\bar{\xi}_{m})(E-\xi_{q})).$ (34)

For propagating modes, $\xi_{m}\in\mathbb{R}$ so that (34) vanishes when
$\xi_{m}=-\xi_{q}$ while for $\xi_{m}=\xi_{q}$,

$\displaystyle(\psi_{m},\sigma_{3}\psi_{m})_{y}=\frac{2n-(E-\xi_{m})^{2}}{2n+(E-\xi_{m})^{2}}=\frac{\xi_{m}}{E}\not=0.$

For evanescent mode $i\xi_{m}\in\mathbb{R}$, then (34) vanishes when
$\xi_{m}=\xi_{q}$. When $\xi_{m}=-\xi_{q}$, we find that $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{e}e^{i\xi_{q}x}=1$, and hence,

$\displaystyle(\psi_{m},\sigma_{3}\psi_{q})_{y}=\frac{2n-E^{2}-i\epsilon_{m}E\sqrt{2n-E^{2}}}{2n}\not=0.$

Here we used

$\displaystyle 2n-(E+\xi_{m})^{2}=2(2n-E^{2}-\xi_{m}E),\quad
c_{m}^{-2}=c_{q}^{-2}=2n+|E-\xi_{m}|^{2}=4n=(c_{m}c_{q})^{-1}.$

Define $\widehat{m}=\begin{cases}m\quad\xi_{m}\in\mathbb{R}\\\ -m\quad
i\xi_{m}\in\mathbb{R}\end{cases}$. Then from the above,
$(\psi_{m},\sigma_{3}\psi_{m^{\prime}})_{y}$ is a nonzero constant independent
of $x$ when $m^{\prime}=\widehat{m}$ while it vanishes otherwise. We observe
that

$\displaystyle D_{x}(\psi,\sigma_{3}\psi)_{y}$
$\displaystyle=D_{x}(\displaystyle\sum_{m,m^{\prime}}(\beta_{m^{\prime}}\psi_{m^{\prime}},\sigma_{3}\beta_{m}\psi_{m})_{y})=\displaystyle\sum_{m}D_{x}(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}\beta_{m})(\psi_{\widehat{m}},\sigma_{3}\psi_{m})_{y}.$ (35)

Substituting (33) in $(H+V-E)\psi=0$ yields

$\displaystyle\displaystyle\sum_{m^{\prime}}(D_{x}\beta_{m^{\prime}})\sigma_{3}\psi_{m^{\prime}}+\beta_{m^{\prime}}V\psi_{m^{\prime}}=0.$
(36)

Multiplying the above by $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
and integrating over $y$ yields

$\displaystyle\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}D_{x}\beta_{m}=-(\beta_{\widehat{m}}\psi_{\widehat{m}},V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}})_{y}.$

Similarly, taking the conjugate of (36) and multiplying by
$\beta_{m}\psi_{m}$, we observe

$\displaystyle\beta_{m}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{D}=-(V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}},\beta_{m}\psi_{m})_{y}.$

Using (35), we finally obtain

$\displaystyle D_{x}(\psi,\sigma_{3}\psi)_{y}$
$\displaystyle=(V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}},\displaystyle\sum_{m}\beta_{m}\psi_{m})_{y}-(\displaystyle\sum_{m}\beta_{\widehat{m}}\psi_{\widehat{m}},V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}})_{y}$
$\displaystyle=(V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}},\displaystyle\sum_{m}\beta_{m}\psi_{m})_{y}-(\displaystyle\sum_{m}\beta_{m}\psi_{m},V\displaystyle\sum_{m^{\prime}}\beta_{m^{\prime}}\psi_{m^{\prime}})_{y}=0.$

The last identity holds since $V$ is Hermitian. This concludes our derivation
that $\tilde{j}_{\psi}(x)$ is independent of $x$ and hence current is
conserved.

#### Numerical verification of current conservation.

As an example, we consider the potential $V=V_{1}\chi_{I}$ in (30), with
support in $I=[0,100]$. We use 10 levels in Alg.2 ($L=10$) to compute the
eigenfunctions of the slab with $n_{x}=12$ and $n_{y}=100$ in each leaf. In
the convergence study in (Fig.7), empirically we saw that for these
parameters, the error in a leaf of a corresponding size (length $\approx 0.1$)
was almost machine precision. Since the merging algorithm is purely algebraic,
we expect similar errors (up to accumulated rounding errors) in the
computation of the final TR matrix and the subsequent solution of $\rho$.

As an example, we consider the case where the energy $E$ is $1.8.$ Since for
this value $2<E^{2}<4$, $\xi_{m}$ is real-valued for three modes
$\psi\equiv\psi_{m}^{V}=\psi_{m,{\rm in}}+\psi_{m,{\rm out}}$ solution of
$(H+V-E)\psi^{V}_{m}=0$ corresponding to $m=(0,-1)$, $m=(1,1)$, and
$m=(1,-1)$. Fig.8 shows the variations in the current about their average
$\tilde{j}_{\psi}(x)-\overline{\tilde{j}_{\psi}(\cdot)}$ for each of these
propagating modes. We observe that current is indeed preserved up to numerical
errors of order $10^{-13}$.

The average values of the currents corresponding to (a), (b), and (c) in the
figure are given approximately by $4.66\,10^{-3}$, $3.06\,10^{-6}$, and
$-0.62$, respectively. These numbers may be explained as follows. The current
of mode $(1,-1)$ is close to the value it would have when $V_{1}=0$ since this
mode is only moderately affected by the potential $V_{1}$. Since the latter
has oscillations in the $x$ variable that concentrate near
$\xi_{(0,-1)}-\xi_{(1,1)}$, it effectively couples the modes $(0,-1)$ and
$(1,1)$. As a result, significant backscattering occurs and the resulting
currents associated to modes $(0,-1)$ and $(1,1)$ are correspondingly smaller.
See discussion in the subsequent section 6.2 for additional detail on
transmission and reflection coefficients.

(a) $\psi_{\rm in}=e^{i\xi_{(0,-1)}x}\varphi_{(0,-1)}$

(b) $\psi_{\rm in}=e^{i\xi_{(1,1)}x}\varphi_{(1,1)}$

(c) $\psi_{\rm in}=e^{i\xi_{(1,-1)}x}\varphi_{(1,-1)}$

Figure 8: Spatial fluctuations of the currents
$\tilde{j}_{\psi}(x)-\overline{\tilde{j}_{\psi}(\cdot)}$ for different
incoming condition, mean current $\overline{\tilde{j}_{\psi}(\cdot)}$ is shown
in the legend.

#### Computation of line conductivity.

The interface conductivity $\sigma_{I}$ introduced in (4) is a physical
observable that models asymmetric transport along the interface $y\sim 0$.
Observing current for modes at energy $E$ and at the interface point $x_{0}$,
we consider the specific conductivity

$\displaystyle 2\pi\sigma_{I}(E,x_{0})={\rm Tr}\ 2\pi
i[H_{V},P]\chi^{\prime}(H_{V})$ (37)

with $P^{\prime}=\delta(x-x_{0})$ and $\chi^{\prime}(H)=\delta(H-E)$ formally.
Computations of general conductivities as in (4) may be obtained as a
superposition of such conductivities $2\pi\sigma_{I}(E,x_{0})$.

In the unperturbed case with $V=0$, the generalized eigenfunctions of $H$ are
given by $\psi_{m}(x;\xi)=e^{ix\xi}\phi_{m}(\xi)$ in (7). Then
$\chi^{\prime}(H)$ can be represented spectrally by

$\displaystyle\chi^{\prime}(H)=\mathcal{F}^{-1}\chi^{\prime}(\hat{H}(\xi))\mathcal{F}=\displaystyle\sum_{m}\mathcal{F}^{-1}\chi^{\prime}(E_{m}(\xi))\Pi_{m}(\xi)\mathcal{F},$

with $\mathcal{F}$ Fourier transform in the first variable and
$\Pi_{m}(\xi)=\phi_{m}(\xi)\otimes\phi_{m}(\xi)$ the rank-one projector for
the $m$th branch at wavenumber $\xi$ corresponding to an energy
$E_{m}(\xi)=\epsilon_{m}\sqrt{\xi^{2}+2n}$. With this, the Schwartz kernel of
the operator $\chi^{\prime}(H)$ is given by

$\displaystyle
K_{\chi^{\prime}(H)}(x,x^{\prime})=\displaystyle\sum_{m}\displaystyle\int_{\mathbb{R}}\dfrac{d\xi}{2\pi}\chi^{\prime}(E_{m}(\xi))\psi_{m}(x;\xi)\otimes\psi_{m}(x^{\prime};\xi).$

For the Dirac equation, $i[H,P]=P^{\prime}(x)\sigma_{3}$, so that for
$\chi^{\prime}(H)=\delta(H-E)$,

$\displaystyle{\rm Tr}\ 2\pi i[H,P]\chi^{\prime}(H)$
$\displaystyle=\displaystyle\sum_{m}\displaystyle\int_{\mathbb{R}}d\xi\chi^{\prime}(E_{m}(\xi))\big{(}\psi_{m}(x_{0};\xi),\sigma_{3}\psi_{m}(x_{0};\xi)\big{)}_{y}$
$\displaystyle=\displaystyle\sum_{m}\Big{|}\dfrac{\partial{\xi_{m}}}{\partial{E}}\Big{|}\big{(}\psi_{m}(x_{0};\xi_{m}),\sigma_{3}\psi_{m}(x_{0};\xi_{m})\big{)}_{y},$

where we have defined the points $\xi_{m}$ such that $E_{m}(\xi_{m})=E$. They
are given by $\xi_{m}=\epsilon_{m}\sqrt{E^{2}-2n}$ and are real-valued as the
above sum over $m$ involves only propagating modes since
$(\phi_{m},\sigma_{3}\phi_{m})=0$ for evanescent modes. Here we have
implicitly used the results obtained in, e.g., [9, 34] that the trace in (37)
may be computed as the integral along the diagonal of the Schwartz kernel of
the operator $2\pi i[H_{V},P]\chi^{\prime}(H_{V})$.

Now, when $V\not=0$ is compactly supported, the Fourier transform may formally
be replaced by a generalized Fourier transform, where on each branch,
$e^{ix\xi}\phi_{m}(\xi)$ is replaced by $\psi^{V}_{m}(x;\xi)$; see [52] for
similar completeness results for the Schrödinger equation.

Here, $\psi_{m}^{V}(x;\xi_{m})$ is the solution in the kernel of $H_{V}-E$
with incoming condition given by $e^{i\xi_{m}x}\phi_{m}(\xi_{m})$ for
$\xi_{m}$ solving $E=E_{m}(\xi)$, i.e. $\xi_{m}=\epsilon_{m}\sqrt{E^{2}-2n}$.

We thus formally obtain the following expression for the conductivity:

$\displaystyle
2\pi\sigma_{I}=\displaystyle\sum_{m}\Big{|}\dfrac{\partial{\xi_{m}}}{\partial{E}}\Big{|}\big{(}\psi_{m}^{V}(x_{0};\xi_{m}),\sigma_{3}\psi^{V}_{m}(x_{0};\xi_{m})\big{)}_{y},\quad\Big{|}\dfrac{\partial{\xi_{m}}}{\partial{E}}\Big{|}=\dfrac{E}{\sqrt{E^{2}-2n}}.$
(38)

The above expression relates the conductivity describing asymmetric transport
to the generalized eigenvectors that our algorithm may estimate robustly and
accurately.

For simplicity of presentation, we decompose the conductivity as follows

$\displaystyle 2\pi\sigma_{I}=\displaystyle\sum_{m}j_{m}(x),\qquad
j_{m}(x):=\Big{|}\dfrac{\partial{\xi_{m}}}{\partial{E}}\Big{|}\big{(}\psi_{m}^{V}(x,\cdot;\xi_{m}),\sigma_{3}\psi^{V}_{m}(x,\cdot;\xi_{m})\big{)}_{y}.$
(39)

Due to the current conservation shown above, $j_{m}(x)$ is in fact independent
of $x$.

#### Numerical Example: Conductivity dependence.

The computation of the conductivity for a given energy $E$ involves all
generalized eigenfunctions of the operator $H+V-E$. We consider a potential
that couples the propagating modes of the unperturbed system given by
$V=V_{0}\chi_{[-1,1]}$, where

$\displaystyle V_{0}$
$\displaystyle=\exp(-y^{2})y\cos(x(\xi_{(0,-1)}-\xi_{(1,1)}))+\exp(-y^{2})y\cos(x(\xi_{(0,-1)}-\xi_{(1,-1)}))$
$\displaystyle+\exp(-y^{2})\cos(x(\xi_{(1,-1)}-\xi_{(1,1)})).$ (40)

Note that in this construction, $V_{0}$ and hence $V$ depend on $E$ implicitly
since the wavenumbers $\xi_{m}$ also do. We then apply Alg.2 with four levels
to calculate the generalized eigenfunctions $\psi_{m}^{V}$ in the interval
$[-1,1]$ with $L=4$ level binary merging. In each leaf, we have $n_{x}=10$,
$n_{y}=100$ to achieve sufficient accuracy. The conductivity is evaluated at
$x_{0}=0$ and $2<E^{2}<4$ ($E>0$). Then we have two branches $m=(0,-1)$ and
$m=(1,\pm 1)$ crossing $E$ with the first branch crossing at $\xi=-E$ and the
second branch crossing at $\xi_{1\pm}=\pm\sqrt{E^{2}-2}$.

In Fig.9, we show $j_{m}$ defined in (39) for different energy levels and
different scalings $\lambda$ of the (energy dependent) potential (where the
perturbation $V$ is replaced by $\lambda V$). We observe that the individual
terms in the sum (38) depend on $E$ and $V$. The change of the current of a
particular eigenfunction given some fixed incoming waves is due to the change
of scattering coefficients when varying $E$ and $V$. However, as a topological
invariant, as expected $\sum_{m}j_{m}$ is close to $-1$ with high accuracy (13
correct digits in this example).

(a) $j_{(0,-1)}$

(b) $j_{(1,1)}$

(c) $j_{(1,-1)}$

(d) $\displaystyle\sum_{m}j_{m}+1$

Figure 9:  Current and conductivity dependence on $E$ and scale $\lambda$ of
$V$ (currents of generalized eigenvectors of $(H+\lambda
V_{0}\chi_{[-1,1]}-E)$).

We also consider the setting of a fixed energy $E_{0}=1.8$ with
$V_{0}=V_{0}(E_{0}=1.8)$ in (40) with wavenumbers $\xi_{m}=\xi_{m}(E_{0})$. We
present in Fig.10 the currents as a function of $E$ close to $E_{0}$ and
$\lambda$ associated to the generalized eigenvectors of $(H+\lambda
V(E_{0}=1.8)-E)$. As above, we observe quantization of the line conductivity
while its individual components (continuously) depend on the energy and
potential fluctuations.

(a) $j_{(0,-1)}$

(b) $j_{(1,1)}$

(c) $j_{(1,-1)}$

(d) $\displaystyle\sum_{m}j_{m}+1$

Figure 10: Current and conductivity dependence on $E$ and scale $\lambda$ of
$V$ (current of generalized eigenvectors of $(H+\lambda
V_{0}(E_{0}=1.8)\chi_{[-1,1]}-E)$).

###  6.2 Computation of far-field scattering matrices

#### Line conductivity and scattering matrix.

Given a potential $V$ compactly supported on $I=[x_{L},x_{R}]$ in the variable
$x$, the far-field scattering matrix may be extracted from the computed TR
matrix by only considering propagating modes. Note from (39) that
$\sqrt{\frac{E}{|\xi_{m}|}}\psi_{m}$ carries normalized current equal to $\pm
1$.

For example, when $E=1.8$, the scattering matrix is the $3\times 3$ submatrix
of the TR matrix related to the propagating modes $(0,-1)$, $(1,1)$ and
$(1,-1)$.

In general, assuming $k$ propagating modes travelling right (with positive
current) and $k+1$ modes travelling left (with negative current so that the
sum of all currents is $2\pi\sigma_{I}=-1$), the scattering matrix $S$ is
defined by,

$\displaystyle S=\left(\begin{array}[]{ll}R_{+}&T_{-}\\\
T_{+}&R_{-}\end{array}\right)$ (43)

where $R_{+}$ is the $k\times(k+1)$ matrix of reflection of the modes from
left to left, $T_{+}$ the $(k+1)\times(k+1)$ transmission matrix of the same
modes to the right (positive values of $y),R_{-}$ the $(k+1)\times k$ matrix
of reflection of the modes from right to right and $T_{-}$ the $k\times k$
matrix of transmission of those modes to the left. Then, from the definition
of the TR matrix, we deduce that

$\displaystyle\begin{cases}R_{+,pn}=e^{ix_{L}\xi_{(n,1)}-ix_{R}\xi_{(p,-1)}}\sqrt{\Big{|}\frac{\xi_{p}}{\xi_{n}}\Big{|}}M_{(p,-1),(n,1)}\\\
T_{+,pn}=e^{ix_{L}\xi_{(n,1)}-ix_{R}\xi_{p+}}\sqrt{\Big{|}\frac{\xi_{p}}{\xi_{n}}\Big{|}}M_{(p,1),(n,1)}\\\
R_{-,pn}=e^{ix_{L}\xi_{(n,-1)}-ix_{R}\xi_{(p,1)}}\sqrt{\Big{|}\frac{\xi_{p}}{\xi_{n}}\Big{|}}M_{(p,1),(n,-1)}\\\
T_{-,pn}=e^{ix_{L}\xi_{(n,-1)}-ix_{R}\xi_{(p,-1)}}\sqrt{\Big{|}\frac{\xi_{p}}{\xi_{n}}\Big{|}}M_{(p,-1),(n,-1)}.\\\
\end{cases}$

The scattering matrix is constrained by the quantization of the asymmetric
transport $2\pi\sigma_{I}=-1$. Substituting the decomposition (25),
$\psi=\sum_{m}\alpha_{m}\phi_{m}$ into the individual term of the conductivity
in (39), we obtain

$\displaystyle
j_{m}(x_{0})=\Big{|}\dfrac{\partial{\xi_{m}}}{\partial{E}}\Big{|}\displaystyle\sum_{q,q^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}\alpha_{q^{\prime};m}(x_{0})(\phi_{q}^{\prime},\sigma_{3}\phi_{q}),$ (44)

here $q=(p,\varepsilon_{q})$, $m=(n,\varepsilon_{m})$, and $\alpha_{q;m}$
denotes the Fourier coefficients of $\psi_{m}^{V}$ projected onto the
$\\{\psi\\}$ basis. Since evanescent modes vanish outside of the support of
$V$, while we proved current remains constant in $x_{0}$, we may choose
$|x_{0}|\to\infty$ in (44) to limit the summation to propagating modes only.

More precisely, when $\varepsilon_{m}>0$, the incoming condition is a right
travelling modes. We take $x_{0}\to\infty$ in (44) and also substitute
$(\phi_{q}^{\prime},\sigma_{3}\phi_{q})=\delta_{q}(q^{\prime})\frac{\xi_{m}}{E}$
to get,

$\displaystyle j_{m}(x_{0})$
$\displaystyle=\frac{E}{|\xi_{m}|}\displaystyle\sum_{q}|\alpha_{q;m}|^{2}(x_{0})\frac{\xi_{q}}{E}\
=\
\displaystyle\sum_{\varepsilon_{q}>0}|\alpha_{q;m}|^{2}(x_{0})\frac{\xi_{q}}{|\xi_{m}|}.$

The second equation is due to incoming condition $\psi^{V}_{{\rm
in};m}=\psi_{m}$ which is solely right travelling. Still for
$\varepsilon_{m}>0$, we have

$\displaystyle|\alpha_{q;m}|^{2}(x_{0})\frac{\xi_{q}}{\xi_{m}}=|M_{q,m}|^{2}\Big{|}\frac{\xi_{q}}{\xi_{m}}\Big{|}=|T_{+,pn}|^{2},$

so that

$\displaystyle j_{m}=\sum_{p}|T_{+,pn}|^{2}.$ (45)

For $\varepsilon_{m}<0$, we consider $x_{0}\to-\infty$ instead since each
$j_{m}(x)$ is independent of $x$ and obtain that

$\displaystyle j_{m}=-\sum_{p}|T_{-,pn}|^{2}.$ (46)

Summing up the above calculations, we thus deduce that the conductivity in
(39) is given by

$\displaystyle 2\pi\sigma_{I}={\rm tr}\ T_{+}^{*}T_{+}-{\rm tr}\
T_{-}^{*}T_{-}.$

#### Scattering matrix of a slab of length $l$.

We now investigate the properties of the scattering matrix beyond the above
quantized constraint, and in particular how it depends on the potential $V$.

For computational purpose, instead of the binary merging proposed in Alg.2, we
merge the leaves sequentially. More precisely, for $i=1\cdots N$, let $M_{i}$
be the sequence of leaf TR matrices on the interval $I_{i}$. We start from the
first leaf TR matrix $M_{1}$ and merge it with the adjacent $M_{2}$ to get a
scattering matrix $I_{1}\cup I_{2}$. Iterating the process until $k=N-1$ to
get the TR matrix of $I_{1}\cup\cdots\cup I_{k}$, we merge it with $M_{k+1}$
to get the TR matrix for $I_{1}\cup\cdots\cup I_{k}\cup I_{k+1}$. This
modified merging procedure is slightly more convenient for showing the
relationship between the length of the interval and the corresponding
scattering matrix.

#### Topologically non-trivial back-scattering.

Fig.6.2 displays the scattering matrix for energy $E=1.8$ (so that
$2<E^{2}<4$) corresponding to the three propagating modes $(0,-1)$, $(1,-1)$,
and $(1,1)$; recall Fig. 1. The parameters used in the simulation are
identical to those used for the current calculation, i.e. for computing the
results in numerical verification of current conservation in section 6.1.

The main observation is that in the presence of sufficiently strong
fluctuations $V$ (i.e., when the support of $V$ is sufficiently large), the
modes $(0,-1)$ and $(1,1)$ are almost fully back-scattered while the mode
$(1,-1)$ is almost fully transmitted. The absence of backscattering is
therefore not topologically protected. Only the global asymmetry
$2\pi\sigma_{I}=-1$ is, which is consistent with the analysis in [13] in the
presence of random fluctuations $V$.

We note that the oscillations in the scattering coefficient of the $(1,-1)$
mode, while small, are not numerical artifacts. These oscillations, which
depend on the choice of $V$, are significantly larger than accuracy of the
algorithm as exemplified by current conservation in, e.g., Fig.8. Although we
chose a potential $V=V_{0}\chi_{[0,l]}$ to maximize coupling between the
$(0,-1)$ and $(1,1)$ mode, this does not guarantee that the remaining modes
are exactly protected.

(a) $(0,-1)\to(0,-1)$

(b) $(1,1)\to(0,-1)$

(d) $(0,-1)\to(1,1)$ (e) $(1,1)\to(1,1)$ (f) $(1,-1)\to(1,1)$ (g)
$(0,-1)\to(1,-1)$ (h) $(1,1)\to(1,-1)$ (i) $(1,-1)\to(1,-1)$ Figure 11:
Entries of scattering matrix against length of slab $l$,
$V=V_{1}\chi_{[0,l]}$, $E=1.8$

The above strong backscattering is in sharp contrast with the setting of
energies $0<E^{2}<2$. In that case, the only allowed propagating mode is
$(0,-1)$. Back-scattering is then obviously absent since no mode carries a
positive current. We thus obtain a back-scatter-free setting for a combination
of topological ($\sigma_{I}=1)$ and energetic ($E^{2}<2$) reasons.

Consider the case $E=1.2$ with a potential $V=V_{4}\chi_{[0,l]}$ where

$\displaystyle V_{4}(x,y)=\exp(-y^{2})\cos(Ex)$

. In Fig.6.2, we find the scattering coefficient $(0,-1)\to(0,-1)$, which is a
$1\times 1$ (transmission) matrix. The coefficient oscillates as a periodic
function of $l$ with absolute value necessarily equal to $1$. However, the
potential does have an influence on the phase of this mode; see Fig.6.2(b)).

(a) $|S_{(1,-1),(1,-1)}|-1$

(b) $arg(S_{(1,-1),(1,-1)})$

Figure 12: Scattering coefficient $(1,-1)\to(1,-1)$ against length of slab
$l$, $V=V_{4}\chi_{[0,l]}$, $E=1.2$

#### Coupling generating oscillatory transmission.

We now consider the potential

$\displaystyle V=V_{2}\chi_{[0,l]},\quad
V2=\exp\left(-y^{2}\right)y\cos\left(\left(\xi_{(0,-1)}-\xi_{(1,-1)}\right)x\right),$

which couples the two left travelling modes $(0,-1)$ and $(1,-1)$. The
computational configuration is the same as for $V=V_{1}\chi_{[0,l]}$ above. In
Fig.6.2, we show the entries of the scattering matrix as a function of the
length of the slab $l$.

Recall that for $V=V_{1}\chi_{[0,l]}$, increasing the length of the scattering
interval increases the coupling between the modes $(0,-1)$ and $(1,1)$, which
are eventually fully backscattered.

In contrast, when $V=V_{2}\chi_{[0,l]}$, the modes with negative current
$(0,-1)$ and $(1,-1)$ are efficiently coupled, which results in the periodic
pattern observed in Fig.6.2.

(a) $(0,-1)\to(0,-1)$

(b) $(1,1)\to(0,-1)$

(d) $(0,-1)\to(1,1)$ (e) $(1,1)\to(1,1)$ (f) $(1,-1)\to(1,1)$ (g)
$(0,-1)\to(1,-1)$ (h) $(1,1)\to(1,-1)$ (i) $(1,-1)\to(1,-1)$ Figure 13:
Entries of scattering matrix against length of slab $l$,
$V=V_{2}\chi_{[0,l]}$, $E=1.8$.

#### Scattering matrix at higher energy.

Our algorithm allows us to compute generalized eigenfunctions of arbitrary
energy level and matrix-valued (Hermitian) perturbations. Consider $E=2.2$,
which corresponds to $5$ propagating modes, i.e., $(0,-1)$, $(1,\pm)$,
$(2,\pm)$, and a perturbation of the form $V=V_{3}\chi_{[-1,1]}$ with

$\displaystyle
V_{3}=\begin{pmatrix}\exp(-y^{2})\cos((\xi_{(1,-1)}-\xi_{(1,1)})x)&0\\\
0&\exp(-y^{2})\cos((\xi_{(2,-1)}-\xi_{(2,1)})x)\end{pmatrix},$

which couples $(0,-1)$ with $(2,\pm)$ and $(1,-1)$ with $(1,1)$ modes. Since
$V_{3}$ is an even function, there is no coupling between modes with adjacent
energy, e.g. $0$ (or $(2,\pm 1)$) with $(1,\pm 1)$. In Fig.6.2 we show the
absolute value of the scattering matrix associated to these $5$ modes.
Firstly, we see that computed coefficients that are supposed to vanish between
the non-coupled modes are indeed so up to machine precision, (see (b), (c),
(f), (i), (j), (k), (n), (o), (q), (r), (v), (w)). Secondly, all modes except
$(0,-1)$ are significantly back-scattered when the slab is long enough, while
the $(0,-1)$ mode strongly propagates, although it admits moderate coupling
with $2\pm$ modes.

(a) $(0,-1)\to(0,-1)$

(b) $(1,1)\to(0,-1)$

(c) $(1,-1)\to(0,-1)$

(d) $(2,1)\to(0,-1)$

(f) $(0,-1)\to(1,1)$ (g) $(1,1)\to(1,1)$ (h) $(1,-1)\to(1,1)$ (i)
$(2,1)\to(1,1)$ (j) $(2,-1)\to(1,1)$ (k) $(0,-1)\to(1,-1)$ (l)
$(1,1)\to(1,-1)$ (m) $(1,-1)\to(1,-1)$ (n) $(2,1)\to(1,-1)$ (o)
$(2,-1)\to(1,-1)$ (p) $(0,-1)\to(2,1)$ (q) $(1,1)\to(2,1)$ (r)
$(1,-1)\to(2,1)$ (s) $(2,1)\to(2,1)$ (t) $(2,-1)\to(2,1)$ (u)
$(0,-1)\to(2,-1)$ (v) $(1,1)\to(2,-1)$ (w) $(1,-1)\to(2,-1)$ (x)
$(2,1)\to(2,-1)$ (y) $(2,-1)\to(2,-1)$ Figure 14: Entries of scattering matrix
against length of slab $l$, $V=V_{3}\chi_{[0,l]}$, $E=2.2$.

We next compute the currents $j_{m}$ from the scattering matrix using (45) and
(46). In Fig.6.2, we see that the currents vanish for the $(1,\pm 1)$ and
$(2,1)$ modes due to (essentially) complete back-scattering while for $(0,-1)$
and $(2,-1)$ mode, they oscillate as the slab length increases. These
calculations also confirm the asymmetric transport properties of the
topological insulator.

(a) $j_{(0,-1)}$

(b) $j_{(1,1)}$

(d) $j_{(2,1)}$ (e) $j_{(2,-1)}$ (f) $\displaystyle\sum_{m}j_{m}$ Figure 15:
Currents from scattering matrix, $j_{m}$ against length of slab $l$,
$V=V_{3}\chi_{[0,l]}$, $E=2.2$.

###  6.3 Existence of Localized Modes

In section 2.3, we argued that $I+V\mathcal{G}$ with $\mathcal{G}=(H-E)^{-1}$
(with outgoing conditions) was invertible as soon $-1$ was not an eigenvalue
of the compact operator $V\mathcal{G}$. We now show that $-1$ may indeed be in
the spectrum of $\mathcal{G}$ for a scalar perturbation $V(x)$ independent of
$y$ and compactly supported in $x$. Specifically, we assume that
$V(x,y)=V_{0}\chi_{I}(x)$, where $V_{0}$ is a fixed constant and $I=[0,l]$ is
an interval. We now identify values of $(E,l,V_{0})$ such that $-1$ is in the
point spectrum of $V\mathcal{G}$.

Assume the existence of $n\in\mathbb{N}$ such that $E^{2}-2n<0$ while
$(E-V_{0})^{2}-2n>0$. Thus, $\psi_{n,\pm 1}$ is evanescent outside of $[0,l]$
and propagating inside $[0,l]$. Denote

$\displaystyle\phi_{\pm}=c_{(n,\pm
1)}\begin{pmatrix}{\mathfrak{a}}\varphi_{n}\\\ (E-\xi_{(n,\pm
1)})\varphi_{n}\end{pmatrix}=c_{\pm}\begin{pmatrix}\sqrt{2n}\varphi_{n-1}\\\
(E\mp\xi)\varphi_{n}\end{pmatrix},$ (47)

where $c_{\pm}^{-2}=2n+|E\mp\xi|^{2}$ and $\xi=\sqrt{E^{2}-2n}\in i\mathbb{R}$
and

$\displaystyle\phi_{\pm}^{\prime}=c_{(n,\pm
1)}^{\prime}\begin{pmatrix}{\mathfrak{a}}\varphi_{n}\\\ (E-\xi_{(n,\pm
1)}^{\prime})\varphi_{n}\end{pmatrix}=c_{\pm}^{\prime}\begin{pmatrix}\sqrt{2n}\varphi_{n-1}\\\
(E^{\prime}\mp\xi^{\prime})\varphi_{n}\end{pmatrix},$ (48)

where $(c_{\pm}^{\prime})^{-2}=2n+|E^{\prime}\mp\xi^{\prime}|^{2}$,
$\xi^{\prime}=\sqrt{E^{\prime 2}-2n}\in\mathbb{R}^{+}$, and
$E^{\prime}=E-V_{0}$.

We seek a solution $(H+V-E)\psi=0$ given by $\psi=e^{i\xi(x-l)}\phi_{+}$ for
$x\geq l$, which vanishes at $+\infty$. At $x=l$, we have by continuity

$\displaystyle\psi(l,y)=a_{+}\phi^{\prime}_{+}(y)+a_{-}\phi^{\prime}_{-}(y),\qquad\begin{pmatrix}a_{+}\\\
a_{-}\end{pmatrix}=(\phi_{+}^{\prime},\phi_{-}^{\prime})^{-1}\phi_{+}.$

Here, $(\phi_{+}^{\prime},\phi_{-}^{\prime})$ is a $2\times 2$ matrix with
columns $\phi^{\prime}_{\pm}$. Solving the above equation over $x\in[0,l]$
yields

$\displaystyle\psi(0,y)=a_{+}e^{-i\xi^{\prime}l}\phi^{\prime}_{+}(y)+a_{-}e^{i\xi^{\prime}l}\phi^{\prime}_{-}(y).$

The function $\psi\in L^{2}(\mathbb{R}^{2};\mathbb{C}^{2})$ is normalizable
only if $\psi$ is proportional to $e^{-i\xi x}\phi_{-}$ for $x<0$, which in
turn is equivalent to

$\displaystyle b_{+}a_{+}e^{-i\xi^{\prime}l}+b_{-}a_{-}e^{i\xi^{\prime}l}=0,$

where $(b_{+},b_{-})$ is the first row of
$(\phi_{+},\phi_{-})^{-1}(\phi_{+}^{\prime},\phi_{-}^{\prime})$.

Simple but tedious calculations show that

$\displaystyle\frac{b_{+}a_{+}}{b_{-}a_{-}}=-\frac{(E-E^{\prime}+\xi+\xi^{\prime})(E-E^{\prime}-\xi-\xi^{\prime})}{(E-E^{\prime}-\xi+\xi^{\prime})(E-E^{\prime}+\xi-\xi^{\prime})}.$

Recalling $\xi\in i\mathbb{R}$ and $\xi^{\prime}\in\mathbb{R}$, we have
$|\frac{b_{+}a_{+}}{b_{-}a_{-}}|=1$ and thus need:

$\displaystyle{\rm
Arg}(\frac{(E-E^{\prime}+\xi+\xi^{\prime})(E-E^{\prime}-\xi-\xi^{\prime})}{(E-E^{\prime}-\xi+\xi^{\prime})(E-E^{\prime}+\xi-\xi^{\prime})})=2\xi^{\prime}l+2k\pi,\quad
k=0,\pm 1,\pm 2,\cdots.$

Any solution $(E,V_{0},l)$ of the above equation provides a localized model
associated to the eigenvalue $E$ of the operator $H+V$.

As an illustration, we consider the case with $E=1.8$, $\xi^{\prime}=1$ and
$n=2$ while $k=0$. We then find the numerical approximation $l\approx 1.5422$.

In the discretization of Alg.1, we let $(n_{x},n_{y})=(60,20)$ to resolve
local modes related to $\phi_{(2,\pm 1)}$. At step 3, we found that
$I+\hat{V}\hat{G}$ was not stably invertible. Numerically, by applying
singular value decomposition on $I+\hat{V}\hat{G}$, we found that the
condition number, equal to the ratio of largest and smallest eigenvalues, is
$3.619\times 10^{15}$. The smallest two singular values are $0.5431$ and
$4.88\times 10^{-16}$. Let $\hat{\rho}$ be the eigenvector associated with the
latter, which forms the null space of $I+\hat{V}\hat{G}$. As a verification,
we compute that $|(I+\hat{V}\hat{G})\hat{\rho}|_{l_{2}}\approx 0$
($2.029\times 10^{-15}$). In Fig.16, we present the corresponding localized
mode $\psi=\int G\hat{\rho}$ solution of $(H+V-E)\psi\approx 0$.

(a) $Re(\psi_{1})$

(b) $Re(\psi_{2})$

(c) $Im(\psi_{1})$

(d) $Im(\psi_{2})$

Figure 16: Localized mode computed by $\psi=\int G\hat{\rho}$.

##  7 Conclusions

This paper presented a fast and accurate algorithm to compute generalized
eigenfunctions of two-dimensional operators displaying a domain wall that
separates two insulating half-spaces. The main assumption in our construction
is that the Green’s function of the unperturbed operator may be described
reasonably explicitly. Eigenfunctions of perturbed operators are then computed
using an integral formulation for a density $\rho$ supported on the same
domain as the perturbation. This allows us to construct generalized
eigenfunctions with support in the whole of $\mathbb{R}^{2}$ from a compactly
supported density.

Our main application was a Dirac operator with linear domain wall, which finds
applications in the analysis of topological phases of matter. Potential
generalizations of the algorithm include block-diagonal systems of Dirac
operators, which also appear naturally in the study of topological insulators
[9, 13], as well as three dimensional generalizations to high-order
topological insulators [34].

We showed a number of applications of the algorithm to quantify the transport
properties of Dirac operators with domain wall. In particular, we deri=ed an
expression for a quantized interface conductivity describing asymmetric
transport in terms of the computed generalized eigenfunctions. We were then
able to confirm numerically that the conductivity was indeed integral valued
up to very high accuracy. We also computed the full far-field scattering
matrix generated by the perturbations. This matrix took the form of a $2n+1$
times $2n+1$ matrix for all energies $E$ such that $2n<E^{2}<2n+2$. Such
matrices may be computed for energies $E$ not belonging to the point spectrum
of the perturbed operator $H+V$. Note that such a point spectrum is
necessarily embedded in the (absolutely) continuous spectrum $\sigma_{\rm
ac}(H+V)=\mathbb{R}$, a fact that would not be allowed for Schrödinger or
Dirac operators in the absence of a domain wall. As future research, we plan
to use the algorithm to analyze the structure of the point spectrum and
possibly resonances for Dirac and related operators.

## Acknowledgment.

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

## References

  * [1] Edward Witten.  Three lectures on topological phases of matter.  Nuovo Cimento Rivista Serie, 39(7):313–370, 2016.
  * [2] Roderich Moessner and Joel E Moore.  Topological Phases of Matter.  Cambridge University Press, 2021.
  * [3] Bogdan Andrei Bernevig.  Topological Insulators and Topological Superconductors.  Princeton University Press, 2013.
  * [4] Emil Prodan and Hermann Schulz-Baldes.  Bulk and boundary invariants for complex topological insulators.  Springer verlag, Berlin, 2016.
  * [5] Masatoshi Sato and Yoichi Ando.  Topological superconductors: a review.  Reports on Progress in Physics, 80(7):076501, 2017.
  * [6] Grigory E. Volovik.  Nonlinear Phenomena in Condensed Matter: Universe in a Helium Droplet.  Springer, 1989.
  * [7] Michel Fruchart and David Carpentier.  An introduction to topological insulators.  Comptes Rendus Physique, 14(9-10):779–815, 2013.
  * [8] JP Lee-Thorp, MI Weinstein, and Y Zhu.  Elliptic operators with honeycomb symmetry: Dirac points, edge states and applications to photonic graphene.  arXiv preprint arXiv:1710.03389, 2017.
  * [9] Guillaume Bal.  Continuous bulk and interface description of topological insulators.  Journal of Mathematical Physics, 60(8):081506, 2019.
  * [10] Alexis Drouot and MI Weinstein.  Edge states and the valley hall effect.  Advances in Mathematics, 368:107142, 2020.
  * [11] Guillaume Bal.  Topological invariants for interface modes.  To appear in Comm. P.D.E., 2022.
  * [12] Jean-Pierre Fouque, Josselin Garnier, George Papanicolaou, and Knut Solna.  Wave propagation and time reversal in randomly layered media, volume 56.  Springer Science & Business Media, 2007.
  * [13] Guillaume Bal.  Topological protection of perturbed edge states.  Communications in Mathematical Sciences, 17(1):193–225, 2019.
  * [14] Ling Lu, John D. Joannopoulos, and Marin Soljačić.  Topological photonics.  Nature Photonics, 8(11):821–829, 2014.
  * [15] Pierre Delplace, JB Marston, and Antoine Venaille.  Topological origin of equatorial waves.  Science, 358(6366):1075–1077, 2017.
  * [16] Anton Souslov, Kinjal Dasbiswas, Michel Fruchart, Suriyanarayanan Vaikuntanathan, and Vincenzo Vitelli.  Topological waves in fluids with odd viscosity.  Physical review letters, 122(12):128001, 2019.
  * [17] Kendall Atkinson and Weimin Han.  Numerical solution of fredholm integral equations of the second kind.  In Theoretical Numerical Analysis, pages 473–549. Springer, 2009\.
  * [18] Leslie Greengard and Vladimir Rokhlin.  On the numerical solution of two-point boundary value problems.  Communications on pure and applied mathematics, 44(4):419–452, 1991\.
  * [19] N. Beams, A. Gillman, and R. Hewett.  A parallel implementation of a high order accurate solution technique for variable coefficient helmholtz problems.  Computers and Mathematics with Applications, 79:996–1011, 2020\.
  * [20] S. Hao and P.G. Martinsson.  A direct solver for elliptic pdes in three dimensions based on hierarchical merging of poincaré–steklov operators.  J. Comput. Appl. Math, 308:419––434, 2016.
  * [21] A. Gillman, P. Young, and P.G. Martinsson.  A direct solver with o(n) complexity for integral equations on one-dimensional domains.  Frontiers of Mathematics in China, 7(2):217–247, 2012.
  * [22] Adrianna Gillman and P.G. Martinsson.  An o(n) algorithm for constructing the solution operator to 2d elliptic boundary value problems in the absence of body loads.  Advances in Computational Mathematics, 40:773–796, 2014.
  * [23] Adrianna Gillman, A. Barnett, and P.G. Martinsson.  A spectrally accurate direct solution technique for frequency-domain scattering problems with variable media.  BIT Numerical Mathematics, 55:141–170, 2015.
  * [24] D. Fortunato, N. Hale, and A. Townsend.  The ultraspherical spectral element method.  Journal of Computational Physics, 436, 2021.
  * [25] Martinsson P.G.  A fast direct solver for a class of elliptic partial differential equations.  J. Sci. Comput., 38:316–330, 2009.
  * [26] P.G. Martinsson.  A direct solver for variable coefficient elliptic pdes discretized via a composite spectral collocation method.  Journal of Computational Physics, 242:460–479, 2013.
  * [27] E. Corona, P.G. Martinsson, and D. Zorin.  An o(n) direct solver for integral equations in the plane.  Advances in Computational Harmonic Analysis, 38:284–317, 2015.
  * [28] V. Minden, K. Ho, A. Damle, and L. Ying.  A recursive skeletonization factorization based on strong admissibility.  SIAM Multiscale Modeling and Simulation, 15, 2017.
  * [29] Lifeng Li.  Note on the s-matrix propagation algorithm.  JOSA A, 20(4):655–660, 2003.
  * [30] Fu Kexiang, Wang Zhiheng, Zhang Dayue, Zhang Jing, and Zhang Qizhi.  A modal theory and recursion rtcm algorithm for gratings of deep grooves and arbitrary profile.  Science in China Series A: Mathematics, 42(6):636–645, 1999.
  * [31] Tom G Mackay and Akhlesh Lakhtakia.  The transfer-matrix method in electromagnetics and optics.  Synthesis Lectures on Electromagnetics, 1(1):1–126, 2020.
  * [32] Anh-Luan Phan and Dai-Nam Le.  Electronic transport in two-dimensional strained dirac materials under multi-step fermi velocity barrier: transfer matrix method for supersymmetric systems.  The European Physical Journal B, 94(8):1–16, 2021.
  * [33] A. Barnett, Bradley Nelson, and J. Matthew Mahoney.  High-order boundary integral equation solution of high frequency wave scattering from obstacles in an unbounded linearly stratified medium.  Journal of Computational Physics, 297:407–426, 2015.
  * [34] Guillaume Bal.  Topological charge conservation for continuous insulators.  arXiv preprint arXiv:2106.08480, 2021.
  * [35] Alexis Drouot.  Microlocal analysis of the bulk-edge correspondence.  Communications in Mathematical Physics, 383(3):2069–2112, 2021\.
  * [36] Peter Elbau and Gian-Michele Graf.  Equality of bulk and edge hall conductance revisited.  Communications in mathematical physics, 229(3):415–432, 2002.
  * [37] A Elgart, GM Graf, and JH Schenker.  Equality of the bulk and edge hall conductances in a mobility gap.  Communications in mathematical physics, 259(1):185–221, 2005.
  * [38] Hermann Schulz-Baldes, Johannes Kellendonk, and Thomas Richter.  Simultaneous quantization of edge and bulk hall conductivity.  Journal of Physics A: Mathematical and General, 33(2):L27, 2000\.
  * [39] Guillaume Bal and Daniel Massatt.  Multiscale invariants of Floquet topological insulators.  Multiscale Modeling & Simulation, 20(1):493–523, 2022.
  * [40] Solomon Quinn and Guillaume Bal.  Approximations of interface topological invariants.  arXiv preprint arXiv:2112.02686, 2021.
  * [41] Matthew J Colbrook, Andrew Horning, Kyle Thicke, and Alexander B Watson.  Computing spectral properties of topological insulators without artificial truncation or supercell approximation.  arXiv preprint arXiv:2112.03942, 2021.
  * [42] Guillaume Bal, Simon Becker, Alexis Drouot, Clotilde Fermanian Kammerer, Jianfeng Lu, and Alexander Watson.  Edge state dynamics along curved interfaces.  arXiv preprint arXiv:2106.00729, 2021.
  * [43] Guillaume Bal, Simon Becker, and Alexis Drouot.  Magnetic slowdown of topological edge states.  arXiv preprint arXiv:2201.07133, 2022.
  * [44] Pipi Hu, Peng Xie, and Yi Zhu.  Traveling edge states in massive dirac equations along slowly varying edges.  arXiv preprint arXiv:2202.13653, 2022.
  * [45] Lawrence C Evans.  Partial differential equations, volume 19.  American Mathematical Soc., 2010.
  * [46] Anne Berthier and Vladimir Georgescu.  On the point spectrum of dirac operators.  Journal of functional analysis, 71(2):309–338, 1987.
  * [47] Michael Reed and Barry Simon.  Methods of modern mathematical physics, IV: Analysis of operators.  Academic press, 1978.
  * [48] Michael Reed and Barry Simon.  Methods of modern mathematical physics, III: Scattering Theory.  Elsevier, 1979.
  * [49] Michael Reed and Barry Simon.  Methods of modern mathematical physics, I: Functional analysis.  Elsevier, 1972.
  * [50] Rainer Kress, V Maz’ya, and V Kozlov.  Linear integral equations, volume 82.  Springer, 1989.
  * [51] J-.Y. Lee and L. Greengard.  A fast adaptive numerical method for stiff two-point boundary value problems.  SIAM J. Sci. Comput., 18:403–429, 1997.
  * [52] Barry Simon.  Schrödinger semigroups.  Bulletin of the American Mathematical Society, 7(3):447–526, 1982\.