Majorana edge modes in number-conserving models with long-range interactions

TL;DR

This paper numerically investigates number-conserving long-range interacting models, providing evidence for Majorana edge modes in finite systems without mean-field assumptions, and introduces a correlation-based method to analyze their structure.

Contribution

It demonstrates the existence of Majorana edge modes in number-conserving, long-range interacting models and develops a new correlation-based approach to study their spatial structure.

Findings

Majorana edge modes appear in open chains with long-range interactions.

Energy splitting between parity sectors vanishes with increasing system size.

Nonlocal, parity-dependent edge correlations are observed.

Abstract

Topological superconductors are believed to host exotic quasiparticle excitations known as Majorana zero-modes, with much of the evidence based on BCS mean-field theory. The direct application of mean-field arguments is tenuous in finite, isolated systems relevant in some experiments. Here, we numerically study fermion number-conserving models with long-range interactions, which under periodic boundary conditions exhibit robust topological and non-topological superconductivity, tuned by the strength of interaction [1]. We find evidence that, on the topological side, Majorana edge modes appear in open chains, manifesting as the vanishing of the energy splitting between odd- and even-parity ground states with increasing system size. Additionally, off-diagonal two-point correlation functions show nonlocal, parity-dependent edge effects consistent with Majorana phenomenology. We develop a correlation-based method revealing the spatial structure of Majorana modes in this fully interacting many-body setting.

Figures

Equations

$${\mathcal{H}}_{\text{RGK}}=\sum _k{\epsilon }_k{c}_k^{†}{c}_k-\frac{16g}{L}\sum _{k,k^{\prime }}{\eta }_k{\eta }_{k^{\prime }}{c}_k^{†}{c}_{-k}^{†}{c}_{-k^{\prime }}{c}_{k^{\prime }},$$

$${\mathcal{H}}_{\text{RGK}}=-\sum _{i=1}^L\sum _{l=1}^2(t_l{c}_i^{†}{c}_{i+l}+h.c.)-\frac{4g}{L}{I}_{\varphi }^{†}{I}_{\varphi },$$

$$I_{\varphi }=2i\sum _k{\eta }_k{c}_k{c}_{-k}=\sum _{i>j}^L\eta (i-j)c_i{c}_j.$$

$$\eta (r)=\frac{8r(-1{)}^r}{\pi (1-4r^2)}.$$

$${|N+1⟩}_o$$

$${|N⟩}_e$$

$$2$$

$$0$$

$$d^{†}$$

$${|N⟩}_e$$

$${|N+1⟩}_o$$

$${|N+1⟩}_o\propto c_{k=0}^{†}{P}_N\cdot \prod _k{e}^{{\lambda }_k{\tilde{c}}_k^{†}{\tilde{c}}_{-k}^{†}}|0⟩$$

$${|N⟩}_e\propto {\tilde{c}}_{\ell ,1}^{†}{|N-1⟩}_o.$$

$$M_{jk}^o$$

$$E_0(L,N)=F(L,N)+(-1{)}^{N}G(L,N),$$

$$E_0^{(k)}(N)=2^k(-1{)}^{\frac{k}{2}}G(L,N)(-1{)}^N+O({\kappa }^{k}F)$$

$$E_0^{(2)}(N)=E_0(N+1)+E_0(N-1)-2E_0(N),$$

$$\begin{array}{r}{E}_0^{(4)}(N)=E_0(N+2)-4E_0(N+1)+6E_0(N)\\ -4E_0(N-1)+E_0(N-2)\end{array}$$

$$H_{int}=-\frac{4g}{L}\sum _{i>j,m>\ell }^L\eta (i-j)\eta (\ell -m)c_j^{†}{c}_i^{†}{c}_m{c}_{\ell }.$$

$$G(L,N)\approx -\frac{4g}{L}\sum _{i>j}^L\eta (i-j)P_{ij}^{\ast }\sum _{m>\ell }^L\eta (\ell -m)\Delta P_{m\ell }.$$

$$G(L,N)\sim g{P}_{loc}{P}_{nl}/L.$$

Full text

††preprint: APS/123-QED

In the last few decades, the physics of Majorana fermions has become an important subfield of condensed matter physics. Majorana fermions were first constructed as real solutions to the relativistic Dirac equation [2]. Formally, any complex fermion can be written as a superposition of two Majorana fermions. However, only in special cases, when a system exhibits topological order, can spatially unpaired Majorana fermions exist. These operators are self-conjugate and commute with the Hamiltonian. Hence, they represent zero-energy quasiparticles and are referred to as Majorana zero-modes (MZMs). MZMs have garnered significant theoretical and experimental interest in recent years due to their potential applications in fault-tolerant quantum computing and topologically protected quantum information storage [3, 4]. MZMs possess unique properties stemming from their non-Abelian exchange statistics and topological protection that make them robust against local perturbations and decoherence. Their potential realization in condensed matter systems such as topological superconductors has sparked intense experimental [5, 6, 7, 8, 9] and theoretical efforts [10, 11, 12, 13, 14, 15] aimed at their detection and manipulation.

The canonical example of a system supporting MZMs is the Kitaev chain [16]. Kitaev’s model is a mean-field description of a $1$D tight-binding chain with p-wave superconducting pairing. It supports a topological phase characterized by the presence of unpaired MZMs localized near the chain’s edges. Kitaev’s model forms the basis for understanding more realistic models of topological states in superconducting wires, with spinful electrons, noise and disorder [17]. However, being a mean-field model, it lacks particle number conservation, which should be an important consideration for describing finite, isolated segments of superconductors or superfluids. Significant effort has thus been dedicated to understanding systems that support MZMs and conserve particle number [18, 19, 20, 21, 22, 23, 24, 25].

