# Anisotropic Topological Anderson Transitions in Chiral Symmetry Classes

Zhenyu Xiao  International Center for Quantum Materials, Peking University,
Beijing 100871, China  Kohei Kawabata  Department of Physics, Princeton
University, Princeton, New Jersey 08544, USA  Institute for Solid State
Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan  Xunlong Luo
Science and Technology on Surface Physics and Chemistry Laboratory, Mianyang
621907, China  Tomi Ohtsuki  Physics Division, Sophia University, Chiyoda-ku,
Tokyo 102-8554, Japan  Ryuichi Shindou<EMAIL_ADDRESS>International
Center for Quantum Materials, Peking University, Beijing 100871, China

###### Abstract

We study quantum phase transitions of three-dimensional disordered systems in
the chiral classes (AIII and BDI) with and without weak topological indices.
We show that the systems with a nontrivial weak topological index universally
exhibit an emergent thermodynamic phase where wave functions are delocalized
along one spatial direction but exponentially localized in the other two
spatial directions, which we call the quasi-localized phase. Our extensive
numerical study clarifies that the critical exponent of the Anderson
transition between the metallic and quasi-localized phases, as well as that
between the quasi-localized and localized phases, are different from that with
no weak topological index, signaling the new universality classes induced by
topology. The quasi-localized phase and concomitant topological Anderson
transition manifest themselves in the anisotropic transport phenomena of
disordered weak topological insulators and nodal-line semimetals, which
exhibit the metallic behavior in one direction but the insulating behavior in
the other directions.

