Majorana zero modes in Y-shape interacting Kitaev wires

Abstract

Motivated by the recent experimental realization of minimal Kitaev chains using quantum dots, we investigate the Majorana zero modes (MZM) in Y-shape Kitaev wires. We solve the associated Kitaev models analytically at the sweet spot (t_h = Δ) and derive the exact form of MZM wave-functions in this geometry. We observed exotic multi-site MZMs located near the junction center, on the nearby edge sites of each leg. This result is important for potential braiding of Majoranas and the performance of Y-junctions made from arrays of quantum dots. Furthermore, we study the stability of local (single-site) and multi-site MZMs modes in the presence of Coulomb repulsion, using density matrix renormalization group theory. Our local density-of-states calculation shows that these multi-site MZMs are as equally topologically protected as the single-site MZMs when in the presence of Coulomb repulsion or when away from the sweet-spot.

Introduction

Majorana zero modes (MZMs) are charge-neutral non-Abelian quasiparticles^1,2. They have attracted much interest because of their potential application in fault-tolerant topological quantum computing^2,3,4,5. The occurrence of zero-bias peaks in tunneling spectroscopy is one of the experimental signatures of MZMs⁶. A promising platform to realize MZMs are semiconductor nanowires proximitized to superconductors⁷, and magnet-superconductor hybrid systems^8,9,10,11, where MZMs are expected to develop at both ends of the wire. However, in these materials the so-called “sweet spot” in couplings, where Kitaev theory is primarily developed, is difficult to reach. Recently, the realization of MZMs was also proposed in quantum-dot-superconductor linear arrays^12,13,14. These quantum-dots systems^{15,16,17,18,19,20} are expected to reduce the problem of random-disorder potential, as compared to the proximitized semiconductor nanowires where the effect of disorder is strong¹² and may create false signals of MZMs in tunneling spectra²¹. Interestingly, the experimental realization of a minimal Kitaev chain has been demonstrated using two quantum dots coupled through a short superconducting-semiconductor hybrid (InSb nanowire)²². Remarkably, in this experiment²², two localized MZMs were observed in tunneling conductance measurements at the sweet spot t_h = Δ. Very recently the experimental realization of a three-site Kitaev chain geometry in a quantum-dot system has been reported²³, opening the doors towards the realization of longer Kitaev chains using many coupled quantum dots.

In topological quantum computation, it is required to move and perform braiding operations of the MZMs^24,25. A strict 1D geometry is not sufficient to perform such braiding operations, because in 1D the MZMs can fuse during their exchange process^26,27. To realize non-Abelian statistics (or braiding), T and Y-shaped wires geometry have been proposed^{24,28,29,30,31,32}. It has been shown that the MZMs in T-shape nanowires can be transformed under exchange, similarly to 2D p + ip superconducting systems, displaying non-Abelian statistics²⁴. Braiding-based gates with MZMs using quantum dot arrays were recently proposed³².

This paper focuses on finding MZM modes near the junction of interacting Y-shaped quantum dots near the sweet-spot (t_h=Δ), important for braiding and for demonstrating the non-Abelian statistics of MZMs. However, in the context of proximity-induced semiconductor nanowires (in the limit of Δ ≪ t_h), there are only a few studies related to ground states of T-shaped or Y-shaped wires^33,34,35. The sub-gap properties of a three-terminal Josephson junction (composed of effective spinless p-wave superconductors), joined into a T-shaped normal-metallic region, has been studied using the scattering matrix approach in the non-interacting limit³⁴. They found that depending upon the superconducting phase of each arm, the Majorana zero mode extended into the metallic region of either all three legs or two legs of the T-shape wire. In Ref. ³⁶ the existence of a zero-energy parafermionic mode was proposed in multi-legged star junctions of quantum Hall states. However, in these studies the precise form of the Majorana wave functions and the effect of Coulomb interactions were not addressed^33,34. The repulsive Coulomb interaction is expected to suppress the pairing-induced bulk-gap and can affect the stability of Majorana modes^13,37.

Motivated by the above described recent progress in the realization of minimal Kitaev chains using quantum dots^22,23, we study the Majorana zero modes in Y-shaped interacting Kitaev wires. We address the exact form of Majorana wavefunctions near the junction and its dependency on the superconducting phases at each wire. Assuming each arm of the Y-shape wire can take different values for the superconducting (SC) phase (see Fig. 1a), we solve exactly the Kitaev model for Y-shapes wires working at the sweet spot t_h = Δ. Remarkably, in terms of Majorana operators, we are able to write four independent commuting Hamiltonians, consisting of the three arms (I, II, III) and one central region (IV), as shown in Fig. 1b. This allow us to diagonalize the full system independently for each region and solve exactly at the sweet spot. We find the expected three single-site localized MZMs on the edge sites of the Y-shape wire. Surprisingly, depending on the SC phase values, we also find exotic multi-site MZMs near the central region, not bounded to be on only one site as in the Kitaev model³⁸. Instead they are located on two or three edge sites of the different legs of the Y-shape wire. Note that this is conceptually totally different from the standard exponential decay of the Majorana wave function away from the sweet spot or under non-ideal conditions. The multi-site nature of the MZMs unveiled here is an intrinsic, previously unknown, property of the Y-geometry system at the sweet spot.

Fig. 1: Schematic representation of the Y-shape Kitaev wires at the sweet-spot.