Most of the considered cases included only short-range interactions, which can only sustain finite-range superconducting correlations in 1D. In contrast, long-range interactions can stabilize phase coherence over long distances even in $1$D. This has been demonstrated by Ortiz et al [1], using integrable long-range interacting Richardson-Gaudin-Kitaev (RGK) models on a ring. They convincingly demonstrated the existence of topological and non-topological superconducting phases by probing ground state fermion parity switches between periodic and antiperiodic boundary conditions. However, the integrability breaks down for open boundary conditions and hence cannot be used to test for the existence of edge modes.

In this work, we fill the knowledge gap by studying RGK models on open segments using DMRG (ITensor). First, we examine the odd-even fermion parity energy gap, $\Delta E$, in the open chain as a function of system size $L$. In the trivial phase, this difference saturates to a finite value. In the topological phase, it decays to zero as a power law $\sim 1/L$. The latter behavior can be attributed to the presence of gapless edge modes. Second, we find that single-fermion two-point correlation functions clearly distinguish the topological phase from the trivial phase. We present a correlation-based method to extract the real-space profiles of excitations that correspond to the parity switches. We find that in the topological case, these excitations have mixed particle-hole character, in close analogy to their mean-field counterparts. That is in contrast to the trivial case, where parity-switching excitations have predominantly particle character, as one would expect in the case of single-particle excitations in the BEC limit of a superconductor. Finally, we calculate observables that can be associated with the parities of Majorana edge modes and find that they are quantized, in direct correspondence with mean-field MZMs. All of these findings clarify the nature of MZMs in fully interacting, number-conserving models.

Richardson-Gaudin-Kitaev Chain. We build upon earlier work by Ortiz et al  [1] that demonstrated the transition between topological and trivial superconductivity in long-range-interacting Richardson-Gaudin-Kitaev (RGK) chains. The Hamiltonian defining the RGK chain is

[TABLE]

where $\varepsilon_{k}=-2t_{1}\cos{k}-2t_{2}\cos{2k}$ is the single-particle dispersion, $g>0$ is the interaction strength, and $c_{k}$ are annihilation operators for spinless particles with momentum $k$. The special form of the interaction term $\eta_{k}=\sin{(k/2)}\sqrt{t_{1}+4t_{2}\cos^{2}{(k/2)}}$ makes the model integrable, allowing for exact numerical solution for large system sizes both for periodic and antiperiodic boundary conditions. General boundary conditions are encoded by the phase $\phi$, $c_{i+L}=e^{i\phi}c_{i}$, with $\phi=0$ corresponding to periodic and $\phi=\pi$ anti-periodic. In real space, the Hamiltonian takes the form

[TABLE]

where

[TABLE]

For concreteness, we will study the case $t_{1}\neq 0$, $t_{2}=0$ (nearest neighbor hopping), and set $t_{1}=1$. Then, the position space pairing potential $\eta(m)$ can be determined in closed form via Fourier transform, which in the limit $L\rightarrow\infty$ is

[TABLE]

Such long-range interactions can stabilize topological phases that tend to be fragile in systems with only short-range interactions [26, 27]

Signatures of topological superconductivity and Majorana zero-modes. Within mean-field descriptions of topological superconductors, one can easily compute bulk topological invariants and directly verify the existence of MZMs. However, the situation is more subtle for interacting systems. In the closed chain, the topological phase transition can be identified from ground state fermion parity switches. The trivial phase is associated with the even ground state parity: all electrons prefer to form Cooper pairs and adding an extra particle tends to increase the energy. In the topological phase, on the contrary, the ground state parity is odd for periodic boundary conditions. Notably, changing the boundary conditions from periodic to anti-periodic has no effect on the ground state parity of the trivial phase (even), while switching it to even in the topological case. In the mean-field description of Read and Green [28], the phases track the occupancy of the zero-momentum fermion mode (empty in the trivial phase and occupied in the topological phase). This mode only exists for periodic boundary conditions and thus the spectral distinction between the trivial and topological phases disappears for anti-periodic boundary conditions.

In number-conserving models, it is also possible to track fermion parity switches under twisted boundary conditions [1]. However, care needs to be taken to eliminate the trivial density-dependent contribution to the energy, which can mask the effect. This can be done by considering the “symmetrized” odd-even ground state energy difference $\Delta E=E_{0}^{\textnormal{odd}}(N,L)-E_{0}^{\textnormal{even}}(N,L)$. For even $N$, we define $E_{0}^{\textnormal{odd}}(N,L):=\frac{1}{2}[E_{0}(N+1,L)+E_{0}(N-1,L)]$ and $E_{0}^{\textnormal{even}}(N,L):=E_{0}(N,L)$. For odd $N$, we swap the definitions of $E_{0}^{\textnormal{odd}}(N,L)$ and $E_{0}^{\textnormal{even}}(N,L)$. The averaging over $N\pm 1$ offsets the trivial (chemical potential) shift, allowing a comparison of energies with the same average density regardless of the parity sector. Ortiz et al [1] studied $\Delta E$ under periodic and anti-periodic boundary conditions to identify the topological phase transition in these models; we verify these findings in SM..1 as a consistency check.