Introduction—The last decades have seen remarkable discoveries of topological
materials [1, 2, 3]. The interplay of disorder and topology leads to new types
of quantum phase transitions, including the quantum Hall plateau transitions
[4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. The universality classes of the disorder-
driven metal-insulator transitions, known as the Anderson transitions, are
characterized by the critical exponents and scaling functions, which are
commonly believed to be determined solely by symmetry and spatial dimensions
[14]. Many theories investigated whether topology can change the universality
classes of the Anderson transitions [15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Still, the role of topology in the
Anderson transitions has been elusive.

Prime examples of three-dimensional (3D) topological materials include nodal-
line semimetals characterized by the weak topological invariant [35, 36, 37,
38]. Several recent experiments realized nodal-line semimetals in solid states
[39, 40, 41], as well as synthetic materials of ultracold atoms [42] and
photonic [43, 44] and phononic [45] systems. Despite the significant interest
in the physics of nodal-line semimetals [46, 47, 48, 49, 50], their unique
transport signatures have remained largely unexplored.

(a) topological models

(b) nontopological models

Figure 1: Phase diagrams of 3D disordered Hamiltonians in the chiral symmetry
classes (a) with and (b) without the weak topological index $\nu_{z}$. The
critical exponents $\nu$ and localization lengths $\xi_{z},\xi_{\perp}\\!\
(\perp=x,y)$ along different directions are shown for different phases. The
nontrivial critical exponents $\nu=0.82\pm 0.02$ and $\nu=1.09\pm 0.03$ are
obtained for class BDI.

In this Letter, we elucidate that the weak topological indices induce a novel
thermodynamic phase in 3D disordered systems, including topological nodal-line
semimetals, in the chiral classes. There, 3D wave functions are delocalized
along one spatial direction and exponentially localized along the other two
spatial directions—quasi-localized phase [Fig. 1(a)]. From extensive numerical
calculations, we evaluate correlation-length critical exponents of the
Anderson transitions among the metallic, quasi-localized, and localized phases
[Fig. 1(a)] and find that they are distinct from the critical exponent in
topologically trivial systems [Fig. 1(b)], signaling new universality classes
induced by the topological indices. Notably, our quasi-localized phase and
concomitant topological Anderson transition are of direct experimental
relevance in the anisotropic transport that exhibits the metallic behavior in
one direction but the insulating behavior in the other directions. While such
anisotropic transport has played an important role in condensed matter physics
[51, 52, 53, 54, 55, 56, 57, 58], our results provide its new universal
mechanism induced by the interplay of disorder and topology.

Lyapunov exponents and topological indices—We study disorder-induced quantum
phase transitions of 3D chiral-symmetric Hamiltonians $\mathcal{H}$. The
localization properties along the $\mu$ direction ($\mu=x,y,z$) are
efficiently captured by the Lyapunov exponents (LEs) along the $\mu$ direction
in the limit $L\rightarrow\infty$, which are eigenvalues of [59, 60]

$\lim_{L_{\mu}\rightarrow\infty}\log{(M^{\dagger}M)^{\frac{1}{2L_{\mu}}}}.$
(1)

Here, $M\equiv M_{L_{\mu}}M_{L_{\mu}-1}\cdots M_{1}$ is the product of
transfer matrices along the $\mu$ direction. The smallest positive LE gives
the inverse of the localization length along the $\mu$ direction [61]. In the
limit $L\to\infty$, the LEs of ${\cal H}$ form several continuous spectra
[62]. If the spectra do not include zero, the wave function is localized along
the $\mu$ direction. By contrast, if the spectra include zero, the
localization length diverges, which means the delocalization of the wave
function. The finite (infinite) localization length leads to the vanishing
(nonvanishing) conductance in the same direction, as shown in the Supplemental
Material [34].

Symmetries of Hamiltonians give constraints on the spectrum of the LEs. For
example, because of Hermiticity of ${\cal H}$, the LEs come in opposite-sign
pairs. Moreover, in the presence of chiral symmetry, ${\cal H}$ can be brought
into the block off-diagonal structure,

$\displaystyle{\cal H}=\begin{pmatrix}0&h\\\ h^{{\dagger}}&0\\\
\end{pmatrix},$ (2)

where the off-diagonal part $h$ is assumed to be a square matrix. Because of
chiral symmetry, the LEs of ${\cal H}$ reduce to the LEs of $h$ and
$h^{{\dagger}}$, which come in opposite-sign pairs, as shown in the
Supplemental Material [34]. Consequently, we only need to calculate the
product of the transfer matrices of $h$.

We demonstrate that a weak topological index $\nu_{\mu}$ imposes another
constraint on the spectrum of the LEs and plays a vital role in the emergence
of the quasi-localized phase in disordered chiral-symmetric systems. To
introduce $\nu_{\mu}$ along the $\mu$ direction in the presence of disorder,
let us insert a magnetic flux $\phi_{\mu}$ through a closed loop along the
$\mu$ direction. Then, the weak topological index $\nu_{\mu}$ is given by the
winding of $\det h\left(\phi_{\mu}\right)$ in Eq. (2) under an adiabatic
insertion of a unit flux [63, 64, 65]:

$\displaystyle\nu_{\mu}\equiv\frac{\rm
i}{L^{2}}\int^{2\pi}_{0}\frac{d\phi_{\mu}}{2\pi}\partial_{\phi_{\mu}}{\rm
Tr}\big{[}\log\big{[}h(\phi_{\mu})]\big{]}\,,$ (3)

where $L^{2}$ is the system size within the two directions perpendicular to
the $\mu$ direction. Here, $\nu_{\mu}$ is not necessarily quantized and takes
an arbitrary real number. Notably, the weak topological index $\nu_{\mu}$ and
LEs of $h$ are related to each other by [34, 66]

$\displaystyle\nu_{\mu}=\frac{1}{2L^{2}}(N_{+,\mu}-N_{-,\mu})\,,$ (4)

where $N_{+,\mu}$ and $N_{-,\mu}$ are the numbers of positive and negative LEs
of $h$ along the $\mu$ direction, respectively.

Suppose ${\cal H}$ has a mobility gap around $E=0$ and its zero-energy state
is characterized by the weak topological indices $\nu_{x}=\nu_{y}=0$,
$\nu_{z}=1$. From Eq. (4), a finite gap exists between the smallest positive
LE and the largest negative LE such that $N_{+,z}-N_{-,z}=2L^{2}$. By
contrast, when disorder is strong enough, the zero-energy state is in a
topologically-trivial localized phase with $N_{+,z}=N_{-,z}$. Between the two
localized phases, $L^{2}$ positive LEs of $h$ cross zero, and $\nu_{z}$
continuously changes from $1$ to $0$ with respect to the disorder strength,
where the localization length $\xi_{z}$ along the $z$ direction always
diverges. Within this finite range with divergent $\xi_{z}$, the zero-energy
state undergoes the Anderson transitions along the $x$ and $y$ directions, and
thus a quasi-localized phase with divergent $\xi_{z}$ and finite $\xi_{x}$ and
$\xi_{y}$ emerges. Below, we clarify its nature, obtain the critical exponents
of the Anderson transitions among the metallic, quasi-localized, and localized
phases, and demonstrate the existence of new universality classes.

Model—As a prototypical example, we study a two-orbital tight-binding model on
a 3D cubic lattice [50]

$\displaystyle{\cal H}$
$\displaystyle=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left\\{\epsilon_{\bm{r}}c^{\dagger}_{\bm{r}}\sigma_{z}c_{\bm{r}}+\left[\sum_{\mu=x,y}\left(t_{\perp}c^{\dagger}_{\bm{r}+\bm{e_{\mu}}}\sigma_{z}c_{\bm{r}}\right)\right.\right.$
$\displaystyle\left.\left.-{\rm
i}t_{\|}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{y}c_{\bm{r}}+t_{\|}^{\prime}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{z}c_{\bm{r}}+\text{H.c.}\right]\right\\}.$
(5)

Here, $c_{\bm{r}}$ is a two-component annihilation operator at the cubic
lattice site ${\bm{r}}$, $\sigma_{\mu}$ $(\mu=x,y,z)$ are Pauli matrices,
$t_{\perp}$, $t_{\|}$, $t_{\|}^{\prime}$ are real-valued parameters, and
$\epsilon_{\bm{r}}$ is a random potential that distributes uniformly in
$[-W/2,W/2]$. We assume $t_{\perp},t_{\|}>0$ for simplicity. This Hamiltonian
respects time-reversal symmetry ${\cal H}={\cal H}^{*}$ and chiral symmetry
${\cal H}=-\sigma_{x}{\cal H}\sigma_{x}$, and hence belongs to class BDI [67,
14, 3]. In addition, the ensemble of Hamiltonians is statistically invariant
under the combination of time reversal and reflection with respect to the $xy$
plane, which requires $N_{+,x}=N_{-,x}$, $N_{+,y}=N_{-,y}$ and
$\nu_{x}=\nu_{y}=0$, as shown in the Supplemental Material [34], while
$\nu_{z}$ can be nonzero. In the clean limit, the Hamiltonian has an energy
gap around $E=0$ with $\nu_{z}=1$ for $4t_{\perp}<2|t_{\|}^{\prime}|$. For
$4t_{\perp}>2|t_{\|}^{\prime}|$, by contrast, the zero-energy state forms a
nodal line in momentum space, resulting in $0<\nu_{z}<1$. In the following, we
focus on the nodal-line-semimetal phase for $t_{\|}=t_{\|}^{\prime}=1/2$,
$t_{\perp}=1$ and study the Anderson transitions of the zero modes along all
the directions. Still, we stress that the weak topological invariant
$\nu_{\mu}$, rather than a nodal line itself, is the main ingredient for the
quasi-localized phase.

Figure 2: $2L^{2}$ Lyapunov exponents (LEs) of the right-upper part $h$ of the
3D nodal-line-semimetal model ${\cal H}$ along the $z$ direction with the
quasi-1D geometry $L\times L\times L_{z}$ ($L=18$, $L_{z}=2\times 10^{6}$),
plotted as a function of the disorder strength $W$. The color scale stands for
the density $\rho(\gamma)$ of the LEs with the normalization
$\int\rho(\gamma)d\gamma=1$. The LEs of ${\cal H}$ are composed of the LEs of
$h$ and $h^{{\dagger}}$. Inset: the largest LE $\gamma_{\rm max}(W,L)$ among
the smaller $L^{2}$ LEs as a function of $W$ in the limit
$L\rightarrow\infty$, obtained by a finite-size scaling fit. The error bars
are smaller than the marks. The plot crosses zero linearly at
$W^{(z)}_{c}=29.45\pm 0.05$; $\xi_{z}\sim(W-W_{c}^{(z)})^{-\nu^{\prime}}$ with
$\nu^{\prime}=1$ for $W>W_{c}^{(z)}$.

Localization length $\xi_{z}$—Figure 2 shows the distribution of LEs $\gamma$
of $h$ in Eq. (2) for the nodal-line-semimetal model ${\cal H}$ in Eq. (5)
along the $z$ direction in the quasi-1D geometry $L\times L\times L_{z}$. The
distribution consists of two separate spectra, each of which contains $L^{2}$
LEs. The upper spectrum is always $\gamma=+\infty$ [34] and irrelevant to the
Anderson transitions. For $W\leq W_{c}^{(z)}\approx 29$, the lower spectrum
includes zero $\gamma=0$. Every positive LE in the lower spectrum for
$W<W_{c}^{(z)}$ crosses zero when we increase $W$. At each crossing point,
$N_{-,z}$ changes by one. For $L\rightarrow\infty$, the crossing points become
dense and $\nu_{z}=1-N_{-,z}/L^{2}$ changes continuously with $W$. For
$W>W_{c}^{(z)}$, all the LEs in the lower spectrum are negative (i.e.,
$N_{-,z}=L^{2}$), and the system is in a localized phase with no weak
topological index $\nu_{z}=0$. At $W=W_{c}^{(z)}$, the maximal LE in the lower
spectrum crosses zero. Notably, $W_{c}^{(z)}$ for $L\rightarrow\infty$ cannot
be determined by fitting $\xi_{z}/L$ with a standard scaling function [e.g.,
see Eq. (7)] because $\xi_{z}$ with finite $L$ diverges at some
$W<W_{c}^{(z)}$. Instead, we map the non-Hermitian matrix $h$ into a well-
localized Hermitian matrix by a similarity transformation [34, 68], where the
localization length obeys a scaling form in the strong disorder limit [69].
Then, we obtain the scaling form of the largest LE $\gamma_{\rm max}(W,L)$,

$\gamma_{\rm max}(W,L)=a/L+\gamma_{\rm max}(W,L=\infty)\,.$ (6)

We numerically verify this scaling and determine the critical disorder
strength $W^{(z)}_{c}=29.45\pm 0.05$ (inset of Fig. 2).

Figure 3: Normalized localization length $\Lambda_{x}\equiv\xi_{x}/L$ along
the $x$ direction as a function of the disorder strength $W$ in the nodal-
line-semimetal model in Eq. (5) with the quasi-1D geometry $L\times L\times
L_{x}$. The black points are the raw data with the error bars. The solid lines
for different $L$ and the dashed vertical line $W_{c}^{(x)}$ with the error
bars are the results of the fitting according to Eq. (7) with $n=3$. The
dashed line $W_{c}^{(z)}$ is evaluated by the fitting of the Lyapunov exponent
along the $z$ direction by Eq. (6). Inset: single-parameter scaling function
of $\Lambda_{x}$. $\Lambda_{\text{corrected}}$ is $\Lambda_{x}$ subtracted by
a contribution of the irrelevant scaling variable $c$ in Eq. (7), and $\phi$
is the relevant scaling variable.

Localization length $\xi_{x}$, $\xi_{y}$—The statistical symmetries mentioned
above require LEs of $h$ along the $x$ and $y$ directions to come in opposite-
sign pairs. Thus, the localization length $\xi_{x}$ along the $x$ direction is
always finite in the quasi-1D geometry with finite $L$. As shown in Fig. 3,
the normalized localization length $\Lambda_{x}(W,L)\equiv\xi_{x}(W,L)/L$
shows scale-invariant behavior at a certain disorder strength $W^{(x)}_{c}$
well below $W=W^{(z)}_{c}$, indicating a quantum phase transition at
$W=W^{(x)}_{c}<W^{(z)}_{c}$. To determine $W^{(x)}_{c}$ and the critical
exponent $\nu$, we use a finite-size scaling function and its polynomial
expansion [70, 61]. The scaling function for $\Lambda_{x}(W,L)$ is Taylor-
expanded with respect to the relevant scaling variable $\phi(w)$ and the least
irrelevant scaling variable $c$ up to the $n$th order and first order,
respectively,

$\displaystyle\Lambda_{x}(W,L)=\sum^{n}_{i=0}\sum^{1}_{j=0}a_{i,j}\big{(}\phi(w)L^{1/\nu}\big{)}^{i}\big{(}cL^{-y}\big{)}^{j},$
(7)

with $w\equiv(W-W^{(x)}_{c})/W^{(x)}_{c}$ and the scaling dimension $-y$
$(<0)$ of the least irrelevant scaling variable around a saddle-point fixed
point. The fitting is carried out by the $\chi^{2}$ fitting method, and the
confidence error bars for the optimal parameters are determined by the Monte
Carlo method, as detailed in the Supplemental Material [34].

The first row in Table 1 shows the fitting results, where
$W^{(x)}_{c}=27.24\pm 0.05$ is significantly smaller than
$W^{(z)}_{c}=29.45\pm 0.05$ and the critical exponent at $W^{(x)}_{c}$ is
evaluated as $\nu=0.82\pm 0.02$. The two different critical disorder strengths
illustrate the emergence of the three distinct phases as a function of the
disorder strength $W$ [Fig. 1(a)]. For $W<W^{(x)}_{c}$, the localization
lengths diverge along all directions (metallic phase). For $W>W^{(z)}_{c}$,
the localization lengths are finite along all directions (Anderson insulator
phase). For $W^{(x)}_{c}<W<W^{(z)}_{c}$, the localization lengths are finite
along the $x$ and $y$ directions but diverge along the $z$ direction (quasi-
localized phase), and $\nu_{z}$ continuously changes as $W$ changes. Our
extensive numerical calculations show that the quasi-localized phase with
divergent $\xi_{z}$ but finite $\xi_{x},\xi_{y}$ universally appears between
metallic and localized phases in different models with nonzero $\nu_{z}$, as
shown in the Supplemental Material [34]. The consistent critical exponent at
$W=W^{(x)}_{c}$ was also obtained in Ref. [50], while a different critical
exponent was obtained in Ref. [71] even in the same class. In this Letter, we
elucidate that this difference originates from the emergence of the quasi-
localized phase, which was not identified previously.

Table 1: Critical disorder strength $W^{(\mu)}_{c}$ and critical exponent $\nu$ for the 3D chiral classes, obtained by the polynomial fitting of the normalized localization length $\Lambda_{\mu}\equiv\xi_{\mu}/L$ along the $\mu$ direction ($\mu=x,y,z$) around critical points of different models with the quasi-one-dimensional geometry $L\times L\times L_{\mu}$. In the column “Topo”, “$\surd$” shows the nonzero weak topological index $\nu_{z}$ around the critical point, and “$\times$” shows zero topological indices in all the directions. The square brackets denote the 95% confidence interval.  Class | Topo |  $\mu$ | $W^{(\mu)}_{c}$ | $\nu$
---|---|---|---|---
BDI | $\surd$ | $x$ | 27.241[27.194,27.303] | 0.820[0.783,0.846]
AIII | $\surd$ | $x$ | 9.143[9.125,9.168] | 0.824[0.776,0.862]
BDI | $\times$ | $z$ | 23.220[23.167,23.293] | 1.089[1.005,1.128]
BDI | $\times$ | $x$ | 23.170[23.098,23.279] | 1.042[0.943,1.099]
AIII | $\times$ | $z$ | 8.091[8.074,8.096] | 1.024[0.973,1.070]

Quasi-localized phase—Now, we clarify the nature of the quasi-localized phase
induced by the weak topological index $\nu_{\mu}$. Let
$\Phi({\bm{r}})=\langle{\bm{r}}|\Phi\rangle$ be a normalized wave function.
The wave function interacts with an effective disorder potential $V_{\rm
eff}=\langle\Phi|V|\Phi\rangle=\sum_{\bm{r}}V({\bm{r}})|\Phi({\bm{r}})|^{2}$,
whose strength is given by $\langle V_{\rm eff}^{2}\rangle=W^{2}P_{2}$ with
the inverse participation ratio
$P_{2}\equiv\sum_{\bm{r}}|\Phi({\bm{r}})|^{4}$. Here, $\langle...\rangle$
denotes the disorder average: $\langle
V({\bm{r}})V({\bm{r}}^{\prime})\rangle=W^{2}\delta_{{\bm{r}},{\bm{r}}^{\prime}}$.
As long as $W$ is finite, the following argument is applicable to general
$V(\bm{r})$, including the box disorder in $[-W/2,W/2]$ used for the numerical
calculations. Let us introduce the integrated weight of the wave function in
the $z$th layer by $|\phi(z)|^{2}=\sum_{x,y}|\Phi({\bm{r}})|^{2}$ and also the
one-dimensional inverse participation ratio
$P_{2}^{z}\equiv\sum_{z}|\phi(z)|^{4}$. $P_{2}^{x},P_{2}^{y}$ can be defined
in the same manner. $P_{2}^{\mu}$ measures the localization property of
$\Phi({\bm{r}})$ along the $\mu$ direction, giving an upper bound of $P_{2}$:
$P_{2}\leq P_{2}^{\mu}\,(\mu=x,y,z)$ [34]. If the wave function is extended
along the $z$ direction (i.e., $P_{2}^{z}\sim L_{z}^{-1}$ [14]), $P_{2}$ and
$\braket{V_{\rm eff}^{2}}$ should vanish for $L_{z}\rightarrow\infty$, and
$\Phi({\bm{r}})$ must be extended along all the directions. If $P_{2}^{z}$ is
finite even for $L_{z}\rightarrow\infty$, by contrast, $P_{2}^{x}$ and
$P_{2}^{y}$ should also be finite for $L_{x},L_{y}\rightarrow\infty$.
Otherwise, $\Phi({\bm{r}})$ is extended within all the directions, which
contradicts finite $P_{2}^{z}$. In the intermediate phase discussed above, we
find that $\xi_{x}$ is finite but $\xi_{z}$ diverges. While finite $\xi_{x}$
means finite $P_{2}$ and $P_{2}^{z}$, divergent $\xi_{z}$ with finite
$P_{2}^{z}$ means that the wave function $\Phi({\bm{r}})$ must be quasi-
localized along the $z$ direction. Thus, the wave function in the intermediate
phase is localized within the $xy$ plane and delocalized only along the $z$
direction—quasi-localized phase. Here, $\Phi({\bm{r}})$ along the $z$
direction shares the same localization properties as wave functions of 1D
chiral-symmetric systems at a topological phase transition, where the 1D
topological index changes [14, 63, 64, 72, 73, 74]. The 3D system in the
intermediate phase is effectively decoupled into 1D wires because of finite
$\xi_{x,y}$.

The emergence of the quasi-localized phase in 3D systems is a consequence of
finite $P_{2}^{\rm 1D}$ at the topological phase transition of 1D chiral-
symmetric systems. Generally, when a $d^{\prime}$-dimensional wave function
$\Phi({\bm{R}})$ in ${\bm{R}}\equiv({\bm{r}},{\bm{s}})$ with
${\bm{r}}=(r_{1},\cdots,r_{d})$ and ${\bm{s}}=(s_{1},\cdots,s_{d^{\prime}-d})$
($d<d^{\prime}$) is made out of coupled $d$-dimensional wave functions
$\psi({\bm{r}})$ at a critical point, $\Phi({\bm{R}})$ is more extended than
$\psi({\bm{r}})$ along the ${\bm{r}}$ direction because of the interlayer
coupling [34]. Thus, the effective disorder strength for the
$d^{\prime}$-dimensional wave function $\Phi({\bm{R}})$ is bounded by the
$d$-dimensional inverse participation ratio $P^{\psi(\bm{r})}_{2}$ of
$\psi(\bm{r})$. When the wave function $\psi({\bm{r}})$ has finite
$P^{\psi({\bm{r}})}_{2}$ at the critical point, the effective disorder
strength can be finite, and $\Phi({\bm{R}})$ can be either extended or
localized within the ${\bm{s}}$ direction. On the other hand, when
$P^{\psi({\bm{r}})}_{2}$ is zero at the critical point, e.g., 2D critical wave
functions at the quantum Hall plateau transition, the effective disorder
strength is zero, and the $d^{\prime}$-dimensional wave function should be
always extended in both ${\bm{r}}$ and ${\bm{s}}$ directions. Notably, the 1D
topological phase transitions in all the three chiral classes are
characterized by finite $P_{2}$ [14]. In the following, we demonstrate the
quasi-localized phases also in the 3D chiral unitary class, which is
consistent with the above argument.

Model without time-reversal symmetry—We add a time-reversal-breaking but
chiral-symmetric disorder $\Delta{\cal H}$ to the model ${\cal H}$ in Eq. (5):

$\displaystyle{\cal H}_{1}={\cal H}+\Delta{\cal H},\quad\Delta{\cal
H}=\sum_{\bm{r}}\epsilon^{\prime}_{\bm{r}}c^{\dagger}_{\bm{r}}\sigma_{y}c_{\bm{r}},$
(8)

with the random potentials $\epsilon^{\prime}_{\bm{r}}$, where
$(\epsilon_{\bm{r}},\epsilon^{\prime}_{\bm{r}})=(V_{\bm{r}}\cos\theta_{\bm{r}},V_{\bm{r}}\sin\theta_{\bm{r}})$,
and $\theta_{\bm{r}}$ and $V_{\bm{r}}$ distribute uniformly in the range of
$[0,2\pi)$ and $[0,W]$, respectively. This model only respects chiral symmetry
and belongs to class AIII, in which the weak topological indices are defined
in the same manner. It shows a similar phase diagram as in the previous model
in class BDI with $W^{(z)}_{c}=9.8\pm 0.1$ and $W^{(x)}_{c}=9.14\pm 0.01$ (see
Fig. 1 and Table 1). The critical exponents are the same as those in the
models in class BDI, which suggests possible super-universality in 3D systems
in the chiral classes with the topological indices.

Models with trivial topological indices—To further clarify the role of the
topological indices, we also study a topologically trivial model in class BDI
with statistical symmetries. The statistical symmetry of time reversal
combined with reflection with respect to the $xz$ or $yz$ plane makes all
three topological indices vanish, as shown in the Supplemental Material [34].
In addition, LEs of $h$ along any direction come in opposite-sign pairs, and
the localization lengths along the $x$ and $y$ directions are the same. On
increasing the disorder strength, the model undergoes the Anderson transition,
where the normalized localization lengths $\Lambda_{x}$ and $\Lambda_{z}$
along the $x$ and $z$ directions both show scale-invariant behaviors. The
critical disorder strengths and critical exponents determined from
$\Lambda_{x}$ and $\Lambda_{z}$ are consistent with each other (see Table 1),
which suggests that the scale-invariant behavior of $\Lambda_{x}$ and
$\Lambda_{z}$ comes from the same quantum phase transition [Fig. 1(b)]. The
evaluated critical exponent $\nu=1.089[1.005,1.128]$ is different from $\nu$
at $W=W^{(x)}_{c}$ of the topological model, and consistent with $\nu$ of the
topologically trivial models in Ref. [71]. We also evaluate the critical
exponent in the chiral unitary class without weak topological indices as
$\nu=1.024[0.973,1.070]$ (Table 1), which is different from $\nu$ of the
topological models in the same symmetry class and consistent with Refs. [71,
33].

Summary and discussion—In this Letter, we show that in 3D systems in the
chiral classes, the weak topological indices induce a disorder-driven quasi-
localized phase where wave functions are delocalized only along one direction
and localized along the other two directions. The critical exponents of the
Anderson transitions among metal, quasi-localized, and localized phases are
all different (Fig. 1). We believe that these conclusions hold also in the
chiral symplectic class (class CII). Our quasi-localized phase leads to the
anisotropic transport phenomena of topological nodal-line semimetals [36, 37,
38, 39, 40, 41, 42, 43, 44, 45], where the conductance along the direction
with the divergent localization length takes finite values with larger
fluctuations, while it vanishes along the other directions in the
thermodynamic limit, as shown in the Supplemental Material [34]. The quasi-
localized phase may potentially find practical applications such as quantum
devices that control the direction of currents.

Our results are also relevant to non-Hermitian physics [75, 76, 77], where the
interplay between disorder and dissipation has recently acquired renewed
interest. In fact, all the disorder-driven phases and phase transitions in
this Letter are characterized by the LEs of the off-diagonal part $h$ in Eq.
(2), which can be considered as a non-Hermitian Hamiltonian [33]. Anisotropy
of ${\cal H}$ corresponds to nonreciprocity of $h$ and leads to transport
phenomena unique to open systems.

3D chiral-symmetric systems also host a strong topological index [1, 2, 3]. By
a similar numerical study, we find that the strong index does not lead to the
quasi-localized phases, not influencing the universality classes of the
Anderson transitions [78]. It also remains to be explored whether the quasi-
localized phase appears and whether the topological indices change the
universality classes of the Anderson transitions in 2D systems, as well as
nodal-line semimetals protected by spatial symmetry.

Acknowledgement— Z.X. thanks Zhida Song, Lingxian Kong, and Yeyang Zhang for
fruitful discussions. Z.X. and R.S. were supported by the National Basic
Research Programs of China (No. 2019YFA0308401) and by National Natural
Science Foundation of China (No. 11674011 and No. 12074008). K.K. was
supported by JSPS Overseas Research Fellowship, and Grant No. GBMF8685 from
the Gordon and Betty Moore Foundation toward the Princeton theory program.
X.L. was supported by National Natural Science Foundation of China of Grant
No. 12105253. T.O. was supported by JSPS KAKENHI Grant 19H00658 and 22H05114.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New Method for High-Accuracy Determination of the Fine-Structure Constant Based on Quantized Hall Resistance, Phys. Rev. Lett. 45, 494 (1980).
  * Chalker and Coddington [1988] J. Chalker and P. Coddington, Percolation, quantum tunnelling and the integer Hall effect, J. Phys. C 21, 2665 (1988).
  * Pruisken [1988] A. M. M. Pruisken, Universal Singularities in the Integral Quantum Hall Effect, Phys. Rev. Lett. 61, 1297 (1988).
  * Huckestein [1995] B. Huckestein, Scaling theory of the integer quantum Hall effect, Rev. Mod. Phys. 67, 357 (1995).
  * Bhaseen  _et al._ [2000] M. J. Bhaseen, I. I. Kogan, O. A. Soloviev, N. Taniguchi, and A. M. Tsvelik, Towards a field theory of the plateau transitions in the integer quantum Hall effect, Nucl. Phys. B 580, 688 (2000).
  * Slevin and Ohtsuki [2009] K. Slevin and T. Ohtsuki, Critical exponent for the quantum Hall transition, Phys. Rev. B 80, 041304 (2009).
  * Prodan  _et al._ [2010] E. Prodan, T. L. Hughes, and B. A. Bernevig, Entanglement Spectrum of a Disordered Topological Chern Insulator, Phys. Rev. Lett. 105, 115501 (2010).
  * Zhu  _et al._ [2019] Q. Zhu, P. Wu, R. N. Bhatt, and X. Wan, Localization-length exponent in two models of quantum Hall plateau transitions, Phys. Rev. B 99, 024205 (2019).
  * Puschmann  _et al._ [2019] M. Puschmann, P. Cain, M. Schreiber, and T. Vojta, Integer quantum Hall transition on a tight-binding lattice, Phys. Rev. B 99, 121301 (2019).
  * Dresselhaus  _et al._ [2022] E. J. Dresselhaus, B. Sbierski, and I. A. Gruzberg, Scaling Collapse of Longitudinal Conductance near the Integer Quantum Hall Transition, Phys. Rev. Lett. 129, 026801 (2022).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Anderson transitions, Rev. Mod. Phys. 80, 1355 (2008).
  * Asada  _et al._ [2002] Y. Asada, K. Slevin, and T. Ohtsuki, Anderson Transition in Two-Dimensional Systems with Spin-Orbit Coupling, Phys. Rev. Lett. 89, 256601 (2002).
  * Asada  _et al._ [2005] Y. Asada, K. Slevin, and T. Ohtsuki, Anderson Transition in the Three Dimensional Symplectic Universality Class, J. Phys. Soc. Jpn. 74, 238 (2005).
  * Onoda  _et al._ [2007] M. Onoda, Y. Avishai, and N. Nagaosa, Localization in a Quantum Spin Hall System, Phys. Rev. Lett. 98, 076802 (2007).
  * Obuse  _et al._ [2007] H. Obuse, A. Furusaki, S. Ryu, and C. Mudry, Two-dimensional spin-filtered chiral network model for the ${\mathbb{Z}}_{2}$ quantum spin-Hall effect, Phys. Rev. B 76, 075301 (2007).
  * Ryu  _et al._ [2007] S. Ryu, C. Mudry, H. Obuse, and A. Furusaki, ${\mathbb{Z}}_{2}$ Topological Term, the Global Anomaly, and the Two-Dimensional Symplectic Symmetry Class of Anderson Localization, Phys. Rev. Lett. 99, 116601 (2007).
  * Nomura  _et al._ [2007] K. Nomura, M. Koshino, and S. Ryu, Topological Delocalization of Two-Dimensional Massless Dirac Fermions, Phys. Rev. Lett. 99, 146806 (2007).
  * Mirlin  _et al._ [2010] A. D. Mirlin, F. Evers, I. V. Gornyi, and P. M. Ostrovsky, Anderson Transitions: Criticality, Symmetries, and Topologies, Int. J. Mod. Phys. B 24, 1577 (2010).
  * König  _et al._ [2012] E. J. König, P. M. Ostrovsky, I. V. Protopopov, and A. D. Mirlin, Metal-insulator transition in two-dimensional random fermion systems of chiral symmetry classes, Phys. Rev. B 85, 195130 (2012).
  * Ringel  _et al._ [2012] Z. Ringel, Y. E. Kraus, and A. Stern, Strong side of weak topological insulators, Phys. Rev. B 86, 045102 (2012).
  * Fu and Kane [2012] L. Fu and C. L. Kane, Topology, Delocalization via Average Symmetry and the Symplectic Anderson Transition, Phys. Rev. Lett. 109, 246605 (2012).
  * Slevin and Ohtsuki [2016] K. Slevin and T. Ohtsuki, Estimate of the Critical Exponent of the Anderson Transition in the Three and Four-Dimensional Unitary Universality Classes, J. Phys. Soc. Jpn. 85, 104712 (2016).
  * Roy  _et al._ [2017] B. Roy, Y. Alavirad, and J. D. Sau, Global Phase Diagram of a Three-Dimensional Dirty Topological Superconductor, Phys. Rev. Lett. 118, 227002 (2017).
  * Luo  _et al._ [2018] X. Luo, B. Xu, T. Ohtsuki, and R. Shindou, Quantum multicriticality in disordered Weyl semimetals, Phys. Rev. B 97, 045129 (2018).
  * Yoshioka  _et al._ [2018] N. Yoshioka, Y. Akagi, and H. Katsura, Learning disordered topological phases by statistical recovery of symmetry, Phys. Rev. B 97, 205110 (2018).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, R. Queiroz, R. Ilan, B. A. Bernevig, and A. Stern, Delocalization Transition of a Disordered Axion Insulator, Phys. Rev. Lett. 127, 016602 (2021).
  * Son and Raghu [2021] J. H. Son and S. Raghu, Three-dimensional network model for strong topological insulator transitions, Phys. Rev. B 104, 125142 (2021).
  * Pan  _et al._ [2021] Z. Pan, T. Wang, T. Ohtsuki, and R. Shindou, Renormalization group analysis of Dirac fermions with a random mass, Phys. Rev. B 104, 174205 (2021).
  * Wang  _et al._ [2021a] T. Wang, Z. Pan, T. Ohtsuki, I. A. Gruzberg, and R. Shindou, Multicriticality of two-dimensional class-D disordered topological superconductors, Phys. Rev. B 104, 184201 (2021a).
  * Luo  _et al._ [2022] X. Luo, Z. Xiao, K. Kawabata, T. Ohtsuki, and R. Shindou, Unifying the Anderson transitions in Hermitian and non-Hermitian systems, Phys. Rev. Research 4, L022035 (2022).
  * [34] See Supplemental Material for summary of known critical exponents in topological insulators, inverse participation ratio, transfer matrix analyses, a relation between weak topological indices and Lyapunov exponents, statistical symmetry, details of finite-size scaling analyses, quasi-localized phase in other topological models, Anderson transition in non-topological models, and detailed numerical results of conductance in metal, quasi-localized and Anderson localized phases, which includes Refs. [79, 80, 81, 82, 83, 84, 85].
  * Burkov  _et al._ [2011] A. A. Burkov, M. D. Hook, and L. Balents, Topological nodal semimetals, Phys. Rev. B 84, 235126 (2011).
  * Schnyder and Brydon [2015] A. P. Schnyder and P. M. R. Brydon, Topological surface states in nodal superconductors, J. Phys.: Condens. Matter 27, 243201 (2015).
  * Fang  _et al._ [2016] C. Fang, H. Weng, X. Dai, and Z. Fang, Topological nodal line semimetals, Chinese Phys. B 25, 117106 (2016).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Bian  _et al._ [2016] G. Bian, T.-R. Chang, R. Sankar, S.-Y. Xu, H. Zheng, T. Neupert, C.-K. Chiu, S.-M. Huang, G. Chang, I. Belopolski, D. S. Sanchez, M. Neupane, N. Alidoust, C. Liu, B. Wang, C.-C. Lee, H.-T. Jeng, C. Zhang, Z. Yuan, S. Jia, A. Bansil, F. Chou, H. Lin, and M. Z. Hasan, Topological nodal-line fermions in spin-orbit metal PbTaSe2, Nat. Commun. 7, 10556 (2016).
  * Schoop  _et al._ [2016] L. M. Schoop, M. N. Ali, C. Straßer, A. Topp, A. Varykhalov, D. Marchenko, V. Duppel, S. S. P. Parkin, B. V. Lotsch, and C. R. Ast, Dirac cone protected by non-symmorphic symmetry and three-dimensional Dirac line node in ZrSiS, Nat. Commun. 7, 11696 (2016).
  * Chen  _et al._ [2022] C. Chen, X.-T. Zeng, Z. Chen, Y. X. Zhao, X.-L. Sheng, and S. A. Yang, Second-Order Real Nodal-Line Semimetal in Three-Dimensional Graphdiyne, Phys. Rev. Lett. 128, 026405 (2022).
  * Song  _et al._ [2019] B. Song, C. He, S. Niu, L. Zhang, Z. Ren, X.-J. Liu, and G.-B. Jo, Observation of nodal-line semimetal with ultracold fermions in an optical lattice, Nat. Phys. 15, 911 (2019).
  * Gao  _et al._ [2018] W. Gao, B. Yang, Biao Tremain, H. Liu, Q. Guo, L. Xia, A. P. Hibbins, and S. Zhang, Experimental observation of photonic nodal line degeneracies in metacrystals, Nat. Commun. 9, 950 (2018).
  * Xia  _et al._ [2019] L. Xia, Q. Guo, B. Yang, J. Han, C.-X. Liu, W. Zhang, and S. Zhang, Observation of Hourglass Nodal Lines in Photonics, Phys. Rev. Lett. 122, 103903 (2019).
  * Deng  _et al._ [2019] W. Deng, J. Lu, F. Li, X. Huang, M. Yan, J. Ma, and Z. Liu, Nodal rings and drumhead surface states in phononic crystals, Nat. Commun. 10, 1769 (2019).
  * Nandkishore [2016] R. Nandkishore, Weyl and Dirac loop superconductors, Phys. Rev. B 93, 020506 (2016).
  * Sur and Nandkishore [2016] S. Sur and R. Nandkishore, Instabilities of Weyl loop semimetals, New J. Phys. 18, 115006 (2016).
  * Syzranov and Skinner [2017] S. V. Syzranov and B. Skinner, Electron transport in nodal-line semimetals, Phys. Rev. B 96, 161105 (2017).
  * Gonçalves  _et al._ [2020] M. Gonçalves, P. Ribeiro, E. V. Castro, and M. A. N. Araújo, Disorder-Driven Multifractality Transition in Weyl Nodal Loops, Phys. Rev. Lett. 124, 136405 (2020).
  * Luo  _et al._ [2020] X. Luo, B. Xu, T. Ohtsuki, and R. Shindou, Critical behavior of Anderson transitions in three-dimensional orthogonal classes with particle-hole symmetries, Phys. Rev. B 101, 020202 (2020).
  * Cohen  _et al._ [1974] M. J. Cohen, L. B. Coleman, A. F. Garito, and A. J. Heeger, Electrical conductivity of tetrathiofulvalinium tetracyanoquinodimethan $({\rm TTF})({\rm TCNQ})$, Phys. Rev. B 10, 1298 (1974).
  * Lilly  _et al._ [1999] M. P. Lilly, K. B. Cooper, J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Evidence for an Anisotropic State of Two-Dimensional Electrons in High Landau Levels, Phys. Rev. Lett. 82, 394 (1999).
  * Du  _et al._ [1999] R. R. Du, D. C. Tsui, H. L. Stormer, L. N. Pfeiffer, K. W. Baldwin, and K. W. West, Strongly anisotropic transport in higher two-dimensional Landau levels, Solid State Commun. 109, 389 (1999).
  * Ando  _et al._ [2002] Y. Ando, K. Segawa, S. Komiya, and A. N. Lavrov, Electrical Resistivity Anisotropy from Self-Organized One Dimensionality in High-Temperature Superconductors, Phys. Rev. Lett. 88, 137005 (2002).
  * Borzi  _et al._ [2007] R. A. Borzi, S. A. Grigera, J. Farrell, R. S. Perry, S. J. S. Lister, S. L. Lee, D. A. Tennant, Y. Maeno, and A. P. Mackenzie, Formation of a Nematic Fluid at High Fields in Sr3Ru2O7, Science 315, 214 (2007).
  * Hinkov  _et al._ [2008] V. Hinkov, D. Haug, B. Fauqué, P. Bourges, Y. Sidis, A. Ivanov, C. Bernhard, C. T. Lin, and B. Keimer, Electronic Liquid Crystal State in the High-Temperature Superconductor YBa2Cu3O6.45, Science 319, 597 (2008).
  * Chu  _et al._ [2010] J.-H. Chu, J. G. Analytis, K. D. Greve, P. L. McMahon, Z. Islam, Y. Yamamoto, and I. R. Fisher, In-Plane Resistivity Anisotropy in an Underdoped Iron Arsenide Superconductor, Science 329, 824 (2010).
  * Fradkin  _et al._ [2010] E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Nematic Fermi Fluids in Condensed Matter Physics, Annu. Rev. Condens. Matter Phys. 1, 153 (2010).
  * MacKinnon and Kramer [1981] A. MacKinnon and B. Kramer, One-Parameter Scaling of Localization Length and Conductance in Disordered Systems, Phys. Rev. Lett. 47, 1546 (1981).
  * MacKinnon and Kramer [1983] A. MacKinnon and B. Kramer, The scaling theory of electrons in disordered solids: Additional numerical results, Z. Physik B 53, 1 (1983).
  * Slevin and Ohtsuki [2014] K. Slevin and T. Ohtsuki, Critical exponent for the Anderson transition in the three-dimensional orthogonal universality class, New J. Phys. 16, 015012 (2014).
  * Markos [1995] P. Markos, Phenomenological theory of the metal-insulator transition, J. Phys.: Condens. Matter 7, 8361 (1995).
  * Mondragon-Shem  _et al._ [2014] I. Mondragon-Shem, T. L. Hughes, J. Song, and E. Prodan, Topological Criticality in the Chiral-Symmetric AIII Class at Strong Disorder, Phys. Rev. Lett. 113, 046802 (2014).
  * Altland  _et al._ [2014] A. Altland, D. Bagrets, L. Fritz, A. Kamenev, and H. Schmiedt, Quantum Criticality of Quasi-One-Dimensional Topological Anderson Insulators, Phys. Rev. Lett. 112, 206602 (2014).
  * Claes and Hughes [2020] J. Claes and T. L. Hughes, Disorder driven phase transitions in weak AIII topological insulators, Phys. Rev. B 101, 224201 (2020).
  * Molinari [2003] L. Molinari, Spectral duality and distribution of exponents for transfer matrices of block-tridiagonal Hamiltonians, J. Phys. A 36, 4081 (2003).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Localization Transitions in Non-Hermitian Quantum Mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * Asada  _et al._ [2004] Y. Asada, K. Slevin, and T. Ohtsuki, Numerical estimation of the $\beta$ function in two-dimensional systems with spin-orbit coupling, Phys. Rev. B 70, 035115 (2004).
  * Slevin and Ohtsuki [1999] K. Slevin and T. Ohtsuki, Corrections to Scaling at the Anderson Transition, Phys. Rev. Lett. 82, 382 (1999).
  * Wang  _et al._ [2021b] T. Wang, T. Ohtsuki, and R. Shindou, Universality classes of the Anderson transition in the three-dimensional symmetry classes AIII, BDI, C, D, and CI, Phys. Rev. B 104, 014206 (2021b).
  * Balents and Fisher [1997] L. Balents and M. P. A. Fisher, Delocalization transition via supersymmetry in one dimension, Phys. Rev. B 56, 12970 (1997).
  * Mathur [1997] H. Mathur, Feynman’s propagator applied to network models of localization, Phys. Rev. B 56, 15794 (1997).
  * Brouwer  _et al._ [1998] P. W. Brouwer, C. Mudry, B. D. Simons, and A. Altland, Delocalization in Coupled One-Dimensional Chains, Phys. Rev. Lett. 81, 862 (1998).
  * Konotop  _et al._ [2016] V. V. Konotop, J. Yang, and D. A. Zezyulin, Nonlinear waves in $\mathcal{PT}$-symmetric systems, Rev. Mod. Phys. 88, 035002 (2016).
  * Feng  _et al._ [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian photonics based on parity-time symmetry, Nat. Photon. 11, 752 (2017).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and $\mathcal{PT}$ symmetry, Nat. Phys. 14, 11 (2018).
  * [78] Z. Xiao, K. Kawabata, X. Luo, T. Ohtsuki, and R. Shindou, unpublished.
  * Crisanti  _et al._ [1993] A. Crisanti, G. Paladin, and A. Vulpiani,  _Products of Random Matrices_ (Springer, Berlin, Heidelberg, 1993).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and Topology in Non-Hermitian Physics, Phys. Rev. X 9, 041015 (2019).
  * Goldhirsch  _et al._ [1987] I. Goldhirsch, P.-L. Sulem, and S. A. Orszag, Stability and Lyapunov stability of dynamical systems: A differential approach and a numerical method, Physica D 27, 311 (1987).
  * Fulga  _et al._ [2012] I. C. Fulga, A. R. Akhmerov, J. Tworzydło, B. Béri, and C. W. J. Beenakker, Thermal metal-insulator transition in a helical topological superconductor, Phys. Rev. B 86, 054505 (2012).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Topological Insulators in Three Dimensions, Phys. Rev. Lett. 98, 106803 (2007).
  * Pendry  _et al._ [1990] J. B. Pendry, A. MacKinnon, and A. B. Pretre, Maximal fluctuations—A new phenomenon in disordered systems, Physica A 168, 400 (1990).
  * Kramer  _et al._ [2005] B. Kramer, T. Ohtsuki, and S. Kettemann, Random network models and quantum phase transitions in two dimensions, Phys. Rep. 417, 211 (2005).

##  I Supplemental Material for
“Topological Anderson Transitions in Chiral Symmetry Classes”

This Supplemental Material is organized as follows. In Sec. I.1, we summarize
known critical exponents between metal and topological-insulator phases and
between metal and ordinary-insulator phases. In Sec. I.2, we introduce the
inverse participation ratio along different directions and prove that coupling
among low-dimensional systems makes wave functions more extended even in the
small coupling limit. In Sec. I.3, we review the polynomial fitting of the
finite-size scaling function and details of Table I in the main text. In Sec.
I.4, we introduce the transfer matrix method and explain properties of
transfer matrices for chiral-symmetric Hamiltonians. In Sec. I.5, we summarize
a relation between weak topological indices and distributions of Lyapunov
exponents (LEs). In Sec. I.6, we show how statistical symmetries require the
weak topological indices to be zero and LEs of the non-Hermitian matrix (right
upper part of a chiral-symmetric Hamiltonian) to come in opposite-sign pairs.
In Sec. I.7, we summarize a scaling form for the maximal and minimal LEs
within a continuum spectrum and show numerical fittings based on this scaling
form. In Secs. I.8 and I.9, we provide detailed numerical studies of the
criticality in chiral-symmetric models with and without non-trivial
topological indices, respectively. In Sec. I.10, we provide detailed numerical
results of the two-terminal conductance of the model with non-trivial
topological indices. The conductance shows the anisotropic transport behavior
in the quasi-localized phase.

###### Contents

  1. I Supplemental Material for  “Topological Anderson Transitions in Chiral Symmetry Classes”
    1. I.1 Summary of known critical exponents
    2. I.2 Inverse participation ratio
      1. I.2.1 Inverse participation ratio along different directions
      2. I.2.2 Wave-function hybridization and inverse participation ratio
    3. I.3 Polynomial fitting
    4. I.4 Transfer matrix, Lyapunov exponents, and localization length
      1. I.4.1 Transfer matrix of a chiral-symmetric Hamiltonian
      2. I.4.2 Transfer matrix of the nodal-line semimetal model
    5. I.5 Weak topological indices and Lyapunov exponents
      1. I.5.1 Winding number in the clean limit
    6. I.6 Statistical symmetry
      1. I.6.1 Statistical symmetry, Lyapunov exponents, and one-dimensional winding number
      2. I.6.2 Statistical symmetry in the nodal-line semimetal model
    7. I.7 Finite-size scaling of Lyapunov exponents
      1. I.7.1 Numerical fitting
    8. I.8 Quasi-localized phases in chiral-symmetric models with weak topological indices
      1. I.8.1 Nodal-line semimetal in class AIII
      2. I.8.2 Nodal-line semimetals in class BDI
      3. I.8.3 Weak topological insulators and ordinary insulators in class BDI
    9. I.9 Anderson transitions in chiral-symmetric models with no weak topological indices
      1. I.9.1 Three-dimensional non-topological models in classes BDI and AIII
      2. I.9.2 Another three-dimensional non-topological model in class AIII
    10. I.10 Conductance and weak topological indices

###  I.1 Summary of known critical exponents

Table 2 summarizes known results of critical exponents of the Anderson
transitions between metal and topological-insulator phases and between metal
and trivial-insulator phases in the same spatial dimensions and symmetry
classes. For each symmetry class and spatial dimensions, the evaluated
critical exponents for the two types of the Anderson transitions are
consistent with each other. However, the critical exponents between
topological-insulator and trivial-insulator phases, as well as those between
topological-semimetal and diffusive-metal phases, can be different from the
ones in Table 2. In our work, we focus on the Anderson transitions in the 3D
chiral classes and demonstrate the different critical exponents due to the
weak topological indices.

Table 2: Correlation-length critical exponents of the Anderson transitions between metal and topological-insulator phases and those between metal and trivial-insulator phases in the same symmetry classes and spatial dimensions.  |  Class |  Metal-topological-insulator |  Metal-trivial-insulator
---|---|---|---
2D | AII |  $2.74\pm 0.12$ [18] |  $2.73\pm 0.02$11195% confidence interval. [15]
| D |  $1.371[1.311,1.437]$11footnotemark: 1 [32] |  $1.348[1.279,1.402]$11footnotemark: 1 [32]
| DIII |  $1.5\pm 0.1$ [28] |  $1.5\pm 0.1$ [28]
3D | A |  $1.34[1.23,1.53]$11footnotemark: 1 [27]222Layered Chern insulator. |  $1.443[1.437,1.449]$11footnotemark: 1 [25]
|  |  $1.42\pm 0.12$ [29]333Axion insulator. |
| AII |  $1.311\pm 0.033$11footnotemark: 1 [30] 444Reference [30] obtained the critical exponent between the metal and topological-insulator phases in a three-dimensional network model belonging to symmetry class AII. The critical exponent is different from the one between the metal and trivial-insulator phases in the same symmetry class and spatial dimensions obtained by the SU(2) model [16]. However, a more careful error analysis is needed, because system sizes in Ref. [30] may not be large enough ($\leq 10$), and the difference between the two exponents is small. |  $1.375\pm 0.016$11footnotemark: 1 [16]
| DIII |  $0.85\pm 0.05$ [26] |  $0.903[0.896,0.908]$11footnotemark: 1 [33]

###  I.2 Inverse participation ratio

####  I.2.1 Inverse participation ratio along different directions

The inverse participation ratio $P_{2}$ measures localization properties of a
wave function $\Phi(\bm{r})$ in $d$ dimension, defined by

$P_{2}\equiv\sum_{\bm{r}}|\Phi(\bm{r})|^{4}$ (9)

with the normalization condition $\sum_{\bm{r}}|\Phi(\bm{r})|^{2}=1$. We have
$P_{2}\leq 1$, where the equality holds only when $\Phi(\bm{r})$ is fully
localized at one lattice site. $P_{2}$ in extended and localized phases show
the different scaling relations with the system size $L$,

$P_{2}\sim\begin{cases}L^{-d}&[\Phi(\bm{r})\text{ is extended}]\,,\\\ {\rm
constant}<1&[\Phi(\bm{r})\text{ is localized}]\end{cases}$ (10)

for $L\to\infty$. The one-dimensional inverse participation ratio
$P^{\mu}_{2}$ ($\mu=x,y,z$) measures localization properties of $\Phi(x,y,z)$
along the $\mu$ direction. The integrated weight of the wave function at $z$,
$|\phi(z)|^{2}\equiv\sum_{x,y}|\Phi(x,y,z)|^{2}$, is regarded as the squared
one-dimensional normalized wave function along the $z$ direction and describes
how the three-dimensional wave function $\Phi(x,y,z)$ is localized along the
$z$ direction. Thus, the inverse participation ratio $P_{2}^{z}$ along the $z$
direction is introduced as

$P_{2}^{z}=\sum_{z}|\phi(z)|^{4}\,\quad,\quad|\phi(z)|^{2}=\sum_{x,y}|\Phi(x,y,z)|^{2}\,.$
(11)

Notably, $P_{2}^{z}$ provides an upper bound of $P_{2}$,

$P_{2}=\sum_{z}|\phi(z)|^{4}\left[\sum_{x,y}\left\lvert\frac{\Phi({\bm{r}})}{\phi(z)}\right\rvert^{4}\right]\leq
P_{2}^{z},$ (12)

where the equality holds only when we have
$\sum_{x,y}\left\lvert\Phi({\bm{r}})/\phi(z)\right\rvert^{4}=1$ for all $z$.
The one-dimensional inverse participation ratio $P_{2}^{\mu}$ along the other
two directions ($\mu=x,y$) is defined in the same manner. In a similar manner,
the $d$-dimensional inverse participation ratio is defined for a normalized
wave function in $d^{\prime}$ dimension ($d<d^{\prime}$). For example, the
following two-dimensional inverse participation ratio $P^{(x,y)}_{2}$ measures
the localization properties of $\Phi(x,y,z)$ within the $xy$ plane,

$P_{2}^{(x,y)}\equiv\sum_{x,y}|\phi(x,y)|^{4}\,,\quad|\phi(x,y)|^{2}=\sum_{z}|\Phi(x,y,z)|^{2}\,,$
(13)

which satisfies

$P_{2}\leq P_{2}^{(x,y)}\,.$ (14)

####  I.2.2 Wave-function hybridization and inverse participation ratio

Suppose that a $d^{\prime}$-dimensional disordered non-interacting Hamiltonian
${\cal H}^{\bm{R}}$ in ${\bm{R}}\equiv({\bm{r}},{\bm{s}})$ with
${\bm{r}}=(r_{1},\cdots,r_{d})$ and ${\bm{s}}=(s_{1},\cdots,s_{d^{\prime}-d})$
($d<d^{\prime}$) consists of $d$-dimensional Hamiltonians ${\cal
H}^{\bm{r}}_{\bm{s}}$ at different ${\bm{s}}$ and coupling ${\cal H}^{\prime}$
among the $d$-dimensional systems. On-site disorder potential $V({\bm{R}})$ is
chosen to distribute uniformly in the range $[-W/2,W/2]$ for all the lattice
sites ${\bm{R}}$. Then, ${\bm{s}}$ can be regarded as different disorder
realizations from the same ensemble for the $d$-dimensional system with the
disorder strength $W$. In this section, we show that even in the small
coupling limit, an eigenstate $\Phi({\bm{R}})$ of ${\cal H}^{\bm{R}}$ with
eigenenergy $E$ is more extended along the ${\bm{r}}$ direction than an
eigenstate $\psi({\bm{r}})$ of ${\cal H}^{\bm{r}}_{\bm{s}}$ with the same
eigenenergy. Here, the small coupling limit means that the maximal eigenvalue
of ${\cal H}^{\prime}$ is much smaller than the mean level spacing of ${\cal
H}^{\bm{r}}_{\bm{s}}$ around $E$.

For ${\cal H}^{\prime}=0$, eigenstates of ${\cal H}^{\bm{R}}$ are given by
eigenstates of ${\cal H}_{\bm{s}}^{\bm{r}}$. In the small coupling limit, we
can treat ${\cal H}^{\prime}$ perturbatively. We introduce an energy window
$[E-\Delta E,E+\Delta E]$ and choose $\Delta E$ to be small enough that each
${\cal H}^{\bm{r}}_{\bm{s}}$ has at most one eigenstate
$\psi_{\bm{s}}({\bm{r}})$ with eigenenergy $E_{\bm{s}}$ in the energy window
and that we have $|E_{\bm{s}}-E|\ll\Delta E$. In the small coupling limit, the
maximal eigenvalue of ${\cal H}^{\prime}$ can be much smaller than $\Delta E$.
Thus, in the lowest order of degenerate perturbation theory, ${\cal
H}^{\prime}$ does not mix unperturbed eigenstates inside the energy window
with those outside the energy window, and an eigenstate $\Phi({\bm{R}})$ of
${\cal H}^{\bm{R}}$ is given by a linear superposition of
$\psi_{\bm{s}}({\bm{r}})$ over different ${\bm{s}}$,

$\Phi({\bm{r}},{\bm{s}})=\sum_{{\bm{s}}^{\prime}}a_{{\bm{s}}^{\prime}}\psi_{\bm{s}^{\prime}}({\bm{r}})\delta_{{\bm{s}},{\bm{s}}^{\prime}}=a_{{\bm{s}}}\psi_{\bm{s}}({\bm{r}})\,,$
(15)

where $\delta_{\bm{s},\bm{s}^{\prime}}$ is the Kronecker delta. Here, we
impose the normalization conditions $\sum_{\bm{s}}|a_{\bm{s}}|^{2}=1$ and
$\sum_{\bm{r}}|\psi_{\bm{s}}({\bm{r}})|^{2}=1$, where we sum only over such
${\bm{s}}$ that ${\cal H}^{\bm{r}}_{\bm{s}}$ has an eigenenergy inside the
window $[E-\Delta E,E+\Delta E]$. $a_{\bm{s}}$ is the ${\bm{s}}$-component of
an eigenstate of an effective Hamiltonian $H^{\bm{R}}$ given as

$\displaystyle({H}^{\bm{R}})_{{\bm{s}},{\bm{s}}^{\prime}}=\sum_{\bm{r}}\sum_{{\bm{r}}^{\prime}}\\!\
\psi^{*}_{\bm{s}}({\bm{r}})\\!\ ({\cal
H}^{\prime})_{{\bm{r}},{\bm{r}}^{\prime}}\\!\
\psi_{{\bm{s}}^{\prime}}({\bm{r}}^{\prime})+E_{s}\delta_{\bm{s},\bm{s}^{\prime}}.$
(16)

In the following, we show that the inverse participation ratio
$P_{2}^{\bm{r}}$ of $\Phi({\bm{R}})$ along the ${\bm{r}}$ direction is always
smaller than the $d$-dimensional inverse participation ratio of
$\psi_{\bm{s}}({\bm{r}})$. The weight $|\phi({\bm{r}}_{0})|^{2}$ of the wave
function $\Phi({\bm{R}})$ on a hyperplane ${\bm{r}}={\bm{r}}_{0}$ is given as

$|\phi({\bm{r}})|^{2}=\sum_{\bm{s}}\left\lvert\Phi({\bm{r}},{\bm{s}})\right\rvert^{2}\,.$
(17)

The inverse participation ratio $P_{2}^{\bm{r}}$ of $\Phi({\bm{R}})$ along the
${\bm{r}}$ direction measures the localization properties of $\Phi({\bm{R}})$
within the ${\bm{r}}$ direction and is given by the sum of the square of the
weight over ${\bm{r}}$,

$\displaystyle P_{2}^{\bm{r}}$
$\displaystyle=\sum_{\bm{r}}|\phi({\bm{r}})|^{4}=\sum_{\bm{r}}\sum_{\bm{s}_{1}}\sum_{\bm{s}_{2}}\left\lvert\Phi({\bm{r}},{\bm{s}_{1}})\right\rvert^{2}\left\lvert\Phi({\bm{r}},{\bm{s}_{2}})\right\rvert^{2}=\sum_{\bm{s}_{1}}\sum_{\bm{s}_{2}}|a_{{\bm{s}_{1}}}|^{2}|a_{{\bm{s}_{2}}}|^{2}\left[\sum_{\bm{r}}|\psi_{\bm{s}_{1}}({\bm{r}})|^{2}|\psi_{\bm{s}_{2}}({\bm{r}})|^{2}\right]$
(18)
$\displaystyle\leq\frac{1}{2}\sum_{\bm{s}_{1}}\sum_{\bm{s}_{2}}|a_{{\bm{s}_{1}}}|^{2}|a_{{\bm{s}_{2}}}|^{2}\sum_{\bm{r}}\left[|\psi_{\bm{s}_{1}}({\bm{r}})|^{4}+|\psi_{\bm{s}_{2}}({\bm{r}})|^{4}\right]=\sum_{\bm{s}}|a_{{\bm{s}}}|^{2}\sum_{\bm{r}}|\psi_{\bm{s}}({\bm{r}})|^{4}\,.$

Here, the equality holds true only when we have
$\psi_{\bm{s}_{1}}({\bm{r}})=\psi_{\bm{s}_{2}}({\bm{r}})$ for all $\bm{r}$,
${\bm{s}}_{1}$, and ${\bm{s}}_{2}$. Notably,
$\sum_{\bm{r}}|\psi_{\bm{s}}({\bm{r}})|^{4}$ is the $d$-dimensional inverse
participation ratio $P_{2}^{\psi_{\bm{s}}(\bm{r})}$ of
$\psi_{\bm{s}}(\bm{r})$, and ${\cal H}^{\bm{r}}_{\bm{s}}$ at different
${\bm{s}}$ belongs to the same ensemble with the same disorder strength $W$.
In the thermodynamic limit ($N_{\bm{r}}\equiv\sum_{\bm{r}}\rightarrow\infty$),
$\Delta E$ goes to zero as the mean level spacing goes to zero, and
$P_{2}^{\psi_{\bm{s}}(\bm{r})}$ at different ${\bm{s}}$ takes the same value
$P_{2}^{\psi_{(}\bm{r})}$. Then, we have

$\displaystyle P^{\bm{r}}_{2}\leq
P^{\psi({\bm{r}})}_{2}\sum_{\bm{s}}|a_{{\bm{s}}}|^{2}=P^{\psi({\bm{r}})}_{2},$
(19)

which proves that within the lowest order in ${\cal H}^{\prime}$, the small
coupling ${\cal H}^{\prime}$ among the $d$-dimensional systems always makes
$d$-dimensional wave functions spatially more extended.

###  I.3 Polynomial fitting

In this section, we present more details about the polynomial fitting [Eq. (7)
in the main text] of the normalized localization length
$\Lambda_{x}(W,L)=\xi_{x}(W,L)/L$ and show details of Table I in the main text
(see Table 3). The scaling function for $\Lambda_{x}(W,L)$ is Taylor-expanded
with respect to the relevant scaling variable $\phi(w)$ and the least
irrelevant scaling variable $\psi(w)$ up to the $n$th order and first order,
respectively,

$\displaystyle\Lambda_{x}(W,L)=\sum^{n}_{i=0}\sum^{1}_{j=0}a_{i,j}\big{(}\phi(w)L^{1/\nu}\big{)}^{i}\big{(}\psi(w)L^{-y}\big{)}^{j},$
(20)

with $w\equiv(W-W^{(x)}_{c})/W^{(x)}_{c}$, and the scaling dimension $-y$
$(<0)$ of the least irrelevant scaling variable around a saddle-point fixed
point. The relevant scaling variable is further expanded around $w=0$ up to
the $m$th order, while only the zeroth-order in $w$ is kept for the irrelevant
scaling variable $\psi(w)$,

$\displaystyle\phi(w)=\sum^{m}_{k=1}b_{k}w^{k},\quad\psi(w)=c.$ (21)

Here, $\\{W^{(x)}_{c},\nu,y,a_{i,j},b_{k},c\\}$ are the fitting parameters. To
avoid the ambiguity in the Taylor expansion of the scaling function, we should
set $a_{0,1}=a_{1,0}=1$. Thus, the number $N_{f}$ of the free parameters in
the fitting is $N_{f}=2(n+1)+m+2$. We minimize $\chi^{2}$ statistics

$\chi^{2}=\sum_{j=1}^{N_{D}}\left(\frac{F_{j}-\Lambda_{j}}{\sigma_{j}}\right)^{2}\,,$
(22)

where $\Lambda_{j}$ and $\sigma_{j}$ are the normalized localization length
and its standard deviation for $(W,L)$ evaluated by the transfer matrix
method, respectively, $F_{j}$ is the value of the polynomial fitting function
for $(W,L)$, and $N_{D}$ is the number of data points. The confidence error
bars for the optimal parameters are determined by the fittings for 1000 sets
of synthetic data for $\Lambda_{x}(W,L)$. The synthetic data are generated
according to a standard deviation from the transfer matrix calculation.

Table 3: Polynomial fitting results of the normalized localization length $\Lambda_{\mu}\equiv\xi_{\mu}/L$ along the $\mu$ direction ($\mu=x,y,z$) around critical points of different models with the quasi-one-dimensional geometry $L\times L\times L_{\mu}$. “$\surd$” in the column “topology” shows that the weak topological index $\nu_{z}$ is non-zero around the critical point, and “$\times$” shows that all the weak topological indices always vanish around the critical point. We show the critical disorder strength $W^{(\mu)}_{c}$, critical exponent $\nu$, scaling dimension $-y$ of the least irrelevant scaling variable, critical localization length $\Lambda_{c}$, the goodness of fitting (GOF), and Taylor-expansion order of $(m,n)$ in Eqs. (20) and (21). The square brackets denote the 95% confidence interval. Note that $\Lambda_{c}$’s here are critical values in the presence of anisotropic spatial geometry and take non-universal values. symmetry class | topology | direction | $m$ | $n$ | GOF | $W^{(\mu)}_{c}$ | $\nu$ | $y$ | $\Lambda_{c}$
---|---|---|---|---|---|---|---|---|---
BDI | $\surd$ | $\mu=x$ | 2 | 3 | 0.15 | 27.241[27.194,27.303] | 0.820[0.783,0.846] | 2.584[2.175,2.955] | 0.134[0.130,0.137]
BDI | $\surd$ | $\mu=x$ | 3 | 3 | 0.14 | 27.243[27.192,27.301] | 0.820[0.787,0.848] | 2.574[2.212,2.947] | 0.134[0.130,0.138]
AIII | $\surd$ | $\mu=x$ | 3 | 3 | 0.47 | 9.143[9.125,9.168] | 0.824[0.776,0.862] | 2.157[1.727,2.519] | 0.225[0.213,0.232]
BDI | $\times$ | $\mu=z$ | 2 | 3 | 0.19 | 23.220[23.167,23.293] | 1.089[1.005,1.128] | 1.926[1.074,3.034] | 0.374[0.352,0.385]
BDI | $\times$ | $\mu=z$ | 3 | 3 | 0.18 | 23.223[23.138,23.409] | 1.088[0.991,1.141] | 1.906[0.604,3.677] | 0.373[0.302,0.389]
BDI | $\times$ | $\mu=x$ | 2 | 3 | 0.23 | 23.170[23.098,23.279] | 1.042[0.943,1.099] | 1.591[0.889,2.543] | 0.281[0.254,0.293]
BDI | $\times$ | $\mu=x$ | 3 | 3 | 0.31 | 23.167[23.101,23.310] | 1.039[0.937,1.100] | 1.607[0.753,2.425] | 0.281[0.239,0.292]
AIII | $\times$ | $\mu=z$ | 2 | 3 | 0.20 | 8.091[8.074,8.096] | 1.024[0.973,1.070] | 0.470[0.450,1.481] | 0.650[0.639,0.706]

###  I.4 Transfer matrix, Lyapunov exponents, and localization length

The transfer matrix method solves an eigenvalue problem of a non-interacting
disordered Hamiltonian $\cal{H}$ recursively. This method is efficient for
obtaining the localization length along one spatial direction, which we call
the $\mu$ direction in the following. In this formulation, the Hamiltonian is
decomposed into a layer structure along the $\mu$ direction,

$\mathcal{H}_{i,j}=H_{i}\delta_{i,j}+V_{i,i+1}\delta_{i,j-1}+V_{i,i-1}\delta_{i,j+1}\,,$
(23)

where $i,j=1,2,\cdots,L_{\mu}$ are indices of the layers, $H_{i}$ is a block
of matrix elements within the $i$th layer, and $V_{i,i\pm 1}$ is a block of
matrix elements between the $i$th layer and the $(i\pm 1)$th layer. The
decomposition assumes that matrix elements appear only between the nearest
neighboring layers or within each layer. In the presence of next-nearest
hopping, one can redefine two neighboring layers as one layer. Let $H_{i}$,
$V_{i,i\pm 1}$ be $m$ by $m$ matrices and
$(\cdots,A_{i-1},A_{i},A_{i+1},\cdots)^{T}$ be an eigenvector of ${\cal H}$
for an eigenenergy $E$:

$H_{i}A_{i}+V_{i,i-1}A_{i-1}+V_{i,i+1}A_{i+1}=EA_{i}.$ (24)

For simplicity, suppose that the disorder terms are present only in the
diagonal matrix elements and that $V_{i,i-1}=V_{+}$ and $V_{i,i+1}=V_{-}$ are
free from disorder. The eigenvectors are solved layer by layer recursively by
a transfer matrix $M_{i}$,

$\left(\begin{matrix}A_{i+1}\\\
A_{i}\end{matrix}\right)=M_{i}\left(\begin{matrix}A_{i}\\\
A_{i-1}\end{matrix}\right),\\!\ \\!\
\,M_{i}\equiv\left(\begin{matrix}-V_{-}^{-1}(H_{i}-E)&-V_{-}^{-1}V_{+}\\\
1_{m\times m}&0_{m\times m}\\\ \end{matrix}\right)\,.$ (25)

The product of the transfer matrices, $M=M_{L_{\mu}}M_{L_{\mu}-1}\cdots
M_{1}$, relates the components of the eigenvector at the $(L_{\mu}+1)$th and
$L_{\mu}$th layers with the components at the first and zeroth layer.
According to Oseledec’s theorem [79], the matrix

$P(E)=\lim_{L_{\mu}\rightarrow\infty}\ln{(M^{\dagger}M)^{\frac{1}{2L_{\mu}}}}=-\lim_{L_{\mu}\rightarrow\infty}\ln{(M^{-1}M^{-1\dagger})^{\frac{1}{2L_{\mu}}}}$
(26)

well converges in the limit $L_{\mu}\to\infty$. Eigenvalues of $P(E)$ are
known as Lyapunov exponents (LEs). If $\mathcal{H}$ is Hermitian, LEs come in
opposite-sign pairs [79]. The inverse of the smallest positive or the largest
negative LE corresponds to the localization length $\xi_{\mu}$ along the $\mu$
direction.

####  I.4.1 Transfer matrix of a chiral-symmetric Hamiltonian

Suppose that a $2n\times 2n$ Hermitian Hamiltonian ${\cal H}$ satisfies chiral
symmetry $\mathcal{C}{\cal H}\mathcal{C}^{-1}=-{\cal H}$ with a chiral
operator $\mathcal{C}$ satisfying $\mathcal{C}^{2}=1$. Eigenvalues of
$\mathcal{C}$ are $\pm 1$, the numbers of which are assumed to be the same.
Then, the unitary matrix ${\cal C}$ is diagonalized as
$\mathcal{C}=\sum_{i=1}^{n}\ket{v_{i}}\bra{v_{i}}-\sum_{i=1}^{n}\ket{u_{i}}\bra{u_{i}}$.
Here, $\ket{v_{1}},\cdots,\ket{v_{n}}$ and $\ket{u_{1}},\cdots,\ket{u_{n}}$
are eigenvectors of $\mathcal{C}$ with eigenvalues $+1$ and $-1$,
respectively. Because of chiral symmetry, we have $\braket{v_{i}}{{\cal
H}}{v_{j}}=\braket{u_{i}}{{\cal H}}{u_{j}}=0$. Thus, the $2n\times 2n$ matrix
${\cal H}$ is decomposed into two $n\times n$ matrices $h$ and $h^{\prime}$ in
the off-diagonal parts,

$\displaystyle{\cal H}=\left(\begin{array}[]{cc}0&h\\\ h^{\prime}&0\\\
\end{array}\right),$ (29)

with

$\displaystyle(h)_{i,j}=\braket{v_{i}}{{\cal
H}}{u_{j}},\quad(h^{\prime})_{i,j}=\braket{u_{i}}{{H}}{v_{j}},$ (30)

satisfying ${h}^{\prime}={h}^{\dagger}$.

Equation (30) does not determine ${h}$ uniquely up to $n\times n$ unitary
transformations, $h\rightarrow{\cal V}^{\dagger}{h}{\cal U}$, where the
unitary transformations ${\cal V}$ and ${\cal U}$ change bases among the
$n$-fold degenerate eigenstates of $\mathcal{C}$. Nonetheless, with a certain
choice of the bases for $\ket{v_{1}},\cdots,\ket{v_{n}}$ and
$\ket{u_{1}},\cdots,\ket{u_{n}}$, any Hermitian Hamiltonian with chiral
symmetry can be decomposed into the off-diagonal form as Eq. (29). The off-
diagonal parts thus introduced are non-Hermitian matrices, in general.

If the chiral operator ${\cal C}$ is diagonal with respect to the layer index
and its matrix elements do not depend on the layer index, the Hamiltonian
$H_{i}$ within the $i$th layer and the hopping matrix $V_{\pm}$ between the
$i$th layer and the $(i\pm 1)$th layers also take the block off-diagonal
structure,

$H_{i}=\begin{pmatrix}0&\tilde{h}_{i}\\\
\tilde{h}_{i}^{\dagger}&0\end{pmatrix},\quad V_{+}=\begin{pmatrix}0&v_{+}\\\
v_{-}^{\dagger}&0\end{pmatrix},\quad V_{-}=\begin{pmatrix}0&v_{-}\\\
v_{+}^{\dagger}&0\end{pmatrix},$ (31)

where $v_{\pm}$ are free from the disorder and independent of the layer index.
The transfer matrix $M^{(\cal H)}_{i}$ of ${\cal H}$ for zero energy $E=0$
reads

$\displaystyle M^{(\cal H)}_{i}$
$\displaystyle=\begin{pmatrix}-\begin{pmatrix}0&v_{-}\\\
v_{+}^{\dagger}&0\end{pmatrix}^{-1}\begin{pmatrix}0&\tilde{h}_{i}\\\
\tilde{h}_{i}^{\dagger}&0\end{pmatrix}&-\begin{pmatrix}0&v_{-}\\\
v_{+}^{\dagger}&0\end{pmatrix}^{-1}\begin{pmatrix}0&v_{+}\\\
v_{-}^{\dagger}&0\end{pmatrix}\\\ \begin{matrix}1&0\\\ 0&1\\\
\end{matrix}&\begin{matrix}0&0\\\ 0&0\\\ \end{matrix}\end{pmatrix}$ (32)
$\displaystyle=\begin{pmatrix}-v_{+}^{\dagger-1}\tilde{h}_{i}^{\dagger}&0&-v_{+}^{\dagger-1}v_{-}^{\dagger}&0\\\
0&-v_{-}^{-1}\tilde{h}_{i}&0&-v_{-}^{-1}v_{+}\\\ 1&0&0&0\\\ 0&1&0&0\\\
\end{pmatrix}\,.$

With a proper unitary transformation ${\cal U}$, the transfer matrix
$M_{i}^{({\cal H})}$ is block-diagonalized into $M_{i}$ and $M^{\prime}_{i}$,

$\displaystyle{\cal U}^{\dagger}M_{i}^{(\cal H)}{\cal
U}=\begin{pmatrix}M_{i}^{\prime}&0\\\ 0&M_{i}\\\ \end{pmatrix},\\!\ \\!\
M_{i}^{\prime}\equiv\begin{pmatrix}-v_{+}^{\dagger-1}\tilde{h}_{i}^{\dagger}&-v_{+}^{\dagger-1}v_{-}^{\dagger}\\\
1&0\\\ \end{pmatrix},\\!\ \\!\
M_{i}\equiv\begin{pmatrix}-v_{-}^{-1}\tilde{h}_{i}&-v_{-}^{-1}v_{+}\\\ 1&0\\\
\end{pmatrix}.$ (33)

Notably, $M_{i}$ and $M^{\prime}_{i}$ are the transfer matrices of the right-
upper part $h$ and left-lower part $h^{\dagger}$ of the Hamiltonian ${\cal H}$
in Eq. (29), respectively. In the canonical basis of Eqs. (29) and (31), $h$
and $h^{\prime}=h^{{\dagger}}$ are decomposed into the layer structure along
the $\mu$ direction,

$\displaystyle h_{i,j}$
$\displaystyle=\tilde{h}_{i}\delta_{i,j}+v_{+}\delta_{i,j-1}+v_{-}\delta_{i,j+1},$
$\displaystyle h^{\prime}_{i,j}$
$\displaystyle=\tilde{h}^{\dagger}_{i}\delta_{i,j}+v^{\dagger}_{-}\delta_{i,j-1}+v^{\dagger}_{+}\delta_{i,j+1},$

where $i,j=1,\cdots,L_{\mu}$ are the indices of layers. From Eqs. (24) and
(25), we obtain $M_{i}$ of $h$ and $M^{\prime}_{i}$ of $h^{\prime}$ as in Eq.
(33). Note that $M_{i}^{\prime}$ is equivalent to $(M_{i}^{\dagger})^{-1}$
under a certain transformation,

$S\equiv\begin{pmatrix}0&-v_{-}^{\dagger}\\\ v_{+}^{\dagger}&0\\\
\end{pmatrix},\quad
SM_{i}^{\prime}S^{-1}=(M_{i}^{\dagger})^{-1}=\begin{pmatrix}0&-v_{-}^{\dagger}v_{+}^{\dagger-1}\\\
1&-\tilde{h}_{i}^{\dagger}v_{+}^{\dagger-1}\end{pmatrix}\,.$ (34)

Thus, the LEs obtained by the product of $SM_{i}^{\prime}S^{-1}$ have signs
opposite to the LEs obtained by the product of $M_{i}$. The non-singular
similarity transformation $S$ does not change LEs. Thereby, the LEs of $h$ and
the LEs of $h^{\prime}=h^{\dagger}$ come in opposite-sign pairs. The LEs of
${\cal H}$ are the sum of the LEs of $h$ and the LEs of $h^{\prime}$.

####  I.4.2 Transfer matrix of the nodal-line semimetal model

The Hamiltonian of the nodal-line semimetal model reads

${\cal
H}=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left\\{(\Delta+\epsilon_{\bm{r}})c^{\dagger}_{\bm{r}}\sigma_{z}c_{\bm{r}}+\left[\sum_{\mu=x,y}\left(t_{\perp}c^{\dagger}_{\bm{r}+\bm{e_{\mu}}}\sigma_{z}c_{\bm{r}}\right)-{\rm
i}t_{\|}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{y}c_{\bm{r}}+t_{\|}^{\prime}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{z}c_{\bm{r}}+\text{H.c.}\right]\right\\}\,,$
(35)

where $c_{\bm{r}}$ is a two-component annihilation operator on the cubic
lattice site ${\bm{r}}$, $\sigma_{\mu}(\mu=x,y,z)$ are the Pauli matrices,
$\Delta$, $t_{\perp}$, $t_{\|}$, and $t_{\|}^{\prime}$ are real-valued
parameters, $\epsilon_{\bm{r}}$ is a random potential that distributes
uniformly in $[-W/2,W/2]$, and ${\bm{e}}_{x}=(1,0,0)$, ${\bm{e}}_{y}=(0,1,0)$,
and ${\bm{e}}_{z}=(0,0,1)$ are the unit vectors. Note that the Hamiltonian in
Eq. (35) reduces to Eq. (5) in the main text for $\Delta=0$. Depending on
$\Delta$ and the other parameters, Eq. (35) describes an ordinary insulator,
topological insulator, and nodal-line semimetal [see Eq. (64)]. ${\cal H}$
satisfies time-reversal symmetry ${\cal H}={\cal H}^{*}$ and chiral symmetry
${\cal H}=-\mathcal{C}^{\dagger}{\cal H}^{\dagger}\mathcal{C}$ with a unitary
operator
$\mathcal{C}_{{\bm{r},\bm{r}^{\prime}}}=\delta_{{\bm{r},\bm{r}^{\prime}}}\sigma_{x}$
with $\mathcal{C}^{*}=\mathcal{C}$, and thus belongs to the chiral orthogonal
class (class BDI).

The chiral operator ${\cal C}$ has eigenvalues $+1$ and $-1$. Since ${\cal C}$
is diagonal with respect to the lattice site, eigenvectors of ${\cal C}$ can
be labelled by the cubic-lattice site ${\bm{s}}\equiv(s_{x},s_{y},s_{z})$:

$\displaystyle\langle{\bm{r}}|v_{\bm{s}}\rangle=\delta_{{\bm{r}},{\bm{s}}}\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
1\\\
\end{array}\right),\quad\langle{\bm{r}}|u_{\bm{s}}\rangle=\delta_{{\bm{r}},{\bm{s}}}\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-1\\\ \end{array}\right),{\color[rgb]{0,0,1}}$ (40)

satisfying ${\cal C}|v_{\bm{s}}\rangle=|v_{\bm{s}}\rangle$ and ${\cal
C}|u_{\bm{s}}\rangle=-|u_{\bm{s}}\rangle$. Following Eq. (30), we construct
the right-upper part $h$ of $\cal H$ on the same cubic lattice as,

$\displaystyle(h)_{{\bm{s}},{\bm{s}}}$
$\displaystyle=\braket{v_{\bm{s}}}{{\cal
H}}{u_{\bm{s}}}=\Delta+\epsilon_{\bm{s}},$ (41)
$\displaystyle(h)_{{\bm{s}}+{\bm{e}}_{\mu},{\bm{s}}}=(h)_{{\bm{s}},{\bm{s}}+{\bm{e}}_{\mu}}$
$\displaystyle=\braket{v_{{\bm{s}}+{\bm{e}}_{\mu}}}{{\cal
H}}{u_{\bm{s}}}=\braket{v_{\bm{s}}}{{\cal
H}}{u_{{\bm{s}}+{\bm{e}}_{\mu}}}=t_{\perp},$ (42)
$\displaystyle(h)_{{\bm{s}}+{\bm{e}}_{z},{\bm{s}}}$
$\displaystyle=\braket{v_{{\bm{s}}+{\bm{e}}_{z}}}{{\cal
H}}{u_{\bm{s}}}=t_{\|}+t_{\|}^{\prime},$ (43)
$\displaystyle(h)_{{\bm{s}},{\bm{s}}+{\bm{e}}_{z}}$
$\displaystyle=\braket{v_{\bm{s}}}{{\cal
H}}{u_{{\bm{s}}+{\bm{e}}_{z}}}=-t_{\|}+t_{\|}^{\prime},$ (44)

for $\mu=x,y$. All the other matrix elements of $h$ are zero. Notably, $h$ can
be regarded as a single-orbital tight-binding model,

$\displaystyle h$
$\displaystyle=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left[(\Delta+\epsilon_{\bm{r}})f^{\dagger}_{\bm{r}}f_{\bm{r}}+\sum_{\mu=x,y}\left(t_{\perp}f^{\dagger}_{{\bm{r}}+{\bm{e}}_{\mu}}f_{\bm{r}}+\text{H.c.}\right)+(t^{\prime}_{\|}+t_{\|})f^{\dagger}_{{\bm{r}}+{\bm{e}}_{z}}f_{\bm{r}}+(t^{\prime}_{\|}-t_{\|})f^{\dagger}_{\bm{r}}f_{{\bm{r}}+{\bm{e}}_{z}}\right]\,,$
(45)

where $f_{\bm{r}}$ and $f^{\dagger}_{\bm{r}}$ are annihilation and creation
operators at site ${\bm{r}}$. While $h$ respects $h=h^{*}$, we have $h\neq
h^{\dagger}$ for $t_{\|}\neq 0$. Hence, $h$ generally belongs to the non-
Hermitian symmetry class AI [80, 33].

The transfer matrix of $h$ along the $z$ direction is given by

$M_{i}=\left(\begin{matrix}-\frac{1}{t^{\prime}_{\|}-t_{\|}}\tilde{h}_{i}&-\frac{t^{\prime}_{\|}+t_{\|}}{t^{\prime}_{\|}-t_{\|}}1_{m\times
m}\\\ 1_{m\times m}&0_{m\times m}\\\ \end{matrix}\right)\,,$ (46)

where $m=L^{2}$ is the degrees of freedom in each layer and the quasi-1D
geometry ($L\times L\times L_{z}$, $L_{z}\gg L$) is considered.
$\tilde{h}_{i}$ is the Hamiltonian within the $i$th layer, which has
$\Delta+\epsilon_{\bm{s}}$ in its diagonal elements and $t_{\perp}$ in its
nearest-neighbor hopping.

For $t^{\prime}_{\|}-t_{\|}=0$, $M_{i}$ is singular, and $m$ eigenvalues of
$\frac{1}{L_{z}}\ln M\equiv\frac{1}{L_{z}}\ln(M_{L_{z}}M_{L_{z}-1}\cdots
M_{1})$ diverge to $\infty$. In fact, $M^{-1}_{i}$ has zero eigenvalues with
multiplicity $m$ for $t^{\prime}_{\|}=t_{\|}$,

$M_{i}^{-1}=\begin{pmatrix}0_{m\times m}&1_{m\times m}\\\
-\frac{t^{\prime}_{\|}-t_{\|}}{t^{\prime}_{\|}+t_{\|}}1_{m\times
m}&-\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{i}\\\
\end{pmatrix}\,\rightarrow\begin{pmatrix}0_{m\times m}&1_{m\times m}\\\
0_{m\times m}&-\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{i}\\\
\end{pmatrix}\quad(t^{\prime}_{\|}-t_{\|}\to 0),$ (47)

where $M^{-1}=M^{-1}_{1}\cdots M^{-1}_{L_{z}-1}M^{-1}_{L_{z}}$ has zero
eigenvalue with multiplicity at least $m$. Therefore, $m$ eigenvalues of
$\frac{1}{2L_{z}}\ln(M^{-1})^{\dagger}M^{-1}$ always diverge to $-\infty$,
while $m$ eigenvalues of $\frac{1}{2L_{z}}\ln M^{\dagger}M$ always diverge to
$+\infty$. The other $m$ finite-valued LEs of $h$ are determined from the
following product:

$\displaystyle p$
$\displaystyle\equiv-\lim_{L_{z}\rightarrow\infty}\frac{1}{2L_{z}}\ln{\left(\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{1}\cdots\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{L_{z}}\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{L_{z}}^{\dagger}\cdots\frac{1}{t^{\prime}_{\|}+t_{\|}}\tilde{h}_{1}^{\dagger}\right)},$
$\displaystyle=\lim_{L_{z}\rightarrow\infty}\frac{1}{2L_{z}}\ln{\left((t^{\prime}_{\|}+t_{\|})\tilde{h}_{1}^{-1\dagger}\cdots(t^{\prime}_{\|}+t_{\|})\tilde{h}_{L_{z}}^{-1\dagger}(t^{\prime}_{\|}+t_{\|})\tilde{h}_{L_{z}}^{-1}\cdots(t^{\prime}_{\|}+t_{\|})\tilde{h}_{1}^{-1}\right)}.$
(48)

###  I.5 Weak topological indices and Lyapunov exponents

We summarize a relationship between the weak topological indices of chiral-
symmetric Hamiltonians and the numbers of positive and negative LEs of its
right-upper part $h$ in the canonical basis in Eq. (29) [66]. Consider a
chiral-symmetric Hamiltonian ${\cal H}(\phi)$, in which a magnetic flux $\phi$
is inserted through a closed loop along the $\mu$ direction. Similarly to
${\cal H}$ in Eq. (29), ${\cal H}(\phi)$ takes a block off-diagonal structure
in a basis where the chiral operator is diagonal,

$\displaystyle{\cal H}(\phi)=\begin{pmatrix}0&h(\phi)\\\
h(\phi)^{\dagger}&0\\\ \end{pmatrix}.$ (49)

The right-upper block $h(\phi)$ is decomposed into a layer structure along the
$\mu$ direction,

$h(\phi)=\begin{pmatrix}\tilde{h}_{1}&v_{-}&0&\cdots&0&\frac{1}{z}v_{+}\\\
v_{+}&\tilde{h}_{2}&v_{-}&\cdots&\cdots&0\\\
0&v_{+}&\tilde{h}_{3}&v_{-}&\cdots&0\\\
\vdots&\vdots&\ddots&\ddots&\ddots&\vdots\\\
zv_{-}&0&\cdots&0&v_{+}&\tilde{h}_{L_{\mu}}\\\ \end{pmatrix}\,,$ (50)

with $z=e^{{\rm i}\phi}$. We assume that the hopping appears only between the
nearest neighboring layers or within each layer. In the presence of next-
nearest neighbor hopping, we can redefine two neighboring layers as one layer.
$v_{-}$, $v_{+}$, and $\tilde{h}_{i}(i=1,2,\cdots,L_{\mu})$ are $m\times m$
matrices, where $m$ is the degrees of freedom of $h$ in each layer. The 1D
winding number $w_{\mu}$ along the $\mu$ direction is defined in terms of
$h(\phi)$,

$\displaystyle w_{\mu}\equiv{\rm
i}\int^{2\pi}_{0}\frac{d\phi}{2\pi}\partial_{\phi}{\rm
Tr}\big{[}\log\big{[}h(\phi)]\big{]}\,.$ (51)

The winding number $w_{\mu}$ is given by the contour integral

$\displaystyle w_{\mu}$ $\displaystyle={\rm
i}\oint_{|z|=1}\frac{dz}{2\pi}\partial_{z}{\rm
Tr}\big{[}\log\big{[}h(z)]\big{]}$ (52) $\displaystyle={\rm
i}\oint_{|z|=1}\frac{dz}{2\pi}\partial_{z}\log\left[\det\left[h(z)\right]\right].$

Here, $\det\left[h(z)\right]$ is a polynomial function in terms of $z$ with
the lowest order $z^{-m}$ and the highest order $z^{m}$ if we have
$\det{v_{-}}\neq 0$ and $\det{v_{+}}\neq 0$. For simplicity, we assume
$\det{v_{-}}\neq 0$ and $\det{v_{+}}\neq 0$ while the following argument can
be generalized to other cases. Then, $z^{m}\det\left[h(z)\right]$ is an
analytic function of $z$, and $w_{\mu}$ is related to the number of zeros of
$z^{m}\det\left[h(z)\right]$ within the circle $|z|=1$,

$\displaystyle w_{\mu}$ $\displaystyle={\rm
i}\oint_{|z|=1}\frac{dz}{2\pi}\partial_{z}\log\left[z^{-m}z^{m}\det\left[h(z)\right]\right]$
(53) $\displaystyle=m+{\rm
i}\oint_{|z|=1}\frac{dz}{2\pi}\partial_{z}\log\left[z^{m}\det\left[h(z)\right]\right]$
$\displaystyle=m+{\rm
i}\oint_{|z|=1}\frac{dz}{2\pi}\frac{\partial_{z}\left[z^{m}\det\left[h(z)\right]\right]}{z^{m}\det\left[h(z)\right]}$
$\displaystyle=m-Z,$

where $Z$ is the weighted number of the zeros, and the residue theorem is used
in the last equality. $z=0$ should not be a zero of
$z^{m}\det\left[h(z)\right]$, since the lowest order of
$\det\left[h(z)\right]$ is $z^{-m}$. Thus, $Z$ is equal to the weighted number
of the zeros of $\det\left[h(z)\right]$ in the disk $|z|<1$.

The number of the zeros of $\det\left[h(z)\right]$ is determined by the LEs of
$h$ along the $\mu$ direction. For $\det\left[h(z)\right]=0$, $h(z)$ has a
zero mode. The presence of the zero modes is given by the transfer matrices of
$h(z)$ along the $\mu$ direction. The transfer matrices for each layer are
given by,

$M_{i}(z)=\begin{cases}\begin{pmatrix}-v_{-}^{-1}\tilde{h}_{1}&-\frac{1}{z}v_{-}^{-1}v_{+}\\\
1&0\\\ \end{pmatrix}&(i=1),\\\
\begin{pmatrix}-v_{-}^{-1}\tilde{h}_{i}&-v_{-}^{-1}v_{+}\\\ 1&0\\\
\end{pmatrix}=M_{i}&(i=2,3,\cdots,L_{\mu}-1),\\\
\begin{pmatrix}-\frac{1}{z}v_{-}^{-1}\tilde{h}_{L_{\mu}}&-\frac{1}{z}v_{-}^{-1}v_{+}\\\
1&0\\\ \end{pmatrix}&(i=L_{\mu}).\\\ \end{cases}$ (54)

Here, $M_{1}(z)$ and $M_{L_{\mu}}(z)$ satisfy

$M_{1}(z)M_{L_{\mu}}(z)=\begin{pmatrix}\frac{1}{z}v_{-}^{-1}\tilde{h}_{1}v_{-}^{-1}\tilde{h}_{L_{\mu}}-\frac{1}{z}v_{-}^{-1}v_{+}&\frac{1}{z}v_{-}^{-1}H_{1}v_{-}^{-1}v_{+}\\\
-\frac{1}{z}v_{-}^{-1}\tilde{h}_{L_{\mu}}&-\frac{1}{z}v_{-}^{-1}v_{+}\\\
\end{pmatrix}=\frac{1}{z}M_{1}(z=1)M_{L_{\mu}}(z=1)\,.$ (55)

Suppose that $(A_{1},A_{2},\cdots,A_{L_{\mu}})^{T}$ is a zero mode of $h(z)$
under the periodic boundary conditions. Then, we have

$\begin{pmatrix}A_{2}\\\ A_{1}\end{pmatrix}=M_{1}(z)\begin{pmatrix}A_{1}\\\
A_{L_{\mu}}\end{pmatrix},\quad\begin{pmatrix}A_{3}\\\
A_{2}\end{pmatrix}=M_{2}\begin{pmatrix}A_{2}\\\
A_{1}\end{pmatrix},\quad\cdots,\quad\begin{pmatrix}A_{1}\\\
A_{L_{\mu}}\end{pmatrix}=M_{L_{\mu}}(z)\begin{pmatrix}A_{L_{\mu}}\\\
A_{L_{\mu}-1}\end{pmatrix},$ (56)

and hence

$\displaystyle\begin{pmatrix}A_{L_{\mu}}\\\ A_{L_{\mu}-1}\end{pmatrix}$
$\displaystyle=M_{L_{\mu}-1}\cdots
M_{2}M_{1}(z)M_{L_{\mu}}(z)\begin{pmatrix}A_{L_{\mu}}\\\
A_{L_{\mu}-1}\end{pmatrix}\equiv\frac{1}{z}M\begin{pmatrix}A_{L_{\mu}}\\\
A_{L_{\mu}-1}\end{pmatrix},$ (57)

where $M\equiv M_{L_{\mu}-1}\cdots M_{2}M_{1}(z=1)M_{L_{\mu}}(z=1)$ is the
product of the transfer matrices without the magnetic flux ($\phi=0$, $z=1$).
Since $M_{1}(z=1)$ and $M_{L_{\mu}}(z=1)$ are statistically equivalent to
$M_{L_{\mu}-1}$, $M_{L_{\mu}-2}$, $\cdots$, $M_{3}$, and $M_{2}$, the
eigenvalues of $\frac{1}{L_{\mu}}\ln M$ in the limit
$L_{\mu}\rightarrow\infty$ are characterized by the LEs of $h$.

If $h(z)$ has a zero mode for a complex value $z$, $M$ has an eigenvalue of
$z$ from Eq. (57), and vice versa. Thus, the number $Z$ of the zeros of
$\det\left[h(z)\right]$ within the disk $|z|<1$ is equivalent to the number of
eigenvalues of $M$ whose absolute values are smaller than $1$. The product $M$
of $L_{\mu}$ random matrices has eigenvalues $e^{\alpha_{j}+{\rm
i}\beta_{j}}(j=1,2,\cdots,2m;\alpha_{j},\beta_{j}\in\mathbb{R})$, and
generally, $\alpha_{j}$ grows linearly with $L_{\mu}$, satisfying

$\gamma_{j}=\lim_{L_{\mu}\rightarrow\infty}\frac{\alpha_{j}}{L_{\mu}}\,,$ (58)

with the LE $\gamma_{j}$ of $M$ [66, 79, 81]. Thus, $Z$ is also the same as
the number $N_{-}$ of the negative LEs of $h$:

$Z=N_{-}.\,$ (59)

In terms of Eq. (53), the 1D winding number is given by

$\displaystyle w_{\mu}$ $\displaystyle=m-N_{-}=\frac{1}{2}(N_{+}-N_{-}),$ (60)

where $N_{+}$ and $N_{-}$ are the numbers of positive and negative LEs,
satisfying $N_{+}+N_{-}=2m$. The weak topological index along the $\mu$
direction is the 1D winding number normalized by the degrees of freedom of $h$
in each layer,

$\displaystyle\nu_{\mu}=\frac{1}{m}w_{\mu}=\frac{1}{2m}(N_{+}-N_{-}).$ (61)

Notably, if a LE is exactly zero, the localization length $\xi_{\mu}$ along
the $\mu$ direction diverges and the winding number $w_{\mu}$ is ill defined.
In the quasi-1D geometry of a 3D disordered Hamiltonian $(L\times L\times
L_{\mu},L_{\mu}\gg L)$, $L^{2}$ LEs are distributed within a finite range and
form a continuous spectrum in the thermodynamic limit $L\rightarrow\infty$.
When the spectrum crosses zero with changing $W$, $w_{\mu}$ also changes from
an integer to another integer, $\nu_{\mu}$ continuously changes with $W$, and
the localization length always diverges.

####  I.5.1 Winding number in the clean limit

When a $d$-dimensional system has translation invariance in all the
$d$-dimensional coordinates, Eq. (51) reduces to

$\displaystyle w_{\mu}={\rm i}\int_{0}^{2\pi}\frac{dk_{\mu}}{2\pi}\\!\
\partial_{k_{\mu}}{\rm Tr}\big{[}\log\left[h({\bm{k}})\right]\big{]}\,,$ (62)

with the momentum ${\bm{k}}\equiv(k_{1},k_{2},\cdots,k_{d-1},k_{d})$ and the
Bloch Hamiltonian $h({\bm{k}})$. The trace includes the sum over momenta along
the directions complementary to the $\mu$ direction. For example, the Bloch
Hamiltonian for the 3D nodal-line semimetal model is given by the two-by-two
matrix

${\cal H}(\bm{k})=[\Delta+2t_{\perp}\left(\cos k_{x}+\cos
k_{y}\right)+2t_{\|}^{\prime}\cos k_{z}]\sigma_{z}-2t_{\|}\sin
k_{z}\sigma_{y}.$ (63)

In the canonical basis where the chiral operator is diagonal, the matrix takes
the block off-diagonal structure with

$h({\bm{k}})=2{\rm i}\\!\ t_{\|}\sin k_{z}-[\Delta+2t_{\perp}\left(\cos
k_{x}+\cos k_{y}\right)+2t_{\|}^{\prime}\cos k_{z}].$ (64)

The complex number $h({\bm{k}})$ winds around zero when $k_{z}$ changes from
$0$ to $2\pi$. For $|\Delta+2t_{\perp}\left(\cos k_{x}+\cos
k_{y}\right)|>2|t_{\|}^{\prime}|$, ${\cal H}$ has an energy gap around $E=0$
and $h({\bm{k}})$ winds around zero clockwise for all $k_{x}$ and $k_{y}$ in
the first Brillouin zone, leading to $w_{z}=L^{2}$ and $\nu_{z}=1$. Here,
$L^{2}$ is the system size within the $xy$ plane. For
$|\Delta+2t_{\perp}\left(\cos k_{x}+\cos k_{y}\right)|<2|t_{\|}^{\prime}|$,
${\cal H}$ has a gap at $E=0$, but $h({\bm{k}})$ does not wind around zero for
any $k_{x}$ and $k_{y}$, leading to $w_{z}=\nu_{z}=0$. For
$2|t_{\|}^{\prime}|-4t_{\perp}<\Delta<2|t_{\|}^{\prime}|+4t_{\perp}$, zero
modes of ${\cal H}$ form a nodal ring in momentum space, and the winding
number is $+1$ ($0$) for the wave numbers $k_{x},k_{y}$ inside (outside) the
nodal ring, leading to $0<\nu_{z}<1$.

###  I.6 Statistical symmetry

An ensemble of disordered Hamiltonians, as a whole, can be invariant under a
symmetry operation even if each disorder realization in the ensemble breaks
the symmetry. Such symmetry of the ensemble is dubbed statistical symmetry
[82]. Statistical symmetry does not influence the symmetry class of the
Hamiltonians since it is not a symmetry of each disordered Hamiltonian. An
example of statistical symmetry is translation symmetry in 3D weak topological
insulators with disorder [83].

Statistical symmetry can make the ensemble averages of physical observables or
topological indices be zero. A prime example is zero Hall conductance due to
statistical time-reversal symmetry. Suppose that ${\cal H}_{\alpha}$ is a
time-reversal-breaking Hamiltonian in an ensemble with statistical time-
reversal symmetry and has a finite Hall conductance $\sigma_{xy}^{\alpha}$. A
time-reversed counterpart ${\cal H}_{\bar{\alpha}}$ of ${\cal H}_{\alpha}$
exists in the same ensemble and has the opposite value $-\sigma_{xy}^{\alpha}$
of the Hall conductance. Such an ensemble has zero Hall conductance on
average,

$\langle\sigma_{xy}\rangle=\frac{1}{2N_{\rm
sample}}\sum_{\alpha}\left(\sigma_{xy}^{\alpha}-\sigma_{xy}^{\alpha}\right)=0\,.$
(65)

####  I.6.1 Statistical symmetry, Lyapunov exponents, and one-dimensional
winding number

In Sec. I.5, the 1D winding numbers and weak topological indices are defined
in terms of the right-upper part $h$ of the chiral-symmetric Hamiltonian
${\cal H}$ in the canonical basis [see Eq. (29)]. Now, we introduce
statistical symmetry of $h$ and show that it requires LEs of $h$ to come in
opposite-sign pairs. Statistical symmetry also makes the 1D winding numbers
and weak topological indices be zero as a whole. Suppose that an ensemble of
$h$ with different disorder realizations, $\big{\\{}\\!\ h\\!\ \big{|}\\!\
\epsilon_{\bm{r}}\in[-W/2,W/2]\big{\\}}$, is symmetric under transposition of
$h$ together with a certain unitary transformation ${\cal U}$:

$\displaystyle\Big{\\{}\\!\ h\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}}=\Big{\\{}\\!\
h^{\prime}\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}},\quad{\rm with}\quad
h^{\prime}\equiv{\cal U}h^{T}{\cal U}^{\dagger}.$ (66)

