# Non-Bloch band theory of generalized eigenvalue problems

Kazuki Yokomizo  Department of Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo, 113-0033, Japan  Taiki Yoda  Department of Physics, Tokyo
Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8551, Japan
Yuto Ashida  Department of Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo, 113-0033, Japan  Institute for Physics of Intelligence, The
University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan

###### Abstract

Waves in a variety of fields in physics, such as mechanics, optics,
spintronics, and nonlinear systems, obey generalized eigenvalue equations. To
study non-Hermitian physics of those systems, in this work, we construct a
non-Bloch band theory of generalized eigenvalue problems. Specifically, we
show that eigenvalues of a transfer matrix lead to a certain condition imposed
on the generalized Brillouin zone, which allows us to develop a theory to
calculate the continuum bands. As a concrete example, we examine the non-
Hermitian skin effect of photonic crystals composed of chiral metamaterials by
invoking our theoretical framework. When the medium has circularly polarized
eigenmodes, we find that each eigenmode localizes at either of edges depending
on whether it is left- or right-circularly polarized. In contrast, when the
medium only has linearly polarized eigenmodes, every eigenmode localizes to
the edge of the same side independent of its polarization. We demonstrate that
the localization lengths of those eigenmodes can be determined from the chiral
parameters and eigenfrequencies of the photonic crystal.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Non-Hermitian physics has been gaining growing interest over the past few
years, in which one can find rich phenomena with no counterparts in Hermitian
regimes [1, 2]. Non-Hermitian systems belong to a class of nonequilibrium
systems which can be described by certain non-Hermitian operators, whose
eigenvalues characterize the dynamical behaviors. A prominent example is the
non-Hermitian skin effect, in which numerous eigenstates are localized at
boundaries of non-Hermitian systems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 17, 18, 19, 20]. Accordingly, it must be necessary to modify the
conventional Bloch band theory in such a way that the continuum energy bands
reproduce the energy levels of non-Hermitian systems with open boundary
conditions. The resulting theory is called the non-Bloch band theory, which
has been successful to calculate the continuum energy bands in the presence of
the non-Hermitian skin effect [4, 21, 22, 23, 24, 25, 26].

On another front, there exist several physical platforms that can be described
by generalized eigenvalue equations in, e.g., mechanics [27], optics [28, 29],
spintronics [30, 31], and nonlinear systems [32]. Notably, even if all
operators appearing in generalized eigenvalue equations are Hermitian, the
corresponding eigenvalues can be complex-valued [33, 34]. Indeed, it has been
proposed that bosonic Bogoliubov-de Gennes systems exhibit the non-Hermitian
skin effect in the regime where the pseudo-Hermiticity is broken [35, 36].
Aside from this example, several non-Hermitian phenomena in various systems
described by generalized eigenvalue equations have been investigated [37, 38,
39]. It remains, however, unclear how and even whether the non-Bloch band
theory can be extended to generalized eigenvalue problems.

In our work, we develop a non-Bloch band theory to study non-Hermitian waves
in a one-dimensional (1D) periodic medium described by a generalized
eigenvalue equation. In the context of the non-Bloch band theory, the
continuum bands can be obtained by determining the generalized Brillouin zone
for the complex Bloch wavenumber. We show that the eigenvalues of a transfer
matrix of the medium lead to a certain condition imposed on the generalized
Brillouin zone. The generalized theory proposed here allows us to explore the
non-Hermitian skin effect in photonic crystals composed of chiral
metamaterials. We confirm that the continuum bands obtained from the
generalized Brillouin zone reproduce the eigenvalues of finite-size systems.
In particular, when the medium has circularly polarized eigenmodes, we find
that each eigenmode can localize at either of edges depending on whether it is
left- or right-circularly polarized. This finding suggests the interesting
possibility of controlling which side excitation modes localize at simply by
changing a polarization of incident light while keeping the medium unchanged.
In contrast, when the medium only has linearly polarized eigenmodes, our
analysis shows that every eigenmode localizes at one edge of the system.

Before concluding this section, we shall compare the present work with the
previous studies [25, 40]. The latter has focused on dielectric media with
optical loss and investigated electromagnetic waves in the two-dimensional
plane where the dielectric tensor is anisotropic. There, one of the key
findings was that the combination of the anisotropy and the loss gives rise to
the non-Hermitian skin effect, and the resulting localization depends on
wavenumbers of a specified direction. In photonic crystals composed of the
dielectric media, for instance, the non-Bloch band theory has revealed that
all the continuum bands share the common generalized Brillouin zone.
Meanwhile, in the present work, we investigate standing waves of
electromagnetic waves forming along the stacking direction of photonic
crystals composed of chiral metamaterials. Here, the key ingredient of the
occurrence of the non-Hermitian skin effect is chirality. Indeed, the
localization lengths of the skin modes depend on chirality parameters and
eigenfrequencies of the photonic crystals. In addition, it is remarkable that
the multiple generalized Brillouin zones appear depending on the continuum
bands, making contrast to the previous studies.

The rest of this paper is organized as follows. We introduce a generalized
eigenvalue equation, explain the construction of a transfer matrix, and show
the condition for the generalized Brillouin zone in Sec. II. We next study the
non-Hermitian skin effect of linearly and circularly polarized eigenmodes of
the photonic crystals by invoking the non-Bloch band theory in Sec. III.
Finally, in Sec. IV, we summarize the results and comment on the perspective
of this work.

##  II Non-Bloch band theory

We develop a theory to calculate the continuum bands in 1D continuous periodic
models described by a generalized eigenvalue equation. Our framework is based
on a transfer matrix method, where the eigenvalues of a transfer matrix lead
to the condition for the generalized Brillouin zone.

###  II.1 Generalized eigenvalue equation

We study the waves in a 1D spatially periodic medium which can be modeled by a
generalized eigenvalue equation. We denote a lattice constant by $a$ and take
the convention of the time dependence of the wavefunction to be $e^{-i\omega
t}$. Furthermore, we assume that the wavefunction is written as a
multicomponent vector:
$\Psi\left(z\right)=\left(\psi_{1}\left(z\right),\dots,\psi_{2n}\left(z\right)\right)^{\rm
T}$ with $n\in{\mathbb{N}}$ and $z\in{\mathbb{R}}$. The generalized eigenvalue
equation governing our setup then reads

$\frac{d}{dz}\Psi\left(z\right)=i\frac{\omega}{c}A\left(z\right)\Psi\left(z\right),$
(1)

where $c$ is a positive constant, and $A\left(z\right)$ is a $2n\times 2n$
matrix satisfying $A\left(z+a\right)=A\left(z\right)$. The form of Eq. (1)
naturally appears by rewriting the Maxwell’s equations, as demonstrated later.

We aim to investigate a method to calculate the continuum bands, including the
asymptotic eigenvalues of the system with open boundary conditions in the
limit of a large system size. To this end, we focus on the plane-wave
expansion of the wavefunction given by

$\Psi\left(z\right)=\sum_{n}\tilde{\Psi}\left(k+\frac{2n\pi}{a}\right)\exp\left[i\left(k+\frac{2n\pi}{a}\right)z\right],$
(2)

where $k$ represents the Bloch wavenumber. The generalized eigenvalue equation
can be then rewritten in the form of a secular equation as follows:

$\displaystyle\left(k+\frac{2n\pi}{a}\right)\tilde{\Psi}\left(k+\frac{2n\pi}{a}\right)$
$\displaystyle-\omega\sum_{n^{\prime}}\tilde{A}_{n-n^{\prime}}\tilde{\Psi}\left(k+\frac{2n^{\prime}\pi}{a}\right)=0,$
(3)

where $\tilde{A}_{n}$ is a Fourier coefficient of $A\left(z\right)$. The
localized nature due to the non-Hermitian skin effect can be taken into
account by considering a complex-valued Bloch wavenumber $k$ [4, 21].
Therefore, it is necessary to determine the generalized Brillouin zone to
obtain the continuum bands by solving Eq. (3).

###  II.2 Transfer matrix

We next utilize a transfer matrix to develop the non-Bloch band theory which
is applicable to the system described by Eq. (1). Specifically, we can define
the transfer matrix as follows:

$\Psi\left(z+a\right)=T\Psi\left(z\right),$ (4)

where the wavefunction in a unit cell is transferred to that in the next unit
cell by the transfer matrix $T$, i.e., the transfer matrix plays the same role
as a translation operator. The eigenvalues of the transfer matrix construct a
set of $e^{ika}~{}\left(k\in{\mathbb{C}}\right)$, which gives the generalized
Brillouin zone. We note that, when there is a translation symmetry,
eigenvalues of a transfer matrix belong to a set of
$e^{ika}~{}\left(k\in{\mathbb{R}}\right)$ [41, 42].

The explicit form of the transfer matrix in our setup can be obtained by

$T={\cal P}\exp\left(i\frac{\omega}{c}\int_{0}^{a}dzA\left(z\right)\right),$
(5)

where ${\cal P}$ is a path-ordered product. Importantly, we find that the
generalized Brilloun zone is formed by the eigenvalues of the transfer matrix
satisfying

$\left|\rho_{n}\right|=\left|\rho_{n+1}\right|,$ (6)

where we number $2n$ eigenvalues so as to satisfy

$\left|\rho_{1}\right|\leq\cdots\leq\left|\rho_{2n}\right|.$ (7)

The trajectories of $\beta=re^{i{\rm Re}\left(k\right)a}$ with
$r\equiv\left|\rho_{n}\right|=\left|\rho_{n+1}\right|$ give the generalized
Brillouin zone. Thereby, one can calculate the continuum bands by solving Eq.
(3) with the complex-valued Bloch wavenumber. We note that Eq. (6) can be
proved by using tight-binding models as shown in Appendix A, and it has also
been recently proved in continuous models [26]. In Sec. III, we shall
demonstrate that the calculated continuum bands indeed reproduces the discrete
eigenvalues of photonic crystals.

##  III Photonic crystal

Our main interest in this section lies in a chiral metamaterial which exhibits
the magnetoelectric coupling via subwavelength structures. We then apply the
non-Bloch band theory developed in Sec. II to photonic crystals composed of
chiral metamaterials, which allows for investigating the non-Hermitian skin
effect in the medium with linearly or circularly polarized eigenmodes.

###  III.1 Setup

Figure 1: Schematic figures of the photonic crystal. (a) Unit cell of the
photonic crystal. The constitutive parameters and thickness of the layer
$j~{}\left(j=1,2\right)$ are set to be
$\left(\varepsilon_{j},\xi_{j},\zeta_{j},\mu_{j}\right)$ and $a_{j}$,
respectively. (b) Finite-size photonic crystal with boundaries. We impose the
perfect electric conductor (PEC) boundary condition on both ends of the
system.

We consider the photonic crystal in which two media exhibiting the
magnetoelectric coupling are alternately stacked in a periodic manner, where
the layer 1 and 2 have the thickness $a_{1}$ and $a_{2}$, respectively [see
Fig. 1(a)]. We then investigate the eigenmodes of the photonic crystal
surrounded by perfect electric conductors (PECs) without external excitation
[see Fig. 1(b)]. We note that the $x$ and $y$ components of the electric
fields vanish at the ends of the system. Let suppose that the electromagnetic
waves form standing waves along the $z$ direction, and the polarization is
perpendicular to the stacking direction. The electromagnetic waves in the
photonic crystal can be described by the Maxwell’s equations

$\displaystyle\left\\{\begin{array}[]{l}{\bm{\nabla}}\times{\bm{E}}=i\omega{\bm{B}},\vspace{5pt}\\\
{\bm{\nabla}}\times{\bm{H}}=-i\omega{\bm{D}},\end{array}\right.$ (10)

associated with the constitutive relation at each layer

$\displaystyle\left(\begin{array}[]{c}{\bm{D}\vspace{5pt}}\\\
{\bm{B}}\end{array}\right)=\left(\begin{array}[]{cc}\varepsilon_{0}\hat{\varepsilon}_{j}&\hat{\xi}_{j}/c\vspace{5pt}\\\
\hat{\zeta}_{j}/c&\mu_{0}\hat{\mu}_{j}\end{array}\right)\left(\begin{array}[]{c}{\bm{E}\vspace{5pt}}\\\
{\bm{H}}\end{array}\right),~{}\left(j=1,2\right),$ (17)

where $\varepsilon_{0}$, $\mu_{0}$, and $c$ are the vacuum permittivity, the
vacuum permeability, and the vacuum speed of light, respectively. Furthermore,
$\hat{\varepsilon}_{j}$ and $\hat{\mu}_{j}$ are the relative permittivity and
permeability tensors, respectively, and $\hat{\xi}_{j}$ and $\hat{\zeta}_{j}$
are the chirality tensors expressing the degrees of the magnetoelectric
coupling. One can then rewrite the Maxwell’s equations as follows:

$\displaystyle\frac{d}{dz}\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\varepsilon_{0}}E_{y}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{y}\left(z\right)\end{array}\right)=i\frac{\omega}{c}A\left(z\right)\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\varepsilon_{0}}E_{y}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{y}\left(z\right)\end{array}\right),$ (26)

