You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
author:
- 'Yutaka Akagi$^1$[^1], Hosho Katsura$^1$, and Tohru Koma$^2$'
title: 'A New Numerical Method for $\mathbb{Z}_2$ Topological Insulators with Strong Disorder'
---

Since the discovery of $\mathbb{Z}_2$ topological insulators by Kane and Mele [@Kane05], many methods have been proposed to compute the $\mathbb{Z}_2$ index. In particular, Fu and Kane found that the calculation of the $\mathbb{Z}_2$ index can be considerably simplified in a system with inversion symmetry [@Fu07]. However, for disordered systems, a numerical determination of the $\mathbb{ Z}_2$ index is still very challenging. Roughly speaking, three types of numerical methods have been proposed so far for this problem:

$\bullet$ The first one was proposed by Kane and Mele themselves [@Kane05], and later some modifications were introduced. Basically, for a given model, the $\mathbb{ Z}_2$ index is defined by a certain Pfaffian with twisted boundary conditions [@Niu85]. The methods were applied to class AII models [@Schnyder08; @Kitaev09; @Ryu10] in two and three dimensions with or without a certain inversion symmetry [@Essin07; @Guo10; @Leung12].

$\bullet$ The second one is based on a scattering matrix theory [@Fulga12]. The $\mathbb{ Z}_2$ indices which are defined by the scattering matrices were numerically computed for models in the classes, AII and DIII, in two and three dimensions [@Fulga12; @Sbierski14].

$\bullet$ The third one was proposed by Loring and Hastings [@Loring10]. The $\mathbb{ Z}_2$ indices are defined by Bott indices which are introduced as an obstruction to approximating almost commuting matrices. For models in the class AII in two and three dimensions, the $\mathbb{ Z}_2$ indices were numerically computed [@Loring10; @Loring15]. For systems with chiral symmetry, Loring and Schulz-Baldes [@Loring17] proposed a numerical method to obtain the values of Bott indices.

In this paper, we propose an alternative method to numerically calculate the $\mathbb{ Z}_2$ indices for disordered topological insulators in arbitrary dimensions. The method is based on the index formulae which were derived in Refs. \[15, 16\]. and has the following two advantages: (i) There is no need to take an average of the $\mathbb{ Z}_2$ index over random variables in a model. (ii) The $\mathbb{ Z}_2$ index is determined by the discrete spectrum of a certain compact operator with a supersymmetric structure. The latter makes it possible to numerically determine the $\mathbb{Z}_2$ index highly efficiently.

Our method is demonstrated for Bernevig-Hughes-Zhang (BHZ) [@Bernevig06; @Yamakage12; @Yamakage13] and Wilson-Dirac [@Wilson74; @Qi08; @Kobayashi13] models whose topological phases are characterized by a $\mathbb{ Z}_2$ index of the class AII in two and three dimensions, respectively. In consequence, the method enables us to determine all of the values of the $\mathbb{ Z}_2$ indices of the strong and weak topological insulators, and the normal insulator phases in the phase diagrams. These values of the $\mathbb{ Z}_2$ indices completely coincide with the predictions in previous studies using a reliable transfer-matrix method [@Yamakage12; @Yamakage13; @Kobayashi13].

To begin with, we introduce two Dirac operators as [@Katsura16a; @Katsura16b] $$\mathcal{D}_{\bm a}({\bm x}) := \frac{x_1 + i x_2 - (a_1 + i a_2)}{| x_1 + i x_2 - (a_1 + i a_2) |}$$ in two dimension (2D), and $$D_{\bm a}({\bm x}) := \frac{1}{|{\bm x} - {\bm a}|} ({\bm x} - {\bm a}) \cdot \mbox{\boldmath $\sigma$}$$ in three dimensions (3D), where ${\bm x} =$ is the position operator and ${\bm a} =$ is a vector for $d=2,3$ \[see Fig. \[Fig:PD\_BHZ\](a)\]. The three-component vector $\mbox{\boldmath $\sigma$}$ is defined by $\mbox{\boldmath $\sigma$}=\mbox{$(\sigma_1, \sigma_2, \sigma_3)$}$ whose components are given by Pauli matrices $\sigma_i$, each of which acts on an auxiliary Hilbert space $\mathbb{C}^2$. The whole Hilbert space is given by the tensor product of the auxiliary $\mathbb{C}^2$ and the original Hilbert space for the Hamiltonian of tight-binding models which we will consider shortly.

