Controllable Majorana vortex states in iron-based superconducting nanowires

Abstract To reveal the non-Abelian braiding statistics of Majorana zero modes (MZMs), it is crucial to design a Majorana platform, in which MZMs can be easily manipulated in a broad topological nontrivial parameter space. This is also an essential step to confirm their existence. In this study, we propose an iron-based superconducting nanowire system with Majorana vortex states to satisfy desirable conditions. This system has a radius-induced topological phase transition, giving a lower bound for the nanowire radius. In the topological phase, the iron-based superconducting nanowires have only one pair of MZMs over a wide range of radii, chemical potential and external magnetic fields. The wave function of MZMs has a sizable distribution at the side edge of the nanowires. This property enables the control of the interaction of MZMs in neighboring vortex nanowires and paves the way for Majorana fusion and braiding.

In Sec. II of the main text, we model the topological surface states in iron-based superconducting nanowires using a TI Hamiltonian which can be written as [1] H TI (k) =Aσ xŝ · k +σ z (M − Bk 2 ) where we adopted long-wavelength approximation in the plane of k = (k x , k y ) for subsequent vortex phase transition occurs at ±µ c = ±A |E g,Z /B| ≈ ±30.6 meV [3], that is, the MZMs exist in the range from −µ c to µ c , which could be verified by the topological invariant ν.
While the radius-induced topological phase transition in thin nanowires has been discussed in Sec. II in the main text.
In the rest of this section, we will apply Bessel expansion on the simplified Fe(Se,Te) model to calculate the topological invariant ν, eigen-energies, and eigen-wavefunctions. We need to adopt k-space for z-direction (the direction of vortex line) when solve the topological invariant, while use full real space for eigen-states. Taking the former for example, we rewrite the TI bands part in the Fe(Se,Te) model Eq. (1) into real space cylindrical coordinate system as Inducing the SC term∆e iϕ , we obtain the entire Hamiltonian H S describing an iron-based SC vortex system. To simplify the following calculations, we could reduce the angular dimension using the continuous rotational symmetry [Ĵ z , H cyl ] = 0 in this cylindrical nanowire, with the z-component of total angular momentumĴ z =L z + τ z (ŝ z − 1)/2 where the orbital part L z = −i ∂ ϕ andτ are the Pauli matrices in the particle-hole space. The wavefunctions take the formÛ ϕ} with the total magnetic quantum number j ∈ Z fulfilling the monodromy of wavefunctions and Ψ (j) (r) is a column vector independent of ϕ. So we could take the transformÛ (j) ϕ to reduce ϕ and get an effective Hamiltonian of r and k z for any certain j.
The radial differential operators in H where |J The topological region of this system can be obtained by regarding the superconducting vortex line as a quasi-1D system with particle-hole symmetry, and the corresponding Z 2 topological invariant can be calculated as [4] where H which is enough to isolate two MZMs at opposite ends of the vortex line. We used Kwant code [5] to construct the Hamiltonian and the PFAPACK library [6] to calculate the Pfaffian.

FINITE SIZE EFFECT OF TI'S LATERAL SURFACE STATES
In the radius-induced topological phase transition of the iron-based superconducting nanowire, the gap of the TI lateral surface states plays a key role and we will drive it in this section.
Let us force on the TI bands' Hamiltonian Eq. (1) at the band inversion point k z = π/c where M = M − 4B 3 /c 2 and M B > 0. In polar coordinate system, we could use the rotational symmetry and rewrite the Hamiltonian as with j ∈ Z + 1/2 now, and divide it into two blocks In the j = 1/2 sector, We define the new wavefunctionΨ(r) = √ rΨ(r) so that 2π |Ψ(r)| 2 dr = 1. Then, the corresponding Hamiltonian becomes When the radius r 0 is very large, we anticipate there is a zero-energy edge stateΨ 0 with 1/r → 0, and thus the Hamiltonian retains Assuming φ(r) ∝ e λ(r−r 0 ) , it becomes When η = −1, it has two roots with positive real part Hence, φ(r) ∝ e λ + (r−r 0 ) − e λ − (r−r 0 ) is the radius wavefunction of the edge states.
When the r 0 is small, the perturbation term can no longer be ignored. For simplicity, we use 1/r ≈ 1/r 0 for this exponentially decayed edge stateΨ 0 . Then we have  [7]. In the space of these two states and their hole part when the SC vortex is induced, we will find the size effect induced topological phase transition.
In addition, note that this energy gap and the phase transition is not due to the overlap between the TI surface states of the opposite surfaces. The energy gap due to the surface states overlapping is directly related to the surface states decay length l 0 , which is determined by the bulk gap. According to the Fe(Se,Te) parameters, the band dispersion is W-shaped and the minimum bulk energy gap is at finite k m ≈ M/B, as shown in Fig. 1. We could estimate the Fermi velocity v F = 2B · k m , the bulk gap E g,b = A · k m , and the surface states decay length l 0 = v F /E g,b = 2B/A ≈ 2 (nm) [8]. However, the topological phase transition in the SC nanowire is at r c ≈ 2ξ ≈ 10 (nm). This means that when the phase transition occurs, the surface state peaks at the nanowire edge and decays to zero before reaching r = 0, as shown in Fig. 2. Therefore, the phase transition here is not attributable to the overlap of the TI surface states.