Even when a bulk topological invariant can be defined, it may not guarantee the existence of gapless edge modes in the open chain [29, 30]. Instead, to directly probe the presence of MZMs in the open system, we analyze the symmetrized odd-even energy gap and the ground state single-particle two-point correlation matrix. As we show, these reveal qualitative differences between the trivial and topological phases.

Open RGK chain: Spectrum. A telltale sign of gapless edge modes is the behavior of the odd-even excitation gap, $\Delta E$, as a function of system size. In the topological phase, the energy cost to change parity sectors is expected to vanish in the thermodynamic limit. In contrast, in the trivial phase, no low-energy edge quasiparticles exist, and $\Delta E$ should remain finite in the thermodynamic limit.

In Fig. 1, we present numerical results for the topological and trivial phases. We choose the reference particle number as $N=L/4$ for each system size; it can be even or odd. We find that in the topological phase, for $g=1$, $|\Delta E|$ decreases as $\sim 1/L$, consistent with the presence of gapless edge modes. We note, however, that this $\Delta E$ does not have an interpretation as the quasiparticle energy cost; instead, the decay of $\Delta E$ should be attributed to a smooth ‘charging energy’ $E_{0}(N)\sim N^{2}/L$. Indeed, the energy cost to go from $N$ to symmetrized $N\pm 1$ states is always positive, regardless of whether $N$ is even or odd, consistent with the ground state energy being a nearly smooth function of $N$ with vanishing discrete curvature $E_{0}(N+1)-2E_{0}(N)+E_{0}(N-1)$ for large $L=4N$. There is, however, a much smaller underlying $(-1)^{N}$ energy oscillation, which can be associated with the MZM energy splitting. It can be extracted by taking higher discrete derivatives of the energy (SM..6), and also found to scale as $1/L$.

On the other hand, in the trivial phase, the splitting $\left|\Delta E\right|$ does not decay to zero, with the odd parity state always being higher in energy. This is consistent with the finite energy cost associated with an unpaired quasiparticle, as expected in a trivial superconductor [31].

Open RGK chain: Correlation matrix. We next examine the single-particle density matrix $M^{ij}=\left<N\right|c_{i}^{\dagger}c_{j}\left|N\right>$, which, as we will show, can directly reveal the structure of MZMs. In the ground state of a non-interacting system, the eigenmodes (and eigenvalues) of $M$ coincide with the single-particle eigenstates (and their occupation). In general, the eigenvalues are bounded between 0 and 1. Tracking changes in $M^{ij}$ between ground states at different $N$ can thus shed light on the structure of the lowest energy quasiparticle excitation connecting these states. In atomic simulators, $M^{ij}$ can be experimentally accessed via quenches [32].

In Figure 2, $M^{ij}$ is plotted for a system of size $L=32$ in four cases: $N=8$ and $N=9$, for $g=1$ and $g=3$. For the trivial case ($g=3$), there is no obvious difference between $N=8$ and $N=9$. In contrast, for the topological case, the correlators involving creation/annihilation operations at the opposite edges are large and differ in sign for $N=8$ and $N=9$. This behavior remains robust as the system size increases (comparing $N=L/4,L/4+1$).

The finiteness of these inter-edge correlators at large $N$ suggests the presence of a topologically protected edge mode.

The eigen-spectra of these matrices in the trivial and topological phases are shown in Fig. 3 (a) and (b). We also define a new matrix $\Delta M\equiv M_{o}-M_{e}$. It is constructed out of $M_{e}\equiv M_{8}$, and $M_{o}=[M_{7}+M_{9}]/2$ in the topological phase and simply $M_{o}=M_{9}$ in the trivial phase. The different choices of $M_{o}$ in the two phases lead to the smallest rank of $\Delta M$. We observe that i) the spectra for even and odd parity states are doubly degenerate except the odd parity states have additional, isolated and nearly exact [math] and $1$ eigenvalues, independent of the phase; ii) the topological phase is characterized by a spectrum spanning the full range from 0 to 1, unlike the trivial phase where all the eigenvalues (except the eigenvalue 1 for the odd parity state) lie between [math] and $0.35$; iii) $\Delta M$ is a low-rank matrix. In the trivial phase, all eigenvalues of $\Delta M$ are close to [math], except one, which is $\approx 1$. In the topological phase, there are two large eigenvalues $\approx\pm 1/2$, with eigenvectors whose weight is peaked near the edges.

The pronounced eigenvalue doubling is evidence of Cooper pairing. As we show in SM..3, a number-projected BCS wavefunction [33], $\psi\sim P_{N}e^{A_{ij}c^{\dagger}_{i}c^{\dagger}_{j}}\left|\Omega\right>$, has this feature for even $N$. The skew-symmetric pairing matrix $A_{ij}$ can be diagonalized into blocks $\lambda_{l}(-i\sigma^{y})$ associated with the $l^{\text{th}}$ pair of orthonormal modes. The correlation matrix is diagonalized by the same orthonormal modes, but with eigenvalues $\lambda_{l}^{2}/(\lambda_{l}^{2}+1)\equiv n_{l}$, corresponding to the probability of a pair being present in the condensate. The difference in the span of eigenvalues between the two phases agrees with a mean-field treatment for a periodic system. The full span in the topological phase can be related with the Zak phase acquiring a nontrivial winding (see SM..3 for details).

