Number-conserving interacting fermion models with exact topological superconducting ground states

TL;DR

This paper introduces a method to construct number-conserving Hamiltonians with exact topological superconducting ground states, including models with Majorana modes, advancing understanding and potential experimental realization.

Contribution

The authors develop a general technique to create number-conserving Hamiltonians that exactly reproduce chosen BCS-type states, including topological superconductors like Kitaev wires and 2D p+ip models.

Findings

Constructed 1D double wire model with Majorana-like ground states

Developed 2D p+ip superconducting model with vortex states

Provided analytic expressions for topologically degenerate ground states

Abstract

We present a method to construct number-conserving Hamiltonians whose ground states exactly reproduce an arbitrarily chosen BCS-type mean-field state. Such parent Hamiltonians can be constructed not only for the usual $s$-wave BCS state, but also for more exotic states of this form, including the ground states of Kitaev wires and 2D topological superconductors. This method leads to infinite families of locally-interacting fermion models with exact topological superconducting ground states. After explaining the general technique, we apply this method to construct two specific classes of models. The first one is a one-dimensional double wire lattice model with Majorana-like degenerate ground states. The second one is a two-dimensional $p_x+ip_y$ superconducting model, where we also obtain analytic expressions for topologically degenerate ground states in the presence of vortices. Our models may provide a deeper conceptual understanding of how Majorana zero modes could emerge in condensed matter systems, as well as inspire novel routes to realize them in experiment.

Figures

Equations

$$K_{\mathrm{mf}}$$

$$|G⟩=\text{sideset}{}^{\prime }{\prod }_p{\alpha }_p{\alpha }_{\bar{p}}|0⟩,$$

$${\hat{A}}_{p{p}^{\prime }}=S_p^{†}{\alpha }_{p^{\prime }}+S_{p^{\prime }}^{†}{\alpha }_p.$$

$$\hat{H}=\sum _{p_1{p}_2{p}_3{p}_4}{H}_{p_1{p}_2;p_3{p}_4}{\hat{A}}_{p_1{p}_2}^{†}{\hat{A}}_{p_3{p}_4},$$

$${\hat{A}}_{p{p}^{\prime }}=S_p^{†}{C}_{p^{\prime }}+S_{p^{\prime }}^{†}{C}_p,$$

$$|G^e⟩=\exp \left(\sum _{i<j}{c}_i^{†}{c}_j^{†}\right)|0⟩,|G^o⟩={\tilde{c}}_{k=0}^{†}|G^e⟩,$$

$${\alpha }_{k\sigma }$$

$$H_{p_1{p}_2;p_3{p}_4}=\frac{16}{L^2}\sum _{j=1}^L\sin (k_{1}j)\sin (k_{2}j)\sin (k_{3}j)\sin (k_{4}j)$$

$$\times \left[p{\delta }_{{\sigma }_1{\sigma }_2{\sigma }_3{\sigma }_4}+q{\delta }_{{\sigma }_1{\sigma }_2}{\delta }_{{\sigma }_3{\sigma }_4}-r({\delta }_{{\sigma }_1{\sigma }_3}{\delta }_{{\sigma }_2{\sigma }_4}+{\delta }_{{\sigma }_1{\sigma }_4}{\delta }_{{\sigma }_2{\sigma }_3})\right],$$

$$\hat{H}$$

$$|G_N^{ee}⟩={\left[\sum _{i<j}(a_i^{†}{a}_j^{†}+b_i^{†}{b}_j^{†})\right]}^{\frac{N}{2}}|0⟩,|G_N^{oo}⟩={\tilde{a}}_0^{†}{\tilde{b}}_0^{†}|G_{N-2}^{ee}⟩.$$

$$\hat{K}={\int }_S\left[\frac{\nabla {\psi }_z^{†}\cdot \nabla {\psi }_z}{2m}-(\Delta {\psi }_z{\partial }_{\bar{z}}{\psi }_z+H.c.)+\mu {\psi }_z{\psi }_z^{†}\right]d^{2}z,$$

$${\hat{K}}_{\mathrm{bulk}}$$

$${\hat{K}}_{\mathrm{bound}}$$

$$|G^e⟩=\exp \left[\frac{1}{2}{\int }_{S}g(z,z^{\prime }){\psi }_z^{†}{\psi }_{z^{\prime }}^{†}{d}^{2}z{d}^2{z}^{\prime }\right]|0⟩,$$

$$2{\partial }_{\bar{z}}g(z,z^{\prime })={\delta }^2(z-z^{\prime }),$$