Here, we assume that the unitary transformation ${\cal U}$ is diagonal in a
spatial coordinate $r_{\mu}$ and is independent of $r_{\mu}$ while it can be
non-diagonal in the other coordinates ${\bm{s}}$,

$\displaystyle({\cal
U})_{{\bm{r}},{\bm{r}}^{\prime}}=\delta_{r_{\mu},r^{\prime}_{\mu}}(u)_{{\bm{s}},{\bm{s}}^{\prime}},$
(67)

with ${\bm{r}}\equiv({\bm{s}},r_{\mu})$ and
${\bm{r}}^{\prime}\equiv({\bm{s}}^{\prime},r^{\prime}_{\mu})$. Then, LEs of
$h$ along the $\mu$ direction come in opposite-sign pairs. The 1D winding
number and weak topological index of $h$ along the $\mu$ direction vanish from
Eqs. (60) and (61).

To see this, we decompose $h$ into a layer structure along the $\mu$
direction,

$(h)_{i,j}=\tilde{h}_{i}\delta_{i,j}+v_{+}\delta_{i,j-1}+v_{-}\delta_{i,j+1}\,,$
(68)

with $i,j=1,\cdots,L_{\mu}$. $\tilde{h}_{i}$ is a block of matrix elements of
$h$ within the $i$th layer, and $v_{\pm}$ is a block of matrix elements of $h$
between the $i$th layer and $(i\mp 1)$th layer. $v_{\pm}$ are free from
disorder and independent of the layer index. If the degree of freedom in each
layer of $h$ is $m$, $\tilde{h}_{i}$ and $v_{\pm}$ are $m\times m$ matrices.
Similarly, $h^{\prime}\equiv{\cal U}h^{T}{\cal U}^{\dagger}$ is also
decomposed into a layer structure along the $\mu$ direction,