a ϕ₁, ϕ₂, and ϕ₃ are the phases of the p-wave superconductor in each arm of the Y-shape wire. Arrows denote the directions of the pairing terms Δ and site index j in each wire. b Pictorial representation of the Majorana zero modes in the Y-shape Kitaev wire unveiled here, at t_h = Δ. There are always three localized MZMs γ (purple color) at each end of the Y-geometry. Depending upon ϕ₁, ϕ₂ and ϕ₃ for each arm, near the wire’s junction we observed only two situations: either (i) one multi-site MZM (χ) or (ii) two multi-site MZMs (χ) and one single-site MZM (γ). Green color MZMs form ordinary fermions at each bond.

Furthermore, we perform much needed unbiased simulations of the Majorana zero modes in Y-shaped interacting Kitaev wires, using density matrix renormalization group (DMRG)^39,40. In the non-interacting limit, we find peaks in the site-dependent local density-of-states (LDOS) at the exact locations predicted by our analytical calculations. In order to compare the stability of single-site vs. multi-site MZMs, we examine the electron and hole components of the LDOS separately^41,42, against the increase in repulsive Coulomb interactions. Interestingly, the LDOS(ω, j) calculations indicate the multi-site MZMs are as equally stable as the single-site MZMs residing at the ends of the Y-shaped wire.

The Hamiltonian for the Y-shaped Kitaev model, with superconducting phases ϕ₁, ϕ₂, and ϕ₃ at each arm, can be divided into four different parts. The Hamiltonian for each leg can be written as:

$${H}^{I}=\mathop{\sum }\limits_{j=1}^{l-1}\left(-{t}_{h}\mathop{c}\nolimits_{j}^{\,\,{\dagger} }{c}_{j+1}+{e}^{i{\phi }_{1}}{{\Delta }}{c}_{j}{c}_{j+1}+H.c.\right),$$

(1)

$${H}^{II}=\mathop{\sum }\limits_{j=l+2}^{2l}\left(-{t}_{h}\mathop{c}\nolimits_{j}^{\,\,{\dagger} }{c}_{j+1}+{e}^{i{\phi }_{2}}{{\Delta }}{c}_{j}{c}_{j+1}+H.c.\right),$$

(2)

$${H}^{III}=\mathop{\sum }\limits_{j=2l+2}^{3l}\left(-{t}_{h}\mathop{c}\nolimits_{j}^{\,\,{\dagger} }{c}_{j+1}+{e}^{i{\phi }_{3}}{{\Delta }}{c}_{j}{c}_{j+1}+H.c.\right).$$

(3)

Moreover, the Hamiltonian for the central site l + 1 joining each leg edge site can be written as:

$$\begin{array}{l}{H}^{IV}=\left(-{t}_{h}\mathop{c}\nolimits_{l}^{\,\,{\dagger} }{{c}_{l+1}}+{e}^{i{\phi }_{1}}{{\Delta }}{c}_{l}{c}_{l+1}+H.c.\right)\\ \left(-{t}_{h}\mathop{c}\nolimits_{l+1}^{\,\,{\dagger} }{{c}_{l+2}+{e}^{i{\phi }_{2}}}{{\Delta }}{c}_{l+1}{c}_{l+2}+H.c.\right)\\ \left(-{t}_{h}\mathop{c}\nolimits_{l+1}^{\,\,{\dagger} }{{c}_{2l+2}}+{e}^{i{\phi }_{3}}{{\Delta }}{c}_{l+1}{c}_{2l+2}+H.c.\right).\end{array}$$

(4)

Results

Description of analytical method

We noticed that the Y-shape Kitaev model can be solved analytically at the sweet spot t_h = Δ and V = 0, and we will focus on three different sets of superconducting phases: (i) ϕ₁ = π, ϕ₂ = 0, and ϕ₃ = 0, (ii) ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0, and (iii) ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2. First, we divide the system Hamiltonian into four different parts in terms of spinless fermionic operators. Then, we rewrite the system Hamiltonian in terms of Majorana operators, using the transformation ${c}_{j}=\frac{1}{\sqrt{2}}{e}^{-i\frac{{\phi }_{k}}{2}}\left({\gamma }_{j}^{A}+i{\gamma }_{j}^{B}\right)$¹, where k refers to each of the legs. In terms of the Majorana operators, remarkably the system can be written as four independent commuting Hamiltonians: (1) the three independent 1D wires (I, II, III) (see Fig. 1b) and (2) the central region (IV), consisting only of five Majorana operators (two from the central site and three from the near edge sites of each leg), as shown in Fig. 1b. As already expressed, this procedure allows us to solve the Y-shape Kitaev model exactly at the sweet spot t_h = Δ and V = 0, for any values of the SC phases of each arm (for the exact solution of 1D Kitaev chain with π-junction, see Supplementary Note II). For all three cases of SC phases discussed above, we find there is one Majorana zero mode at the outer edge (single-site) of each arm (purple color in Fig. 1b), as expected intuitively. Our main finding is that depending on the phase values of each arm, in addition to the above mentioned outer edge MZMs, we find at the central region either (i) only one multi-site MZM (χ) or (ii) two multi-site MZMs (χ) accompanied with one single-site MZM (γ). These multi-site MZMs (χ) are located at the edge sites (near the central region) of each arm in the crucial limit of t_h=Δ that allows for the analytic solution. These exotic multi-site MZMs could be realized in quantum dot experiments, by changing the phases of each arm in a Y-shape geometry of quantum dots arrays.

