Majorana zero modes and long range edge correlation in interacting Kitaev chains: analytic solutions and density-matrix-renormalization-group study.

We study Kitaev model in one-dimension with open boundary condition by using exact analytic methods for non-interacting system at zero chemical potential as well as in the symmetric case of Δ = t, and by using density-matrix-renormalization-group method for interacting system with nearest neighbor repulsion interaction. We suggest and examine an edge correlation function of Majorana fermions to characterize the long range order in the topological superconducting states and study the phase diagram of the interating Kitaev chain.

Introduction

Majorana[1] zero mode (MZM) has attracted a lot of attention in the recent years[2-4], which may emerge as a novel excitation in some topological condensed matter systems. MZMs obey non-Abelian statistics and have potential application to build robust qubits against decoherence in quantum computation[5,6]. The emergence of MZMs has been theoretically proposed in a number of condensed matter systems, including chiral p-wave superconductors[7,8], ν = 5/2 fractional quantum Hall system[9], the interface between a topological insulator and an s-wave superconductor[10], proximity-induced superconductor for spin-orbit coupled nanowires[11,12], spin-orbit coupled semiconductor with externally applied Zeeman field[13-15], and ferromagnetic atoms in proximity to superconductors[16,17]. There also exist various experimental efforts to realize and detect MZMs in these proposed systems[18-27].

Among these candidates, the one-dimensional (1D) systems are of special theoretical interest for possible generalization to interacting systems. The interaction may change properties drastically in 1D systems. The Fermi liquid description of the interacting Fermi gas usually works in 2D or 3D. However, it breaks down in 1D and the systems become Luttinger liquids. Fortunately, there have been a number of many-body techniques suitable to study various 1D problems[28], which make the generalization of the MZMs in 1D models accessible. On the other hand, the interaction will modify topological systems violently, e.g. the non-interacting classification of fermionic systems[29-31] will “collapse” and there exists a continuous path connecting trivial and topological phases in 1D[32].

Kitaev chain[7] is a prototype of 1D systems possessing MZMs at the two edges. The non-interacting Kitaev model was initially solved in a ring with periodic boundary condition. The edge state was then proposed to exhibit MZM. The model has been generalized to interacting case with nearest neighboring repulsive interaction. The interacting Kitaev model does not have analytic solutions in general cases except for a set of specially tuned parameters[33,34]. The model can also be studied by numerical methods[34-36]. In general, interacting effects on MZMs have been investigated in various systems, e.g. nanowires[35,37-40], multiband nanowires[41], helical liquids[42], two-leg ladders[43], Josephson junctions[44], Abrikosov vortex lattice[45] and topological insulator/superconductor heterostructure[46]. The interplay of disorder and interaction has also been analyzed[47,48]. The MZM is stable against weak perturbations including the interaction and disorder. However, the generic interaction effect remains an open question, although lots of efforts have been made, which includes the exact solution[49], topological classification[32,50], entanglement entropy investigation[51], many-body MZM operator[52,53], super-symmetry approaches[34,54-56] and parafermion edge zero mode[57-61].

In this paper, we shall first study non-interacting Kitaev chain of length L with open boundary condition by using an analytic method, which is accessible at the symmetric points with zero chemical potential and equal pairing and hopping amplitudes, Δ = t. We propose a correlation function of the two Majorana operators as a long range order parameter to describe non-trivial topological state with edge MZMs and calculate the long range correlation function explicitly. We then study Kitaev model with nearst neighboring repulsion interaction in open boundary condition by using density matrix renormalization group (DMRG) method. We show that the qualitative feature of the long range correlation remain unchanged in the interacting systems provided that the system is in the topological non-trivial phase. The phase diagram in the interacting model will also be discussed. This work is a generalization of our previous work on exact solution for interacting Kitaev chain at symmetric point[62]. The exact solution can be obtained only at special point and we have to resort to numerical methods for generic parameters. In this paper we explore the phase diagram with a generic chemical potential μ, and demonstrate that the edge correlation is not only valid in the non-interacting system but also in more generic interacting systems.

This paper is organized as follows. In Section 2, the model Hamiltonian are presented and Majorana fermion representation is introduced. In Section 3, we study non-interacting models by using analytic solutions. A single-particle correlation function is introduced and its edge component is used to describe the topological order. In Section 4, numerical DMRG analysis is carried out to study interacting systems. Section 5 is devoted to discussions.

Model