The simplest prescription to create a charged excitation is to unbind a pair of modes with eigenvalue $\lambda_{*}$, and set the occupation of one mode in the pair to $1$, and the other to [math]. Up to $1/N$ corrections, this would leave the spectrum of $M$ unchanged, except for the pair of eigenvalues $n_{*}$ converting into 1 and 0. Consequently, with the same accuracy, $\Delta M$ would be a rank-2 matrix with eigenvalues $1-n_{*}$ and $-n_{*}$. In the trivial phase, the data are well explained with the choice $\lambda_{*}\approx 0$, both for open and periodic boundary conditions, indicating that the excitation is essentially particle-like. In the topological phase, the $\pm 1/2$ eigenvalues of $\Delta M$ imply $\lambda_{*}\approx 1$. This indeed works well for periodic boundary conditions where a pair with eigenvalues $n_{*}\approx 1/2$ are converted into 1 and 0 upon addition of a particle (see SM..3). However, for open boundary conditions in the topological phase, there is no evidence of such an unbinding of a selected pair with $\lambda_{*}\approx 1$ (see spectra of $M_{7},M_{8},M_{9}$ in Fig. 3 (a)), motivating a more general ansatz.

Particle-hole symmetric ansatz. Inspired by mean-field considerations, we conjecture the following relations between ground states in the topological phase with open boundary conditions,

[TABLE]

with normalization $\sum_{j}\left|\alpha_{j}\right|^{2}=\sum_{j}\left|\beta_{j}\right|^{2}=1$. In the mean-field case, as well as in the number-projected Kitaev wave function [20], the excitation connecting states of different parities is an equal weight mixed particle-hole mode. In correspondence with [20], we anticipate $\alpha_{j}=\alpha_{L+1-j}$ and $\beta_{j}=-\beta_{L+1-j}$, which are significant only near the edges. Note that we interpret the ‘hole’ contribution in this number-conserving setting as a fermionic annihilation operator acting on top of Cooper pair creation, $C^{\dagger}=\sum_{N}\left|N+2\right>\left<N\right|$, on a state $\left|N\right>_{e}$.

Imposing normalization on the states in Eqs. (5) yields

[TABLE]

where we introduced the anomalous correlation matrix $P^{jk}_{e/o}=\langle N+2|c^{\dagger}_{j}c^{\dagger}_{k}|N\rangle$ for $N$ even/odd, and defined $M\equiv(M_{e}+M_{o})/2$ and $\Delta M\equiv M_{o}-M_{e}$ (and analogously for $P$) in anticipation that $\Delta M$ and $\Delta P$ isolate the inter-edge correlations associated with the topological phase in the large $N$ limit. Einstein summation over repeated indices is implied. Given the observed low-rank of $\Delta M$ and the eigenvalues $\pm 1/2$, it is natural to consider the ansatz $\Delta M^{jk}=(\beta_{j}\beta_{k}-\alpha_{j}\alpha_{k})\tilde{M}$ with $\tilde{M}\approx 1/2$. Eq. (6) is then satisfied by choosing $\Delta P^{jk}=(\beta_{j}\alpha_{k}-\beta_{k}\alpha_{j})\tilde{P}$, with $\tilde{P}=1/2$. For number-projected Kitaev wavefunctions, this symmetry between regular and anomalous correlators is present [20]. More generally, there is a physical principle why we can expect $\tilde{M}=\tilde{P}$, and thus $\tilde{M}=1/2$, in agreement with the numerically obtained eigenvalues of $\Delta M$. A non-vanishing nonlocal correlator $\Delta M^{1,L}=\left<N\right|{c^{\dagger}_{1}c_{L}}\left|N\right>$ seems to allow instantaneous teleportation of a particle from one edge of the system to another, while remaining in the ground state [34, 35]. To avoid this, it is necessary that the state $c^{\dagger}_{1}\left|N\right>$ have equal overlap with both $c^{\dagger}_{L}\left|N\right>$ and $c_{L}\left|N+2\right>$, which requires that inter-edge correlators in $M_{ij}$ and $P_{ij}$ (represented by $\Delta M$ and $\Delta P$) are identical in amplitude, giving $\tilde{M}=\tilde{P}$.

MZM operators. The above analysis points to the existence of a fermionic mode whose occupation number switches between even- and odd-parity states. The creation operator associated with this mode is given by

[TABLE]

which acts as $\left|N+1\right>_{o}=d^{\dagger}\left|N\right>_{e}$ and $\left|N\right>_{e}=d\left|N+1\right>_{o}$, as can be seen from (5). The first term preferentially creates extra density at the edges, while the second term first adds an extra Cooper pair, leading to a $2/L$ change in density uniformly throughout the system, and then removes density near the edges. Given the equal amplitude of particle and hole contributions at the edges, the net effect is an extra particle smeared uniformly through the system. Thus, even though $d$ nominally looks like a local operator, its effect is global. Moreover, the change in density $\sim 1/L$ vanishes in the thermodynamic limit and thus does not lead to locally observable consequences, as should be the case in the topological phase.