where $A\left(z\right)$ is represented by

$\displaystyle
A\left(z\right)=i\left(\begin{array}[]{cc}\sigma_{y}\hat{\zeta}_{j}&\sigma_{y}\hat{\mu}_{j}\vspace{5pt}\\\
-\sigma_{y}\hat{\varepsilon}_{j}&-\sigma_{y}\hat{\xi}_{j}\end{array}\right),$
(29)

when $z$ is in the layer $j$. We note that $\sigma_{0}$ and
$\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$ denote a $2\times 2$ identity
matrix and the Pauli matrices, respectively.

Throughout our work, we study the case where the relative permittivity and
permeability tensors satisfy $\hat{\varepsilon}_{j}=\varepsilon_{j}\sigma_{0}$
and $\hat{\mu}_{j}=\mu_{j}\sigma_{0}$, respectively, and
$i\hat{\xi}_{j}\sigma_{y}$ and $-i\sigma_{y}\hat{\zeta}_{j}$ are
simultaneously diagonalizable. Let $P_{j}$ denote the matrix diagonalizing
$i\hat{\xi}_{j}\sigma_{y}$ and $-i\sigma_{y}\hat{\zeta}_{j}$. The linear
transformations of the electric and magnetic fields at each layer,

$\displaystyle\begin{array}[]{l}\left(\begin{array}[]{c}\tilde{E}_{+}\left(z\right)\vspace{5pt}\\\
\tilde{E}_{-}\left(z\right)\end{array}\right)=P_{j}^{-1}\left(\begin{array}[]{c}E_{x}\left(z\right)\vspace{5pt}\\\
E_{y}\left(z\right)\end{array}\right),\vspace{5pt}\\\
\left(\begin{array}[]{c}\tilde{H}_{+}\left(z\right)\vspace{5pt}\\\
\tilde{H}_{-}\left(z\right)\end{array}\right)=P_{j}^{-1}\sigma_{y}\left(\begin{array}[]{c}H_{x}\left(z\right)\vspace{5pt}\\\
H_{y}\left(z\right)\end{array}\right),\end{array}$ (40)

then allow us to derive the block diagonalization form of Eq. (26) as follows:

$\displaystyle\frac{d}{dz}\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}\tilde{E}_{+}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}\tilde{H}_{+}\left(z\right)\vspace{5pt}\\\
\sqrt{\varepsilon_{0}}\tilde{E}_{-}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}\tilde{H}_{-}\left(z\right)\end{array}\right)=i\frac{\omega}{c}\left(\begin{array}[]{cc}\tilde{A}_{+}\left(z\right)&O\vspace{5pt}\\\
O&\tilde{A}_{-}\left(z\right)\end{array}\right)\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}\tilde{E}_{+}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}\tilde{H}_{+}\left(z\right)\vspace{5pt}\\\
\sqrt{\varepsilon_{0}}\tilde{E}_{-}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}\tilde{H}_{-}\left(z\right)\end{array}\right),$ (51)

where $\tilde{A}_{\pm}\left(z\right)$ with $z$ being in the layer $j$ are
represented by

$\displaystyle\tilde{A}_{\pm}\left(z\right)=\left(\begin{array}[]{cc}-\zeta_{j,\pm}&i\mu_{j}\vspace{5pt}\\\
-i\varepsilon_{j}&-\xi_{j,\pm}\end{array}\right).$ (54)

Here, $\xi_{j,\pm}$ and $\zeta_{j,\pm}$ are the eigenvalues of
$i\hat{\xi}_{j}\sigma_{y}$ and $-i\sigma_{y}\hat{\zeta}_{j}$, respectively. It
is worthwhile to mention that the eigenvectors of $i\hat{\xi}_{j}\sigma_{y}$
or, equivalently, $-i\sigma_{y}\hat{\zeta}_{j}$ determine polarizations of the
eigenmodes. For the sake of later convenience, we denote the eigenvalues of
Eq. (54) by $\lambda_{j,\pm},\bar{\lambda}_{j,\pm}$.

###  III.2 Transfer matrix and generalized Brillouin zone

According to Eq. (51), the transfer matrix of the photonic crystal also has
the block diagonalization form, which is given by

$\displaystyle T=\left(\begin{array}[]{cc}T_{+}&O\vspace{5pt}\\\
O&T_{-}\end{array}\right),$ (57)

where

$T_{\pm}={\cal
P}\exp\left(i\frac{\omega}{c}\int_{0}^{a}dz\tilde{A}_{\pm}\left(z\right)\right).$
(58)

To obtain the generalized Brillouin zone, we calculate the determinants of the
block transfer matrices $T_{\pm}$. Let $U_{j,\pm}~{}\left(j=1,2\right)$ denote
the matrices diagonalizing $T_{\pm}$. Equation (58) can be then written as

$\displaystyle
T_{\pm}=U_{2,\pm}\left(\begin{array}[]{cc}e^{i\omega\lambda_{2,\pm}a_{2}/c}&0\vspace{5pt}\\\
0&e^{i\omega\bar{\lambda}_{2,\pm}a_{2}/c}\end{array}\right)U_{2,\pm}^{-1}$
(61) $\displaystyle\times
U_{1,\pm}\left(\begin{array}[]{cc}e^{i\omega\lambda_{1,\pm}a_{1}/c}&0\vspace{5pt}\\\
0&e^{i\omega\bar{\lambda}_{1,\pm}a_{1}/c}\end{array}\right)U_{1,\pm}^{-1},$
(64)

and one can obtain

$\det
T_{\pm}=\exp\left[i\frac{\omega}{c}\sum_{j=1}^{2}\left(\lambda_{j,\pm}+\bar{\lambda}_{j,\pm}\right)a_{j}\right].$
(65)

We note that, in this case, the following facts enable us to straightforwardly
calculate the generalized Brillouin zone. First, each of the eigenvalues
$\rho_{1,\pm},\rho_{2,\pm}$ of the block transfer matrices $T_{\pm}$ induces
the condition for the generalized Brillouin zone, which means that

$\left|\rho_{1,\pm}\right|=\left|\rho_{2,\pm}\right|.$ (66)

Second, the product of the eigenvalues is equal to the determinant of the
block transfer matrix. Therefore, one can get the generalized Brillouuin zones
as closed loops formed by $r_{\pm}e^{i{\rm Re}\left(k\right)a}$, where
$r_{\pm}$ are given by

$r_{\pm}=\sqrt{\left|\det T_{\pm}\right|}.$ (67)

###  III.3 Linearly polarized eigenmodes

We first investigate the case where the eigenmodes of the photonic crystal are
linearly polarized. The chirality tensors are then given by

$\displaystyle\hat{\xi}_{j}=\left(\begin{array}[]{cc}0&\xi_{xy,j}\vspace{5pt}\\\
\xi_{yx,j}&0\end{array}\right),~{}\hat{\zeta}_{j}=\left(\begin{array}[]{cc}0&\zeta_{xy,j}\vspace{5pt}\\\
\zeta_{yx,j}&0\end{array}\right)$ (72)

with $j=1,2$. One can easily check that the eigenvectors of
$i\hat{\xi}_{j}\sigma_{y}$ and $-i\sigma_{y}\hat{\zeta}_{j}$ are indeed
linearly polarized, i.e., they are proportional to $\left(1,0\right)^{\rm T}$
and $\left(0,1\right)^{\rm T}$. The corresponding eigenvalues are given by
$\left(\xi_{j,+},\xi_{j,-}\right)=\left(-\xi_{xy,j},\xi_{yx,j}\right)$ and
$\left(\zeta_{j,+},\zeta_{j,-}\right)=\left(-\zeta_{yx,j},\zeta_{xy,j}\right)$.
We focus on the governing equation obtained from one block of Eq. (51), which
can be written as

$\displaystyle\frac{d}{dz}\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{y}\left(z\right)\end{array}\right)=i\frac{\omega}{c}A\left(z\right)\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{x}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{y}\left(z\right)\end{array}\right),$ (77)

where $A\left(z\right)$ is represented by

$\displaystyle
A\left(z\right)=\left(\begin{array}[]{cc}\zeta_{yx,j}&\mu_{j}\vspace{5pt}\\\
\varepsilon_{j}&\xi_{xy,j}\end{array}\right),$ (80)

when $z$ is in the layer $j$. We note that Eq. (77) and the other block of Eq.
(51) describe the linearly polarized eigenmodes whose polarization vectors lie
in the $x$ and $y$ axis, respectively. Clearly, these eigenmodes are
independent of each other.

As shown in Sec. III.2, the generalized Brillouin zone can be obtained by
$re^{i{\rm Re}\left(k\right)a}$, where $r$ is given by

$r=\exp\left[-\frac{1}{2c}{\rm
Im}\left(\omega\sum_{j=1}^{2}\left(\xi_{xy,j}+\zeta_{yx,j}\right)a_{j}\right)\right].$
(81)

We note that Eq. (81) reproduces the result obtained in Ref. [39] when
$\xi_{xy,j}=\zeta_{yx,j}$. One can infer from Eq. (81) that not only the
chirality parameters but also the eigenfrequencies contributes to the
localization of the non-Hermitian skin effect. Meanwhile, the non-Hermitian
skin effect disappears when the system becomes reciprocal, i.e.,
$\hat{\xi}=-\hat{\zeta}^{\rm T}$.

Equation (81) allows us to determine the generalized Brillouin zones [see Fig.
2(b)], which deviate from the conventional Brillouin zone. We then calculate
the continuum bands from the generalized Brillouin zones [red lines in Fig.
2(a)] and confirm that the continuum bands reproduce the eigenvalues of a
finite-size system [red dots in Fig. 2(c)]. We note that the red continuum
bands are distinct from the ones obtained from the real-valued Bloch
wavenumber [black curves in Fig. 2(a)]. Finally, Fig. 2(d) shows that the
eigenmodes are localized due to the non-Hermitian skin effect, and this
localization occurs at only one side of the system.

Figure 2: Non-Hermitian skin effect of the linearly polarized eigenmodes. (a)
Continuum bands calculated from the complex-valued Bloch wavenumber (red) and
real-valued Bloch wavenumber (black). (b) Generalized Brillouin zones (color
loops) and conventional Brillouin zone (black dashed circle). (c) Eigenvalues
of a finite-size system under the perfect electric conductor (PEC) boundary
condition (red) and a periodic boundary condition (black). The system size is
set to be $20a$. (d) Spatial profiles of all the eigenmodes included in the
band 1. We set the parameters as
$\varepsilon_{1}=\mu_{1}=1,\xi_{xy,1}=\zeta_{yz,1}=0,\varepsilon_{2}=4,\mu_{2}=1,\xi_{xy,2}=-0.3i,\zeta_{yx,2}=0.1$,
and $a_{1}=0.75a$.

###  III.4 Circularly polarized eigenmodes

We next investigate the case where the eigenmodes of the photonic crystal are
circularly polarized. The chirality tensors are given by

$\hat{\xi}_{j}=\xi\sigma_{0},~{}\hat{\zeta}_{j}=\zeta_{j}\sigma_{0},~{}\left(j=1,2\right),$
(82)

for which the eigenvectors of $i\hat{\xi}_{j}\sigma_{y}$ and
$-i\sigma_{y}\hat{\zeta}_{j}$ are proportional to $\left(1,\pm i\right)^{\rm
T}$. The corresponding eigenvalues are given by
$\left(\xi_{j,+},\xi_{j,-}\right)=\left(-i\xi_{j},i\xi_{j}\right)$ and
$\left(\zeta_{j,+},\zeta_{j,-}\right)=\left(i\zeta_{j},-i\zeta_{j}\right)$.
The governing equations then read

