key: cord-1012809-yls7jz3p
authors: Mercado, Alejandro; Sahoo, Sharmistha; Franz, M.
title: High-temperature Majorana zero modes
date: 2021-06-23
journal: Phys Rev Lett
DOI: 10.1103/physrevlett.128.137002
sha: cbd90eedc750cfe6ba7d42a18880f054cbe1e9d6
doc_id: 1012809
cord_uid: yls7jz3p

We employ general arguments and numerical simulations to show that unpaired Majorana zero modes can occur in cores of Abrikosov vortices at the interface between a three-dimensional topological insulator, such as Bi$_2$Se$_3$, and a twisted bilayer of high-$T_c$ cuprate superconductor, such as Bi$_2$Sr$_2$CaCu$_2$O$_{8+delta}$. When the twist angle is close to $45^{rm o}$ the latter forms a fully gapped topological superconductor up to temperatures approaching its native $T_csimeq 90$K. Majorana zero modes in these structures will persist up to unprecedented high temperatures and, depending on the quality of the interface, may be protected by gaps with larger magnitudes than other candidate systems.

Introduction -Unpaired Majorana zero modes (MZMs) continue to capture the imagination of researchers due to their promise in new technologies but also as an outstanding intellectual challenge in the physics of quantum matter. An important component in the ongoing worldwide effort is design and fabrication of new platforms where MZMs occur robustly under a wide range of conditions and at accessible temperatures [1] [2] [3] [4] [5] [6] . In this Letter we discuss a novel platform based on high-T c cuprate superconductors where MZMs could persist up to the native critical temperature of these materials, which can be as high as 90K in Bi-based superconductors such as Bi 2 Sr 2 CaCu 2 O 8+δ (Bi2212). The proposed platform is a variation on the Fu-Kane paradigm [7] which takes advantage of the anomalous nature of the surface state in the 3D strong topological insulator (STI) proximitized by an ordinary superconductor. The Fu-Kane superconductor may have already been realized in the FeSe family of materials [8] [9] [10] [11] although many questions remain, owing primarily to the complexities associated with the unambiguous MZM detection which remains a key outstanding challenge [12] .

An appealing route towards realizing MZMs at higher temperatures is to employ high-T c cuprate superconductors. However, the SC order parameter in all cuprates is known to have a d x 2 −y 2 symmetry [13] which implies Dirac nodes in the low-energy quasiparticle excitation spectrum. As a result MZMs in the proposed setups that involve cuprates [14] [15] [16] [17] would occur on the background of a sea of gapless Dirac excitations. Unambiguous detection and technological applications, however, demand MZMs protected by a gap, preferably with a large amplitude. In this Letter we propose a strategy to circumvent these inherent limitations of cuprates by exploiting the recently discussed novel physics in twisted bilayers. According to the theoretical work [18, 19] a bilayer assembled from two monolayers of high-T c cuprate spontaneously breaks time reversal T and forms a fully gapped superconductor when the twist angle θ is close to 45 o . Bi2212 has been demonstrated to superconduct at high temperature in its monolayer form [20] and is therefore a natural candidate material for such a bilayer. The Tbroken state of the bilayer is topological with Chern number C = 2 or 4 (depending on system parameters) and can be regarded as a d x 2 −y 2 + id xy superconductor, or d + id for short. Lacking significant spin-orbit coupling (SOC), however, the system will not by itself host an MZM in the core of an Abrikosov vortex [21] .

Our idea, illustrated in Fig. 1 , is to use the twisted cuprate bilayer to proximitize a surface of an STI, which supplies the requisite SOC. Owing to the T -breaking nature of the d + id phase and several other factors discussed below, the mathematical description of the system differs significanly from the canonical Fu-Kane construction [7] and the existence of MZMs is not a priori obvious. The two also belong to different symmetry classes under the Teo-Kane topological classification [22] . Nevertheless we demonstrate here that a simple model describing the structure in Fig. 1 hosts an unpaired MZM in the core of an Abrikosov vortex. Given the topological stability of isolated MZMs we argue that their existence is a robust feature of an STI surface proximitized with a d + id SC order, insensitive to microscopic details.