Next, we define the $\mathbb{Z}_2$ index for an infinite-volume system which is a tight-binding model on a square $\mathbb{Z}^2$ or cubic lattice $\mathbb{Z}^3$ [@comment]. Let $P_{\rm F}$ be the projection operator onto the states below the Fermi energy $E_{\rm F}$. The difference of two projections is defined as [@Avron94b; @Katsura16b] $$\begin{aligned}
A := \begin{cases}
    P_{\rm F} - \mathcal{D}_{\bm a}^* P_{\rm F} \mathcal{D}_{\bm a} & {\rm in \: \: 2D} \\
    P_{\rm F} - D_{\bm a} P_{\rm F}D_{\bm a} & {\rm in \: \: 3D},
  \end{cases}
\label{eq:A_op}\end{aligned}$$ where $\mathcal{D}_{\bm a}^*$ is the adjoint of the Dirac operator $\mathcal{D}_{\bm a}$. Then, the $\mathbb{Z}_2$ index $\nu$ is defined as $$\begin{aligned}
  \nu = {\rm dim} \; {\rm ker}\; (A - 1) \; \; {\rm modulo} \; 2.
  \label{eq:Z2-index}\end{aligned}$$ When the Fermi energy $E_{\rm F}$ lies in a spectral gap or a mobility gap, the $\mathbb{Z}_2$ index is known to be well defined [@Katsura16b]. In the following, we will consider such situations.

Now we describe our numerical scheme for calculating the $\mathbb{Z}_2$ index. Let $\Lambda$ and $\Omega$ be two finite regions satisfying $\Omega\subset\Lambda\subset\mathbb{R}^d$ and ${\bm a}\in\Omega$ as in Fig. \[Fig:PD\_BHZ\]. First, we approximate the Fermi projection $P_{\rm F}$ in the infinite volume by the corresponding Fermi projection in the finite region $\Lambda$. More precisely, the approximate one is given by $$\begin{aligned}
P_{\rm F}^{(\Lambda)}:= \sum_{E_n \le E_{\rm F}} \ket{n} \bra{n},\end{aligned}$$ where $\ket{n}$ are eigenstates of the tight-binding Hamiltonian $\mathcal{H}$ on $\Lambda$ which we consider, and we have denoted by $E_n$ the corresponding eigenvalues. To avoid the presence of gapless edge/surface states, we impose periodic boundary conditions for the Hamiltonian in practical numerical calculations. For the operator $A$ in Eq. (\[eq:A\_op\]), we replace the Fermi projection $P_{\rm F}$ by $P_{\rm F}^{(\Lambda)}$, and write $A^{(\Lambda)}$ for the approximate one. Further, we restrict the operator $A^{(\Lambda)}$ to the subregion $\Omega$ as $$\label{AOmegaLambda}
A_{\Omega}^{(\Lambda)} := \chi_\Omega A^{(\Lambda)} \chi_\Omega,$$ where $\chi_\Omega$ is the characteristic function of the subregion $\Omega$.

