Topological phase diagrams of exactly solvable non-Hermitian interacting Kitaev chains

TL;DR

This paper presents exact and numerical phase diagrams for non-Hermitian interacting Kitaev chains, revealing how non-Hermiticity affects topological phases and degeneracies in complex quantum systems.

Contribution

It provides the first comprehensive phase diagrams for non-Hermitian interacting Kitaev chains, including exact phase boundaries and exploration beyond solvable regimes.

Findings

Hermitian phases can disappear with increased non-Hermitian effects

Non-Hermitian topological degeneracy can persist beyond solvable regimes

The work offers a detailed characterization of non-Hermitian topology in interacting systems

Abstract

Many-body interactions give rise to the appearance of exotic phases in Hermitian physics. Despite their importance, many-body effects remain an open problem in non-Hermitian physics due to the complexity of treating many-body interactions. Here, we present a family of exact and numerical phase diagrams for non-Hermitian interacting Kitaev chains. In particular, we establish the exact phase boundaries for the dimerized Kitaev-Hubbard chain with complex-valued Hubbard interactions. Our results reveal that some of the Hermitian phases disappear as non-Hermiticty is enhanced. Based on our analytical findings, we explore the regime of the model that goes beyond the solvable regime, revealing regimes where non-Hermitian topological degeneracy remains. The combination of our exact and numerical phase diagrams provides an extensive description of a family of non-Hermitian interacting models. Our results provide a stepping stone toward characterizing non-Hermitian topology in realistic interacting quantum many-body systems.

Equations

$$\mathcal{H}=-\sum _{j=1}^{L-1}\left[t_j\left(c_j^{†}{c}_{j+1}+c_{j+1}^{†}{c}_j\right)+{\Delta }_j\left(c_j^{†}{c}_{j+1}^{†}+c_{j+1}{c}_j\right)\right]$$

$$+\sum _{j=1}^{L-1}(U_j-\mathrm{i}{\delta }_j)\left(2n_j-1\right)\left(2n_{j+1}-1\right)-\mu \sum _{j=1}^L(n_j-\frac{1}{2}),$$

$${\mathcal{O}}_j$$

$$\mathcal{H}$$

$$+\sum _j{\tilde{U}}_j\left[f_{j+1}^{†}{f}_j+f_j^{†}{f}_{j+1}-f_j^{†}{f}_{j+1}^{†}-f_{j+1}{f}_j\right],$$

$$\mathcal{H}$$

$$-\sum _k(t+\tilde{U})\left[w_k{f}_{-kA}^{†}{f}_{kB}^{†}+w_k^{\ast }{f}_{kB}{f}_{-kA}\right],$$

$$\frac{{\Lambda }_k^2}{4}$$

$$\frac{U}{t}=\sqrt{\frac{{\delta }^2}{t^2}-\frac{(1\pm \eta {)}^2}{(1\mp \eta {)}^2}},$$

$$\mathcal{H}=\mathrm{i}\sum _j\left[t_j{{\rm Y}}_j^B{{\rm Y}}_{j+1}^A+{\tilde{U}}_j{{\rm Y}}_j^A{{\rm Y}}_{j+1}^B\right].$$

$${\mathcal{H}}_{\mathrm{I}}$$

$${\mathcal{H}}_{\mathrm{II}}$$

$${{\rm Y}}_j^A$$

$${{\rm Y}}_j^B$$

$$a_{\mathrm{I},j}=-{\left(\frac{U(1+\eta )}{t(1-\eta )}\right)}^{j-1},a_{\mathrm{II},j}=-{\left(\frac{t(1+\eta )}{U(1-\eta )}\right)}^{j-1}.$$

$$\Xi (n)=\sum _{\alpha }|⟨{\Psi }_{\alpha }|c_n|GS⟩{|}^2+|⟨{\Psi }_{\alpha }|c_n^{†}|GS⟩{|}^2$$

$$\mathcal{H}=\sum _j\left[-t_j{\sigma }_j^x{\sigma }_{j+1}^x+{\tilde{U}}_j{\sigma }_j^z{\sigma }_{j+1}^z\right],$$

$$\mathcal{H}=\sum _j\left[-t_j{\sigma }_j^x{\sigma }_{j+1}^x+{\tilde{U}}_j{\sigma }_j^y{\sigma }_{j+1}^y\right].$$