Compared to the STI proximitized using a conventional s-wave SC our proposed setup features several key advantages including potentially larger excitation gaps and higher critical temperatures. In addition, due to its extremely short coherence length ξ, the cuprate bilayer will exhibit very few Caroli-de Gennes-Matricon states [23] in the vortex core thus enhancing the MZM protection from quasiparticle poisoning and decoherence effects.

The model -The low-energy theory of the STI surface proximitized with SC order is given by the Bogoliubov-de Gennes (BdG) Hamiltonian

Here h k = v(σ y k x − σ x k y ) describes the single massless Dirac fermion on the STI surface [24] and h T k = σ y h * −k σ y is its time-reversed counterpart. v is the surface state velocity, σ α are Pauli matrices in spin space and µ denotes the chemical potential. The d + id order parameter is defined by

and is, in this representation, proportional to the unit matrix in spin space. ∆ and ∆ denote the d x 2 −y 2 and d xy components, respectively. The energy eigenvalues can be written as

It is important to note that in view of Eq. (2) the d + id order parameter yields a fully gapped spectrum only when µ = 0. As we shall see having to work at nonzero chemical potential is one reason why solving for the zero mode analytically becomes more difficult than in the corresponding s-wave case where ∆ k is a constant and one can work at µ = 0. We model a vortex in the SC order parameter by writing the Hamiltonian (1) in the real-space representation

where p = −i∇ is the momentum operator and we performed a uniform π/2 rotation around σ z for future convenience. In order to enable analytical progress we focus here on the pure chiral limit ∆ = ∆ which, for a single vortex placed at the origin, yields a system with full rotation symmetry. The gap operator can be constructed as described in Ref. [25] and is given bŷ

Here a 0 is the lattice constant of the underlying lattice model, p ± = p x ± ip y , {a, b} = 1 2 (ab + ba) and θ(r) denotes the phase of the space dependent gap function ∆(r). The double symmetrization is required to ensure that the resulting operator is hermitian and the last term is needed to maintain gauge invariance [25] .

Vortex core zero mode solution -The order parameter in the presence of a single vortex is given by

where (r, ϕ) represent the polar coordinates of vector r and integer n encodes vorticity. For the sake of simplicity we assume ∆ 0 (r) to be constant for all r, although in reality the amplitude is required to vanish at the origin.

To find the zero-energy eigenstate Ψ 0 (r) bound to the vortex we follow the strategy outlined in Ref. [26] . The Hamiltonian (4) obeys the antiunitary particle-hole symmetry ΞHΞ = −H with Ξ = τ y σ y K, where τ α are Pauli matrices in the Nambu space and K denotes complex conjugation. The p-h symmetry dictates that for every eigenstate of H with energy E there exists an eigenstate Ψ = ΞΨ with energy −E. A non-degenerate eigenstate at zero energy must therefore be also an eigenstate of Ξ, that is ΞΨ 0 = λΨ 0 . Furthermore, because Ξ 2 = 1 it follows that |λ| = 1, and hence λ = e iβ is a pure phase. When seeking the zero mode we can always perform a global U(1) rotation Ψ 0 → e −iβ/2 Ψ 0 fixing λ = 1, which we will assume to be the case henceforth.

If we denote Ψ 0 = (u, v) T with u, v two-component spinors in spin space, it is easy to see that the zero-mode condition ΞΨ 0 = Ψ 0 implies that v = iσ y u * . Thus, Ψ 0 is fully determined in terms of only two independent complex components u = (u ↑ , u ↓ ) T which obey the following equation

We solve this equation by passing to polar coordinates and adopting the ansatz [26] u(r, ϕ) = e ilϕ e −iπ/4 χ ↑ (r) e i(π/4+ϕ) χ ↓ (r) (8) with l integer and χ σ (r) assumed real. Recalling that p ± = e ±iϕ (−i∂ r ± ∂ ϕ /r) Eq. (7) becomes

