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.

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. where † c c ( ) j j is fermion annihilation (creation) operator on site j, = † n c c j j 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 → φ c e c j i j . 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

Model
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 Z f 2 symmetry, which is defined as is the total fermion number, and it is obvious that . Z f 2 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 real

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 λ j 1 and λ j 2 and is given by where B is a L × L real matrix, With the help of SVD, B = UΛV T , where Λ is a real diagonal matrix, U and V are real orthogonal matrices, H μ can be diagonalized as follows, 2 are the complex fermion operators.
In the weak pairing region, μ < 2t, we find that (See Appendix A for details) the smallest singular value Λ k is nonzero given by k L 0 and the corresponding matrix elements (1 ) k vL v 2 1/2 0 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 = μ L = μ/2, resulting in Λ ko = 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 is Λ k0 given in Eq. 10, approaches to zero with the exponential factor e −Lln(2t/μ) 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 read , namely, it coincides to its antiparticle. Using Eq. A2a, we have 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 jl Hermitian. Especially, the edge component of G jl is given when j = 1 and l = L, Note that the correlation function G jl 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 1L in the thermodynamic limit reflects the existence of edge modes.
The edge correlation function G 1L is easy to calculate in the case of Δ = t and U = 0, and is given for the ground state |0〉 by As proved by Lieb et al. 63 , this summation is of order of O(1/L). When μ < 2t, The nonvanishing value of G 1L 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 1L for the topological degenerate state |1〉. Thus, for a generic ground state GS , the edge correlation function in the thermodynamic limit is given by L L 1 2 Note that the nonzero contribution U V k Lk 1 0 0 comes from the Majorana zero mode k 0 . Other modes mainly distribute in the bulk and the contributions to G 1L 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 1L . In the quantum critical region, L z 1 with critical exponent z = 1. Now we would like to examine the behavior of G ij 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 1L ≠ ±G L1 . 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.
Therefore, we propose to use the edge correlation function G 1L 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 Here the normalization factors are given by The wave vector k I 's are given by the following equation, sin , I I and k II 's are determined by II 0 with v determined by (32) For this k II 0 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 k II 0 mode vanishes in the thermodynamic limit, Λ = .

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 ij . 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 .

Phase diagrams.
The TSC phase is detected by the two-fold degenerate ground states with opposite fermion number parity Z f 2 and CAT phase is the two-fold degenerate ground states with opposite particle-hole symmetry Z p 2 . In contrast, the two ground states of CDW and ICDW phase have the same Z f 2 . In practice, we compute the matrix elements for Z f 2 or Z p 2 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 ij . We also compute the single-particle correlation function G ij 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 ij serves an efficient measurement for edge states and thereby the topological order.
Edge correlation function G 1L . The nonvanishing edge correlation function G 1L characterizes the topological order. We fix Δ = t and study G 1L as a function of μ and U. The result is plotted in Fig. 7. The value of G 1L 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.
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 〈 | | 〉 n 0 0 j 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.

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 ij and find that there exists a long-range edge correlation G 1L 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 1L 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 1L ∝ w z with z = 1 near the critical point, where w = Δ, μ c − μ, etc. is a control parameter that drives the system from a topologically nontrivial phase to a topologically trivial phase.

Appendix A Exact diagonalization of non-interacting Kitaev chains with Δ = t
In this appendix, we provide details in exact diagonalization of the matrix B in Eq. 8. We write the matrix B in the SVD form 33 , T where the matrix Λ = Λ k is diagonal. The matrices U and V are orthogonal transformations The energy spectra of the Hamiltonian H are given by the singular values of the matrix B. We note the orthogonal matrices U and V diagonalize BB T and B B T , respectively The singular values Λ k are the non-negative square roots of the eigenvalues of BB T . Similar diagonalization was found by Lieb et al. in the study of Heisenberg-Ising model 63 . The orthogonal matrices U and V are found to be The graphical solution is shown in the Fig. A1. For μ ≥ 2t, there are L real roots, including all the normal modes. For μ < 2t, there are L − 1 real roots and one complex root Then for this special mode we have