Under the above gap assumption, the operator $A$ in Eq. (\[eq:A\_op\]) is compact even in the infinite volume limit. Therefore, the spectrum of $A$ is discrete and has no accumulation point except for zero. This implies that an eigenstate of $A$ is localized if the corresponding eigenvalue is not equal to zero. Let $\lambda\ne 0$ be an eigenvalue of $A$. From these observations, it is clear that the eigenvalue $\lambda$ can be approximated by an eigenvalue $\lambda'$ of the approximate operator $A_\Omega^{(\Lambda)}$ if the subregion $\Omega$ is sufficiently large. In addition to this, the cutoff function $\chi_\Omega$ in Eq. (\[AOmegaLambda\]) enables us to remove the boundary effects due to the finite size of the region $\Lambda$.

As the first demonstration, we compute the $\mathbb{ Z}_2$ index of the BHZ model [@Bernevig06; @Yamakage12; @Yamakage13] with disorder on a square lattice $\mathbb{Z}^2$. The Hamiltonian is written as $$\begin{aligned}
    \mathcal{H}_{\rm 2D} \! = \! \sum_{\bm x} \! \sum_{k=1,2} \! \bigl[ c^{\dagger}_{\bm x} t_k c_{{\bm x}
+{\bm e}_k} \! + \!{\rm h.c.}  \bigr]
                \! + \! \sum_{\bm x} c^{\dagger}_{\bm x} (\tau_0 \! \otimes \! \epsilon_{\bm x}) c_{\bm x},
\label{Eq:BHZ}\end{aligned}$$ where $c_{\bm x}=$ , and $c_{{\bm x} i \alpha}$ is the annihilation operator of an electron with spin $\alpha$ and orbital $i$ at site ${\bm x}$. The hopping matrices, $t_1$ and $t_2$, are given by $t_1= g_1\tau_0 \otimes s_3 -\frac{i}{2} g_2 \tau_0 \otimes s_2 + \frac{i}{2} g_3 \tau_2 \otimes s_3$, and $t_2= g_1\tau_0 \otimes s_3 -\frac{i}{2} g_2 \tau_3 \otimes s_1 - \frac{i}{2} g_3 \tau_2 \otimes s_0$, where $\tau_k$ and $s_k$ are Pauli matrices for the orbital and the spin, respectively. Here, $g_1, g_2$ and $g_3$ are real parameters. The on-site potential $\epsilon_{\bm x}$ is given by $\epsilon_{\bm x}=$ , where $W_{\bm x}^{+}$ and $W_{\bm x}^{-}$ are a random potential whose distribution is uniform in the interval with a positive parameter $W$. ${\bm e}_1$ and ${\bm e}_2$ are the unit vectors in the $x$ and $y$ directions, respectively. As is well known [@Bernevig06; @Yamakage12; @Yamakage13], this Hamiltonian (\[Eq:BHZ\]) belongs to the symmetry class AII. We set the Fermi energy which is located at the center of the energy gap.

In the following, we write $\lambda_i$ for the [*i*]{}-th eigenvalue of the operator $A_{\Omega}^{(\Lambda)}$ of Eq. (\[AOmegaLambda\]) in descending order including multiplicity [@comment2]. We note that the spectrum of $A_{\Omega}^{(\Lambda)}$ is included in the interval $[-1,1]$.

Before showing our numerical results, we abbreviate the topological and the ordinary insulator phases as TI and OI, respectively, in the phase diagram [@Yamakage12; @Yamakage13], and write $\nu$ for the value of the $\mathbb{ Z}_2$ index. Figure \[Fig:PD\_BHZ\](c) and (d) show, respectively, $\lambda_1$ and as a function of the mass $\Delta$ and the strength $W$ of disorder. In the region TI, our numerical results are satisfactory because and . Actually, these imply $\nu=1$, i.e., the phase is topological as we expected. In the region OI, $\lambda_1$ is significantly different from 1, and thus $\nu=0$. These results show that the $\mathbb{ Z}_2$ index enables us to distinguish TI from OI. In OI phase, one notices $\lambda_1\simeq\lambda_2\lesssim 0.2$ as seen in Fig. \[Fig:PD\_BHZ\](d). This degeneracy is nothing but a consequence of the two symmetries [@Katsura16a; @Katsura16b], the time-reversal symmetry of the Hamiltonian and the supersymmetric structure of the operator $A$. This kind of degeneracy is very useful to determine the $\mathbb{ Z}_2$ index. In the region $W \lesssim 6$, our results are in good agreement with the previous results [@Yamakage13] which were obtained by using a transfer-matrix method. In the region with a sufficiently large $W$, Anderson localization is expected to occur. For the intermediate critical region between these two regions, our method is not under control because the existence of a significantly nonvanishing spectral or mobility gap cannot be expected. In fact, our numerical results in this region show large fluctuations in both $\lambda_1$ and $\lambda_1-\lambda_2$.