In terms of Majorana operators, the Hamiltonian of each leg is independent of the SC phases (ϕ₁, ϕ₂, ϕ₃). Using the transformations ${c}_{j}^{I}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{1}/2}\left({\gamma }_{A,j}^{I}+i{\gamma }_{B,j}^{I}\right)$, ${c}_{j}^{II}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{2}/2}\left({\gamma }_{A,j}^{II}+i{\gamma }_{B,j}^{II}\right)$ and ${c}_{j}^{III}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{3}/2}\left({\gamma }_{A,j}^{III}+i{\gamma }_{B,j}^{III}\right)$, the Hamiltonian for the three legs H^I, H^II and H^III can be written in terms of Majorana operators as:

$${H}^{I}=-2i{{\Delta }}\mathop{\sum }\limits_{j=1}^{l-1}\left({\gamma }_{A,j+1}^{I}{\gamma }_{B,j}^{I}\right),$$

(5)

$${H}^{II}=-2i{{\Delta }}\mathop{\sum }\limits_{j=l+2}^{2l}\left({\gamma }_{A,j+1}^{II}{\gamma }_{B,j}^{II}\right),$$

(6)

$${H}^{III}=-2i{{\Delta }}\mathop{\sum }\limits_{j=2l+2}^{3l}\left({\gamma }_{A,j+1}^{III}{\gamma }_{B,j}^{III}\right),$$

(7)

where we have used the sweet-spot property t_h = Δ. In these equations, the Majorana operators ${\gamma }_{A,1}^{I}$, ${\gamma }_{B,2l+1}^{II}$, and ${\gamma }_{B,3l+1}^{III}$ are absent²⁴, and commute with these Hamiltonians, which indicates the presence of three single-site MZMs at the edge sites of the Y-shape Kitaev wire (see Fig. 1b), similarly as in the original Kitaev chain exact solution. For the central region, the Hamiltonian H^IV depends upon the SC phases. To illustrate the physics unveiled here, we solve H^IV for three different cases, in order to understand the nature of the central MZMs.

The case ϕ ₁ = π, ϕ ₂ = 0, and ϕ ₃ = 0

To explain the exact solution of Majorana wave functions in Y-shape geometries, we start with phase values ϕ₁ = π, ϕ₂ = 0, and ϕ₃ = 0 on arms I, II, and III, respectively. For these phase values, the pairing term at each arm preserves the rotational symmetry of the system, as under 120° rotation around the central sites l + 1, the system Hamiltonian remains invariant due to the π phase in arm I (see Supplementary Note I, for more detail)^43,44.

For the central sites, using the relations ${c}_{l}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{1}/2}\left({\gamma }_{A,l}^{I}+i{\gamma }_{B,l}^{I}\right)$, ${c}_{l+1}=\frac{1}{\sqrt{2}}\left({\gamma }_{A,l+1}^{IV}+i{\gamma }_{B,l+1}^{IV}\right)$, ${c}_{l+2}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{2}/2}\left({\gamma }_{A,l+2}^{II}+i{\gamma }_{B,l+2}^{II}\right)$, and ${c}_{2l+2}=\frac{1}{\sqrt{2}}{e}^{-i{\phi }_{3}/2}\left({\gamma }_{A,2l+2}^{III}+i{\gamma }_{B,2l+2}^{III}\right)$, the sector H^IV (with ϕ₁ = π, ϕ₂ = 0, and ϕ₃ = 0) can be transformed in terms of Majorana operators as:

$${H}^{IV}=-2i{{\Delta }}\left({\gamma }_{B,l}^{I}+{\gamma }_{A,l+2}^{II}+{\gamma }_{A,2l+2}^{III}\right){\gamma }_{B,l+1}^{IV}.$$

(8)

Note that in Eq. (8) the Majorana operator ${\gamma }_{A,l+1}^{IV}$ at the central site l + 1 is absent, signaling the presence of a single-site MZM at site j = 16 (for a L = 46 sites system). Next, we write Eq. (8) in terms of a 4 × 4 matrix in the basis of $\left({\gamma }_{B,l}^{I},{\gamma }_{A,l+2}^{II},{\gamma }_{A,2l+2}^{III},{\gamma }_{B,l+1}^{IV}\right)$ and obtained four eigenvalues ($-\sqrt{3},\sqrt{3},0,0$). The last two eigenvalues (e₃, e₄) are zero denoting MZMs, and the associated eigenvectors ${\chi }_{3}=-\frac{1}{\sqrt{2}}{\gamma }_{B,l}^{I}+\frac{1}{\sqrt{2}}{\gamma }_{A,2l+2}^{III}$ and ${\chi }_{4}=-\frac{1}{\sqrt{6}}{\gamma }_{B,l}^{I}+\sqrt{\frac{2}{3}}{\gamma }_{A,l+2}^{II}-\frac{1}{\sqrt{6}}{\gamma }_{A,2l+2}^{III}$ emerge, with properties ${{\chi }^{{\dagger} }}_{3}={\chi }_{3}^{{\dagger} }$ and ${{\chi }^{{\dagger} }}_{4}={\chi }_{4}^{{\dagger} }$. In addition, χ₃ and χ₄ also commute with H^IV, and its amplitude splits into multiple sites (i.e. not bounded at only one site); these properties confirm the presence of two multi-site MZMs in the central region. The MZM χ₃ is located on two sites j = 15 and j = 32 with equal amplitude, whereas the MZM χ₄ is distributed on three sites j = 15, 17 and 32 with different amplitudes (see Fig. 2a). Note that this is unrelated to exponential decaying wave functions, as it occurs away from the sweet spot. There is no “tail” in these multi-site wave functions but instead the core is fully spread into a relatively small number of sites. The diagonalized Hamiltonian for the central site can be written as