$${\sigma }_i^x$$

$${\sigma }_i^y$$

$${\sigma }_i^z$$

$$\mathcal{H}$$

$$+\sum _j{\tilde{U}}_j\left[f_{j+1}^{†}{f}_j+f_j^{†}{f}_{j+1}-f_j^{†}{f}_{j+1}^{†}-f_{j+1}{f}_j\right].$$

$$\mathcal{H}$$

$$-\sum _k(t+\tilde{U})\left[w_k{f}_{-kA}^{†}{f}_{kB}^{†}+w_k^{\ast }{f}_{kB}{f}_{-kA}\right].$$

$$z_k$$

$$w_k$$

$$\mathcal{H}$$

$$\frac{{\Lambda }_k^2}{4}$$

$$\frac{U}{t}=\sqrt{\frac{{\delta }^2}{t^2}-\frac{(1\pm \eta {)}^2}{(1\mp \eta {)}^2}},$$

Full text

Introduction.

Many-body interactions play a crucial role in Hermitian quantum systems. The emergent correlation effects in these systems give rise to a variety of collective phenomena, such as spontaneous symmetry breaking Koma and Tasaki (1994); Belyaev (2006); Brauner (2010); Dong et al. (2017), phase transitions Santos et al. (2016); Heyl (2018); Carollo et al. (2020); Serwatka et al. (2023), and the emergence of fractionalized quasiparticles Zhang et al. (2020); Hashisaka et al. (2021); Kaskela and Lado (2021); Wouters et al. (2022). Understanding these rich phenomena often requires the combination of analytical and numerical techniques due to the scarcity of exact solutions for many-body models. Nonetheless, especially in one-dimensional systems, analytical solutions in specific regimes are attainable Gangadharaiah et al. (2011); Katsura et al. (2015); Ezawa (2017); Wang et al. (2017); Miao et al. (2017); Zvyagin (2022). Away from these parameter regimes, employing various numerical methods Lin (1990); Vidal (2004); Schollwöck (2005); Kotliar et al. (2006); Gangadharaiah et al. (2011); Stoudenmire et al. (2011); Silvi et al. (2019); Tuovinen (2021) allows unveiling the underlying physics of many-body systems in generic scenarios.

The presence of losses and dissipation in real systems provide natural platforms realizing non-Hermitian models Esaki et al. (2011); Malzard et al. (2015); Gong et al. (2018); Yao and Wang (2018); Shen et al. (2018); Kawabata et al. (2019); Yokomizo and Murakami (2019); Zhou and Lee (2019); Borgnia et al. (2020); Ashida et al. (2020). Non-Hermitian quantum models have risen as a new paradigm to manipulate and interpret various emergent phenomena Sayyad et al. (2022a); Okuma and Sato (2023). Here, Non-Hermiticity emerges as the effective description Prosen (2008); Lieu et al. (2020); Sayyad et al. (2021); Jr et al. (2022); Yang et al. (2022); Talkington and Claassen (2022); Starchl and Sieberer (2022); Gneiting et al. (2022) of out-of-equilibrium and open quantum systems, e.g., in superconducting qubits Chen et al. (2021); Abbasi et al. (2022); Chen et al. (2022), giving rise to various phenomena absent in the Hermitian counterpart. Paradigmatic examples are occurrence of various non-Hermitian degeneracies Leykam et al. (2017); Bergholtz et al. (2021); Sayyad and Kunst (2022); Sayyad et al. (2022b); Sayyad (2022) and non-Hermitian (bulk) skin states Song et al. (2019); Zhang et al. (2022). Both of these phenomenologies are mainly explored in effectively single-particle non-Hermitian Hamiltonians while non-Hermitian many-body effects have remained relatively unexplored, partly due to limitations of numerical methodsFreund (1992); Freund et al. (1993); Guo et al. (2022); Carden (2011); Zhang and Dai (2016); Chen et al. (2022). In particular, recent efforts have addressed one-dimensional non-Hermitian fermionic Fukui and Kawakami (1998); Buča et al. (2020); Zhang and Song (2021); Nakagawa et al. (2021); Yoshida and Katsura (2022); Hyart and Lado (2022) and bosonic Yamamoto et al. (2022); Wang et al. (2022) Hubbard models. Here, the non-Hermiticity is incorporated by having nonreciprocal hopping Fukui and Kawakami (1998), complex Hopping Zhang and Song (2021), or complex Hubbard interaction Buča et al. (2020); Nakagawa et al. (2021); Yoshida and Katsura (2022). The latter form of non-Hermiticity provides effective descriptions for experiments on open quantum systems with two-body loss Lewenstein et al. (2007); Zhang et al. (2016); Gross and Bloch (2017); Rosso et al. (2022).