$(h^{\prime})_{i,j}=u\tilde{h}_{i}^{T}u^{\dagger}\delta_{i,j}+uv^{T}_{-}u^{\dagger}\delta_{i,j-1}+uv^{T}_{+}u^{\dagger}\delta_{i,j+1}\,.$
(69)

From Eqs. (24) and (25), the transfer matrices of $h$ and $h^{\prime}$ are
obtained as

$\displaystyle
M_{i}=\left(\begin{array}[]{cc}-(v_{-})^{-1}\tilde{h}_{i}&-(v_{-})^{-1}v_{+}\\\
1&0\\\ \end{array}\right),\ \ \
M^{\prime}_{i}=\left(\begin{array}[]{cc}-u(v^{T}_{+})^{-1}\tilde{h}^{T}_{i}u^{\dagger}&-u(v^{T}_{+})^{-1}v^{T}_{-}u^{\dagger}\\\
1&0\\\ \end{array}\right).$ (74)

The two matrices are related to each other by the following symmetry,

$\displaystyle S{M^{\prime}_{i}}^{T}S^{-1}=M^{-1}_{i},$ (75)

with

$\displaystyle S\equiv\left(\begin{array}[]{cc}0&-(v_{-})^{-1}u^{T}\\\
(v_{+})^{-1}u^{T}&0\\\ \end{array}\right).$ (78)

