Engineering of many-body Majorana states in a topological insulator/s-wave superconductor heterostructure

We study a vortex chain in a thin film of a topological insulator with proximity-induced superconductivity—a promising platform to realize Majorana zero modes (MZMs)—by modeling it as a two-leg Majorana ladder. While each pair of MZMs hybridizes through vortex tunneling, we hereby show that MZMs can be stabilized on the ends of the ladder with the presence of tilted external magnetic field and four-Majorana interaction. Furthermore, a fruitful phase diagram is obtained by controlling the direction of magnetic field and the thickness of the sample. We reveal many-body Majorana states and interaction-induced topological phase transitions and also identify trivial-superconducting and commensurate/incommensurate charge-density-wave states in the phase diagram.

superconducting gap (1 meV), this hybridization completely destroys MZMs. To save MZMs in this experimental setup, first we consider a 1D dense vortex array in the thin TI film and tune the chemical potential right at the Dirac point of the surface modes so that additional chiral symmetry is preserved. The symmetry suppresses the hybridization of MZMs on the same surface to zero. Hence, the interaction of four Majoranas becomes leading order [21][22][23] , so many-body Majorana wavefunctions have to be considered for the full characterization of the system's quantum phases. Furthermore, when the vortex array is tilted by a magnetic field, the Majorana interaction assists a MZM to appear on the end of the vortex array. Such a many-body effect, though it was less investigated previously, not only provides additional degrees of freedom to engineer MZMs but also open an avenue to study interacting topological physics.
In this report, we propose a possible realization of a one-dimensional vortex array in a superconducting TI film device that can be represented by a tilted ladder model of many Majorana fermions associated with the Fu-Kane model 5 , as shown in Fig. 1(b). In this system, various intravotex and intervortex couplings between Majorana fermions are tunable with the control of the chemical potential and the vortex's incline angle by an external magnetic field. Performing the density-matrix-renormalizaion-group (DMRG) calculations [24][25][26][27] , we obtain the many-body ground state of the system and present interacting phase diagrams as a function of these Majorana couplings. The presence of Majorana interaction enlarges the topological region of the Majorana ladder in the phase diagram; it leads to a MZM localized on the end of the ladder even in the presence of the Majorana hybridization.