![(Color online). (a) The location of ${\bm a}$ is chosen to be the center of the finite lattice $\Lambda$ of linear size $L$. (b) The subregion $\Omega$ is chosen so that its center is the same as ${\bm a}$ and it does not overlap with the boundary of $\Lambda$. (c) and (d), respectively, show $\lambda_1$ and in the BHZ model as a function of the mass $\Delta$ and the disorder strength $W$. The obtained value $\nu$ of the $\mathbb{Z}_2$ index is indicated in both OI and TI phases. The values of the parameters used are , , , and the system size is . The curves of the phase boundaries with the dots are plotted by using the results in Ref. \[19\]. []{data-label="Fig:PD_BHZ"}](Fig1.eps){width="1.0\columnwidth"}

The second example is the Wilson-Dirac model [@Kobayashi13] with disorder on the cubic lattice $\mathbb{Z}^3$. The Hamiltonian is written as $$\begin{aligned}
    \mathcal{H}_{\rm 3D} = \mathcal{H}_{\sf 0} + \mathcal{H}_{\sf hop} + \mathcal{H}_{\sf dis},
\label{Eq:WD}\end{aligned}$$ where $$\begin{aligned}
    \mathcal{H}_{\sf 0} \! \! \! \! \! &=& \! \! \! \! \! \sum_{\bm x} \! \! \sum_{k=1,2,3} \! \bigl[\frac{i t}{2} c^{\dagger}_{{\bm x}+{\bm e}_k} {\alpha}_k c_{\bm x} 
                         \! - \! \frac{m_2}{2} c^{\dagger}_{{\bm x}+{\bm e}_k} \beta c_{\bm x} \! + \! {\rm h.c.} \bigr] \notag\\ 
                         &+& \! \! \! \! \! (m+3m_2) \sum_{\bm x} c^{\dagger}_{\bm x} \beta c_{\bm x}, \\
    \mathcal{H}_{\sf hop} \! \! \! \! \! &=& \! \! \! \! \! \sum_{\bm x} \sum_{k=1,2,3} \bigl[ t_0 c^{\dagger}_{{\bm x}+{\bm e}_k} c_{\bm x} + {\rm h.c.} \bigr], \\
    \mathcal{H}_{\sf dis} \! \! \! \! \! &=& \! \! \! \! \! \sum_{\bm x} v_{\bm x} c^{\dagger}_{\bm x} c_{\bm x}.\end{aligned}$$ Here, $c_{\bm x}=$ , the vector ${\bm e}_k$ is the unit vector in $k=x,y,z$ direction, and $\alpha_k =$ , $\beta =$ ; $v_{\bm x}$ is the on-site random potential whose distribution is uniform in the interval with a positive parameter $W$. This Hamiltonian (\[Eq:WD\]) belongs to the symmetry class AII for or . In the following, we will treat the case with .

![(Color online). (a) and (b), respectively, show $\lambda_1$ and in the Wilson-Dirac-type model as a function of the mass parameter $m_0/m_2$ and the disorder strength $W/m_2$. The numerical values $\nu$ of the $\mathbb{ Z}_2$ index are indicated in the phases WTI, STI and OI. The parameters used are , , , and the system size is . The curves of the phase boundaries with the dots are plotted by using the results in Ref. \[22\].[]{data-label="Fig:PD_Wilson-Dirac"}](Fig2.eps){width="1.0\columnwidth"}