In this Letter, combining exact analytical results and numerical calculations, we establish the phase diagram of a family of non-Hermitian Kitaev chains Rainis and Loss (2012); Lieu (2019); Sayyad et al. (2021); Sakaguchi et al. (2022). Our interacting model consists of a complex-valued many-body interaction that may host Majorana modes, in particular, realizable in an array of Josephson junctions Hassler and Schuricht (2012). Our results reveal that, depending on the relative couplings, non-Hermitian interacting phases with topological degeneracies emerge in the system. We show how increasing the non-Hermiticity parameters affect some of the many-body topological phases. We furthermore present that the topological degeneracies of the model remain in the non-analytically solvable regime by numerically solving the interacting problem.

Model.

The Hamiltonian for the non-Hermitian dimerized Kitaev-Hubbard chain, schematically shown in Fig. 1, reads

[TABLE]

where $L$ denotes the length of the chain, $c^{\dagger}_{j}\leavevmode\nobreak\ (c_{j})$ creates (annihilates) an spinless fermion at site $j$ associated with the fermion density $n_{j}=c^{\dagger}_{j}c_{j}$. Here $\mu$ adjusts the onsite energy, and $t_{j}$, $\Delta_{j}$ and $U_{j}$ are, respectively, real-valued site-dependent hopping amplitude, superconducting pairing amplitude, and Hubbard interaction. The dimerized parameter ${\cal O}_{j}\in\{t_{j},\Delta_{j},U_{j},\delta_{j}\}$ for $1\leq j\leq L$ reads

[TABLE]

where $\eta$ is the real-valued dimerization parameter and ${\cal O}\in\{t,\Delta,U,\delta\}$ stands for site-independent parameters; see also Fig. 1.

The Hamiltonian in Eq. (1) at $\mu=0$ remains invariant under $c_{j}\rightarrow(-1)^{j}c^{\dagger}_{j}$ which enforces the charge conjugation symmetry. We note that respecting this symmetry ensures that eigenvalues of the Hamiltonian come in complex-conjugate pairs making this model not directly realizable in open quantum systems Lieu et al. (2020), which can be resolved by a negative imaginary shift of all eigenvalues Joglekar and Harter (2018). Nevertheless, preserving the charge conjugation symmetry in Hermitian models ($\delta=0$) at the symmetric point $\Delta=t$ allows calculating exact phase diagrams for arbitrary $\eta$ Gangadharaiah et al. (2011); Katsura et al. (2015); Ezawa (2017); Wang et al. (2017); Miao et al. (2017). In the following, we present that there exists an analytical solution for the non-Hermitian model when the charge conjugation symmetry is respected, i.e., at $\mu=0$ and $\Delta=t$. When $\Delta\neq t$ or $\mu\neq 0$, we compute the topological phase diagram numerically.

Exact solution at $\mu=0$ and $\Delta=t$.

To obtain the exact phase diagram of the model Hamiltonian at $\mu=0$, we employ two Jordan-Wigner transformations and one spin rotation; see the Supplemental materials (SM) for details SuppMat . This procedure map our initial non-Hermitian interacting Hamiltonian into a non-Hermitian quadratic fermionic model given by

[TABLE]

which in the momentum-space casts

[TABLE]

where $\tilde{U}=U-\mathrm{i}\delta$. Diagonalizing this Hamiltonian, we obtain the energy spectrum of this four-band system given by

[TABLE]

where $\tilde{U}^{2}=U^{2}-\delta^{2}-2\mathrm{i}\delta U$. As eigenvalues of our system appear in complex-conjugate pairs, due to the charge-conjugation, zero modes in $\Lambda_{k}$ emerge when $\mathop{\mathrm{Re}}[\Lambda_{k}]=0$ at