Results
Experimental setup of a two-leg Majorana ladder. We start from the Fu-Kane heterostructure 5 , which is a 3D strong TI thin film on the top of an s-wave type-II superconductor. In this thin film, both top and bottom TI surfaces exhibit effective time-reversal-symmetric p ± ip superconductivity, via the superconducting proximity effect. Experimentally this setup has been demonstrated in Bi 2 Te 3 thin films grown on a NbSe 2 substrate 14,17 , as shown in Fig. 1(a). With an external magnetic field turned on, vortices are generated on the TI surfaces and each end of the vortices hosts a MZM [28][29][30] . However, the induced superconducting gap on the naked (top) surface is much smaller than the bottom surface in contact with the superconductor 15,17,20 . Furthermore, the MZMs at the two ends of the vortex can tunnel through the vortex line and then hybridize, such that they do not possess zero energy. For this purpose, we consider a tilted magnetic field, which can effectively enlarge the distance between the MZM on the top and bottom surfaces, and weaken the hybridization.
Inspired by the one-dimensional vortex chain with a tilted magnetic field in the copper oxide thin films 31 , we consider a strongly anisotropic vortex array, which turns out to be a one-dimensional stripe along a certain direction determined by an external magnetic field, as shown in Fig. 1(b). With the tilted fields, the MZMs (red dots) at the top and bottom surfaces oppositely shift and form a tilted two-leg ladder, as shown in Fig. 1(c). On the same surface, the wavefunction of the MZM may overlap with its nearest neighbors and contribute to intra-leg hopping γ γ + t j j 1 1 and λ λ + t j j 2 1 for the top and bottom surfaces, respectively. The thickness of the TI determines the coupling between the top and bottom Majoranas along the same vortex line, γ λ t s j j . As the magnetic field is titled, the hybridization between γ j+1 and λ j becomes non-negligible, resulting in γ λ + t d j j 1 . In addition to the single-particle hopping, there exists interaction among the MZMs. Assuming that the tilted angle θ is small enough such that the distance between γ j and λ j is less than that between γ j+1 and λ j , at the leading order we can have interaction The tight-biding model for the Majorana ladder, describing top-chain (γ) and bottom-chain (λ) Majorana fermions coupled by intra-leg tunnelings t 1 γ i γ i+1 and t 2 λ i λ i+1 , inter-leg tunnelings t s γ j λ j and t d λ j γ j+1 , and four-Majorana interaction Uγ j λ j λ j+1 γ j+1 . where L − 1 in the first summation indicates the open boundary condition. One feasible way to control t 1,2 is to adjust the spatial distance between vortices, which can be artificially tuned via the magnitudes of magnetic fields; at the same time, however, the four-Majorana interaction is weaken. To keep the interaction strength, one needs to tune the chemical potential at the surface Dirac point to preserve additional chiral symmetry. The Majorana hybridization on the surface, which is forbidden by the symmetry, vanishes, and the Majorana interaction, which preserves the symmetry, survives [21][22][23] .
On the other hand, in the noninteracting limit, U = 0, the topology of the Majorana ladder is determined by t d /t s . The topological region where MZM reside on the vortex cores can be exactly determined at Here we drop the constant energy shift after the basis transformation. We also consider a grand canonical ensemble such that the filling of fermions changes with the on-site term − s j j . The spinless fermionic Hamiltonian (2) has a similar structure to an interacting Kitaev chain 10,32,33 . However, we should emphasize that the realistic system (a Majorana ladder in a vortex chain) described by our model is fundamentally different from that in the previous study. Our proposed heterostructure provides a very different mechanism of tuning the model parameters, enabling the exploration of a wider phase diagram. For example, the last term, which has the form of nearest-neighbor electronic interaction, is actually determined by the overlap between four Majorana fermions in a plaquette γ λ λ γ + + j j j j 1 1 . Therefore the strength U is related to the sample thickness and distance between two vortices on the same surfaces and can hence be fine tuned (compared with the hardly tunable electronic interaction in solid). To capture the salient physics, below we consider non-negative tight-binding parameters and interaction strength, i.e. ≥ t 0 s d 1,2, , and U ≥ 0, the same intra-leg tunnelings on the top and bottom surfaces t 1 = t 2 , and the inter-leg hopping t s = 1 as the energy unit. Moreover, we are interested in the phases of the entire ladder, which should not be sensitive to the boundary condition, so we neglect the boundary term + U n n 2 ( ) L 1 in Eq. (2) in following calculations. Before getting into the details of the system's phase diagram, we briefly point out that the original Hamiltonian o f E q . ( 1 ) h a s a n o t h e r e q u i v a l e n t f o r m i n t e r m s o f P a u l i m a t r i c e s σ x , y, z a s σ σ σ σσ σσ σ σ = ∑ − − + + + . This Hamiltonian describes a spin chain with transverse Zeeman field (t s ), anisotropic Dzyaloshinskii-Moriya interaction (t 1, 2 ) 34, 35 , and anisotropic exchange interaction (t d and U). Our proposed heterostructure may thus find applications as a test bed to such interesting spin systems.
Phase diagram. The Hamiltonian Eq. (2) is effectively an 1D fermion chain. To study the many-body physics, we implement the DMRG method to perform numerical simulation, and investigate the ground state phase diagram. We compute the energy gap defined as the difference of the ground state energy in the even parity (P = 1) and odd parity (P = −1) , the difference in paired entanglement spectra δε (δε = 0 indicates two-fold degeneracy of entanglement spectrum), charge structure factor S(q) (which indicate the strength of charge density wave with momentum q) and filling n of fermions (see Method: Physical Quantities for the details). There are several distinct phases such as trivial superconducting phase (TvSC), topological Majorana zero mode (MZM), incommensurate charge-density-wave liquid (IDW) and commensurate charge-density-wave insulators (CDWI) depending on the system parameters. We summarize the phase diagram as a function of t d and U and at a variety of t 1, 2 in Fig. 2.
First let us simply consider the t 1,2 = 0 case as the chemical potential is adjusted at the surface Dirac node, i.e. the intra-leg hopping vanish, in Fig. 2(a) 21 . In this case, the vortices on the same TI surfaces far separate in space. The noninteracting MZMs exist as t d > t s . As the interaction is slightly turned on, we found that the ladder is still in the topological phase. Since the bulk-boundary correspondence still holds in interacting class D, a many-body MZM is localized on each end of the ladder for non-zero U 23 . This many-body MZM is adiabatically connected to the single-particle MZM without the interaction. In Fig. 3(a-d), we show details of four physical quantities vs U at fixed t d = 1.2t s to identify the topological phase. The ground state energy in the even parity (P = +1) and odd parity (P = −1) sectors are doubly degenerate, so ΔE = 0 in Fig. 3(a). Furthermore, in Fig. 3(b) δε = 0, which indicates double degeneracy in the entanglement spectra, leads to the topological phase extended from U = 0. Upon increasing U, on the other hand, the smoothly decreasing filling n and the absence of featured peaks in S(q) show no indication of other physical phases.
The MZM phase region is extended as interaction U increases until the phase transition U c 2 [red circles in Fig. 2(a)] to the incommensurate charge-density-wave liquid (IDW). To show the IDW region, we can see the double degeneracy between the even-and odd-parity ground states is clearly lifted 32 (even if the energy gap is quite small) in Fig. 3(a). Meanwhile, as shown in Fig. 3(b) the double degeneracy of entanglement spectrum disappears, i.e. δε > 0. The charge structure factor S(q) shows peaks at the incommensurate wave vector at ≅ q k 2 F , where k F is the Fermi vector. An example can be seen Fig. 4(a) that the charge structure factors S(q) at different interaction strength as t 1,2 = 0 and t d = 1.2t s . In this regime, filling n still decreases smoothly upon increasing U, and the Fermi vector k F as well as the peak locations of S(q) move towards to a larger q. This charge 2k F instability of IDW state is also reminiscent of a similar feature of a Luttinger liquid 36,37 . In Fig. 2(a), the red circles describe the phase boundary between MZM and IDW.
As U increases across the other phase boundary [blue triangles in Fig. 2(a) or blue line in Fig. 3(c)], the system opens a gap and a CDWI is detected. The dominant peak occurs at q = π and the CDW order parameter survives in the thermodynamic limit. At this moment, the filling approaches .  n 0 5 or half-filling. The ground state is parity odd (P = −1) and ΔE ≠ 0. In a classical analogy, electrons are loaded on every other lattice site. The blue triangles (U CDW ) depict the phase boundary between IDW and CDWI. By DMRG, we can distinguish the distinct phases and pin out the phase boundary by observing variations in ΔE, δε and S(π). The phase transition between  IDW and CDW was also discovered theoretically by bond entropy method 38 . Experimentally, the appearance of CDWI can be measured by Coulomb drag 39 or by thermodynamics method 40 .
Next we turn to the t d < t s regime, which physically corresponds to small tilted angle θ in Fig. 1(b). We used t d = 0.8t s as demonstration presented in Fig. 3(e-h). In the noninteracting limit (U = 0), the system is a trivial superconductor (TvSC) with a finite gap, because the Majorana hybridization t s between the top and bottom TI surfaces destroys the topological phase. The ground state is parity even (P = 1) and the entanglement spectrum shows no paired degeneracy, so both δε ∆ ≠ E, 0 at small U in Fig. 3(e,f). However, at a sufficiently strong (but not too strong) interaction strength, the ladder undergoes the topological phase transition at U c 1 , and MZMs emerge at each end of the ladder. The ground state has double degeneracy and the entanglement spectrum appears in pair, i.e. δε ∆ = = E 0. Back to the phase diagram Fig. 2(a), we can clearly see that the topological state is adiabatically connected to the MZM in the t d > t s regime. This MZM is driven by finite interactions, as an interaction-induced topological state. The black squares in Fig. 2(a) describe the phase boundary between TvSC to MZM (our DMRG calculation shows a weak finite-size effect on the phase boundary). Upon increasing interaction strength, the MZM region is enlarged, implying that a moderate interaction stabilizes the topological MZM, even with less tilted magnetic fields. In the large-U side, the ground states are still characterized as the IDW and CDWI, similar to the observation in the t d > t s regime.
Next we move to consider ≠ t 0 1,2 as the chemical potential is not located at the surface Dirac node. In reality, TI materials with chemical potential exactly at the Dirac node have not been discovered, so the intra-leg tunneling between the MZMs is inevitable. Therefore it is important to investigate how MZM responds to finite t 1,2 . The phase diagrams in Fig. 2(b) and (c) consider finite values of t 1,2 . The influence of t 1,2 is remarkable in the interacting Majorana ladder. In Fig. 2(b) using t 1,2 = 0.1t s , it is obvious to see that, compared to (a), where t 1,2 = 0, the MZM regime shrinks. However, this phase still extends to a finite range whereas the IDW regime is enlarged. There exists a critical t d to harbor the interaction-driven MZM, which is ∼ . t t 0 5 d c s . Below this point, the TvSC phase directly turns to the IDW state, and the MZM disappears. Both TvSC and IDW show trivial behavior in the entanglement spectra. To pin out the boundary, indicated by green stars in Fig. 2(b), we examine the gap magnitudes and observe S(q). At TvSC, ∆ ≠ E 0 and no featured peak in S(q), whereas at IDW, ∆  E 0 and S(q) has a peak located at q = 2k F , as shown in Fig. 4(b). We summarize the variation of the physics observables at t d = 0.4t s and t 1,2 = 0.1t s and at variety of U in Fig. 3(i-l). At U = 0, there is an energy gap. The gap ΔE decreases to a small but finite value as U increases to U SL c and remains small in , ΔE rapidly increases and S(π) jumps to a finite value. In the whole range δε ≠ 0, so no MZM exists.
For further stronger intra-leg tunneling, the region of many-body MZMs becomes even smaller. Figure 2(c) shows the case of = .
s . Therefore, the presence of intra-leg tunneling t 1, 2 will corrupt the stability of the many-body MZM. We have numerically estimated that, as = . ⪆ t t t 0 5 s 1 2 , the system no longer supports the interaction-driven MZM. This implies that the chemical potential has to be tuned to close to the surface Dirac node to suppress the intra-leg tunnelings. In experiment, the magnitudes of magnetic fields are required to be appropriately tuned, such that the vortices are away from each other to lower the intra-leg hybridization but close enough to strengthen the four-Majorana interaction.