$${\hat{H}}_{\mathrm{bulk}}={\int }_{S}W(z_1,z_2;z_3,z_4)A_{z_1,z_2}^{†}{A}_{z_3,z_4}\prod _{j=1}^4{d}^2{z}_j,$$

$$\hat{H}={\hat{H}}_{\mathrm{bulk}}+{\hat{K}}_{\mathrm{bound}}.$$

$$\hat{H}$$

$$g_{12,34}(z,z^{\prime })$$

$${\psi }_{12,34}(z_1,z_2\dots z_N)=\mathrm{Pf}\{g_{12,34}(z_i,z_j)\},$$

$$H_{\mathrm{Kitaev}}=\sum _{j}t(-c_j^{†}{c}_{j+1}+c_j{c}_{j+1}+H.c.)-\mu \hat{N}.$$

$${\alpha }_k=\sum _{j=1}^L(u_j^k{c}_j+v_j^k{c}_j^{†}).$$

$$(E_k+\mu )u_j^k$$

$$(E_k-\mu )v_j^k$$

$$u_0^k=v_0^k,u_{L+1}^k=-v_{L+1}^k.$$

$$u_j^k=\lambda (k)e^{i(kj-{\theta }_k)}-\lambda (-k)e^{-i(kj-{\theta }_k)},v_j^k=e^{i(kj-{\theta }_k)}-e^{-i(kj-{\theta }_k)}$$

$$k(L+1)=2{\theta }_k+m\pi ,m\in \mathbf{Z},\tan {\theta }_k=\frac{2\Delta \sin k}{E_k+\mu +2t\cos k},$$

$${\alpha }_k=e^{i{\theta }_k}(i\cos {\theta }_k{c}_{-k}^{†}-\sin {\theta }_k{c}_k)-(k\to -k).$$

$$E_k=2t,{\theta }_k=\frac{k}{2},k=\frac{m\pi }{L},m=1\dots (L-1),$$

$$c_k|G^e⟩=\frac{1}{\sqrt{L}}\sum _{j=1}^L{e}^{-ikj}{c}_j\exp \{\sum _{i<j^{\prime }}{c}_i^{†}{c}_{j^{\prime }}^{†}\}|0⟩=\frac{1}{\sqrt{L}}\sum _{j,j^{\prime }}{e}^{-ikj}\mathrm{sgn}(j^{\prime }-j)c_{j^{\prime }}^{†}|G^e⟩=[i\cot \frac{k}{2}{c}_{-k}^{†}-\frac{1+(-1{)}^m}{1-e^{ik}}{c}_{k=0}^{†}]|G^e⟩,$$

$$[e^{i\frac{k}{2}}\sin \frac{k}{2}{c}_k-(k\to -k)]|G^e⟩=[e^{i\frac{k}{2}}i\cos \frac{k}{2}{c}_{-k}^{†}-(k\to -k)]|G^e⟩,$$

Full text

Introduction.

Topological superconductors have become an active research area in condensed matter and cold atom physics Kitaev ; TPSC-RMP . From a fundamental viewpoint, they provide examples of topological phases that have been classified systematically TPorder1 ; TPorder2 ; TPorder3 ; TPorder4 ; TPorder5 ; TPorder6 . At a practical level, the non-Abelian statistics MR1991 ; Nayak1996 ; Ivanov of Majorana zero modes and the robustness of degenerate ground states against local perturbations have made topological superconductors components of promising architectures for fault-tolerant quantum computation TPQC-RMP . Experimental signatures Exp1 ; Exp2 ; Exp3 ; Exp4 of Majorana zero modes and topological superconductivity in solid state systems call for more realistic theoretical descriptions of the relevant physics, for example the effect of interactions.

Most of the theoretical research in this area has begun with non-interacting mean-field Hamiltonians with an effective $p$-wave pairing term, from which one obtains topologically-protected degenerate ground states and Majorana zero modes that emerge as effective quasiparticles. In nature, however, quantum systems typically contain sizable interacting terms that challenge the validity of mean-field theory. It is therefore important to understand better the criteria for the aforementioned topological phenomena to persist in the interacting case. Moreover, the mean-field approximation breaks the number conservation, which obscures the connection to realistic, number-conserving systems.

Understanding the interplay of topology, number conservation, and interactions is challenging because the interactions usually prevent exact solution by analytic or numerical techniques. In one dimension, special tools are available, and progress has been made using bosonization Bos1 ; Bos2 ; Bos3 and numerical methods [density-matrix renormalization group (DMRG)] PZoller . Exactly solvable models in this area are still rare InteractingKitaev , and the three number-conserving Majorana models RGK ; SDiehl ; NLang that have been proposed are one-dimensional (1D). Having new families of exactly solvable models with realistic local interactions, especially in higher dimensions, will therefore shed light on the characterization of topological phenomena and Majorana zero modes in intrinsically interacting and number-conserving systems. In addition, these results provide new Hamiltonians that can be used to experimentally realize topological states.