Since $u$ and $v_{\pm}$ in $S$ are independent of the layer index, the same
symmetry holds between the products of the transfer matrices, $M\equiv
M_{L_{\mu}}\cdots M_{1}$ and $M^{\prime}\equiv M^{\prime}_{L_{\mu}}\cdots
M^{\prime}_{1}$,

$\displaystyle S{M^{\prime}}^{T}S^{-1}=M^{-1}.$ (79)

Since a non-singular similarity transformation $S$ does not change LEs,
eigenvalues of
$P^{\prime}(0)\equiv\lim_{L_{\mu\rightarrow\infty}}\frac{1}{2L_{\mu}}\ln({M^{\prime}}^{\dagger}M^{\prime})$
are opposite to eigenvalues of
$P(0)\equiv\lim_{L_{\mu\rightarrow\infty}}\frac{1}{2L_{\mu}}\ln(M^{\dagger}M)$.
Since $h$ and $h^{\prime}$ are in the same ensemble, $M_{i}$ and
$M^{\prime}_{j}$ ($i,j=1,\cdots,L_{\mu}$) are random matrices with the same
possibility distribution. Thus, according to Oseledec’s theorem [79], the
eigenvalues of $P(0)$ and $P^{\prime}(0)$ converge to the same values in the
limit $L_{\mu}\rightarrow\infty$. Then, the eigenvalues of $P(0)$, as well as
the LEs of $h$, must come in opposite-sign pairs in the limit
$L_{\mu}\rightarrow\infty$. Because of $N_{+}=N_{-}$ for the LEs along the
$\mu$ direction, $\nu_{\mu}$ and $w_{\mu}$ of $h$ vanish from Eqs. (60) and
(61).