$${H}^{IV}=2\sqrt{3}{{\Delta }}\left({{\chi }_{2}}^{{\dagger} }{\chi }_{2}-\frac{1}{2}\right),$$

(9)

where ${\chi }_{2}=\frac{i}{\sqrt{6}}{\gamma }_{B,l}^{I}+\frac{i}{\sqrt{6}}{\gamma }_{A,l+2}^{II}+\frac{i}{\sqrt{6}}{\gamma }_{A,2l+2}^{III}+\frac{1}{\sqrt{2}}{\gamma }_{B,l+1}^{IV}$ is an ordinary fermion.

Fig. 2: Majorana zero modes in a Y-shape Kitaev wire, for the case ϕ_1=π, ϕ₂₌₀, and ϕ₃₌₀, at the sweet-spot (t_h = Δ = 1).

a There are a total of six MZMs. For V = 0, three single-site MZMs (γ₁, γ₂, and γ₃) at the edge sites j = 1, 31, 46 and in addition one single-site MZM (γ₄) at the central site j = 16. Near this central region, there are two additional multi-site Majoranas: χ₃ (distributed on sites j = 15 and 32) and χ₄ (distributed on sites j = 15, 17, and 32), which leads to a spectral weight 2/3 (compared to the single-site localized MZMs). At robust repulsion V = 2 the Majoranas remain exponentially localized over a few sites. The electron parts of the LDOS(ω) vs. ω are shown for sites (b) j = 16, (d) j = 17, and (f) j = 31 and hole parts of the LDOS(ω) vs. ω are shown for sites (c) j = 16, (e) j = 17, and (g) j = 31. For V≤2, the electron and hole parts of the LDOS(ω) are almost equal and shows peaks at ω = 0, indicating the presence of stable MZMs on these sites.

The other remaining Hamiltonians for each arm can be diagonalized using fermionic operators ${d}_{k,j}=\frac{1}{\sqrt{2}}\left({\gamma }_{B,j}^{k}+i{\gamma }_{A,j+1}^{k}\right)$, where k = I, II, or III. The diagonalized Hamiltonian H^I, H^II, and H^III are:

$${H}^{I}=2{{\Delta }}\mathop{\sum }\limits_{j=1}^{l-1}\left({d}_{I,j}^{\,{\dagger} }{d}_{I,j}-\frac{1}{2}\right),$$

(10)

$${H}^{II}=2{{\Delta }}\mathop{\sum }\limits_{j=l+2}^{2l}\left({d}_{II,j}^{\,{\dagger} }{d}_{II,j}-\frac{1}{2}\right),$$

(11)

$${H}^{III}=2{{\Delta }}\mathop{\sum }\limits_{j=2l+2}^{3l}\left({d}_{III,j}^{\,{\dagger} }{d}_{III,j}-\frac{1}{2}\right).$$

(12)

In conclusion, our analytical calculations predict a total of six MZMs, a surprisingly large number, for the case of ϕ₁ = π, ϕ₂ = 0 and ϕ₃ = 0. These six MZMs lead to eight fold-degeneracy in the ground state of the system, after including the result of Majorana fusion⁴, which we also find is fully consistent with our numerical Lanczos calculations (see Supplementary Note III, for more details).

Next, we analyze the stability of MZMs in the presence of a nearest-neighbor interaction H_I = Vn_jn_j+1, where n_j is the electronic density at site j. We calculate the LDOS, using DMRG for a Y-shaped geometry with system size L = 46. As shown in Fig. 2a, the site dependent LDOS(ω = 0, j) shows sharp peaks for the edge sites j = 1, 31, and 46, indicating three single-site localized MZMs at each edge of the arms. At the central site l + 1, there is a sharp peak with same height as for the edge sites, showing the presence of a single-site MZM γ₄ at site j = 16, as already discussed. Interestingly, there are three other peaks in that LDOS(ω = 0, j) on sites j = 15, 17, and 32. From the expressions of χ₃ and χ₄ presented in the previous page, it can be deduced that the height should be exactly 2/3 for sites j = 15, 17, and 32 as compared to the edge sites. These three peaks signal the presence of two multi-site MZMs, distributed over three central sites (j = 15, 17, and 32).

Increasing the interaction strength to V = 2, the Majoranas are no longer strictly localized at a single site j. Now the MZMs are decaying over a few more sites and consequently the peak height of the LDOS(ω = 0, j) decreases (Fig. 2a). This shows that the Majorana zero modes are topologically protected against moderate values of the Coulomb interaction.

To compare the topological protection against V, for single-site and multi-site MZMs, we calculate the electron and hole parts of LDOS(ω, j), separately for the edge and central sites (right panel of Fig. 2). In Figs. 2b, c, we show the electron and hole part of LDOS(ω, j) for the central site j = 16. With increasing V, the peak height of the electron and hole portions of the LDOS(ω) at ω = 0 decrease to the same values, showing the preservation of its MZM nature (γ = γ^†)⁴¹. Due to the rotational symmetry of the system, sites near the center j = 15, 17, and 32 are equivalent and they behave very similarly increasing V. In Figs. 1d, e, we show the electron and hole parts of LDOS(ω, j) for site j = 17. As discussed previously, the two multi-site MZMs χ₃ and χ₄ are distributed on sites (j = 15 and 32) and (j = 15, 17, and 32), with total amplitude 2/3 on each site (j = 15, 17, and 32), which leads to a spectral weight 2/3 (compared to the single-site MZMs with weight 1) in the LDOS(ω, j) for site j = 17 (also for j = 15 and 32 at V = 0). Figure 2f, g shows the electron and hole part of LDOS(ω, j) for the edge site j = 31 with increasing V (note: sites j = 1 and 46 are equivalent). The rate of decrease in peak height in electron and hole part of LDOS(ω, j), for multi-site MZMs at site j = 17 and single-site MZM at edge site j = 31 are almost identical when increasing V (in the range V≤2). The LDOS(ω, j) of the local MZM at the central site j = 16 decreases with a slightly faster rate because it develops a finite overlap with χ₃ and χ₄, with increasing V.