Without loss of generality, we consider a chain of spinless fermions with open boundary condition. The Hamiltonian of such an interacting Kitaev chain iswhere is fermion annihilation (creation) operator on site j, is the fermion number operator, t is the hopping matrix element, and Δ is the p-wave superconducting pairing potential induced by the proximity effect, μ is the chemical potential controlling the electron density, and U is the nearest neighbor interaction. One can always choose Δ real and non-negative by the global transformation . Similarly, one can study the case of t ≥ 0 and μ ≥ 0 only, since the parameter transformations t → −t and μ → −μ can be realized by by the gauge transformation and particle-hole cojugation respectively. Note that all these transformations will keep other parameters unchanged. In this paper, we only consider repulsive nearest neighbor interaction with U ≥ 0. When U = 0, this model will reduce to the usual (non-interacting) Kitaev chain[7].

The Hamiltonian has the fermion number parity symmetry, which is defined aswhere is the total fermion number, and it is obvious that and . conserves in the whole parameter space. In the presence of the pairing potential Δ, the total fermion number is not conserved but only conserved modulo 2.

Majorana fermion representation

We shall use the Majorana fermion representation to investigate the interacting Kitaev chain. Following Katsura et al.[33], we split one complex fermion operator into two Majorana fermion operators

The Majorana fermion operators are realand satisfy the anticommutation relationswhere a, b = 1, 2. In the Majorana fermion representation, the Hamiltonian of the interacting Kitaev chain becomes

Non-interacting Kitaev chains

In this section, we consider the non-interacting Kitaev chains with open boundary condition and discuss the relations among the topological degeneracy, the Majorana zero mode, and the edge correlation functions. We shall use analytic method to exactly solve the two non-interacting cases with Δ = t, U = 0 and μ = 0, U = 0 by the singular value decomposition (SVD) in Majorana fermion representation.

Non-interacting chains with Δ = t

In this case, the transition between the topological superconductor and the trivial superconductor can be studied by tuning the chemical potential μ. The non-interacting Hamiltonian H is quadratic in and and is given bywhere B is a L × L real matrix,

With the help of SVD, B = UΛV , where Λ is a real diagonal matrix, U and V are real orthogonal matrices, H can be diagonalized as follows,where Λ ≥ 0 are singular values of the matrix B, and are the complex fermion operators.

In the weak pairing region, μ < 2t, we find that (See Appendix A for details) the smallest singular value Λ is nonzero given byand the corresponding matrix elements where is the normalization factor, and v is a positive real number determined by Eq. A12.

It is worth noting that a similar model has been solved by Katsura et al.[33] using SVD. In their case, the chemical potential is half of the bulk’s value at edge, μ 1 = μ  = μ/2, resulting in Λ = 0.

Topological degeneracy and the edge mode

It is well known that there exist two topologically distinct phases in the non-interacting Kitaev chain model[7,64,65]. For strong pairing μ > 2t, the system is in the trivial superconducting state, while for weak pairing μ < 2t, the system is in the topological superconducting state.

In the trivial superconducting state, the energy spectrum is gapped and the ground state is non-degenerate. However, in the topological superconductor, the energy gap between the ground state |0〉 and the first excited state is Λ given in Eq. 10, approaches to zero with the exponential factor e − in the large L limit. Thus, the k 0-mode is a zero mode and the topological superconductor has two-fold degenerate ground states in thermodynamic limit. In other words, it is a gapped system with two-fold topological degeneracy.

Now we shall check that the first excited state |1〉 is an edge mode. It is a single particle (hole) excited state. The particle and hole parts of the wavefunction readandrespectively, where Eqs 11 and A2a have been used in the derivation. It is easy to see that this zero mode has a complex wave vector k 0 = π + iv and the wavefunction is well localized at edges with localization length v −1 as demonstrated in Fig. 1.

Figure 1

The particle wavefunction for the k 0-mode with L = 500 and v = 0.2.

Now we would like to examine that the k 0 mode is indeed a Majorana mode, say, , namely, it coincides to its antiparticle. Using Eq. A2a, we have

By Eq. 11, we find that

So that there exists one Majorana mode with at the edge j = 1 and another Majorana mode with at the edge j = L.

Fermion number parity and edge correlation function

There are two characterizing features for topological ordered systems, (base-manifold dependent) ground state degeneracy and gapless edge states.

We note the ground state |0〉 and the excited state |1〉 have opposite fermion number parity

In the thermodynamic limit, the first excited |1〉 is degenerate with the ground state |0〉.

We define the following single-particle correlation function at two sites j and l,where the imaginary i is introduced to make G Hermitian. Especially, the edge component of G is given when j = 1 and l = L,

