Abstract
The nonuniform \(\mathbb {Z}_{2}\) symmetric Kitaev chain, comprising alternating topological and normal regions, hosts localized states known as edge-zero modes (EZMs) at its interfaces. These EZMs can pair to form qubits that are resilient to quantum decoherence, a feature expected to extend to higher symmetric chains, i.e., parafermion chains. However, finite-size effects may impact this ideal picture. Diagnosing these effects requires first a thorough understanding of the low-energy physics where EZMs may emerge. Previous studies have largely focused on uniform chains, with nonuniform cases inferred from these results. While recent work [Narozhny, Sci. Rep. 7, 1447 (2017)] provides an insightful analytical solution for a nonuniform \(\mathbb {Z}_{2}\) chain with two topological regions separated by a normal one, its complexity limits its applicability to chains with more regions or higher symmetries. Here, we present a new approach based on decimating the highest-energy terms, facilitating the scalable analysis of \(\mathbb {Z}_{n}\) chains with any number of regions. We provide analytical results for both \(\mathbb {Z}_{2}\) and \(\mathbb {Z}_{3}\) chains, supported by numerical findings, and identify the critical lengths necessary to preserve well-separated EZMs.
Similar content being viewed by others
Introduction
Quantum computation has long grappled with the challenge of quantum decoherence1,2,3,4,5,6,7,8,9,10. Among various efforts, it has been demonstrated that topological systems such as Kitaev chain with \(\mathbb {Z}_{2}\) symmetry11,12,13,14,15,16 and its generalizations, i.e., parafermion chains with \(\mathbb {Z}_{n}\) (\(n\ge 3\)) symmetry17,18,19,20,21,22,23,24,25,26, can help tackle this issue (for convenience, we will henceforth refer to the Kitaev chain as the \(\mathbb {Z}_{2}\)-symmetric parafermion chain). These models exhibit distinct phases-gapped normal (N) and topological (T) separated by a critical point. Notably, within the T phase, certain decoherence-protected degrees of freedom manifest at the edges known as edge zero modes (EZMs). These EZMs, which commute with the Hamiltonian and display non-Abelian statistics27,28, are resilient to local perturbations, making them suitable for constructing robust qudits (for \(\mathbb {Z}_{2}\) symmetry, the term ’qubit’ is commonly used). However, the need for quantum computation goes beyond a single qudit. To reach this goal we may use many copies of these chains in T phase, e.g., \(n_T\), while making sure that the EZMs from the end of different regions are far enough that cannot hybridize and disappear. The common practice is linking \(n_T\) T-regions via \(n_N\) (\(= n_T - 1\)) N-regions (Fig. 1) which act as isolaters. This would result in \(2 \times n_T\) EZMs at the boundaries and a \(n^{n_T}\)-fold topological ground sate degeneracy.
(a) Schematic of a \(\mathbb {Z}_{2}\) symmetric chain featuring two topological regions (T1, T2) separated by a normal one (N1). Green and blue areas refer to T and N regions, respectively. For \(\mathbb {Z}_{3}\), substitute \(\gamma\)’s with \(\eta\)’s. (b) Four expected EZMs are highlighted in red.
We know that the above picture is idealized, as in reality, due to finite-size effects, there would be some splitting among the \(n^{n_T}\) states. If the splitting is significant, it may cause hybridization and dephasing, making the quantum computation less fault-tolerant. Analytical results studying the finite-size effects on the low-energy states of such systems, where EZMs emerge and play a key role in quantum computation, are rare, especially for \(\mathbb {Z}_{3}\) and higher-symmetry parafermion chains with \(N_T > 1\). This scarcity is largely due to the challenges of directly deriving the low-energy Hamiltonian of nonuniform chains. In Ref.29, an analytical study of the low-energy behavior for a \(\mathbb {Z}_{2}\) chain with \(n_T = 2\) was presented; however, the approach is complex and becomes impractically challenging for cases with \(n_T > 2\) or for chains with \(\mathbb {Z}_{n}\) symmetry. Thus, a new method that simplifies this process would be highly valuable for advancing studies in this area.
In this work, we show that by implementing the Decimation of the Highest-Energy term (DHE) method30,31,32,33, one not only recovers the result of29 but also extends it to cases with \(n_T \ge 3\). This method also enhances our understanding of parafermion chains with arbitrary \(n_T\), particularly for \(\mathbb {Z}_{3}\) symmetry, an area less explored compared to \(n_T = 1\). Our findings highlight the versatility of this approach across various symmetries. To validate our results using the DHE method, we employ numerical techniques such as Exact Diagonalization and Density Matrix Renormalization Group (DMRG)34,35,36,37,38,39. Finally, our analysis of the ground state energy splitting underscores the importance of region lengths in ensuring well-separated EZMs for optimal performance.
The paper is organized as follows: We first discuss the nonuniform \(\mathbb {Z}_{2}\) model, followed by an analysis of the \(\mathbb {Z}_{3}\) model. Next, we examine the ground-state energy splitting due to finite-size effects. Finally, we present our concluding remarks.
\(\mathbb {Z}_{2}\) Symmetry
The Hamiltonian of a \(\mathbb {Z}_{2}\) chain with \(n_T\) number of T regions reads
where
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.
Applying this method to Eq. (2) and neglecting constant terms yields
where
Labeling the \(\gamma ^{a}\)’s and \(\gamma ^{b}\)’s in Eq. (3) according to their positions at the interfaces (Fig. 1(b)) provides
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
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.
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
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
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.
\(\mathbb {Z}_{3}\) Symmetry
In this section, we aim to extend our findings to the \(\mathbb {Z}_{3}\) model, with certain adaptations. We begin by considering the Hamiltonian of the \(\mathbb {Z}_{3}\) chain, given by
where now
In this case, the couplings J and f are real and non-negative, and \(\alpha = e^{-i\phi }\) and \(\hat{\alpha } = e^{-i\hat{\phi }}\) (with \(\phi , \hat{\phi } \in [0, \pi /3)\)) represent the chirality factors17,18,19,20,25. The operators \(\eta _i\) are generalizations of Majorana fermions with the properties \(\eta ^3 = I\), \(\eta ^{\dagger } = \eta ^2\), and \(\eta _i \eta _{j>i} = \omega \eta _j \eta _i\), where \(\omega = e^{2\pi i / 3}\). The parameters \(n_T\), \(n_N = n_T - 1\), \(L_{Ti}\), \(L_{Ni}\), and L retain their previous meanings.
The DHE method is applicable to this model, as demonstrated for a single-region chain in the T phase in Ref.33. In Sec. II of Supplemental Material, we show its applicability to the N phase. By following procedures analogous to those used for the \(\mathbb {Z}_{2}\) symmetry, we obtain
where \(\eta\)’s refer to the EZMs residing at the interfaces, and
It is important to note that for certain angles, the \(\Delta\)’s do not diverge during the process. The effective chirality factors are now given by \(\beta _{Ti} = e^{-i\theta _{Ti}}\) and \(\beta _{Ni} = e^{-i\theta _{Ni}}\), where
When \(n_T = 1\) and, for example, \(\phi = \hat{\phi } = 0\), we find that \(H^{\text {eff}}_{\mathbb {Z}_{3}} \approx \Delta _{T1} (\bar{\omega } \eta ^{\dagger }_{1} \eta _{2} + h.c.)\), with \(\Delta _{T1} = f_{T1} \left( \frac{2 f_{T1}}{3J_{T1}} \right) ^{(L_{T1}-1)}\). This result is in agreement with Ref.20 (particularly Eq. 15, although note that due to different model definitions, the factor of 2 is absent there). It is also worth noting that, at \(\phi = \hat{\phi } = 0\), the EZMs display weak characteristics, whereas at \(\phi = \hat{\phi } = \pi /6\), the model exhibits super-integrability, leading to strong EZMs17,18,19,20.
Our investigations for \(n_T > 1\) revealed results similar to those for the \(\mathbb {Z}_{2}\) case. Specifically, we identified some localized configurations, the number of which follows Eq. (6), as well as some delocalized cases. However, our focus here remains on the localized configurations, as before. We will review two examples for \(n_T = 3\), which give rise to the configurations shown in Fig. 4a and b. We do not examine \(n_T = 2\) separately, as we expect the patterns to resemble those in Fig. 3 by setting \(L_{N2} = L_{T3} = 0\), effectively eliminating EZMs numbered 5 and 6. Thus, these examples are representative for both \(n_T = 2\) and \(n_T = 3\).
Here, we consider two sets of parameters yielding configurations similar to those shown in Fig. 4a and b, respectively. In the first set, we have \(L_{T1} = L_{T2} = L_{N2} = 6\), \(L_{N1} = L_{T3} = 4\), \(J_{Ti} = f_{Ni} = 1; \forall i\), \(f_{T1} = f_{T2} = J_{N2} = 0.1\), \(f_{T3} = 0.2\), and \(J_{N1} = 0.05\). We also assume \(\phi _{Ti} = \hat{\phi }_{Ti} = \pi /6; \forall i\) and \(\phi _{Ni} = \hat{\phi }_{Ni} = 0; \forall i\). In this setup, we have
(a) \(\Delta _{N2} = 8.78\times 10^{-9} \ll \Delta _{N1} = 6.17\times 10^{-8} \ll \Delta _{T1} = \Delta _{T2} = 4.87 \times 10^{-7} \ll \Delta _{T3} = 1.04\times 10^{-3}\).
Applying the DHE method and following the procedure used for \(\mathbb {Z}_{2}\), based on the competition among the \(\Delta\)’s, it is evident that \(\eta _5\) and \(\eta _6\) pair first, followed by \(\eta _3\) and \(\eta _4\), and finally \(\eta _1\) and \(\eta _2\). Thus,
where \(\epsilon ^{(a)}_{12} = \Delta _{T1},~\epsilon ^{(a)}_{34} = \Delta _{T2}\), and \(\epsilon ^{(a)}_{56} = \Delta _{T3}\).
For the second set, we maintained the same parameters as in the example above but set \(J_{N1} = 0.3\) and applied the DHE method. We obtain \(\Delta _{N2} = 8.78 \times 10^{-9}\). Here, we observe that
(b) \(\Delta _{N2} \ll \Delta _{T1} = \Delta _{T2} = 4.87 \times 10^{-7} \ll \Delta _{N1} = 4.80\times 10^{-4} < \Delta _{T3} = 1.04\times 10^{-3}\)
and thus
where \(\epsilon ^{(b)}_{14} = \frac{2}{3} \frac{\Delta _{T1}\Delta _{T2}}{\Delta _{N1}} = 3.30\times 10^{-10}\), \(\epsilon ^{(b)}_{23} = \Delta _{N1}\), and \(\epsilon ^{(b)}_{56} = \Delta _{T3}\).
Unfortunately, standard methods for diagonalizing \(H_{\mathbb {Z}_{2}}\) cannot be applied to \(H_{\mathbb {Z}_{3}}\)17,40. Thus, alternative approaches are needed to numerically identify coupled EZMs. In Sec. III of Supplemental Material, we demonstrate that the Green’s function \(G_{k,n} = |\langle {\eta ^{a\dagger }_{k} \eta ^{b}_n} \rangle |\) effectively serves this purpose. Specifically, long-range correlations reveal which EZMs are coupled. Numerical results for the examples in Fig. 6, obtained using DMRG (for details Sec. IV of Supplemental Material), align well with our expectations. Furthermore, the plots in Fig. 7, obtained using DMRG, confirm nearly 27-fold degenerate ground states with a finite gap to excited states in the aforementioned examples. These findings suggest that the conclusions drawn from \(\mathbb {Z}_{2}\) symmetry can be extended to the \(\mathbb {Z}_{3}\) case.
Heatmap illustrating correlation functions in a \(\mathbb {Z}_{3}\) parafermion chain with \(n_T = 3\), obtained using the DMRG method. Red circles highlight long-range correlations among coupled EZMs. Panel (a) corresponds to a parameter set where \(\Delta _{N2} \ll \Delta _{N1} \ll \Delta _{T1} = \Delta _{T2} \ll \Delta _{T3}\), leading to the configuration shown in Fig. 4a. Panel (b) corresponds to a parameter set where \(\Delta _{N2} \ll \Delta _{T1} = \Delta _{T2} \ll \Delta _{N1} < \Delta _{T3}\), leading to the configuration shown in Fig. 4b. The specific parameters are \(L_{T1} = L_{T2} = L_{N2} = 6\), \(L_{N1} = L_{T3} = 4\), \(J_{Ti} = f_{Ni} = 1\), \(f_{T1} = f_{T2} = J_{N2} = 0.1\), \(f_{T3} = 0.2\), and \(J_{N1} = 0.05\) for panel (a), while \(J_{N1} = 0.3\) for panel (b). The chirality factors were set to \(\phi _{Ti} = \hat{\phi }_{Ti} = \pi /6\) and \(\phi _{Ni} = \hat{\phi }_{Ni} = 0\). As in the \(\mathbb {Z}_{2}\) examples, tuning the system parameters alters the coupling patterns of the EZMs. More details are provided in the main text.
Panels (a) and (b) show the energy spectrum, highlighting a clear gap between the many-body ground state manifold and the first excited states for the parameter sets (a) and (b) in the \(\mathbb {Z}_3\) examples, respectively. These results were obtained using the DMRG method. Panel (a) corresponds to the parameter set where \(\Delta _{N2} \ll \Delta _{N1} \ll \Delta _{T1} = \Delta _{T2} \ll \Delta _{T3}\), while panel (b) corresponds to the parameter set where \(\Delta _{N2} \ll \Delta _{T1} = \Delta _{T2} \ll \Delta _{N1} < \Delta _{T3}\). The presence of a finite energy gap confirms the stability of the ground state manifold, which is essential for maintaining well-separated EZMs. This stability is crucial for ensuring the robustness of qudits in topological quantum computation.
DMRG results for the GES in a \(\mathbb {Z}2\) chain with \(n_T = 2\) are shown in panel (a), where varying \(L_{N1}\) while keeping other parameters fixed leads to \(\text {GES} \propto m^{L_{N1}}\), analogous to \(\Delta _{N1}\). Panel (b) further confirms that \(m = J_{N1}/f_{N1}\). Similar analyses for a \(\mathbb {Z}_{3}\) chain validate the relation \(\text {GES} \propto m^L\). Panels (c) and (d) depict the dependence of m on chirality in the T and N phases, respectively. The dashed lines represent theoretical predictions from Eq. (12), while the data points (circles, stars, and squares) correspond to numerical DMRG results. Chirality has a significant effect on m at small J, f ratios but becomes negligible as the ratio increases, leaving m primarily determined by J and f. These findings underscore the crucial role of critical lengths in maintaining well-separated EZMs.
Ground state energy splitting analysis
The concept of critical lengths is essential for maintaining well-separated EZMs and ensuring the robustness of the topological phase against finite-size effects. If the system size is smaller than the critical length, interactions between EZMs lift degeneracies in the ground-state manifold (and excited states), resulting in ground-state energy splitting (GES) and EZM hybridization, thereby compromising their topological protection. This has direct implications for quantum computation, where well-separated zero modes are crucial for preserving logical qubits against decoherence. Ensuring that the system size exceeds the critical length significantly enhances stability. Here, we demonstrate that the values of \(\Delta\) reflect the GES induced by finite-size effects in each region, allowing us to define critical lengths accordingly.
In \(\mathbb {Z}_{2}\) chains, when \(n_T = 1\), we confirm the Kitaev result11 from Eq. (5), affirming the link between \(\Delta _{T1}\) and GES. It’s crucial to ensure that \(L_{T1} \gg \xi _{T1}\) to maintain well-separated EZMs, defining the critical length as \(\xi _{T1} = 1/|\ln (f_{T1}/J_{T1})|\). For \(n_T = 2\), our numerical results from DMRG, depicted in Fig. 8(a, b), show that GES \(\propto (\frac{J_{N1}}{f_{N1}})^{L_{N1}}\), while ensuring \(L_{Tn} \gg \xi _{Tn}\) for \(n=1,2\), with all other parameters fixed. This procedure can generally be repeated for other regions, yielding similar outcomes. Consequently, critical lengths can be defined for each region as:
where \(q = \frac{f}{J}\) (\(q = \frac{J}{f}\)) for T (N) regions. Thus, to guarantee well-separated EZMs, it’s necessary to have \(L_{Ti} \gg \xi _{Ti}\) and \(L_{Ni} \gg \xi _{Ni}\) for all i.
In \(\mathbb {Z}_{3}\) chains, the GES analysis obtained using DMRG, depicted in Fig. 8c,d, mirrors the findings in the \(\mathbb {Z}_{2}\) case, exhibiting an anticipated trend for small ratios of J and f, albeit with deviations for larger ones. Once more, ensuring well-separated EZMs necessitates \(L_{Ti} \gg \xi _{Ti}\) and \(L_{Ni} \gg \xi _{Ni}\) for all i, utilizing Eq. (16) with \(q = \frac{2\cos \phi }{(2\cos \phi )^2-1} \times \frac{f}{J}\) (\(q = \frac{2\cos \hat{\phi }}{(2\cos \hat{\phi })^2-1} \times \frac{J}{f}\)) for T (N) regions.
As easily verified, the condition \(L \gg \xi\) is met for all regions in the examples studied throughout this paper.
Conclusion
In this study, we introduced a novel approach based on the decimation of the highest energy term to derive the low-energy Hamiltonian for \(\mathbb {Z}_{2}\) and \(\mathbb {Z}_{3}\) symmetric parafermion chains, encompassing a range of topological and normal regions. By applying the DHE method, we confirmed previously established results for the \(n_T = 1, 2\) \(\mathbb {Z}_{2}\) chains and \(n_T = 1\) \(\mathbb {Z}_{3}\) chains. We extended our analysis to the \(n_T = 3\) case for both \(\mathbb {Z}_{2}\) and \(\mathbb {Z}_{3}\) symmetries, validating the analytical results with numerical computations. Furthermore, our investigation of the ground state energy splitting provided a crucial metric for ensuring well-separated EZMs in finite-size systems, which was vital for their stability. Overall, the method we presented proved highly versatile and can be readily adapted to investigate higher symmetries, paving the way for further exploration in the field.
Data Availability
All data generated or analyzed during this study are included in this published article [and its supplementary information files].
References
Nielsen, M. A. & Chuang, I. L. Quantum Computation and Quantum Information (Cambridge University Press, 2010).
Rieffel, E. G. & Polak, W. H. Quantum Computing: A Gentle Introduction (MIT press, 2011).
DiVincenzo, D. P. & Loss, D. Quantum computers and quantum coherence. J. Magn. Magn. Mater. 200, 202–218 (1999).
Knill, E. & Laflamme, R. Theory of quantum error-correcting codes. Phys. Rev. A 55, 900 (1997).
Lidar, D. A., Chuang, I. L. & Whaley, K. B. Decoherence-free subspaces for quantum computation. Phys. Rev. Lett. 81, 2594 (1998).
Zurek, W. H. Decoherence and the transition from quantum to classical. Phys. Today 44, 36–44 (1991).
Schlosshauer, M. Quantum decoherence. Phys. Rep. 831, 1–57 (2019).
Chuang, I. L., Laflamme, R., Shor, P. W. & Zurek, W. H. Quantum computers, factoring, and decoherence. Science 270, 1633–1635 (1995).
Shor, P. W. Fault-tolerant quantum computation. In Proceedings of 37th Conference on Foundations of Computer Science (ed. Shor, P. W.) 56–65 (IEEE, 1996).
Shor, P. W. Scheme for reducing decoherence in quantum computer memory. Phys. Rev. A 52, R2493 (1995).
Kitaev, A. Y. Unpaired majorana fermions in quantum wires. Phys. Usp. 44, 131 (2001).
Leumer, N., Marganska, M., Muralidharan, B. & Grifoni, M. Exact eigenvectors and eigenvalues of the finite kitaev chain and its topological properties. J. Phys. Condens. Matter 32, 445502 (2020).
Rančić, M. J. Exactly solving the kitaev chain and generating majorana-zero-modes out of noisy qubits. Sci. Rep. 12, 19882 (2022).
Sung, K. J., Rančić, M. J., Lanes, O. T. & Bronn, N. T. Simulating majorana zero modes on a noisy quantum processor. Quantum Sci. Technol. 8, 025010 (2023).
Miao, J.-J., Jin, H.-K., Zhang, F.-C. & Zhou, Y. Majorana zero modes and long range edge correlation in interacting kitaev chains: Analytic solutions and density-matrix-renormalization-group study. Sci. Rep. 8, 488 (2018).
Burnell, F., Shnirman, A. & Oreg, Y. Measuring fermion parity correlations and relaxation rates in one-dimensional topological superconducting wires. Phys. Rev. B 88, 224507 (2013).
Fendley, P. Parafermionic edge zero modes in zn-invariant spin chains. J. Stat. Mech. Theory Exp. 2012, P11020 (2012).
Fendley, P. Free parafermions. J. Phys. A Math. Theor. 47, 075001 (2014).
Alicea, J. & Fendley, P. Topological phases with parafermions: Theory and blueprints. Annu. Rev. Condens. Matter Phys. 7, 119–139 (2016).
Jermyn, A. S., Mong, R. S., Alicea, J. & Fendley, P. Stability of zero modes in parafermion chains. Phys. Rev. B 90, 165106 (2014).
Iemini, F., Mora, C. & Mazza, L. Topological phases of parafermions: A model with exactly solvable ground states. Phys. Rev. Lett. 118, 170402 (2017).
Vaezi, A. Fractional topological superconductor with fractionalized majorana fermions. Phys. Rev. B 87, 035132 (2013).
Vaezi, M.-S. & Vaezi, A. Numerical observation of parafermion zero modes and their stability in 2d topological states. Preprint at arXiv:1706.01192 (2017).
Mahyaeh, I. Study of the Phase Diagram of Zn Symmetric Chains (Stockholm University, 2020).
Zhuang, Y., Changlani, H. J., Tubman, N. M. & Hughes, T. L. Phase diagram of the z 3 parafermionic chain with chiral interactions. Phys. Rev. B 92, 035154 (2015).
Wu, L.-A. & Lidar, D. Qubits as parafermions. J. Math. Phys. 43, 4506–4525 (2002).
Nayak, C., Simon, S. H., Stern, A., Freedman, M. & Das Sarma, S. Non-abelian anyons and topological quantum computation. Rev. Mod. Phys. 80, 1083–1159 (2008).
Rao, S. Introduction to abelian and non-abelian anyons. Topology and condensed matter physics 399–437 (2017).
Narozhny, B. Majorana fermions in the nonuniform ising-kitaev chain: exact solution. Sci. Rep. 7, 1447 (2017).
Dasgupta, C. & Ma, S.-K. Low-temperature properties of the random heisenberg antiferromagnetic chain. Phys. Rev. B 22, 1305 (1980).
Fisher, D. S. Critical behavior of random transverse-field ising spin chains. Phys. Rev. B 51, 6411 (1995).
Shivamoggi, V. B. Majorana Fermions and Dirac Edge States in Topological Phases (University of California, 2011).
Benhemou, A., Angkhanawin, T., Adams, C. S., Browne, D. E. & Pachos, J. K. Universality of z 3 parafermions via edge-mode interaction and quantum simulation of topological space evolution with rydberg atoms. Phys. Rev. Res. 5, 023076 (2023).
White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 69, 2863 (1992).
White, S. R. Density-matrix algorithms for quantum renormalization groups. Phys. Rev. B 48, 10345 (1993).
Schollwöck, U. The density-matrix renormalization group. Rev. Mod. Phys. 77, 259 (2005).
Schollwöck, U. The density-matrix renormalization group in the age of matrix product states. Ann. Phys. 326, 96–192 (2011).
Verstraete, F. et al. Density matrix renormalization group, 30 years on. Nat. Rev. Phys. 5, 273–276 (2023).
Vaezi, M.-S. & Vaezi, A. Entanglement distance between quantum states and its implications for a density-matrix renormalization group study of degenerate ground states. Phys. Rev. B 96, 165129 (2017).
Cobanera, E. & Ortiz, G. Fock parafermions and self-dual representations of the braid group. Phys. Rev. A 89, 012328 (2014).
Acknowledgements
We are grateful to Abolhassan Vaezi and Reza Asgari for their valuable and insightful discussions.
Author information
Authors and Affiliations
Contributions
M.-S.V. conceived the project, devised the methodology, conducted the investigation, analyzed the data, generated the visualizations, and drafted the original manuscript. M.M.N. contributed to the investigation, visualizations, and assisted in revising the draft. M.-S.V. supervised the project.
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Nasiri Fatmehsari, M.M., Vaezi, MS. Low-energy physics and finite-size effects in nonuniform parafermion chains. Sci Rep 15, 11446 (2025). https://doi.org/10.1038/s41598-025-94657-z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-025-94657-z