As discussed in Sec. I.4, the LEs of $h$ and the LEs of $h^{\dagger}$
generally come in opposite-sign pairs. Thus, statistical Hermitian-conjugation
symmetry of $h$ also requires the LEs of $h$ to come in the opposite-sign
pairs. Suppose that an ensemble of $h$ with different disorder realizations,
$\big{\\{}\\!\ h\\!\ \big{|}\\!\ \epsilon_{\bm{r}}\in[-W/2,W/2]\big{\\}}$, is
symmetric under Hermitian conjugation together with a certain unitary
transformation ${\cal U}$ defined in Eq. (67):

$\displaystyle\Big{\\{}\\!\ h\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}}=\Big{\\{}\\!\
h^{\prime}\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}},\quad{\rm with}\quad
h^{\prime}\equiv{\cal U}h^{\dagger}{\cal U}^{\dagger}.$ (80)

Then, the LEs of $h$ along the $\mu$ direction should come in opposite-sign
pairs, leading to $\nu_{\mu}=w_{\mu}=0$.

It is also notable that statistical symmetry of $h$ leads to statistical
symmetry of ${\cal H}$. If an ensemble of $h$ is invariant under the
combination of transposition and a unitary operation in Eq. (66), the
corresponding ensemble of ${\cal H}$ is invariant under the combination of
time reversal and a unitary operation:

$\Big{\\{}\\!\ {\cal H}=\begin{pmatrix}0&h\\\ h^{\dagger}&0\end{pmatrix}\\!\
\Big{|}\\!\ \epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}}=\Big{\\{}\\!\
{\cal H}^{\prime}=\begin{pmatrix}0&h^{\prime}\\\
h^{\prime\dagger}&0\end{pmatrix}\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}},\quad{\cal
H}^{\prime}\equiv\begin{pmatrix}0&{\cal U}\\\ {\cal U}&0\end{pmatrix}{\cal
H}^{*}\begin{pmatrix}0&{\cal U}^{\dagger}\\\ {\cal
U}^{\dagger}&0\end{pmatrix}.$ (81)

If an ensemble of $h$ is invariant under the combination of Hermitian
conjugation and a unitary transformation in Eq. (80), the corresponding
ensemble of ${\cal H}$ is invariant under the following unitary operation:

$\Big{\\{}\\!\ {\cal H}=\begin{pmatrix}0&h\\\ h^{\dagger}&0\end{pmatrix}\\!\
\Big{|}\\!\ \epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}}=\Big{\\{}\\!\
{\cal H}^{\prime}=\begin{pmatrix}0&h^{\prime}\\\
h^{\prime\dagger}&0\end{pmatrix}\\!\ \Big{|}\\!\
\epsilon_{\bm{r}}\in\big{[}-W/2,W/2\big{]}\Big{\\}},\quad{\cal
H}^{\prime}\equiv\begin{pmatrix}0&{\cal U}\\\ {\cal U}&0\end{pmatrix}{\cal
H}\begin{pmatrix}0&{\cal U}^{\dagger}\\\ {\cal U}^{\dagger}&0\end{pmatrix}.$
(82)

####  I.6.2 Statistical symmetry in the nodal-line semimetal model

We show that the 3D nodal-line semimetal model has statistical transposition
symmetries for the $\mu=x,y$ directions. The nodal-line semimetal model takes
a block off-diagonal structure in the canonical basis, where the upper-right
part $h$ is given by Eq. (45). Transposition of $h$ exchanges
$t^{\prime}_{\|}+t_{\|}$ and $t^{\prime}_{\|}-t_{\|}$ in Eq. (45). Since the
disorder potential $\epsilon_{\bm{r}}$ is statistically equivalent for
different lattice points ${\bm{r}}$, we can introduce a mirror operation with
respect to the $xy$ plane as a unitary transformation of Eqs. (66) and (67),

$\displaystyle{\cal
U}_{(r_{x},r_{y},r_{z}|r^{\prime}_{x},r^{\prime}_{y},r^{\prime}_{z})}=\delta_{r_{x},r^{\prime}_{x}}\delta_{r_{y},r^{\prime}_{y}}\delta_{r_{z},-r^{\prime}_{z}}.$
(83)

This mirror operation exchanges $t^{\prime}_{\|}+t_{\|}$ and
$t^{\prime}_{\|}-t_{\|}$ as well as $\epsilon_{r_{x},r_{y},r_{z}}$ and
$\epsilon_{r_{x},r_{y},-r_{z}}$ while $\epsilon_{{\bm{r}}}$ is statistically
equivalent for different ${\bm{r}}$. Thus, an ensemble for $h$ defined by Eq.
(45) is statistically invariant under transposition with the unitary
transformation. The symmetry of Eqs. (66) and (83) requires the LEs of $h$
along the $x$ ($y$) direction to come in opposite-sign pairs, leading to
$w_{x(y)}=\nu_{x(y)}=0$.

###  I.7 Finite-size scaling of Lyapunov exponents

The LEs of a chiral-symmetric Hamiltonian ${\cal H}$ are the sum of LEs of the
right-upper part $h$ of ${\cal H}$ in the canonical basis and their opposite-
sign exponents [see Eq. (29)]. In the quasi-one-dimensional (quasi-1D)
geometry, the LEs of ${\cal H}$, as well as the LEs of $h$, comprise continuum
spectra for the limit $L\to\infty$ [62]. In the nodal-line semimetal model
${\cal H}$ with $t_{\|}=t_{\|}^{\prime}$, $m=L^{2}$ LEs of diverge to
$+\infty$, and the other $m$ LEs of $h$ form a finite spectrum around
$\gamma=0$ (see also Fig. 2 in the main text). For $t_{\|}\neq
t_{\|}^{\prime}$, on the other hand, all the $2m=2L^{2}$ LEs of $h$ comprise
either one or two continuum spectra around $\gamma=0$, depending on the
disorder strength (see Figs. 9(c) and 9(d)).

In the following discussion, we focus on the case with the two continuum LE
spectra. Generalization to the other cases is straightforward. For finite $L$,
the transfer matrix study of the non-Hermitian Hamiltonian $h$ in the quasi-1D
geometry gives a discrete set of $2m$ LEs,

$\displaystyle\big{\\{}\gamma^{(1)}_{\rm min}(W,L),\cdots,\gamma^{(1)}_{\rm
max}(W,L),\gamma^{(2)}_{\rm min}(W,L),\cdots,\gamma^{(2)}_{\rm
max}(W,L)\big{\\}},$ (84)

with $\gamma^{(1)}_{\rm min}(W,L)<\cdots<\gamma^{(1)}_{\rm
max}(W,L)<\gamma^{(2)}_{\rm min}(W,L)<\cdots<\gamma^{(2)}_{\rm max}(W,L)$. In
the limit $L\to\infty$, all the $m$ LEs from $\gamma^{(j)}_{\rm min}(W,L)$ to
$\gamma^{(j)}_{\rm max}(W,L)$ ($j=1,2$) form a continuum spectrum that ranges
from $\gamma^{(j)}_{\rm min}(W)$ to $\gamma^{(j)}_{\rm max}(W)$, satisfying

$\displaystyle\lim_{L\rightarrow\infty}\gamma^{(j)}_{\rm
min}(W,L)\equiv\gamma^{(j)}_{\rm min}(W),\ \
\lim_{L\rightarrow\infty}\gamma^{(j)}_{\rm max}(W,L)\equiv\gamma^{(j)}_{\rm
max}(W),$ (85)

with $j=1,2$. For some disorder strength, a finite gap
$2\Delta\equiv\gamma^{(2)}_{\rm min}(W)-\gamma^{(1)}_{\rm max}(W)$ exists
between the two continuum LEs spectra (see Figs. 7, 9(c), and 9(d)).

When the gap $2\Delta$ is much larger than $L^{-1}$, $\gamma^{(1)}_{\rm
max}(W,L)$ and $\gamma^{(2)}_{\rm min}(W,L)$ can be fitted well by the
following scaling functions:

$\displaystyle\gamma^{(1)}_{\rm max}(W,L)=-\frac{a}{L}+\gamma^{(1)}_{\rm
max}(W).$ (86)

and

$\displaystyle\gamma^{(2)}_{\rm min}(W,L)=\frac{a}{L}+\gamma^{(2)}_{\rm
min}(W),$ (87)

Notably, this scaling holds irrespective of whether their limits
$\gamma^{(2)}_{\rm min}(W)$ and $\gamma^{(1)}_{\rm max}(W)$ are close to zero.
To see this, we first note that LEs of $h$ with different $t_{\|}$ and
$t^{\prime}_{\|}$ are related by an imaginary gauge transformation along the
$z$ direction [68]. Let the imaginary gauge transformation with an imaginary
gauge ${\rm i}g$ act on $h$ by

$h\rightarrow h_{g}\equiv V_{g}hV_{g}^{-1},$ (88)

where $V_{g}$ is a diagonal matrix whose diagonal element takes $e^{jg}$ in
the $j$th layer along the $z$ direction:

$V_{g}=\begin{pmatrix}e^{g}1_{m\times m}&0&0&\cdots&0\\\ 0&e^{2g}1_{m\times
m}&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&\cdots&\cdots&\cdots&e^{L_{z}g}1_{m\times m}\\\ \end{pmatrix}\,.$ (89)

The transfer matrix $M_{i}(g)$ of $h_{g}\equiv V_{g}hV_{g}^{-1}$ along the $z$
direction is obtained from Eq. (46) as

$M_{i}(g)=\left(\begin{matrix}-e^{g}\frac{1}{t^{\prime}_{\|}-t_{\|}}\tilde{h}_{i}&-e^{2g}\frac{t^{\prime}_{\|}+t_{\|}}{t^{\prime}_{\|}-t_{\|}}1_{m\times
m}\\\ 1_{m\times m}&0_{m\times m}\\\ \end{matrix}\right)\,,$ (90)

and satisfies

$M_{i}(g)=e^{g}S\left(\begin{matrix}-\frac{1}{t^{\prime}_{\|}-t_{\|}}\tilde{h}_{i}&-\frac{t^{\prime}_{\|}+t_{\|}}{t^{\prime}_{\|}-t_{\|}}1_{m\times
m}\\\ 1_{m\times m}&0_{m\times m}\\\
\end{matrix}\right)S^{-1}=e^{g}SM_{i}S^{-1},\quad
S=\left(\begin{matrix}1_{m\times m}&0_{m\times m}\\\ 0_{m\times
m}&e^{-g}1_{m\times m}\\\ \end{matrix}\right)\,.$ (91)

Thus, the LEs of $h_{g}$ are obtained from the LEs of $h$ in Eq. (84),

$\big{\\{}\gamma^{(1)}_{\rm min}(W,L)+g,\cdots,\gamma^{(1)}_{\rm
max}(W,L)+g,\gamma^{(2)}_{\rm min}(W,L)+g,\cdots,\gamma^{(2)}_{\rm
max}(W,L)+g\big{\\}}.$ (92)

Similarly, the LEs of $h^{\dagger}_{g}$ differ from the LEs of $h^{\dagger}$
by $-g$. Suppose that the gap between $\gamma^{(1)}_{\rm max}(W)$ and
$\gamma^{(2)}_{\rm min}(W)$ is much larger than $L^{-1}$, and choose the
imaginary gauge $g$ in such a way that a midpoint of the gap comes around
zero,