Note that the correlation function G is a block of single-particle(hole) density of matrix, which can be generalized to interacting systems and reflects the site-distribution of single-particle component in a many-particle wavefunction. As long as the bulk is uniform, the finite value of G 1 in the thermodynamic limit reflects the existence of edge modes.

The edge correlation function G 1 is easy to calculate in the case of Δ = t and U = 0, and is given for the ground state |0〉 by

When μ ≥ 2t,

As proved by Lieb et al.[63], this summation is of order of O(1/L). When μ < 2t,

The nonvanishing value of G 1 for μ < 2t in the thermodynamic limit reflects the topological order in the topological superconductor state. In this topological phase, we can also calculate edge correlation function G 1 for the topological degenerate state |1〉.

Thus, for a generic ground state , the edge correlation function in the thermodynamic limit is given by

Note that the nonzero contribution comes from the Majorana zero mode k 0. Other modes mainly distribute in the bulk and the contributions to G 1 is of order of O(1/L), which is neglectable in the thermodynamic limit. At the quantum critical point μ = 2t, we have v = 0 and the wave vector of the Majorana zero mode becomes real k 0 = π. The k 0-mode is no longer localized at edges but merges into the bulk, resulting in vanishing edge correlation function G 1. In the quantum critical region,with critical exponent z = 1.

Now we would like to examine the behavior of G inside the bulk, which can be done numerically. Two topologically distinct examples are investigated and shown in Figs 2 and 3 respectively. The first example is given by Δ = t, μ = 3t, U = 0, which is in the topologically trivial phase, where a peak appears at short range with i~j while long range correlation is absent. The second example is given by Δ = t, μ = t, U = 0, which is in the nontrivial topological superconductor phase. There exhibits a long range peak at i = 1 and j = L, and long range correlation is still absent inside the bulk. We note the edge correlation is not symmetric or antisymmetric, i.e. G 1 ≠ ±G . Hence there is no peak at i = L and j = 1. If we use parameters with t < 0, the peak will appear at i = L and j = 1. So it is a matter of choice. The point is there is a edge correlation function corresponding to the Majorana zero mode.

Figure 2

Correlation function |G | for a topologically trivial state, Δ = t, μ = 3t, U = 0.

Figure 3

Correlation function |G | for a topologically nontrivial state, Δ = t, μ = 0, U = 0.

Therefore, we propose to use the edge correlation function G 1 to characterize the topological order and emerged edge states. We shall examine this for the non-interacting systems with different parameters in the next subsection and for the interacting systems in the next section.

Non-interacting chains with μ = 0

In this subsection, we utilize non-interacting Kitaev chains with μ = 0 to study how topological order will vanish as the superconducting gap Δ approaches zero. The Hamiltonian now reads

We are able to diagonalize the Hamiltonian H Δ by SVD as before. There exist two kinds of modes in this situation. For the first kind of modes, the two orthogonal matrices U and V are found to be

The second kind of modes is given by

Here the normalization factors are given byand

Corresponding singular values are given by

The wave vector k ’s are given by the following equation,and k ’s are determined by

Besides L − 1 real k ’s, there exists a single complex k in the second kind modes,with v determined by

For this mode we have

Then the normalization factor can be written explicitly,and the singular value reads

It is easy to see that the singular value of mode vanishes in the thermodynamic limit,

The (single particle) wavefunction of this zero mode is given bywhich has nonzero value only near the edge in the thermodynamic limit. Similarly, one can verify that at edges. Hence the -mode is the Majorana zero mode localized at edges. When Δ → 0, the wave vector of the zero mode becomes real and the Majorana zero mode is no longer localized at edges. This is consistent with the condition for the boundary Majorana fermion argued by Kitaev[7], i.e. the presence of an arbitrary small superconducting gap Δ.

Now we compute the edge correlation function G 1 for the ground state |0〉,and for the topological degenerate state ,

For small but finite Δ, we havewith critical exponent z = 1. Thus the edge correlation function vanished as Δ → 0.

Interacting Kitaev chains: DMRG analysis

In this section, we shall study interacting Kitaev chains by carrying out DMRG calculations in the language of matrix product states[66] with various model parameters in Hamiltonian 1 and system size up to L = 140. We compute the energy of low lying states, local particle density, as well as the single-particle correlation function G .

Phase diagrams

Figure 4 displays the phase diagram at Δ = t obtained from the combination of exact solutions and DMRG calculations. As a function of μ and U, there are five distinct phases, trivial superconductor (SC), topological superconductor (TSC), commensurate charge density wave (CDW), incommensurate charge density wave (ICDW) and Shrödinger-cat-like state (CAT). The five different phases are separated from each other by critical lines. Such a phase diagram is consistent with previous studies[33-35] except the CAT states at μ = 0 obtained by exact solution[62].

