End-to-end correlations in the Kitaev chain

The interdependence between long range correlations and topological signatures in fermionic arrays is examined. End-to-end correlations, in particular classical correlations, maintain a characteristic pattern in the presence of delocalized excitations and this behavior can be used as an operational criterion to identify Majorana fermions in one-dimensional systems. The study discusses how to obtain the chain eigenstates in tensor-state representation together with the proposed assessment of correlations. Outstandingly, the final result can be written as a simple analytical expression that underlines the link with the system's topological phases.


I. INTRODUCTION
Majorana fermions [1,2] are described by a highly versatile formalism that provides conceptual-as well as technical-tools, so much so that Majorana particles can be found in solid state systems incarnated as collective phenomena. One particular scenario where this identification takes place is the Kitaev chain [3], a onedimensional fermionic system that displays a fundamental relation between Majorana excitations and topology, allowing the use of geometrical arguments to establish a connection between the presence of Majorana fermions in the open chain and the parity of the periodic chain, providing in this way simple and operational conditions to allocate Majorana particles, specifically, superconductivity, which ensures a gap in the bulk, and an odd number of Dirac points in the band structure of the noninteracting system. This approach has been rather successful in giving a qualitative characterization of the Majorana chain, and the discovery has been attracting a lot of attention over the past decades due to its potential applications in quantum information, prompting experimental verification in state-of-the-art setups, usually in the form of zero-bias conductance peaks on the edges of one-dimensional structures, as for instance in [4][5][6] to mention only the most recent studies.
Even though the Kitaev chain is well understood in terms of its topological structure, a more quantitative description in terms of the state's mean values is clearly of interest. Such a description is necessary to thoroughly characterize the behavior of the state's observables in the topological phase. This characterization can then be used on other systems where lack of integrability does not allow a direct identification of Majorana excitations [7]. The fundamental observation is that since the excitations supported by Majorana fermions are highly delocalized, it is reasonable to expect they enhance the correlations between the end sites of the chain. If that is indeed the case, experimental verification could be improved if it were possible to simultaneously test electron density on both of the chain ends. The study of edge correlations in fermion chains has been addressed in references [8,9] for specific cases that permit analytical progress [10,11].
Here the analysis covers the whole range of parameters and is conceptually exact, albeit with a numerical component. This approach allows to find a generalized expression for the correlations that complements the results derived analytically.
The fact that the Kitaev chain is an integrable model does not preclude the need for numerical analysis. Complications arise because the chain's Hilbert space grows as 2 N , being N the number of sites, and many features manifest exclusively in the thermodynamic limit. These complications can be circumvented by the use of tensor state techniques. Such techniques can be implemented in different ways. One such way is Density Matrix Renormalization Group (DMRG) [7,8,12], which minimizes the energy over the space of eigenstates of the chain's local density matrices. Another approach is to use a tensor representation as a variational network in an abstract way. This is characteristic of the method known as Matrix Product States (MPS). The way tensor state techniques are applied here is different from these two, and is more in accordance with the updating protocols introduced in reference [13]. A family of methods based on such protocols is known as Time Evolving Block Decimation (TEBD). However, the path followed in this report differs from this denomination. First, time-evolving-or step-integration is not incorporated, and second, the implementation is exact, so that numerical approximations like splittings of operator exponentials, which are huge error-contributors to TEBD, are not employed whatsoever. The techniques applied here to fermion chains have been first developed in the context of bosonic arrays in [14,15], although with some key differences, the most important one being the inclusion of pairing terms integrally in the current formalism, which is possible thanks to the decomposition of fermionic operators in terms of Majorana operators. Another aspect that contrasts with other works is that the tensor state formulation is carried completely in the fermionic Fock-space, without the extra work of incorporating the so called Wigner-Jordan transformation to reformulate the problem in terms of a spin chain, as seems to be frequent in DMRG applications to fermion systems. The Kitaev chain, also known as the Majorana chain, is described by the following Hamiltonian [3] Constant w is the next-neighbor hopping intensity, while µ is the chemical potential, which relates to the total number of fermions in the wire. Parameter ∆ is the intensity of the pairing and is known as the superconducting gap. Creation and annihilation operators follow fermionic anticommutation rules {ĉ j ,ĉ k } = 0 and {ĉ j ,ĉ † k } = δ k j . Open or periodic boundary conditions are enforced by takingĉ N +1 = 0 orĉ N +1 =ĉ 1 respectively. The model describes a quantum wire in proximity to the surface of a p-wave superconductor [1][2][3], and also a spin-polarized 1-D superconductor, since only one spin component is being considered and the pairing term mixes modes with opposite crystal momentum, as can be shown by switching to a momentum basis. Independently of boundary conditions, Hamiltonian (1) commutes with the parity operatorΠ The symmetry associated to this operator is not spatial, instead, it is a parity associated to number of particles.
Let us now introduce the Majorana operators (MOs) corresponding to site j: A key feature of these MOs is that they are hermitian, γ † k =γ k . It can be shown that the anticommutation relations are given by {γ k ,γ j } = 2δ j k , so thatγ 2 k = 1. Since there are two Majorana modes for every site, the total number of modes doubles. Equations (3) and (4) can be inverted and the result can be used to write Hamiltonian (1) asĤ When w = ∆ = 0, this Hamiltonian becomes diagonal in the Fock basis with a non-degenerate spectrum and a ground state that depends on the sign of µ. From (5) it can be seen that such a particular case corresponds to a chain where MOs from the same site pair up. This behavior is the generic signature of the trivial phase. In contrast, when w = |∆| > 0 and µ = 0, the pairing takes place between Majorana operators from neighbor sites, as can be seen in the transformed Hamiltonian A key feature of this expression is that it lacks bothγ 1 andγ 2N , which are unpaired. From these one can build an uncoupled mode, satisfying {f N ,f † N } = 1. This is a highly delocalized fermionic mode, having equitable contributions on the edges. The simplest physical operator that can be in this way built isf † Nf N , which consequently commutes with the Hamiltonian. The Hilbert space associated to such a term contains two states, one occupied and one empty. As the Hamiltonian does not include terms that could operate on neither of these, it is energetically equivalent to have many-body eigenstates with or without a particle on the aforementioned mode, i.e., the energy cost associated to this mode is zero, being that the reason whyf N is known as a Zero Mode whileγ 1 andγ 2N are known as Majorana Zero Modes (MZMs) [1] or Edge Modes. As a consequence, the whole spectrum of (6) becomes twofold degenerate. The fact that neitherγ 1 norγ 2N appear in the Hamiltonian implies that they do not pick up oscillatory phases in the Heisenberg picture, which makes them robust against this kind of decoherence mechanism. From the previous arguments it can be deduced that (6) commutes with the symmetry operator It can be noticed thatQ R is unitary and its eigenvalues are 1 for a filled mode and −1 for an unfilled one. Unlikê Π,Q R determines a symmetry only for a specific set of parameters. In spite of the spectrum being degenerate in this case, it is possible to build ground states |ψ G of (6) that are also eigenstates of (8), so that Performing the same calculation with any normalized state that is not an eigenstate ofQ R would result in a lower value, so that the maximum is linked to eigenstates of Hamiltonians having unpaired MOs completely localized at the ends. A totally analogous case would be obtained if focus is made on the point w = −|∆| < 0 and µ = 0, yieldingĤ This time it isγ 2 andγ 2N −1 which do not appear in the Hamiltonian, thus giving rise to the symmetry operator which follows a relation analogous to (9). According to reference [3], there are unpaired MOs, or Majorana fermions, over the region in parameter space surrounding the particular cases studied above as long as the gap of the equivalent periodic chain does not vanish. Following this argument it can be shown that unpaired MOs prevail in the region defined by 2|w| < |µ|. In the general case such operators only become completely unpaired in the thermodynamic limit, although with exponential convergence, and they are not completely, but still highly, localized at the ends.
The above observations suggest that in order to scan for unpaired MOs it is useful to exploit their relation with the state's local-symmetries. Let us therefore define the operatorQ It is reasonable to expect that the mean values ofQ, which actually measure end-to-end single-particle hopping, determine the degree of localization of MOs at the chain ends, taking extreme values when they are completely localized and vanishing when there is none. In order to test this conjecture the following measure is proposed The numerical calculation of this expression is not always efficient because long-range correlations are involved, hence a practical approach is desirable. Next section focuses on presenting a way in which the eigenstates of (1) as well as Z can be effectively calculated.