where∆ l,l = e −ilϕ∆ e −il ϕ . Because∆ ∼ e i(n+2)ϕ the requirement that∆ l,l+1 and∆ l+1,l are ϕ-independent implies 2l + 1 = n + 2. It follows that a non-degenerate zero-energy solution can exist only when l = (n + 1)/2, that is, l = 1 for a vortex and l = 0 for an antivortex. Working out the required∆ l,l operators and passing to a dimensionless coordinate ρ = rk F with k F = µ/v Eqs. (9) become, for the n = 1 vortex, with δ = ∆ 0 µa 2 0 /v 2 the dimensionless gap amplitude. It is straightforward to obtain exact numerical solutions of Eqs. (10) for an arbitrary gap parameter δ; an example is given in Fig. 2 . We also show in Supplementary Material (SM) that a simple analytic form

solves Eqs. (10) in the asymptotic region at large ρ for small δ. As seen in Fig. 2 χ ∞ (ρ) actually turns out to be a very good approximation to the full solution for all distances ρ. Similar results can be established for an antivortex (n = −1) and are given in SM along with additional details of the calculation.

Lattice model results -To confirm the existence of the zero mode in a broader context we performed numerical simulations within a minimal lattice model. These do not depend on the assumption of pure chiral limit and in addition provide useful information on the entire vortex spectrum, not just the zero mode.

Although strictly speaking the STI surface cannot be described by a T -invariant 2D lattice model [27, 28] a workaround exists for situations, such as in the presence of an Abrikosov vortex, where T is explicitly broken. As demonstrated in Refs. [29, 30] a simple lattice model can be constructed that faithfully describes the low-energy 1 d+id -wave physics of the STI surface in the presence of various perturbations, such as magnetic or SC coating. The Hamiltonian is given by

where 

and the lattice BdG Hamiltonian H latt k is then constructed following Eq. (1).

To implement the Abrikosov vortex we express H latt k in the real-space representation. We consider a square sample containing L × L lattice sites with open boundary conditions and the vortex placed at the center. In this representation the BdG Hamiltonian becomes a hermitian matrix of size 4L 2 which we diagonalize numerically for system sizes up to L = 59. More details on the model and its real-space implementation in the presence of vortices can be found in Ref. [31] . Fig. 3 shows the local density of states (LDOS)

where the sum is over positive energy eigenvalues n of H latt and Ψ n (r) = (u n (r), v n (r)) T are the corresponding eigenstates. When evaluated at the vortex center ρ(ω, r) shows a pronounced peak near zero energy, Fig. 3a , associated with a pair of eigenvalues ± 0 close to zero energy, separated by a gap from the rest of the spectrum. We identify these as linear combinations of MZMs bound to the vortex core and the egde, respectively. Indeed when decomposed into Ξ eigenstates using the procedure described e.g. in Ref. [32] the corresponding wavefunction amplitudes are localized at the vortex core and the edge, respectively (insets Fig. 3 ). The small energy splitting is due to the finite size of our sample and should decay exponentially with L. This is indeed observed in Fig. 3b where ∆E = 2 0 is plotted on a semilog scale.

Outlook -Our results establish the existence of unpaired MZMs in cores of Abrikosov vortices at the STI surface with an induced order parameter of d x 2 −y 2 + id xy symmetry. The latter may occur in twisted bilayers of high-T c cuprate superconductors such as Bi2212 [18] . For the sake of analytical tractability the MZMs were shown to exist in the very simplest of models. Nevertheless their topological nature ensures stability against arbitrary perturbations that do not close the gap.

While the mechanism leading to the zero mode is ultimately similar to the classic Fu-Kane construction [7] it is worth noting that the two are not adiabatically connected and the existence of the MZM in the d ± id case is by no means obvious. From the standpoint of the general classification of defects [22] the two systems belong to different symmetry classes: the former respects time reversal T (in addition to the p-h symmetry Ξ) while in the latter T is broken. As discussed in SM the T breaking in the present case implies that MZM wavefunctions in the vortex and the antivortex are not simply related to one another, in stark contrast to the Fu-Kane setup. In addition the fact that the d + id gap function vanishes at k = 0 implies the necessity to work at a non-zero chemical potential µ which in turn brings about the oscillatory nature of the MZM wavefunctions highlighted in Fig. 2 .