Similarly, we abbreviate the weak, strong topological, the ordinary insulator and the diffusive metal phases as WTI, STI, OI, and M, respectively [@Kobayashi13], and write $\nu$ for the value of the $\mathbb{ Z}_2$ index [@comment3]. Figure \[Fig:PD\_Wilson-Dirac\] shows $\lambda_1$ and as a function of the mass parameter $m_0/m_2$ and the strength $W/m_2$ of disorder. In the region OI, $\nu=0$ because . In the region STI, and , and hence $\nu=1$. In the region WTI, but $\lambda_3$ is significantly different from 1 \[see Fig. \[Fig:cross-section\](b)\], and hence $\nu=0$ because the multiplicity of the eigenvalue is two. As for the region M, we cannot expect our method to be effective because the spectral or mobility gap is expected to vanish if the metallic character of the spectrum is true. To summarize, as seen in Fig. \[Fig:PD\_Wilson-Dirac\], our numerical results for the $\mathbb{ Z}_2$ index are in good agreement with the predictions of Ref. \[22\]. In particular, the phase boundaries between WTI and STI, and between STI and OI are considerably sharp. Moreover, although our method is expected to be useless in the metallic phase, there do not appear fluctuations like those in the two-dimensional case.

Although the $\mathbb{Z}_2$ index vanishes in both the OI and WTI phases, there is a definite difference between them as follows: OI phase yields no eigenvalue while WTI phase yields the eigenvalue which is doubly degenerate. If the eigenvalue is related to surface states, one can expect, from our numerical results, that the multiplicity of coincides with the number of Dirac cones which appear on the surface of the system. In fact, it was numerically observed in a generalized Kane-Mele model on a diamond lattice that WTI (STI) phase exhibits two surface Dirac cones (odd number of Dirac cones) [@Fu07b].

![(Color online). (a) $W/m_2$ and (b) $m_0/m_2$ dependences of the eigenvalues $\lambda_i$, $i=1,2,\ldots,5$, of the operator $A_{\Omega}^{(\Lambda)}$ for the system size . (a) and (b) correspond to cross-sections of Fig. \[Fig:PD\_Wilson-Dirac\] for and , respectively. The inset in (a) shows $\lambda_1$ and $\lambda_2$ as a function of the inverse volume of the system ($1/L^3$) for fixed parameters, and . The green (yellow) line is a fit to the numerical data of $\lambda_1$ ($\lambda_2$). []{data-label="Fig:cross-section"}](Fig3.eps){width="0.8\columnwidth"}

Figure \[Fig:cross-section\] shows $W/m_2$ and $m_0/m_2$ dependences of the eigenvalues $\lambda_i$, $i=1,2,\ldots,5$, of the operator $A_{\Omega}^{(\Lambda)}$. According to Ref. \[22\], STI emerges for and for . One can see that $\lambda_1$ is significantly different from $\lambda_2$, and that $\lambda_2$ and $\lambda_3$ are degenerate, and similarly, $\lambda_4$ and $\lambda_5$ are degenerate. As mentioned in the case of two dimensions, this even degeneracy is a consequence of the two symmetries [@Katsura16a; @Katsura16b], the time-reversal symmetry of the Hamiltonian and the supersymmetric structure of the operator $A$. For the region in Fig. \[Fig:cross-section\](a), the diffusive metallic phase appears. In this region, one can see that the above double degeneracy in the spectrum of $A_{\Omega}^{(\Lambda)}$ is lifted due to the vanishing of the spectral or mobility gap in the metallic phase.