We can now relate the operator $d$ to the enigmatic MZM operators. Due to the bimodal support of $d$ and the special relationship between the vectors $\alpha$ and $\beta$, the operators $\Gamma^{\dagger}_{L}=d^{\dagger}+dC^{\dagger}=1/\sqrt{2}\sum_{j}(\alpha_{j}-\beta_{j})(c_{j}^{\dagger}+c_{j}C^{\dagger})$ and $\Gamma^{\dagger}_{R}=i(d^{\dagger}-dC^{\dagger})=i/\sqrt{2}\sum_{j}(\alpha_{j}+\beta_{j})(c_{j}^{\dagger}-c_{j}C^{\dagger})$ appear to be localized to a single edge each. Thus, the combinations $\alpha_{j}\pm\beta_{j}$ can be thought of as MZM wavefunctions (shown as green and red lines in Fig. 3d). However, just as in the case of the operator $d$, the effect of $\Gamma_{L,R}$ is global due to the presence of the Cooper pair operator, increasing the density in the system uniformly by $1/L$.

We finally note that the canonical self-conjugate Majorana operators $\gamma_{L}=1/\sqrt{2}\sum_{j}(\alpha_{j}-\beta_{j})(c_{j}^{\dagger}+c_{j})$, $\gamma_{R}=i/\sqrt{2}\sum_{j}(\alpha_{j}+\beta_{j})(c_{j}^{\dagger}-c_{j})$ retain significance beyond mean-field theory. We find that the MZM parity operator $\left<N\right|i\gamma_{R}\gamma_{L}\left|N\right>$ is nearly quantized when evaluated over fixed number states as $\approx 0.5\times(-1)^{N}$ [e.g., $95\%$ ($99.9\%$) accuracy for $L=48$ and $N$ even (odd)]. The value is only half-integer because it misses the anomalous contribution of the same magnitude.

Conclusion and outlook. In this work, we have investigated the manifestations of MZMs in a number-conserving setting. While we focused on the interacting Richardson-Gaudin-Kitaev chain, we expect that our results are more general. We found that the topological phase, unlike the trivial phase, is highly sensitive to the boundary conditions: the sign of the odd-even parity energy gap changes between periodic and antiperiodic boundary conditions, and vanishes with increasing system size for open boundary conditions. Even more remarkably, the single-particle correlation matrix acquires non-local parity-dependent inter-edge contributions that remain large as the system size increases. From the correlation matrix, we reconstructed the non-local fermionic mode that converts between different number states and can be related to the Majorana modes. The sensitivity to boundary conditions and the existence of nonlocal boundary modes is most likely related to the Cooper pair wavefunction becoming infinitely-ranged (only algebraically decaying) in topological superconductors [28]. Both this exotic behavior of the pair wavefunction and the MZM modes may be accessible by measuring the correlation matrix in cold atom gas experiments. We leave detailed study of the pair wavefunction, more general interactions cases, and braiding for future work.

Acknowledgments: We acknowledge discussion with S. K. Yip and E.-A. Kim. JTM would like to thank the University of Chicago REU program. IM and KA acknowledge support by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

Supplemental Material

.1 Spectra of closed RGK Chain

Consider slowly changing the boundary conditions from periodic to antiperiodic. The topological invariant $Q$ can be defined as the difference in effective ground state fermion parities between chains with periodic and antiperiodic boundary conditions [1]. In particular, an odd number of gap closings corresponds to a nontrivial $Q$, signaling that the superconductor has nontrivial topology. We illustrate this behavior in Fig. 4.

Fig. 5 shows the odd-even excitation gap $\Delta E$ (defined in the main text) as a function of system size $L$ for periodic and antiperiodic boundary conditions, for $g=1$ and $g=3$. For $g=1$, in the topological phase, $\Delta E$ does not decay to zero. Under periodic boundary conditions, the system favors the odd-parity ground state, regardless of the parity of $N$, even though there are some oscillations as a function of $N$, with the odd-parity state being somewhat more favored for odd $N$ relative to even $N$. For antiperiodic boundary conditions, however, the favored ground state flips — the system now prefers the even-parity ground state. This boundary-condition–sensitive parity inversion can be viewed as a signature of nontrivial topology.

On the other hand, in the trivial phase realized for larger $g$, the even-parity sector is favored for any $N$ and any boundary conditions. This aligns with conventional superconducting behavior and highlights the distinction between the phases.

.2 Evidence for BEC in trivial phase

In Fig. 6, we observe a dramatic change in $\langle c_{i}^{\dagger}c_{i+1}\rangle$ across the phase boundary: the correlator drops by two orders of magnitude between $g=1$ and $g=3.$ This correlator measures kinetic energy.

There are two possible interpretations of this effect: the formation of small “heavy” Cooper pairs or macroscopic phase separation. In the former case, the state would remain a (BEC) superconductor, while in the latter it would be an insulator. We note that the model with only short-range attraction exhibits phase separation at strong coupling (the Ising regime of the corresponding XXZ spin model [20]). However, in long-range interacting RGK models, where interaction alternates in space between attraction and repulsion (Eq. 4), the phase separation can be averted in favor of BEC. Indeed, the preference for even fermion parity states that we found both for closed and open boundary conditions for $g=3$ support the Cooper pair BEC scenario. So does the significant superfluid stiffness, which manifests in the energy dependence on the boundary conditions under phase twist (Fig 4b) [36]. The kinetic energy in this case is determined by the hopping of Cooper pairs, which for the nearest-neighbor Cooper pairs can be estimated as $t^{2}/U$, where $U\sim g$ is the nearest-neighbor attraction  [37].

.3 Number-projected BCS ansatz

The limited change in the spectrum of the correlation matrix $M$ between different parity states in the trivial phase hints at an underlying mean-field description, whereby the change in the correlation matrix occurs due to the occupation/emptying of a single mode. Following BCS, we propose even and odd particle number wavefunctions