In this letter we take a bottom-up approach to these fundamental issues: Starting from a general BCS-type mean-field ground state $\ket{G}$, we show how to construct number-conserving parent Hamiltonians that have $\ket{G}$ as a ground state (with no approximation). This construction enables us to realize the physics of Majorana zero modes in interacting number-conserving systems in an exact manner. Following the general construction, we build specific models including a 1D Majorana double wire and a 2D $p_{x}+ip_{y}$ topological superconductor, and obtain analytic expressions for the degenerate ground states in the presence of edges and vortices.

General construction.

Suppose we have an effective mean-field Hamiltonian in some BCS-like theory in any dimension with or without spin

[TABLE]

where $E_{p}=\sqrt{\xi_{p}^{2}+|\Delta_{p}|^{2}}$, $p$ indexes the single-particle states (including momentum, spin, or any other quantities necessary), $\bar{p}$ denotes the time reversed state of $p$, and $\alpha_{p}=u_{p}a_{p}-v_{p}a_{\bar{p}}^{\dagger}$ are Bogoliubov quasi-particle operators with $|u_{p}|^{2}+|v_{p}|^{2}=1$ and $v_{p}/u_{p}=-(E_{p}-\xi_{p})/\Delta_{p}^{*}$. The BCS-like ground state of $K_{\text{mf}}$ is (up to normalization)

[TABLE]

where the prime means each pair $p\bar{p}$ appears exactly once.

Our goal is to construct a number-conserving Hamiltonian whose ground state is $|G\rangle$. To do this, we first separate each $\alpha_{p}$ into creation and annihilation parts $\alpha_{p}\equiv C_{p}-S_{p}^{\dagger}$ with $C_{p}=u_{p}a_{p}$ and $S_{p}^{\dagger}=v_{p}a_{\bar{p}}^{\dagger}$, and define

[TABLE]

From $\alpha_{p}|G\rangle=0$, we know that $\hat{A}_{pp^{\prime}}|G\rangle=0$, and thus a parent Hamiltonian for $|G\rangle$ can be constructed by

[TABLE]

where the matrix $H_{p_{1}p_{2};p_{3}p_{4}}$ is required to be Hermitian. This construction suffices for $\ket{G}$ to be a zero-energy eigenstate of $\hat{H}$. To ensure it is a ground state, we require that the matrix $H_{p_{1}p_{2};p_{3}p_{4}}$ is positive-definite. Notice that $\hat{H}$ conserves total particle number ${\hat{N}}=\sum_{p}a_{p}^{\dagger}a_{p}$ since ${\hat{A}}_{pp^{\prime}}$ can be rewritten as

[TABLE]

which follows because $S_{p}^{\dagger}S_{p^{\prime}}^{\dagger}+S_{p^{\prime}}^{\dagger}S_{p}^{\dagger}$ vanishes by fermionic antisymmetry. Then the ground state of $\hat{H}$ with a definite particle number $N$ is simply given by the projection of $|G\rangle$ to the $N$-particle subspace $|G_{N}\rangle=\hat{P}_{N}|G\rangle$.

The double wire model.

As a first specific example of this construction, we construct one-dimensional models that reproduce the ground states of Kitaev’s 1D wire Kitaev , $\hat{H}_{\mathrm{Kitaev}}=\sum_{j}(-tc_{j}^{\dagger}c_{j+1}+\Delta c_{j}c_{j+1}+\mathrm{H.c.})-\mu\hat{N}$, where $t$, $\mu$, and $\Delta$ denote the hopping amplitude, the chemical potential and the superconducting gap, respectively. The Kitaev’s model has a special point at $\mu=0,~{}\Delta=t$ where the Hamiltonian can be rewritten as a sum of mutually commuting local operators and the spectrum is non-dispersing Kitaev . To simplify our calculation, we focus on this special point (though our construction protocol is general and could be applied to other points as well Suppl ) where the Hamiltonian has doubly degenerate ground states given by (up to normalization)

[TABLE]

where $\tilde{c}^{\dagger}_{k}=\frac{1}{\sqrt{L}}\sum^{L}_{j=1}e^{ikj}c^{\dagger}_{j}$ Suppl . The superscripts $e,o$ denote even and odd fermion parity, respectively.