[TABLE]

which is obtained at $k=\pm\pi/2$. In the Hermitian limit ($\delta=0$), Eq. (6) reproduces the Hermitian results Gangadharaiah et al. (2011); Katsura et al. (2015); Ezawa (2017); Wang et al. (2017); Miao et al. (2017). Eq. (6) shows the boundaries between various phases in our system, presented in yellow dashed lines in Fig. 3. These boundaries delineate phases with fourfold, twofold, and no degeneracies, respectively, shown in red, green, and blue in Fig. 2. The order of degeneracy is determined by $\chi=\sum_{i=0}^{L}\exp[-\lambda|\varepsilon_{i}-\varepsilon_{0}|]$. This quantity measures the degeneracy between first $\varepsilon_{0}$ and the $i$the smallest eigenvalue $\varepsilon_{i}$ within the energy resolution of $1/\lambda(=0.05)$. Here the eigenvalues are calculated using the numerical exact diagonalization method. The difference between the numerical boundaries and the exact solution should be attributed to the finite size effect. In the thermodynamic limit ($L=\infty$), one can recover the exact phase boundaries; see the SM SuppMat . The topological superconducting phase residing in the boundaries surrounding $U=0$ Kitaev (2001) shrinks as non-Hermiticty increases. At the critical value $\delta=t$, this topological phase fades away, resulting in the mixing of other two-fold degenerate phases. It is also worth noting that the exact phase boundaries in Eq. (6) are associated with transitions between non-Hermitian spectra with different types of the gap in the non-Hermitian effective model in Eq. (4); see also the SM SuppMat .

Many-body Majorana edge modes.

To identify Majorana modes in our model, in the next step, we rewrite the Hamiltonian in Eq. (3) in terms of Majorana fermions ${\Upsilon}_{j}^{A}=f^{\dagger}_{j}+f_{j}$ and ${\Upsilon}_{j}^{B}=\mathrm{i}(f^{\dagger}_{j}-f_{j})$.The Hamiltonian then reads

[TABLE]

We note that interactions and hopping amplitudes between the same sublattices or within each unit cell vanish. Hence, Eq. (7) can be decoupled into two independent non-interacting Kitaev chains with length $L/2$ such that ${\cal H}={\cal H}_{\rm I}+{\cal H}_{\rm II}$ where

[TABLE]

with $\Phi^{A}_{{\rm I},j}=\Upsilon_{2j-1}^{A}$, $\Phi^{B}_{{\rm I},j}=\Upsilon^{B}_{2j}$, $\Phi^{B}_{{\rm II},j}=\Upsilon^{B}_{2j-1}$, and $\Phi^{B}_{{\rm II},j}=\Upsilon^{A}_{2j}$. Introducing Majorana particles from the electron operators as $\gamma_{j}^{A}=c^{\dagger}_{j}+c_{j}$ and $\gamma_{j}^{B}=\mathrm{i}(c^{\dagger}_{j}-c_{j})$, one can show that ${\Upsilon}$ operators are products of $\gamma$ operators such that Goldstein and Chamon (2012); Yang and Feldman (2014); Kells (2015); Miao et al. (2017); McGinley et al. (2017); Ezawa (2017)

[TABLE]

[TABLE]

These relations keep the Majorana anti-commutation relations unchanged, i.e., $\{\Upsilon^{\alpha}_{i},\Upsilon^{\beta}_{j}\}=\{\gamma^{\alpha}_{i},\gamma^{\beta}_{j}\}=2\delta_{i,j}\delta^{\alpha,\beta}$. We note that $\Upsilon$ operators, which are comprised of odd (even) number of Majorana fermions ($\gamma$s), belong to subsystem ${\rm I}\leavevmode\nobreak\ ({\rm II})$ described by ${\cal H}_{\rm I}\leavevmode\nobreak\ ({\cal H}_{\rm II})$.