The case ϕ ₁ = 0, ϕ ₂ = 0, and ϕ ₃ = 0

Let us consider now the same Y-shape geometry but with the same phase ϕ = 0 on each arm. Surprisingly, for the ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0 case, the pairing term in the Hamiltonian breaks rotational symmetry, as after a 120° anti-clockwise rotation around the central sites l + 1, the pairing term in leg I changes its sign (becomes negative due fermionic anticommutations) (see also Supplementary Note I, for details). Arms II and III have reflection symmetry around the central site l + 1. Using the same transformations as previously discussed, the H^I, H^II and H^III of each arm in terms of Majorana operators can be written in similar form as described using Eqs. (5), (6), (7). Again, in these equations, the Majorana operators ${\gamma }_{A,1}^{I}$, ${\gamma }_{B,2l+1}^{II}$, and ${\gamma }_{B,3l+1}^{III}$ are absent, indicating the presence of three edge single-site MZMs at sites j = 1, 31, and 46 of a Y-junction with a total of 46 sites. The Hamiltonian for the central region H^IV can be transformed as:

$${H}^{IV}=-2i{{\Delta }}\left[{\gamma }_{A,l+1}^{IV}{\gamma }_{B,l}^{I}+\left({\gamma }_{A,l+2}^{II}+{\gamma }_{A,2l+2}^{III}\right){\gamma }_{B,l+1}^{IV}\right].$$

In the above equation, H^IV has a reflection symmetry $\left[{\gamma }_{A,l+2}^{II}\leftrightarrow {\gamma }_{A,2l+2}^{III}\right]$. Thus, defining the operators

$$\begin{array}{lll}{R}_{1}&=&\frac{1}{\sqrt{2}}\left({\gamma }_{A,l+2}^{II}+{\gamma }_{A,2l+2}^{III}\right),\\ {R}_{2}&=&\frac{1}{\sqrt{2}}\left({\gamma }_{A,l+2}^{II}-{\gamma }_{A,2l+2}^{III}\right),\end{array}$$

(13)

the Hamiltonian H^IV further simplifies as

$${H}^{IV}=-2i{{\Delta }}{\gamma }_{A,l+1}^{IV}{\gamma }_{B,l}^{I}-2\sqrt{2}i{{\Delta }}{R}_{1}{\gamma }_{B,l+1}^{IV}.$$

(14)

In Eq. (14) the operator R₂ is absent and has the properties ${{R}_{2}}={R}_{2}^{\,{\dagger} }$, and [H^IV, R₂] = 0, indicating R₂ is a Majorana zero mode. The Majorana zero mode ${R}_{2}=\frac{1}{\sqrt{2}}\left({\gamma }_{A,l+2}^{II}-{\gamma }_{A,2l+2}^{III}\right)$ is equally distributed on sites 17 (II leg) and 32 (III leg) (see Fig. 3a), showing that R₂ is indeed a multi-site MZM. Next, we write Eq. (14) in terms of a 4 × 4 matrix in the basis of $\left[{\gamma }_{A,l+1}^{IV},{R}_{1},{\gamma }_{B,l+1}^{IV},{\gamma }_{B,l}^{I}\right]$ and obtain four eigenvalues ($-\sqrt{2},\sqrt{2},-1,1$). The central region diagonal Hamiltonian can be written in terms of ordinary fermions χ₂ and χ₄ as

$${H}^{IV}=2\sqrt{2}{{\Delta }}\left({{\chi }_{2}}^{{\dagger} }{\chi }_{2}-\frac{1}{2}\right)+2{{\Delta }}\left({{\chi }_{4}}^{{\dagger} }{\chi }_{4}-\frac{1}{2}\right),$$

(15)

with ${\chi }_{2}=\frac{1}{\sqrt{2}}\left(i{R}_{1}+{\gamma }_{B,l+1}^{IV}\right)$ and ${\chi }_{4}=\frac{1}{\sqrt{2}}\left(i{\gamma }_{B,l+1}^{IV}+{\gamma }_{B,l}^{I}\right)$ (see Supplementary Note IV, for more detail). The diagonalized Hamiltonian H^I, H^II, and H^III takes similar form as Eqs. (10), (11), (12), in terms of the ordinary d_k,j fermionic operators. In summary, our analytical calculation finds a total of four MZMs (three single-site at the edge sites and one multi-site MZM near the center region). These four MZMs results in fourfold degeneracy in the ground state of the system, which is also consistent with our full-diagonalization numerical results.

Fig. 3: Majorana zero modes in a Y-shape Kitaev wire, with ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0, and at t_h = Δ = 1.