$\displaystyle\frac{d}{dz}\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{\pm}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{\pm}\left(z\right)\end{array}\right)=i\frac{\omega}{c}A_{\pm}\left(z\right)\left(\begin{array}[]{c}\sqrt{\varepsilon_{0}}E_{\pm}\left(z\right)\vspace{5pt}\\\
\sqrt{\mu_{0}}H_{\pm}\left(z\right)\end{array}\right),$ (87)

where $A_{\pm}\left(z\right)$ is represented by

$\displaystyle A_{\pm}\left(z\right)=\left(\begin{array}[]{cc}\mp
i\zeta_{j}&i\mu_{j}\vspace{5pt}\\\ -i\varepsilon_{j}&\pm
i\xi_{j}\end{array}\right),$ (90)

when $z$ is in the layer $j$. We note that
$\left(E_{+}\left(z\right),H_{+}\left(z\right)\right)$ and
$\left(E_{-}\left(z\right),H_{-}\left(z\right)\right)$ correspond to the left-
and right-circularly polarized eigenmodes, respectively. In contrast to the
linearly polarized case discussed in Sec. III.3, these circularly polarized
eigenmodes are not independent of each other, which means that the
corresponding eigenvalues are degenerate. Equation (67) allows us to calculate
the generalized Brillouin zones of the left- $\left(+\right)$ and right-
$\left(-\right)$ circularly polarized eigenmodes, where $r_{\pm}$ reads

$r_{+}=\exp\left[-\frac{1}{2c}{\rm
Re}\left(\omega\sum_{j=1}^{2}\left(\xi_{j}-\zeta_{j}\right)a_{j}\right)\right],~{}r_{-}=\frac{1}{r_{+}}.$
(91)

We show the generalized Brillouin zones of the left- and right-circularly
polarized eigenmodes in Figs. 3(b1) and (b2), respectively, and the
corresponding continuum bands are shown in Fig. 3(a). One can confirm that the
eigenvalues of the two circularly polarized eigenmodes are degenerate, while
the generalized Brillouin zones are distinct from each other. Figures. 3(a)
and (c) indicate that the continuum bands indeed reproduce the eigenvalues of
a finite-size system. Furthermore, as shown in Fig. 3(d), the left-circularly
polarized eigenmode (red) is localized at the right boundary, while the right-
circularly polarized eigenmode (blue) is localized at the opposite boundary.
We emphasize that these localization behaviors are in a stark contrast to the
case where the medium only has linearly polarized eigenmodes; in the latter,
all the eigenmodes localize to the edge of the same side. Moreover, in the
case of the circularly polarized eigenmodes, the non-Hermitian skin effect
persists even if the system is reciprocal, i.e., $\hat{\xi}=-\hat{\zeta}^{\rm
T}$. We speculate that these qualitative differences originate from the
underlying symmetries [7].

Figure 3: Non-Hermitian skin effect of the circularly polarized eigenmodes.
(a) Continuum bands calculated from the complex-valued Bloch wavenumber (red)
and real-valued Bloch wavenumber (black). (b) Generalized Brillouin zones
(color loops) and conventional Brillouin zone (black dashed loop). The
generalized Brillouin zones of the left- and right-circularly polarized
eigenmodes are shown in (b1) and (b2), respectively. (c) Eigenvalues of a
finite-size system under the perfect electric conductor (PEC) boundary
condition (red) and a periodic boundary condition (black). The system size is
set to be $20a$. (d) Spatial profiles of the left- (red) and right- (blue)
circularly polarized eigenmodes with $\omega a/2\pi c=0.2523$. We set the
parameters as
$\varepsilon_{1}=\mu_{1}=1,\xi_{xy,1}=\zeta_{yz,1}=0,\varepsilon_{2}=4,\mu_{2}=1,\xi_{2}=0.1,\zeta_{2}=-0.1$,
and $a_{1}=0.75a$.

##  IV Summary and Discussions

In this work, we develop the non-Bloch band theory of non-Hermitian systems
described by the generalized eigenvalue equation. One can calculate the
generalized Brillouin zone from the eigenvalues of the transfer matrix.
Notably, the theory allows us to investigate non-Hermitian physics of
multicomponent wavefunctions exhibiting the non-Hermitian skin effect, in
contrast to Refs. [25, 26] which have studied the non-Hermitian skin effect of
single-component wavefunctions. Indeed, we demonstrate that the theory is
applicable to electromagnetic waves in the photonic crystals composed of the
chiral metamaterials. As a result, we find that the localization properties of
the skin modes can qualitatively change depending on the polarization of the
eigenmodes. Notably, when the medium has circularly polarized eigenmodes, the
skin modes can localize at either of edges depending on whether it is left- or
right-circularly polarized. While we have focused on the case where the medium
has linearly or circularly polarized eigenmodes in this work, one can study
the non-Hermitian skin effect of elliptically polarized eigenmodes by setting
the chirality tensors to be

$\displaystyle\hat{\xi}_{j}=\left(\begin{array}[]{cc}\xi_{xx,j}&0\vspace{5pt}\\\
0&\xi_{yy,j}\end{array}\right),~{}\hat{\zeta}_{j}=\left(\begin{array}[]{cc}\zeta_{xx,j}&0\vspace{5pt}\\\
0&\zeta_{yy,j}\end{array}\right)$ (96)

satisfying $\xi_{xx,j}\zeta_{xx,j}=\xi_{yy,j}\zeta_{yy,j}$ with $j=1,2$.

It merits further study to analyze how the eigenmodes of the photonic crystal
with open boundary conditions can be actually excited by incident light. We
point out that the results revealed in this work suggest the possibility of
controlling the localization properties of excitation modes simply by changing
a polarization of incident light while keeping the medium unchanged. For
example, in the photonic crystal investigated in Sec. III.4, we expect that a
modulation of incident light from a linear polarization to a circular
polarization would change the localization positions of the excitation modes
from both ends to one end in steady-state regime. Meanwhile, it is also
intriguing to study dynamical behaviors of those excited modes in transient
regimes.

###### Acknowledgements.

We are grateful to Shuichi Murakami, Ryo Okugawa, Ryo Takahashi, and Tsuneya
Yoshida for variable discussion. K. Y. was supported by JSPS KAKENHI (Grants
No. JP21J01409 and No. JP23K13027). Y. A. acknowledges support from the Japan
Society for the Promotion of Science (JSPS) through Grant No. JP19K23424 and
from JST FOREST Program (Grant No. JPMJFR222U, Japan) and JST CREST (Grant No.
JPMJCR2312, Japan).

##  Appendix A Transfer matrix in tight-binding models

In this appendix, we discuss a method to construct a transfer matrix in 1D
non-Hermitian tight-binding systems on the lines of Ref. [42]. We derive the
condition for the generalized Brillouin zone in terms of the eigenvalues of
the transfer matrix. One can indeed see that the band structure calculated
from the generalized Brillouin zone reproduces the energy eigenvalues under
open boundary conditions.

###  A.1 Setup

We start from 1D non-Hermitian tight-binding systems which are described by

$H=\sum_{i}\sum_{j=-N}^{N}\sum_{\mu,\nu=1}^{q}t_{j,\mu\nu}c_{i+j,\mu}^{\dagger}c_{i,\nu},$
(97)

where $c_{i,\nu}$ is an annihilation operator of a particle at sublattice
$\nu$ in the $i$th unit cell, and $t_{j,\mu\nu}$ is not equal to
$t_{-j,\nu\mu}^{\ast}$. We assume that the particles hop to the $N$th nearest
unit cell. In Eq. (97), sublattice sites in a unit cell are connected to up to
the $2N$th-nearest-neighbor unit cells via the hopping amplitudes. Meanwhile,
we can enlarge a unit cell so that sublattice sites in a given cell are
connected to only the nearest-neighbor cells [41, 42]. The enlarged unit cell
is called a supercell including $s\geq qN$ degrees of freedom, the definition
of which is not unique. Let ${\bm{c}}_{n}$ denote a vector of the annihilation
operators of the particles in the $n$th supercell. Equation (97) can be then
rewritten as the reduced form including only the nearest-neighbor hopping
amplitudes:

$H=\sum_{n}\left({\bm{c}}^{\dagger}_{n}J_{L}{\bm{c}}_{n+1}+{\bm{c}}_{n}^{\dagger}M{\bm{c}}_{n}+{\bm{c}}_{n+1}^{\dagger}J_{R}^{\dagger}{\bm{c}}_{n}\right),$
(98)

where $J_{L}$ and $J_{R}$ are hopping matrices, and $M$ is an onsite matrix.
The reduced form of the real-space eigenequation reads

$J_{L}{\bm{\Psi}}_{n+1}+M{\bm{\Psi}}_{n}+J_{R}^{\dagger}{\bm{\Psi}}_{n+1}=E{\bm{\Psi}}_{n},$
(99)

where ${\bm{\Psi}}_{n}$ is a wavefunction for the $n$th supercell. We note
that the systems become Hermitian when $J_{L}=L_{R}$ and $M$ is a Hermitian
matrix.

In the following, we focus on the systems with the hopping matrices satisfying

$J_{L}^{2}=J_{R}^{2}=O$ (100)

and

${\rm rank}\left(J_{L}\right)={\rm rank}\left(J_{R}\right)=r,$ (101)

for the sake of simplicity. We note that taking a sufficiently large supercell
ensures Eq. (100). Indeed, if $J_{L}=J_{R}$, the physical interpretation of
Eq. (100) is that, in a given supercell, no sublattice sites are connected to
both the left and right adjacent supercells. Meanwhile, Eq. (101) means that,
in a given supercell, $r$ sublattice sites are connected to the adjacent
supercell. In this situation, we derive the transfer matrix of the systems
from Eq. (99). To this end, we utilize the singular value decompositions of
the hopping matrices given by