The inset of Fig. \[Fig:cross-section\](a) shows system-size dependences of $\lambda_1$ and $\lambda_2$ for fixed parameters and . As the system size increases, $\lambda_1$ and $\lambda_2$ converge to $1$ and about $0.7$, respectively. Thus in the infinite-volume limit, we can clearly conclude that the $\mathbb{ Z}_2$ index $\nu$ is unity. One can perform a similar analysis and extract the eigenvalues of $A$ in the infinite-volume limit for other values of the parameters.

As seen in Fig. \[Fig:cross-section\](b), the double degeneracy of appears in the region of WTI where the corresponding parameter satisfies . Clearly, one can see that $\lambda_3$ is separated from , and hence the $\mathbb{ Z}_2$ index $\nu$ is equal to zero. Thus, if the first and the second eigenvalues, $\lambda_1$ and $\lambda_2$, are degenerate, information about the third eigenvalue $\lambda_3$ is useful to determine the $\mathbb{ Z}_2$ index.

Summary: We have presented our new method for numerically calculating the $\mathbb{ Z}_2$ indices of topological insulators with strong disorder. Our method has the following two advantages: (i) There is no need to take an average of the $\mathbb{ Z}_2$ index over random variables in a model. (ii) The $\mathbb{ Z}_2$ index is determined by the discrete spectrum of a certain compact operator with a supersymmetric structure. These properties make it possible to numerically determine the $\mathbb{Z}_2$ index highly efficiently. In order to check the effectiveness of our method, we have demonstrated that all of the numerical values of the $\mathbb{ Z}_2$ indices completely coincide with the predictions in previous studies using a reliable transfer-matrix method [@Yamakage12; @Yamakage13; @Kobayashi13] for the two-dimensional Bernevig-Hughes-Zhang and the three-dimensional Wilson-Dirac models. Thus, the strong topological insulator phases can be characterized by the $\mathbb{ Z}_2$ index in the index formulae [@Katsura16a; @Katsura16b], and can be clearly distinguished from other phases. We believe that the good agreement between the two different approaches is one of the steps toward the understanding of the nature of $\mathbb{Z}_2$ topological insulators although we cannot definitely compare our method with other approaches mentioned at the beginning of the present paper. Finally, we remark that the generalization of our method to models in other symmetry classes in arbitrary dimensions is straightforward.

[101]{}

C. L. Kane and E. J. Mele, *${{ Z}}_2$ Topological Order and the Quantum Spin ${{\rm H}}$all Effect*, Phys. Rev. Lett. **95**, 146802 (2005).

L. Fu and C. L. Kane, *Topological insulators with inversion symmetry*, Phys. Rev. B **76**, 045302 (2007).

Q. Niu, D. J. Thouless, and Y.-S. Wu, *Quantized ${{\rm H}}$all conductance as a topological invariant*, Phys. Rev. B **31**, 3372 (1985).

A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, *Classification of topological insulators and superconductors in three spatial dimensions*, Phys. Rev. B **78**, 195125 (2008).

A. Kitaev, *Periodic table for topological insulators and superconductors*, AIP Conference Proceedings **1134**, 22 (2009).

S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, *Topological insulators and superconductors: tenfold way and dimensional hierarchy*, New J. Phys. **12**, 065010 (2010).

A. M. Essin and J. E. Moore, *Topological insulators beyond the ${{\rm B}}$rillouin zone via ${{\rm C}}$hern parity*, Phys. Rev. B **76**, 165307 (2007).

H.-M. Guo, *Topological invariant in three-dimensional band insulators with disorder*, Phys. Rev. B **82**, 115122 (2010).

B. Leung and E. Prodan, *Effect of strong disorder in a three-dimensional topological insulator: Phase diagram and maps of the ${{ Z}}_2$ invariant*, Phys. Rev. B **85**, 205136 (2012).

I. C. Fulga, F. Hassler, and A. R. Akhmerov, *Scattering theory of topological insulators and superconductors*, Phys. Rev. B **85**, 165409 (2012).