a In this case, there are a total of four MZMs. Three are localized at the end sites j = 1, 31, and 46. Near the central region, there is one extra multi-site MZM χ₄. The multi-site MZM R₂ is equally distributed among the central sites j = 17 and 32, leading to an spectral weight exactly 1/2 (compared to the localized MZMs with weight 1) as shown in the LDOS(ω = 0, j) on these central sites. For V = 2, the MZMs are spread over a few sites compared to the V = 0 case. The electron part of the LDOS(ω) vs. ω are shown for sites (c) j = 1, (e) j = 17, and (g) j = 31 and the hole part of the LDOS(ω) vs. ω are shown for sites (b) j = 1, (d) j = 17, and (f) j = 31. For V≤2 the electron and hole parts of the LDOS(ω) are almost equal and shows peaks at ω = 0, indicating the presence of stable MZMs on sites j = 1, 17, 31, 32, and 46. Note the weight at j = 17 is half that of sites j = 1 and j = 31.

Figure 3 shows DMRG results for L = 46 sites at t_h = Δ and for different values of the Coulomb interaction V. At V = 0, the site dependent LDOS(ω = 0, j) shows sharp localized peaks at the edge sites j = 1, 31, and 46, indicating three single-site MZMs on those edge sites, as expected. Near the center, the LDOS(ω = 0, j) displays two peaks at sites j = 17 and 32, with height 1/2 compared to the edge sites, suggesting the presence of a multi-site MZM. With increase in interaction (V = 2), these MZMs remain exponentially localized over a few sites (Fig. 3a). To compare the stability of single-site and multi-site MZMs, we calculate the electron and hole parts of LDOS(ω, j) for different values V. The peak values for LDOS^e(ω) (Fig. 3b, f) and LDOS^h(ω) (Fig. 3c, g) at ω = 0, for edge sites j = 1 and 31, decrease to the same values with increase in V. This shows that the characteristic features of the single-site MZM remain and the spectral weight of electron and hole part of LDOS(ω, j) are equal⁴¹, at moderate values of V≤2.

Interestingly, the spectral weight of electron and hole of LDOS(ω) for site j = 17, takes value half (compared to the single-site MZMs on edge sites) at V = 0. This is because the multi-site MZM ${R}_{2}=\frac{1}{\sqrt{2}}\left({\gamma }_{A,17}^{II}-{\gamma }_{A,32}^{III}\right)$ is equally distributed at sites j = 17 and 32. Increasing the repulsion strength V, the peak values of LDOS^e(ω) (Fig. 3d) and LDOS^h(ω) (Fig. 3e) are reduced (but still take the same values). The rate of decrease in peak height for single- and multi-site MZMs are almost the same, which suggests these Majorana modes are equally topologically protected against V.

Figure 4 shows the numerical Bogoliuvov-de Gennes (BdG)^1,2 results for L = 46 sites at V = 0 and for different values of the pairing amplitude Δ, to analyze what occurs away from the sweet spot. At Δ = t_h = 1, the three edge MZMs are single-site localized on one end site of each leg, and do not decay exponentially, whereas the only multi-site MZM (R₂), is localized at sites j = 17 and 32 (with equal peak height in LDOS(ω, j)). At Δ/t_h = 0.6 the four MZMs are quite stable and exponentially decay to a few sites. On the other hand, for even smaller values of Δ/t_h, such as 0.2, the edge sites MZMs decay exponentially over many sites⁴⁵. Note that the multi-site Majorana also decays exponentially over many sites of arms II and III. For that reason, for Δ = 0.2, the multi-site Majorana (R₂) overlaps with the edge MZMs (γ₂ and γ₂) on arms II and III (see Fig 4). The overlap of central MZM to edge MZMs could give rise to interesting features in electron transport experiments³⁵.

Fig. 4: Majorana zero modes in the Y-shape Kitaev wire for different values of Δ/t_h with ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0.

The LDOS(ω = 0, j) vs site j for different values of Δ/t_h. For Δ/t_h = 1.0 and 0.6, the three single-site MZMs are almost localized at the end sites j = 1, 31, and 46, while near the central region, the multi-site Majorana χ₅ is located on sites j = 17 and 32. For Δ/t_h = 0.2, the central MZM (R₂) overlaps with the edge MZMs (γ₂ and γ₃).

The case ϕ ₁ = 0, ϕ ₂ = 0, and ϕ ₃ = π/2

Finally, we consider the Y-shape Kitaev wires with phases ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2. This limit is also equivalent to two perpendicular Kitaev chains with phase difference of π/2 (T-shape wire)^24,46. The Hamiltonian for three legs H^I, H^II, and H^III takes the same form as Eqs. (5), (6), (7), in terms of Majorana operators. As expected, in these equations the Majorana operators ${\gamma }_{A,l}^{I}$, ${\gamma }_{B,2l+1}^{II}$, and ${\gamma }_{B,3l+1}^{III}$ are absent, indicating the presence of three single-site MZMs at the edge sites of the Y-shape Kitaev wire (see Fig. 5). The central region, H^IV, in terms of Majorana operators becomes:

$$\begin{array}{ll}{H}^{IV}=&-\sqrt{2}i{{\Delta }}{\gamma }_{A,2l+2}^{III}\left({\gamma }_{A,l+1}^{IV}+{\gamma }_{B,l+1}^{IV}\right)\\ &-2i{{\Delta }}\left({\gamma }_{A,l+1}^{IV}{\gamma }_{B,l}^{I}+{\gamma }_{A,l+2}^{II}{\gamma }_{B,l+1}^{IV}\right).\end{array}$$

(16)