Because Bi2212 is known to superconduct up to 90K in its monolayer form our proposed construction could open new directions in the pursuit of Majorana fermions at high temperatures. An important measure of the potential usefulness in applications is the size of the minigap ∆ min protecting the Majorana zero mode from thermal excitations. In our proposed setup the minigap is given by the level spacing δE between the Caroli-de Gennes-Matricon vortex core states. In SM we estimate δE 1/ρ(µ)πR 2 v where ρ(ω) = ω/2π 2 v 2 is the TI surface density of states and R v denotes the vortex core radius, related to the coherence length ξ. We show that owing to the small DOS in the TI surface and the short cuprate coherence length ξ the minigap can be large such that only a handfull of bound states appear inside the vortex core. Obtaining a reliable estimate for the minigap size, however, is challenging because it depends sensitively on the quality of the TI/SC interface as well as the minimum SC gap in the twisted bilayer. The latter in turn depends on details of the interlayer coupling at nonzero twist as discussed in recent works [33] [34] [35] [36] [37] . For these reasons quantitatively accurate estimate of the minigap and other relevant parameters will likely require aditional input from experimental studies specifically aimed at targeting the above uncertainties.

Early work on interfacing topological insulators with (untwisted) Bi2212 already produced encouraging results. Mechanically bonded tunnel junctions [38] showed evidence of proximity-induced superconductivity at temperatures up to 80K with an induced SC gap of 13 meV in Bi 2 Se 3 (and 5 meV in Bi 2 Te 3 ), a good fraction of the Bi2212 maximum gap (∼ 45 meV). Similar results were obtained by growing Bi 2 Se 3 films on Bi2212 crystals [39] where angle-resolved photoemission spectroscopy revealed an induced pairing gap ∼ 15 meV. (We note however, that a subsequent study in a similar setting failed to detect a measurable gap [40] ). These results underscore the urgent need for further studies of STI/cuprate heterostructures. We also hope that our predictions will motivate efforts aimed at fabricating twisted cuprate bilayers and interfacing them with materials that exhibit strong spin-orbit coupling, in search for new platforms capable of supporting MZMs at elevated temperatures.

Approximate analytic solution By appealing to the Bessel function identities (∂ x ± m/x)J m (x) = ±J m∓1 (x) it is easy to see that the solution to Eqs. (10) in the gapless limit δ = 0 takes the form

For x 1 the Bessel functions behave as

which indicates that χ 0 (ρ) is not normalizable, just as one would expect of a scattering state in the absence of a gap. Now consider the effect of turning on small nonzero δ. We seek the solution in the form f (ρ)χ 0 (ρ) where f is a slowly-varying envelope function. In the long-distance limit ρ 1 the dominant contribution to the pairing operator in Eqs. (10) will then come from ∂ 2 ρ acting on χ 0 (ρ) while the other terms will be suppressed by powers of 1/ρ. Furthermore, given the asymptotic form Eq. (S2) we may approximate in this limit ∂ 2 ρ J m (ρ) ≈ −J m (ρ) + o(ρ −1 ). For small δ and large ρ Eq. (10) can thus be approximated as

The solution is an exponentially decaying oscillatory function χ ∞ (ρ) given in Eq. (11) . This wavefunction is normalizable for any δ > 0 and represents a legitimate zero-energy eigenstate bound to a vortex defect in a gapped system. Perhaps surprisingly we find that χ ∞ (ρ) is actually very close to the exact solution obtained by numerical integration of full Eqs. (10) at all distances.

For an antivortex we consider Eq. (9) with n = −1 and l = 0. Working out the expression for the required pairing operators∆ 0,1 and∆ 1,0 it becomes

The full solution, once again, must be obtained by numerical integration; an example is shown in Fig. S1 . The long-distance asymptotic solution can be obtain by the same argument as above and reads

As for the vortex we find that when the dimensionless gap δ is small this asymptotic solution turns out to approximate the exact numerical solution very well at all distances (see Fig. S1 for an explicit comparison).