The quadratic Hamiltonian in Eq. (7) may host two types of boundary modes ($\cal Q$). These boundary modes are constructed from linear combinations of $\Phi^{A/B}_{\rm I/II}$ operators, i.e., ${\cal Q}^{\alpha}_{\beta}=\sum_{j\geq 0}{a}_{\beta,j}\Phi^{\alpha}_{\beta,j}$ with $\alpha=A,B$ and $\beta={\rm I,II}$. As ${\cal Q}$ consists of higher-order multiple Majorana fermions, it is dubbed ”many-body Majorana operator” McGinley et al. (2017). In the Hermitian limit, ${\cal Q}$ operators are conserved, $[{\cal H},{\cal Q}]=0$, and using the iteration procedure, one can determine the coefficients ($a$) as Fendley (2012); Ezawa (2017); McGinley et al. (2017)

[TABLE]

In non-Hermitian systems, the operator $\cal O$ is conserved if it satisfies $[{\cal H}_{R},{\cal O}]=\{{\cal H}_{I},{\cal O}\}=0$ where ${\cal H}={\cal H}_{R}+\mathrm{i}{\cal H}_{I}$ Sayyad (2022). We note that based on the structure of ${\cal H}_{\rm I/II}$ in Eqs. (8) and (9), $\{{\cal H}_{I},{\cal O}\}=0$ is by construction satisfied and fulfilling $[{\cal H}_{R},{\cal O}]=0$ results in obtaining Eq. (12). Hence, the boundary modes in our non-Hermitian system are continuously ($\delta\rightarrow 0$) connected to the zero-energy boundary modes.

The Majorana boundary mode ${\cal Q}_{\rm I}$ consists of odd numbers of higher-order Majorana operators ($\gamma$) is fermionic and satisfies $\{{\cal Z}_{2},{\cal Q}_{\rm I}\}=0$, in the infinite chain limit Fendley (2012). However, the ${\cal Q}_{\rm II}$ mode comprises even numbers of higher-order Majorana operators ($\gamma$) is bosonic as $[{\cal Z}_{2},{\cal Q}_{\rm II}]=0$ Katsura et al. (2015); McGinley et al. (2017); Ezawa (2017). Regions with fourfold degeneracies in Fig. 2 host both $({\cal Q}_{\rm I},{\cal Q}_{\rm II})$. The topological superconducting phase enclosing $U=0$ merely hosts ${\cal Q}_{\rm I}$, in agreement with Hermitian noninteracting intuition Kitaev (2001). The other two twofold degenerate phases accommodate ${\cal Q}_{\rm II}$.

Beyond exact solutions.

Let us now present the phase diagram of our system away from the integrable regime $\Delta=t$ and $\mu=0$. First, we consider $\Delta/t\in\{0.5,1.5\}$ and plot the associated phase diagrams with $\delta\in\{0.0,1.0\}$ in Fig. 3 for a finite size system with $L=16$. Similar to the phase diagram of the system at $\Delta=t$, we witness phases with fourfold, twofold, and no degeneracies. Comparing the exact phase boundaries at $\Delta=t$, in yellow dashed lines, with boundaries of regions with $2n$fold degeneracies, in red or green, we identify deformation of the phase boundaries toward $U<0$ ($>0$) for $\Delta<t\leavevmode\nobreak\ (>t)$.

At finite chemical potential and $\delta=0$, all fourfold degeneracies are lifted, and merely phases with twofold degeneracies remain in the phase diagram; see Fig. 4. The Hermitian phase boundary $U=-0.5t$ at $\eta=0$ is consistent with previous calculations on the Kitaev-Hubbard chain Mahyaeh and Ardonne (2020); see panels (a) at $L=16$, (c) in the thermodynamics limit and the SM SuppMat . Witnessing merely twofold degeneracies in phases with $U>1$ in the Hermitian limit also persists as non-Hermiticty is increased; see Fig. 4(b) at $L=16$. While no portion of the topological superconducting phase, region encircling $U=0$, is present in (b) obtained at $\mu=0.25$ and $L=16$, extrapolating the phase diagram using different system sizes, shown in (d), reveals the survival of this phase at $L=\infty$.

Experimental measurement.