Discussion
With proper strength of the Majorana interaction, the topological region is tremendously enlarged; by tilting a small angle of magnetic field MZMs appear on the ladder ends, even in the presence of the Majorana hybridization between the top and bottom surfaces. As shown in Fig. 2 the intra-leg tunneling t 1,2 of Majorana Fermions on the same surface shrinks the topological region. To enlarge the region, t 1,2 can be tuned to zero by adjusting the chemical potential right at the surface Dirac node. For the recent experiment of the heterostructure on Bi 2 Se 3 thin films 15  to expect the MZM on the end of vortex array of the naked surface. Our current proposal directly solves one of the major difficulties of the Fu-Kane model: usually the TI film has to be thin enough to induce the superconductivity gap on the naked surface, but such a thin film can lead to the Majorana hybridization, which destroys MZMs. Tilting the magnetic fields can both reduce the hybridization and enhance the interaction and hence rescue MZMs.
To explore other many-body phases, such as IDW and CDWI, one needs other TI materials to provide larger values of U/t s . It is interesting to see the transition between IDW and MZM, which only occurs with Majorana interactions. Such a topological phase transition is beyond the single-particle picture. The IDW state sharing similarities with a Luttinger liquid could be identified using the Coulomb drag measurement 39 .
Although the physics of many-body MZM and its topology has been discussed extensively [41][42][43][44][45][46] , promising platforms for such systems are barely found in the literature. In this report, we have designed a realizable experimental setup to investigate interaction effects on topological states.

