(1)
where
$$\begin{aligned} H_{Ti}= & i J_{Ti}\sum ^{{\varvec{\textrm{r}}}_{Ti} + L_{Ti}-1}_{j = {\varvec{\textrm{r}}}_{Ti} +1} {\gamma ^{b}_{j} \gamma ^{a}_{j+1}} + if_{Ti} \sum ^{{\varvec{\textrm{r}}}_{Ti} + L_{Ti}}_{j = {\varvec{\textrm{r}}}_{Ti}+1} {\gamma ^{a}_{j} \gamma ^{b}_{j}},\nonumber \\ H_{Ni}= & i J_{Ni}\sum ^{{\varvec{\textrm{r}}}_{Ni} + L_{Ni}}_{j = {\varvec{\textrm{r}}}_{Ni}} {\gamma ^{b}_{j} \gamma ^{a}_{j+1}} + if_{Ni} \sum ^{{\varvec{\textrm{r}}}_{Ni} + L_{Ni}}_{j = {\varvec{\textrm{r}}}_{Ni}+1} {\gamma ^{a}_{j} \gamma ^{b}_{j}}. \end{aligned}$$
(2)
Here, \(\gamma\)’s are the Majorana fermion operators satisfying the relations \(\gamma = \gamma ^{\dag }\), and \(\{\gamma ^{\alpha }_i,\gamma ^{\beta }_j\} = 2 \delta _{ij}\delta _{\alpha \beta }\). The non-negative couplings of the n-th T (N) region with \(L_{Ti}\) (\(L_{Ni}\)) number of sites satisfy the relation \(f_{Ti}/J_{Ti} < 1\) (\(f_{Ni}/J_{Ni} > 1\)). For convenience, we have absorbed the interaction terms between T and N regions into the N regions. The vectors \({\varvec{\textrm{r}}}_{Ti}\) and \({\varvec{\textrm{r}}}_{Ni}\) also refer to real space starting point of i-th T and N regions (see Fig. 1a), respectively.
As stated earlier we aim to study low-energy physics where EZMs appear and hold particular significance for quantum computation by employing the DHE method. This method is an iterative process that involves eliminating the most energetic term within the Hamiltonian and replacing it with effective longer-range interactions using second-order perturbation theory at each iteration. For instance, within the T region where \(J > f\), decimating J-let’s say between \(\gamma ^{b}_{1}\) and \(\gamma ^{a}_{2}\) -leads to an effective interaction \(f_{\text {eff}} = \frac{f^2}{J}\) between \(\tilde{\gamma }^{a}_{1}\) and \(\tilde{\gamma }^{b}_{2}\), which can be approximated by \(\gamma ^{a}_{1}\) and \(\gamma ^{b}_{2}\), respectively (Fig. 2). Smaller f/J ratios yield better approximations. Conversely, within the N region where \(f > J\), one should decimate f terms, resulting in an effective interaction \(J_{\text {eff}} = \frac{J^2}{f}\). The details of this method applied to a single T and N regions are available in Ref.32.
Schematic illustration of the DHE method. Decimating the J couplings \(\gamma ^{b}_{1}\) and \(\gamma ^{a}_{2}\) in panel (a) leads to the configuration of (b).
Applying this method to Eq. (2) and neglecting constant terms yields
$$\begin{aligned} H_{\mathbb {Z}_{2}, Ti}\approx & i \Delta _{Ti} \gamma ^{a}_{{\varvec{\textrm{r}}}_{Ti} +1} \gamma ^{b}_{{\varvec{\textrm{r}}}_{Ti}+L_{Ti}}, \nonumber \\ H_{\mathbb {Z}_{2}, Ni}\approx & i \Delta _{Ni} \gamma ^{b}_{{\varvec{\textrm{r}}}_{Ni}} \gamma ^{a}_{{\varvec{\textrm{r}}}_{Ni}+L_{Ni}+1}, \end{aligned}$$
(3)
where
$$\begin{aligned} \Delta _{Ti} = f_{Ti}(\frac{f_{Ti}}{J_{Ti}})^{L_{Ti-1}},~ \Delta _{Ni} = J_{Ni}(\frac{J_{Ni}}{f_{Ni}})^{L_{Ni}}. \end{aligned}$$
(4)
Labeling the \(\gamma ^{a}\)’s and \(\gamma ^{b}\)’s in Eq. (3) according to their positions at the interfaces (Fig. 1(b)) provides
$$\begin{aligned} H^{\text {eff}}_{\mathbb {Z}_{2}} \approx i \sum ^{n_T}_{i=1} \Delta _{Ti} \gamma _{2i-1} \gamma _{2i} + i \sum ^{n_N}_{i=1} \Delta _{Ni} \gamma _{2i} \gamma _{2i+1}. \end{aligned}$$
(5)
Interestingly, this outcome resembles the famous Kitaev chain. One can easily confirm that for \(n_T = 1\), the result for a single-region Kitaev chain11 is retrieved, i.e., \(H^{\text {eff}}_{\mathbb {Z}_{2}} = i \epsilon _{1} \gamma _{1} \gamma _{2}\), where \(\epsilon _{1} = \Delta _{T1}\) is the ground state energy splitting. The associated non-local quasi-fermion operators then reads \(c_{12} \approx (\gamma _{1} + i \gamma _{2})/\sqrt{2}\).
Unfortunately, for \(n_T > 1\), the situation is not as trivial as in the single-region case. This becomes evident from Eq. (5), where the competition among \(\Delta\)’s plays a major role in determining which EZMs couple to form a quasi-fermion and which energy levels correspond to the lowest energies. However, as we will demonstrate, applying the DHE method to Eq. (5) can still be beneficial. To better understand this, let us first review the case \(n_T = 2\), for which an analytical result is available29. Here, different scenarios are possible. Consider, for example, two extreme cases:
- (a)
$$\Delta _{T1}, \Delta _{T2} \gg \Delta _{N1}$$
.
- (b)$$\Delta _{T1}, \Delta _{T2} \ll \Delta _{N1}.$$
In case (a), if we further assume that \(\Delta _{T1} < \Delta _{T2}\), then using the DHE method, we conclude that \(\gamma _3\) couples with \(\gamma _4\) with a coupling \(\epsilon _{34} = \Delta _{T2}\). By eliminating this pair (as part of the DHE procedure), we are left with \(\gamma _1\) and \(\gamma _2\), which are now compelled to couple with \(\epsilon _{12} = \Delta _{T1}\). That is, \(H^{(a)}_{\mathbb {Z}_{2},n_T=2} = i (\epsilon _{12} \gamma _{1} \gamma _{2} + \epsilon _{34} \gamma _{3} \gamma _{4})\). The quasi-fermions are \(c_{12} = (\gamma _1 + i \gamma _2)/\sqrt{2}\) and \(c_{34} = (\gamma _3 + i \gamma _4)/\sqrt{2}\), with the corresponding low-energy states (LES) given by \(|n_{12}, n_{34} \rangle\). Evidently, \(\epsilon _{12} < \epsilon _{34}\). Note that if we had assumed \(\Delta _{T1} > \Delta _{T2}\), then we would have instead \(\epsilon _{12} > \epsilon _{34}\). Case (b) is more interesting. Here, \(\gamma _2\) and \(\gamma _3\) are coupled first with \(\epsilon _{23} = \Delta _{N1}\). Continuing the DHE procedure and eliminating this pair, we are left with \(\gamma _1\) and \(\gamma _4\), which now interact with an effective strength \(\Delta _{\text {eff}} = \frac{\Delta _{T1} \Delta _{T2}}{\Delta _{N1}}\). Therefore, they couple with \(\epsilon _{14} = \Delta _{\text {eff}}\), leading to \(H^{(b)}_{\mathbb {Z}_{2},n_T=2} = i (\epsilon _{14} \gamma _{1} \gamma _{4} + \epsilon _{23} \gamma _{2} \gamma _{3})\). The quasi-fermions are \(c_{14} = (\gamma _1 + i \gamma _4)/\sqrt{2}\) and \(c_{23} = (\gamma _2 + i \gamma _3)/\sqrt{2}\), with the corresponding LES given by \(|n_{14}, n_{23} \rangle\). These results are in exact agreement with Ref.29.
In cases (a) and (b), the EZMs are localized at the interfaces, according to Fig. 3a and b, respectively. However, one might ask what happens if the \(\Delta\)’s are comparable. Let us assume, for the moment, that \(\Delta _{T1} \simeq \Delta _{T2} \simeq \Delta _{N1}\). Unfortunately, in such cases, the DHE method does not work well, as its main assumption is that some couplings are significantly smaller or larger than the rest (or at least than the neighboring couplings at each step). As shown in Ref.29, when the \(\Delta\) values are comparable, the zero modes (ZMs) become ”delocalized”, meaning they are spread across multiple sites. For instance, one ZM may appear as a weighted combination of \(\gamma _1\) and \(\gamma _3\). However, the primary goal of topological quantum computation is to leverage the localized nature of EZMs, which ensures robustness against errors. Delocalization, on the other hand, compromises this robustness. Therefore, in this study, we focus solely on cases with localized ZMs. We observed that the number of localized configurations, denoted by \(n_C\), for \(n_T = 2\) is two. Our investigations further revealed that, in general, this pattern follows the Catalan numbers
$$\begin{aligned} n_C(n_T) = \frac{(2n_T)!}{(n_T + 1)! n_T!}. \end{aligned}$$
(6)
The central insight leading to this result is that the lines connecting different pairs of EZMs must not intersect. For instance, with \(n_T = 3\), we identified five possible localized configurations, as illustrated in Fig. 4.
The panels (a) and (b) show two ways of linking EZMs in a \(\mathbb {Z}_{2}\) or \(\mathbb {Z}_{3}\) chain with \(n_T = 2\) to form localized configurations.
The panels (a–e) show various ways of linking EZMs in a \(\mathbb {Z}_{2}\) or \(\mathbb {Z}_{3}\) chain with \(n_T = 3\) to form localized configurations.
Now, let us investigate our method for \(n_T = 3\) and test it against numerical results. We will go over two examples that illustrate the method for the rest. Let us begin with the first set of parameters leading to the configuration shown in Fig. 4(d): \(L_{T_i} = L_{N_i} = 10; \forall i\), \(J_{T_i} = J_{N_i} = 1; \forall i\), \(f_{T1} = f_{T2} = f_{T3} = 0.1\), and \(f_{N1} = 6.75\), \(f_{N2} = 5\). Trivially, in this case,
(d) \(\Delta _{T1} = \Delta _{T2} = \Delta _{T3} = 1\times 10^{-10} \ll \Delta _{N1} = 5.09\times 10^{-9} \ll \Delta _{N2} = 1.02\times 10^{-7}\).
The approach involves utilizing the DHE method once again as follows. Given that \(\Delta _{N2}\) significantly surpasses the other couplings, we infer that \(\gamma _4\) and \(\gamma _5\) primarily engage with a strength of \(\epsilon ^{(d)}_{45} = \Delta _{N2} = 1.02 \times 10^{-7}\). Continuing the process by decimating this pair, we find that \(\gamma _3\) and \(\gamma _6\) interact with an effective strength of \(\Delta _{\text {eff}} = \frac{\Delta _{T2} \Delta _{T3}}{\Delta _{N2}} = 0.98 \times 10^{-13}\). Upon comparing the remaining couplings (\(\Delta _{\text {eff}} \ll \Delta _{T1} \ll \Delta _{N1}\)) and identifying the maximum, we observe that \(\gamma _2\) and \(\gamma _3\) interact with a strength of \(\epsilon ^{(d)}_{23} = \Delta _{N1} = 5.09 \times 10^{-9}\). Subsequently, after decimating this pair, we are left with \(\gamma _1\) and \(\gamma _6\), which are compelled to interact with a strength of \(\epsilon ^{(d)}_{16} = \frac{\Delta {T1} \Delta _{\text {eff}}}{\epsilon ^{(d)}_{23}} = \frac{\Delta _{T1} \Delta _{T2} \Delta _{T3}}{\Delta _{N1} \Delta _{N2}} = 1.91 \times 10^{-15}\). Thus, we obtain
$$\begin{aligned} H^{(d)}_{\mathbb {Z}_{2},n_T=3} = i \epsilon ^{(d)}_{16} \gamma _{1}\gamma _{6} + i \epsilon ^{(d)}_{23} \gamma _{2}\gamma _{3} + i \epsilon ^{(d)}_{45} \gamma _{4}\gamma _{5}, \end{aligned}$$
(7)
leading to the association of the state \(| {\psi ^{(d)}} \rangle = | {n_{16},n_{23},n_{45}} \rangle\). The numerical results obtained using the Exact Diagonalization method (see Sec. I of Supplemental Material for details) show that \(\epsilon ^{(d)}_{16} = 2.13 \times 10^{-15} \ll \epsilon ^{(d)}_{23} = 4.93 \times 10^{-9} \ll \epsilon ^{(d)}_{45} = 0.97 \times 10^{-7}\), in excellent agreement with our findings. In Fig. 5(a), we further present the numerical results for the three lowest single-body wave functions, each exhibiting the anticipated behavior.
In panels (a) and (b), the green stars, red squares, and blue circles represent the single-body wave functions of the first, second, and third excited states, respectively. The parameters used for these systems, with \(L_{Ti} = L_{Ni} = 10\) and \(J_{Ti} = J_{Ni} = 1\), are: for panel (a), \(f_{T1} = f_{T2} = f_{T3} = 0.1\), \(f_{N1} = 6.75\), and \(f_{N2} = 5\); and for panel (b), \(f_{T1} = f_{T3} = 0.05\), \(f_{T2} = 0.3\), and \(f_{N1} = f_{N2} = 6.75\). The parameters in panels (a) and (b) are chosen such that the EZMs couple (and form a quasi-fermion) according to configurations (d) and (e) of Fig. 4, respectively. As evident, tuning the system parameters alters the coupling patterns of the EZMs.
The second set of parameters, which leads to the configuration shown in Fig. 4e, was obtained by keeping all previous parameters unchanged while setting \(f_{T1} = f_{T3} = 0.05\), \(f_{T2} = 0.3\), and \(f_{N1} = f_{N2} = 6.75\). We find that
(e) \(\Delta _{T1} = \Delta _{T3} = 9.77\times 10^{-14} \ll \Delta _{N1} = \Delta _{N2} = 5.09\times 10^{-9} \ll \Delta _{T2} = 5.90\times 10^{-6}\).
By applying the DHE technique as before, it is evident that \(\gamma _3\) and \(\gamma _4\) primarily interact with a magnitude of \(\epsilon ^{(e)}_{34} = \Delta _{T2} = 5.90 \times 10^{-6}\). Upon decimating this pair, we observe \(\gamma _2\) and \(\gamma _5\) interacting with \(\Delta _{\text {eff}} = \frac{\Delta _{N1} \Delta _{N2}}{\Delta _{T2}} = 4.39\times 10^{-12}\). Given that \(\Delta _{T1},~\Delta _{T3} \ll \Delta _{\text {eff}}\), it follows that \(\gamma _2\) and \(\gamma _5\) should pair with strength \(\epsilon ^{(e)}_{25} = \Delta _{\text {eff}} = 4.39\times 10^{-12}\). Ultimately, \(\gamma _{1}\) and \(\gamma _{6}\) interact with \(\epsilon ^{(e)}_{16} = \frac{\Delta _{T1} \Delta _{T3}}{\epsilon ^{(e)}_{25}} = \frac{\Delta _{T1} \Delta _{T2} \Delta _{T3}}{\Delta _{N1}\Delta _{N2}} = 2.17\times 10^{-15}\). This confirms
$$\begin{aligned} H^{(e)}_{\mathbb {Z}_{2},n_T=3} = i \epsilon ^{(e)}_{16} \gamma _{1}\gamma _{6} + i \epsilon ^{(e)}_{25} \gamma _{2}\gamma _{5} + i \epsilon ^{(e)}_{34} \gamma _{3}\gamma _{4}, \end{aligned}$$
(8)
resulting in the state \(| {\psi ^{(e)}} \rangle = | {n_{16},n_{25},n_{34}} \rangle\). The numerical results obtained using the Exact Diagonalization method indicate that \(\epsilon ^{(e)}_{16} = 2.24\times 10^{-15} \ll \epsilon ^{(e)}_{25} = 4.29\times 10^{-12} \ll \epsilon ^{(e)}_{34} = 5.27\times 10^{-6}\) which aligns very well with our discoveries. Fig. 5b further shows the numerical results for the three lowest single-body wave functions, each displaying the expected characteristics.
It is worth noting that the above configurations can arise from different parameter sets. For instance, in configuration (d), choosing \(f_{N1} = 5\) and \(f_{N2} = 6.75\) yields the same structural pattern. This flexibility extends beyond this example and holds across various cases (not explicitly mentioned in this manuscript), simplifying parameter selection in certain situations.
Now, we try to highlight key aspects relevant to broader configurations. For \(n_T > 3\), the process remains similar to smaller \(n_T\) cases: we systematically identify the largest \(\Delta\) and follow the steps outlined in the DHE method. However, questions may arise regarding the scalability of our approach, e.g., the number of localized configurations. Using Eq. 6 and applying Stirling’s approximation, we find that the number of localized configurations grows as \(n_C(n_T) \sim 4^{n_T}/n_T^{3/2}\). While this exhibits exponential growth, the polynomial factor \(n_T^{3/2}\) moderates the scaling, making it slightly slower than pure exponential growth. Thus, engineering system parameters to target specific configurations becomes increasingly difficult as \(n_T\) increases and requires careful tuning.
Nevertheless, the DHE method remains applicable for determining which configuration emerges, provided two key conditions are met: (1) the lengths \(L_{Ti}\) and \(L_{Ni}\) are sufficiently large to suppress unintended hybridization between EZMs, and (2) a clear hierarchy exists among the couplings \({\Delta _{Ti}, \Delta _{Ni}}\) to ensure unambiguous iterative decimation steps. When these conditions hold, the method scales polynomially with \(n_T\), as each decimation step reduces the problem size by a constant factor. The polynomial scaling arises from the recursive nature of the DHE method, where the number of steps grows linearly with \(n_T\), and each step involves simple algebraic updates to effective couplings. Consequently, the method can accommodate chains with arbitrary \(n_T\), as the iterative process applies to each successive EZM pair. However, if the neighboring couplings become comparable at any step, the perturbative approach breaks down, leading to delocalized zero modes. Thus, the method remains effective for \(n_T > 3\) as long as a clear coupling hierarchy is maintained, ensuring EZM localization and the validity of the DHE method.
Before analyzing the \(\mathbb {Z}_{3}\) model, a critical question remains: Do these different localized configurations impact topological quantum computation? Addressing this requires rigorous investigation, as identifying configurations that optimize qubit stability is crucial for practical applications. This will be a primary focus of our future work, which aims to determine the configurations that best preserve EZM robustness. Our present study establishes a solid foundation for these investigations by developing a systematic framework for analyzing the low-energy physics of nonuniform chains. To the best of our knowledge, no prior work has provided such a detailed and applicable derivation, which we believe will be essential for advancing the understanding of finite-size effects and the configuration landscape in parafermion systems.
The principles outlined here extend naturally to the \(\mathbb {Z}_{3}\) case, demonstrating the versatility of our approach.