We now address the signatures of the zero modes from the experimental point of view. In a tunneling experiment with a local probe Drost et al. (2017); Kempkes et al. (2019, 2018); Huda et al. (2020); Dvir et al. (2023), the conductance at zero bias $G(\omega=0,n)$ depends on the probability of extracting (injecting) an electron in site $n$ at energy $\omega$ as $\text{dI/dV}(n,\omega)\sim\sum_{\alpha}|\langle GS|c^{(\dagger)}_{n}|\Psi_{\alpha}\rangle|^{2}\delta(\omega-E_{\alpha}+E_{GS})$, where $|GS\rangle$ is the ground state and $H|\Psi_{\alpha}\rangle=E_{\alpha}|\Psi_{\alpha}\rangle$ are the many-body excited states of the system. In the presence of topological degeneracy, the ground state presents zero mode excitations that distinguish different ground states. In that scenario, the zero bias conductance at site $n$ can be written as $\text{dI/dV}(n,\omega=0)\sim\Xi(n)$ where

[TABLE]

which directly images the probability of a local excitation between the ground state and its degenerate manifold. Here, $\alpha$ runs over the ground state manifold. In particular, $\Xi(n)$ allows one to directly observe the emergence of topological zero modes associated with the topological degeneracy of the ground state. With the previous quantity, the emergence of topological zero modes associated to the topological degeneracy of the non-Hermitian model can be directly imaged. We show in Fig. 5 the local correlator computed for the interacting non-Hermitian model for different system sizes. As the system becomes larger, an edge excitation emerges in the model, which in the thermodynamic limit leads to decoupled modes between the two edges; see also the SM SuppMat . The previous quantity has been directly imaged in the realization of the current model in the topological phase with $\delta=U=0$ and $\eta\neq 0$ Drost et al. (2017), and minimal a chain with $\Delta\neq 0$Dvir et al. (2023).

Conclusion.

To summarize, we have presented a family of non-Hermitian interacting models featuring different classes of topological degeneracies. While non-interacting non-Hermitian models can be studied with conventional methodologies, the inclusion of many-body interactions renders exploring non-Hermitian systems greatly challenging. Our manuscript establishes a family of solvable interacting non-Hermitian models, providing ideal systems for benchmarking methodologies to treat interacting non-Hermitian models. Besides showing the emergence of different topological phases in the solvable limit, we provided the many-body operators accounting for the topological degeneracy of the model. We showed how non-Hermiticity modifies topological many-body models, substantially impacting the topological phases of Hermitian systems. We benchmarked our analytical construction with exact numerical calculations of the full many-body system, demonstrating that even for finite systems, the emergence of topological modes and topological degeneracies can be observed. Our results establish a versatile family of models featuring interacting topology, providing a starting point for higher dimensional solvable interacting models.

Acknowledgement.

S.S. thanks Vittorio Paeno and Abolhassan Vaezi for helpful discussions. J.L.L. acknowledges the computational resources provided by the Aalto Science-IT project, and the financial support from the Academy of Finland Projects No. 331342 and No. 336243 and the Jane and Aatos Erkko Foundation.

Appendix A Mapping the non-Hermitian interacting Kitaev chain into a noninteracting Hamiltonian

Following Refs. Ezawa (2017); Wang et al. (2017), we now present the generalization of mapping our non-Hermitian interacting model into a non-interacting Hamiltonian when $t=\Delta$ and $\mu=0$.

In the first step, we represent the fermion operators in terms of the spin operators, described by Pauli matrices, using the Jordan-Wigner transformation Sela et al. (2011); Katsura et al. (2015); Miao et al. (2017); Mahyaeh and Ardonne (2020); Liu et al. (2021); Bi et al. (2021); Zvyagin (2022) given by $c_{i}=K_{i}\sigma^{-}_{i}$, $c_{i}^{\dagger}=\sigma^{+}_{i}K^{+}_{i}$, and $2n_{i}-1=\sigma_{i}^{z}$ with $\sigma^{\pm}=\sigma^{x}\pm\mathrm{i}\sigma^{y}$ and $K_{i}=\prod_{j<i}(-\sigma^{z}_{j})$. As a result, the spin representation of the Hamiltonian in Eq. (1) casts

[TABLE]

where $\tilde{U}_{j}=U_{j}-\mathrm{i}\delta_{j}$. As the XY spin model is exactly solvable Schultz et al. (1964), one can transform the above XZ spin model, which may not be exactly solvable, in Eq. (S1) into an XY model. For this purpose, we perform a $\pi/2$ spin rotation on all spin operators around the $x$ axis using $R=\exp[-\mathrm{i}\pi\sum_{j}\sigma^{x}_{j}/2]$ Miao et al. (2017). The subsequent Hamiltonian yields