Equation (16) can be written as a 5 × 5 matrix in the basis $\left[{\gamma }_{A,l+2}^{II},{\gamma }_{A,2l+2}^{III},{\gamma }_{B,l}^{I},{\gamma }_{A,l+1}^{IV},{\gamma }_{B,l+1}^{IV}\right]$. After diagonalizing H^IV, we obtained five eigenvalues $\left(-\sqrt{2},\sqrt{2},-1,1,0\right)$ (see Supplementary Note V, for more details). The last eigenvalue e₅ = 0 and its eigenvector ${\chi }_{5}=-\frac{1}{2}{\gamma }_{A,l+2}^{II}+\frac{1}{\sqrt{2}}{\gamma }_{A,2l+2}^{III}+\frac{1}{2}{\gamma }_{B,l}^{I}$ has the property ${{\chi }_{5}}={\chi }_{5}^{\,{\dagger} }$, and $[{H}^{IV},{\chi }_{5}^{}]=0$, confirming the presence of a spread MZM near the junction. The MZM χ₅ is distributed on three sites j = 15, 17, and 32, showing the multi-site nature of χ₅ (see Fig. 5). The central region H^IV can be written in a diagonal form as:

$${H}^{IV}=2\sqrt{2}{{\Delta }}\left({{\chi }_{2}}^{{\dagger} }{\chi }_{2}-\frac{1}{2}\right)+2{{\Delta }}\left({{\chi }_{4}}^{{\dagger} }{\chi }_{4}-\frac{1}{2}\right),$$

(17)

where ${\chi }_{2}=\frac{i}{2\sqrt{2}}{\gamma }_{A,l+2}^{II}+\frac{i}{2}{\gamma }_{A,2l+2}^{III}-\frac{i}{2\sqrt{2}}{\gamma }_{B,l}^{I}+\frac{1}{2}{\gamma }_{A,l+1}^{IV}+\frac{1}{2}{\gamma }_{B,l+1}^{IV}$ and ${\chi }_{4}=\frac{i}{2}{\gamma }_{A,l+2}^{II}+\frac{i}{2}{\gamma }_{B,l}^{I}-\frac{1}{2}{\gamma }_{A,l+1}^{IV}+\frac{1}{2}{\gamma }_{B,l+1}^{IV}$. Note that the diagonalized system Hamiltonian, in case of the SC phase (ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2) and (ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0) take a similar form (see Eqs. (15) and (17)), which lead to the same energy spectrum for both cases, although the multi-site Majoranas wavefunctions are quite different for these two cases. The remaining Hamiltonians H^I, H^II, and H^III, after diagonalizing in terms of the ordinary d_k,j fermionic operators, take similar forms as Eqs. (10), (11), (12)). In conclusion, our analytical calculations find a total of 4 MZMs. The three single-site MZMs are located at edge sites in their natural positions, while a multi-site MZM χ₅ is situated near the central region. These four MZMs results in fourfold degeneracy in the ground state of the system.

Fig. 5: Majorana zero modes in the Y-shape Kitaev wire with ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2, and at t_h = Δ = 1.

The LDOS(ω = 0, j) shows four Majorana zero modes: (i) three single-site MZMs at the end sites j = 1, 31, and 46, while (ii) near the central region, there is one multi-site Majorana χ₅. For V = 0, the MZM χ₅ is distributed among the central sites j = 15, 17, and 32, with spectral weight 1/4, 1/4, and 1/2, respectively (compared to single-site MZMs with weight 1). For V = 2 the MZMs are spread over more sites compared to the V = 0 case.

In Fig. 5, we present the DMRG calculations with ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2 for different values of V using a system size L = 46. Similarly to the previous cases, the LDOS(ω = 0, j) shows sharp peaks for the edge sites j = 1, 31, and 46, indicating three single-site localized edge MZMs. Interestingly, near the center the LDOS(ω = 0, j) shows three peaks with heights 1/4, 1/4, and 1/2 (compared to the edge sites) at sites j = 15, 17, and 32, respectively. These peaks in the LDOS(ω = 0, j) indicate the presence of a multi-site MZM near the central site. The peak height can be explained by the special form of the multi-site MZM wavefunction ${\chi }_{5}=-\frac{1}{2}{\gamma }_{A,l+2}^{II}+\frac{1}{2}{\gamma }_{B,l}^{I}+\frac{1}{\sqrt{2}}{\gamma }_{A,2l+2}^{III}$, showing that χ₅ is distributed on sites j = 15, 17, and 32 with amplitudes exactly 1/4, 1/4, and 1/2, respectively. Increasing the repulsion V, the peak height of LDOS(ω = 0, j) decreases for different sites and these MZMs become exponentially localized over a few sites (Fig. 5). We find that the peak height of the electron and hole portions of LDOS(ω = 0, j) take the same values even for V = 2, indicating the MZMs are quite stable against repulsive interaction for V≤2.

Discussion

In this publication, we studied the Y-shaped interacting Kitaev chains using analytical and DMRG methods for different superconducting phases at each arm. The Y-shape quantum wires are potentially important to perform braiding experiments in topological quantum computations. We predict the unexpected presence of multi-site Majorana zero modes, in coexistence with the standard single-site Majoranas of the Kitaev model. Based on our analytical and DMRG results: (i) For ϕ₁ = π, ϕ₂ = 0, and ϕ₃ = 0, we predict a total of six MZMs, a large number. There are three single-site MZMs on each edge sites, one single-site MZM at the central site, and also there are two multi-site MZMs near the central region, which results in three peaks at sites l (arm I), l + 1 (arm II), and 2l + 2 (arm III), with heights 2/3 (as compared to the edge sites with height 1), in the site-dependent LDOS calculations. (ii) For ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = 0, we find a total of four MZMs, three single-site localized at the edge sites and one multi-site MZM near the central site. The latter is equally distributed at sites l + 2 (arm II) and 2l + 2 (arm III), leading to two peaks with height 1/2 (compared to the edge sites) as unveiled by the site-dependent LDOS. (iii) For ϕ₁ = 0, ϕ₂ = 0, and ϕ₃ = π/2, we find a total of four MZMs as well: three single-site localized at the edge sites and one multi-site MZM near the center. This multi-site MZM is distributed on sites l (arm I), l + 2 (arm II), and 2l + 2 (arm III), leading to three peaks with heights 1/4, 1/4, and 1/2 in the LDOS calculation.