Method
Noninteracting Majorana ladder. To determine the topological phase, we solve the ladder Hamiltonian , the noninteracting Majorana Hamiltonian in momentum basis is given by The topology of the Majorana ladder can be characterized by the Pfaffian of the Hamiltonian at k = 0 and π Hence, the topological region is located at > t t d s in this non-interacting system, irrespective of values of t 1 and t 2 . We also expect t 1 and t 2 small enough to keep the system insulating. In the topological phase, the MZMs reside in the vortex cores whereas they vanish in the trivial phase as the magnetic field goes through the TI without tilting. By changing the tilted angle of the magnetic fields, one can manipulate the ratio of t d /t s to trigger a topological transition between trivial and topological phases.
Computational Methods. With a finite interaction in Eq. (2), exact characterization of the ground state is beyond the single-particle picture. Although one can still perform the Hartree-Fock approximation to decouple the interaction term as with χ = + † c c j j j 1 , the mean-field approach neglects the quantum fluctuations and may not accurately capture the ground state properties in one dimension. To study the many-body physics, we implement the density matrix renormalization groups method (DMRG) to the Hamiltonian of Eq. (2). The DMRG method has been shown to be an efficient numerical algorithm to describe one-dimensional correlated systems [24][25][26][27] and has also been successfully applied on interacting systems with Majorana fermions 32,33,47,48 . The approximated ground-state wave function as well as the entanglement spectrum can be easily obtained via iterative numerical renormalization. We set the number of states kept per block up to m = 120 and compare three different sizes L = 200, 400 and 600 to examine finite-size effects. Furthermore, we keep the truncation errors less than 10 −8 . Although the particle number is not conserved here, in the DMRG calculations, we can still employ parity = − ∑ P ( 1) n j j as a good quantum number to label the quantum state. Thus, the ground state energies E 0 (P) and the wave functions ψ P ( ) 0 associated with the specific parity sector P are accessible.
Physical Quantities. The first signature we use to identify the MZMs is a zero energy gap between the lowest even-parity and odd-parity states, where the subscript l represents partially tracing out the degrees of freedom of the right block. The topological phase has two-fold degeneracies of the entire entanglement spectrum. Rather than observing the entanglement spectrum, throughout the main context, we compute the unitless δε defined as to distinguish the topological from trivial phases 49 . The first summation is over the ground states in two parity sectors. In the topological phase, both the ground state and the entanglement spectra are doubly degenerate, so all the paired entanglement spectrum difference ε ε − + ( ) n P n P 1 vanish and δε = 0. This property is robust even in the presence of interaction 33 and easily implemented with numerical simulation.
In the large-U limit, the system is a commensurate charge-density-wave (CDW) insulator, which is topologically trivial since neither the entanglement spectrum nor the ground-state energy shows double degeneracy. A CDW state has electrons residing on every other lattice sites (in a classical picture) to lower the interaction energy + Un n 4 j j 1 , and can hence be characterized by the structure factor of charge-charge correlations The CDW order parameter can be defined as = →∞ O S q L lim ( )/ L CDW . In the thermodynamic limit, a finite O CDW implies the existence of the long-range CDW ordering. At q = π, the CDW ordering is commensurate, labeled as CDWI in the phase diagrams in Fig. 2. For other values of q, it is incommensurate; in the main context, we have q = 2k F instability in the incommensurate charge-density-wave liquid state, where k F is the Fermi wave vector. In the spinless chain, half-filling = .
n 0 5 corresponds to Figure 4(a) and (b) show the (unnormalized) charge structure factor S(q) for (a) t 1,2 = 0, t d = 1.2t s (b) t 1,2 = 0.1t s , t d = 0.4t s . In both figures, at relative weak interaction strength, labeled by the black squares, no peaks are observed. They are located in the MZM and TvSC states in (a) and (b), respectively. At moderate interaction, S(q) develops peaks located roughly at q = 2k F . Upon increasing U, the filling n decreases and approaches to 0.5, and k F moves toward to π/2. This feature reveals the 2k F charge instability and characterizes the IDW state. At q = π, for (a) U = 1.7t s and for (b) U = t s , the ground state is a CDWI, and consistently, .  n 0 5 at this moment.

Estimation of Majorana coupling and interaction.
To estimate the strengths of the physical parameters, we consider the thin film of topological insulator Bi 2 Se 3 on the top of the superconductor NbSe 2 , which is an experimental realization 15 50 . This hybridization leads to non-zero energy Majorana fermions residing on the vortices.
We can estimate the values of t 1 and t 2 based on the parameters of the superconductivity, since in the absence of the superconductivity the Majoranas are delocalized on the top and bottom layers, regardless of superconductor proximity effect 20 , we use the NbSe 2 superconducting gap G SC ~ 1 meV. The estimated values of the intra-leg tunnelings are given by where the distance between two Majoranas d v on the same surface is about 50 nm and the decay length of Majorana hybridization strength on the topological insulator λ M is close to the London penetration depth of the superconductor (40 nm) when the depth is smaller than ν G / F SC 28 . By comparing with t d and t s , t 1 and t 2 can be neglected.
The interaction U for two Majoranas on the top and two Majoranas on the bottom comes from the Coulomb interaction of two electrons (holes), each of which is the overlap between two Majorana wavefunctions; hence, the strength of the interaction U of four Majorana α 1,2,3,4 has been written in the density function ρ of electron (hole) 21