Majorana fermions in the nonuniform Ising-Kitaev chain: exact solution

A quantum computer based on Majorana qubits would contain a large number of zero-energy Majorana states. This system can be modelled as a connected network of the Ising-Kitaev chains alternating the “trivial” and “topological” regions, with the zero-energy Majorana fermions localized at their interfaces. The low-energy sector of the theory describing such a network can be formulated in terms of leading-order couplings between the Majorana zero modes. I consider a minimal model exhibiting effective couplings between four Majorana zero modes – the nonuniform Ising-Kitaev chain, containing two “topological” regions separated by a “trivial” region. Solving the model exactly, I show that for generic values of the model parameters the four zero modes are localized at the four interface points of the chain. In the special case where additional inversion symmetry is present, the Majorana zero modes are “delocalized” between two interface points. In both cases, the low-energy sector of the theory can be formulated in terms of the localized Majorana fermions, but the couplings between some of them are independent of their respective separations: the exact solution does not support the “nearest-neighbor” form of the effective low-energy Hamiltonian.

Physicists have been fascinated with Majorana fermions ever since their discovery 1, 2 in 1937, when Ettore Majorana found a completely real (i.e. not containing complex coefficients) representation of the Dirac equation. The solutions of the Majorana equation describe neutral fermions -particles that obey the Fermi statistics, but at the same time are their own antiparticles. Whether they exist in nature as elementary particles is still an open question. It has been hypothesized that neutrinos might be Majorana fermions. This hypothesis could be experimentally confirmed by observation of an elusive process known as the neutrinoless double beta decay 3 , which is the focus of considerable experimental efforts.
Recently, it has become possible to imitate the ideas of the relativistic field theory in solids. Following the success of graphene research, further novel materials have been identified as Dirac 4 and Weyl 5, 6 semimetals. These materials exhibit a finite number of band crossings at the Fermi level (the so-called Dirac and Weyl points) 7,8 . To a good approximation, low-energy excitations near these points are characterized by the (quasi)-relativistic spectrum allowing one to observe phenomena previously belonging to the realm of high energy physics, such as the Bell-Adler-Jackiw chiral anomaly 9-12 . At the same time, signatures of Majorana fermions were found in nanowires with proximity-induced superconductivity [13][14][15][16][17][18][19][20] . While the physics of such systems is rather complex, the effective low-energy Hamiltonian describing the nanowire is essentially that of the one-dimensional (1D) p-wave superconductor, i.e. the continuous limit of the Kitaev model [21][22][23][24] . The 1D lattice model proposed by Kitaev 21 exhibits a quantum phase transition between two gapful (massive) phases, known as the "trivial" and "topological". The "trivial" phase is characterized by a single non-degenerate ground state, while the "topological" phase possesses a ground state that is nearly doubly degenerate: for any finite-size, open chain the difference between the energies of the lowest-lying excited state and the ground state is exponentially small, ∝ −  L exp( / ) 0 , in the length of the chain L (here  0 is a certain correlation length defined below). In the thermodynamic limit, the energy difference vanishes and the ground state becomes truly degenerate. This is a manifestation of a well known theorem in statistical physics 25 : spontaneous symmetry breaking and the corresponding vacuum degeneracy may only occur in the thermodynamic limit. The lowest-lying excitation in the "topological" phase of the Kitaev model is a fermion with a wave function that Scientific RepoRts | 7: 1447 | DOI: 10.1038/s41598-017-01413-z is nonzero (with exponential accuracy) only near the two edges of the chain 21 . This fermion can be described in terms of two Majorana fermions, one at each edge 21,26 . It is these objects, known as the "Majorana bound states" 13,20,24 or "Majorana zero modes" 14,17,19 , that have been observed.
Arguably the main driving force behind the pursuit of Majorana zero modes in solids is the possibility of applications to quantum computing 21,27,28 . The basic building block of a quantum computer, the qubit, can be realized as a coupled system of four spatially separated zero-energy Majorana states 26,29,30 . It is expected that a Majorana qubit would have a rather long coherence time due to its topological nature 29,31 . Quantum computer can then be envisioned as a connected network of such qubits. Certain logical operations in such a computer can also be performed topologically by means of braiding (or adiabatic interchange) of Majorana fermions 29,32 .
Alternatively, one can search for Majorana fermions in manifestly discreet systems [33][34][35][36] . For instance, one may engineer the Majorana bound states using Josephson qubits 37,38 to build an artificial spin chain 33,34 that is designed to be an experimental realization of the 1D quantum Ising model [39][40][41][42][43] . The quantum Ising chain with open boundary conditions is formally equivalent to the Kitaev chain 21,28,34,42,43 (note, that the two models do not enjoy the same level of the topological protection 21,44 ). This equivalency is based on the Jordan-Wigner transformation 45 that is commonly used in 1D theories to express the spin-1/2 operators in terms of creation and annihilation operators of spinless fermions 40 . In fact, the original solution 43 of the 1D quantum Ising model was based on the consequent application of the Jordan-Wigner transformation and the Bogolyubov transformation 46 , mapping the model onto a system of free fermions 42 . The simplicity of the resulting physical picture may be deceptive, since both the Jordan-Wigner and Bogolyubov transformations are nonlocal 44 . Although the original Hamiltonian contains only nearest-neighbor couplings, the model may develop long-range correlations. In fact, the ground state of the open-ended chain is characterized by the "end-to-end" correlation function 43 that vanishes in the "trivial" phase (in the thermodynamic limit), but remains finite in the "topological" phase. This result can be interpreted in terms of a nonlocal fermion operator that is a linear combination of the Jordan-Wigner fermions at both ends of the chain. The lowest excited state of the open-ended chain in the "topological" phase (i.e. the state that is nearly degenerate with the ground state) possesses a similar structure. The wave function of this state decays exponentially away from the chain ends and hence can be represented as a linear combination of the two states localized at either end of the chain. The existence of such edge states has been known for a very long time 42 , but they were not interpreted in terms of Majorana fermions and related to the quantum information theory before the work of Kitaev 21 .
A quantum computer based on Majorana qubits would contain a large number of Majorana zero modes. Whether the device will be built using the nanowires [22][23][24] or the artificial spin chains 33,34 , one can envision the effective model of the system as a connected network of the Ising-Kitaev chains alternating the "trivial" and "topological" regions, with the zero-energy Majorana fermions localized at their interfaces 27 . The low-energy sector of such a theory can be formulated in terms of leading-order couplings between the Majorana zero modes 27,30,47,48 . These couplings are often chosen based on physical intuition. Given the nonlocal relation between the Majorana zero modes and the Kitaev (or Jordan-Wigner) fermions, it is desirable to test that intuition against a rigorous solution of a representative microscopic model. This is the principle goal of the present work.
In this paper I consider a minimal model exhibiting effective couplings between Majorana zero modes -the nonuniform Ising-Kitaev chain, containing two "topological" regions separated by a "trivial" region. Based on the common intuition, one would expect that this model possesses four Majorana zero modes, each localized at one of the four interface points of the chain 27, 47 (i.e. the two chain ends and two edges of the "trivial" region). I present the exact solution of the model and identify the region of model parameters where the above expectation is indeed fulfilled. However, the exact solution also exhibits situations where the intuitive expectation is not fulfilled. In particular, inversion symmetry (in the case where the two "topological" regions are identical) leads to "delocalization" of the Majorana zero modes between two interface points. While one can use a basis rotation to express the low-energy sector of the theory in terms of four localized Majorana operators, the corresponding states will no longer be the eigenstates of the model. The low-energy Hamiltonian will then contain effective couplings between some of these modes that are independent of their respective separations. I also demonstrate that the symmetric case in not the only situation exhibiting the "delocalization" of the Majorana bound states. As an example, I show that the "delocalization" may also occur in the variant of the model, where one of the chain ends is coupled to one of the intermediate sites forming a T-junction (or a Y-junction), see Fig. 1.