[TABLE]

Here, $\sum_{ij}A_{ij}c^{\dagger}_{i}c^{\dagger}_{j}$ is the Cooper-pair creation operator (with $A_{ij}$ a skew-symmetric matrix), and the projection operator $P_{N}$ for $N$ even simply selects the state with $N/2$ Cooper pairs. This wavefunction is a number projected mean-field BCS state. The operators $\tilde{c}_{l,1},\tilde{c}_{l,2}$ represent orthonormal modes in the canonical basis of $A_{ij}$, with eigenvalues $\lambda_{l}>0$. The odd parity state can be obtained by breaking one, say the $\ell^{\text{th}}$ ‘component’ of the Cooper pair wavefunction and fixing the occupation probability of the modes $\tilde{c}^{\dagger}_{l=\ell,1}$ and $\tilde{c}^{\dagger}_{l=\ell,2}$ to $1$ and [math], respectively. We note here that on a ring with periodic boundary conditions there are two special modes, $k=0$ and $k=\pi$, which lack time-reversed partners, and should be treated separately [1, 24].

The wavefunction $\left|N\right>_{e}$ yields a correlation matrix $M$ with a fully doubly degenerate spectrum and eigenvalues $\lambda_{l}^{2}/(1+\lambda_{l}^{2})$ (in the large $N,L$ limit). The derivation is trivial for the non-projected BCS wavefunction; then, since that the operator $c^{\dagger}_{i}c_{j}$ does not change the number sector, the result must match the most probable number projected states (ones with $N\approx\langle N\rangle$), up to $\mathcal{O}\left(1/N\right)$ corrections to the eigenvalues.

In the odd parity state, the spectrum should remain nearly identical, except for the doubly-degenerate eigenvalues associated with the $\ell^{\text{th}}$ pair of states being replaced by eigenvalues [math] and $1$. $\Delta M$ is then a low-rank matrix with non-zero eigenvalues $-\lambda_{\ell}^{2}/(1+\lambda_{\ell}^{2}),1/(1+\lambda_{\ell}^{2})$. Within this approximation, the orthonormal modes are identical for $M_{e},M_{o}$ and $\Delta M$.

We find that the above ansatz works well in the trivial phase for both open and periodic boundary conditions. In the odd parity state, the mode with $k\approx 0$, which was empty in the even $N$ state, becomes occupied (appearing with an eigenvalue $1$ in the correlation spectrum).

The spectrum of $M$ for the topological phase with periodic boundary conditions can be explained equally well with this ansatz, with the modification that the $k=0$ mode remains unpaired and occupied in both the odd and even parity states. That makes the odd-parity state of the form

[TABLE]

a good reference in this case. The Cooper pair operator includes all momenta except [math] and $\pi$, which do not have time-reversed partners. This odd-parity ground state does not have any excited quasiparticles. Instead, it is the even-parity state that has one Cooper pair unbound, with the wavefunction

[TABLE]

In Fig. 7, we show how the spectrum changes for $L=32$ between $N=8$ and $9$ with periodic boundary conditions. In the odd parity state ($N=9$), there are two distinct modes at $n=0$ and $1$, which upon inspection are found to be $k=\pi$ and [math] modes, respectively. Also, this state has a pair of modes near $n\approx 0.5$. These are a pair of modes with intermediate values of momenta, $(k_{*},-k_{*})$. They seem to vanish in the $N=8$ state. However, in reality, they split up, ending near $n=0$ and 1, as expected from our ansatz for the topological state with periodic boundary conditions.

In contrast, the topological state with open boundary conditions, which is the focus of the main text, presents a significant challenge. As can be seen from Fig. 3a, the difference between odd and even $N$ states is subtle, with the spectra changing very gradually with $N$. Moreover, the eigenmodes of $\Delta M$ with eigenvalues $\pm 1/2$ are not eigenstates of $M_{e},M_{o}$, which is also inconsistent with the simple projected BCS ansatz. For certain, there is no pair unbinding near the Fermi level (near occupancy $n\sim 0.5$). The notable difference between the odd and even cases is that in the odd-number state, there are two distinct unpaired modes with occupancies very close to [math] and $1$. We discuss the nature of these modes in Sec. .4. It turns out that the latter (fully occupied) mode is a direct descendant of the $k=0$ mode in a closed system, while the empty mode descends from the $k=\pi$ mode (see Fig. 8).

Zak phase.—