We consider a double wire geometry that has two parallel one-dimensional chains, with fermion creation operators on each chain given by $a^{\dagger}_{j}\equiv c^{\dagger}_{j,1},~{}b^{\dagger}_{j}\equiv c^{\dagger}_{j,2}$, respectively. Our aim is to construct a number-conserving lattice model on these two wires whose ground states are direct products of the Kitaev ground states on each wire, projected to fixed total particle number $\ket{G_{N}}={\hat{P}}_{N}(\ket{G_{A}}\otimes\ket{G_{B}})$. We will show that the resulting Hamiltonian also leads to Majorana-like edge modes and robust ground state degeneracy.

The direct product of Kitaev ground states $|G_{A}\rangle\otimes|G_{B}\rangle$ is annihilated by Bogoliubov operators Suppl

[TABLE]

where $\sigma=1,2$, and the quasimomentum $k$ is quantized with open boundary condition $k=\frac{m\pi}{L},~{}m=1\ldots(L-1)$ 111Notice that with open boundary condition $k=m\pi/L$, the $\tilde{c}_{k\sigma},\tilde{c}^{\dagger}_{k^{\prime}\sigma^{\prime}}$ no longer satisfy the canonical anti-commutation relations, but $\{\alpha_{k\sigma},\alpha^{\dagger}_{k^{\prime}\sigma^{\prime}}\}=\delta_{kk^{\prime}}\delta_{\sigma\sigma^{\prime}}$..

Having identified the $\alpha_{k\sigma}$, and therefore the $C_{k\sigma}$ and $S_{k\sigma}^{\dagger}$, Eq. (4) gives a family of parent Hamiltonians. By choosing

[TABLE]

where $p_{j}=(k_{j},\sigma_{j})$ and $p,q,r$ are arbitrary real coefficients, one obtains Hamiltonians ${\hat{H}}$ with interactions that are local in real space Suppl ,

[TABLE]

where $J_{||,j}=b_{j}a_{j}-b_{j+1}a_{j+1},$ $J_{=,j}=a_{j+1}a_{j}-b_{j+1}b_{j},~{}J_{\times,j}=b_{j+1}a_{j}-b_{j}a_{j+1}$ and $\alpha,\beta,\gamma$ are real numbers determined by $p,q,r$ Suppl with constraint $\alpha+\beta+\gamma=0$. In order for the matrix Eq. (8) to be positive-definite, the $\alpha,\beta,\gamma$ should satisfy $\alpha<0,\beta<t,\gamma<t$, which, combined with $\alpha+\beta+\gamma=0$, gives the triangle region that is shown in Fig. 1. The center of the triangle $\alpha=-t,~{}\beta=\gamma=t/2$ reproduces the model in Ref. SDiehl .

Since $\hat{H}$ preserves total particle number $N=N_{A}+N_{B}$ and single wire fermion parity $P^{A,B}=(-1)^{\hat{N}_{A,B}}$, ground states in each $N$-particle sector are doubly degenerate. For example, if $N$ is even, we have (up to normalization)

[TABLE]

This degeneracy is topologically protected in the sense that all local perturbations in the bulk, even including the ones that violate single wire parity, take the form of an identity matrix when projected to the ground state subspace. For example, using the same arguments as in Refs. NLang ; SDiehl , we explicitly find that the energy splitting $\Delta E$ due to perturbation $a_{j}^{\dagger}b_{j}+\mathrm{H.c.}$ scales as $\Delta E\sim e^{-j/l_{0}}$ (assuming $j<L/2$) for some finite length scale $l_{0}$.

The 2D $p_{x}+ip_{y}$ model.

Majorana zero modes are expected to appear at certain boundaries and in the cores of vortices in $p_{x}+ip_{y}$ superconductors. We now show that Eq. (4) can be used to construct number-conserving parent Hamiltonians for $p_{x}+ip_{y}$ topological superconductors that share the same ground state as the mean-field model

[TABLE]

where $S$ denotes an arbitrary region in the 2D plane with complex coordinates $z=x+iy$, $\partial_{z}=(\partial_{x}-i\partial_{y})/2$, $\partial_{\bar{z}}=(\partial_{x}+i\partial_{y})/2$, $d^{2}z=dxdy$, and $\psi_{z}$ is the fermionic annihilation operator at position $z$. The term $\Delta\psi_{z}\partial_{\bar{z}}\psi_{z}+\mathrm{H.c.}$ characterizes chiral $p$-wave pairing, and $\mu\psi_{z}\psi_{z}^{\dagger}$ is the chemical potential term (which differs from the usual convention by a constant). Although in principle we can construct number-conserving parent Hamiltonians for all values of $(m,\Delta,\mu)$, in the following we only consider a special point $\mu=\frac{m\Delta^{2}}{2}$ and use natural units $2m=m\Delta=1$ for simplicity. With an integration by parts, $\hat{K}$ can be separated into a bulk Hamiltonian and a boundary term $\hat{K}=\hat{K}_{\mathrm{bulk}}+\hat{K}_{\mathrm{bound}}$, with