$\displaystyle\left\\{\begin{array}[]{l}J_{L}=X\Xi_{L}Y^{\dagger},\vspace{5pt}\\\
J_{R}=V\Xi_{R}W^{\dagger},\end{array}\right.$ (104)

where $\Xi_{L}$ and $\Xi_{R}$ are $r\times r$ diagonal matrices including the
singular values of $J_{L}$ and $J_{R}$, respectively, and $X,Y,V$, and $W$ are
$s\times r$ matrices satisfying

$X^{\dagger}X=Y^{\dagger}Y=V^{\dagger}V=W^{\dagger}W=\mbox{1}\mbox{l},$ (105)

and

$X^{\dagger}Y=V^{\dagger}W=O.$ (106)

It is worth noting that Eq. (106) is always ensured because of Eq. (100). We
can then rewrite Eq. (99) as

${\bm{\Psi}}_{n}=GX\Xi_{L}{\bm{\beta}}_{n+1}+GW\Xi_{R}{\bm{\alpha}}_{n-1},$
(107)

where ${\bm{\alpha}}_{n}=V^{\dagger}{\bm{\Psi}}_{n}$,
${\bm{\beta}}_{n}=Y^{\dagger}{\bm{\Psi}}_{n}$, and $G=\left(E-M\right)^{-1}$.
For the sake of convenience, let $G_{ab}$ denote $B^{\dagger}GA$ for $s\times
r$ matrices $A,B$. Equation. (107) leads to the recursion equation for
${\bm{\alpha}}_{n}$ and ${\bm{\beta}}_{n}$ given by

$\displaystyle\left(\begin{array}[]{c}{\bm{\beta}}_{n+1}\vspace{5pt}\\\
{\bm{\alpha}}_{n}\end{array}\right)=T\left(\begin{array}[]{c}{\bm{\beta}}_{n}\vspace{5pt}\\\
{\bm{\alpha}}_{n-1}\end{array}\right),$ (112)

where $T$ is the transfer matrix explicitly written as

$\displaystyle
T=\left(\begin{array}[]{cc}\Xi_{L}^{-1}G_{xy}^{-1}&-\Xi_{L}G_{xy}^{-1}G_{wy}\Xi_{R}\vspace{5pt}\\\
G_{xv}G_{xy}^{-1}&\left(G_{wv}-G_{xv}G_{xy}^{-1}G_{wy}\right)\Xi_{R}\end{array}\right).$
(115)

We note that the size of the transfer matrix is independent of the choice of a
supercell.

###  A.2 Condition for the generalized Brillouin zone

We show that the condition for the generalized Brillouin zone can be obtained
from the eigenvalues of the transfer matrix. For the sake of simplicity, we
suppose that the transfer matrix is a diagonalizable matrix, for which the
eigenvalue problem is written as

$T{\bm{\varphi}}^{\left(l\right)}=\rho_{l}{\bm{\varphi}}^{\left(l\right)}$
(116)

with $l=1,\cdots,2r$, where we number the eigenvalues so as to satisfy
$\left|\rho_{1}\right|\leq\cdots\leq\left|\rho_{2r}\right|$. In the following,
We focus on a finite-size system with open boundary conditions expressed by

${\bm{\Psi}}_{0}={\bm{\Psi}}_{L+1}=0,$ (117)

where $L$ is a system size. One can immediately obtain the boundary equation
from Eq. (112) as follows:

$\displaystyle\left(\begin{array}[]{c}{\bm{0}\vspace{5pt}}\\\
{\bm{\alpha}}_{L}\end{array}\right)=T^{L}\left(\begin{array}[]{c}{\bm{\beta}}_{1}\vspace{5pt}\\\
{\bm{0}}\end{array}\right).$ (122)

We note that one can take ${\bm{\alpha}}_{L}$ and ${\bm{\beta}}_{1}$ to be
arbitrary vectors. To solve Eq. (122), we expand the vectors included in this
equation in terms of the eigenvectors of the transfer matrix, which leads to

$\displaystyle\left(\begin{array}[]{c}{\bm{\beta}}_{1}\vspace{5pt}\\\
{\bm{0}}\end{array}\right)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{2r}a_{l}{\bm{\varphi}}^{\left(l\right)},$ (125)
$\displaystyle\left(\begin{array}[]{c}{\bm{0}\vspace{5pt}}\\\
{\bm{\alpha}}_{L}\end{array}\right)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{2r}a_{l}\left(\rho_{l}\right)^{L}{\bm{\varphi}}^{\left(l\right)}.$
(128)

Let ${\cal P}_{\alpha}=\left(O,\mbox{1}\mbox{l}\right)$ and ${\cal
P}_{\beta}=\left(\mbox{1}\mbox{l},O\right)$ denote $r\times 2r$ matrices.
Acting ${\cal P}_{\alpha}$ and ${\cal P}_{\beta}$ on Eqs. (125) and (128),
respectively, we get

$\displaystyle\begin{array}[]{l}\displaystyle\sum_{l=1}^{2r}a_{l}{\cal
P}_{\alpha}{\bm{\varphi}}^{\left(l\right)}={\bm{0}},\vspace{5pt}\\\
\displaystyle\sum_{l=1}^{2r}a_{l}\left(\rho_{l}\right)^{L}{\cal
P}_{\beta}{\bm{\varphi}}^{\left(l\right)}={\bm{0}}.\end{array}$ (131)

Equation (131) is a set of algebraic equations for $a_{1},\cdots,a_{2r}$, and
we can recast it to the form of a matrix equation as follows:

$\displaystyle\left(R_{1}{\bm{\varphi}}^{\left(1\right)}~{}\cdots~{}R_{2r}{\bm{\varphi}}^{\left(2r\right)}\right)\left(\begin{array}[]{c}a_{1}\vspace{5pt}\\\
\vdots\vspace{5pt}\\\ a_{2r}\end{array}\right)={\bm{0}},$ (135)

where $R_{l}$ are $2r\times 2r$ matrices defined by

$\displaystyle
R_{l}=\left(\begin{array}[]{cc}\left(\rho_{l}\right)^{L}\mbox{1}\mbox{l}&O\vspace{5pt}\\\
O&\mbox{1}\mbox{l}\end{array}\right).$ (138)

The condition that Eq. (135) has a nontrivial solution finally reads

$\displaystyle\left|\begin{array}[]{ccc}\varphi_{1}^{\left(1\right)}\left(\rho_{1}\right)^{L}&\cdots&\varphi_{1}^{\left(2r\right)}\left(\rho_{2r}\right)^{L}\vspace{5pt}\\\
\vdots&\vdots&\vdots\vspace{5pt}\\\
\varphi_{r}^{\left(1\right)}\left(\rho_{1}\right)^{L}&\cdots&\varphi_{r}^{\left(2r\right)}\left(\rho_{2r}\right)^{L}\vspace{5pt}\\\
\varphi_{r+1}^{\left(1\right)}&\cdots&\varphi_{r+1}^{\left(2r\right)}\vspace{5pt}\\\
\vdots&\vdots&\vdots\vspace{5pt}\\\
\varphi_{2r}^{\left(1\right)}&\cdots&\varphi_{2r}^{\left(2r\right)}\end{array}\right|=0.$
(145)

We remark that one can get the eigenenergies of a finite-size system with
edges by solving Eq. (145), since this equation is defined only by the
eigenvalues and eigenvectors of the transfer matrix.

To study the condition for the generalized Brillouin zone, we expand Eq. (145)
as follows:

$\sum_{P,P^{\prime}}\sum_{Q,Q^{\prime}}F\left(\varphi_{i\in Q}^{\left(j\in
P\right)},\varphi_{i^{\prime}\in Q^{\prime}}^{\left(j^{\prime}\in
P^{\prime}\right)}\right)\prod_{i\in P}\left(\rho_{i}\right)^{L}=0,$ (146)

where $P$ and $P^{\prime}$ denote two disjoint subsets of the set
$\left\\{1,\cdots,2r\right\\}$ with $r$ elements,
$Q=\left\\{1,\cdots,r\right\\}$, and
$Q^{\prime}=\left\\{r+1,\cdots,2r\right\\}$. While it is difficult to get the
exact solutions of Eq. (146) in general, one can investigate the asymptotic
behavior of the eigenenergies in $L\rightarrow\infty$. The key observation is
that, when

$\left|\rho_{r}\right|=\left|\rho_{r+1}\right|$ (147)

is satisfied, the dominant contributions to the left-hand side of Eq. (146)
are the term including $\left(\rho_{r}\rho_{r+2}\cdots\rho_{2r}\right)^{L}$
and the one including $\left(\rho_{r+1}\rho_{r+2}\cdots\rho_{2r}\right)^{L}$
in the limit of a large system size. As a result, the asymptotic form of Eq.
(146) can be obtained by

$\left(\frac{\rho_{r}}{\rho_{r+1}}\right)^{L}=-\frac{F\left(\varphi_{i\in
Q}^{\left(j\in P_{0}\right)},\varphi_{i^{\prime}\in
Q^{\prime}}^{\left(j^{\prime}\in
P_{0}^{\prime}\right)}\right)}{F\left(\varphi_{i\in Q}^{\left(j\in
P_{1}\right)},\varphi_{i^{\prime}\in Q^{\prime}}^{\left(j^{\prime}\in
P_{1}^{\prime}\right)}\right)},$ (148)

where
$P_{0}=\left\\{r+1,\cdots,2r\right\\},P_{0}^{\prime}=\left\\{1,\cdots,r\right\\},P_{1}=\left\\{r,r+2,\cdots,2r\right\\}$,
and $P_{1}^{\prime}=\left\\{1,\cdots,r-1,r\right\\}$. Thus, one can get the
continuum energy bands, since the change of the relative phase between
$\rho_{r}$ and $\rho_{r+1}$ produces dense sets of the eigenenergies.
Accordingly, the sets of the eigenvalues of the transfer matrix satisfying Eq.
(147) form the generalized Brillouin zone.

###  A.3 Non-Hermitian Su-Schrieffer-Heeger model

We show that Eq. (147) allows us to calculate the generalized Brillouin zone
and the continuum energy bands by using the non-Hermitian Su-Schrieffer-Heeger
model with asymmetric next-nearest-neighbor hopping amplitudes [see Fig.
4(a)]. The real-space Hamiltonian of this system reads

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{n}\left[t_{1}\left(c_{n,\alpha}^{\dagger}c_{n,\beta}+c_{n,\beta}^{\dagger}c_{n,\alpha}\right)\right.$
$\displaystyle+t_{2}\left(c_{n,\beta}^{\dagger}c_{n+1,\alpha}+c_{n+1,\alpha}^{\dagger}c_{n,\beta}\right)$
$\displaystyle\left.+\left(t_{3}+\frac{\gamma}{2}\right)c_{n,\alpha}^{\dagger}c_{n+1,\beta}+\left(t_{3}-\frac{\gamma}{2}\right)c_{n+1,\beta}^{\dagger}c_{n,\alpha}\right],$

where all the parameters take positive values, and we assume $t_{3}>\gamma/2$
for simplicity. We note that Eq. (LABEL:eqappa24) has the form of Eq. (97)
with $N=1$ and $q=2$.

To calculate the transfer matrix of the system in the manner explained in Sec.
A.1, we take a supercell including four sublattices [see Fig. 4(b)]. One can
then obtain the reduced Hamiltonian in the form of Eq. (98) and also the
reduced real-space eigenequation for the wavefunction
${\bm{\Psi}}_{n}=\left(\Psi_{n,{\rm A}},\Psi_{n,{\rm B}},\Psi_{n,{\rm
C}},\Psi_{n,{\rm D}}\right)$ in the form of Eq. (99). Here, the hopping
matrices and onsite matrix are given by

$\displaystyle J_{R}^{\dagger}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&0&0&t_{2}\vspace{5pt}\\\
0&0&t_{3}-\gamma/2&0\vspace{5pt}\\\ 0&0&0&0\vspace{5pt}\\\
0&0&0&0\end{array}\right),$ (154) $\displaystyle J_{L}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&0&0&0\vspace{5pt}\\\
0&0&0&0\vspace{5pt}\\\ 0&t_{3}+\gamma/2&0&0\vspace{5pt}\\\
t_{2}&0&0&0\end{array}\right),$ (159) $\displaystyle M$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&t_{1}&0&t_{3}+\gamma/2\vspace{5pt}\\\
t_{1}&0&t_{2}&0\vspace{5pt}\\\ 0&t_{2}&0&t_{1}\vspace{5pt}\\\
t_{3}-\gamma/2&0&t_{1}&0\end{array}\right).$ (164)

We note that the rank of both the hopping matrices is $2$, which means that
the size of the transfer matrix becomes $4$. Let
$\rho_{r}~{}\left(r=1,\dots,4\right)$ denote the eigenvalues of the transfer
matrix. The condition for the generalized Brillouin zone can be obtained as
follows:

$\left|\rho_{2}\right|=\left|\rho_{3}\right|,$ (165)

where $\left|\rho_{1}\right|\leq\cdots\leq\left|\rho_{4}\right|$ is satisfied.

We calculate the generalized Brilloun zone [Fig. 4(c)] and the continuum
energy bands [Fig. 4(d)] based on Eq. (165). We indeed confirm that the
continuum energy bands reproduces the eigenenergies of a finite-size system
with open boundary conditions.

Figure 4: Schematic figures, generalized Brillouin zones, and eigenenergies of
the non-Hermitian Su-Schrieffer-Heeger model. (a) Unit cells represented by
the dashed boxes. (b) Supercell of our choice represented by the box. (c)
Generalized Brillouin zones. (d) Continuum energy bands (red) and
eigenenergies of a finite-size system with open boundary conditions (black).
The system size is set to be $L=50$. We set the system parameters as
$t_{1}=1,t_{2}=0.7,t_{3}=1.2$, and $\gamma=4/3$.