II. REDUCTION OF THE KITAEV CHAIN BY A SERIES OF UNITARY TRANSFORMATIONS
Hamiltonian (5) has the following general structurê where the coefficients A kl form a real antisymmetric matrix A kl = −A lk . Following Kitaev [3],Ĥ can be diagonalized by an unitary transformation that reduces it toĤ Theζ's are MO that can be expressed as linear combi- MatrixŴ is such that it transformÂ in the following The single-body energies ǫ k are real-and-positive whileŴ is a real orthonormal matrix satisfyingŴŴ † = I. This factorization is a particular case of a procedure known as Schur decomposition [16]. The diagonalized Hamiltonian can be written in terms of standard fermionic modeŝ This can be checked by replacingζ 2k−1 =f k +f † k and (15). The system eigenenergies are given by in such a way that n k = {0, 1}. The system's ground state corresponds to the case where all n k = 0, therefore If there are one or more single-body vanishing-energies ǫ k = 0, the spectrum becomes degenerate, because there is no energy difference between occupied and unoccupied zero modes.
In order to obtain the eigenstates an approach similar to that in reference [17] is adopted. First, Hamiltonian (18) is written aŝ The W j,k 's are the components of matrixŴ , as defined by (17). An unitary transformation acting on two consecutive MOs can be defined aŝ The effect of this transformation on the Hamiltonian is calculated throughĤ . This can be performed by applying the same transformation on every operator composingĤ. The operation over a linear combination of consecutive MOs is written aŝ where W ′ j = W j cos θ − W j−1 sin θ and W ′ j−1 = W j−1 cos θ+W j sin θ. The contribution ofγ j can be taken out by choosing in such a way that W ′ j = 0. If W j−1 = 0 and W j = 0, then θ = π 2 . If both W j−1 and W j are zero, then θ = 0 is enforced. In any other case the angle is well defined because the W 's are real. It is practical to choose the angle θ so that sign(sin θ) = sign(W j ) and sign(cos θ) = sign(W j−1 ). In this way W ′ j−1 = W 2 j−1 + W 2 j , leaving a positive coefficient. In order to highlight the dependence of θ with respect to W j and W j−1 , θ j is used from now on. When this operation is applied on the whole Hamiltonian, the mechanism can be described as an overall action on the diagonal modes: As a result,γ 2N vanishes fromζ 1 and the vertically aligned coefficients are in some way affected. The process continues by applying another transformation aimed at canceling the next component, which generates a similar effect on the stack of coefficientŝ The process is repeated, removing one component in each step and advancing towardγ 1Û The last transformation eliminatesγ 2 and leaves onlyγ 1 multiplied by W ′ 1,1 = W 2 1,1 + W 2 1,2 + ... + W 2 1,2N = 1. Because all the transformations are unitary, the orthogonality of the coefficients must be preserved. Hence, if onlyγ 1 remains in the top row, there cannot beγ 1 -terms in the rest of the stack. The last operation then leaves the following arrangementγ A similar series of operations can be devised to reduce the second row, this time avoiding any transformation involvingγ 1 in order to keep the first mode folded. The process can be repeated with the same intended effect in each step. However, the last transformation brings up an additional issue. Let us notice that before the last operation the stack of components looks likê is aimed at folding the antepenultimate row, however, due to the orthonormality of the original modes, it folds the last row too. However, there is no guarantee that the resulting coefficient is positive, since the transformation only takes care of the sign of the coefficients of the row being folded. Consequently, after the folding is finished, there are two possible states of the stackγ In the first case, when all the coefficients are positive, the transformed Hamiltonian becomes The eigenstates of this Hamiltonian can be identified as occupation states, and the corresponding eigenergies are given by (19). The vacuum |0 , or state without fermions, is simultaneously the system's ground state.
With regard to the second case, let us first point out that fermionic operators are given in terms of MO by the relation It can be seen that a negative sign inγ 2N induces an particle-hole transformation,ĉ † N ⇔ĉ N , in such a way that after passing to the fermionic basis the reduced Hamiltonian becomes The eigenstates of this reduced Hamiltonian can be built as in (25), but taking into account that the ground state is not the vacuum but the state with one fermion in the N -th site Likewise, the expression for the associated eigenenergy is (19). To obtain the eigenstates of the original Kitaev chain, |ψ l , one applies the transformations in reverse order over the states (25) or (27), depending on the result of the folding. The operation can be written as Both |ψ l and |ϕ l are eigenstates corresponding to E l , because unitary transformations do not change eigenvalues. Using (3) and (4) it can be shown that the transformations that appear in (28) are given bŷ andÛ Transformations with j even operate only on site j 2 and in matrix form they can be written aŝ This matrix is written with respect to occupation states with the order |0 , |1 . Transformations with j odd operate non-trivially on consecutive sites j−1 2 and j+1 2 through the following matrix representation In this case the basis order is |00 , |01 , |10 , |11 (the first position for site j−1 2 and the second for site j+1 2 ). The reducibility of the matrix underlines the fact that the Hamiltonian commutes with the parity operator (2) and therefore the eigenvectors inhabit spaces with even or odd parity. Because these matrices only mix states with the same parity, |ψ l and |ϕ l have equal parity.
In order to obtain the eigenstates of the Kitaev chain, first matrixŴ is gotten using standard numerical routines. The entries of this matrix are then used to get the folding angles θ j,k from Eq. (23). These angles are then employed to build the transformations composing expression (28) in matricial form. Since such transformations involve neighbor sites only, expression (28) can be computed using the updating protocols described in [13], so that the final result is expressed in tensorial representation. A detailed description of how to incorporate tensor product tecniques particularly on this problem is given in appendix A.