[TABLE]

To diagonalize this spin Hamiltonian, we introduce a second Jordan-Wigner transformation Wang et al. (2017); Ezawa (2017); Ding and Zhong (2021); Wada et al. (2021) given by

[TABLE]

Using the above relations, we rewrite the XY spin model as a fermionic quadratic model, which reads

[TABLE]

As the bulk-boundary correspondence is not violated in our non-Hermitian system, searching for the presence of zero modes can be done in the momentum space ($k$) using the Fourier transformation $f_{k}=\sum_{j}f_{j}\exp[\mathrm{i}kja]$ with $a(=1)$ be the lattice constant.

The quadratic Hamiltonian in the momentum space on a bipartite unit cell with sublattices $(A,B)$ casts

[TABLE]

Here $\tilde{U}=U-\mathrm{i}\delta$ and $z_{k}$ and $w_{k}$ read

[TABLE]

The non-Hermitian Hamiltonian in Eq. (S7) can be written in a matrix form using the vector operator ${\cal C}_{k}^{\dagger}=(f^{\dagger}_{kA},f^{\dagger}_{kB},f_{-kA},f_{-kB})$ as

[TABLE]

where $a_{1}=-t+\tilde{U}$ and $a_{2}=-(t+\tilde{U})$. Expressing ${\cal H}=\frac{1}{2}\sum_{k}\Lambda_{k}{\cal C}^{\dagger}_{k}{\cal C}_{k}$, four eigenvalues $\Lambda_{k}$ then read

[TABLE]

where $\tilde{U}^{2}=U^{2}-\delta^{2}-2\mathrm{i}\delta U$. The gap closure in $\Lambda_{k}$ occurs when we impose $\mathop{\mathrm{Re}}[\Lambda_{k}]=0$ resulting in

[TABLE]

which is obtained at $k=\pm\pi/2$. We emphasize that due to the charge conjugation symmetry enforcing $\mathop{\mathrm{Re}}[\Lambda_{k}]=0$ ensures $\mathop{\mathrm{Im}}[\Lambda_{k}]=0$.

While Eq. (S12) determines phase boundaries, we note that characterizing the properties of each phase cannot be addressed from the spectra of the Hamiltonian in Eq. (S10). This is because the topological character of various phases may not be preserved after performing nonlocal (Wigner) transformations McGinley et al. (2017); Ezawa (2017). Nevertheless, as non-Hermiticity does not violate the bulk-boundary correspondence, we can deduce some pieces of information regarding zero modes from effective noninteracting Hamiltonians.

For this purpose, we start with presenting the phase boundaries in Eq. (S12) as a function of $(\eta,U,\delta)$ in Fig. S1. The two-dimensional intersection of this phase boundaries at $\delta=$ $0.0$ (a), $0.5$ (b), $1.0$ (c), and $1.5$ (d) is shown in Fig. S2. The Hermitian phase diagram in Fig. S2(a) comprises seven phases including the topological superconducting (TSC), (topological) superconducting dimer ((top-)SCD), (topological) single-electron dimer ((top-)SED), charge density wave (CDW), and Schrödinger cat (CAT) phases Ezawa (2017). The CAT phase is a superposition of two superconducting states Miao et al. (2017); Ezawa (2017).

To explore the zero modes in the phase diagram, we now look at the spectra of the system with an open boundary condition at $\mu=0$. We first consider the Hermitian limit with $\delta=0.0$ and at $\eta\pm 0.2$ along the cuts shown by dashed-dotted lines in Fig. S2. The associated spectra for $\eta=-0.2(0.2)$ are shown in Fig. S3. Comparing two pairs of smallest absolute values of eigenvalues in (a, b), shown in red and blue, we realize that the CDW, CAT, and TSC phases are twofold degenerate. The spectra also exhibit fourfold degeneracies in top-SCD and top-SED; see panel (a). However, these fourfold degeneracies are lifted for $\eta>0$ in the SCD and SED phases; see panel (b).