## References

  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Kunst  _et al._ [2018] F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Phys. Rev. Lett. 121, 026808 (2018).
  * Yao and Wang [2018] S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Borgnia  _et al._ [2020] D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Phys. Rev. Lett. 124, 056802 (2020).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * Zhang  _et al._ [2020] K. Zhang, Z. Yang, and C. Fang, Phys. Rev. Lett. 125, 126402 (2020).
  * Yi and Yang [2020] Y. Yi and Z. Yang, Phys. Rev. Lett. 125, 186802 (2020).
  * Brandenbourger  _et al._ [2019] M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Nat. Commun. 10, 4608 (2019).
  * Xiao  _et al._ [2020] L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Nat. Phys. 16, 761 (2020).
  * Weidemann  _et al._ [2020] S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Science 368, 311 (2020).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Nat. Phys. 16, 747 (2020).
  * Hofmann  _et al._ [2020] T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši, C. H. Lee, A. Bilušić, R. Thomale, and T. Neupert, Phys. Rev. Research 2, 023265 (2020).
  * Ghatak  _et al._ [2020] A. Ghatak, M. Brandenbourger, J. van Wezel, and C. Coulais, Proc. Nat. Ac. Sc. USA 117, 29561 (2020).
  * Chen  _et al._ [2021] Y. Chen, X. Li, C. Scheibner, V. Vitelli, and G. Huang, Nat. Commun. 12, 5935 (2021).
  * Zhang  _et al._ [2021] L. Zhang, Y. Yang, Y. Ge, Y.-J. Guan, Q. Chen, Q. Yan, F. Chen, R. Xi, Y. Li, D. Jia, _et al._ , Nat. Commun. 12, 6297 (2021).
  * Wang  _et al._ [2022] W. Wang, X. Wang, and G. Ma, Nature 608, 50 (2022).
  * Liang  _et al._ [2022] Q. Liang, D. Xie, Z. Dong, H. Li, H. Li, B. Gadway, W. Yi, and B. Yan, Phys. Rev. Lett. 129, 070401 (2022).
  * Gu  _et al._ [2022] Z. Gu, H. Gao, H. Xue, J. Li, Z. Su, and J. Zhu, Nature Communications 13, 7668 (2022).
  * Yokomizo and Murakami [2019] K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123, 066404 (2019).
  * Yokomizo and Murakami [2020] K. Yokomizo and S. Murakami, Prog. Theor. Exp. Phys. 2020, 12A102 (2020).
  * Kawabata  _et al._ [2020] K. Kawabata, N. Okuma, and M. Sato, Phys. Rev. B 101, 195147 (2020).
  * Yang  _et al._ [2020] Z. Yang, K. Zhang, C. Fang, and J. Hu, Phys. Rev. Lett. 125, 226402 (2020).
  * Yokomizo  _et al._ [2022] K. Yokomizo, T. Yoda, and S. Murakami, Phys. Rev. Research 4, 023089 (2022).
  * [26] Y.-M. Hu, Y.-Q. Huang, W.-T. Xue, and Z. Wang, arXiv:2310.08572 .
  * Inoue and Murakami [2019] T. Inoue and S. Murakami, Phys. Rev. B 99, 195443 (2019).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, Phys. Rev. Lett. 100, 013904 (2008).
  * Raman and Fan [2010] A. Raman and S. Fan, Phys. Rev. Lett. 104, 087401 (2010).
  * Shindou  _et al._ [2013a] R. Shindou, J.-i. Ohe, R. Matsumoto, S. Murakami, and E. Saitoh, Phys. Rev. B 87, 174402 (2013a).
  * Shindou  _et al._ [2013b] R. Shindou, R. Matsumoto, S. Murakami, and J.-i. Ohe, Phys. Rev. B 87, 174427 (2013b).
  * [32] T. Isobe, T. Yoshida, and Y. Hatsugai, arXiv:2310.12577 .
  * Isobe  _et al._ [2021] T. Isobe, T. Yoshida, and Y. Hatsugai, Phys. Rev. B 104, L121105 (2021).
  * Isobe  _et al._ [2023] T. Isobe, T. Yoshida, and Y. Hatsugai, Nanophotonics 12, 2273 (2023).
  * McDonald  _et al._ [2018] A. McDonald, T. Pereg-Barnea, and A. A. Clerk, Phys. Rev. X 8, 041031 (2018).
  * Yokomizo and Murakami [2021] K. Yokomizo and S. Murakami, Phys. Rev. B 103, 165123 (2021).
  * Rosa and Ruzzene [2020] M. I. Rosa and M. Ruzzene, New J. Phys. 22, 053004 (2020).
  * Braghini  _et al._ [2021] D. Braghini, L. G. Villani, M. I. Rosa, and J. R. de F Arruda, J. Phys. D 54, 285302 (2021).
  * Yan  _et al._ [2021] Q. Yan, H. Chen, and Y. Yang, Prog. Electromagn. Res. 172, 33 (2021).
  * [40] T. Yoda, Y. Moritake, K. Takata, K. Yokomizo, S. Murakami, and M. Notomi, arXiv:2303.05185 .
  * Dwivedi and Chua [2016] V. Dwivedi and V. Chua, Phys. Rev. B 93, 134304 (2016).
  * Kunst and Dwivedi [2019] F. K. Kunst and V. Dwivedi, Phys. Rev. B 99, 245116 (2019).


This concludes the construction of the decoder.

To show correctness, we need to show $\tilde c_1 \in \hat c_1 + B_1$.

We now show

\begin{equation*}
\tilde c_1 = \hat c_1 + \partial_2 (\hat s'_2(V_{00}) + \hat w_2(V_{00}))
\end{equation*}

by checking the equality on $\tilde c_1(E_{*0})$, $\tilde c_1(E_{0*})$, $\tilde c_1(E_{*1})$, $\tilde c_1(E_{1*})$ individually.

Check $\tilde c_1(E_{*0})$.

\begin{equation*}
\hat c_1(E_{*0}) + H_B^\leftarrow \hat s'_2(V_{00}) + H_B^\leftarrow \hat w_2(V_{00})
= \hat c_1(V_{00}, E_{*0}) + H_B^\leftarrow \hat s'_2(V_{00}) + 0
= s_1(V_{00}, E_{*0})
= \tilde c_1(E_{*0}).
\end{equation*}

Check $\tilde c_1(E_{0*})$.

\begin{equation*}
\hat c_1(E_{0*}) + H_A^\uparrow \hat s'_2(V_{00}) + H_A^\uparrow \hat w_2(V_{00})
= \hat c_1(V_{00}, E_{0*}) + H_A^\uparrow \hat s'_2(V_{00}) + 0
= s_1(V_{00}, E_{0*})
= \tilde c_1(E_{0*}).
\end{equation*}

Check $\tilde c_1(E_{*1})$.

\begin{align*}
\hat c_1(E_{*1}) + H_B^\rightarrow \hat s'_2(V_{00}) + H_B^\rightarrow \hat w_2(V_{00})
&= \hat c_1(V_{01}, E_{*1}) + H_B^\rightarrow \hat s'_2(V_{00}) + H_B^\rightarrow \hat w_2(V_{00}) \\
&= s_1(V_{01}, E_{*1}) + H_B^\rightarrow \hat s'_2(V_{01}) + H_B^\rightarrow \hat s'_2(V_{00}) + H_B^\rightarrow \hat w_2(V_{00}) \\
&= s_1(V_{01}, E_{*1}) + H_B^\rightarrow \hat t'_2(E_{0*}) + H_B^\rightarrow \hat w_2(V_{00}) \\
&= s_1(V_{01}, E_{*1}) + H_B^\rightarrow \tilde t'_2(E_{0*}) + H_B^\rightarrow \hat w_2(V_{01}) \\
&= s_1(V_{01}, E_{*1}) + H_B^\rightarrow \tilde t'_2(E_{0*}) + 0 \\
&= \tilde c_1(E_{*1}).
\end{align*}

Check $\tilde c_1(E_{1*})$.

\begin{align*}
\hat c_1(E_{1*}) + H_A^\downarrow \hat s'_2(V_{00}) + H_A^\downarrow \hat w_2(V_{00})
&= \hat c_1(V_{10}, E_{1*}) + H_A^\downarrow \hat s'_2(V_{00}) + H_A^\downarrow \hat w_2(V_{00}) \\
&= s_1(V_{10}, E_{1*}) + H_A^\downarrow \hat s'_2(V_{10}) + H_A^\downarrow \hat s'_2(V_{00}) + H_A^\downarrow \hat w_2(V_{00}) \\
&= s_1(V_{10}, E_{1*}) + H_A^\downarrow \hat t'_2(E_{*0}) + H_A^\downarrow \hat w_2(V_{00}) \\
&= s_1(V_{10}, E_{1*}) + H_A^\downarrow \tilde t'_2(E_{*0}) + H_A^\downarrow \hat w_2(V_{10}) \\
&= s_1(V_{10}, E_{1*}) + H_A^\downarrow \tilde t'_2(E_{*0}) + 0 \\
&= \tilde c_1(E_{1*}).
\end{align*}

Linear Time:

It is clear that the construction at each step takes linear time.
For the steps besides Step4 where we apply the co-decoder assumption,
we only use simple operations such as linear map or inverse of a linear map.
And Step4 takes linear time by assumption.

§ OPTIMAL ROBUST TENSOR CODES

This section shows that random codes have linear robustness with high probability.
We improve on the result in [11, 12]
by using the idea of puncturing and a new counting argument.

Recall that $C_A$ and $C_B$ are linear codes of length $n_a$ and $n_b$ with rate $\rho_a$ and $\rho_b$. For simplicity we assume that $n_a=n_b=\Delta$.
An $s$-punctured code of $C_A$ is obtained  by first choosing $s$ coordinates $I_a \subset [n_a]$ then consider the codewords of $C_A$ restricted to $[n_a]-I_a$.
Notice that $\norm{c}_{[n_a] \times [n_b]} = \abs{c}$ is identical to the Hamming weight, so we will mostly use $\abs{c}$ in this section to simplify the notation.

We now recall <Ref>, which is the main theorem to prove in this section.

The proof follows a similar strategy in [11, 12], where the key is to show that a codeword with a small weight $|c| < w = \Theta(\Delta^2)$ is “structured”, i.e. $c$ is only supported on a few columns and a few rows, with high probability.

To show codewords with small weights are “structured”, we show all non-zero codewords in a random code have some column or row with large weight with high probability. Because the punctured code of a random code is still roughly a random code, the same property also applies to its punctured codes.
Now, since we assumed the codeword $c\in \tensorcode(C_A, C_B)$ has small weight, we can remove a few columns and rows with large weights, such that the rest have small weight in all columns and rows.
We then apply the property of the punctured code above which implies the rest is $0$, so $c$ is only supported on those removed columns and rows, i.e. $c$ is “structured”.

When $c$ is “structured”, one can then find $c_a$ supported on the few columns and $c_b$ supported on the few rows.
This means the cancellation in $c_a + c_b$ could only happen in the intersection of those columns and rows which is small.
Since each column of $c_a$ is a codeword, when the distance is large, $\abs{c_a} \ge d_1 \norm{c_a}_{[n_b]} = \Theta(\Delta) \norm{c_a}_{[n_b]}$.
This implies codewords with small weight satisfy the inequality for robustness $\abs{c} = \abs{c_a} + \abs{c_b} - \textnormal{ small number of cancellations } \ge \Theta(\Delta) (\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]})$.

When $c$ has large weight $|c| \ge w = \Theta(\Delta^2)$,
because $\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]} \le 2 \Delta$,
the inequality for robustness $\abs{c} \ge d_2 (\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]})$ is easily satisfied
by setting $d_2 = w/(2\Delta) = \Theta(\Delta)$.

We now go through the steps carefully below.
We first state two lemmas, show the theorem, then prove the lemmas.
The first lemma says each non-zero codeword has at least one row or column with large weight (which implies codewords with small weight are “structured”).
The second lemma says “structured” codewords satisfy robustness.

Fix $\rho_a, \rho_b \in (0,1)$, let $\sigma \in (0, 1), \tau \in (0, (1-\sigma)/2)$ satisfy

\begin{equation*}
2 h(\sigma) + 2 (1-\sigma)h(\frac{\tau}{1-\sigma}) < \frac{3}{4} \frac{(1 - \sigma - \rho_a)(1 - \sigma - \rho_b)}{1-\sigma}.
\end{equation*}

Let $C_A, C_B$ be random codes sampled from the uniform distribution with length $\Delta$ and dimensions $k_a=\rho_a \Delta, k_b=\rho_b \Delta$.
Then as $\Delta$ goes to infinity, with probability tending to $1$, the following holds:
for any $s = \sigma \Delta$-punctured code $C_A', C_B'$,
all non-zero codewords in $\tensorcode(C_A', C_B')$ have at least one row or one column with weight $\ge t = \tau \Delta$.
In other words, if a codeword in $\tensorcode(C_A', C_B')$ has all its rows and columns with weight $< t$, then the codeword is $0$.