Our work is motivated by recent progress on quantum dot experiments, where the MZM modes appears at the sweet spot t_h = Δ. Consequently, working in this idealized limit is no longer an abstract idea addressed by theory but it has an associated reality in actual experiments. In this limit for the one-dimensional quantum-dots chain, the edge sites MZMs are fully localized at one site, to be compared to the semiconducting nanowires where Δ ≪ t_h, and the MZMs decay over many sites. For this reason, the Y-shape array of quantum-dots can provide the opportunity to perform braiding experiments using only a few sites, compared to the semiconducting nanowires that need a larger system. As expressed before, we found the surprising result that some MZMs are multi-site near the junction of the Y geometry. This is not the canonical exponential decay when away from the sweet spot or when adding correlations, but the MZM is spread over a small number of sites as in a “box” with sharp boundaries. We found the exotic result that some sites contain 1/2 of a Majorana, some 1/4 of a Majorana, and others 2/3 of a Majorana. This conclusion is in agreement with exact analytical results at the sweet spot and confirmed numerically.

The knowledge of the multi-site Majorana wave-function shape is essential when we exchange the MZMs near the junction in such quantum dot systems. Braiding requires that the Majorana wave functions do not overlap. Near the central region of Y-shape quantum wire, the junction can be also made by three dots mutually coupled to each other in a triangular geometry⁴⁷ (instead of three wires coupled to one central quantum dots). Interestingly, we find that the multi-site Majorana modes will still appear near the junction (with different form). Readers are referred to the Supplementary Note VI, for more detail.

Furthermore, we compare the stability of single- and multi-site MZMs against the repulsive interaction V, by calculating the electron and hole part of LDOS(ω, j) separately. Our DMRG results shows the single-site edge MZMs and multi-site MZMs are equally stable, as the peak values of LDOS(ω, j) reduce with similar rate, for moderate values of repulsive interaction. We also checked the stability of MZMs away from the sweet spot. For Δ/t_h not too different from 1, the MZMs are quite stable and almost localized on the same sites as in the sweet spot. For the smaller values of Δ, as in semiconducting nanowires, the central MZMs decays exponentially at each center over many sites and overlap with the exponentially decaying edge MZMs.

We believe that our finding of multi-site MZMs, and its stability against the Coulomb repulsion and deviation from sweet-spot, will be useful to build fully functional Y-shape junction made from array of quantum-dots^22,23. The multi-site MZMs should be observed in quantum-dots experiments, close to the sweet spots, using just seven quantum dots in a Y-shape geometry in the tunneling-conductance measurements²². In this paper, we primarily focused on finding the physical location of the MZMs in a Y-shape geometry Kitaev chain in the Δ = t_h limit realizable in quantum dots. In the near future, it will be also interesting to study MZMs in the X-shaped Kitaev wire, and analyze the effect of disorder and temperature on these systems. A recent study shows the X-shape wire is also quite important for the braiding process in quantum wires^48,49.

Methods

DMRG method

In order to solve numerically the Y-shaped Kitaev Hamiltonian and measure observables, we have used the density matrix renormalization group (DMRG) method^50,51 with DMRG++³⁹. We performed our DMRG calculations within the two-site DMRG approach, for a system size L = 46 sites and employing m = 1500 states, with truncation error ≤10⁻¹⁰.

Local density-of-states

We have calculated the local density-of-states LDOS(ω, j) as a function of frequency ω and site j, via the Krylov-space correction vector DMRG; for a technical review see⁴⁰. The electron part of the LDOS(ω, j) is⁴¹:

$$LDO{S}^{e}(\omega ,j)=\frac{1}{\pi }Im\left[ \left\langle {\psi }_{0}\left\vert {c}_{j}^{\,{\dagger} }\frac{1}{\omega +H-({E}_{g}-i\eta )}{c}_{j}\right\vert {\psi }_{0} \right\rangle \right],$$

(18)

and the hole part of LDOS(ω, j) is⁴¹:

$$LDO{S}^{h}(\omega )=\frac{-1}{\pi }Im\left[ \left\langle {\psi }_{0}\left\vert {c}_{j}\frac{1}{\omega -H+({E}_{g}-i\eta )}{c}_{j}^{\,{\dagger} }\right\vert {\psi }_{0} \right\rangle \right],$$

(19)

where c_i is the fermionic annihilation operator while ${c}_{j}^{{\dagger} }$ is the creation operator, and E_g is the ground state energy. We use as broadening parameter η = 0.1 as in previous studies^26,52. The total local density-of-states is defined as LDOS(ω, j)= LDOS^e(ω, j) + LDOS^h(ω, j). For the Majorana zero mode, it is expected that the peak values of LDOS^e(ω, j) and LDOS^h(ω, j) be at or very close to ω = 0.