Results
The nonuniform Ising-Kitaev chain. The open-ended, nonuniform quantum Ising chain is described by In this paper I focus on the minimal model supporting effective couplings between Majorana zero modes choosing the applied field h n to be piece-wise uniform (see Fig. 2  In this case, the chain is split into three regions such that the two "topological" regions (of the length N 1 and N 2 = L − N 1 − M) are separated by the "trivial" region of the length M. Since physical properties of the model are determined by the ratios of the applied fields to the exchange coupling J, it is convenient to factor out the exchange constant J introducing the parameters The first term in Eq. (5) describes the excitation spectrum of the model in terms of free fermion operators η k . The two remaining terms yield the ground state energy. For an arbitrary choice of λ n , the energies  k can be found The chain contains two "topological" regions (the blue line, with N 1 sites, and the purple line, with N 2 sites) and one "trivial" region (the red line, with M sites). The T-junction is located at the site N 0 . The peaks represent the absolute values of the real-space amplitudes of the Majorana zero modes calculated for N 1 = M = 20, N 2 = 10, and N 0 = 11. The model exhibits four Majorana zero modes (corresponding to the nearly four-fold degeneracy of the ground state). Two of them are localized at a single interface point each: the dark blue at the site N 1 and the green at the T-junction. Note, that this amplitude is spread over only two (out of three) branches at the junction. The remaining two zero modes are delocalized between two interface points, the sites 1 and N 1 + M. One of them is illustrated by the red peaks.

Figure 2.
The nonuniform Ising-Kitaev chain split into two "topological" (dark green) and one "trivial" (red) regions. The first "topological" region is characterized by the parameter λ 1 < 1 and occupies the left part of the chain, ⩽ ⩽ n N 1 1 . The next M sites are occupied by the "trivial" phase with λ 2 > 1. The remaining portion of the chain of the length = − − N L N M 2 1 is occupied by the second "topological" region with λ 3 < 1.
numerically with arbitrary precision. In the specific case (3), the model can also be solved analytically. Below I present the results of the analytic solution and compare them to the numerical results.
Nearly degenerate ground states. The Ising-Kitaev chain split into two "topological" and one "trivial" region possesses two single-particle excitations (hereafter denoted by k = 1, 2) that are nearly degenerate with the ground state. As long as the parameters λ i are not too close to unity and the sizes of the three regions are not too small, such that the three quantities λ N 1 2 1 , λ N 3 2 2 , and λ − M 2 2 are exponentially small, the energies and the wavefunctions of these states can be found analytically. Already the leading-order expression shows excellent agreement with the exact numerical diagonalization of the model as illustrated in Figs 3 and 4. The visible discrepancy between the analytic and numerical results for λ 1 ~ 1 is to be expected: there the above parameters cease being exponentially small and the approximate analytic expressions become invalid.
Without specifying the relation between the three exponentially small parameters, even the leading-order expressions for the two eigenvalues 1,2  are rather cumbersome. Therefore, here I focus on two representative limiting cases (the symbolic expression λ ( ) N  hereafter denotes the omitted subleading, exponentially small terms).  Asymmetric chain. If the two "topological" regions of the chain are not symmetric, then compact expressions for the energies 1,2  can be found under following assumptions.
, the two nearly zero-energy states are determined by the two "topological" regions of the chain, independently of the size of the "trivial" region. The first excited state has the energy while the energy of the second excited state is These results are illustrated in Fig. 3 by the dashed lines to the right of the vertical grid line (marking the end of the above parameter region λ λ ). Vanishing of the energies (6) and (7) at the point λ 1 = λ 3 = 1 is the artifact of the approximation. As the parameters λ i approach unity, the approximate expressions reported here become invalid (while it is possible to write down exact expressions for 1,2  that are valid also near the critical point, their algebraic complexity renders them practically useless).
, the larger eigenvalue  2 is determined by the "trivial" region of the chain while the energy of the lowest excited state is determined by the two "topological" regions combined These results are illustrated in Fig. 3 by the dashed lines on the left side of the vertical grid line.
Symmetric chain. In the symmetric case, λ λ = N N 1 3 1 2 , the two energies (7) and (6) coincide. In this case, one has to consider the subleading terms neglected so far, as the eigenvalues of the finite-size chain (2) are never truly degenerate.
, the resulting energies are given by This result is illustrated in Fig. 4 to the right of the vertical grid line (the exponentially small difference between the two energies (10) is indistinguishable on the scale of the plot).
• "Strong barrier". In the limit λ , the leading behavior of the energy eigenvalues is given by Eqs (6) and (7). The corresponding amplitudes α 1(2)n and β 1(2)n can also be written in compact from, again retaining only the leading exponential terms. The first excited state (6) is characterized by the amplitudes where the symbolic expression  λ ( ) N denoting the subleading terms is omitted in some lines for brevity and Hence with exponential accuracy, the lowest-energy excitation of the model can be described by the single fermion, η 1 , confined to the second "topological" region of the chain, cf. Eq. (12). The amplitudes (15) are illustrated in the two top panels in Fig. 5. The spread of the localized Majorana states over several lattice sites exhibited by Eq. (15) is a generic feature 21 that can be seen also in the continuum limit 32 . The second eigenvalue (7) is characterized by the amplitudes The corresponding excitation η 2 is confined to the first "topological" region of the chain. The amplitudes (17) are illustrated in the two bottom panels in Fig. 5. The results (15) and (17)  the two lowest-energy excitations of the model behave similarly to those of the two independent "topological" regions. In particular, they exhibit four nearly zero-energy Majorana fermions localized at the edges of the "topological" regions.
• "Weak barrier". In the limit, λ λ λ −  ,  The amplitude β 1n is identical with Eq. (15), but the amplitude α 1n has changed. In the first "topological" region of the chain, it behaves as the corresponding amplitude of the second excited state (17) of the strong barrier case. Moreover, there is a nonzero probability to find this quasiparticle also at the interface between the "trivial" and the second "topological" regions, see Fig. 5, i.e. the corresponding Majorana fermion is essentially "delocalized" between the two points.
0 999 1,1 and α = . 0 107 1,31 . Whether this feature survives the thermodynamic limit depends on what happen to the value of the product of a small and large parameters λ λ N M 1 2 1 in the limiting procedure. The second eigenvalue (7) is characterized by the amplitudes , 1 , , 1 , ( ), 1 . In contrast to the strong barrier case, the wave-function of lowest-energy fermion η 1 is now mostly spread between the two outer edges of the chain, with a small weight at the interface between the "trivial" and the second "topological" region. The second excitation η 2 is mostly confined to the edges of the "trivial" region, with the small weight at the beginning of the chain, see Fig. 5.

Symmetric chain.
• The two lowest-energy excitations of the symmetric chain with the strong barrier, , are characterized by the amplitudes, see Fig. 6 (the symbolic expression  λ ( ) N denoting the subleading terms is omitted for brevity) In this case, the excitations η † 1,2 are no longer confined to one of the two "topological" regions of the chain, but are spread symmetrically over both of them. • In the weak barrier case, the amplitudes α n and β n are still described by Eqs (18) and (19).
The exponential decay of the amplitudes (15)- (20) can be described in terms of a correlation length, which is specific to each of the three regions of the chain In experiments on discreet systems [33][34][35][36] , the realistic values of λ i might not be extreme and hence the correlation lengthes (21) might not be very small. In such case, even the localized Majorana fermions are spread over several lattice sites as illustrated in Figs 5 and 6. Effective low energy theory. Applications to quantum computation 27 involve adiabatic manipulations of the Majorana zero modes. This means that any external perturbation applied to the system should be slow enough to avoid exciting higher-energy gapped states. The remaining low-energy sector of the theory consists of the ground state |GS〉 and the nearly degenerate excitations that can be interpreted in terms of Majorana zero modes.
For the specific model considered in this paper, the low-energy sector contains four states where the last state is the two-particle excitation. These four states can be further split into two groups of mutually orthogonal states, belonging to the two parity sectors of the model where the total fermion number is either even or odd. Projecting the Hamiltonian (2) onto either of the above sectors, one finds the effective low-energy theory. In the one-fermion (odd) sector, the effective Hamiltonian has the simplest form in the basis of the Majorana fermions γ α 1(2) ( ) and γ β 1(2) ( ) .
• In the asymmetric chain with the strong barrier, i.e. in the limit λ , the localized Majorana fermions describe the exact eigenstates of the model. Hence, the projected Hamiltonian in the basis the four Majorana states (counted from left to right) has the block-diagonal structure , the Majorana amplitudes (20) are not localized at single interface points in the chain. One can still represent the effective Hamiltonian in the basis of the four localized Majorana fermions. However, now these objects are no longer associated with the exact eigenstates and hence additional couplings appear. Introducing a short-hand notation for the eigenvalues (10) 1 (2) sym I find the following Hamiltonian (using an obvious basis rotation) The Hamiltonian (25) is no longer block-diagonal: the two pairs of the Majorana zero modes, γ 2 , γ 3 and γ 1 , γ 4 , are now coupled. The effective coupling of both pairs is the same (given by δ) despite the large difference in their respective separation. This result is the consequence of the nonlocal nature of the eigenstates of the model. • In the weak barrier case, the low-energy excitations are no longer confined to the "topological" regions of the chain. In order to take into account the "delocalization" of the amplitudes α 1(2)n , see Eqs (18) and (19), I denote the ratio of the peak values of α 1(2)n as κ λ λ ≈  s c 1/ 3 1 N M 1 2 1 (up to exponentially small corrections) and find the low-energy Hamiltonian (in the same basis of γ i counted from left to right) In contrast to Eqs (23) and (25), the leading terms in Eq. (26) are confined to the antidiagonal. This results reflects the structure of the single-fermion excitations (18) and (19): the lowest-energy excitation (18) is mostly spread between the two outer edges of the chain (and is described by γ 1 , γ 4 ), while the fermion (19) -to the edges of the "trivial" region (corresponding to γ 2 , γ 3 ). Neglecting weak "delocalization" completely (i.e. in the limit κ → 0), the Hamiltonian (26) can be made block-diagonal by renumbering the localized Majorana operators.

Kitaev chain with a T-junction.
Consider now a modified model where one of the chain ends is coupled to an intermediate site forming a T-junction (sometimes also referred to as a Y-junction), see Fig. 7. Such a model can be easily formulated in the fermion language by adding another coupling term to the Hamiltonian (2) The Hamiltonian (27) is quadratic and can be diagonalized by the same exact method as the model (2). Due to the more complicated topology and larger number of parameters, the model (27) exhibits many more parameter regimes than the chain (2). A comprehensive discussion of all of these regimes will be presented elsewhere. In this paper, I focus on a single parameter regime, where the nearly zero-energy eigenstates are described by the amplitudes α n that are delocalized between two interface points similarly to Eq. (20). In particular, I consider the regime, where the junction site N 0 is approximately in the middle of its "topological" region, which is twice as long as the second "topological" region. Choosing the equal parameters λ 1 = λ 3 < 1 describing the "topological" regions, I achieve a configuration that is somewhat analogous to the symmetric chain. is occupied by the second "topological" region with λ 3 < 1. The T-junction is located at the site N 0 . In this paper, I report results for the case where the T-junction is approximately in the middle of the first "topological" region of the chain. The results of the exact numerical diagonalization of the model (27) in the chosen parameter regime (with N 1 = 20, M = 20, N 2 = 10, and N 0 = 11) are presented in Fig. 8 by the red dots and illustrated in Fig. 1. The solution is characterized by the amplitudes α 1,2 sym that are delocalized between the edge point of the chain and one of the borders of the "trivial" region. The amplitude α 2 sym appears to be in perfect agreement with Eq. (20), while the amplitude α 1 sym shows a barely perceptible deviation (see the solid curves in the two left panels in Fig. 8). Now the amplitudes β 1,2 sym are no longer delocalized. The amplitude β 1 sym is localized at the edge of the "trivial" region and is perfectly described by Eq. (17), see the top right panel in Fig. 8. The amplitude β 2 sym is localized at the T-junction point, see the bottom right panel in Fig. 8, and can be described analytically by The asymmetry of this amplitude -β 2 sym is spread over only two out of the three branches of the junction -is related to the chirality of the junction 50 described by the explicitly asymmetric Hamiltonian (27), as well as to the time reversal symmetry 52 of the Hamiltonian (27).

Discussion
In this paper I presented the exact analytic solution of the nonuniform Ising-Kitaev chain with open boundary conditions. The motivation for this work was two-fold: (i) I was motivated by the proposal 33, 34 of experimental realization of zero-energy Majorana states in an artificial spin chain engineered using Josephson qubits; such a system would be discreet and, given current technological limitations, contain not too many qubits; (ii) I wanted to reach a better understanding of the the effective coupling between the Majorana zero modes in networked systems used as paradigmatic examples of possible applications to quantum computing 27 , in particular, the common "nearest-neighbor" form of the effective low-energy Hamiltonian 30,47,48 . The model solved in this paper represents the first step in reaching these goals. In the case, where the two "topological" regions of the chain are separated by a "trivial" region, the exact analytic form of the eigenvalues and eigenvectors can be found. Moreover, the resulting low-energy theory does not support the "nearest-neighbor" approach 30,47,48 .
Indeed, in the generic parameter regime with the "strong barrier", the single-fermion sector of the effective low-energy theory contains two single-fermion states, each confined to its own "topological" region as if these regions were disconnected. Similarly to the original Kitaev model 21 , each of these states can be interpreted in terms of two Majorana states localized at the edges of the "topological" regions, as expected 27,30,47,48 . In the basis of the localized Majorana fermions, the effective Hamiltonian has the block-diagonal form (23). Here the two pairs of Majorana fermions γ 1 , γ 2 and γ 3 , γ 4 form the two single-fermion eigenstates. The orthogonality of the single-fermion eigenstates leads to the absence of any coupling between γ 2 and γ 3 , that is typically included in the "nearest-neighbor" approach 30,47,48 .
In contrast, in the specially fine-tuned case of the symmetric chain the single-fermion eigenstates are equally spread between the two "topological" regions. Consequently, the corresponding Majorana fermions are localized not at one, but at two separate interface points. While the resulting low-energy Hamiltonian can of course be represented in the above basis of localized Majorana states, the latter are no longer related to the exact eigenstates of the model. As a result, the Hamiltonian (24) exhibits additional couplings between the two pairs γ 2 , γ 3 and γ 1 , γ 4 . These couplings are identical despite the large difference in the separation between γ 2 , γ 3 and γ 1 , γ 4 . Again, this contradicts the "nearest-neighbor" approach 30,47,48 , where the coupling between γ 1 , γ 4 is not included due to the larger separation (as compared to other pairs of Majorana fermions). Now, in the case of the weak barrier the amplitudes (18) and (19) also exhibit "delocalization" between two interface points (although to a significantly lesser degree). Here, the dominant terms in the effective Hamiltonian (26) are the couplings between γ 2 , γ 3 and γ 1 , γ 4 , while the couplings γ 1 , γ 2 and γ 3 , γ 4 appear with the typically small coefficient κ (due to weak "delocalization"). As a result, the Hamiltonian (26) is also incompatible with the "nearest-neighbor" form 30,47,48 .
Finally, "delocalization" of the Majorana zero modes between two separated interface points is not an artifact 53 of the model (2). In particular, the modified model (27) also exhibits this "delocalization" (again, requiring some fine-tuning).
The results of this paper are relevant for experimentalists designing small systems hosting multiple zero-energy Majorana states [33][34][35][36] . In particular, in systems involving relatively few Josephson qubits with conservative parameter values the spreading of the Majorana zero modes over a few qubits and their "delocalization" between two well separated points are generic effects 53 that need to be taken into account while interpreting the experimental data and especially when planning any kind of manipulation of the Majorana bound states by some external bias. "Delocalization" of the Majorana zero modes and the corresponding separation-independent effective coupling between some pairs of the localized Majorana fermions could be observed in experimental nanowire samples in the presence of additional symmetries incorporated into the system design.

Methods
An arbitrary quadratic Hamiltonian, including Eqs (2) and (27)   In the particular case of the Hamiltonian (2), the resulting diagonal form is given by Eq. (5). The outlined diagonalization procedure is applicable to any quadratic Hamiltonian in any dimensionality and allows for an efficient numerical calculation of the eigenvalues k  and the coefficients α kj and β kj with arbitrary precision. However, analytic solution is manageable only in a few relatively simple cases. Fortunately, the model (2) with the specific choice (3) of the applied fields is one of them. The exact single-particle energies k  of this model can be expressed as λ λ λ λ λ λ = + + ϑ = + + ϑ = + + ϑ 1 2 cos 1 2 cos 1 2 cos ,  The latter equalities in Eq. (34) provide additional constraints on ϑ i , which guarantee the uniqueness of the solution.
Consider now the energy (34) as a function of ϑ i , regardless of which values of ϑ i are allowed by Eq. (35). For real ϑ i , this function exhibits a minimum at π ϑ = i . The minimum value of the energy gives a reasonable lower bound for the bulk gap Hence, any subgap states including nearly zero-energy states are described by complex solutions to Eq. (35).