Previously, we said the robustness parameter can be thought of as the higher dimensional generalization of linear distance.
Here, <Ref> provides an alternative generalization that could also be useful when trying to study higher dimensional tensor codes.
Note that in one-dimension, these two definitions coincide.

[Modification of <cit.> or <cit.>] lemmaRobustLemma
Suppose $C_A$ and $C_B$ have distance $d_1$.
If $c \in \tensorcode(C_A, C_B)$ is supported on $I_a \times [n_b] \cup [n_a] \times I_b$
and $|I_a|, |I_b| < d_1$,
then there exists $c_a \in C_A \otimes \F2^{n_b}$ supported on $[n_a] \times I_b$ and $c_b \in \F2^{n_a} \otimes C_B$ supported on $I_a \times [n_b]$
such that $c = c_a + c_b$.

Furthermore, if $|I_a|, |I_b| < d_1/2$, we have

\begin{equation*}
\abs{c} \ge \frac{d_1}{2}(\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]}).
\end{equation*}

We first show linear distance. The <Ref> in the theorem statement implies $h(\delta_1) < 1 - \rho_a$ and $h(\delta_1) < 1 - \rho_b$ because $h(\delta_1) < 2 h(\delta_1/2)$ for $\delta_1 < 1/2$.
By Gilbert-Varshamov bound, we know $C_A, C_B$ have distance $d_1 = \delta_1 \Delta$ with probability tending to $1$.

Next, we prove that $(C_A,C_B)$ is $\delta_2 \Delta$-robust when $C_A, C_B$ have distance $d_1$ and the condition in <Ref> holds.
More precisely, because the inequality in <Ref> holds for $\sigma = \delta_1/2$ and $\tau = 4 \delta_2/\delta_1$,
we have that for all $s = \sigma \Delta$-punctured codes $C_A', C_B'$, if a codeword in $\tensorcode(C_A', C_B')$ has all columns and rows with weight $< t = \tau \Delta$,
then the codeword is $0$.
Note that the same holds for all punctured codes $C_A', C_B'$ with puncture $\le s$.

Given $c \in \tensorcode(C_A, C_B)$, we consider the following two cases. First, assume that $c$ has small weight $|c| < d_1 t / 2$. We can remove $d_1/2-1=s-1$ columns $I_a$ and $s-1$ rows $I_b$ such that $c_{([n_a]-I_a)\times([n_b]-I_b)}$ has weight less than $t$ for each column and each row.
By <Ref>,
$c_{([n_a]-I_a)\times([n_b]-I_b)}$ has to be $0$.
So $c$ is supported on at most $s-1$ columns and $s-1$ rows.
Now, because $s - 1 < s = d_1/2$, by <Ref>

\begin{equation*}
\abs{c} \ge \frac{d_1}{2} (\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]}).
\end{equation*}

Second, suppose $c$ has large weight $|c| \ge d_1 t / 2$.
Because $c \in \tensorcode(C_A, C_B)$, we can always write $c = c_a + c_b$.
For any choice of $c_a, c_b$, we have $\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]} \le 2 \Delta$.
Together with the fact that $c$ has large weight, we have

\begin{equation*}
\abs{c} \ge \frac{d_1 t}{4 \Delta} (\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]}).
\end{equation*}

Therefore, the code pair $(C_A, C_B)$ is $\min(\frac{d_1}{2}, \frac{d_1 t}{4 \Delta}) = \frac{d_1 t}{4 \Delta} = \delta_2 \Delta$ robust.

We now prove the first lemma.
The proof is similar to the proof of the Gilbert-Varshamov bound where we count the number of bad events and take the union bound.
Here, we organize the bad events by the rank of the matrices.
We first count the number of bad matrices with rank $r$.
Then, we compute the probability for a rank $r$ matrix to be a codeword. Finally, we take the union bound.

We again first state the two claims, show the lemma, then prove the claims.

[Number of bad matrices is small]
Let $\cM(\Delta, r, t)$ be the set of matrices in $\F2^{\Delta \times \Delta}$ with rank $r$ where each row and each column has weight $< t$. Assume $t \le \Delta/2$.

\begin{equation*}
\abs{\cM(\Delta, r, t)} \le 2^{2r\Delta h(\frac{t}{\Delta}) + 2\Delta h(\frac{r}{\Delta})}.
\end{equation*}

[A bad matrix is likely not a codeword]
Let $c \in \F2^{\Delta \times \Delta}$ be a matrix of rank $r$.

\begin{equation*}
\cP(\Delta, r, k_a, k_b) \coloneqq \Pr[c \in \Sigma(C_A, C_B)] \le 512 (r+1) 2^{-\frac{3}{4} \frac{(\Delta-k_a) (\Delta-k_b)}{\Delta} r},
\end{equation*}

where $C_A, C_B$ are uniformly sampled from linear codes of length $\Delta$ and dimension $k_a, k_b$.

To show the lemma, we apply the union bound over the bad events.
First, there are $\binom{\Delta}{s}^2$ possible $s$-punctured code pairs $C_A', C_B'$ when puncturing both $C_A$ and $C_B$.
Next, there are $\abs{\cM(\Delta-s, r, t)}$ many rank $r$ matrices where all rows and columns have weight $<t$.
Finally, each of those rank $r$ matrix is a codeword with probability $\le \cP(\Delta-s, r, k_a, k_b)$ which we explain below.
Recall that $C_A$ is a uniform distribution over dimension $k_a$,
so $C_A'$ is a mixture of uniform distribution over dimensions $k'_a \le k_a$.
Suppose $C_A'$($C_B'$) has dimension $k'_a$($k'_b$) with probability $p_{k'_a}$($q_{k'_b}$),
then a rank $r$ matrix is a codeword with probability $\le \sum_{k'_a=0}^{k_a} \sum_{k'_b=0}^{k_b} p_{k'_a} q_{k'_b} \cP(\Delta-s, r, k'_a, k'_b)$.
Because smaller dimension leads to smaller probability of becoming a codeword, $\cP(\Delta-s, r, k'_a, k'_b) \le \cP(\Delta-s, r, k_a, k_b)$,
the probability for a rank $r$ matrix to be a codeword $\le \cP(\Delta-s, r, k_a, k_b)$ as claimed.

Therefore, the probability where the main condition fails to hold (there are some $s$-punctured code pair $C'_A, C'_B$ with non-zero codeword in $\tensorcode(C_A', C_B')$
where all its rows and columns has weight $< t$) is

\begin{align*}
&\le \binom{\Delta}{s}^2 \sum_{r=1}^{\Delta-s} \abs{\cM(\Delta-s, r, t)}\, \cP(\Delta-s, r, k_a, k_b) \\
&\le 2^{2 \Delta h(\frac{s}{\Delta})} \cdot
\sum_{r=1}^{\Delta-s} 2^{2r(\Delta-s) h(\frac{t}{\Delta-s}) + 2(\Delta-s) h(\frac{r}{\Delta-s})} \cdot
512 (r+1) 2^{-\frac{3}{4} \frac{(\Delta-s-k_a) (\Delta-s-k_b)}{\Delta-s} r}
\end{align*}

where the last inequality uses $\binom{\Delta}{s} \le 2^{\Delta h(\frac{s}{\Delta})}$ and the results from <Ref> and <Ref>.

Finally, we find the requirements for the value to approach $0$ as $\Delta$ goes to infinity.
This boils down to two conditions.
The first is

\begin{equation*}
2 \frac{\Delta - s}{\Delta} h(\frac{t}{\Delta-s}) < \frac{3}{4} \frac{(\Delta - s - k_a)(\Delta - s - k_b)}{\Delta (\Delta - s)}
\end{equation*}

so that large $r$ has exponentially decaying contribution.
The second is

\begin{equation*}
2h(\frac{s}{\Delta}) + 2\frac{\Delta-s}{\Delta} h(\frac{t}{\Delta-s}) < \frac{3}{4} \frac{(\Delta - s - k_a)(\Delta - s - k_b)}{\Delta (\Delta - s)}
\end{equation*}

so that the contribution for $r=1$ approaches $0$.
It is clear that we can reduce the two inequalities to just the second one.
By rewriting the variables $s = \sigma \Delta, t = \tau \Delta, k_a = \rho_a \Delta, k_b = \rho_b \Delta$, we see that if

\begin{equation*}
2 h(\sigma) + 2 (1-\sigma)h(\frac{\tau}{1-\sigma}) < \frac{3}{4} \frac{(1 - \sigma - \rho_a)(1 - \sigma - \rho_b)}{1-\sigma}
\end{equation*}

as $\Delta$ goes to infinity, with probability tending to $1$,
for any $s = \sigma \Delta$-punctured code $C_A', C_B'$,
all non-zero codewords in $\tensorcode(C_A', C_B')$ have at least one row or one column with weight $\ge t = \tau \Delta$.

The idea is that knowing the entries of specific $r$ rows and $r$ columns of a rank $r$ matrix $c$ is sufficient to determine rest of the entries.
So to upper bound the number of rank $r$ matrices $c$ where each row and each column has weight $< t$,
we suffice to upper bound the number of  configuration supported on $r$ rows and $r$ columns,
$c_{I_a \times [n_b] \cup [n_a] \times I_b}$,
such that each row and each column has weight $< t$,
where $I_a$ and $I_b$ are the indices of some $r$ rows and $r$ columns.
We discuss the details below.

Given $c \in \cM(\Delta, r, t)$ of rank $r$ matrices where each row and each column has weight $< t$.
Because $c$ has rank $r$, there exists a $r \times r$ submatrix with rank $r$, i.e. full rank.
Let the $r$ rows and $r$ columns be indexed by $I_a \subset [n_a]$ and $I_b \subset [n_b]$.
We now claim that knowing the value of $c$ on $I_a \times [n_b] \cup [n_a] \times I_b$ is enough to uniquely determine $c$ through the following identity:

\begin{equation*}
c_{i_a, i_b} = c_{I_a \times i_b} c_{I_a \times I_b}^{-1} c_{i_a \times I_b} = \sum_{j_a \in I_a, j_b \in I_b} c_{j_a, i_b} (c_{I_a \times I_b})^{-1}_{j_a, j_b} c_{i_a, j_b}
\end{equation*}

where $c_{I_a \times I_b}^{-1}$ is the matrix inverse.

Because $c$ has rank $r$, the vector $c_{(i_a \cup I_a) \times i_b}$ is a linear combination of the vectors $c_{(i_a \cup I_a) \times j_b}$ for $j_b \in I_b$.
Say $c_{(i_a \cup I_a) \times i_b} = \sum_{j_b \in I_b} v_{j_b} c_{(i_a \cup I_a) \times j_b}$.
Using the coordinates of $I_a$, $c_{I_a \times i_b} = \sum_{j_b \in I_b} v_{j_b} c_{I_a \times j_b}$, one has $v = c_{I_a \times i_b} c_{I_a \times I_b}^{-1}$.
When plug $v$ back to $c_{i_a, i_b} = \sum_{j_b \in I_b} v_{j_b} c_{i_a, j_b}$, we obtain the desired identity.

(0,0) rectangle (6,-6);
(2,0) rectangle (2,-6);
(3,0) rectangle (3,-6);
(0,-2) rectangle (6,-2);
(0,-3) rectangle (6,-3);

(2.5,-2.5)node$c_{i_a, i_b}$;
(1,-2.5)node$c_{i_a \times I_b}$;
(2.5,-1)node$c_{I_a \times i_b}$;
(1,-1)node$c_{I_a \times I_b}$;

(0,0) –node[auto]$I_b$ (2,0) –node[auto]$i_b$ (3,0);
(0,0) –node[auto,swap]$I_a$ (0,-2) –node[auto,swap]$i_a$ (0,-3);