The span of the eigenvalues $\lambda_{l}$ differs in the two phases: in the topological phase, the eigenvalues $n_{l}=\lambda^{2}_{l}/(\lambda^{2}_{l}+1)$ span the entire range from [math] to $1$ smoothly in both even and odd parity states. This is in correspondence with mean-field expectations. Under periodic boundary conditions, the eigenmodes are pairs with momenta $(k,-k)$. The eigenvalues associated with these modes $\lambda_{k}$ are $v_{k}/u_{k}$, where the coefficients $v_{k},u_{k}$ are the usual Bogoluibov coefficients that appear in the BCS wavefunction $\psi\sim\prod\limits_{k}\left(u_{k}+v_{k}c^{\dagger}_{k}c^{\dagger}_{-k}\right)\left|0\right>$. For example, from the mean-field solution of Kiteav model, $\lambda_{k}=\frac{v_{k}}{u_{k}}=\frac{E_{k}-\xi_{k}}{E_{k}+\xi_{k}}=\frac{\Delta_{k}^{2}}{(E_{k}+\xi_{k})^{2}}$, where $\xi_{k}=-2t\cos k-\mu,\Delta_{k}=\Delta\sin k$ and $E_{k}=\sqrt{\Delta^{2}_{k}+\xi^{2}_{k}}$. The topological phase is obtained for $-2t<\mu<2t$, that is, while the chemical potential lies inside the band. Also, $\Delta_{k}=0$ at the edges of the band (at $k=0,$ $\pi$). This implies $E_{k}=\left|\xi_{k}\right|$ at the edges of the band as well. Clearly, in the trivial phase, the demoninator $E_{k}+\xi_{k}$ is always finite for $\mu<-2t$, as $\xi_{k}>0$ for all $k\in(0,\pi)$. Thus, $\lambda_{k}$ is always finite and the range of eigenvalues in the correlation matrix does not approach $1$ (in the even $N$ state). In the topological phase, the denominator must vanish at one edge of the band (and faster than the numerator), while the numerator vanishes at the other edge. This implies that the eigenvalues $\left|\lambda_{k}\right|$ span $(0,\infty)$ and the eigenvalues of the correlation matrix span $(0,1)$.

The full span of the correlation matrix spectrum correlates with the non-trivial winding of the Zak phase of the single-particle wavefunction $\tan^{-1}\left(\Delta_{k}/\xi_{k}\right)$. The Zak phase winds from [math] to $2\pi$ as momentum runs from $k=-\pi$ to $k=\pi$ in the topological phase and does not wind in the trivial phase.

.4 The special eigenvalue-1 mode in the topological phase with odd $N$

Even though the simple projected BCS ansatz does not fully capture the odd-parity state in the topological phase, a careful examination of the eigenstates of $M_{o}$ reveals the appearance of an unpaired fully populated fermion mode. The spatial structure of this mode is similar to the fermion mode in [24], which converts between even and odd fermion number Kitaev ground states, $\left|\Psi^{N+1}_{o}\right>\propto a_{1}^{\dagger}\left|\Psi_{e}^{N}\right>$. They found that $a_{1}$ is an equal superposition of all site modes. Indeed, this is close to what we find numerically for the eigenvector with the eigenvalue $\approx 1$ in Fig. 8.

It is easy to prove that if $\left|\Psi^{N+1}_{o}\right>\propto a_{1}^{\dagger}\left|\Psi_{e}^{N}\right>$, there must be an eigenvalue $1$ in the correlation matrix. First let’s normalize, $\left|\Psi^{N+1}_{o}\right>=a_{1}^{\dagger}\left|\Psi_{e}^{N}\right>/\sqrt{C}$, where $C=\left<\Psi_{e}^{N}\right|a_{1}a_{1}^{\dagger}|\Psi_{e}^{N}\rangle$. Now, expand site fermions $c_{j}$ in terms of orthogonal modes $a_{1},\ a_{2},...,a_{L}$, where all modes except $a_{1}$ are chosen arbitrarily, $c_{j}=A_{jn}a_{n}$, and $A$ is an orthogonal matrix. Then

[TABLE]

where $G_{nm}=\left<\Psi^{N}_{e}\right|a_{1}a_{1}^{\dagger}a^{\dagger}_{n}a_{m}\left|\Psi^{N}_{e}\right>/C$ is a symmetric matrix of size $L-1\times L-1$. The second term therefore is rank-deficient, spanned by the vectors orthogonal to $A^{-1}_{1j}$. But $A^{-1}_{1j}$ are nothing but $v_{j}$ in the expansion $a_{1}=v_{j}c_{j}$. Thus, we have proven that the mode $a_{1}^{\dagger}$ indeed has an associated eigenvalue $=1$ in the correlation matrix.

We have verified that the mode $a_{1}^{\dagger}=v_{j}c_{j}^{\dagger}$ converts between even and odd parity states by comparing $\left<N\right|a_{1}c_{i}^{\dagger}c_{j}a_{1}^{\dagger}\left|N\right>/\left<N\right|a_{1}a_{1}^{\dagger}\left|N\right>$ with $\left<N+1\right|c_{i}^{\dagger}c_{j}\left|N+1\right>$ (see Fig. 9). However, it cannot be considered a Majorana mode. It is easy to check that the weight $C$ is low: $\left<N\right|a_{1}a_{1}^{\dagger}\left|N\right>\sim O(1/L)$, which also implies that $\left<N\right|a_{1}^{\dagger}a_{1}\left|N\right>\sim 1-O(1/L)$. We stress here that the ansatz of Eq. (5) conjectures an equality that identifies a normalized quasiparticle creation operator in contrast to the action of $a_{1}$ on $\left|N\right>_{e}$. We have verified numerically that for $g=1$, $v^{T}M_{e}v\approx 0.97$. Since $v$ is normalized, the $a_{1}$ mode is already $97\%$ occupied in the even state. Acting on the even state with $a_{1}^{\dagger}$ promotes the occupancy to 100%. The essentially uniform change in the density induced by $a^{\dagger}_{1}$ is consistent with the expectation that the states with odd and even parities are locally indistinguishable in the thermodynamic limit.

The low efficiency of this operator is due to the fact that for the internal sites $i$ away from the edges, the expectation value $\left<\Psi^{N+1}_{o}\right|c^{\dagger}_{i}\left|\Psi^{N}_{e}\right>$ vanishes rapidly, which can be seen explicitly in the example of the Kitaev wave function [20]. In Fig. 10, we plot $\left<N\right|a_{1}c^{\dagger}_{i}\left|N\right>/\sqrt{C}\approx\left<N+1\right|c^{\dagger}_{i}\left|N\right>$ which shows the same qualitative behavior.