$\displaystyle\gamma^{(1)}_{\rm max}(W)+g<0<\gamma^{(2)}_{\rm min}(W)+g.\ \ $
(93)

Then, a zero mode of a chiral-symmetric Hamiltonian ${\cal H}_{g}$ that has
$h_{g}$ and $h^{\dagger}_{g}$ in the off-diagonal blocks is in the Anderson
insulator phase, and its localization length $\xi_{z}$ is much shorter than
$L$. Depending on $g$, the localization length is given by either

$\displaystyle\frac{1}{\xi_{z}(W,L)}=-(\gamma^{(1)}_{\rm max}(W,L)+g),$ (94)

or

$\displaystyle\frac{1}{\xi_{z}(W,L)}=\gamma^{(2)}_{\rm min}(W,L)+g.$ (95)

In the Anderson insulator phase, a finite-size scaling of the normalized
localization length $\Lambda(L)\equiv\xi(L)/L$ is described by a function of
the single parameter $\Lambda$ (i.e., single-parameter scaling) [69]:

$\frac{d\ln\Lambda}{d\ln L}=\beta({\Lambda}).$ (96)

For small $\Lambda$, we have $\beta(\Lambda)\to-1$ and
$\ln\Lambda\simeq\ln\xi(L=\infty)-\ln L$. When the localization length $\xi$
is much shorter than $L$, one may expand the $\beta$ function in small
$\Lambda$, $\beta(\Lambda)=-1+a\Lambda+{\cal O}(\Lambda^{2})$, and retain the
zeroth and first order in $\Lambda$,

$\frac{d\ln\Lambda}{d\ln L}=-1+a\Lambda.$ (97)

This differential equation may be solved by an integration in a domain of
$[L,L_{0}]$ with $L\ll L_{0}$,

$\displaystyle\frac{L}{\xi}\equiv\frac{1}{\Lambda}=a+\frac{L}{L_{0}}\Big{(}\frac{1}{\Lambda_{0}}-a\Big{)}\equiv
a+\frac{L}{L_{0}}\Big{(}\frac{L_{0}}{\xi_{0}}-a\Big{)},$ (98)

with $\xi_{0}\equiv\xi(L=L_{0})$. When $L_{0}$ goes to infinity, $\xi_{0}$
converges to finite $\xi(L=\infty)\equiv\xi(\infty)$. Thus, we obtain the
lowest-order finite-size scaling form of the quasi-1D localization length,

$\displaystyle\frac{1}{\xi(L)}=\frac{a}{L}+\frac{1}{\xi(\infty)}.$ (99)

Now that the localization length along the $z$ direction is much shorter than
$L$ in Eqs. (94) and (95), we may use the same scaling function not only for
$\xi_{z}(W,L)$ but also for $\gamma^{(2)}_{\rm min}(W,L)$ and
$\gamma^{(1)}_{\rm max}(W,L)$. In fact, the scaling forms of Eqs. (86) and
(87) work well for the numerical fittings. Note that the coefficient $a$ in
Eq. (99) takes a non-universal value in general (see the fitting values in
Table 8). To obtain the scaling form for $\gamma^{(1)}_{\rm min}(W,L)$, let us
choose large positive $g$ and make all the LEs of $h_{g}$ be positive,

$\displaystyle 0<\gamma^{(1)}_{\rm min}(W,L)+g.$ (100)

In the chiral-symmetric Hamiltonian ${\cal H}_{g}$ that has such $h_{g}$ and
its Hermitian conjugate $h_{g}^{{\dagger}}$, $E=0$ is in the weak topological
insulator phase $(\nu_{z}=1)$ and the localization length $\xi_{z}$ is given
by $\gamma^{(1)}_{\rm min}(W,L)+g$. If we assume that the finite-size scaling
of $\xi_{z}$ in the weak topological insulator phase is also described by the
single parameter scaling function of $\Lambda_{z}\equiv\xi_{z}/L$, we also
obtain the scaling function for $\gamma^{(1)}_{\rm min}(W,L)$ as

$\displaystyle\gamma^{(1)}_{\rm min}(W,L)=\frac{a}{L}+\gamma^{(1)}_{\rm
min}(W).$ (101)

This scaling form also works well for the numerical data of $\gamma^{(1)}_{\rm
min}(W,L)$.

Figure 4: Schematic pictures of Lyapunov exponents (LEs) of the nodal-line
semimetal model ${\cal H}$ and LEs of ${\cal H}_{g}$ generated by an imaginary
gauge transformation. The LEs of ${\cal H}$ and ${\cal H}_{g}$ are the sum of
LEs of its right-upper non-Hermitian Hamiltonian $h$ and their Hermitian
conjugate $h^{{\dagger}}$ in Eq. (29). The LEs calculated in the quasi-1D
geometry ($L\times L\times L_{z}$, $L_{z}\gg L$) comprise two continuum LEs
spectra in the limit $L\rightarrow\infty$. The grey-shaded regions in the left
(right) sides of the vertical axis denote the continuous spectra formed by LEs
of $h$ ($h^{\dagger}$). The dotted horizontal line denotes zero $\gamma=0$.
The LEs of ${\cal H}$ and ${\cal H}_{g}$, as a whole, are symmetric around
zero. The smallest positive or the largest negative LE (marked in bold)
corresponds to the inverse of the quasi-1D localization length. In (a), $g$ is
chosen such that finite $\gamma^{(2)}_{\rm min}(W)+g$ corresponds to the
inverse of the localization length of ${\cal H}_{g}$. In (b), $g$ is chosen
such that $\gamma^{(1)}_{\rm min}(W)+g$ corresponds to the inverse of the
localization length.

####  I.7.1 Numerical fitting

To show the validity of the scaling forms in Eqs. (86) and (87), we use the
standard $\chi^{2}$ fitting method to fit the data of $\gamma^{(i)}_{\rm
min/max}(W,L)$ ($i=1,2$) with larger $L$. For fixed $W$, we minimize the
following $\chi^{2}$ function with respect to $a$ and $\gamma^{(i)}_{\rm
min/max}(W)$ in Eqs. (86) and (87),

$\chi^{2}=\sum^{D}_{j=1}\left(\frac{F_{j}-\gamma_{\rm
min/max}^{(i)}(W,L)}{\sigma_{j}}\right)^{2},$ (102)

where $j$ specifies the data point of $\gamma^{(i)}_{\rm min/max}(W,L)$ with
different $L$. The number $D$ of the data points is typically $7$ ($14\leq
L\leq 28$) and $9$ ($18\leq L\leq 34$). $F_{j}$ is the fitted value from Eqs.
(86) and (87) for different $L$ specified by $j$. $\sigma_{j}$ is the standard
deviation of $\gamma_{\rm min/max}^{(i)}(W,L)$ estimated from the transfer
matrix calculation [61]. The finite-size scaling fit works well in the nodal-
line semimetal models with or without time-reversal symmetry that are studied
in this work. Some fitting results for the nodal-line semimetal models
($\Delta=0$, $t_{\|}=t_{\|}^{\prime}=1/2$, $t_{\perp}=1$) that are studied in
the main text are shown in Fig. 5 (a,c) and Table 4. In Table 4, the Monte
Carlo method is used to generate pseudo-data sets and evaluate the 95%
confidence interval of the fitted values of $``a"$ and $``\gamma_{\rm
min/max}^{(1)}(W)"$.

Using $\gamma^{(1)}_{\rm max}(W)$ with the confidence interval, we determine
the critical disorder strength $W_{c}^{(z)}$ between the non-localized region
and the localized phase. For example, the fitting results in Table 4 show
$\gamma_{\text{max}}=0.0020\\!\ \\!\ [0.0017,0.0022]$ for $W=29.5$ and
$\gamma_{\text{max}}=-0.0042\\!\ \\!\ [-0.0044,0.0040]$ for $W=29.4$ in the
nodal-line semimetal model in symmetry class BDI. The results suggest that
$W_{c}^{(z)}$ must be between $29.4$ and $29.5$ with the 95% confidence:
$W_{c}^{(z)}=29.45\\!\ \\!\ [29.4,29.5]$. $W^{(z)}_{c}$ for the nodal-line
semimetal models in symmetry classes BDI and AIII are summarized in Table 5.
Table 5 also shows $W^{(x)}_{c}$ in the same nodal-line semimetal models for
the same sets of parameters. Comparisons between $W_{c}^{(z)}$ and
$W_{c}^{(x)}$ illustrate that $W^{(x)}_{c}<W^{(z)}_{c}$ and
$|W^{(z)}_{c}-W^{(x)}_{c}|$ is around 10% of $W^{(x)}_{c}$. This concludes the
presence of the quasi-localized phase in these nodal-line semimetal models.
Similarly, comparisons of $W_{c}^{(z)}$ and $W_{c}^{(x)}$ in other types of
nodal-line semimetal models (Tables 7 and 9) also suggest the presence of the
quasi-localized phases (see below).

Table 4: Finite-size scaling analysis of $\gamma^{(1)}_{\rm max}(W,L)$ for several disorder strength $W$ around $W=W_{c}^{(z)}$ in the nodal-line semimetal models ($\Delta=0$, $t_{\|}=t_{\|}^{\prime}=1/2$, $t_{\perp}=1$) with time-reversal symmetry (symmetry class BDI) and without time-reversal symmetry (symmetry class AIII). The square brackets are the 95% confidence intervals determined by the Monte Carlo analyses. symmetry class | $W$ | $L$ | $\gamma_{\rm max}(W,L=\infty)$ | $a$ | GOF
---|---|---|---|---|---
BDI |  29.4 |  18 - 34 |  0.0020 [0.0017,0.0022] |  -3.890 [-3.893,-3.888] |  0.58
BDI | 29.5 | 18 - 34 | -0.0042 [-0.0044,-0.0040] | -3.689 [3.691,-3.687] | 0.14
AIII | 9.70 | 14 - 28 | 0.0143 [0.0128,0.0156] | -1.921 [-1.948,-1.892] | 0.02
AIII | 9.80 | 14 - 28 | 0.0015 [0.0000,0.0029] | -1.929 [-1.957,-1.899] | 0.50
AIII | 9.90 | 14 - 28 | -0.0114 [-0.0129,-0.0099] | -1.931 [-1.961,-1.904] | 0.22
Table 5: Comparison between the critical disorder strength $W_{c}^{(z)}$ in the $z$ direction (weak topological index $\nu_{z}\neq 0$) and the critical disorder strength $W_{c}^{(x)}$ in the $x$ direction (weak topological index $\nu_{x}=0$) for the nodal-line semimetal models in symmetry classes BDI and AIII ($\Delta=0$, $t_{\|}=t_{\|}^{\prime}=1/2$, $t_{\perp}=1$). The square brackets denote the 95% intervals. The confidence intervals of $W_{c}^{(x)}$ are determined by the 95% confidence intervals of $\gamma_{\rm max}(W)$.  symmetry class | $W_{c}^{(x)}$ | $W_{c}^{(z)}$
---|---|---
BDI |  27.24[27.19,27.30] |  29.45[29.4,29.5]
AIII | 9.14[9.12,9.17] | 9.8[9.7,9.9]

(a) class BDI

(b) class BDI

(c) class AIII

(d) class AIII

Figure 5: (a,c) $\gamma^{(1)}_{\text{max}}(W,L)$ as a function of the system
size $L$ for the different disorder strength $W$ in the nodal-line semimetal
models ($\Delta=0$, $t_{\|}=t_{\|}^{\prime}=1/2$, $t_{\perp}=1$) in (a)
symmetry class BDI and (c) symmetry class AIII. The solid lines are the
fitting curves from Eq. (87). A cross-section of the fitting curve at $1/L=0$
determines
$\gamma^{(1)}_{\text{max}}(W)\equiv\lim_{L\rightarrow\infty}\gamma^{(1)}_{\text{max}}(W,L)$.
(b,d) $\gamma^{(1)}_{\rm max}(W)$ as a function of $W$ around $W=W_{c}^{(z)}$
in the nodal-line semimetal models in (b) symmetry class BDI and (d) symmetry
class AIII. Insets of (b,d): distributions of the Lyapunov exponents (LEs) of
the right-upper part $h$ of the nodal-line semimetal model ${\cal H}$ as a
function of $W$ in the larger range of $\gamma$ and $W$. The LEs of the nodal-
line semimetal models are the sum of the LEs of $h$ and their opposite-sign
exponents.

###  I.8 Quasi-localized phases in chiral-symmetric models with weak
topological indices

In Sec. I.7, we describe the finite-size scaling analysis of the LEs along the
$z$ direction in the nodal-line models with the weak topological index
$\nu_{z}\neq 0$. The analysis enables determinations of the phase boundary
$W^{(z)}_{c}$ of the non-localized region. The non-localized region comprises
the metal and quasi-localized phases. In the nodal-line models with
$\nu_{x}=\nu_{y}=0$, the phase transition between the metal and quasi-
localized phases is characterized by the localization properties along the $x$
or $y$ direction. In this section, we discuss the localization properties
along the $x$ direction in the chiral-symmetric models with $\nu_{z}\neq 0$
and $\nu_{x}=0$. We demonstrate the presence of the quasi-localized phases
inside the non-localized region for all the models.

####  I.8.1 Nodal-line semimetal in class AIII

We discuss a nodal-line semimetal in Eq. (35) with the time-reversal-breaking
disorder. The Hamiltonian ${\cal H}_{1}$ has the two types of random
potentials,

${\cal
H}_{1}=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left\\{(\Delta+\epsilon_{\bm{r}})c^{\dagger}_{\bm{r}}\sigma_{y}c_{\bm{r}}+\epsilon_{\bm{r}}^{\prime}c^{\dagger}_{\bm{r}}\sigma_{z}c_{\bm{r}}+\left[\sum_{\mu=x,y}\left(t_{\perp}c^{\dagger}_{\bm{r}+\bm{e_{\mu}}}\sigma_{z}c_{\bm{r}}\right)-{\rm
i}t_{\|}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{y}c_{\bm{r}}+t_{\|}^{\prime}c^{\dagger}_{\bm{r}+\bm{e_{z}}}\sigma_{z}c_{\bm{r}}+\text{H.c.}\right]\right\\},$
(103)

where the random potential $\epsilon_{\bm{r}}$ ($\epsilon_{\bm{r}}^{\prime}$)
respects (breaks) time-reversal symmetry and distributes uniformly in
$\epsilon_{\bm{r}}^{2}+\epsilon_{\bm{r}}^{\prime 2}\leq W^{2}$. The parameters
are chosen to be $\Delta=0$, $t_{\|}=t_{\|}^{\prime}=1/2$, and $t_{\perp}=1$.
The Hamiltonian ${\cal H}_{1}$ only satisfies chiral symmetry ${\cal
H}_{1}=-\sigma_{x}{\cal H}_{1}\sigma_{x}$, and hence belongs to class AIII.

According to Eq. (30), the chiral-symmetric Hamiltonian is decomposed into the
block-off diagonal structure in the basis that diagonalizes the chiral
operator ${\cal C}\equiv\sigma_{x}$. The right-upper part $h_{1}$ of ${\cal
H}_{1}$ in this basis is given by

$h_{1}=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left[(\Delta+\epsilon_{\bm{r}}+{\rm
i}\epsilon_{\bm{r}}^{\prime})f^{\dagger}_{\bm{r}}f_{\bm{r}}+\sum_{\mu=x,y}\left(t_{\perp}f^{\dagger}_{{\bm{r}}+{\bm{e}}_{\mu}}f_{\bm{r}}+\text{H.c.}\right)+(t^{\prime}_{\|}+t_{\|})f^{\dagger}_{{\bm{r}}+{\bm{e}}_{z}}f_{\bm{r}}+(t^{\prime}_{\|}-t_{\|})f^{\dagger}_{\bm{r}}f_{{\bm{r}}+{\bm{e}}_{z}}\right]\,.$
(104)

Transposition exchanges $t^{\prime}_{\|}+t_{\|}$ and $t^{\prime}_{\|}-t_{\|}$.
Thus, as a unitary transformation in Eqs. (66) and (67), we can consider the
mirror operation with respect to the $xy$ plane as in Eq. (83). Since both
$\epsilon_{\bm{r}}$ and $\epsilon_{\bm{r}}^{\prime}$ are statistically
equivalent for different lattice points ${\bm{r}}$, an ensemble of $h_{1}$
defined in Eq. (104) is statistically invariant under the combination of
transposition and the mirror operation. The symmetry in Eqs. (66) and (83)
requires the LEs of $h_{1}$ along the $x$ and $y$ directions to come in
opposite-sign pairs, leading to $\nu_{x}=\nu_{y}=0$.

We calculate the localization length $\xi_{x}$ of $h_{1}$ along the $x$
direction in the quasi-1D geometry ($L_{x}\times L\times L$, $L_{x}\gg L$).
The normalized localization length $\Lambda_{x}\equiv\xi_{x}/L$ shows scale-
invariant behavior around $W_{c}^{(x)}=9.14\pm 0.01$ (Fig. 6). Fitting by the
polynomial expansion of the finite-size scaling functions [see Eqs. (20) and
(21)], we determine the critical disorder strength $W^{(x)}_{c}$ and the
critical exponent (see Table 3). Figure 6 shows the normalized localization
length for different $W$ and $L$ together with the fitting curves.

In Sec. I.7, we use the finite-size scaling of LEs to obtain the critical
disorder strength $W^{(z)}_{c}=9.8\\!\ [9.7,9.9]$. For $W<W_{c}^{(z)}$, the
localization length along the $z$ direction diverges. $W^{(x)}_{c}$ is well
within the non-localized region, $W^{(x)}_{c}<W^{(z)}_{c}=9.8\\!\ [9.7,9.9]$,
demonstrating the presence of the quasi-localized phase in the nodal-line
semimetal model without time-reversal symmetry.

Figure 6:  Normalized localization length $\Lambda_{x}\equiv\xi_{x}/L$ along
the $x$ direction as a function of the disorder strength $W$ in the nodal-line
semimetal model in class AIII [Eq. (35) with $\Delta=0$,
$t_{\|}=t_{\|}^{\prime}=1/2$, $t_{\perp}=1$]. $\xi_{x}$ is calculated in the
quasi-1D geometry ($L\times L\times L_{x}$). The black points are the raw data
with the error bars. The solid lines for different $L$ and the dashed vertical
line $W_{c}^{(x)}\approx 9.14$ with the error bars are the results of the
fitting according to Eqs. (20) and (21) with $(m,n)=(2,3)$. The dashed line
$W_{c}^{(z)}\approx 9.8$ is evaluated by the fitting of the Lyapunov exponents
along the $z$ direction by Eq. (86).

####  I.8.2 Nodal-line semimetals in class BDI

We also study the localization properties along the $x$ direction in other
nodal-line semimetal models in class BDI: i) nodal-line semimetal model in Eq.
(35) with the different parameters $\Delta=0$, $t_{\|}=\sinh
g,t_{\|}^{\prime}=\cosh g$ ($g=0.22,1$), $t_{\perp}=1$ and ii) another nodal-
line semimetal model ${\cal H}_{2}$ with an extended Fermi line running across
the Brillouin zone,

${\cal
H}_{2}=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left\\{(\Delta+\epsilon_{\bm{r}})c^{\dagger}_{\bm{r}}\sigma_{z}c_{\bm{r}}+\left[t_{\perp}c^{\dagger}_{\bm{r}+\bm{e_{x}}}\sigma_{z}c_{\bm{r}}+\sum_{\mu=y,z}\left(-{\rm
i}t_{\|}c^{\dagger}_{\bm{r}+\bm{e_{\mu}}}\sigma_{y}c_{\bm{r}}+t_{\|}^{\prime}c^{\dagger}_{\bm{r}+\bm{e_{\mu}}}\sigma_{z}c_{\bm{r}}\right)+\text{H.c.}\right]\right\\},$
(105)