[TABLE]

where $\partial S$ denotes the boundary of $S$. Since $\hat{K}_{\mathrm{bulk}}$ is by construction positive-definite, the ground states of $\hat{K}$ should be annihilated by the operator $\alpha_{z}\equiv 2\partial_{\bar{z}}\psi_{z}-\psi_{z}^{\dagger}$ for all $z\in S$ in order to minimize $\hat{K}_{\mathrm{bulk}}$ (we will account for the boundary term momentarily). The ground states with even fermion parity can in general be constructed as

[TABLE]

where the two-particle wave function $g(z,z^{\prime})$ satisfies $g(z,z^{\prime})=-g(z^{\prime},z)$ and

[TABLE]

which guarantees that $\alpha_{z}|G^{e}\rangle=0$. To simultaneously minimize the boundary term $\hat{K}_{\mathrm{bound}}$, the function $g(z,z^{\prime})$ should satisfy certain boundary conditions that depend on the geometry of the region $S$, which we will discuss later.

To find a parent Hamiltonian for the mean-field ground state $|G^{e}\rangle$, we again follow our general construction given in Eqs. (4) and (5) where we identify $C_{z}=2\partial_{\bar{z}}\psi_{z}$ and $S^{\dagger}_{z}=\psi^{\dagger}_{z}$, leading to

[TABLE]

where $W(z_{1},z_{2};z_{3},z_{4})$ is a positive-definite Hermitian matrix. We further restrict ourselves to Hamiltonians describing short-ranged interactions, i.e. $W(z_{1},z_{2};z_{3},z_{4})$ tends to zero sufficiently fast when the distance between any two points $|z_{i}-z_{j}|$ becomes large. Furthermore, the boundary term $\hat{K}_{\mathrm{bound}}$ in Eq. (The 2D $p_{x}+ip_{y}$ model.) should be added into $\hat{H}_{\mathrm{bulk}}$ to uniquely pick out the same set of ground states as $\hat{K}$

[TABLE]

The new interacting Hamiltonian $\hat{H}$ harbors the topological $p_{x}+ip_{y}$ ground state. It is number-conserving because both $\hat{H}_{\mathrm{bulk}}$ and $\hat{K}_{\mathrm{bound}}$ preserve total particle number.

As a specific example, we choose $W(z_{1},z_{2};z_{3},z_{4})=w(z_{1}-z_{3})\delta^{2}(z_{1}-z_{2})\delta^{2}(z_{3}-z_{4})/4$ in Eq. (15) where $w(-z)=[w(z)]^{*}$ and $w(0)=1$, after rearranging terms we get

[TABLE]

The function $w(z-z^{\prime})$ is required to be positive-definite in the matrix sense (i.e. $\tilde{w}_{\mathbf{k}}=\int w(\mathbf{r})e^{i\mathbf{k}\cdot\mathbf{r}}d^{2}\mathbf{r}$ should be real and positive for all $\mathbf{k}$), and should decay sufficiently fast as $|z-z^{\prime}|$ becomes large.

Degenerate ground states with vortices.

It is known that vortices in the mean-field model Eq. (11) have localized Majorana zero-modes NRead ; Ivanov giving rise to topologically-protected degeneracy and non-Abelian statistics. It would be interesting to see whether these important properties survive in our number-conserving model, as these properties are crucial for the realization of topological quantum computation TPQC-RMP . As a first step, we need to obtain the ground state wavefunctions, a calculation whose results we present in this section. One remarkable feature of our model is that the analytic expressions of degenerate ground states could be exactly obtained even when there are an arbitrary number of vortices in the 2D plane. These explicit expressions could give us deeper insight into the topological properties of the ground states and provide a platform to study non-Abelian statistics in number-conserving interacting models in an exact manner.

We consider the geometry shown in Fig. 2, with $2M$ vortices lying in an unbounded 2D plane located at $\eta_{1},\eta_{2},\ldots,\eta_{2M}$, respectively. We assume that the core of each vortex is localized inside a radius $r_{0}$ much smaller than the minimal distance between any two vortices. We use the gauge convention in which the superconducting order parameter is the same everywhere [i.e. still consider the same Hamiltonians in Eqs. (11) and (15)] while fermion fields are anti-periodic around each vortex $\psi^{\dagger}_{\theta+2\pi}=-\psi^{\dagger}_{\theta}$.