(3,.7)node[rotate=90]$\left.\rule{0pt}{\dimexpr 3cm}\right\}$;
(-.7,-3)node$\left\{\rule{0pt}{\dimexpr 3cm}\right.$;

(0,0) – (6,0) – (6,-2) – (2,-2) – (2,-6) – (0,-6) – cycle;
in-6, -5.5, ..., 6
(, 0) – +(6, -6);

Given $c$ of rank $r$. If $c_{I_a \times I_b}$ has rank $r$, then $c$ is uniquely determined by the value on the shaded region $c_{I_a \times [n_b] \cup [n_a] \times I_b}$.

So to bound the size of $\cM(\Delta, r, t)$, we suffice to bound the number of choices of $I_a$ and $I_b$ and the number of choices of $c_{I_a \times [n_b] \cup [n_a] \times I_b}$.

We first bound the number of choices of $I_a$ and $I_b$. $I_a$ and $I_b$ each has size $r$ and is a subset of $[\Delta]$, so the number of choices

\begin{equation*}
\le \binom{\Delta}{r}^2 \le 2^{2\Delta h(\frac{r}{\Delta})}
\end{equation*}

where we use the inequality $\binom{\Delta}{r} \le 2^{\Delta h(\frac{r}{\Delta})}$.

Next, we bound the number of choices of $c_{I_a \times [n_b] \cup [n_a] \times I_b}$ given $I_a$ and $I_b$.
It is enough to find $r$ columns and $r$ rows of weight $< t$,
so the number of choices

\begin{equation*}
\le \left(\sum_{i=0}^{t-1} \binom{\Delta}{i}\right)^{2r} \le 2^{2r\Delta h(\frac{t}{\Delta})}
\end{equation*}

where we use the inequality $\sum_{i=0}^{t-1} \binom{\Delta}{i} \le 2^{\Delta h(\frac{t}{\Delta})}$ when $t \le \Delta/2$.

Multiply the two, we have

\begin{equation*}
\abs{\cM(\Delta, r, t)} \le 2^{2r\Delta h(\frac{t}{\Delta}) + 2\Delta h(\frac{r}{\Delta})}.
\end{equation*}

Let $H_A\colon \F2^\Delta \to \F2^{m_a}$ and $H_B\colon \F2^\Delta \to \F2^{m_b}$ be the parity check matrices of $C_A$ and $C_B$ where $m_a = \Delta - k_a$ and $m_b = \Delta - k_b$.
Using the chain rule we have

\begin{equation*}
\Pr[c \in \tensorcode(C_A, C_B)] = \Pr[(H_A \otimes H_B) c = 0] = \sum_{r'=0}^{\min(r,m_a)} \Pr[\rk(H_A c) = r'] \cdot \Pr[H_B (H_A c) = 0 | \rk(H_A c) = r']
\end{equation*}

where $\rk$ is the rank and $H_B (H_A c)$ means we first apply $H_A$ to $c$ then apply $H_B$ to $H_A c$.
So what remains is to bound $\Pr[\rk(H_A c) = r']$ and $\Pr[H_B (H_A c) = 0 | \rk(H_A c) = r']$.

To do so, we introduce the Gaussian binomial coefficients.
The Gaussian binomial coefficients is $\binom{x}{y}_2 \coloneqq \frac{[x]_2!}{[y]_2![x-y]_2!}$ where $[x]_2 \coloneqq 2^x - 1$ and $[x]_2! \coloneqq \prod_{i=1}^x [i]_2$.
It is known that $\binom{x}{y}_2$ is the number of $y$-dimensional vector subspace of an $x$-dimensional vector space over $\F2$.
This is the primary use of the Gaussian binomial coefficients.

Using the knowledge of Gaussian binomial coefficient, we first bound $\Pr[\rk(H_A c) = r']$.
The idea is to consider all the rank $r'$ subspaces of $\F2^{m_a}$ and compute the probability that $H_A c$ is contained in one of those subspaces.
The number of the rank $r'$ subspace is simply $\binom{m_a}{r'}_2$.
Let $V$ be the column space of $c$, i.e. the vector space spanned by the column vectors of $c$.
Because $c$ has rank $r$, $\dim V = r$.
Now, we consider a rank $r'$ subspace $V'$ of $\F2^{m_a}$ and ask what is the probability that $H_A V \subset V'$, i.e.
$V \subset H_A^{-1} V' = \set{v \in \F2^\Delta: H_A v \in V'}$.
Because $H_A$ is full rank, $\dim H_A^{-1} V' = r' + k_a$.
So the probability is $\binom{r'+k_a}{r}_2 \big/ \binom{\Delta}{r}_2$, i.e. the number of rank $r$ subspace in $H_A^{-1} V'$ divided by the number of rank $r$ subspace in $\F2^\Delta$.
Overall we have,

\begin{equation*}
\Pr[\rk(H_A c) = r'] \le \binom{m_a}{r'}_2 \binom{r'+k_a}{r}_2 \bigg/ \binom{\Delta}{r}_2.
\end{equation*}

(Note that this is an inequality because $H_A V$ could have rank less than $r'$.)

Now we bound $\Pr[H_B (H_A c) = 0 | \rk(H_A c) = r']$.
Let $W$ be the row space of $H_A c$ and $W'$ be the kernel $H_B^{-1}(0)$.
We know $\dim W = r'$, because $\rk(H_A c) = r'$, and $H_B^{-1}(0) = k_b$, because $H_B$ has full rank.

\begin{equation*}
\Pr[H_B (H_A c) = 0 | \rk(H_A c) = r'] = \binom{k_b}{r'}_2 \bigg/ \binom{\Delta}{r'}_2.
\end{equation*}

Combine the above results, we have

\begin{equation*}
\Pr[c \in \tensorcode(C_A, C_B)]
\le \sum_{r'=0}^{\min(r, m_a)} \left( \binom{m_a}{r'}_2 \binom{r'+k_a}{r}_2 \middle/ \binom{\Delta}{r}_2 \right) \cdot \left( \binom{k_b}{r'}_2 \middle/ \binom{\Delta}{r'}_2 \right).
\end{equation*}

What is left to be done is to simplify the formula.

We first use the following bound of the Gaussian binomial coefficient.

\begin{equation*}
2^{y(x-y)} \le \binom{x}{y}_2 \le 8 \cdot 2^{y(x-y)}
\end{equation*}

When $x=y$ or $y=0$, the inequality holds trivially.
Otherwise, we will show when $z > 0$

\begin{equation*}
\frac{1}{4} \cdot 2^{z(z+1)/2} \le [z]_2! \le \frac{1}{2} \cdot 2^{z(z+1)/2}
\end{equation*}

which implies the desired result.
By definition
$[z]_2! = \prod_{i=1}^z (2^i-1) = 2^{z(z+1)/2} \prod_{i=1}^z (1-2^{-i})$.
So we suffice to bound $\prod_{i=1}^z (1-2^{-i})$.
For the upper bound, it is clear that $\prod_{i=1}^z (1-2^{-i}) \le 1/2$.
For the lower bound, we have

\begin{equation*}
1/4 = (1-2^{-1})(1-2^{-1}) \le (1-2^{-1})(1-2^{-2})(1-2^{-2})
\le ... \le (1-2^{-z}) \prod_{i=1}^z (1-2^{-i}) \le \prod_{i=1}^z (1-2^{-i})
\end{equation*}

where we apply $(1-2^{-i}) \le (1-2^{-i-1})^2$ iteratively.

The bound on the Gaussian binomial coefficient implies
[Note that $512$ can be improved to $8$ by bounding $\binom{r'+k_a}{r}_2 / \binom{\Delta}{r}_2 \le 2^{-r(\Delta - r' - k_a)}$
and $\binom{k_b}{r'}_2 / \binom{\Delta}{r'}_2 \le 2^{-r'(\Delta - k_b)}$.]

\begin{equation*}
\Pr[c \in \tensorcode(C_A, C_B)]
\le 512 \sum_{r'=0}^{\min(r, m_a)} 2^{-(m_a-r')(r-r')-m_b r'}.
\end{equation*}

Finally, we show $-(m_a-r')(r-r')-m_b r' \le -\frac{3}{4} \frac{m_a m_b}{\Delta} r$ for $0 \le r' \le m_a$.
We first rewrite $-(m_a-r')(r-r')-m_b r' = (m_a-r')(m_b-(r-r')) - m_a m_b$.
Then divide both sides by $m_a m_b$.
If $m_b - (r-r') \ge 0$ we have

\begin{align*}
(1-\frac{r'}{m_a})(1-\frac{r-r'}{m_b}) - 1
&\le (1-\frac{r'}{\Delta})(1-\frac{r-r'}{\Delta}) - 1 \\
&= -\frac{r}{\Delta} + \frac{r'(r-r')}{\Delta^2} \\
&\le -\frac{r}{\Delta} + \frac{\frac{1}{4}r^2}{\Delta^2} \\
&\le -\frac{3}{4} \frac{r}{\Delta}.
\end{align*}

Note the first inequality uses $m_a - r' \ge 0$ and $m_b - (r-r') \ge 0$.
If $m_b - (r-r') < 0$ we have

\begin{equation*}
(1-\frac{r'}{m_a})(1-\frac{r-r'}{m_b}) - 1
\le - 1
\le -\frac{3}{4} \frac{r}{\Delta}.
\end{equation*}

\begin{equation*}
\Pr[c \in \tensorcode(C_A, C_B)] \le 512 (r+1) 2^{-\frac{3}{4} \frac{m_a m_b}{\Delta} r}.
\end{equation*}

We now study the second lemma. We first recall the statement.

The intuition is that one can construct $c_a$ using $c_{([n_a] - I_a) \times I_b}$ by requiring each column to be a codeword in $C_A$.
Because each column is only missing $|I_a| < d_1$ values,
which is less than the distance,
one can uniquely recover the codeword.
This intuition is correct, but requires more steps to make it rigorous as we will do below.

Let $H_A' = H_A|_{I_a}\colon \F2^{I_a} \to \F2^{m_a}$ be the restriction of $H_A$.
Because $C_A$ has distance $d_1$ and $|I_a| < d_1$, $H_A'$ is injective.
Since $H_A'$ is an injective linear map, there exists (not unique) a left inverse linear map $J_A'\colon \F2^{m_a} \to \F2^{I_a}$, such that

\begin{equation*}
J_A' H_A' = I_{I_a}
\end{equation*}

where $I_{I_a}$ is the identity for $\F2^{I_a}$.
Similarly, we can define $H_B'$ and $J_B'$.

Using $J_A'$ we can explicitly write down a recovery map for the erased codeword.

For any codeword $v \in C_A \subset \F2^{n_a}$,
one can recover $v$ from $v_{[n_a]-I_a}$ through the following linear map

\begin{equation*}
v = (J_A' H_A|_{[n_a]-I_a} \Vert I_{[n_a] - I_a}) v_{[n_a]-I_a}
\end{equation*}

where $\Vert$ is concatenation.

Notice that $J_A' H_A|_{[n_a]-I_a} \colon \F2^{[n_a] - I_a} \to \F2^{I_a}$ and $I_{[n_a] - I_a} \colon \F2^{[n_a] - I_a} \to \F2^{[n_a] - I_a}$, so
$J_A' H_A|_{[n_a]-I_a} \Vert I_{[n_a] - I_a} \colon \F2^{[n_a] - I_a} \to \F2^{[n_a]}$.

Because $v \in C_A$ is a codeword, we have $H_A v = 0$,
so $H_A' v_{I_a} = H_A|_{[n_a]-I_a} v_{[n_a]-I_a}$.
Multiply both side by $J_A'$, we obtain $J_A' H_A' v_{I_a} = v_{I_a} = J_A' H_A|_{[n_a]-I_a} v_{[n_a]-I_a}$ which is the desired result.

This suggests one can recover $c_a$ from $c_{([n_a]-I_a) \times I_b}$
by mapping each column supported on $[n_a]-I_a$ to the full codeword supported on $[n_a]$.
However, we need to first show each $c_{([n_a]-I_a) \times i_b}$ is a truncation of a codeword in $C_A$ which we will do now.
Because $H_A H_B c = 0$, each column of $H_B c$ is a codeword in $C_A$.
Because $J_B' H_B' = I_{I_b}$, $c_{([n_a]-I_a) \times I_b} = J_B' H_B' c_{([n_a]-I_a) \times I_b} = (J_B' (H_B c))|_{([n_a]-I_a) \times [m_b]}$,
indeed each column of $c_{([n_a]-I_a) \times I_b}$ is a truncation of codeword.

This leads to the following choice of $c_a$ and $c_b$

\begin{equation*}
c_a = ((J_A' H_A|_{[n_a] - I_a} \Vert I_{[n_a] - I_a}) c_{([n_a]-I_a) \times I_b}) \big\Vert 0_{[n_a] \times ([n_b] - I_b)}
\end{equation*}

\begin{equation*}
c_b = ((J_B' H_B|_{[n_b] - I_b} \Vert I_{[n_b] - I_b}) c_{I_a \times ([n_b]-I_b)}) \big\Vert 0_{([n_a] - I_a) \times [n_b]}.
% \footnote{It is also possible to combine the discussion with the formula and directly writes
%   \begin{equation*}
%     c_a = (J_B' \Vert 0_{[n_b] - I_b}) (J_A' H_A|_{[n_a] - I_a} \Vert I_{[n_a] - I_a}) H_B|_{I_b} c_{([n_a]-I_a) \times I_b}.
%   \end{equation*}
% }
\end{equation*}

By construction $c_a$ is supported on $[n_a] \times I_b$ and $c_b$ is supported on $I_a \times [n_b]$.
And from the discussion we know $c_a \in C_A \otimes \F2^{n_b}$ and $c_b \in \F2^{n_a} \otimes C_B$.
What is left to be shown is $c = c_a + c_b$.

By construction we know $c$ agrees with $c_a + c_b$ on everywhere outside of $I_a \otimes I_b$,
that is $c' \coloneqq c + c_a + c_b$ is supported on $I_a \otimes I_b$.
Now, $H_B c'$ is supported on $I_a \otimes [m_b]$ and each column is a codeword in $C_A$ (because $H_A (H_B c') = 0$).
Because the distance of $C_A$ is $d_1$ and $|I_a|<d_1$, we have $H_B c' = 0$.
Similarly, $c'$ is supported on $I_a \otimes I_b$ and each row is a codeword in $C_B$ (because $H_B c' = 0$).
Because the distance of $C_B$ is $d_1$ and $|I_b|<d_1$, we have $c' = 0$.
This concludes the first half of the lemma.

The second half is simple.
Since each column of $c_a$ is a codeword, it is either $0$ or has weight $\ge d_1$.
For those with weight $\ge d_1$, after removing the indices in $I_a$ it still have weight $> d_1/2$.
That is

\begin{equation*}
\abs{c_{([n_a] - I_a) \times I_b}} \ge \frac{d_1}{2} \norm{c_a}_{[n_b]}.
\end{equation*}

\begin{equation*}
\abs{c_{I_a \times ([n_b] - I_b)}} \ge \frac{d_1}{2} \norm{c_b}_{[n_a]}.
\end{equation*}

\begin{equation*}
\abs{c} \ge \frac{d_1}{2} (\norm{c_a}_{[n_b]} + \norm{c_b}_{[n_a]}).
\end{equation*}

[1]
Irit Dinur, Shai Evra, Ron Livne, Alexander Lubotzky, and Shahar Mozes.
Locally testable codes with constant rate, distance, and locality.
arXiv preprint arXiv:2111.04808, 2021.

[2]
Daniel Gottesman.
Fault-tolerant quantum computation with constant overhead.
arXiv preprint arXiv:1310.2984, 2013.

[3]
A Yu Kitaev.
Fault-tolerant quantum computation by anyons.
Annals of Physics, 303(1):2–30, 2003.

[4]
Jean-Pierre Tillich and Gilles Zémor.
Quantum ldpc codes with positive rate and minimum distance
proportional to the square root of the blocklength.
IEEE Transactions on Information Theory, 60(2):1193–1202,

[5]
Michael H Freedman, David A Meyer, and Feng Luo.
Z2-systolic freedom and quantum codes.
In Mathematics of quantum computation, pages 303–338. Chapman
and Hall/CRC, 2002.

[6]
Shai Evra, Tali Kaufman, and Gilles Zémor.
Decodable quantum ldpc codes beyond the square root distance barrier
using high dimensional expanders.
arXiv preprint arXiv:2004.07935, 2020.

[7]
Tali Kaufman and Ran J Tessler.
New cosystolic expanders from tensors imply explicit quantum ldpc
codes with $\Omega(\sqrt{n}\log^k n)$ distance.
arXiv preprint arXiv:2008.09495, 2020.

[8]
Matthew B Hastings, Jeongwan Haah, and Ryan O'Donnell.
Fiber bundle codes: breaking the n 1/2 polylog (n) barrier for
quantum ldpc codes.
In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory
of Computing, pages 1276–1288, 2021.

[9]
Pavel Panteleev and Gleb Kalachev.
Quantum ldpc codes with almost linear minimum distance.
IEEE Transactions on Information Theory, 2021.

[10]
Nikolas P Breuckmann and Jens N Eberhardt.
Balanced product quantum codes.
IEEE Transactions on Information Theory, 67(10):6653–6674,

[11]
Pavel Panteleev and Gleb Kalachev.
Asymptotically good quantum and locally testable classical ldpc
arXiv preprint arXiv:2111.03654, 2021.

[12]
Anthony Leverrier and Gilles Zémor.
Quantum tanner codes.
arXiv preprint arXiv:2202.13641, 2022.

[13]
Anthony Leverrier, Vivien Londe, and Gilles Zémor.
Towards local testability for quantum coding.
Quantum, 6:661, 2022.

[14]
Ting-Chun Lin and Min-Hsiu Hsieh.
$c^3$-local testable codes from lossless expanders.
arXiv preprint arXiv:2201.11369, 2022.

[15]
Oded Goldreich.
Short locally testable codes and proofs: A survey in two parts.
In Property testing, pages 65–104. Springer, 2010.

[16]
Nikolas P Breuckmann and Jens Niklas Eberhardt.
Quantum low-density parity-check codes.
PRX Quantum, 2(4):040101, 2021.

[17]
Eric Dennis, Alexei Kitaev, Andrew Landahl, and John Preskill.
Topological quantum memory.
Journal of Mathematical Physics, 43(9):4452–4505, 2002.

[18]
Nicolas Delfosse and Naomi H Nickerson.
Almost-linear time decoding algorithm for topological codes.
Quantum, 5:595, 2021.

[19]
Guillaume Duclos-Cianci and David Poulin.
Fast decoders for topological quantum codes.
Phys. Rev. Lett., 104:050504, Feb 2010.

[20]
Pavel Panteleev and Gleb Kalachev.
Degenerate quantum ldpc codes with good finite length performance,

[21]
Michael Sipser and Daniel A Spielman.
Expander codes.
IEEE transactions on Information Theory, 42(6):1710–1722,

[22]
Anthony Leverrier, Jean-Pierre Tillich, and Gilles Zémor.
Quantum expander codes.
In 2015 IEEE 56th Annual Symposium on Foundations of Computer
Science, pages 810–824. IEEE, 2015.

[23]
Ting-Chun Lin and Min-Hsiu Hsieh.
Good quantum ldpc codes with linear time decoder from lossless
arXiv preprint arXiv:2203.03581, 2022.

[24]
Nathan Linial* and Roy Meshulam*.
Homological connectivity of random 2-complexes.
Combinatorica, 26(4):475–487, 2006.

[25]
Mikhail Gromov.
Singularities, expanders and topology of maps. part 2: From
combinatorics to topology via algebraic isoperimetry.
Geometric and Functional Analysis, 20(2):416–526, 2010.

[26]
Fernando Granha Jeronimo, Shashank Srivastava, and Madhur Tulsiani.
Near-linear time decoding of ta-shma’s codes via splittable
In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory
of Computing, pages 1527–1536, 2021.

[27]
Tali Kaufman and Izhar Oppenheim.
High order random walks: Beyond spectral gap.
Combinatorica, pages 1–37, 2020.

[28]
Nima Anari, Kuikui Liu, Shayan Oveis Gharan, and Cynthia Vinzant.
Log-concave polynomials ii: high-dimensional walks and an fpras for
counting bases of a matroid.
In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory
of Computing, pages 1–12, 2019.

[29]
Vedat Levi Alev and Lap Chi Lau.
Improved analysis of higher order random walks and applications.
arXiv preprint arXiv:2001.02827, 2020.

[30]
Nima Anari, Kuikui Liu, and Shayan Oveis Gharan.
Spectral independence in high-dimensional expanders and applications
to the hardcore model.
arXiv preprint arXiv:2001.00303, 2020.

[31]
Yotam Dikstein, Irit Dinur, Yuval Filmus, and Prahladh Harsha.
Boolean function analysis on high-dimensional expanders.
In Approximation, Randomization, and Combinatorial Optimization.
Algorithms and Techniques (APPROX/RANDOM 2018). Schloss
Dagstuhl-Leibniz-Zentrum fuer Informatik, 2018.

[32]
Mitali Bafna, Max Hopkins, Tali Kaufman, and Shachar Lovett.
Hypercontractivity on high dimensional expanders.
arXiv preprint arXiv:2111.09444, 2021.

[33]
Tom Gur, Noam Lifshitz, and Siqi Liu.
Hypercontractivity on high dimensional expanders.
arXiv preprint arXiv:2111.09375, 2021.

[34]
Irit Dinur and Tali Kaufman.
High dimensional expanders imply agreement expanders.
In 2017 IEEE 58th Annual Symposium on Foundations of Computer
Science (FOCS), pages 974–985. IEEE, 2017.

[35]
Yotam Dikstein and Irit Dinur.
Agreement testing theorems on layered set systems.
In 2019 IEEE 60th Annual Symposium on Foundations of Computer
Science (FOCS), pages 1495–1524. IEEE, 2019.

[36]
Irit Dinur, Yuval Filmus, Prahladh Harsha, and Madhur Tulsiani.
Explicit sos lower bounds from high-dimensional expanders.
arXiv preprint arXiv:2009.05218, 2020.

[37]
Max Hopkins and Ting-Chun Lin.
Explicit lower bounds against $\omega(n)$-rounds of sum-of-squares.
arXiv preprint arXiv:2204.11469, 2022.

[38]
Dorit Aharonov and Lior Eldar.
Quantum locally testable codes.
SIAM Journal on Computing, 44(5):1230–1262, 2015.

[39]
Dorit Aharonov, Itai Arad, and Thomas Vidick.
Guest column: the quantum pcp conjecture.
Acm sigact news, 44(2):47–79, 2013.

[40]
Shlomo Hoory, Nathan Linial, and Avi Wigderson.
Expander graphs and their applications.
Bulletin of the American Mathematical Society, 43(4):439–561,

[41]
A Philips, R Lubotsky, and P Sarnak.
Ramanujan graphs.
Combinatorica, 8(3):261–277, 1988.

[42]
Noga Alon and Fan RK Chung.
Explicit construction of linear sized tolerant networks.
Discrete Mathematics, 72(1-3):15–19, 1988.

[43]
Eli Ben-Sasson and Madhu Sudan.
Robust locally testable codes and products of codes.
In Approximation, Randomization, and Combinatorial Optimization.
Algorithms and Techniques, pages 286–297. Springer, 2004.

[44]
Irit Dinur, Madhu Sudan, and Avi Wigderson.
Robust local testability of tensor products of ldpc codes.
In Approximation, Randomization, and Combinatorial Optimization.
Algorithms and Techniques, pages 304–315. Springer, 2006.

[45]
Allen Hatcher.
Algebraic Topology.
Cambridge University Press, 2002.

[46]
R Tanner.
A recursive approach to low complexity codes.
IEEE Transactions on information theory, 27(5):533–547, 1981.

[47]
Tali Kaufman, David Kazhdan, and Alexander Lubotzky.
Ramanujan complexes and bounded degree topological expanders.
In 2014 IEEE 55th Annual Symposium on Foundations of Computer
Science, pages 484–493. IEEE, 2014.

[48]
Shai Evra and Tali Kaufman.
Bounded degree cosystolic expanders of every dimension.
In Proceedings of the forty-eighth annual ACM symposium on
Theory of Computing, pages 36–48, 2016.

[49]
Alexander Lubotzky.
Ramanujan complexes and high dimensional expanders.
Japanese Journal of Mathematics, 9(2):137–169, 2014.

[50]
Daniel A Spielman.
Linear-time encodable and decodable error-correcting codes.
IEEE Transactions on Information Theory, 42(6):1723–1731,