ANALYTIC SOLUTION OF THE ZEEMAN FIELD INDUCED EDGE MZMS
In this section, we are going to solve the MZM localized at the boundary between superconducting topological surface and Zeeman field analytically, and estimate the energy gap between MZMs and excited states by perturbation to show the quantization of excited energies.
in polar coordinate system. The Zeeman field is mapped intoV Z (r) = −V zŝz Θ(r 0 − r) exists in the center of the disk with a radius r 0 corresponding to the radius of the original nanowire.
Outside the Zeeman field is the s-wave SC vortex region with the order parameter simplified as∆ (r)e iϕ = −iŝ y ∆ 0 e iϕ Θ(r − r 0 ) since it's negligible in the center for a large Zeeman field V z ∆ 0 . In summary, the Bogoliubov-de Gennes Hamiltonian of this disk could be written as Again, taking advantage of the rotational symmetry [Ĵ z , H d ] = 0 and using the method introduced in Appendix , we could get the Hamiltonian H (j) d (r, ∂ r ) for certain total magnetic quantum number j. Next we are going to look for the zero-energy wavefunction Ψ 0 (r) = ψ e↑ ψ e↓ ψ h↑ ψ h↓ T at j = 0 in the Zeeman field and SC region respectively, which satisfies the Schrödinger equations Firstly, let us consider the central Zeeman field region r < r 0 . Since the SC term vanishes, the Hamiltonian H with θ = 2 arctan[(V z + µ)/(V z − µ)] indicating the spin polarization direction and ρ 1 = where ξ is the coherence length of SC. For ρ 1 r 1, the approximate formula for m-order Bessel function of imaginary argument I m (ρ 1 r) ≈ e ρ 1 r / √ 2πρ 1 r is independent of m. For V z µ, θ → π/2 so that cos(θ/2) = sin(θ/2). Then the zero-energy wavefunction on the Zeeman field side could be written as for the Zeeman field of sufficient large strength V 2 z ∆ 2 0 + µ 2 and big size r 0 ξ 0 . C 1 , C 2 are two undetermined constants.
Secondly, we will concentrate on the SC region where r > r 0 without Zeeman field. For with While for spin-down components, we could find another zero-energy solution ψ e↓ ψ h↓ . It diverges at r = 0 but is still valid in this central Zeeman split model. Now, matching the wavefunctions at the boundary r = r 0 with ones in the central region, we get the coefficients C 2 = −iC 1 in Eq. (22), and the zero-energy state's wavefunction in SC region for µ = 0 case is Then, we could pursuit for the low-energy excited energies with m = 0 for µ = 0. The m = 0 part in Hamiltonian Eq. (20) can be view as a perturbation and the excited energies could be calculated as E j = Ψ 0 |H (j) |Ψ 0 / Ψ 0 |Ψ 0 . The zero-energy wavefunction has been solved in µ = 0, but we could simplify it before the integration. Note that, leaving from r = r 0 , the wavefunction in Zeeman field region is almost exponentially decay with length ρ −1 1 ≈ ξ 0 ∆ 0 /V z , while in SC region the decay length is approximate to ξ 0 . Since ∆ 0 V z , the main part of integration is in the SC side, and we could ignore the Zeeman field part of the wavefunction in calculation. Furthermore, when r 0 ξ 0 , the inverse proportional factor appeared in spin-down components could be neglected. Under these approximations, we finally get with the dimensionless quantitiesẼ j = E j /∆ 0 andr 0 = r 0 /ξ 0 = r 0 /(πξ). Therefore, the gap between MZMs and the lowest excited energy E 1 , which protects the MZMs' information, is approximately inversely proportional to the radius of Zeeman field r 0 .
For µ = 0 in SC region, refer to the wavefunctions we just solved, we could transform the Hamiltonian in Eq. (20) into the basis , then the solution that converges at infinity could be found via one of the submatrices in the block anti-diagonalized Hamiltonian with ρ 2 = µ/A and J m (N m ) is the Bessel (Neumann) function of n order. C 3 , C 4 are two undetermined constants. Recalling the boundary conditions Ψ 0 (r 0 ) ∝ 1 −i −i 1 T for µ 2 V 2 z − ∆ 2 0 , we could approximate the zero-energy wavefunction as Again, we could calculate the excited energies by perturbation theory for µ = 0 cases.
Here, in the integrand we have used an approximation to turn part of the probability density into a standard exponential shape for simplicity of calculation ( And finally the results of the excited energies are with the dimensionless quantitiesẼ j = E j /∆ 0 ,r 0 = r 0 /ξ 0 = r 0 /(πξ) andμ = µ/∆ 0 . The energy gap E 1 is inversely proportional to the Zeeman field's radius r 0 while approximately inversely proportional to the square of chemical potential µ. Thus, in order to get a large energy gap to reduce the interaction between MZMs and excited states, it is important to control the chemical potential close to the Dirac point in this system and make the cylinder slenderer.

VORTICES IN FINITE SIZE SYSTEMS
In this section, we are going to investigating the free energy of the superconducting nanowire vortex system, and prove that there is indeed a finite range of the external magnetic field that limits the vortices number to one.
Applying a external magnetic field H ext on the superconducting nanowire in the length direction, the first vortex will penetrate the SC when H ext is at the nanowire's lower critical . And further more, we assume that the second vortex will appear at H (n=1) c1 +δH (we use H for magnetic field instead of Hamiltonian in this section). With the change of H ext , the Gibbs free energy of the superconducting system is always continuous, and through that, we could relate H c1 and δH to the vortices' free energies [10].
In general, if there are n vortices penetrating the nanowire, each with a flux Φ, then the Gibbs free energy G n can be written as Here, F n is the Helmholtz free energy of the superconductor containing n vortices and L is the length of nanowire. The energy of n vortices is ∆F , the continuity of Gibbs free energy requires G 0 = G 1 , and this gives the energy of a single vortex At H ext = H (n=1) c1 + δH, from G 1 = G 2 we obtain ∆F 2 − ∆F 1 = (H (n=1) c1 + δH)ΦL/(4π).
Here the energy of two vortices consists of ∆F 2 = 2∆F 1 + F int with F int the interacting term.
Then we have Therefore, there is a finite suitable magnetic field range δH as the positive F int .
Next we will estimate δH in a superconducting nanowire with a radius r 0 [10]. From the Ginzburg-Landau equation and Maxwell equation, we could derive the magnetic field distribution generated by a single vortex h 1 (r), which approximates to the 0-order Hankel function of imaginary argument at ξ < r < λ with λ the penetration depth of magnetic field and ξ the SC coherence length.
h 0 is the characteristic magnetic field strength. And through this, we could calculate the single vortex energy and the interaction energy between two vortices F int = λ 2 h 0 h 1 (r 2 )/2, where h 1 (r 2 ) indicates the magnetic field distribution at the second vortex core generated by the first one. To minimize the repulsion between two vortices as well as the Gibbs free energy G 2 , the distance between the two vortices should be maximized, namely 2r 0 in the cross-section of the nanowire, then Substituting them into Eq. (33) and Eq. (34), and according to the approximate shape of field h 1 (r), we get