Competing Energy Scales in Topological Superconducting Heterostructures

Artificially engineered topological superconductivity has emerged as a viable route to create Majorana modes. In this context, proximity-induced superconductivity in materials with a sizable spin–orbit coupling has been intensively investigated in recent years. Although there is convincing evidence that superconductivity may indeed be induced, it has been difficult to elucidate its topological nature. Here, we engineer an artificial topological superconductor by progressively introducing superconductivity (Nb), strong spin–orbital coupling (Pt), and topological states (Bi2Te3). Through spectroscopic imaging of superconducting vortices within the bare s-wave superconducting Nb and within proximitized Pt and Bi2Te3 layers, we detect the emergence of a zero-bias peak that is directly linked to the presence of topological surface states. Our results are rationalized in terms of competing energy trends which are found to impose an upper limit to the size of the minigap separating Majorana and trivial modes, its size being ultimately linked to fundamental materials properties.

Line-spectroscopy data acquired across different vortices and at different magnetic fields confirm the existence of a zero-bias-peak emerging with a relatively sharp transition from the superconducting gap. Consistently with the data reported in Fig.3i, these additional data also confirm that the zero-bias-peak does not split right off the vortex center.

MODELING TOPOLOGICAL SUPERCONDUCTING HETEROSTRUCTURES
To understand the physics of unmasked in the vortices measured in proximity-coupled s-wave superconductors -3D time-reversal invariant topological insulators (STI) system, we begin with the Hamiltonian for the 3D topological insulator phase written in the real-space representation as In Eq. (2), α = 1, 2, 3 to denote the appropriate gamma matrix, which we define here as In the definition of the (4 × 4) Gamma matrices, τ α represents the orbital degree of freedom (A, B) and σ α represents the spin degree of freedom (↑, ↓) corresponding to a complete basis that we define as c k = c kA↑ , c kA↓ , c kB↑ , c kB↓ T . In order to be able to coherently interpret the measured results, we utilize both the lattice and continuum forms of Eq. (2) in which certain parameters within the Hamiltonian will acquire different definitions. In particular, in the continuum description, the momentum dependent and the mass term, whose value will dictate the presence of absence of topological states in Eq. (2), is defined as In the preceding definitions presented in Eq. (3) and Eq. (4), the parameters b, m 0 and v F are all materials dependent parameters and may be adjusted to suit the specific properties of the system in question. In all subsequent calculations, we seek to have a qualitative understanding of the physics rather than more specific quantitative agreement, therefore, we set the combination of v F = 1 without any loss in understanding. In a similar fashion to the situation for the continuum case, when we consider the definition of the lattice model that is associated with Eq. (2), we require a commensurate set of definitions for vectors d α and M . In the lattice description of the model, we define these parameters to be while In Eq. (5) and Eq. (6), we require the addition of an additional parameter, a 0 , that represents the lattice constant which we hereafter set to be unity without any loss of generality. The above model possesses time-reversal symmetry that is represented by T = I 2×2 iσ y K in which K represents the operation of complex conjugation. Regardless of referring to the lattice or continuum model, when the mass parameter, m 0 = 0 in either Eq. (4 or Eq. (6), then the Hamiltonian is that of a gapped insulator. With the use of the continuum model, we may, assuming that the system possesses translational invariance, calculate the resultant energy spectrum as In Eq. (7), we note that there are only two branches in the spectrum due to each of the branches being doubly degenerate. In keeping with the convention in dealing with 3D topological insulators such as Bi 2 Se 3 , as is the case here, it is common to set b > 0 thereby delineating the topological and trivial regimes to occur when m 0 > 0 and m 0 < 0 respectively.
In the subsequent analysis, we will keep with the aforementioned parameter convention when describing topological and trivial phases.
In Eq. (8), A is the vector potential and ∆ is the superconducting pairing magnitude which, for s-wave superconductivity, may be written as ∆( r) = ∆ 0 ( r)I ⊗ iσ y where this term effectively pairs opposite spin quasiparticles in the doubled Hamiltonian in basis Ψp = c p c−p † . Given the experimental structure, we assume that the superconductive pairing is homogeneous in the x − y-plane in the absence of vortices and will vary in the direction perpendicular to the the heterostructure, ∆ 0 (z), as there is only superconductivity on one side of the topological insulator. Thus, the superconductivity should be strongest near the interface and decay as we move to the surface, where the measurements are taken.
Considering previous work, we understand that in the presence of a magnetic field, the superconductor will proliferate vortices on the surface that may be accounted for in our theory by a winding of the superconducting order parameter, ∆ 0 (r) = ∆ 0 e iθ(r , where θ(r) where B(r) is the magnetic field, φ 0 is the flux, and λ is the penetration depth. Assuming that we have a vortex positioned within our system at the origin, we may solve the London equation for the magnetic field where K n (x) are modified Bessel functions of the second kind, and the magnetic flux within a corresponding radius of r is given by In both our continuum and numerical lattice models, we consider two different cases of the above defined penetration depth, λ, in our analysis: the "thin-flux" limit and the "thick-flux" limit [5]. Naturally, in a consideration of magnetic penetration in a superconductor, both λ(z), ∆(z) acquire a position dependence that should be taken care of in a self-consistent fashion. In this work, we instead opt to explore the qualitative physics of the STI heterostructure by ignoring the exact functional position dependence of the penetration depth and magnitude of the superconducting order parameter and instead using a presumed dependence for each quantity. Nonetheless, the two limits refer to distinct yet adiabatically connected regions of parameter space. Specifically, the "thin-flux" limit occurs when the topological insulator layer is sufficiently thin and the magnetic flux that penetrates into the bottom layer of the topological insulator does not have sufficient distance to spread before reaching the top surface. Within the parameters that we have defined for the STI system, the "thin-flux" limit occurs when λ ≈ a. In the "thick-flux" limit, that occurs when λ a, the penetrating flux is dispersed uniformly within the TI film prior to reaching the top surface.
In an effort to connect our theoretical framework to the experimental STI, we introduce a single π-flux tube into the STI and proceed by applying k · p perturbation theory. Within perturbation theory, we treat the momentum in theẑ-direction, p z as a perturbation within the basis of wavefunctions corresponding to the zero energy eigenstates, psi 1 0 and psi 2 0 . As seen previously [5], the gapless Hamiltonian for the flux tube is H f lux = p z σ x which corresponds to the insertion of a quantum spin Hall edge state into the STI. We desire to understand the evolution of the behavior of the zero energy eigenstates as we modify λ to increasingly large numbers compared to the lattice constant thereby efficiently traversing the "thin-flux" to "thick-flux" limits. As the penetration depth is increased, we expect the addition of a mass term that monotonically increases with increasing λ. Therefore, the flux tube Hamiltonian including the effects of λ induced masses m x and m y may be written as The resultant energy spectrum for H f lux is E ± p 2 z + m 2 x + m 2 y corresponding to the expected result of a gapped spectrum at p z = 0 in the presence of broken time-reversal symmetry breaking. Clearly the size of the gap that is induced by the presence of the flux, E ind gap = m 2 x + m 2 y is governed by Eq. (11). In order to determine the ultimate size of the gap we E ind gap , we may apply perturbation theory. Nonetheless, in order to apply perturbation theory, we require a linearized version of Eq. (2 in the presence of a magnetic flux, here written where terms that are linear in momentum are kept In Eq. (13), we have explicitly used the magnetic vector potential A = φ r 2 e (−yx+xŷ) written in the symmetric gauge. Given the symmetries of STI system, particularly in the presence of vortices, we my now rewrite the linear Hamiltonian in Eq. (13) in polar coordinates as In Eq. (14), the two momenta components expressed, P θ and P −θ , may be written as and Utilizing Eq. (14) in conjunction with first-order perturbation theory results in a corresponding first-order approximation for the induced gap resulting from the insertion of flux within the continuum model [6] as To complete the picture of the physics of the STI with the continuum approach, we must now expand our Hamiltonian to be able to include the proximity-induced s-wave superconductivity within the simplified flux-line picture of Eq. (11) as The energy spectrum of Eq. (18) consists of four non-degenerate bands, where the spectrum is gapped so long as |∆| < |E ind gap | and represents the solution within the "thin-flux" regime.
However, the continuum model has limitations in that we cannot easily study the "thickflux" limit due to the spatial dependence of the superconducting magnitude, ∆, as the flux spreads out in the topological film. Therefore, we must instead turn to a numerical model based on the lattice representation of the STI in order to be capable of understanding the physics of the "thick-flux" regime and the interpolation between the "thin-flux" and "thick-flux" regimes. We begin the extension of the theoretical model within the lattice representation where we rewrite Eq. (2) in a slightly different manner in order to emphasize the on-site and lattice contributions in the position representation as with H(m) = m 0 − 3b and H(δ) = bΓ 0 +iA δ· Γ 2 [5]. In Eq. (20), r denotes the lattice position, δ = [±a 0x , ±a 0ŷ , ±a 0ẑ ] describes nearest neighbor hopping between lattice sites, and Γ = Γ 1x + Γ 2ŷ + Γ 3ẑ . As we consider the interface between the topological insulator and a type-II s-wave superconductor, we recognize that there must be an even number of vortices present on the surface as a result of the periodic boundary conditions in our numerical solution [5,7,8]. Furthermore, we expect that the coupling between these disparate materials is sufficiently strong so that the superconducting Cooper pairs may tunnel from the superconductor to the topological insulator thereby imparting a superconducting paring po- the phase, φ(r) is a gauge field, however due to the presence of the proliferated vortices, the gauge is not a pure one. Therefore, we perform a gauge transformation to move the phase to the diagonal terms in the Hamiltonian. As the singular gauge transformation considers each of the two vortices on a surface, then the gauge transformation is considered to be "bipartite" [7]. More specifically, the "bipartite" gauge transformation operates on both the particle and hole part of Eq. (21) and its main purpose here is to avoid an unclear multivalued definition for the line integral of the flux [5,7], which must retain path independence up to an integer multiple of 2π. The full BdG Hamiltonian thus becomes phase. Regarding the underlying superconductor, we set the initial pairing magnitude, ∆ = 1 and, we assume that the penetration depth, λ is a parameter that we may set and that it is dependent on the thickness of the topological insulator in theẑ-direction.
With the system parameters defined, we now examine the numerical results in Fig. (7).
Examination of Fig. (7), demonstrates the results of diagonalization of Eq. (22). The solution grid we employ corresponds to a numerical lattice of size n x = 28 points in the x-direction, n y = 20 points in theŷ-direction and n z = 24 points in theẑ-direction with the lattice constant a = 1 for all of the simulations performed in this work. Furthermore, the superconducting order parameter decays away from the surface of the superconductor linearly dropping to zero within 6 layers of the surface of the superconducting interface. We have examined a range of superconducting magnitude decays as a function of position as we move away from the superconductor-topological insulator interface and find that the results presented in Fig. (7) remain largely unchanged.
In Fig. (7)(a), we plot the energies of the two lowest lying states in the STI heterostructure as the penetration depth, λ, varies in magnitude represented in units of the lattice constant, a. The range of λ considered sweeps the range of the "thin-flux" limit to the "thick-flux" limit and, correspondingly, examine the changes in the discrepancy between the lowest and second lowest energies in the STI heterostructure. Clearly, as the penetration depth is increased we move from the "thin-flux" limit to the "thick-flux" limit, in which the Majorana modes are localized on the surfaces of the STI [9]. Fig. (7)(b) illustrates the wavefunction for the lowest lying state in the STI heterostructure within the "thin-flux" limit where the physics is dominated by the wormhole effect that allows Majorana states on the top surface be delocalized between the surfaces creating a gap that pushes the energy of the lowest lying state towards the second lowest lying energy state for small λ. As we increase λ we observe that the lowest lying energies in the system are exponentially decreasing as they approach the ideal, or zero energy. In the intermediate figures Fig. (7)(c) and (d) corresponding to the points of λ that interpolate between the "thin-flux" regime and the "thick-flux" regime, we observe that the weight of the wavefunction shifts towards the top and bottom surfaces indicating that the states are becoming more localized along the surfaces of at the top and bottom of the topological film. In Fig. (7)(e), we observe for λ = 1 that the modes are now localized on the surfaces with no weight of the wavefunction remaining in the bulk of the topological material.