It is interesting to note that the asymptotic solution χ ∞ (ρ) in Eq. (S5) has the same form as the exact zeromode solution for a vortex in the Fu-Kane model at a non-zero chemical potential [26] . In that model, however, δ = ∆ 0 /µ, reflecting the fact that for an s-wave order parameter the gap is simply ∆ 0 . Another difference is that in the Fu-Kane model vortex and antivortex solutions are simply related by time-reversal symmetry. In the present case time reversal relates vortex in a d+id SC to an antivortex in a d−id SC. Therefore, there is no simple symmetry relation between vortex and antivortex solutions in a d + id state, as can be seen e.g. by comparing Eqs. (S5) and (11) or the corresponding figures. 

According to the classic work of Caroli, de Gennes and Matricon [23] a tower of bound states with energy spacing δE ∼ ∆ 2 /E F is present in the vortex core of a conventional s-wave BCS superconductor. We may expect such finite-energy CdGM states to also exists in the present setup in addition to the Majorana zero mode. Because the topological protection depends on the size of the minigap (defined as the energy of the lowest excited state, δE in our serup) it is important to estimate its magnitude.

As noted in Ref. [23] the vortex core can be regarded as a disk of radius R v = cξ in which the SC order has been locally suppressed. Here ξ = v F /π∆ is the BCS coherence length and c a constant of order one. Correspondingly, the core levels can be viewed as electron states of the underlying normal metal confined inside this disk. In a conventional BCS superconductor the number of such states with energy inside the gap ∆ can be estimated as N CdGM ρ(E F )∆πR 2 v , where ρ(E F ) = m/π 2 is the density of states of the underlying normal metal in two dimensions. Assuming equally spaced core states the level spacing becomes

This result agrees with the detailed calculation of Ref. [23] if one takes c = π/ √ 2 2.22. In conventional superconductors E F typically exceeds ∆ by two or three orders of magnitude making δE too small for individual CdGM states to be experimentally resolved. In NbSe 2 , for example, ∆ 1.1 meV while E F 350 meV. This makes the expected level spacing δE 0.003 meV, too small to be resolved by scanning tunneling spectroscopy (or any other probe).

We may similarly estimate the core level spacing in the proximitized TI surface. The key difference here is the linear Dirac dispersion k = v|k| which yields the normal-state density of states ρ(ω) = ω/2π 2 v 2 . Assuming the chemical potential of the surface electrons µ ∆ TI , where ∆ TI is the proximity-induced gap, we may estimate the number of the core states N CdGM ρ(µ)∆ TI πR 2 v . From this, the level spacing can be de-duced as

where we have taken c = π/ √ 2 as before. We note that the proximity-induced gap ∆ TI has canceled out in this expression, the energy spacing being determined by the density of states ρ(µ) and the coherence length ξ. Because the SC order in the TI is induced by the proximity effect we expect that the vortex core size will be controlled by the coherence length ξ of the superconductor. In a pure d-wave SC the coherence length is momentumdependent, reflecting the momentum-dependent structure of the d-wave gap function. In the d + id phase that we assume occurs in the twisted cuprate bilayer the core size will be controled by the twist-induced minimum gap that we denote ∆ topo (the notation here reflects the fact that ∆ topo protects the topological edge modes). We thus take ξ = v F /π∆ topo as reflected in the last expression of Eq. (S7). For sufficiently small surface chemical potential µ (but still large compared to ∆ TI ) Eq. (S7) indicates that minigap δE can be large. As a concrete example we consider µ = 100meV, ∆ topo = 10meV and v/v F = 1, which yields δE ∼ 4meV. While there is a large uncertainty associated with this estimate it nevertheless suggests that a situation with very few CdGM states inside the induced gap ∆ TI can be readily achieved, in contrast to the conventional SC case discussed above. (Note that ∆ TI is some fraction of ∆ topo , depending on the nature of the TI/SC interface.) The above conclusion is in line with the results of our lattice simulation shown in Fig. 3 which indicates a small number of CdGM states inside the gap, with the majority clustered close to ∆ TI .

Magic angles and current-induced topology in twisted nodal superconductors

Josephson effects in twisted nodal superconductors

Emergent interfacial superconductivity between twisted cuprate superconductors

Doping a moiré mott insulator: A t-j model study of twisted cuprates

Doping phase diagram of a hubbard model for twisted bilayer cuprates