Despite their small contribution, dropping the internal sites from the expression for $a_{1}^{\dagger}$ while keeping the operator normalized would lead to exaggerated amplitudes near the edges, “overfilling” them and causing inhomogeneous density distribution inconsistent with topological indistinguishability of odd and even ground states. If we insist on the operator acting “locally” near the edges, we would need to offset the edge charge creation by charge annihilation. This is the physical origin of the hole contribution in the operator $d$ in Eq. (8)—it creates charge throughout the system and depletes it at the edges. Of course, $d^{\dagger}$ is not truly localized to the edges since it involves the Cooper pair creation operator that affects the density in the whole chain.

.5 Anomalous correlators

In Fig. 11, we plot the anomalous correlation matrix $P^{jk}_{e/o}=\langle N+2|c^{\dagger}_{j}c^{\dagger}_{k}|N\rangle$ for even/odd $N$. It has a structure similar to $M^{jk}_{e/o}$, except it is antisymmetric due to Fermi exchange statistics. Notably, off-diagonal correlators are large, consistent with the no-teleportation argument (see main text). This is made clearer by Fig. 12, where we plot $\Delta M^{jk}$ and $\Delta P^{jk}$, which isolate the off-diagonal structure. We note here that there is a gauge freedom associated with these anomalous correlation matrices due to the lack of $U(1)$ symmetry. In particular, since they are real, the sign can be chosen arbitrarily to make short-range correlations independent of the parity.

Using the anomalous correlators, we can numerically verify the normalization conditions Eqs. 6 and 7. They evaluate to $2.38$ and $-0.03$, respectively. These deviations from the expected values are not surprising, given the non-negligible difference in filling fraction between $N=8$ and $N=10$, as well as finite system size and DMRG truncation effects.

.6 MZM overlap energy

The apparent smoothness of the discrete second derivative of $E_{0}(N)$ in Fig. 1 suggests that $E_{0}(N)$ itself is a smooth function of $N$. However, a more careful analysis reveals that $E_{0}(N)$ is in fact not smooth. There are tiny oscillations with the parity of $N$, which may be anticipated given that for finite $L$, the MZM overlap is finite, albeit small, and the occupation number of the MZM alternates between even and odd $N$. To isolate the parity-sensitive component of $E_{0}(N)$, we write

[TABLE]

where the functions $F$ and $G$ are smooth functions of their arguments, with $|G|\ll|F|$, as revealed by numerics in Fig. 1. The smooth part can be decomposed into Fourier harmonics, each with a characteristic wave number. Taking a discrete derivative of $F$ amounts to multiplication of each mode by its characteristic frequency, effectively suppressing its contribution. In contrast, the effect of each successive discrete derivative on the second term is to amplify it by a factor of 2, yielding

[TABLE]

for even derivatives. The second discrete derivative,

[TABLE]

evaluated at $L=4N$ (Fig. 1) appears smooth. In contrast, the fourth derivative

[TABLE]

(again evaluated at $L=4N$) reveals clear oscillations with $N$ (Fig. 13). The oscillations can be attributed to tiny oscillations in $E_{0}(N)$ originating from occupying and de-occupying the MZM as the parity of $N$ changes. The magnitude of the “overlap energy,” $G(L,N)\sim|E^{(4)}_{0}(N)/2^{4}|$, is indeed quite small even for our limited sizes accessible with DMRG.

Of course, an even better approximation for the overlap energy could be obtained by taking more derivatives (e.g., $E^{(n)}_{0}(N)/2^{n}$), but $E^{(4)}_{0}(N)$ already exhibits clear sign oscillations, with $|E^{(4)}_{0}(N)|$ being nearly smooth.

We can estimate the value of $G$ by isolating the energy difference between the odd- and even-parity states, using our knowledge of the correlation matrices. Since the kinetic energy operator is spatially local, it is unimportant for evaluating this energy difference, which manifests in long-range correlators. We thus focus on the interaction energy

[TABLE]

We seek the topological contribution to the energy difference. For estimation purposes, we assume that the state is approximately Gaussian (mean-field like), and thus can be Wick-decoupled in the normal and anomalous channels. The topological splitting then comes from terms proportional to $M\Delta M$ and $P\Delta P$. The latter contributes to the odd-even splitting as

[TABLE]

For estimation purposes, we use the fact that $\eta(r)$ scales as $(-1)^{r}/r$. Then the first sum scales with the system size as $LP_{loc}$, while the second sum scales as $P_{nl}/L$. The contribution from the regular correlation matrix, $M\Delta M$, is parametrically smaller (by $1/L$) due to the special form of the interaction, and hence can be neglected. Finally,

[TABLE]

The $1/L$ scaling of the odd-even splitting is indeed consistent with the slow decay observed in Fig. 13 (remember, $L=4N$). Quantitatively, using the values from Fig. 11, $P_{loc}\sim P_{nl}\sim 0.1$ for $L=32,\ N=8$, $G(32,8)\sim 10^{-2}/32\sim 3\times 10^{-4}$. This has the correct scale, comparable with the obtained $E^{(4)}_{0}/2^{4}$. We thus conclude that the slow decay of the MZM splitting in the RGK chain can be attributed to the long-range nature of the interactions.