III. RESULTS
Before addressing the study of correlations in the open chain, the numerical approach proposed in the previous section is tested using the spectrum of the periodic chain. The single-body energies are given by Such energies are numerically calculated as a by-product of the Schur decomposition in (17) and then compared against these analytical results. Next, Eq. (19) is used to find the ground state energy. The tensorial representation of the ground state is then obtained following the folding protocol discussed before. The energy of such a ground state is calculated from this tensorial representation. This can be done as N times the energy of two consecutive sites, but only for eigenstates with translational symmetry, which is the case as long as such eigenstates are nondegenerate. This result is then compared with the quantity obtained before as the sum of single-body-energies. The top plot in figure 1 shows the absolute difference between the two estimates. The bottom plot depicts the mean number of particles, showing a behavior consistent with the contribution of a zero mode at |µ| = 2|w|.
Having verified the technique, let us now study correlations in the open chain. If the ground state is nondegenerate, Z in (13) can be determined from whereρ 1N is the reduced density matrix of the chain ends. The computation of such a matrix from a state written in tensorial representation is described in appendix B. MatrixQ comes given bŷ where P is the ground state's parity. Z is found as the saturation value of (34) with respect to N . The results are shown in 3.5x10 with total localization of MOs and vanishing values characterize the trivial phase. One would expect that the non vanishing values of Z provide an assessment of the level of localization of Majorana fermions at the edges. Interestingly, the numerical values taken by Z in table I correspond to rational numbers, which allows to fit the data to the following analytical function The fact that correlations between the edge sites are conditioned by the existence of unpaired MOs is readily noticeable in this elementary formula. As can be appreciated, the relation is given in terms of simple algebraic functions, so that power law coefficients are rational.  correspond to the cases i: ∆ = w and ii: µ = 0 [8]. Both instances display structural coincidence with Z in spite of differences with the correlation measure. This is because terms such asĉ 1ĉN and its conjugate do not contribute to the expression iγ 1γ2N in the thermodynamic limit.
The presented evidence hints that end-to-end correlations are good indicators of the effects generated by edge modes in one dimension. Due to the topological features of these systems, it is reasonable to assume that this result is robust in the presence of disorder or interaction. As a consequence, correlations constitute a useful inspection mechanism whenever a decomposition in terms of diagonal modes is not feasible or a topological analysis of the system's band structure is not practical. Noticeably, the actual correlation measure seems to be relevant. Entanglement, which accounts for quantum correlations, vanishes exponentially as N grows due to mixness developed byρ 1N .