where $\epsilon_{\bm{r}}$ takes real values and distributes uniformly in
$[-W/2,W/2]$, and the parameters are chosen to be $\Delta=0$, $t_{\|}=\cosh
g,t_{\|}^{\prime}=\sinh g$ ($g=0.2$), $t_{\perp}=1$. The Hamiltonian ${\cal
H}_{2}$ satisfies time-reversal symmetry ${\cal H}_{2}={\cal H}^{*}_{2}$ and
chiral symmetry ${\cal H}_{2}=-\sigma_{x}{\cal H}_{2}\sigma_{x}$, and hence
belongs to class BDI. In terms of Eq. (30), the chiral-symmetric Hamiltonian
${\cal H}_{2}$ is decomposed into the block-off diagonal structure in the
canonical basis where $\sigma_{x}$ is diagonalized. The right-upper part
$h_{2}$ of ${\cal H}_{2}$ in this basis is given by

$h_{2}=\sum_{\bm{r}=(r_{x},r_{y},r_{z})}\left\\{(\Delta+\epsilon_{\bm{r}})f^{\dagger}_{\bm{r}}f_{\bm{r}}+\left(t_{\perp}f^{\dagger}_{{\bm{r}}+{\bm{e}}_{x}}f_{\bm{r}}+\text{H.c.}\right)+\sum_{\mu=y,z}\left[(t^{\prime}_{\|}+t_{\|})f^{\dagger}_{{\bm{r}}+{\bm{e}}_{\mu}}f_{\bm{r}}+(t^{\prime}_{\|}-t_{\|})f^{\dagger}_{\bm{r}}f_{{\bm{r}}+{\bm{e}}_{\mu}}\right]\right\\}\,.$
(106)

Transposition exchanges $t^{\prime}_{\|}+t_{\|}$ and $t^{\prime}_{\|}-t_{\|}$.
Thus, as a unitary transformation in Eqs. (66) and (67), we can apply a
$\pi$-rotation around the $x$ axis,

$\displaystyle{\cal
U}_{(r_{x},r_{y},r_{z}|r^{\prime}_{x},r^{\prime}_{y},r^{\prime}_{z})}=\delta_{r_{x},r^{\prime}_{x}}\delta_{r_{y},-r^{\prime}_{y}}\delta_{r_{z},-r^{\prime}_{z}}.$
(107)

Since $\epsilon_{\bm{r}}$ is statistically equivalent for different lattice
points ${\bm{r}}$, an ensemble of $h_{2}$ defined in Eq. (106) is
statistically invariant under the combination of transposition and the
$\pi$-rotation. The symmetry in Eqs. (66) and (107) requires the LEs of
$h_{2}$ along the $x$ direction to come in opposite-sign pairs, giving rise to
$\nu_{x}=0$. The hopping along the $y$ direction and the hopping along the $z$
direction are symmetric in Eq. (106). Thus, after transposition, we can also
apply a mirror operation with respect to the plane with fixed $y+z$:

$\displaystyle{\cal
U}_{(r_{x},r_{y},r_{z}|r^{\prime}_{x},r^{\prime}_{y},r^{\prime}_{z})}=\delta_{r_{x},r^{\prime}_{x}}\delta_{r_{y},-r^{\prime}_{z}}\delta_{r_{z},-r^{\prime}_{y}}.$
(108)

Transposition exchanges $t^{\prime}_{\|}+t_{\|}$ and $t^{\prime}_{\|}-t_{\|}$
and the mirror operation puts them back. Thus, an ensemble of $h_{2}$ is
statistically invariant under the combination of transposition and the mirror
operation. The symmetry in Eqs. (66) and (108) requires the LEs of $h_{2}$
along the $r_{(0,1,-1)}\equiv r_{y}-r_{z}$ direction to come in opposite-sign
pairs, leading to $\nu_{y}=\nu_{z}$.

Figure 7: $2L^{2}$ Lyapunov exponents (LEs) of the right-upper part $h$ of the
nodal-line semimetal model ${\cal H}$ in the canonical basis
($t_{\perp}=3/10,t^{\prime}_{\|}=1,t_{\|}=1/4,\Delta=2$) calculated along the
$z$ direction with the quasi-1D geometry $L\times L\times L_{z}$.
Distributions of the $2L^{2}$ LEs are plotted as a function of the disorder
strength $W$. The LEs of the nodal-line semimetal model are the sum of the
$2L^{2}$ LEs of $h$ and their opposite-sign exponents.

(a) nodal-line semimetal with $\nu_{z}\neq 0$ and $\nu_{x}=\nu_{y}=0$
($g=0.22$)

(b) nodal-line semimetal with $\nu_{z}\neq 0$ and $\nu_{x}=\nu_{y}=0$ ($g=1$)

(c) nodal-line semimetal with with $\nu_{y}=\nu_{z}\neq 0$ and $\nu_{x}=0$

Figure 8:  Normalized localization length $\Lambda_{x}\equiv\xi_{x}/L$ along the $x$ direction as a function of the disorder strength $W$ for different system sizes $L$. (a, b) Nodal-line semimetal model in Eq. (35) with $\Delta=0$, $t_{\|}=\sinh g,t_{\|}^{\prime}=\cosh g$ [(a) $g=0.22$ and (b) $g=1$], $t_{\perp}=1$, and (c) nodal-line semimetal model in Eq. (105) with $\Delta=0$, $t_{\|}=\cosh g,t_{\|}^{\prime}=\sinh g$ ($g=0.2$), $t_{\perp}=1$. The black points are the raw data with the error bars. The solid lines for different $L$ and the dashed vertical line $W_{c}^{(x)}$ with the error bars are the results of the fitting according to Eqs. (20) and (21) with $(m,n)=(2,3)$. The dashed line $W_{c}^{(z)}$ is evaluated by the fitting of the Lyapunov exponents along the $z$ direction by Eqs. (86) and (87). Table 6: Polynomial fitting results of the normalized localization length $\Lambda_{x}\equiv\xi_{x}/L$ in the $x$ direction around its scaling invariant point of the different nodal-line semimetal models ${\cal H}$ [Eq. (35)] and ${\cal H}_{2}$ [Eq. (105)]. Both models are characterized by the parameters $\Delta=0$, $t_{\|}=\sinh g,t_{\|}^{\prime}=\cosh g$, $t_{\perp}=1$. The column “parameter” specifies the value of $g$. The data of $\Lambda_{x}$ in these models are calculated with the quasi-1D geometry ($L\times L\times L_{x}$, $L_{x}\gg L$). The range of the system size $L$, critical disorder strength $W^{(x)}_{c}$, critical exponent $\nu$, scaling dimension $-y$ of the least irrelevant variable, critical localization length $\Lambda_{c}$, the goodness of fitting (GOF), and the fitting order $(m,n)$ in Eqs. (20) and (21) are shown. The square brackets are the 95% confidence error bars determined by the synthetic data. Note that $\Lambda_{c}$ is non-universal because of the spatial anisotropy. model | parameter | $L$ | $m$ | $n$ | GOF | $W^{(x)}_{c}$ | $\nu$ | $y$ | $\Lambda_{c}$
---|---|---|---|---|---|---|---|---|---
${\cal H}$ | $g=0.22$ | 20-34 | 2 | 3 | 0.61 | 21.688[21.647,21.729] | 0.866[0.792,0.916] | 1.840[1.608,2.130] | 0.239[0.232,0.246]
${\cal H}$ | $g=0.22$ | 20-34 | 3 | 3 | 0.58 | 21.730[21.662,21.803] | 0.860[0.679,0.918] | 1.619[1.284,2.031] | 0.232[0.218,0.243]
${\cal H}$ | $g=1$ | 18-28 | 2 | 3 | 0.14 | 41.284[41.209,41.365] | 0.787[0.754,0.818] | 2.256[2.013,2.516] | 0.110[0.108,0.113]
${\cal H}$ | $g=1$ | 18-28 | 3 | 3 | 0.24 | 41.358[41.264,41.463] | 0.785[0.751,0.817] | 2.088[1.816,2.381] | 0.108[0.104,0.111]
${\cal H}_{2}$ | $g=0.2$ | 20-28 | 2 | 3 | 0.32 | 23.214[23.146,23.334] | 0.857[0.712,0.920] | 2.595[1.845,3.300] | 0.223[0.208,0.231]
${\cal H}_{2}$ | $g=0.2$ | 20-28 | 3 | 3 | 0.31 | 23.229[23.131,23.422] | 0.855[0.691,0.921] | 2.513[1.593,3.418] | 0.222[0.197,0.232]
Table 7: Comparison of the two critical disorder strengths, $W^{(x)}_{c}$ and $W^{(z)}_{c}$, in the nodal-line semimetal models ${\cal H}$ [Eq. (35)] and ${\cal H}_{2}$ [Eq. (105)]. The parameters of these models are the same as in Table 6. Both models are characterized by the parameters $\Delta=0$, $t_{\|}=\sinh g,t_{\|}^{\prime}=\cosh g$, $t_{\perp}=1$. The column “parameter” specifies the value of $g$. In the two models, the non-localized regions extend from $W=0$ to $W=W^{(z)}_{c}$, and the quasi-localized phases extend from $W^{(x)}_{c}$ to $W^{(z)}_{c}$. symmetry class |  model |  parameter | $W_{c}^{(x)}$ | $W_{c}^{(z)}$
---|---|---|---|---
BDI | ${\cal H}$ | $g=0.22$ |  21.73[21.66,21.80] |  23.4[23.3,23.5]
BDI | ${\cal H}_{2}$ | $g=1$ | 41.36[41.26,41.46] | 45.3[45.2,45.4]
BDI | ${\cal H}_{2}$ | $g=0.2$ | 23.23[23.13,23.42] | 24.7[24.6,24.8]

For the directions with non-zero weak topological indices [i.e., $z$ direction
in Eq. (45) and $y,z$ directions in Eq. (106)], we calculate the LEs of the
right-upper part [i.e., $h$ in Eq. (45) and $h_{2}$ in Eq. (106)] as a
function of the disorder strength $W$. In all these models, the distributions
of the LEs for large $L$ show the $W$-dependence described in Fig. 7. For
$W=0$, the $2L^{2}$ LEs form a continuum spectrum in the large $L$ limit,
including zero $\gamma=0$. When $W$ increases, the spectrum splits into the
two continuous spectra. For $W>W^{(z)}_{c}$, all the $L^{2}$ LEs in the lower
spectrum become negative. The non-localized region extends from $W=0$ to
$W=W^{(z)}_{c}$, while the Anderson insulator phase appears in
$W>W^{(z)}_{c}$. $W^{(z)}_{c}$ is determined by the finite-size scaling of
$\gamma^{(1)}_{\rm max}$ as in Sec. I.7, which is summarized in Table 7.

For the directions with zero weak topological indices [i.e., $x,y$ directions
in Eq. (45) and $x$ direction in Eq. (106)], we calculate the localization
length $\xi_{x}$ in the quasi-1D geometry ($L_{x}\times L\times L$, $L_{x}\gg
L$). Figure 8 shows the normalized localization length
$\Lambda_{x}\equiv\xi_{x}/L$ as a function of the disorder strength $W$. In
all these models, $\Lambda_{x}$ shows scale-invariant behavior inside the non-
localized region. Fitting $\Lambda_{x}$ around the scale-invariant points by
the scaling functions of Eqs. (20) and (21), we evaluate the critical exponent
$\nu$ and the critical disorder strength $W_{c}^{(x)}$, as summarized in Table
6. The critical disorder strength $W_{c}^{(x)}$ is far below $W_{c}^{(z)}$,
demonstrating the presence of the quasi-localized phases for
$W_{c}^{(x)}<W<W_{c}^{(z)}$ in these three models (Table 7). The evaluated
critical exponents are consistent with the critical exponent
$\nu=0.820[0.783,0.846]$ shown in the main text for the nodal-line semimetal
in class BDI. This consistency suggests that all the phase transitions between
the metal and quasi-localized phases in the 3D models in symmetry class BDI
are of the same nature.

Note also that $\nu$ in the nodal-line semimetal model of Eq. (45) with
$t_{\|}=\cosh g,t_{\|}^{\prime}=\sinh g$ ($g=0.22$) shows the larger error
bars in their fitting results (see Table 6). These larger error bars may stem
from a severe crossover effect. For $t_{\|}^{\prime}=\sinh g=0$, the nodal-
line semimetal model has an extra unitary symmetry ${\cal H}=\sigma_{z}{\cal
H}\sigma_{z}$. The Hamiltonian can be block-diagonalized into two parts, and
each block belongs to orthogonal class. For nonzero but small $g$, this
unitary symmetry is only weakly broken. Thus, the finite-size systems with
smaller $g$ must suffer from a stronger crossover effect.

####  I.8.3 Weak topological insulators and ordinary insulators in class BDI

Zero-energy states of the nodal-line semimetal model in Eq. (35) can be either
in a topological insulator state with $(\nu_{x},\nu_{y},\nu_{z})=(0,0,1)$ or
in an ordinary insulator state with $(\nu_{x},\nu_{y},\nu_{z})=(0,0,0)$,
depending on its tight-binding parameters. For simplicity, we assume
$\Delta,t_{\perp},t_{\|}>0$ in Eq. (35). For
$\Delta+4t_{\perp}<2|t_{\|}^{\prime}|$
($\Delta-4t_{\perp}>2|t_{\|}^{\prime}|$), the zero-energy states of ${\cal H}$
in Eq. (35) are in the topological (ordinary) insulator state in the clean
limit ($W=0$).

In Ref. [50], the localization lengths of Eq. (35) along the $x$ direction
were calculated with the quasi-1D geometry ($L_{x}\times L\times L$, $L_{x}\gg
L$), and the two consecutive disorder-driven phase transitions were identified
for the following set of parameters:

$t_{\perp}=3/10,\quad t^{\prime}_{\|}=1,\quad
t_{\|}=1/4,\quad\Delta=1/2\quad\text{(topological insulator in the clean
limit).}$ (109)

The two phase transitions are i) a transition from the topological insulator
phase to the diffusive metal phase at $W^{(x)}_{c,1}=3.135\\!\ [3.132,3.138]$
and ii) a transition from the diffusive metal phase to the Anderson insulator
phase at $W^{(x)}_{c,2}=11.96\\!\ [11.92,12.02]$, respectively. In addition,
Ref. [50] studied another parameter set,

$\displaystyle t_{\perp}=3/10,\quad t^{\prime}_{\|}=1,\quad
t_{\|}=1/4,\quad\Delta=4\quad\text{(ordinary insulator in the clean limit),}$
(110)

where a disorder-driven phase transition from the 3D ordinary band insulator
phase to the diffusive metal phase [Fig. 9(b)] was found at
$W^{(x)}_{c,3}=4.76\\!\ [4.75,4.77]$. The normalized localization length
$\Lambda_{x}\equiv\xi_{x}/L$ shows scale-invariant behavior at these critical
disorder strengths. From the finite-size scaling analyses, it was clarified
that these phase transitions are universally characterized by the same
critical exponent $\nu=0.82\pm 0.04$ [50], which is consistent with the
evaluations of the disordered nodal-line semimetal models studied in the main
text.

(a) $t_{\perp}=3/10,t^{\prime}_{\|}=1,t_{\|}=1/4,\Delta=1/2$

(b) $t_{\perp}=3/10,t^{\prime}_{\|}=1,t_{\|}=1/4,\Delta=4$

(c) $\gamma-W$ in (a)

(d) $\gamma-W$ in (b)

Figure 9: (a,b) Schematic phase diagrams of ${\cal H}$ in different parameter
sets, (a) weak topological and (b) ordinary insulator sides.
$W_{c,i}^{(x)}(i=1,2,3)$ stand for critical points of the Anderson transitions
in the $x$ or $y$ direction [50]. The weak topological index $\nu_{z}$ is
finite for $W_{c,1}^{(z)}<W<W_{c,2}^{(z)}$ of (a) and for $W^{(z)}_{c,3}<W$ of
(b). The $2L^{2}$ Lyapunov exponents of $h$ in the same parameter sets as (a)
and (b) are shown in (c) and (d), respectively.

For these two sets of parameters, we study the localization length and the
winding number along the $z$ direction. We calculate the LEs of the right-
upper part $h$ of ${\cal H}$ in the canonical basis [i.e., Eq. (45)] with the
quasi-1D geometry ($L\times L\times L_{z}$, $L_{z}\gg L$). For the parameters
in Eq. (109), the LEs show a $W$-dependence described as Fig. 9(c). The
topological insulator is stable under weak disorder. For $W<W^{(z)}_{c,1}$,
all the $2L^{2}$ LEs are positive, and the localization length is finite. The
winding number $w_{z}$ along the $z$ direction is $L^{2}$, giving rise to
$\nu_{z}=1$. The non-localized region appears from $W=W^{(z)}_{c,1}$ to
$W=W^{(z)}_{c,2}$ ($>W^{(z)}_{c,1}$), where a continuous spectrum of LEs
includes zero $\gamma=0$. When $W$ increases from $W^{(z)}_{c,1}$ to
$W^{(z)}_{c,2}$, $L^{2}$ positive LEs cross zero and become negative;
$\nu_{z}$ changes from 1 to 0. For $W>W^{(z)}_{c,2}$, the $L^{2}$ LEs are
positive and the other $L^{2}$ LEs are negative, leading to $\nu_{z}=0$. For
the parameters in Eq. (110), the LEs shows a $W$-dependence described as Fig.
9(d), where the non-localized region appears from $W=W^{(z)}_{c,3}$ to
$W=W^{(z)}_{c,4}$ ($>W^{(z)}_{c,3}$). In the non-localized region, the number
of positive LEs is greater than the number of negative ones, leading to
$\nu_{z}>0$.

From the finite-size scaling analyses of the minimal or maximal LEs by Eqs.
(86), (87), or (101), we determine the phase boundaries of the non-localized
regions as $W^{(z)}_{c,1}=3.08\\!\ [3.07,3.09]$, $W^{(z)}_{c,2}=13.3\\!\
[13.2,13.4]$, and $W^{(z)}_{c,3}=4.56\\!\ [4.55,4.57]$ (see Table 8). From a
comparison of these numbers with $W^{(x)}_{c}$ obtained in Ref. [50] (see
Table 9), we conclude that the quasi-localized phases appear inside the non-
localized regions:

Table 8: Finite-size scaling analyses of $\gamma^{(1)}_{\rm max}(W,L)$ or
$\gamma^{(1)}_{\rm min}(W,L)$ for the several disorder strength $W$ around
$W=W_{c}^{(z)}$ for the disordered topological insulator model with the
parameters in Eq. (109) (shown as “P1” in the “parameter set”), and the
disordered ordinary insulator model with the parameters in Eq. (110) (shown as
“P3” in the “parameter set”). The square brackets are the 95% confidence error
bars determined by the Monte Carlo analyses.

parameter set | $W$ | $L$ | $\gamma_{\rm min}(W,L)$ | $a$ |  GOF
---|---|---|---|---|---
P1 |  3.07 |  24 - 60 |  0.0013 [0.0012,0.0014] |  0.050 [0.046,0.053] | 0.64
P1 | 3.08 | 24 - 60 | -0.0002 [-0.0003,-0.0000] | 0.051 [0.047,0.054] | 0.90
P1 | 3.09 | 24 - 60 | -0.0015 [-0.0016,-0.0014] | 0.048 [0.045,0.051] | 0.46

parameter set | $W$ | $L$ | $\gamma^{(1)}_{\rm max}(W,L)$ | $a$ |  GOF
---|---|---|---|---|---
P1 |  13.2 |  14 - 28 |  0.006 [0.002,0.010] |  -0.967 [-1.048,-0.889] | 0.15
P1 | 13.3 | 14 - 28 | 0.0007 [-0.003,0.004] | -0.954 [-1.027,-0.880] | 0.74
P1 | 13.4 | 14 - 28 | -0.006 [-0.010,-0.002] | -0.951 [-1.026,-0.870] | 0.11
P3 | 4.60 | 14 - 28 | -0.007 [-0.011,-0.0028] | -1.177 [-1.257,-1.097] | 0.74
P3 | 4.62 | 14 - 28 | -0.003 [-0.007,0.001] | -1.091 [-1.163,-1.009] | 0.19
P3 | 4.63 | 14 - 28 | 0.007 [0.002,0.0109] | -1.150 [-1.231,-1.064] | 0.79