B. Sbierski and P. W. Brouwer, *${{ Z}}_2$ phase diagram of three-dimensional disordered topological insulators via a scattering matrix approach*, Phys. Rev. B **89**, 155311 (2014).

T. A. Loring and M. B. Hastings, *Disordered topological insulators via ${C}^*$-algebras*, EPL (Europhysics Lett.) **92**, 67004 (2010).

T. A. Loring, *K-theory and pseudospectra for topological insulators*, Ann. Phys. **356**, 383 (2015).

T. A. Loring and H. Schulz-Baldes, *Finite volume calculation of [K-theory]{} invariants*, Preprint arXiv:1701.07455 (2017).

H. Katsura and T. Koma, *The ${{ Z}}_2$ index of disordered topological insulators with time reversal symmetry*, J. Math. Phys. **57**, 021903 (2016).

H. Katsura and T. Koma, *The Noncommutative Index Theorem and the Periodic Table for Disordered Topological Insulators and Superconductors*, Preprint arXiv:1611.01928 (2016).

B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, *Quantum Spin ${{\rm H}}$all Effect and Topological Phase Transition in ${{\rm H}}$g${{\rm T}}$e Quantum Wells*, Science **314**, 1757 (2006).

A. Yamakage, K. Nomura, K.-I. Imura, and Y. Kuramoto, *${{ Z}}_2$ Topological Anderson Insulator*, J. Phys.: Conf. Ser. **400**, 042070 (2012).

A. Yamakage, K. Nomura, K.-I. Imura, and Y. Kuramoto, *Criticality of the metal-topological insulator transition driven by disorder*, Phys. Rev. B **87**, 205141 (2013).

K. G. Wilson, *Confinement of quarks*, Phys. Rev. D **10**, 2445 (1974).

X.-L. Qi, T. L. Hughes, and S.-C. Zhang, *Topological field theory of time-reversal invariant insulators*, Phys. Rev. B **78**, 195424 (2008).

K. Kobayashi, T. Ohtsuki, and K.-I. Imura, *Disordered Weak and Strong Topological Insulators*, Phys. Rev. Lett. **110**, 236803 (2013).

We note that a tight-binding model on an arbitrary lattice in $d$ dimension can be mapped onto the model on $\mathbb{Z}^d$ with suitably chosen hopping integrals.

J. Avron, R. Seiler, and B. Simon, *The Index of a Pair of Projections*, J. Func. Anal. **120**, 220 (1994).

Here, $A^{(\Lambda)}$ is expressed in terms of $P_{\rm F}$ and $\tilde{\mathcal{D}}_a({\bm x}) :=$ as $A^{(\Lambda)}=P_{\rm F} - {\tilde {\cal D}}^*_a P_{\rm F} {\tilde {\cal D}}_a$, where ${\mathbb I}_{4}$ is the 4-dimensional identity matrix and $\mathcal{D}_a$ should be thought of as a matrix of dimension $|\Lambda|$.

In this case, the operator $A^{(\Lambda)}$ is written as $A^{(\Lambda)} = \tilde{P}_{\rm F} - \tilde{D}_a \tilde{P}_{\rm F} \tilde{D}_a $, where $\tilde{P}_{\rm F} := P_{\rm F}  \otimes {\mathbb I}_{2}$, $\tilde{D}_a({\bm x}) := \sum_{i=1}^3$ , and $D_a^i({\bm x}) :=  (x_i - a_i)/|{\bm x} - {\bm a}|$. Here, ${\mathbb I}_{2}$ and ${\mathbb I}_{4}$ are the 2- and 4-dimensional identity matrices, respectively. Each operator $D_a^i$ should be thought of as a matrix of dimension $|\Lambda|$.

L. Fu, C. L. Kane, and E. J. Mele, *Topological Insulators in Three Dimensions*, Phys. Rev. Lett. **98**, 106803 (2007).

[^1]: E-mail address: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>