IV. CONCLUSIONS
Arguments supporting the suitability of end-to-end correlations as indicators of unpaired Majorana operators in the Kitaev chain are discussed. The proposal is verified implementing a folding protocol in combination with tensor-state representation to numerically find a given correlation criterion. The results can be written as a consistent analytical expression that evidence the connection with Majorana fermions. These findings support the hypothesis that the same approach can be used in systems with additional elements like disorder or interaction. Given the characteristics of the Kitaev chain, it would be interesting to apply similar methodologies to study Berry phases around the degeneracy points where Majorana fermions are completely localized.
It is quite likely the folding mechanism employed here have potential applications besides the Kitaev chain. First, with some modifications it can be adjusted to calculate time evolution or thermodynamic state. Second, it can also be applied to chains with long-range hopping or long-range pairing. However, in the way the method currently works, it can be used only for integrable models, because it is the diagonal modes what is actually folded. It is therefore desirable to develop a more versatile technique with a broader field of application. Nonetheless, the protocol can still be useful if interaction terms are reduced in a mean-field fashion, although it is not known how reliable such an approach is. Similarly, it is possible the method has applications in the study of open quantum systems and the numerical solution of the Lindblad equation.

V. ACKNOWLEDGMENTS
This research has been funded by Vicerrectoría de Investigaciones, Extensión y Proyección Social from Universidad del Atlántico under the project "Simulación numérica de sistemas cuánticos altamente correlacionados".
The reduction scheme presented in the main text can be used to write the state as a product of tensors [13]. The basic principle behind such a representation is the use of Schmidt vectors [18] that emerge in one-dimensional systems to build basis states that support the global quantum state. Schematically, the state of a fermion chain can be represented as in figure 2, using empty circles for unoccupied sites and black circles for sites with one fermion. Let us initially focus on one site of the chain, arbitrarily chosen. The Hilbert space of that site can be expanded using a local basis |k . The elements of such a basis are |0 , to represent an empty site, and |1 , to represent an occupied site. To complement the Hilbert space of the chain, one can consider the Schmidt vectors covering all the sites to the left of |k , plus the Schmidt vectors to the right, as shown in the upper draw of figure 2. As can be seen, such vectors are represented as |µ ⊢ and |ν ⊣ respectively. As these vectors are taken as a basis, the total quantum state is given as a superposition of such vectors, as follows The variables λ µ and λ ν are Schmidt coefficients and as such are real and positive. Although these coefficients can in principle be absorbed in the definition of the Γ ′ s, their inclusion is an integral part of the protocol. The superposition coefficients are stored in the components of tensor Γ k µν . As a result this tensor is in general complex. Notice that the Schmidt vectors are orthogonal If the same expansion is done for each place of the chain, a set of tensors with no apparent connection among them is created and can serve as a representation of |ψ . One positive aspect of this representation is that a local unitary operation like (31) acting on site l has a simple implementation The change involves redefining the tensors as follows To see how coefficients from different sites relate, let us take vector |µ ⊢ and write it as a product of the local basis, |j , and the Schmidt vectors to the left, |ξ ⊢ , as indicated in the middle sketch of figure 2. In addition, let us represent this expansion in the following manner Replacing this expression in (A1) it results In the last expression the chain has been divided as a Schmidt decomposition with Schmidt coefficients λ µ . This implies that the set of vectors must be a set of Schmidt vectors to the right, making µ ⊣ |µ ′ ⊣ = δ µ ′ µ . The Schmidt decomposition of the state can then be written in the familiar form Let us now consider an unitary operation acting on consecutive sites l and l + 1, as for instance transformation (32). The operation can be represented aŝ The resulting expression is no longer an evident Schmidt decomposition but it can be rearranged as one in the next manner. Let us write the operation in parenthesis as In the last step the pairs of indices (ξ, J) and (K, ν) have been replaced by single indices α and β. Notice that grouping indices is essentially a notation change. It resorts to the possibility of joining Hilbert spaces from adjacent sections of the chain. Matrix M α,β has no restrictions apart from normalization. It is in general complex and is not necessarily square. Such a matrix can be written as a product of (less arbitrary) matrices applying a singular value decomposition (SVD) [18] M α,β = Both T α,α ′ and T β ′ ,β (different matrices) are complex and unitary, their rows (or columns) being orthogonal vectors. They are also square matrices. Matrix Λ α ′ ,β ′ is real and diagonal.
Normalization requires λ 2 1 + λ 2 2 + λ 2 3 + ... = 1. All the λ's are positive. In many numerical applications the number of λ's is artificially fixed. Here the number of coefficients is handled dynamically and only those below numerical precision are discarded.
The double sum in (A10) can be reduced to a single sum One can in addition write the labels α and β in terms of the original labels In the last step the components of the T 's have been renamed. Notice that the Γ's in the last sum are different from the ones appearing in the initial state. No emphasized distinction is made in order not to overload the notation, but tensors with µ ′ are all new. Also notice that neither λ ξ nor λ ν have changed. As J and K are integer labels without explicit meaning, it is valid to rename them with their lower-case equivalents j and k.
Introducing the final expression in (A9) giveŝ The state is in "canonical form", i.e., written with respect to the (new) Schmidt vectors of the chain, since the states formed as and are orthogonal and normalized because they are the entries of matrices T α,α ′ and T β ′ ,β respectively. This representation is very convinient to calculate local mean values. Using (A1) it can be shown that the reduced density matrix of a given site iŝ A mean value corresponding to a matrixτ that operates exclusively on that site can be found as One can work in an analogous way in the space of two consecutive positions using the corresponding reduced density matrix. It can be shown that this matrix can be written aŝ Sometimes it is also useful to know how to obtain the state coefficients in the Fock basis in terms of this tensor representation. Such a relation can be derived following the arguments in reference [13], thus yielding These operations can be efficiently computed if the number of Schmidt coefficients involved is not too large.