Figure 4

Phase diagram for the interacting Kitaev chain with Δ = t. SC stands for trivial superconductor, TSC stands for topological superconductor, CDW stands for charge density wave, ICDW stands for incommensurate charge density wave and CAT stands for Shrödinger-cat-like (CAT) state. Data points are obtained within DMRG for different system sizes. Rhombuses denote SC states, circles denote TSC states, up-triangles denote ICDW states, down-triangles denote CDW states and squares denote CAT states.

The TSC phase is detected by the two-fold degenerate ground states with opposite fermion number parity and CAT phase is the two-fold degenerate ground states with opposite particle-hole symmetry . In contrast, the two ground states of CDW and ICDW phase have the same . In practice, we compute the matrix elements for or in the subspace spanned by the two lowest lying states, |0〉 and |1〉, and diagonalize the 2 × 2 matrix to obtain two eigenvalues. The distinction between ICDW and CDW can be made through local particle density and its Fourier transformation. For a CDW state, there exists a single peak at Q = π, while for a ICDW state, there appear two peaks in the Fourier spectrum.

When μ = 0, as U increases, the ground state changes from CAT to TSC and to CDW directly via the critical point U = ± t. When μ > 0, as U increase, the ground state changes from SC to TSC, ICDW and to CDW in the large U limit.

Single-particle correlation function G

We also compute the single-particle correlation function G defined in Eq. 16 for ground states. Similar to exactly solvable systems shown in Figs 2 and 3, long range correlation is absent inside the bulk. When the system is in the TSC phase, there exists a single long range peak at i = 1 and j = L. Figures 5 and 6 demonstrate two TSC states with Δ = t, μ = 0, U = 0.5 t and Δ = t, μ = t, U = 0.5 t respectively. So that G serves an efficient measurement for edge states and thereby the topological order.

Figure 5

Single-particle correlation function G for the TSC ground state with Δ = t, μ = 0 and U = 0.5t. The system size is L = 100.

Figure 6

Single-particle correlation function G for the TSC ground state with Δ = t, μ = t and U = 0.5t. The system size is L = 100.

Edge correlation function G1

The nonvanishing edge correlation function G 1 characterizes the topological order. We fix Δ = t and study G 1 as a function of μ and U. The result is plotted in Fig. 7. The value of G 1 is finite in TSC phase and vanishes in other topologically trivial phases. Thus this order parameter is valid both in the non-interacting and interacting systems to study the topological order.

Figure 7

Ground state edge correlation function G 1 as function of μ and U. Δ = t and the system size is L = 140. Squares denote SC states, circles denote TSC states, up-triagnles denote ICDW states, and down-triangles denote CDW states.

Local density of states

We can distinguish the ICDW and CDW phases by observing their local density distribution and corresponding Fourier spectrum. When the ground state is a CDW, its Fourier spectrum will have a single peak at Q = π; while for a ICDW state there are two peaks.

For various model parameters, we use the DMRG method to obtain the ground state |0〉 and local density for each site j. The Fourier spectrum is obtained by taking fast Fourier transformation of the local density distribution, whose average value has been subtracted. Here we show two typical figures of ICDW and CDW in Fig. 8.

Figure 8

Local density distribution and density spectrum. In up figure the local density of ICDW oscillates nonuniformly and its Fourier spectrum has two peaks near Q = π. In bottom figure the local density of CDW forms a bipartite lattice and the Fourier spectrum has single peaks at Q = π.

Conclusion

In summary, we have studied in this paper the Kitaev chains with open boundary condition by using analytic exact solution method for the non-interacting model and by using DMRG method for the interacting model.

We study a locally defined single-particle correlation function G and find that there exists a long-range edge correlation G 1 in the topologically nontrivial phase which is absent in topologically trival phases, while long range correlation is always absent inside bulk for all the phases. Thus, we propose that G 1 can be used to characterize the topological order in 1 + 1D fermionic systems and use it to describe quantum phase transitions between topologically trivial and nontrivial phases. It is found that G 1 ∝ w with z = 1 near the critical point, where w = Δ, μ − μ, etc. is a control parameter that drives the system from a topologically nontrivial phase to a topologically trivial phase.

Note added

After this and related paper[62] were posted (arXiv:1608.08382 and arXiv:1610.04485), there appeared several followed works, where the extension to disordered[67] and dimerized systems[68-70] were studied.