Further looking at Fig. S2 reveals that the boundaries of the TSC phase shrink toward diminishing this phase as $\delta$ increases. Reducing the regions with the TSC phase results in lifting the separation between pairs of phases, namely (top-SCD, top-SED), (SCD, SED), and (CAT, CDW). As the order of degeneracies in both components of these pairs are identical, we still detect fourfold, twofold, and no degeneracies, respectively, in the green, white, and blue regions in all panels of Fig. S2. This can be seen in the spectra of the system at $\delta=1.0$ along cuts at $\eta=\pm 0.2$ in Fig. S4. Note that the real and imaginary parts of eigenvalues are zero at degenerate points due to the particle-hole symmetry in the system.

Aside from $2n$-fold degeneracies with $n\in\{0,1,2\}$, each non-Hermitian phase in Fig. S2 exhibits a particular non-Hermitian gap in their spectrum. We demonstrate this point by plotting the complex spectra of systems belonging to the brown, blue, lime, and white regions, shown by circle points in Fig. S2(b). In the TSC phases, the spectrum displays a real line gap shown by a red rectangle in Fig. S5 (a). The spectra exhibit the imaginary line gap within the CDW and CAT phases as exemplified in Fig S5(b). The gap becomes the point gap both in phases with fourfold or no degeneracies, e.g., in the top-SED shown in Fig S5(c) and in the SED phase presented in Fig. S5(d).

Appendix B Topological modes in the thermodynamic limit

In the previous sections, we have focused on relatively small-length systems that can be solved with exact diagonalization. In this section, we show that the edge modes associated with topological degeneracies evolve into fully decoupled modes in the thermodynamic limit. In order to surpass the length limitations of exact diagonalization, in this section, we solve the interacting quantum many-body model using a non-Hermitian tensor network formalism dmr ; ITe ; Fishman et al. (2022), targeting the ground state and lowest excited states Hyart and Lado (2022); Chen et al. (2022). We show in Fig. S6 the spatially resolved correlator computed with the non-Hermitian tensor network formalism. It is observed that as the system becomes larger, the local edge modes become more decoupled (Fig. S6(a)), giving rise to fully decoupled modes for large systems; see panels (b) and (c) in Fig. S6.

Appendix C Finite size effects and extrapolation to the thermodynamic limit

In the main text, we have focused on studying the case of finite systems numerically. Due to finite size effects, the degeneracy of the ground state is lifted due to the coupling of the topological excitations. In order to extract the topological degeneracy in the thermodynamic limit, a finite size scaling of the energy splittings can be performed. In this section, we show that with a size scaling, we can recover the exact analytic results in the interacting limit from exact numerical calculations in finite systems. The finite size scaling relies on taking a scaling for the excited state energies as $\epsilon_{\alpha}(L)=a_{\alpha}+b_{\alpha}/L$, so that in the thermodynamic limit $\epsilon_{\alpha}(L=\infty)=a_{\alpha}$. The excited state energies are taken as $\epsilon_{\alpha}=|E_{\alpha}-E_{0}|$, where $E_{\alpha}$ is the complex eigenenergy of the many-body Hamiltonian and $E_{0}$ the ground state energy. By computing the energies for a set of finite-size systems, the coefficients $a,b$ can be extracted, and the energies in the thermodynamic limit are obtained. With the extracted energy differences, the degeneracy of the ground state can be computed as $\chi=\sum_{\alpha}e^{-\lambda\epsilon_{\alpha}}$ with $\lambda$ being the energy smearing.

We present in Fig. S7 how this methodology allows obtaining the degeneracies in the thermodynamic limit. It is observed that for finite-size systems, the phase boundaries are substantially shifted and are size dependent; see panels (a), (b), and (c) in Fig. S7). Using the extrapolation scheme noted above, the correct phase boundaries known from the analytic solution are recovered, as shown in Fig. S7(d).

We further perform a similar extrapolation scheme to the Hermitian phase diagram at $\mu=0.25t$. Here we employ the exact diagonalization methods in Fig. S8 and the tensor-network formalism in Fig. S9. As it is evident, increasing the length of the chain in panels (a)-(c) results in improving the boundaries of twofold degenerate topological phases. In the thermodynamic limit, the extrapolated phase diagrams, presented in Fig. S8(d) and Fig. S9(d), recover the exact phase boundaries at $\mu=0$ and $U\geq-0.5t$. We note that the phase boundary at $\eta=0$ is located at $U=-0.5t$ (vertical white line), in agreement with the previous Hermitian calculations Mahyaeh and Ardonne (2020).