Example
Let us initially consider a chain with no fermions. In the Fock basis the state is given by When this state is split in adjacent parts, the corresponding Schmidt decomposition is trivial: There is one vector to the left, one vector to the right and the only Schmidt coefficient is 1. From this observation the canonical decomposition can be built directly λ = 1, Γ 0 11 = 1, Γ 1 11 = 0.
The same pattern repeats for every place of the chain. Now let us consider the following unitary operation For simplicity let us suppose thatÛ acts on the first two places. The action of this operator on the state iŝ To build a canonical decomposition (the canonical decomposition is not unique), one sees the state as a composition of a local basis plus the Schmidt vectors to the right and left. Taking the local basis of the first site, the state can be written aŝ Vectors |ν 1 = |00...0 and |ν 2 = |10...0 are normalized and orthogonal, therefore, they are valid Schmidt vectors. In this form it is possible to read out the canonical coefficients, finding The superscript [1] is added to emphasize that these coefficients correspond to the first site. With respect to this decomposition, the estate can be visualized in the following manner (with the superscript omitted) U [3] |ψ = λ 1 Γ 0 11 |0 |ν 1 + λ 2 Γ 1 12 |1 |ν 2 .
This case has been sufficiently simple to allow a direct determination of the canonical representation. In other circumstances the protocol presented in the previous section can be used to build a representation in accordance with the original proposal using a systematic approach. To find the reduced density matrix the state is written as a tensor product making explicit reference to the components of each site |ψ = α1,α2,...,αN−1 |α 0 α 1 [1] |α 1 α 2 [2] ... where The superscript in square brackets is used to specify the position in the chain associated to the corresponding tensor. Such a superscript is dropped in the subsequent development to simplify the notation. The reduction can be effectuated by bracketing corresponding spaceŝ The above expression can also be written as a concatenation of index contractions Thus, it all can be written as a single connecting matrix The calculation of M {α1α ′ 1}{αN −1α ′ N −1} unavoidably involves all the tensors in the bulk of the representation and is the numerically heaviest task of the procedure. The resulting expression is a 4 × 4 matrix in the Fock basis of the chain edges.