In the mean-field model Eq. (11) there is a Majorana zero-mode $\gamma_{j}$ with $\gamma_{j}^{\dagger}=\gamma_{j}$ and $[\gamma_{j},\hat{K}]=0$ localized at the $j^{\mathrm{th}}$ vortex. In total we have $2M$ localized Majorana modes $\gamma_{1},\gamma_{2},\ldots,\gamma_{2M}$ which could be combined to $M$ independent fermion operators $a_{j}=(\gamma_{2j-1}+i\gamma_{2j})/2$, leading to $2^{M}$ degenerate mean-field ground states $|n_{1},n_{2},\ldots,n_{M}\rangle$ with $n_{j}=0,1,~{}1\leq j\leq M$. In our number-conserving model defined in Eq. (16), the ground states with $N$ particles are obtained by projecting the mean-field ground states to the $N$-particles sector. Only those mean-field states with fermion parity equal to $(-1)^{N}$ survive this projection, therefore we are left with $2^{M-1}$ fold degeneracy in each sector.

As an example, we consider the $2M=4$ case and assume $N$ to be even. One of the mean-field ground states $|G_{12,34}\rangle$ could be constructed by Eq. (13) with

[TABLE]

where the terms in the bracket is to guarantee that $g_{12,34}(z,z^{\prime})$ is anti-periodic around each vortex, in accordance with our gauge convention. From the identity $2\partial_{\bar{z}}\frac{1}{z-z^{\prime}}=2\pi\delta^{2}(z-z^{\prime})$ it is easy to see that $g_{12,34}(z,z^{\prime})$ satisfies Eq. (14), and it can be checked that the state $|G_{12,34}\rangle$ also minimizes $\hat{K}_{\mathrm{bound}}$ up to some small corrections 2D_TPSC_mf . Applying the projection operator $\hat{P}_{N}$, we get an $N$-particle ground state with multiparticle wave function (up to normalization)

[TABLE]

where Pf denotes the Pfaffian of the anti-symmetric $N\times N$ matrix $g_{12,34}(z_{i},z_{j})$. The form of the wave function given in Eq. (19) is very similar to one of the Moore-Read Pfaffian states with four quasiholes MR1991 ; Nayak1996 , which were constructed to describe the excitations in the $\nu=5/2$ fractional quantum hall effect (FQHE). By permuting the indices $1,2,3,4$ we get two other degenerate ground states with wave functions $\psi_{13,24}(z_{1}\ldots z_{N})$ and $\psi_{14,23}(z_{1}\ldots z_{N})$. However, using the same method in Ref. Nayak1996 we can prove that these three states are linearly dependent and the space spanned by them is actually two dimensional, consistent with our previous argument.

The non-Abelian statistics of the mean-field ground states of the $p_{x}+ip_{y}$ model have been well-studied in Ref. Ivanov . Braiding the $j^{\mathrm{th}}$ and the $(j+1)^{\mathrm{th}}$ vortices adiabatically gives rise to a unitary rotation $\hat{B}_{jj+1}=\exp(\frac{\pi}{4}\gamma_{j+1}\gamma_{j})$ on the ground state subspace. In our number-conserving model constructed in Eq. (16), the process of braiding preserves total particle number, thus we have to recalculate the Berry’s matrix for each $N$-particle sector, which will be the subject of future work. We expect that such calculation could be done using similar methods in Ref. Bonderson2011 , where unitary evolution of the Moore-Read Pfaffian states due to braiding of quasiholes are calculated with the help of the plasma analogy.

Conclusion.

We have constructed infinite families of number-conserving, interacting Hamiltonians with exact BCS-like ground states, with specific models including a 1D Majorana double wire and a 2D $p_{x}+ip_{y}$ topological superconductor. In the $p_{x}+ip_{y}$ model we obtained analytic expressions of degenerate ground states with four vortices, and pointed out their similarity to the Moore-Read Pfaffian states with four quasiholes constructed in the $\nu=5/2$ FQHE context. Our models give us a deeper theoretical understanding of topological phenomena in interacting systems, set a viable framework for building more realistic models of topological superconductors, and may provide useful guidelines for experimental realization of Majorana zero modes.

I The Double wire model

In this section we first present the detailed derivations of Eqs. (6-9) in the main text, and then we discuss some alternative derivations of the parent Hamiltonian.

I.1 Diagonalization of Kitaev’s Hamiltonian in momentum space with open boundary

To obtain the Bogoliubov operators in Eq. (7) and the form of ground states in Eq. (6) in our main text, here we present the momentum space diagonalization of Kitaev’s Hamiltonian with $\Delta=t$

[TABLE]

To this end we search for Bogoliubov eigenmodes defined as

[TABLE]

Being the eigenmodes of $H_{\mathrm{Kitaev}}$ with energy $E_{k}>0$, they satisfy $[\alpha_{k},H_{\mathrm{Kitaev}}]=E_{k}\alpha_{k}$, which gives difference equations on $u^{k}_{j},v^{k}_{j}$

[TABLE]

with boundary conditions

[TABLE]

To solve these equations, we notice that the ansatz solutions

[TABLE]

with $\lambda(k)=-i\frac{2\Delta\sin k}{E_{k}+\mu+2t\cos k}$ and $E_{k}=\sqrt{\mu^{2}+4t\mu\cos k+4t^{2}}$ satisfy Eq. (I.1). The boundary conditions in Eq. (S4) give constraints on the quasi-momentum $k$ and the real parameter $\theta_{k}$

[TABLE]

where $0<k<\pi$ and $0<\theta_{k}<\pi/2$. The Bogoliubov operators are (up to normalization)

[TABLE]

As an aside, we mention that if we replace Eq. (7) in our main text by Eq. (S7) and use the same matrix in Eq. (8), then, still following our general construction, we can get a bigger family of number-conserving, short-range interacting Hamiltonians with ground states $|G^{e,o}\rangle$ depending on $\mu/t$, and this method can be generalized to construct parent Hamiltonians for Kitaev’s ground states at arbitrary points $(t,\Delta,\mu)$ (even including points in the topologically trivial phase).

At the $\mu=0$ point we get especially simple expressions

[TABLE]

which leads to the single wire version of Bogoliubov operators in Eq. (7) after normalization. To verify that the expressions given in Eq. (6) are indeed the ground states of $H_{\mathrm{Kitaev}}$, we show that $|G^{e}\rangle$ and $|G^{o}\rangle$ are annihilated by all $\alpha_{k}$. We have

[TABLE]

where $\mathrm{sgn}(x)=x/|x|$ for $x\neq 0$ and $\mathrm{sgn}(0)=0$. It follows that

[TABLE]

leading to $\alpha_{k}|G^{e}\rangle=0$. Furthermore, it can be easily checked that $\{c^{\dagger}_{k=0},\alpha_{k}\}=0$ for all $k=m\pi/L,~{}1\leq m\leq L-1$, thus $\alpha_{k}|G^{o}\rangle=\alpha_{k}c^{\dagger}_{k=0}|G^{e}\rangle=-c^{\dagger}_{k=0}\alpha_{k}|G^{e}\rangle=0$. We conclude that Eq. (6) indeed gives us the ground states of $H_{\mathrm{Kitaev}}$ at the point $t=\Delta,~{}\mu=0$.

I.2 Detailed derivation of Eq. (9)

To verify that the combination of Eqs. (4) and (8) indeed give the local form of Eq. (9), we first notice that

[TABLE]

Using the completeness and orthonormality of $\sin kj$,

[TABLE]

we have

[TABLE]

where $C_{j\sigma}=c_{j+1,\sigma}-c_{j\sigma},~{}S_{j\sigma}^{\dagger}=-c^{\dagger}_{j+1,\sigma}-c^{\dagger}_{j\sigma}$. The parent Hamiltonian given in Eqs. (4) and (6) can then be expanded in position space (we use the shorthand $\sum_{\boldsymbol{k}}=\sum_{k_{1},k_{2},k_{3},k_{4}}$ and $\sum_{\boldsymbol{\sigma}}=\sum_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}$)

[TABLE]

where $\beta=-(q+r)/2,~{}\gamma=(q-r)/2$ and $t=(p+q-3r)/2$. By expanding the last line of Eq. (S14) we get the form of Eq. (9) in the main text (with $\alpha=r=-\beta-\gamma$).

I.3 Positive region of $\hat{H}$

We now prove that the matrix $H_{p_{1}p_{2};p_{3}p_{4}}$ given in Eq. (8) is positive semi-definite in the triangle region shown in Fig. 1 in the main text. Notice that $H_{p_{1}p_{2};p_{3}p_{4}}=H_{k_{1}k_{2};k_{3}k_{4}}\cdot\lambda_{\sigma_{1}\sigma_{2};\sigma_{3}\sigma_{4}}$ with orbital part $H_{k_{1}k_{2};k_{3}k_{4}}=\frac{16}{L^{2}}\sum^{L}_{j=1}\sin k_{1}j\sin k_{2}j\sin k_{3}j\sin k_{4}j$ and spin part $\lambda_{\sigma_{1}\sigma_{2};\sigma_{3}\sigma_{4}}=p\delta_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}+q\delta_{\sigma_{1}\sigma_{2}}\delta_{\sigma_{3}\sigma_{4}}-r(\delta_{\sigma_{1}\sigma_{3}}\delta_{\sigma_{2}\sigma_{4}}+\delta_{\sigma_{1}\sigma_{4}}\delta_{\sigma_{2}\sigma_{3}})$. The orbital part is always positive, since for any vector $f_{kk^{\prime}}$ we have

[TABLE]

where $f_{j}=\sum_{kk^{\prime}}f_{kk^{\prime}}\sin kj\sin k^{\prime}j$. For the spin part $\lambda_{\sigma_{1}\sigma_{2};\sigma_{3}\sigma_{4}}$, we write it in the matrix form (assume the order $aa,ab,ba,bb$)

[TABLE]

with $\lambda_{1}=\begin{bmatrix}-\alpha&-\alpha\\ -\alpha&-\alpha\end{bmatrix}$ acting on ($ab,ba$) and $\lambda_{2}=\begin{bmatrix}2t-\gamma-\beta&\gamma-\beta\\ -\gamma-\beta&2t-\gamma-\beta\end{bmatrix}$ acting on $(aa,bb)$. Thus $\lambda$ is positive-definite if and only if both $\lambda_{1}$ and $\lambda_{2}$ are positive-definite. The condition that $\lambda_{1}$ is positive-definite gives $-\alpha>0$, while $\lambda_{2}$ is positive-definite gives $|\gamma-\beta|<2t-\gamma-\beta$, which simplifies to $\alpha<0,~{}\beta<t,~{}\gamma<t$, leading to the triangle region in the main text.

I.4 Alternative derivations of the parent Hamiltonian

The derivation of $\hat{H}$ presented above enables us to see how the double wire model follows from our general construction and can be generalized to arbitrary points of Kitaev’s model. However, for the double wire parent Hamiltonian constructed in our main text, simpler derivations exist. Actually, Eq. (S13) gives us annihilators of the double wire ground state $|G_{A}\rangle\otimes|G_{B}\rangle$ in position space. Thus we can directly build the parent Hamiltonian $\hat{H}$ using the real space version of Eqs. (4) and (5) with $C_{j\sigma}=c_{j+1,\sigma}-c_{j\sigma},~{}S_{j\sigma}^{\dagger}=-c^{\dagger}_{j+1,\sigma}-c^{\dagger}_{j\sigma}$ (or equivalently, directly go to the last line of Eq. (S14) without working in momentum space at all), leading to the same Hamiltonian Eq. (9) in our main text. This derivation is a direct generalization of the one given in Ref. SDiehl .

Another simple derivation is based on using an alternative basis of single wire ground states (at $\Delta=t,\mu=0$) AoP2014

[TABLE]

and observing the following properties

[TABLE]

where $\bar{n}_{i}\equiv 1-n_{i}$. With this, it is easy to check that the number-conserving single wire operator $[(c^{\dagger}_{j}c_{j+1}+\mathrm{H.c.})+2(n_{j}-\frac{1}{2})(n_{j+1}-\frac{1}{2})-\frac{1}{2}]$ annihilates $|G^{\eta}\rangle$. To include interwire couplings, we notice that

[TABLE]

where $\hat{U}^{\square}_{j}=n^{a}_{j}n^{a}_{j+1}n^{b}_{j}n^{b}_{j+1},~{}~{}\hat{U}^{(3p)}_{j}=[n^{a}_{j}n^{a}_{j+1}(n^{b}_{j}+n^{b}_{j+1})+(a\leftrightarrow b)]$, and $|G^{\eta_{A},\eta_{B}}\rangle=|G^{\eta_{A}}_{A}\rangle\otimes|G^{\eta_{B}}_{B}\rangle$ is the double wire ground state constructed by direct product of single wire ground states. Therefore

[TABLE]

for $\alpha+\beta+\gamma=0$. It then follows that the Hamiltonian constructed in Eq. (9) in the main text satisfies $\hat{H}|G^{\eta_{A},\eta_{B}}\rangle=0$, i.e. $|G^{\eta_{A},\eta_{B}}\rangle$ is an eigenstate of $\hat{H}$ with zero energy. This method can only tell us that $|G^{\eta_{A},\eta_{B}}\rangle$ is an eigenstate of $\hat{H}$. To find the positive region of $\hat{H}$ (where $|G^{\eta_{A},\eta_{B}}\rangle$ become its ground states), we still have to turn to other means.

In Fig. S1 we draw a pictorial representation of the interaction terms of the parent Hamiltonian Eq. (9), including two types of nonlinear terms: interactions and correlated pair tunnelings.