https://orcid.org/0000-0003-1537-4788
Abstract
Majorana fermions are spin-1/2 neutral particles that are their own antiparticles; they were initially predicted by Ettore Majorana in particle physics but their observation still remains elusive. The concept of Majorana fermions has been borrowed by condensed matter physics, where, unlike particle physics, Majorana fermions emerge as zero-energy quasiparticles that can be engineered by combining electrons and holes and have therefore been called Majorana zero modes. In this review, we provide a pedagogical explanation of the basic properties of Majorana zero modes in unconventional superconductors and their consequences in experimental observables, putting a special emphasis on the initial theoretical discoveries. In particular, we first show that Majorana zero modes are self-conjugated and emerge as a special type of zero-energy surface Andreev bound states at the boundary of unconventional superconductors. We then explore Majorana zero modes in 1D spin-polarized p-wave superconductors, where we address the formation of topological superconductivity and the physical realization in superconductor–semiconductor hybrids. In this part we highlight that Majorana quasiparticles appear as zero-energy edge states, exhibiting charge neutrality, spin-polarization, and spatial nonlocality as unique properties that can already be seen from their energies and wavefunctions. Next, we discuss the analytically obtained Green’s functions of p-wave superconductors and demonstrate that the emergence of Majorana zero modes is always accompanied by the formation of odd-frequency spin-triplet pairing as a unique result of the self-conjugate nature of Majorana zero modes. We finally address the signatures of Majorana zero modes in tunneling spectroscopy, including the anomalous proximity effect, and the phase-biased Josephson effect.
1. Introduction
Majorana fermions were initially introduced by Ettore Majorana in 1937 as real solutions to the Dirac equation, later known as the Majorana equation, which describes spin-1/2 particles that are identical to their own antiparticles [1]. Since their conception, Majorana fermions have been the subject of intense research in nuclear and particle physics but their detection still remains an open problem [2–6]. Despite the challenges at high energies, the Majorana equation has been shown to naturally appear at low energies in the condensed matter realm, where its solutions, the Majorana fermions, appear instead as emergent quasiparticles composed of electrons and holes in systems with superconducting order. As a matter of fact, quasiparticles in superconductors are described by the Bogoliubov–de Gennes (BdG) equation [7–13], where the inherent particle–hole symmetry makes it acquire the same form as the Majorana equation [14–19]. Majorana quasiparticles in condensed matter emerge as edge or boundary states in a special type of superconductor known as a topological superconductor. In two dimensions, the topological superconducting phase can be created in spinless superconductors with px ± ipy order parameters or pair potentials, which then host Majorana quasiparticles propagating along the boundaries or at zero energy bound to a vortex core [20–30]. In one dimension, the topological superconducting phase appears in spinless (or spin-polarized) p-wave superconductors where Majorana quasiparticles form at the edge at zero energy [31]. At this point, we note that the 1D spinless p-wave superconductor is often referred to as the Kitaev model. Moreover, Majorana states in condensed matter are also called Majorana zero modes (MZMs) or Majorana bound states (MBSs) for obvious reasons in order to distinguish them from their high-energy counterparts. MZMs or MBSs are sometimes called Majorana fermions but this might be an abuse of the language given that these quasiparticles are not exactly those originally thought up by Ettore Majorana in high-energy physics. In contrast to the high-energy Majorana fermions, the low-energy MZMs exhibit yet another distinguishing property that is based on the fact that braiding MZMs implements non-Abelian unitary transformations, revealing their non-Abelian statistics [15,26,29,32–35] that can be used for topologically protected quantum computation [36–40]. MZMs are also special because a pair of them defines a nonlocal fermion, where information is stored nonlocally and is thus immune to local sources of decoherence. Therefore, MZMs are not only of fundamental relevance, characterizing topological superconductivity, but they also exhibit interesting properties useful for possible applications in future quantum technologies.
The origin of MZMs as boundary states in topological superconductors also reflects a less discussed property, that they in fact represent a special type of zero-energy surface Andreev bound state (ZESABS) in unconventional superconductors [41,42]. In contrast to conventional superconductors, where the order parameter is isotropic (s-wave) in space and spin-singlet, the order parameter in an unconventional superconductor is often anisotropic in space, with a nodal structure on the Fermi surface [43]. The research on ZESABSs in unconventional superconductors dates back to the 1980s in the context of superfluid 3He [44,45] and later on the pairing symmetry of high-Tc cuprates [41,46]; for a review on cuprates, see Ref. [47]. Interestingly, one of the authors (Y.T.) has shown that ZESABSs produce a zero-bias conductance peak (ZBCP) in high-Tc cuprates [48] and revealed that tunneling spectroscopy in unconventional superconductors is a phase-sensitive probe [49]. An outstanding feature of these studies is that the conductance formula derived in Refs. [41,48,49] reveals that the ZESABSs induce perfect Andreev reflection, which is the effect giving rise to ZBCPs, pioneering the search of MZMs via conductance signatures in topological superconductors; see the next paragraph. The existence of ZBCPs has been confirmed by many experiments of tunneling spectroscopy in high-Tc cuprates [50–57], which also allowed the establishment of spin-singlet d-wave pairing symmetry in cuprates [41,58,59]. Along these lines, ZBCPs have also been observed via tunnel spectroscopy in CeCoIn5 [60] and PuCoGa5 [61], whose pair potential also has spin-singlet d-wave symmetry. Furthermore, ZESABSs have been detected in grain boundary Josephson junctions made of cuprates [62], with signals leading to an enhancement and also a nonmonotonic temperature dependence in the supercurrents. Interestingly, Josephson junctions formed by d-wave superconductors have also been shown to host nonsinusoidal current–phase curves as a unique phenomenon due to the unconventional superconducting pair symmetry [63–67]; see also Ref. [68].
The formation of ZESABSs in unconventional superconductors is a universal phenomenon that is not restricted to superconductors with a pair potential having a nodal structure and spin-singlet symmetry as discussed in the previous paragraph. In fact, ZESABSs have also been shown to appear in spin-triplet p-wave and spin-triplet f-wave superconductors [25,69–74], where their identification was entirely due to the appearance of robust ZBCPs. In the simplest situation, it has been shown that tunneling into the edge of a 1D spin-polarized p-wave superconductor gives rise to a ZBCP quantized to 2e2/h due to perfect reflection induced by an MZM, signaling the topological superconducting nature of such a state. This conclusion has its origin in the conductance formula derived by Y.T. in Refs. [41,48,49] in the context of d-wave superconductors; see also Refs. [75–77]. These findings enabled the potential of MZMs, and in general of ZESABSs, for distinguishing between spin-singlet and spin-triplet superconductors. Motivated by this fact, Y.T. later obtained a charge transport formula in diffusive normal metal–unconventional superconductor junctions [78,79], which revealed the robust nature of ZBCPs and further helped to clarify the impact of ZESABSs. Furthermore, a zero-energy peak has been shown to appear in the local density of states (LDOS) of the diffusive normal metal (DN) when attached to a spin-triplet superconductor [80], a process that has been attributed to the penetration of a ZESABS generated at the interface and that accompanies the leakage of superconducting correlations into DN. This phenomenon, together with the ZBCP, in spin-triplet superconductors was called the anomalous proximity effect [80–82] to distinguish it from the conventional proximity effect by spin-singlet superconductors, which instead is characterized by a zero-energy dip in the LDOS [83,84].
Further insights into the physical origin of the anomalous proximity effect have been obtained by analyzing the symmetries of the proximity-induced superconducting correlations or Cooper pairs in the DN region. In this respect, it has been identified that the formation of ZESABSs in spin-singlet d-wave superconductors accompanies the emergence of odd-frequency spin-singlet p-wave pairs, which, however, cannot penetrate into the DN part [85–87]. In contrast, in spin-triplet p-wave superconductors, the induced pair correlations acquire odd-frequency spin-triplet s-wave symmetry that is robust against impurity scattering and thus penetrate into DN, whose accompanied zero-energy LDOS peak can be detected as a ZBCP in a tunneling experiment [42,85,88]. We note that the odd-frequency pair symmetry mentioned here corresponds to a symmetry where the pair amplitude (or Cooper pair wavefunction) is an odd function under the exchange of time coordinates of two electrons forming the Cooper pair; see Refs. [42,89–94] for relevant reviews. Since the ZESABSs in 1D spin-polarized p-wave superconductors are MZMs, the formation of odd-frequency spin-triplet s-wave pairing corresponds to a strong signature of Majorana physics and topological superconductivity [42,94]. This intriguing relationship has been demonstrated in junctions based on p-wave superconductors [80,95–98], d-wave superconductors [99], topological insulators [100–113], Weyl semimetals [114,115], semiconductors with Rashba spin–orbit coupling [99,116–118], superconductor junctions with quantum anomalous Hall insulator [119], superconducting doped topological insulator [120], Majorana arrays [121,122], Sachdev–Ye–Kitaev setups [123], and interacting MZMs [124]. Moreover, it has been revealed that unconventional superconductors hosting ZESABSs, such as spin-singlet d-wave or spin-triplet p-wave, can be classified as topological superconductors, with topological invariants that distinguish the absence or presence of the anomalous proximity effect [125]. At the moment, there is a relevant experimental report detecting ZBCPs in CoSi2/TiSi2 heterostructures [126,127], which could be an indicator of the anomalous proximity effect and of odd-frequency spin-triplet s-wave superconducting pairing. Thus, the anomalous proximity effect, via zero-energy LDOS peaks, ZBCPs, and odd-frequency spin-triplet s-wave pairing, can serve as a powerful signature of topological superconductivity and MZMs.
The motivation to discover MZMs is, therefore, enormous, including the characterization of a new state of matter and also their potential for future quantum applications. There is, however, a crucial problem that is based on the fact that intrinsic spin-polarized p-wave superconductors, needed for realizing topological superconductivity and MZMs, are rare in nature. Nonetheless, it has been shown that 1D spin-polarized p-wave superconductivity and MZMs can be engineered by combining conventional spin-singlet s-wave superconductivity, spin–orbit coupling, and a strong magnetic field [29,34,128–133], which represents one of the most studied physical realizations of topological superconductivity; for the realization using topological insulators, see Refs. [134–137], while for using chains of magnetic atoms, see Refs. [138–145]. The simplicity of the recipe has inspired an impressive number of theoretical studies addressing the main properties of MZMs and their experimental signatures, including quantized ZBCPs [146–156], nonsinusoidal current–phase curves [157–165], the anomalous proximity effect, and odd-frequency pairing [96,99,116]. On the experimental side, ZBCPs have been mostly pursued [153,166–174] but the unambiguous detection of MZMs still remains challenging; few experiments have also explored the Josephson effect, which show interesting promising signatures [175–178]. One of the main current problems with conductance experiments is that these systems also host ZBCPs at finite magnetic fields due to topologically trivial zero-energy states [147,154,176,179–214], thus making it difficult to detect MZMs clearly. By now it is well understood that, while ZBCPs indeed represent a necessary condition for MZMs, it is not sufficient and our efforts need to go beyond local conductance signatures [199,202,204,206,215–223]. In this regard, recent advances have reported the fabrication of high-quality samples, where bottom-up engineering of a minimal Kitaev chain is being pursued [224–226] and control over the emergent spin-triplet superconductivity seems to be feasible [227–230]. It is therefore reasonable to think that the unambiguous identification of true MZMs is very likely to occur in the near future.
The aim of this review is to provide a pedagogical discussion about MZMs from the viewpoint of ZESABSs in unconventional superconductors, placing special emphasis on the initial theoretical discoveries. We start by introducing the BdG equations in unconventional superconductors in Section 2 and therein by simple arguments we demonstrate the emergence of ZESABSs and MZMs and their relationship. In Section 3 we discuss the properties of MZMs in 1D p-wave superconductors, highlighting their models, topological properties, and a brief discussion on the physical realization in the semiconductor–superconductor system. In Section 4, we investigate the emergent pair correlations in 1D p-wave superconductors by means of the Green’s function, something that we believe is crucial for understanding MZMs but often ignored. In Sections 5, 6, and 7 we discuss experimental signatures of MZMs in the charge conductance, anomalous proximity effect, and Josephson effect.
We also note that this review does not intend to cover all the exciting developments in the field of topological superconductivity and MZMs. There are brilliant reviews covering different aspects of these topics and we would like the readers to refer to those works. For a review on Majorana fermions in nuclear, particle, and solid state physics, we refer to Ref. [231]. For readers interested in physical implementations of p-wave superconductivity and experimental signatures of MZMs, we refer to Refs. [152–156,232–238]. For details on MZMs in topological superconductors and their topological properties, we refer to Refs. [239,240]. For readers interested in MZMs in topological insulators, magnetic chains, and helical liquids, we refer to Refs. [136,137,241], Ref. [242], and Ref. [243], respectively. For readers interested in modeling MZMs in semiconductor–superconductor structures, we refer to Ref. [244]. For details on the potential application of MZMs in topological quantum computing, we refer the readers to Refs. [36–40,245]. For an extensive discussion on the emergent superconducting pair symmetries in unconventional superconductors and the relationship between topology and MZMs in topological superconductors, we refer to Refs. [42,94].
2. Surface Andreev bound states and Majorana zero modes in unconventional superconductors
In this section, we elucidate the formation of surface Andreev bound states (SABSs) and their relation to MZMs in unconventional superconductors. For this purpose, we first briefly discuss the structure of the pair potential in unconventional superconductors, outline the BdG approach to describing quasiparticles in superconducting systems, and finally, show how SABSs and MZMs emerge in unconventional superconductors.
2.1. Order parameter of unconventional superconductors
Superconductors are characterized by a pair potential or order parameter Δ, which can be understood as a macroscopic wavefunction of Cooper pairs whose symmetries determine the type of superconductor. Generally speaking, the pair potential depends on two spatial coordinates and is a matrix in spin space,
where |$\boldsymbol {r},\boldsymbol {r}^{\prime }$| and ↑, ↓ denote spatial coordinates and spins, respectively. Although Eq. (1) corresponds to the general definition of the pair potential, for pedagogical purposes, we first analyze it in a spatially uniform system, which is standard in the literature. In this case, it is useful to define the center of mass and relative coordinates according to |$\boldsymbol {R}=(\boldsymbol {r} + \boldsymbol {r}^{\prime })/2$|, |$\boldsymbol {x}=\boldsymbol {r}-\boldsymbol {r}^{\prime }$|. Hence, |${\bf \Delta }\left(\boldsymbol {r},\boldsymbol {r}^{\prime } \right)= {\bf \Delta }\left(\boldsymbol {R},\boldsymbol {x} \right)$| does not depend on |$\boldsymbol {R}$|. After Fourier transformation with respect to |$\boldsymbol {x}$|, we can define |${\bf \Delta }\left(\boldsymbol {k}\right)$| with
On one hand, for the spin-singlet pair potential, the following relations are satisfied:
On the other hand, for the spin-triplet pair potential, |${\bf \Delta }\left(\boldsymbol {k} \right)$| can be written as
where |$\mathbf {d}\left({\boldsymbol k} \right) = -\mathbf {d}\left(-{\boldsymbol k} \right)$| represents the vector containing three spin-triplet components, σj is the jth Pauli matrix in spin space, and |$\boldsymbol{\sigma }$| the vector of Pauli matrices. More explicitly, for pairing between opposite spins with spin projection Sz = 0, the components of |${\bf \Delta }\left({\boldsymbol k} \right)$| exhibit the relations given by
Similarly, for equal-spin pairing with Sz = ±1, we obtain
From the above discussion, we conclude that the spin structure of the pair potential gives rise to spin-singlet and spin-triplet superconductors. Furthermore, the dependence on the momentum |$\boldsymbol {k}$| gives rise to distinct spatial symmetries, sometimes referred to as parity. For instance, superconductors with a spin-singlet pair potential that is isotropic in space exhibit s-wave symmetry; these superconductors are also known as conventional superconductors [246]. Moreover, it is also possible to find more exotic dependences, such as spin-singlet d-wave superconductors and spin-triplet p-wave superconductors.
For a Cooper pair formed between two electrons near the Fermi surface, the pair potential |$\Delta (\boldsymbol {k})$| discussed here acquires a simple form in momentum space. For the spin-singlet s-wave, spin-triplet p-wave, and spin-singlet d-wave, respectively, they can be written as
where |$(\hat{k}_{x},\hat{k}_{y})={\boldsymbol k}_{F}/k_{F}$| represent a unit vector around the Fermi surface, with |${\boldsymbol k}_{F}$| the Fermi wave vector and |$k_{F}=|{\boldsymbol k}_{F}|$|. These pair potentials are schematically shown in Fig. 1.
Schematic representation of pair potentials with the following symmetries in momentum space: A: s-wave, B: p-wave, and C: d-wave. The blue and pink colors represent opposite signs in momentum space.
In the case of systems that are not uniform in space, Eq. 1 needs to be described using two spatial coordinates. However, the spin structure for the distinct types of superconductors remains. In fact, for spin-singlet superconductors, the components of the pair potential read
Analogously, for spin-triplet superconductors with opposite-spin pairing with Sz = 0, we obtain
while for equal-spin-triplet superconductors with spin projections Sz = ±1, we find
As for spatially uniform systems, the pair potential for nonuniform systems can be compactly written as
where ψs describes the spin-singlet pair potential, while d represents the vector containing the three spin-triplet components.
Before closing this part, we note that the distinct types of pair potentials above can emerge intrinsically but they can also be engineered by combining conventional superconductors and other materials. While we do not specify the particular generating mechanism of superconductivity, we will try to point out the origin of a given type of superconductivity whenever necessary. We refer to Ref. [43] for a more detailed discussion on unconventional superconductors.
2.2. Bogoliubov–de Gennes equations
To model the superconductors described in the previous subsection, it is common to employ the Bogoliubov–de Gennes equations, which are commonly written in Nambu space and given by [41,247]
where
is the BdG Hamiltonian, with
describing the normal state part, while |$\Delta _{\sigma \sigma ^{\prime }}({\boldsymbol r},{\boldsymbol r}^{\prime })$| represents an anisotropic superconducting pair potential, which is common in unconventional superconductors as seen in the previous section. Moreover, μ is the chemical potential measured from the band bottom, |$U({\boldsymbol r})$| is the one-body potential, m is the effective mass, and uσ and vσ the electron and hole components of the wavefunction, namely, Ψ = (u↑, u↓, v↑, v↓)T written in so-called Nambu space. The system of equations offered by the matrix equations (12) is commonly referred to as the Bogoliubov–de Gennes (BdG) equations.
Under certain circumstances, Eq. (12) can be decomposed into two 2 × 2 matrices. For instance, for spin-singlet superconductors with a pair potential described by Eq. (8), the two blocks of Eq. (12) are given by
It is thus straightforward to write the corresponding BdG equations for spin-singlet s-wave superconductors as
where we have used the fact that s-wave pair potentials are local in space and satisfy |$\Delta (\boldsymbol {r},\boldsymbol {r}^{\prime })=\Delta (\boldsymbol {r})\delta ({\boldsymbol r} - {\boldsymbol r}^{\prime })$|.
Similarly, for mixed-spin-triplet superconductors with a pair potential described by Eq. (9), the two blocks of Eq. (12) are given by
Moreover, for equal-spin-triplet superconductors with a pair potential described by Eq. (10), the two blocks of Eq. (12) can be written as
Then, to describe superconductors with a given symmetry in the pair potential and find their energies and wavefunctions, one needs to solve the equations written above, which is not always a simple task. One evident difficulty is the second derivative appearing in the normal Hamiltonian h given by Eq. (14). Thus, to simplify the discussion about the formation of SABSs, it is convenient to also simplify the BdG equations.
2.3. Andreev equations
To further simplify the BdG equations, here we discuss the case where they can be reduced to be a 2 × 2 matrix and apply the so-called Andreev approximation where the chemical potential is the largest energy scale. This approximation is sometimes referred to as the quasiclassical approximation. Thus, to study spin-singlet, spin-triplet with Sz = 0, and spin-triplet with Sz = ±1, respectively, it is sufficient to consider Eqs. (17) or (18), Eqs. (19) or (20), and Eqs. (21) or (22); see Refs. [41,247]. Then, we start from
where |$\tilde{u}\left(\boldsymbol {r} \right)=u_{\uparrow }\left(\boldsymbol {r} \right)$|, |$\tilde{v}\left(\boldsymbol {r} \right)=v_{\downarrow }\left(\boldsymbol {r} \right)$| for spin-singlet and spin-triplet pairing with Sz = 0, and |$\tilde{u}\left(\boldsymbol {r} \right)=u_{\uparrow }\left(\boldsymbol {r} \right)$|, |$\tilde{v}\left(\boldsymbol {r} \right)=v_{\downarrow }\left(\boldsymbol {r} \right)$| for equal-spin-triplet pairing with Sz = ±1. Here, +(−) corresponds to spin-singlet (spin-triplet) superconductors.
In the quasiclassical description, it is common to consider |$\mu \gg |\Delta (\boldsymbol {r},\boldsymbol {r^{\prime }})|$|, which implies that we can consider a Cooper pair formed between two electrons on the Fermi surface and the pair potential is determined by the direction of the quasiparticle on the Fermi surface. The electron and hole wavefunction components |$\tilde{u}(\boldsymbol {r})$| and |$\tilde{v}(\boldsymbol {r})$| can be expressed by the product of the rapidly oscillating term |$e^{i \boldsymbol {k}_{F} \cdot \boldsymbol {r}}$| and the slowly varying terms |$u(\hat{\boldsymbol {k}},\boldsymbol {r})$| and |$v(\hat{\boldsymbol {k}},\boldsymbol {r})$| as follows:
where |$\hat{\boldsymbol {k}}=\boldsymbol {k}_{F}/|\boldsymbol {k}_{F}|$| and |$|\boldsymbol {k}_{F}|=\sqrt{2m\mu /\hbar ^{2}}$|. In the following, we ignore the second derivative of |$\hat{\Psi }(\hat{\boldsymbol {k}},\boldsymbol {r})$|. To find the integral in Eqs. (23), we use the approximation that the size of the Cooper pair (the effective range of the interelectron potential) is smaller than the length scale of the pair potential [247]. Then, for the first equation in Eqs. (23), we obtain
while for the second equation in Eqs. (23), we get
where |${\boldsymbol x}={\boldsymbol r}-{\boldsymbol r}^{\prime }$|, |${\bf R}=({\boldsymbol r}+{\boldsymbol r}^{\prime })/2$|, and |$\Delta (\hat{\boldsymbol {k}},\boldsymbol {r})$| satisfies
which represents an effective pair potential with momentum |$\hat{\boldsymbol {k}}$| at position |${\boldsymbol r}$|. To arrive at the third equality in Eq. (25), we have used |$\Delta ({\boldsymbol x},{\boldsymbol r}-{\boldsymbol x}/2)\approx \Delta ({\boldsymbol x},{\boldsymbol r})-({\boldsymbol x}/2)\cdot \nabla \Delta ({\boldsymbol x},{\boldsymbol r})\approx \Delta ({\boldsymbol x},{\boldsymbol r})$|, a valid approximation when the size of the Cooper pair is much smaller than the coherence length of the pair potential. Then, Eqs. (23) can be collected into a 2 × 2 matrix for |$\Psi (\hat{\boldsymbol {k}},\boldsymbol {r})$| given by Eq. (24), leading to an eigenvalue matrix equation given by
where vF is the Fermi velocity and we have used the relations |$\Delta (-\hat{\boldsymbol {k}},\boldsymbol {r})=\Delta (\hat{\boldsymbol {k}},\boldsymbol {r})$| for spin-singlet superconductor, while |$\Delta (-\hat{\boldsymbol {k}},\boldsymbol {r})=-\Delta (\hat{\boldsymbol {k}},\boldsymbol {r})$| for a spin-triplet superconductor. Equation (28) represents the so-called Andreev equations, which are simpler than the BdG equations because they involve linear terms in the derivatives of the diagonals in contrast to the second derivatives of the BdG equations in Eqs. (12). As we will see below, this will be useful when finding SABSs in unconventional superconductors.
2.4. Surface Andreev bound state and Majorana bound states in unconventional superconductors
In this part, we employ the methodology presented in the previous two subsections and explore the emergence of SABSs in unconventional superconductors. To characterize the formation of SABSs, we need to investigate the surface or boundary of unconventional superconductors using the BdG or Andreev equations given by Eqs. (12) or (28). For this purpose, we consider a 2D unconventional superconductor junction with a flat surface located at x = 0 and superconductor (S) defined for x > 0; see Fig. 2. Thus, the perpendicular direction along y is translationally invariant. A quasiparticle in S thus exhibits a wavefunction that has the form of Eqs. (24), which, taking into account that momentum parallel to the interface is conserved, can be written as
where |$\hat{\boldsymbol {k}}$| is related to the injection angle as
where kFx = kFcosθ, kFy = kFsinθ, and |$k_{F}=|\boldsymbol {k}_{\rm F}|$|. We note that Eq. (30) implies that functions labeled by |$\hat{\boldsymbol {k}}$| can be labeled by θ. Thus, the electron and hole components of the wavefunction become dependent on θ and x, which means that |$\Psi (\hat{\boldsymbol {r}},x)$| can be written as
with −π/2 < θ < π/2. Then, within the quasiclassical approximation, Ψ(θ, x) satisfies the Andreev equations given by Eqs. (28), which, for x ≠ 0, read
where j = + for σ = 1 and θj = θ, while j = − for σ = −1 and θj = π − θ. Also, Δ(θj, x) represents the order parameter of the unconventional superconductor parametrized by the angle θj, and Θ(x) is the Heaviside step function, being Θ(x) = 1 for x > 0, in S. Now, we need to solve the above linear coupled differential equations. Before doing so, we note that, in general, the spatial dependence of the pair potential should be determined self-consistently, which is known to give reduced values of the pair potential near the interface. However, here we assume a spatially constant pair potential inside the superconductor because it was previously shown that SABSs and also MZMs are not influenced by the spatial depletion of the pair potential [86,87].
A semi-infinite 2D NS junction along x with its interface located at x = 0 and an unconventional superconductor for x > 0. The surface of the junction is shown as a thick vertical line, while the electron-like and hole-like quasiparticles in the superconductor are depicted by solid blue and dotted red arrows and feel pair potentials parametrized by θ+ = θ and θ− = π − θ. The direction of the arrow denotes that of the group velocity of the quasiparticles.
Then, the wavefunction Ψ(θ, x) can be written as
where
and
Then, by introducing convenient parameters Γ+ and Γ− given by
the wavefunction in S can be written as
with k± = kFx ± γ±. The wavefunction Ψ(θ, x) has two terms. The first and second terms represent electron-like and hole-like quasiparticle wavefunctions. Since Imγ+ > 0 and Imγ− < 0 are satisfied, neither term has divergence for x → ∞, as expected. For these reasons, a state at the boundary x = 0 is called a SABS, which is the focus of this part.
To fully characterize the wavefunction Ψ(θ, x) in S, we need to find the coefficients c and d, which can be found by imposing boundary conditions for the wavefunction at x = 0 given by
when searching for a nontrivial solution, we find the condition of the existence of the SABS to be given by
where Γ± are given by Eqs. (37). Now, to obtain the energy of the SABS, we plug Γ± from Eqs. (37), using Eqs. (34), into Eq. (40) and solve for E. When |Δ+| = |Δ−|, where Δ± = Δ(θ±), a simple expression is obtained for the SABSs, which are given by [248]
where Δ0 is the maximum value of the pair potential, δϕ = ϕ− − ϕ+, and exp(iϕ±) = Δ±/|Δ±|. Interestingly, when Δ+ = −Δ− is satisfied, i.e. when
we obtain δϕ = π, which, remarkably, gives rise to SABSs with zero energy Eb = 0. Thus, the condition given by Eq. (42) gives rise to zero-energy SABSs, which we will denote as ZESABSs [41,46]. The ZESABSs have attracted a lot of attention and have become one of the central topics in the context of high-Tc cuprates [46]. Furthermore, the ZESABSs have been observed by tunneling spectroscopy, giving rise to the establishment of the pairing symmetry in high-Tc cuprates as spin-singlet d-wave. [41,48]. In what follows we discuss ZESABSs for unconventional superconductors with specific profiles of their order parameters, describing spin-triplet p-wave and spin-singlet d-wave pair potentials.
2.4.1. ZESABSs in d-wave superconductors
In the case of spin-singlet d-wave superconductors, we would like to remind the reader that d-wave pairing can be realized in high-Tc cuprates. In this case, the d-wave pair potential can be expressed by the angle measured from the a-axis of the basal plane of the cuprate since one of the directions of the lobe of d-wave pairing coincides with that of the a-axis. Then, we consider a junction with a misorientation angle α, where α denotes the angle between the normal to the interface and the a-axis of the d-wave superconductor. Thus, the pair potential Δ± felt by a quasiparticle with an injection angle θ is given by
A schematic of the surface of α = 0 and α = π/4 is shown in Fig. 3. We now focus on when ZESABSs appear, taking into account the condition discussed before, Δ+Δ− < 0. (i) For |α| ≤ π/4, a ZESABS is generated for π/4 − |α| < |θ| < π/4 + |α| [41,46,48,49,249]. (ii) For α = 0, ZESABSs are absent for any θ. (iii) On the other hand, at α = ±π/4, a ZESABS appears for all θ except for the nodal direction θ = 0. Interestingly, the ZESABS here has a dispersionless flat band as a function of momentum kFy = kFsin θ [80]. This flat band ZESABS has a topological invariant defined for fixed kFy [125].
Schematic figure showing the surface of a d-wave superconductor with misorientation angle α: (a) α = 0, (b) α = π/4. The blue and pink colors indicate values of opposite signs in the pair potential.
2.4.2. ZESABSs in spin-triplet p-wave superconductors
In the case of spin-triplet p-wave superconductors, we restrict our analysis first to 1D, where θ = 0. In this case, the pair potential is expressed by Δ+ = −Δ− = Δ0, which is the condition for having a ZESABSs [44,45,72]. Hence, the p-wave superconductor hosts a ZESABS at the boundary. Since these junctions are in 1D, the ZESABS can be seen as an edge state. It is also important to remark that for a spin-triplet px-wave pair potential with either Sz = 0 or Sz = ±1 and where Δ(θ±) is given by Δ(θ±) = ±Δ0cos θ, the ZESABS become flat as a function of ky, as in the case discussed for d-wave superconductors in the previous subsection.
2.4.3. Wavefunctions of the ZESABSs
In this part, we analyze the wavefunctions of the ZESABSs found in the previous two subsections. We first focus on spin-singlet d-wave and spin-triplet p-wave superconductors with Sz = 0 and then analyze the case with Sz = ±1. For the d-wave superconductor, we choose α = π/4 in Eq. (43), which implies that the pair potential is then given by Δ(θ+) = −Δ(θ−) = Δ0sin(2θ). This pair potential describes spin-singlet dxy-wave superconductors.
2D spin-singlet d-wave and spin-triplet p-wave superconductors with Sz = 0: In this case, for the d-wave superconductor, we choose α = π/4 in Eq. (43) and the pair potential is then given by Δ(θ+) = −Δ(θ−) = Δ0sin 2θ. This order parameter describes spin-singlet dxy-wave superconductors. For the spin-triplet superconductor, we consider Sz = 0 with opposite-spin pairing. Then, from Eq. (38), the wavefunctions of the ZESABSs at the boundary of the spin-singlet d-wave and spin-triplet p-wave are given [125], respectively, by
where |$\sigma (\bar{\sigma }$|) denotes the spin direction with ↑(↓), ↓(↑); ↑ corresponds to σs = 1; and ↓ corresponds to σs = −1. The subscript (0) in the wavefunction indicates that it corresponds to Sz = 0 spin projection. Moreover, we note that |$\xi _{d}^{-1}= (2 \Delta _{0} \sin \theta )/(\hbar v_{F})$| represents the superconducting coherence length in the spin-singlet d-wave superconductor, |$\xi _{p}^{-1}= \Delta _{0}/(\hbar v_{F})$| the superconducting coherence length in the spin-triplet p-wave superconductor. At this point, we remember that the wavefunctions of the ZESABSs given by Eqs. (44) are solutions of Eq. (23). Thus, their wavefunctions define the so-called BdG transformation, which diagonalizes the BdG equations in terms of new quasiparticles, also known as Bogoliubov quasiparticles or Bogolons. Thus, it is possible to define the annihilation operator of a Bogoliubov quasiparticle associated with the ZESABS as
where |$\psi _{j}^{(i)}$| is the wavefunction of the ZESABS given by Eqs. (44) and |$\alpha _{k\sigma }=(a_{k,\sigma },a^{\dagger }_{-k,\bar{\sigma }})^{\rm T}$| is the Nambu spinor. Thus, plugging in the expressions given by Eqs. (44), the annihilation operator of the Bogoliubov quasiparticle is written as
where the coherence factors |$u^{\vphantom{\dagger }}_{k_{Fy}}(x)$| and |$v^{\vphantom{\dagger }}_{k_{Fy}}(x)$| are those defined by the electron and hole entries of the wavefunction in Eq. (44). It is now simple to identify that the coherence factors of the wavefunctions in Eq. (44) satisfy |$u^{\vphantom{*}}_{k_{Fy}\sigma }(x)=v_{-k_{Fy}\sigma }^{*}(x)$|. As a result, the Bogoliubov operator satisfies the following expression:
This equation reveals an interesting nontrivial relation between creation and annihilation operators in different spin sectors. Even though Eq. (47) already reflects the self-conjugate nature pointed out for Majorana states in the introduction, it still does not represent quasiparticles with Majorana nature since Eq. (47) indicates the relation between different spin species.
2D spin-triplet p-wave superconductors with Sz = ±1: We now turn our attention to 2D superconductors with spin-triplet px-wave pair potentials having Sz = ±1. In this case, the two electrons forming Cooper pairs have the same spin direction. Having only one spin direction implies that Cooper pairs behave as spin-polarized or spinless ones. The wavefunction of the ZESABS is given by
where kFx is determined from |$k_{F}=\sqrt{k_{Fy}^{2}+k_{Fx}^{2}}$|. As before, we now define a Bogoliubov annihilation operator as
where the coherence factors obey |$u^{\vphantom{*}}_{k_{Fy}\sigma }(x)=v_{-k_{Fy}\sigma }^{*}(x)$|. Therefore, the Bogoliubov operator satisfies the following expression:
which reveals that the ZESABSs in a spin-triplet p-wave superconductor with Sz = ±1 exhibit Majorana character [250] and the quasiparticles described by Eq. (50) can be referred to as MZMs or MBSs. Moreover, it is worth noting that, if the time-reversal symmetry is satisfied, MZMs appear as Kramers pairs [250]. For a fully spin-polarized state with one of the spin degrees of freedom quenched, we can ignore the spin index σ and write
This situation corresponds to a single MZM appearing at the surface of a spin-polarized p-wave superconductor for fixed momentum kFy. The existence of a nondegenerate MZM at fixed kFy has been extensively investigated. In particular, it has been shown that in 1D, where kFy = 0, the 1D spin-polarized (or spinless) p-wave superconductor realizes a topological phase that is characterized by the emergence of MZMs that fulfill the relation given by Eq. (51). In the next section, we provide further details on this matter.
2.5. Summary
In this section, we have first discussed the symmetries of the order parameters in unconventional superconductors and then outlined a brief introduction of the BdG formalism, which is useful for describing quasiparticles in superconductors. Furthermore, we have employed these concepts and methods to demonstrate the formation of ZESABSs in spin-singlet d-wave and spin-triplet p-wave superconductors. We highlighted that ZESABSs in spin-polarized p-wave superconductors fulfill the self-conjugate Majorana condition and thus become MZMs, whose wavefunction is exponentially localized at the interface and decays towards the bulk of the superconductor following an oscillatory profile.
3. Majorana zero modes in p-wave superconductors
In the previous section, we discussed that MZMs represent a special family of ZESABSs and emerge at the edge of a p-wave superconductor. In this part, we explore further the properties of such exotic zero-energy states, by investigating the continuum and tight-binding models of spin-polarized 1D p-wave superconductors. Furthermore, we also address their relation to topology and how they can be engineered in semiconductor–superconductor hybrids.
3.1. Topological superconductivity and MZMs in p-wave superconductors
In this part, we show that spin-polarized p-wave superconductors exhibit a topological phase that is characterized by the emergence of MZMs. We further show that this topological superconducting phase is characterized by a topological invariant. To address these points we explore the continuum and tight-binding descriptions of spin-polarized p-wave superconductors.
3.1.1. Continuum model
We start by exploring the continuum model of 1D spin-polarized p-wave superconductors, which follows the discussion presented in Section 2.1 for unconventional superconductors. We note that the momentum is a good quantum number and the BdG Hamiltonian can be written as
where |$\psi _{k}=(c_{k},c_{-k}^{\dagger })^{\rm T}$| is the Nambu spinor and
where ξk = ℏ2k2/2m − μ, Δk = kΔ is the momentum-dependent p-wave order parameter, and μ the chemical potential. The eigenvalues of H(k) are |$E_{k,\pm }=\pm \sqrt{\xi _{k}^{2}+k^{2}|\Delta |^{2}}$|, which form a gapped spectrum as long as μ ≠ 0. Interestingly, at μ = 0 the spectrum becomes gapless at k = 0, a critical point that separates two important regimes. For μ > 0, the system in the normal state (Δ = 0) has a Fermi surface and is thus metallic, while for μ < 0 the system with Δ = 0 is an insulator without a Fermi surface; see Fig. 4(a). This is of particular relevance because it reveals that H(k) exhibits two phases. The insulating phase (μ < 0) is gapped at Δ = 0 and is adiabatically connected to the spectrum of the gapped superconductor when μ < 0 and Δ ≠ 0; see Fig. 4(b). In contrast, the metallic phase μ > 0 is gapless at Δ = 0 but becomes gapped only at Δ ≠ 0; see Fig. 4(c). Hence, the phases with μ > 0 and μ < 0 cannot be adiabatically connected without exhibiting gap closing and reopening. Therefore, these two phases are topologically distinct. It has been shown that the phase with μ > 0 is topological, and we will explore this more below.
(a) Energy as a function of the momentum of the normal part of the p-wave superconductor for μ = 0.5 meV > 0 (black curve) and μ = −0.5 meV < 0 (brown dashed curve). The dotted gray horizontal line marks the bottom of the band given by μ > 0; for μ < 0 the band is shifted to higher energies. For μ > 0 the system is metallic and exhibits an odd number of pairs of Fermi points, denoted by magenta-filled circles. For μ < 0, no states at the Fermi level are found and the system is then insulating. (b), (c) Eigenvalues of the BdG Hamiltonian for a p-wave superconductor given by Eq. (53) as a function of momentum for μ = −0.5 meV < 0 (b) and μ = 0.5 meV > 0 (c), where the orange and blue curves correspond to Δ = 0 and Δ = 0.3 meV ≠ 0, respectively. A finite Δ opens a gap at the Fermi points.
3.1.2. Tight-binding model
The lattice version of the p-wave superconductor can be simply given by [31]
where cj (|$c_{j}^{\dagger }$|) destroys (creates) a spinless fermionic state at site j, μTB denotes the onsite energies in the tight-binding (TB) model, similar to the chemical potential of the previous part, t is the nearest-neighbor hopping, and ΔTB is the p-wave pair potential. We label the chemical potential and the pair potential with the subscript TB to distinguish them from their corresponding continuum counterparts. In this part we are interested in exploring the two phases identified in the previous subsection and also the impact when the system is of finite length.
To show that this lattice model describes the physics discussed for the bulk Hamiltonian in Eq. (53), we consider Eq. (54) to form a closed loop with periodic boundary conditions such that c1 = cN + 1 and the second sum in Eq. (54) runs up to N. The closed-loop Hamiltonian is then given by |$\tilde{H}=\hat{H}-t[c_{N}^{\dagger }c_{N+1}+c_{N+1}^{\dagger }c_{N}]+\Delta _{\rm TB} c_{N} c_{N+1}+\Delta _{\rm TB}^{*} c_{N+1}^{\dagger } c_{N}^{\dagger }$|. We can then perform a lattice Fourier transformation of |$\tilde{H}$| and obtain
where |$\psi _{k}=(c_{k},c_{-k}^{\dagger })$| and
with |${\boldsymbol h}=(0,\tilde{\Delta }_{k},\tilde{\xi }_{k})$| and |${\boldsymbol \tau }=(\tau _{x},\tau _{y},\tau _{z})$|, |$\tilde{\xi }_{k}=-2t{\rm cos}(ka)-\mu _{\rm TB}$|, |$\tilde{\Delta }_{k}=-2\Delta _{\rm TB}{\rm sin}(ka)$|, with a being the lattice spacing. We can now see that Eq. (56) indeed corresponds to the Hamiltonian describing the bulk Kitaev model given by Eq. (54). For this purpose, we obtain the spectrum of |$\tilde{H}(k)$|, which is given by |$E_{k,\pm }=\pm \sqrt{[2t{\rm cos}(ka)+\mu _{\rm TB}]^{2}+4|\Delta _{\rm TB}|^{2}{\rm sin}^{2}(ka)}$|. The energy spectrum of Ek,± is plotted in Fig. 5. Then, by expanding these eigenvalues for k ∼ 0, it is clear that we recover the eigenvalues and Hamiltonian of the continuum model after an appropriate redefinition of the TB parameters. By analogy to what we discussed in the bulk case, here we also obtain critical points when |$\tilde{\xi }_{k}=0$| and |$\tilde{\Delta }_{k}=0$|. The pair potential vanishes when k = 0 and k = ±π/a. At these momenta, the normal dispersion vanishes when μTB = −2t and μTB = +2t, respectively. These two lines μTB = ±2t thus separate two gapped superconducting phases, |μTB| < |2t| and |μTB| > |2t|, which are connected only by closing and reopening the energy spectrum gap. Before we go further we note that for |μTB| < |2t|, the system is metallic with a Fermi surface and thus corresponds to the metallic phase discussed for the bulk Hamiltonian in the previous section.
Having understood the appearance of two phases that are connected via the closing and reopening of an energy gap, we now explore the tight-binding Hamiltonian given by Eq. (54) under open boundary conditions in these two phases. For this purpose, it is convenient to perform a transformation into a new set of operators given by
where the operators |$\gamma _{j}^{A(B)}$| are also known as Majorana operators. These Majorana operators satisfy
and follow from the Bogoliubov operators obtained in Section 2. Any fermionic operator that satisfies the conditions given by Eqs. (58) is a Majorana operator. It is worth noting that the second expression in Eqs. (58) indeed expresses the essence borrowed from Majorana fermions: a particle created by the operator γ† is identical to its antiparticle created by γ. We note that it is always possible to perform the Majorana transformation given by Eqs. (57); any fermion operator can be defined in terms of two new operators. However, it is worth pointing out that sometimes this only complicates the problem and does not lead to any novel physics. In our current case, we will show that such decomposition leads to interesting physics.
Thus, the Hamiltonian given by Eq. (54) in terms of these new operators reads
where μTB, ω− = ΔTB − t, and ω+ = ΔTB + t represent, respectively, the hopping amplitudes between Majorana operators of the same fermionic site j, between the first Majorana operator of site j with the second Majorana operator of site j + 1, and between the second Majorana operator of site j with the first Majorana operator of site j + 1. In Fig. 6(a) we schematically show Eq. (59), where the dashed blue, solid green, and solid red lines denote μTB, ω−, and ω+, respectively. Depending on the parameters μTB, ΔTB, and t, the Hamiltonian given by Eq. (59) exhibits different interesting properties, which we discuss next. In particular, we illustrate the difference between the topological and trivial phases by looking at two special limits.
(a) Energy as a function of the momentum of the normal part of the p-wave superconductor within the TB description at μTB = 0. The gray dotted horizontal lines mark the points at which the energy gap closes when tuning the chemical potential μTB appropriately. Inside the region between the gray dotted lines, i.e. |μTB| < 2t, the system is metallic and exhibits an odd number of pairs of Fermi points, denoted by magenta filled circles. Outside this region, however, the system is insulating. (b) Eigenvalues of the BdG Hamiltonian for a p-wave superconductor given by Eq. (56) as a function of momentum for μTB = 2.2t > 2t, where the orange and blue curves correspond to ΔTB = 0 and ΔTB = 0.2 meV ≠ 0, respectively. (c) Same as in (b) but with μTB = 0 < 2t; a finite ΔTB opens a gap at k = ±π/2.
![(a) The Kitaev chain with N sites in the fermionic and Majorana basis given by Eq. (59), depicted by boxes and blue gradient filled circles. In this case, μTB, ω−, and ω+ denote, respectively, the hopping amplitudes between Majorana operators of the same fermionic site j, between the first Majorana operator of site j with the second Majorana operator of site j + 1, and between the second Majorana operator of site j with the first Majorana operator of site j + 1. (b) Representation of the Kitaev chain in the trivial phase, where Majorana operators of the same physical site are coupled to form a Dirac fermion as in Eq. (60). (c) Kitaev chain in the topological phase, where two unpaired Majorana operators appear at the end of the chain; see Eq. (61). Adapted from Ref. [251].](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig6.jpeg?Expires=1747856791&Signature=IgkmZ63wKnelRWShp7q5TQ2xEQOIuxU7oXBDFfa51hkXn6u8zokIZ1P7~J1XbbTD6x4G7sz35fk~KTCf~uh9kkDQimPlViiyacaefq89rfVRKxKC5tjQz0hqNCgMXiHsMlyuiZyJunVXTzG9TBDVHAN~sXvkeO-eBNZcWfZinTHWUcwrsPfVq83sog3qryF~INehyOEdCunJgQ5u3fannXM4orRWiJGnskDhNT6o3Lxg8faO~l1k6glWxir~yEeBRklsT9NQW~mf52B6s6iy0QEwyWQTY9LeAwkp43DPkL6FS87q3QS3ZB2V1UDJalWt-mIAeUXdOh9lVibDzd4~0w__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
(a) The Kitaev chain with N sites in the fermionic and Majorana basis given by Eq. (59), depicted by boxes and blue gradient filled circles. In this case, μTB, ω−, and ω+ denote, respectively, the hopping amplitudes between Majorana operators of the same fermionic site j, between the first Majorana operator of site j with the second Majorana operator of site j + 1, and between the second Majorana operator of site j with the first Majorana operator of site j + 1. (b) Representation of the Kitaev chain in the trivial phase, where Majorana operators of the same physical site are coupled to form a Dirac fermion as in Eq. (60). (c) Kitaev chain in the topological phase, where two unpaired Majorana operators appear at the end of the chain; see Eq. (61). Adapted from Ref. [251].
Regime with |μTB| > |2t|: In this case, we can consider t = ΔTB = 0 as an example. Thus, the Hamiltonian given by Eq. (59) reads
The sum in this equation tells us that Majorana operators from the same physical site are paired together to form a regular fermion, as depicted in Fig. 6(b). In this case, the ground state is given by all states being empty (μTB < 0) or occupied (μTB > 0). It is therefore clear that there is not nothing important happening in this regime, which is a consequence of the system being topologically trivial.
Regime with |μTB| < |2t|: To inspect this phase, we consider μTB = 0 and ω− = 0. Then, the Hamiltonian given by Eq. (59) acquires the following form:
which clearly reveals that Majorana operators on adjacent lattice sites are coupled but, interestingly, operators |$\gamma _{1}^{A}$| and |$\gamma _{N}^{B}$| do not appear in the description; see Fig. 6(c). Thus, the Majorana operators |$\gamma _{1}^{A}$| and |$\gamma _{N}^{B}$| represent Majorana modes emerging at zero energy and located at the ends of the chain. It is worth noting that a consequence of |$\gamma _{1}^{A}$| and |$\gamma _{N}^{B}$| not appearing in the Hamiltonian is that they do not couple to any Majorana operator and therefore commute with the Hamiltonian, namely, |$[H,\gamma _{1}^{A}]=[H,\gamma _{N}^{B}]=0$|. This implies that Majorana states with |$\gamma _{1}^{A}$| and |$\gamma _{N}^{B}$| can be created without energy cost and, therefore, the ground state remains at the same energy.
Further insights are obtained by defining new fermionic operators for effectively writing Eq. (61) and other fermionic operators to describe the Majorana pair not present in Eq. (61). We start by defining a new set of fermionic operators that involve only Majorana operators included in the Hamiltonian given by Eq. (61),
which leads to a new form of Eq. (61) that reads
This Hamiltonian can also be written as
where the ground state energy is then given by E0 = −t(N − 1) for t > 0.
In relation to the two unpaired Majorana operators located at the ends of the chain, |$\gamma _{1}^{A}$| and |$\gamma _{N}^{B}$|, they can be fused into a fermion operator with the following transformation:
Interestingly this fermion exhibits a very unusual property: it is nonlocal in space because it has contributions from both ends of the chain. Moreover, as already pointed out, creating such a nonlocal fermion requires zero energy because their composing Majorana operators are absent in Eq. (63). At this point, we remember that |$[H,\gamma _{1}^{A}]=[H,\gamma _{N}^{B}]=0$|, which implies that |$Hf\ket {GS}=E_{0}f\ket {GS}$| and |$Hf^{\dagger }\ket {GS}=E_{0}f^{\dagger }\ket {GS}$|, where |$\ket {GS}$| is the ground state with energy E0. Thus, the ground states of H can be spanned by f and f†. This can be visualizing by writing the occupation number operator |$n=f^{\dagger }f=(1+i\gamma _{1}^{A}\gamma _{N}^{B})/2$| and the fermion parity operator |$P=(-1)^{n}=-i\gamma _{1}^{A}\gamma _{N}^{B}$|. Then, |$\bra {GS}P\ket {GS}=-\bra {GS}[2f^{\dagger }f-1]\ket {GS}=\mp 1$| for occupied (n = 1) or empty (n = 0) ground states. Hence, we can use n to label the ground state, which can be seen to have a two-fold degeneracy arising from the two possible occupancies n = 0, 1. Thus, the two ground states are |$\ket {0}$| and |$\ket {1}$|, which satisfy |$f\ket {0}=0$| and |$\ket {1}=f^{\dagger }\ket {0}$|, with |$\braket {1}{0}=0$|, which reveals that the ground state degeneracy of the Kitaev’s model for |2μTB| < 2t is indeed two-fold. The parity state therefore defines a degenerate two-level quantum system, which is the definition of a quantum bit (qubit) and can be used to encode quantum information [31]. Furthermore, we remember that the number operator defining the parity corresponds to a spatially nonlocal fermion with zero energy. This implies that the parity state cannot be accessed with any local measurement on one of the bound states at the end of the chain, which shows that the information in such a qubit is stored nonlocally and is robust against local sources of decoherence [31]. See Refs. [39,40,238,245] for a detailed discussion on how Majorana operators can be used for designing robust qubits.
We have therefore seen that the two phases |μTB| < |2t| and |μTB| > |2t| are different and can be distinguished because the former phase hosts MZMs at the edges. We also note that these two phases are topologically distinct and can be distinguished by means of topological invariants. In the case of 1D topological superconductors, it has been shown that one can define a topological invariant given by M = (−1)ν, where ν represents the number of pairs of Fermi points in the Brillouin zone in the normal state (Δ = 0) spectrum of the Kitaev model under periodic boundary conditions (56). In this formulation, an odd (even) number of pairs of Fermi points signal the emergence of a topological (trivial) phase. Interestingly, for |μTB| < |2t|, the number of pairs of Fermi points is one and therefore M = −1, signaling that the phase with |μTB| < |2t| is topological. It is also important to note that, for |μTB| < |2t|, the system is metallic with a Fermi surface, as for the continuum model with μ > 0. Before ending this part, we note that MZMs in this case can also be characterized by a Z2 topological index also known as the Majorana number M, which can be calculated as [31]
where A represents the Hamiltonian written as an antisymmetric matrix given by
Here, M = +1 corresponds to topologically trivial states, while M = −1 characterizes the topologically nontrivial state. In Kitaev’s 1D chain, the Majorana number M is given by
with Pf(A)2 = det(A). Therefore, Eq. (68) implies that the topological phase occurs for 2t > |μTB|, given ΔTB ≠ 0. We can hence conclude that the phase where MZMs appear is a topological superconducting phase.
3.2. Physical realization using semiconductor–superconductor hybrids
As we have demonstrated in previous subsections, to realize topological superconductivity and MZMs, it is necessary to have a spin-triplet p-wave superconductor. There are in principle two ways to have this exotic superconductor. One is to find materials with intrinsic spin-triplet p-wave pair correlations, which, despite intense efforts, such as in Refs. [252–263], still represents an ongoing challenge. Another possibility is to engineer such unconventional superconductivity in the laboratory. This has been predicted to be achieved by combining common ingredients, which include conventional spin-singlet s-wave superconductors, spin–orbit coupling (SOC), and magnetic fields. These ingredients have been exploited in semiconductors [29,128–133], topological insulators [134–137,264], and chains of magnetic atoms [138–143]; see also Ref. [34]. In this part we discuss the realization of 1D topological superconductivity and MZMs by combining conventional spin-singlet s-wave superconductors, Rashba spin–orbit coupling, and an external magnetic field. In particular, this part focuses on understanding the continuum model and its lattice representation.
3.2.1. Continuum model
We consider a single channel nanowire in one dimension with SOC and Zeeman interaction, where the radius of the wire is small compared to the Fermi wavelength and there is a single 1D occupied mode. Thus, this 1D nanowire is modeled by
where ψ = (ψ↑, ψ↓)T, with ψσ representing the annihilation operator of an electron at position x with spin σ = ↑, ↓, while H0 is the Hamiltonian density that reads
Here, the first term is the kinetic part with px = −iℏ∂x being the momentum operator, m is the effective electron’s mass in the nanowire, and μ is the chemical potential that determines the filling of the nanowire, measured from the bottom of the band. The second term corresponds to the Rashba SOC Hamiltonian with the spin direction along the y-axis, where αR represents the strength of the Rashba SOC. The last term in Eq. (70) is the Zeeman Hamiltonian with |$B=g\mu _{B}\mathcal {B}/2$| due to an external magnetic field |$\mathcal {B}$| applied here along the nanowire x-axis, where g is the wire’s g-factor and μB the Bohr magneton. Note that the Zeeman field is perpendicular to the SOC axis, which, as we will see later, is crucial for achieving topological superconductivity.
At this point, it is important to understand the ingredients of the Hamiltonian H0 describing the Rashba nanowire given by Eq. (70). For this purpose, we find the eigenvalues and eigenvectors of H0. In terms of plane waves |$\Psi _{k}(x)=\phi \, e^{ikx}$|, we get
where ξk = ℏ2k2/(2m) − μ, and |$\gamma _{k}=(i\alpha _{R} k+B)/\sqrt{B^{2}+\alpha _{R}^{2}k^{2}}$|. To visualize the behavior of the eigenvalues, in Fig. 7(a) we present them as a function of momentum k for distinct values of the parameters. In the absence of SOC and Zeeman energy, the energy bands are two superimposed parabolas that are degenerate in spin; see the gray curve in Fig. 7(a). A finite SOC shifts the energy bands for each spin along the momentum axis such that they become centered at ±ksoc = ±mαR/ℏ2 and by an overall energy |$E_\mathrm{soc}=m\alpha _{R}^{2}/2\hbar ^{2}$|; see the red and green dashed curves in Fig. 7(a). A finite Zeeman field, perpendicular to the SOC, lifts the degeneracy at k = 0 and opens a gap of 2B, as shown by the magenta curves and yellow region in Fig. 7(a). Interestingly, when B > μ, only the lowest band labeled by − is partially occupied and the nanowire behaves as if spin-polarized or spinless. In this regime, the system exhibits two Fermi points inside the gap corresponding to counterpropagating states with different spins, which implies that the spin projection is locked to the momentum. For this reason, this regime is also known as helical and the opened gap at k = 0 is sometimes referred to as a helical gap [180,265,266]. The behavior of these Rashba bands, having an odd number of Fermi points, is similar to what is found when describing the normal bands of the Kitaev model, revealing the potential of the Rashba nanowire for topological superconductivity, which is discussed next.
(a) Energy bands versus momentum k in the normal state of a 1D Rashba semiconductor. The gray curve corresponds to spin degenerate bands for free electrons, where the chemical potential is measured from the bottom of the band. The red and green dashed curves correspond to energy bands at finite Rashba SOC, which split in spin but still remains degenerate at k = 0. Magenta curves show the bands labeled by ± at finite Zeeman field when B > μ, which opens a gap k = 0 indicated by the yellow-shaded region. Inside this region, states at the Fermi level (red-filled circles) exhibit different spin projections and propagate in opposite directions, giving rise to the concept of helicity. (b)–(e) Evolution of the energy versus momentum in the superconducting state for different values of the Zeeman field defined in terms of the critical field |$B_{c}=\sqrt{\mu ^{2}+\Delta ^{2}}$|. (f) Zeeman dependence of the two energy gaps that appear in the energy bands of the Rashba superconductor at finite Zeeman field. The shaded gray (red) regions indicate their trivial (topological) origin.
To see the emergence of topological superconductivity, we now incorporate conventional spin-singlet s-wave superconductivity in the nanowire. In principle, the nanowire is coupled to a superconductor and then electrons in the nanowire feel an effective superconducting pair potential due to the so-called superconducting proximity effect where superconducting correlations leak into the nanowire [267,268]. Thus, the proximity effect involves tunneling between these systems, implicitly requiring a good interface between the nanowire and superconductor. At a phenomenological level, we consider that conventional superconductivity in the nanowire is described by
where Δ = eiφΔ is the spin-singlet s-wave pair potential and φ is the superconducting phase. Even though Δ is complex in general, in what follows we consider it to be real. Thus, taking into account the normal and superconducting parts described by Eqs. (69) and (72), the full system Hamiltonian is written as
This Hamiltonian already provides a good start to analyze the impact of SOC, Zeeman field, and superconductivity. However, to see how a spinless p-wave superconductor can be engineered, it is instructive to write the Hamiltonian given by Eq. (73) in the helical basis that diagonalizes |$\mathcal {H}_{0}$| constructed from Eqs. (71) as follows [132,233]:
where ψ± are operators that annihilate states in the upper/lower bands labeled by ± at momentum k with energy εk, ±. Moreover, ϕ± correspond to the normalized wavefunctions given by Eq. (71). We further decompose Eq. (74) into the two spinor components, which read
By introducing Eqs. (75) into Eq. (73), the full Hamiltonian in this new basis is given by
where
represent the distinct pair potentials appearing in the nanowire due to the interplay of Rashba SOC, Zeeman interaction, and spin-singlet s-wave superconductivity. The first line in Eq. (76) represents the normal Rashba nanowire Hamiltonian, which, as expected, is diagonal in this helical basis. The second and third lines, respectively, correspond to the superconducting part containing terms that describe pairing between states of the same (Δ−−, + +(k)) and different (Δ+−(k)) bands. The band index σ = ± plays the role of pseudospin (or spin broadly speaking) and therefore, having pairing between states of the same band σ, can be seen to be similar to having spin-triplet pairing between electrons of the same spin, namely, spin-polarized pairing. Analogously, pairing between states of different σ can be seen to be a sort of spin-singlet type of pairing in the band index. Another important property of the pair potentials in Eqs. (77) is that they exhibit an interesting dependence on momentum. In fact, the intraband pair potentials Δ−−, + +(k) are odd functions in momentum k, which reflects their p-wave pair symmetry, similar to what we discussed in Section 2.1. In contrast, the interband Δ+− is an even function in momentum and thus has s-wave pair symmetry. We thus see that the Hamiltonian for a nanowire with Rashba SOC, Zeeman interaction, and conventional spin-singlet s-wave superconductivity is able to realize spin-polarized p-wave superconductivity as shown by Eq. (76), which, however, in general, coexists with spin-singlet s-wave superconductivity.
To further inspect the emergence of p-wave superconductivity, we explore the eigenvalues of Eq. (76). For this reason, we first write it in Nambu space as
where the BdG Hamiltonian is given by
Then, by a simple diagonalization, we obtain the eigenvalues of HBdG, which are given by
where εk, ± are given by Eq. (71) and |$\Delta _{\sigma \sigma ^{\prime }}$| by Eqs. (77). In Figs. 7(b)–(e) we show the evolution of the eigenvalues given by Eq. (80) as a function of momentum k for different values of the Zeeman field B. At B = 0, there are four bands with a finite energy gap opened by Δ at the Fermi points kF, ±; see Fig. 7(b). We remind ourselves that kF, ± are found from εk, ± = 0 and given by |$k_{F,\pm }=\sqrt{k_{F}^{2}+2k_\mathrm{soc}^{2}\mp \sqrt{(k_{F}^{2}+2k_\mathrm{soc}^{2})^{2}-k_{F}^{4}+k_{Z}^{4}}}$|, with |$k_{F}=\sqrt{2m\mu /\hbar ^{2}}$|, |$k_{Z}=\sqrt{2mB/\hbar ^{2}}$|, and ksoc = mαR/ℏ2. At finite Zeeman fields, the bands lift the spin degeneracy at high energies and also develop distinct energy gaps at low and high Fermi momenta denoted by Δ1(2); see the orange and black double arrows in Fig. 7(c). Since these distinct gaps are formed at the lower band E− from Eq. (80), they can be roughly defined as
It is thus evident that, because kF, ± depend on the SOC, Zeeman field, and chemical potential, the gaps Δ1(2) will necessarily exhibit a different behavior when varying such parameters. On further increasing the Zeeman field, the low momentum gap decreases such that the spectrum becomes gapless at B = Bc in Fig. 7(d), where
defines a critical field. The high momentum gap Δ2, however, remains finite, with a value that can be roughly constant under strong SOC. For B > Bc, the gap reopens such that for sufficiently large fields the lowest gap is Δ2 and the higher band is no longer within the low-energy description; see Fig. 7(e). To further understand the gap closing and reopening, we plot Δ1, 2 in Fig. 7(f) at strong SOC, which clearly shows the discussed behavior and reveals that Δ1 exhibits roughly a linear dependence on B. This can be further confirmed by evaluating the energy bands at k = 0, which gives
and clearly shows the gap closing at B = Bc. Thus, the Rashba nanowire with Zeeman field and superconductivity hosts two distinct gapped phases that are connected only by a gap closing and reopening. The phase with B > Bc only contains the lowest band, where only pairing between states of the same band is possible and, therefore, behaves as if it was effectively spin-polarized or spinless.
It has been shown that the two phases, with B > Bc and B < Bc, are topologically different, and B = Bc defines a topological phase transition into a topological superconducting phase with MZMs [130,131]. As in Kitaev’s model discussed in the previous subsection, the topological invariant that distinguishes these phases is given by M = (−1)ν, where ν is the number of pairs of Fermi points in the normal dispersion, which is shown in Fig. 7(a). Indeed, the system hosts an odd number of pairs of Fermi points when |μ| < B, which is fulfilled by B > Bc. Thus, the regime with B > Bc is expected to be topological. We anticipate that this topological phase contains MZMs, protected by the effective gap Δeff = Min(Δ1, Δ2), at the wire ends. Above a certain field |$B_c^{(2)}$|, the gap Δeff saturates at Δ2 [Fig. 7(f)], which is strongly dependent on the strength of SOC and approximated as [161] |$\Delta _{2}\approx 2\Delta E_\mathrm{soc}/\sqrt{E_\mathrm{soc}(2E_\mathrm{soc}+\sqrt{B^{2}+4E_\mathrm{soc}^{2}})}$|. In this regime, a single band is left in the description and pairing occurs between states of the same band. In the following, we show that for high Zeeman fields B ≫ Bc, the Hamiltonian |$\mathcal {H}$| given by Eq. (76) describes Kitaev’s model.
3.2.2. Strong magnetic fields and mimicking the Kitaev model
To have further insight into the previous discussion, the system has to exhibit p-wave pairing symmetry according to Kitaev’s model. Thus, it is convenient to project the system Hamiltonian onto the lower band −. This is allowed because for reaching the topological phase one needs a strong Zeeman field, then the upper band + (see Fig. 7) can be removed from the low-energy problem. Therefore, we can write
where the superconducting pairing potential, or order parameter,
has p-wave symmetry. Now, we write the previous Hamiltonian in the BdG form,
where
whose energy spectrum is given by
The previous equation is in essence the energy spectrum of a p-wave superconductor, thus being in concordance with Kitaev’s model described in the previous section. Therefore, it is not a coincidence that the model given by the Hamiltonian (73) indeed describes Majorana-like physics when B > Bc.
From the point of view of topology, the infinite, Zeeman-polarized, single-channel semiconducting nanowire in proximity to a conventional s-wave superconductor, described in this section, belongs to the so-called 1D D class [269], which has an invariant ν′ that may be ν′ = 0 (topologically trivial) or ν′ = 1 (nontrivial). This system undergoes a band topological transition from ν′ = 0 to ν′ = 1 when the Zeeman splitting B, perpendicular to the spin–orbit axis, exceeds a critical value |$B_c=\sqrt{\mu _S^2+\Delta ^2}$|, where μS and Δ are the wire’s Fermi energy and induced gap respectively.
3.2.3. Lattice model
Now we inspect a finite superconducting nanowire with SOC, Zeeman field, and spin-singlet s-wave superconductivity. We thus discretize the Hamiltonian given by Eq. (73) into a spinful 1D tight-binding chain that is given by
where cσn annihilates a fermionic state with spin σ at site n, μ is the chemical potential, t = ℏ2/(2ma2) is the hopping, m is the effective mass, and a is the lattice spacing. In Eq. (89), 〈n, n′〉 indicates hopping between nearest-neighbor sites. Moreover, Δ represents the induced spin-singlet s-wave pair potential, |$t^\mathrm{soc}_{\pm 1}=\pm \alpha _\mathrm{R}/(2a)$| is the SOC hopping, where αR is the SOC strength that defines a SOC length as λsoc = ℏ2/(mαR), |$B=g\mu _B \mathcal {B}/2$| is the Zeeman field due to an external magnetic field B, where g is the g-factor, and |$\mathcal {B}$| is the magnetic field along the wire. The length of the nanowire is given by L = Na, where N is the number of lattice sites.
Here we are interested in exploring the impact of the Zeeman field on the spectrum because we saw before that it induces a topological phase transition at B = Bc. We thus perform a diagonalization of Eq. (89) and plot its spectrum as a function of the Zeeman field in Figs. 8(a)–(c). In the absence of any SOC, a finite pair potential induces a gap in the low-energy spectrum at B = 0, which does not host energy levels in agreement with the Anderson theorem [270]. In fact, induced superconductivity in the absence of SOC and Zeeman field remains with spin-singlet symmetry as no mechanism is present to mix spin states. At finite Zeeman field but without SOC, the energy levels split in spin, where the up and down levels shift by ±B/2, which then produces a reduction of the size of the induced superconducting gap. Thus, different spin components cross zero energy at B = Δ, signaling the closing of the induced superconducting gap, marked by the magenta dashed line. The closing of the superconducting gap can be understood to be a result of Zeeman depairing where spin-singlet Cooper pairs are destroyed. As B acquires values larger than Δ, there is a region of Zeeman fields roughly given by Δ < B < Bc where energy levels possess both spin components, which, however, depends on having a finite μ. For strong B, such as for B > Bc, the energy levels become spin-polarized and develop crossings at zero and finite energies.
Low-energy spectrum of a nanowire with a finite length L and spin-singlet s-wave superconductivity as a function of the Zeeman field. (a) shows the spectrum at zero SOC at L = 1.5 μm, where the dashed green line marks B = Δ. In (b), (c) the spectrum at finite SOC is shown for two values of the length of the wire, where a magenta dashed line marks the topological phase transition B = Bc into a topological phase with Majorana quasiparticles. Here, δ1(2) denote the lowest gaps in the trivial (topological) phases. In (d)–(f) the wavefunctions of the lowest state at high Zeeman fields are shown as a function of space. The values of the Zeeman fields are indicated by dotted vertical lines in (a)–(c). Parameters: Δ = 0.5 meV, μ = 0.5 meV, αR = 40 meV nm.
In the presence of SOC, the spectrum exhibits interesting differences, as we can observe in Figs. 8(b), (c). The first feature is that the gap closing now occurs at B = Bc, where Bc is the critical field defined by Eq. (82), which marks the topological phase transition separating the trivial and topological phases of the Rashba superconductor nanowire. This gap closing is more visible for longer wires; as expected, the condition B = Bc in Eq. (82) corresponds to a bulk system. Second, for B > Bc the gap reopens and the finite SOC moves all the Zeeman crossings to higher energies but, remarkably, leaves only two energy levels that oscillate around zero energy. These zero-energy oscillatory levels are separated by an energy gap δ2 from the quasicontinuum, which is strongly dependent on the size of the SOC. As the wire is made longer, the amplitude of the oscillations is reduced such that these levels acquire zero energy. For L → ∞, δ1(2) take the values of Δ1(2) defined by Eqs. (81). The strong dependence on the length of the wire already suggests an interesting property related to their spatial behavior. Further insights on the spatial dependence of these intriguing zero-energy states is obtained from their wavefunctions Ψ = (u↑, u↓, v↑, v↓), which are plotted in Figs. 8(d)–(f) for a fixed strong Zeeman field. While the energy level closest to zero energy without SOC exhibits homogeneous oscillations in space just like a particle in a box, the wavefunctions of the zero-energy states at finite SOC develop an interesting profile. Their wavefunctions are localized at the edges and decay towards the center of the system in an exponentially oscillatory fashion, with a decay length ξM that is shorter for longer wires [147,149,161,271]. It can also be verified that the charge density |uσ|2 − |vσ|2 exhibits homogeneous oscillations for short wires and vanishes for long wires [272]. These zero-energy energy levels define the emergence of MZMs or MBSs.
3.3. Summary
We have shown that 1D spin-polarized p-wave superconductors exhibit a topological phase where MZMs emerge as zero-energy edge states. Moreover, we have pointed out that MZMs also exhibit spatial nonlocality due to their wavefunction localization at the edges of the wire. Finally, we have shown that 1D spin-polarized p-wave superconductivity can be realized by combining conventional spin-singlet s-wave superconductivity, SOC, and a Zeeman field perpendicular to the SOC axis.
4. Green’s functions of p-wave superconductors with MZMs
Having discussed the emergence of MZMs in 1D p-wave superconductors, here we would like to explore how they manifest on the local density of states (LDOS) and what the nature is of their superconducting pair correlations or simply pair amplitudes. To address these two questions, we calculate the Green’s functions of 1D p-wave superconductors, which permit us to directly obtain the LDOS and pair amplitudes. For completeness, here we first focus on the continuum model of 1D p-wave superconductors and inspect the impact of having semi-infinite and finite length systems. Furthermore, for pedagogical purposes, we start by discussing how to obtain the Green’s functions by using scattering states.
4.1. Brief discussion on obtaining the Green’s functions from scattering states
To obtain the Green’s functions, in this review, we follow the scattering formalism initially formulated by McMillan [273]. Within this framework, the retarded Green’s functions can be simply obtained by considering a scattering region and appropriately combining incoming and outgoing waves as
where the subscript t denotes the transpose operation and |$\Psi ^{(j)}_\mathrm{out(in)}$| represents outgoing (incoming) waves for electron- and hole-like quasiparticles; here the labels j = + and j = − denote the electron- and hole-like characters of the wavefunction. In Eqs. (90), |$\Psi ^{(j)}_{m}$| and |$\tilde{\Psi }^{(j)}_{m}$| represent solutions to the BdG equations defined by
where |$\hat{H}(x,x^{\prime })$| and |$\hat{H}^{t}(x^{\prime },x)$| are, respectively, BdG and transposed BdG Hamiltonians whose form depends on the respective problem at hand.
To fully characterize Gr(x, x′, E), it is necessary to find the coefficients αi and βi, which is done by employing the equation of motion [10,274]:
By integrating once (twice) around x = x′, this equation provides the necessary boundary conditions to find αi and βi. Once the retarded Green’s function is fully characterized, we are in a position to obtain the LDOS and pair amplitudes. It is important to point out that in Nambu space the Green’s function is a 2 × 2 matrix, which can be written as
where the diagonal represents the normal electron–electron and hole–hole Green’s functions, while the off-diagonal components correspond to the anomalous electron–hole and hole–electron Green’s functions. It is thus evident that the diagonal and off-diagonal parts of Eq. (93) enable the calculation of the LDOS and superconducting pair amplitudes, respectively. In particular, the LDOS can be calculated as
where the Tr is taken over the other degrees of freedom such as spin.
In relation to the pair amplitudes, if spin is an active degree of freedom, then Fr(x, x′, E) is a matrix and needs to be decomposed in order to know which type of superconducting pair symmetry is favored. Thus, we can define the spin components of the pair amplitudes as |$f^{r}_{\sigma \sigma ^{\prime }}(x,x^{\prime },E)=[F^{r}(x,x^{\prime },E)]_{\sigma \sigma ^{\prime }}$|, and, taking into account the Fermi–Dirac statistics, we find the following condition:
where we have used the advanced Green’s function in passing to negative energies (or frequencies) Ga(x, x′, E) = [Gr(x′, x, E)†]. The condition given by Eq. (95) is extremely relevant because it determines all the possible emergent pair symmetries in the presence of the quantum numbers determined by spins, spatial coordinates, and finite energy. Therefore, an allowed pair symmetry must be antisymmetric under the exchange of all the quantum numbers involved.
To see the possible pair symmetries in a system with active spin, considering the Nambu basis given by |$(c_{\uparrow },c_{\downarrow },c^{\dagger }_{\uparrow },c_{\downarrow }^{\dagger })$|, it is useful to write Fr(x, x′, E) as
where σj is the jth Pauli matrix in spin space and repeated indices in the second term imply summation. Thus, |$f_{0}^{r}$| and |$f_{3}^{r}$| correspond to spin-singlet and mixed spin-triplet pairs with spin projection Sz = 0, written as |$f_{0(3)}^{r}=(f_{\uparrow \downarrow }^{r} \mp f_{\downarrow \uparrow }^{r})/2$|. Similarly, |$f_{1}^{r}=(f_{\downarrow \downarrow }^{r}- f_{\uparrow \uparrow }^{r})/2$| and |$f_{2}^{r}=(f_{\downarrow \downarrow }^{r}+f_{\uparrow \uparrow }^{r})/2i$| give rise to |$if_{2}^{r}\mp f_{1}^{r}$|, which correspond to equal-spin-triplet pairs with spin projection Sz = ±1. We thus see that pair amplitudes can have spin-singlet and spin-triplet symmetries. Moreover, to comply with the antisymmetry condition imposed by Eq. (95), each of the amplitudes |$f_{i}^{r}$| can still be even or odd under the exchange of spatial coordinates and/or energy as long as the total exchange of quantum numbers makes |$f_{i}^{r}$| antisymmetric. Thus, spin-singlet pairs are allowed to have the following symmetry: even-frequency, spin-singlet, even-parity (ESE) or odd-frequency, spin-singlet odd-parity (OSO) [275–277], both pair symmetries consistent with the antisymmetry condition. For the spin-triplet pair amplitudes, they are allowed to have the following symmetries: even-frequency, spin-triplet, odd-parity (ETO) or odd-frequency, spin-triplet even-parity (OTE) [89]. Therefore, when the pair amplitude depends on the frequency, spins, and spatial coordinates of the paired electrons, there are four pair symmetries allowed by Fermi–Dirac statistics; see Refs. [86,87,92,94,248]. We stress that, for spin-polarized p-wave superconductors, the equal-spin triplet component is relevant, implying the appearance of either ETO or OTE pair symmetries, which is the focus of the following subsections.
4.2. Bulk p-wave superconductor
Here we use the recipe discussed in the previous subsection and calculate the Green’s function for a spin-polarized p-wave superconductor in the bulk. Before going further, it is important to write down the spatial dependence of the Hamiltonian for a 1D spin-polarized p-wave superconductor, which, in general, is given by
Now, we denote the electron and hole components of Ψ(x) and |$\tilde{\Psi }(x)$| as u(x) and v(x), respectively. Then, by applying the quasiclassical approximation used in Eq. (25), we obtain
where |$\Delta (\hat{k},x)=\hat{k} \Delta _{0}$| is the p-wave pair potential that satisfies |$\Delta (\hat{k},x)=-\Delta (-\hat{k},x)$|, with |$\hat{k}={\rm sgn}(k)$| and Δ(k) = −Δ(−k). By plugging Eqs. (98) into the BdG equations defined by Eq. (97), and using Eqs. (91), we obtain the incoming and outgoing wavefunctions in the bulk. For μ ≫ Δ0 and μ ≫ |E|, they read
where Γ = Δ0/(E + Ω), Ω is defined in Eq. (36), and
Having found the incoming and outgoing wavefunctions, we can now calculate the Green’s function for a bulk p-wave superconductor without any interface by using Eq. (90). Since there is no interface, the processes of both normal and Andreev reflections are absent, restricting the combination of incoming and outgoing waves in Eq. (90) to occurring only between electron-like (and hole-like) quasiparticles, i.e. α2, 3 = 0 and β2, 3 = 0. To find the rest of the coefficients (α1, 4 and β1, 4), we use the boundary conditions defined by Eqs. (92), which here acquire the following forms:
This way we find α1, 4 and β1, 4 and, after a simple manipulation considering μ ≫ Ω and k± = kF ± γ, we obtain the Green’s function in a bulk spin-polarized p-wave superconductor given by
where η = m/(iℏ2kF), γ = Ω/(ℏvF), vF = 2μ/(ℏkF) is the Fermi velocity, |$k_{F}=\sqrt{2m\mu /\hbar ^{2}}$|, and |$\hat{\tau }_{i}$| represents the ith Pauli matrix in Nambu space. It is now evident that, as discussed in Section 4.1, the retarded Green’s function obtained in Eq. (105) is a 2 × 2 matrix in Nambu space, with its diagonal and off-diagonal elements representing the normal and anomalous Green’s functions:
Thus, the LDOS obtained by using |$G_{0}^{r}$| in Eq. (94) is given by ρ(x, E) = −Im[ηE/(πΩ)], as expected in a bulk 1D spin-polarized p-wave superconductor [96]. In relation to the pair amplitude Fr, we first note that it corresponds to a spin-polarized pair potential, so its spin symmetry is triplet. Second, Fr is an odd function under the exchange of x and x′, taking zero value locally in space at x = x′, a symmetry that reflects its p-wave symmetry. Third, if we define Fa(x, x′, E), which is an advanced part of Fr(x, x′, E), Fr(x, x′, E) satisfies Fr(x, x′, −E) = Fa(x, x′, E), pointing out its even-frequency symmetry. In summary, the pair amplitude Fr exhibits an even-frequency, spin-triplet, p-wave symmetry, an ETO pair symmetry class according to the antisymmetry condition given by Eq. (95) and consistent with the symmetry of the pair potential for the spin-polarized p-wave superconductor considered here; see Eq. (97).
4.3. Semi-infinite p-wave superconductor
In this part, we inspect the Green’s function when having one interface in the 1D spin-polarized p-wave superconductor discussed in the previous section. In particular, we consider a semi-infinite 1D spin-polarized p-wave superconductor with the interface at x = 0, as schematically shown in Fig. 9.
Schematic representation of a semi-infinite 1D spin-polarized p-wave superconductor located for x ≥ 0 and with a interface at x = 0. In the topological phase, it hosts a single MZM at x = 0.
As discussed in Section 4.1, to obtain the Green’s function, we need to obtain the incoming and outgoing waves for |$\hat{H}(x,x^{\prime })$| and |$\hat{H}^{t}(x,x^{\prime })$|, which in this case are given by
where f1(x) = cos (k+x) − cos (k−x), f2(x) = Γ2sin (k+x) − sin (k−x), and f3(x) = sin (k+x) − Γ2sin (k−x). Moreover, we emphasize that Eqs. (107) satisfy the boundary conditions at the interface |$\Psi _\mathrm{in}^{\left(1 \right)}\left(x=0 \right)=\Psi _\mathrm{in}^{\left(2 \right)}\left(x=0 \right)=0$| and |$\tilde{\Psi }_\mathrm{in}^{\left(1\right)}\left(x=0 \right)=\tilde{\Psi }_\mathrm{in}^{\left(2\right)}\left(x=0 \right)=0.$| Thus, by plugging Eqs. (107) into Eqs. (90) and using the boundary condition given by Eq. (104), we obtain the retarded Green’s function [96]
Despite its apparent complexity, this Green’s function reveals interesting features due to the boundary effect. The first feature to note is that |$\hat{G}^{R}$| here has two clear parts, one involving properties of the bulk and another involving properties of the boundary. In fact, the element in square brackets in the second line of Eq. (108), and multiplied by the exponential coefficient, is exactly the bulk Green’s function found in the previous subsection and given by Eq. (105). Thus, the rest of the elements can be clearly interpreted to be an effect coming entirely from the boundary.
To further inspect the impact of the boundary, let us write down the normal and anomalous components of Eq. (108), which read
where the first line of each expression for |$G_{0}^{r}$| and Fr represents the contributions coming from the bulk, while the remaining lines show the contribution due to the boundary. It is now useful to inspect the LDOS, which, using |$G_{0}^{r}$| in Eq. (94), reads
This LDOS expression shows that, when a boundary is present as we have here, the total LDOS contains a bulk contribution (first term) and also a component that depends on space (second, third, and fourth terms in square brackets). For subgap energies, Ω is imaginary as seen in Eqs. (36), which implies that the exponential term in front of the square brackets exhibits an exponential and oscillatory decay for x > 0 that, interestingly, strongly depends on E. Of particular relevance is the behavior close to E = 0, which can clearly induce a divergent profile around the edge component of the LDOS; see the last term in Eq. (110). This property is a manifestation of the emergence of an MZM at the boundary. We note that the LDOS for a semi-infinite spin-singlet s-wave superconductor is given by [96]
which, although partially similar to Eq. (110), does not contain a contribution that produces a zero-energy peak (ZEP) at E = 0. Therefore, a ZEP in the LDOS can be seen to be a manifestation of MZMs at the edge of spin-polarized p-wave superconductors.
When it comes to the pair amplitude Fr in Eq. (109), we also note two contributions arising from the bulk (first line) and boundary (second line). In Section 4.1, we identified that the symmetry of the bulk pair amplitude is ETO and therefore vanishes at x = x′. Interestingly, in Eqs. (109) we see that Fr retains contributions even when x = x′, which are given by
an expression that has strong similarities with the last term in the LDOS of Eq. (110). This pair amplitude is local in space, thus having even parity (s-wave), and decays in an exponentially and oscillatory fashion as x increases. Moreover, it is spin-polarized, so it has spin-triplet symmetry. When it comes to the energy dependence (or frequency), it is evident that it is an odd function, with a divergent profile at E = 0, revealing its odd-frequency symmetry; as imposed by Eq. (95), when E ↔ −E we compare with the advanced counterpart of Fr(x, x, E) in Eq. (111). Thus, the pair amplitude in Eq. (111) has OTE symmetry, which is in agreement with the antisymmetry condition given by Eq. (95) discussed in Section 4.1. For completeness, we also write the pair amplitude in Eq. (111) in Matsubara frequencies by performing an analytic continuation as E → iωn where ωn = 2πκBT(n + 1/2), which reads
where |$\gamma _{n}=\sqrt{\omega _{n}^{2} + \Delta _{0}^{2}}/\hbar v_{F}$|. We thus see that Eq. (112) is an odd function under ωn ↔ −ωn. The OTE symmetry of Fr(x, x, E), with intriguing dependences, reveals the nature of the emergent pair correlations in the presence of MZMs.
4.3.1. Quasiclassical Green’s functions
Further understanding of the relationship between MZMs and OTE pair correlations is obtained by analyzing the quasiclassical Green’s function [278]. For this purpose, it is useful to write the retarded Green’s function found in Eq. (108) as
where |$\hat{\mathcal {G}}^{\alpha \beta }(x,x^{\prime }) \equiv \hat{\mathcal {G}}^{\alpha \beta }(x,x^{\prime },E)$|, with α = ±, β = ±, expresses the envelope of the Gor’kov Green’s function given by
where |$\hat{\tau }_{i}$| is the ith Pauli matrix in Nambu space.
Now, we use |$\hat{\mathcal {G}}^{++}(x \pm 0,x,E)$| and |$\hat{\mathcal {G}}^{--}(x \pm 0,x,E)$| and define a new Green’s function, also referred to as a quasiclassical Green’s function [279], as
where |$\hat{\gamma }_{3}^{\alpha \alpha }$| is obtained from |$\hat{\gamma }_{3}^{++}=\hat{\tau }_{0}$| and |$\hat{\gamma }_{3}^{--}=-\hat{\tau }_{0}$|. Then, we obtain expressions for |$\hat{g}^{++}(x, E)$| and |$\hat{g}^{--}(x, E)$|, given by
where
Before further showing the properties of Eq. (117), we note that the quasiclassical Green’s functions |$\hat{g}^{\alpha \alpha }$| satisfy the normalization condition given by
which is specific to the quasiclassical Green’s functions and can be proved by using Eqs. (116) and (117); see Ref. [279]. Moreover, we note that |$\hat{g}^{++}(x,E)$| and |$\hat{g}^{--}(x,E)$| satisfy the following differential equation given by
where |$\hat{\Delta }_{tr}=i\Delta _{\alpha }\hat{\tau }_{2}$|, with α = ±, Δ± = ±Δ0. Equation (119) is also known as the Eilenberger equation; see Ref. [280] for more details.
Now, we are in the position to discuss how properties of Eqs. (117) can be seen as the quasiclassical counterparts of the expressions found in Eqs. (109) but with a small variation for the normal component. In fact, while g±(x) expresses the normal Green’s function with quasiparticle excitation, since |$\hat{g}^{\alpha \alpha }(x,E)$| is defined multiplied by i in Eq. (115), the real part of g±(x) already represents the normalized LDOS, which by construction is normalized by its value in the normal state. Furthermore, f1 ±(x) and f2 ±(x) represent the pair amplitudes that, as can be seen in Eqs. (117), correspond to odd- and even-frequency symmetries, respectively [86,87]. As a note, we mention that, within the quasiclassical approach, the spatial coordinate x is now the center of mass coordinate. Another feature to note is that g±(x) and f2 ±(x) exhibit contributions from the bulk and boundary: the bulk term is independent of space, while the boundary term depends on space, in a similar manner to that obtained for Eqs. (109).
In Figs. 10, 11, and 12, we show the spatial dependence of the quasiclassical Green’s functions g(x) = g+(x), f2(x) = f2 +(x), and f1(x) = f1 +(x). For numerical purposes, we replace E by E + iδ with an infinitesimal imaginary number δ = 0.001Δ0. When inspecting g(x) and f1, 2(x) at energies above the gap E > Δ0, they develop an oscillatory behavior, as predicted by Eqs. (117). At x = 0, f2(x) = 0, while g(x) and f1(x) exhibit finite values. For energies within the gap, |E| < Δ0, the factor γ becomes imaginary according to Eq. (36), inducing an exponential decay from the interface in the boundary contributions, with a decay length ∼2ξ; see Fig. 11. For the chosen energy, the imaginary parts g(x) and the real part of f1, 2(x) become dominant, taking values of the same order. It is interesting to note that the imaginary part of g(x) decays towards the bulk of the system, a behavior that is also exhibited by the OTE pair amplitude f2(x). In contrast, the real part of the ETO amplitude increases as x takes values deeper inside the superconductor; see Fig. 11. This is of course an expected behavior because the pair potential of the bulk superconductor exhibits ETO pair symmetry.
Spatial dependence of the quasiclassical Green’s functions at E/Δ0 = 1.2, where the different panels correspond to A: g(x), B: f2(x), and C: f1(x). In each panel, the solid (a) and dashed (b) curves show the Re and Im parts, respectively. Parameters: δ = 0.001Δ0, ξ = ℏvF/Δ0.
Spatial dependence of the quasiclassical Green’s functions at E/Δ0 = 0.5, where the different panels correspond to A: g(x), B: f2(x), and C: f1(x). In each panel, the solid (a) and dashed (b) curves show the Re and Im parts, respectively. Parameters: δ = 0.001Δ0, ξ = ℏvF/Δ0.
Spatial dependence of the quasiclassical Green’s functions at E/Δ0 = 0, where the different panels correspond to A: g(x), B: f2(x), and C: f1(x). In each panel, the solid (a) and dashed (b) curves show the Re and Im parts, respectively. Parameters: δ = 0.001Δ0, ξ = ℏvF/Δ0.
To close this part, we analyze the spatial dependence of the quasiclassical Green’s functions at E = 0, which is presented in Fig. 12. As we observe, the spatial dependence of the real part of g(x) and the imaginary part of f1(x) coincide, developing a maximum at x = 0 with huge values, followed by an exponential decay towards the bulk of the system. This peak at the boundary is caused by the presence of an MZM, since the spin-polarized p-wave superconductor that we consider here is in the topological phase. In consequence, this further demonstrates that the emergence of an MZM at the boundary is always accompanied by the generation of OTE pairing and that there is a one-to-one correspondence between them [42,116]. This view is further supported by Figs. 13 and 14, which show the energy dependence of the LDOS Re[g(x)] and the OTE pair amplitude Im[f1(x)] at x = 0. Since the real part of g(x), namely, Re[g(x)], represents the LDOS normalized by its value in the normal state, there is a direct correspondence between the ZESABS and OTE pairing at very low energies. A huge peak forms at E = 0 in both the LDOS and the OTE pair amplitude due to the emergence of an MZM at the boundary, which demonstrates the strong relationship between LDOS and OTE at E = 0. It is worth noting that for finite subgap energies (|E| < Δ0) the LDOS and OTE pairing have the same dependence, while only the LDOS captures the gap edges at |E| = Δ0. Nevertheless, the intriguing signature of both the OTE pairing and LDOS is their ZEP due to an MZM at the boundary. As we enter into the superconductor, the height and width of the ZEP is suppressed as shown in Fig. 14 at x = ξ. This behavior is indeed consistent with the Majorana origin of the ZEP, thus showing that both the LDOS and the OTE pair amplitude at the boundary of a topological superconductor can be used to detect MZMs and their exotic pair symmetry.
![Energy dependence of LDOS and OTE pair amplitude at x = 0, where A: Re[g(x)] and B: Im[f1(x)]. δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig13.jpeg?Expires=1747856791&Signature=ykYvgMOf1vGyuGmOLZuPnzQ3dCnn1Jn-u0XazJZSCdoZwrIftelggQCQz5rlkQYLNQhy4~cDVo8DMJrORdyfUDsMfhk9ZwmBp2CccvGsio5RNBZlQF66jB~G05FMJ4U33OfBzfb1Gfwt4eJCL~8Q7c4rOYrzUzJXvYKn1qoJXtMDbMKTTH4LaLeaPBOXBkblyRDTKVVgH1SwoqzUeeYO3~~OYZS5b2IpTFs9zRGkuedrB8UqrN9-rYeDhZTWIKrj6qW4GmA5XotiheWf3F1XdfuXPXkyVc3lTypCUx3aC8EEY~5Mm3QdkVXchED1ILRR4jDhh5FcostQUqKgurXqRA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Energy dependence of LDOS and OTE pair amplitude at x = 0, where A: Re[g(x)] and B: Im[f1(x)]. δ = 0.001Δ0.
![Energy dependence of LDOS and OTE pair amplitude at x = ξ, where A: Re[g(x)] and B: Im[f1(x)]. δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig14.jpeg?Expires=1747856791&Signature=SKsHJivSCvTuySJxZyYS4hFobTHcwXAZBWdT7uvJur9TjSpPonvdekUAIfM8wvx7ZE1Es58-zz45ODpg4OI75vmajoXXW4bNTL7l0dXDyzq4Dy4fZKdk5zqqyKuy7eCXSG-18cyWHbWxdHhQhhpbS0jHD13Gs2VijGPkPM~j~ula~cLn6vPZVdjeJghpG5bC~k28XJxSga-mIdzXYat1AI0IBrn774cdzrFLzV1ZcEOvgiZPpFl4KHPkeoM81a4RZFEK7idMgS4m6m159lXM4fpJyjzRXgzckLKcq36a-gie0YbSJZYiM6bfato6b0ZoBgtuugmRcMKLJxMTPAnG4Q__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Energy dependence of LDOS and OTE pair amplitude at x = ξ, where A: Re[g(x)] and B: Im[f1(x)]. δ = 0.001Δ0.
4.4. Finite length p-wave superconductor
In this part we explore the Green’s function properties of a spin-polarized p-wave superconductor with finite length L, as schematically shown in Fig. 15. Since MZMs appear at the boundary in the topological phase, we expect that an MZM emerges at each end of the system as depicted by red filled circles in Fig. 15.
Schematic of the finite length Kitaev chain model for 0 ≤ x ≤ L. Two Majorana bound states are located at x = 0 and x = L.
Before constructing the Green’s function for a finite length 1D p-wave superconductor, it is instructive to first obtain the energy of the bound state, following the same spirit of Section 2.4. In this regard, we recall that the total wavefunction of the bound state of the 1D p-wave superconductor described by Eq. (97) for |E| < Δ0 reads
which is composed of two electron-like (proportional to A and B) and two hole-like (proportional to C and D) wavefunctions at momenta ±k± given by Eq. (103). Now, we use hard-wall boundary conditions at both x = 0 and x = L and then write down a system of equations for the total wavefunction in Eq. (120), namely, Ψ(x = 0) = 0 and Ψ(L = 0) = 0. We thus obtain the following system of equations for the coefficients:
A nontrivial solution of a bound state is then obtained by postulating that these four equations should be satisfied for nonzero values of A, B, C, and D, in a similar way as we have done in Section 2.4 when deriving the bound state energy at one boundary. Thus, in the present case, we obtain the bound state solution given by
where δs and δf are defined by
The bound state expression given by Eq. (122) can be further simplified for μ ≫ Δ0, where we can approximate cos δ ∼ cosh(2γBL) and cos δ1 ∼ cos(2kFL), with |$\gamma _{B}={\sqrt{\Delta _{0}^{2}-E_{b}^{2}}}/{\hbar v_{F}}$|. Then, the bound state energy from Eq. (122) reduces to
At this point, we are in a position to draw some conclusions about the bound state energy of a spin-polarized p-wave superconductor of finite length L. The first feature to notice is that there are two solutions, a positive and a negative, which, under general conditions of finite but short L, exhibit an oscillatory behavior with sizeable amplitudes and crossings at zero energy when kFL = nπ with |$n\in \mathbb {Z}$| like in a Fabry–Pérot interferometer. Furthermore, it can be shown that the amplitude of such oscillations increases as the chemical potential increases. Interestingly, for very large L, the two bound state energies acquire vanishing values such that Eb = 0. These zero-energy bound states represent the two MZMs located at the ends x = 0 and x = L of the 1D spin-polarized p-wave superconductor.
Having understood that a finite topological superconductor hosts two MZMs, one at each end, whose energies depend on the system length, now we turn our attention to their impact on the Green’s function. As discussed in Section 4.1, in order to construct the Green’s function, we must first obtain the outgoing and incoming wavefunctions associated with the BdG equations of |$\hat{H}$| and |$\hat{H}^{t}$|; see Eqs. (91). The outgoing and incoming wavefunctions associated with |$\hat{H}$| are given by
which satisfy the boundary condition given by
while the coefficients a and b read
For |$\hat{H}^{t}$|, the outgoing and incoming wavefunctions are given by
which satisfy
Then, by using the incoming and outgoing states given by Eqs. (125) and (127) in Eqs. (90) and (104), we find the retarded Green’s function, which, after some manipulations, reads
where
and |$X_{1}=\hat{\tau }_{0}$|, |$X_{2}=\hat{\tau }_{3}$|, |$Y_{1}=\hat{\tau }_{1}$| and |$Y_{2}=i\hat{\tau }_{2}$|. In the above κ±, γ, and Γ are given by
with Ω given by Eq. (36). At this moment, it is worth highlighting some properties of the retarded Green’s function given by Eq. (128). Even though it has an evident complicated structure, it is possible to note that it depends on the length of the superconductor L in an oscillatory fashion via sines and cosines in the square brackets of the first and second lines. Moreover, the Green’s function depends on the spatial coordinates x and x′ via the coefficients Ci. Lastly, there is an overall energy dependence, for instance, via Γ = Δ0/(E + Ω) and Ω defined in Eq. (36). It is a 2 × 2 matrix as expected, with its diagonals and off-diagonals representing the normal and anomalous Green’s functions; see Eqs. (93). For visualization purposes, it is useful to decompose the total Green’s function into its normal and anomalous parts as
where |$\mathcal {G}(x,x^{\prime },E)$| corresponds to the normal Green’s function, and |$\mathcal {F}_{1}(x,x^{\prime },E)$| and |$\mathcal {F}_{2}(x,x^{\prime },E)$| to the OTE and ETO pair amplitudes, respectively. We now explore the correspondence between MZMs and the OTE pairing and, for this reason, we look at |$\mathcal {G}(x,x,E)$| and |$\mathcal {F}_{1}(x,x,E)$| defined in Eq. (135). They are given by
By plugging in the expressions for Γ and κ± from Eqs. (133) and (134), we obtain more compact expressions given by
In the following, we take these two previous expressions and analyze the corresponding LDOS and pair amplitude from the normal and anomalous parts, respectively. Since the spin symmetry is spin-triplet and we are evaluating locally in space, the pair amplitude has an OTE symmetry. In Figs. 16 and 17, we plot |${\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$|, |${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$| as a function of E at a fixed position close to the edge. In this case, |$\rho (E)=-{\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$| represents the LDOS normalized by its value in the normal bulk state value. For LkF = 1000π, the LDOS |${\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$| has a sharp peak at zero energy and, interestingly, a similar line shape appears for the real part of the OTE pair amplitude |${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|, as shown in Figs. 16(A), (B). The imaginary part |${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$| is also enhanced around E = 0 but it is an odd function of E in contrast to the |${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|; see Fig. 16(C). Away from the resonant condition given by LkF = 1000.5π, the LDOS and the OTE pairing develop a small splitting, which is reflected in the formation of two peaks around zero energy; see Figs. 17(A), (B). This zero-energy splitting is a consequence of the finite spatial overlap of the two MZMs emerging at either side of the system, an effect that can be understood by a simple inspection of Eq. (124). At the position of the bound states around zero energy, the amplitudes of both |${\rm Real}[\mathcal {F}_{1}(x_{b},x_{b},E)]$| and |${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$| are enhanced, as seen in Figs. 17(B) and (C).
![Energy dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at xb = L/1113 and LkF = 1000π: A: $\rho (E)=-{\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$, B: ${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$, and C: ${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig16.jpeg?Expires=1747856791&Signature=KJK2fgmk~g97K5ok1kjgqcJKegc5mRqA~LxFR-Y5aRpmFR3DMkc1I8FNd9lSZY2mNngaT-~NG78Rh7Qp3xhDYLmlxGCGkDGZb5~ZGW50dux9okIIlr0CMak7gr0wxBKlKzttYriL8Z8DvwiiNG7b-Rboi7HORVNRuhhPWs6YQ5uUhCln1rqzrnzpMpsUpYZWUD7r8iIF24RtLV6QLe2JWCQP9oJs3HV5P4ImrwlM2gFT7Zwpxbf-m0TF7o-wcG6doA3GiqYIrr7SDyY4R6FCYblO9yFhJgOzu120I-uc77wPya~wJsySsCxNaMh5oOWRaiG-NEBBNb-VkAfa4koeGA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Energy dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at xb = L/1113 and LkF = 1000π: A: |$\rho (E)=-{\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$|, B: |${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|, and C: |${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.
![Energy dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at xb = L/1113 and LkF = 1000.5π: A: $\rho (E)=-{\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$, B: ${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$, and C: ${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig17.jpeg?Expires=1747856791&Signature=2YwfSa9usfGyLmfPojgCDuQQpYdk4mtRFvCT4p9O4VvDPUcD9iYL7Ce2OlFYSyqRXV2hfUvuqv0va7ED4lwLrDDsC-N3Hlj8iSwjnN8XsWGYdZ9bEqWducBus903xfIbqKTS1Ja9SzUUA08Zh03oaVRGagSXpjxmvqUCHoU~XJ1FlcfuqxatizOLUL56E6-FgblNzrOQwr2L1UfwVJzSMiUxQichIJ88jT4isEUgU9A-NyOcr~yz3y6RvFB6vJ9G4wCB77E-EYXMuozMekbxeMMVpgVfYGXtf1jSfsrdUj3-IrS-WYhzIcSxke3fcBMcV7y8q0J42cZCLVvph48haQ__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Energy dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at xb = L/1113 and LkF = 1000.5π: A: |$\rho (E)=-{\rm Im}[\mathcal {G}(x_{b},x_{b},E)]$|, B: |${\rm Re}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|, and C: |${\rm Im}[\mathcal {F}_{1}(x_{b},x_{b},E)]$|. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.
Having seen the relation between the local OTE pairing and LDOS as a function of energy, now we explore their spatial dependence. For this purpose, in Figs. 18–20, we plot ρ(E), |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| as functions of x for distinct values of the system length. The first observation is that all these quantities have a rapidly oscillatory behavior with a period proportional to the inverse of kF and an exponential decay from the edges with the decay length given by ℏvF/Δ0. The exponential decay reflects the localization of MZMs at the edges of the system, where only very long systems host true MZMs, as indeed seen in Figs. 18, 19, and 20. Another feature is that the spatial dependence of ρ(E), |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| is very similar in all cases, although the Re and Im parts of the pair amplitudes exhibit distinct signs at the edges. Furthermore, we note that the fact that the pair amplitudes exponentially decay from the edges reveals that this type of pairing is a unique property tied to the edge where MZMs form. As a result of this, both |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$| and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| become zero at x = L/2, a common feature of the OTE pair amplitude that is independent of L. As seen from Fig. 18, |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$| and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| become positive (negative) with L/2 > x > 0 (L/2 < x < L) for kFL = 1000. This property is a general profile for kFL = nπ where sin kFL = 0 is satisfied as seen from Eq. (139).
![Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 1000π: A: $\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$, B: ${\rm Re}[\mathcal {F}_{1}(x,x,E)]$, and C: ${\rm Im}[\mathcal {F}_{1}(x,x,E)]$. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig18.jpeg?Expires=1747856791&Signature=X85gC~pGBCbJ890Ru2~UjhQWg341LqwqHOC~eRX84XC3pm0sYyCVlW5xWz0Cy1BZWpas2JUoSD6TrGw7ZAvbHIOq-ZQYl4WiHnPPEP89i093WhNFzVQ09vXbKZy6grVw5k~KRHPv43ReFfj9vKNwvt1HeArnQ7Wc8AL3xk0lDkpnJ9~na6zC-HfR52b1zExpwdvKrelJ1vhuc5NzYZTQ8rAopVapEeRIcutOxCLqLZHa-OeF~n2IkwAF~HoHiN-eg7F4lalqLa4b5EQtytR96Md16QDKSQkTaTQbtSYGTIFr4-PcnMB7Gdfo6jpJ~Moos5xrMeejQsTACikewiUpyg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 1000π: A: |$\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$|, B: |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and C: |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$|. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.
![Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 1000.5π: A: $\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$, B: ${\rm Re}[\mathcal {F}_{1}(x,x,E)]$, and C: ${\rm Im}[\mathcal {F}_{1}(x,x,E)]$. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig19.jpeg?Expires=1747856791&Signature=K5ERRZh43NufbTYT4iUm18Ju61o24YkOOf2bUwET7TUfX6Gvi9MptT62HRo9TFs4lremQGcKBA0a-EkVZ1XbdixIJz5X07ZmO2M0OC6USZ0PC-26fuG35thCiM-EnJivaGiJH~XG2bSITnVvZJyZwbSvp5~aaMzAmqgQenMeh7AdzviUsCt4mlanzNy9ozg7GpaHn7l7EEnXyT3e1eM8fHj5gKbUv7LreplRH5TU6KlKZimwKQ-97YgaSfHx15OrJX0blPo3hK0XQlA34mOvpXBYr9Rjmj2REACTGb91J59-PAIL3Tl5ESHAuDxpxwdZtkeZkukXpGvccqIDj9oMTg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 1000.5π: A: |$\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$|, B: |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and C: |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$|. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.
![Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 3000.5π: A: $\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$, B: ${\rm Re}[\mathcal {F}_{1}(x,x,E)]$, and C: ${\rm Im}[\mathcal {F}_{1}(x,x,E)]$. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig20.jpeg?Expires=1747856791&Signature=BoPw9EpzGPb0ojdKAfFK2BR-rSCBPS1gLAos832PukzwVEzFITXb~qGspmeQf0h0PsTsl7QTMlfBqnTdqMAHrRBRCbEXGhPjPiQIY8c5U4orYCxT4CZ5yGYGEA9Wof4QoP5OaCRMcmBe8K49OyG77RO9HNW4GMdZXPZIFUq0s0tG3NStO3KHEhSMSRai8A1oROPlhefYtP6Po4~heHs2JOxPlVOyKpqOTPZKUeiUaoqWF0qq5vMX81sLDc7-i3TxCtrcwjMM7QETWlMuVUG5lAFD3uVYiqlPnPq7~iUSewO4WdfZRPeuvASbmoVWIuerAC6tWj7D9RcIY8BHT4rxqw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Spatial dependence of the LDOS and the local OTE pair amplitude in a finite 1D spin-polarized p-wave superconductor at E = 0.015Δ0 and LkF = 3000.5π: A: |$\rho (E)=-{\rm Im}[\mathcal {G}(x,x,E)]$|, B: |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and C: |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$|. Parameters: Δ0 = 0.002μ and δ = 0.001Δ0.
On the other hand, for sin (kFL) = 1000.5π, as shown in Fig. 19, the magnitudes of ρ(E), |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| are suppressed as compared to the case with sin (kFL) = 1000π (Fig. 18). This is due to the destructive interference of MZMs localized at the left and right edges of the system. For kFL = 3000.5π, ρ(E), |${\rm Re}[\mathcal {F}_{1}(x,x,E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x,x,E)]$| are strongly localized at x = 0 and L, and these values are significantly suppressed around x = L/2 (Fig. 20).
We next explore the Green’s function components by fixing one spatial coordinate close to the left edge and evolving with the other coordinate as we move towards the bulk of the system. This is particularly useful for obtaining information about nonlocal correlations, which can provide further understanding of MZMs because they exhibit intrinsic spatial nonlocality. In Figs. 21, 22, and 23 we plot the spatial dependence of |${\rm Re}[\mathcal {G}(x_{b}, x^{\prime },E)]$|, |${\rm Im}[\mathcal {F}_{2}(x_{b}, x^{\prime },E)]$|, and |${\rm Im}[\mathcal {F}_{1}(x_{b}, x^{\prime },E)]$| at xb = L/1113 at E = 0.015Δ0 + iδ, with δ = 0.001Δ0. Thus, changing g the values of x′ provides information about the nonlocal properties of these Green’s functions. Before describing the features of the Green’s functions, we note that, at this energy, the real part of |$\mathcal {G}(x_{b}, x^{\prime },E)$| is larger than the imaginary term. On the other hand, for |$\mathcal {F}_{1}(x_{b}, x^{\prime },E)$| and |$\mathcal {F}_{2}(x_{b}, x^{\prime },E)$| the situation is that their imaginary components are larger than their real counterparts. In Fig. 21, we choose kFL = 1000π where the wavefunction of an electron in the normal state vanishes at both x = 0 and x = L. Since we set Δ0/μ = 500, we get L/ξ = π with ξ = ℏvF/Δ0. The first feature that we identify is that |${\rm Re}[\mathcal {G}(x_{b}, x^{\prime },E)]$| and |${\rm Im}[\mathcal {F}_{1}(x_{b}, x^{\prime },E)]$| acquire equally large values, which exceed those of |${\rm Im}[\mathcal {F}_{2}(x_{b}, x^{\prime },E)]$|. Moreover, we observe that |${\rm Re}[\mathcal {G}(x_{b}, x^{\prime },E)]$|, |${\rm Im}[\mathcal {F}_{1}(x_{b}, x^{\prime },E)]$|, and |${\rm Im}[\mathcal {F}_{2}(x_{b}, x^{\prime },E)]$| have an oscillatory behavior, with gradually decaying envelope functions as we move further from the left interface at x = 0. The order of the period of the rapid oscillations is the inverse of kF and their amplitude decays with a decay length determined by ξ. When we move to LkF = 1000.5π, which is slightly different from the condition for kFL discussed in the previous paragraph, the situation changes dramatically, as clearly seen in Fig. 22. We first note that the magnitudes of |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$| and |${\rm Im}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$| are suppressed around x′ ∼ 0 as compared to the case for LkF = 1000 in Fig. 21. Interestingly, the amplitudes of the rapid oscillations are enhanced with the increase of x′ for |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$| and |${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$|. On the other hand, for the OTE pair amplitude |${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$| we see that the amplitude of the rapid oscillations is suppressed. The sensitivity of the L dependence of the Green’s functions shown in Figs. 21 and 22 is a consequence of the hybridization of the wavefunctions of MZMs localized at the left and right edges. We also note that the suppression of |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$| and |${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$| can be seen at other values of kFL, as we show in Fig. 23 for kFL = 3000.5π, which is again a result of the reduction of the spatial overlap between MZMs. The enhancement of the amplitudes of the rapid oscillations shown in Fig. 22 is specific to the 1D p-wave model and never appears in the corresponding s-wave model.
![Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 1000π. Distinct panels correspond to: A: ${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$, B: ${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$, and C: ${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig21.jpeg?Expires=1747856791&Signature=FRcn4dhlq5GDlrRcu1ro9Cmh--hHb4NGttztOZk8qrZuz93QsxDf0uxxVJKzg910ZWo3Ja1BEeobg6qd9PqiEDDjzSLKFhWi26E9eaVhDRd1Knz2MV5AuYDGiDEZV52HnPc0cbR-bGIYXa8eJJeVyKq9HtZyBkfH-mvVgYdIXsNp~-3AhsYqZ~x0L5GofIr5i4lVZC18AcmjkJ6XGJg2mWdTWzOqcnQzezcdCEfI7ZKlwhCpEWbFTvl4afLgAmo0cAHJGBUL1uihQCwK37rOeWd9Bd~op6H8na3SXTe6Tic6N7LmnM6vrZUCkoQPPtW5zvXCiC2lz3A2aCpPUDBSQA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 1000π. Distinct panels correspond to: A: |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$|, B: |${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$|, and C: |${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$|. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.
![Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 1000.5π. Distinct panels correspond to: A: ${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$, B: ${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$, and C: ${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig22.jpeg?Expires=1747856791&Signature=FJ36iUff4PXK4aNVhqga7zOnFMGSQX0SBZcyvXuSfBLPfwSeSOSSq~U76ep8wC4ehgbL4o0px2UuhIG1r4N~1XDFpnddK0D7oPvZmCA2XCAhY39c5rfs0z~LNCxkZb0wkJDskzkFIQRER3LlbKZz95IrXeY3m5WwvBUbE6k6UN1X~~1yBVPJ~oZJGjowMJUbj5jfSed9KG8TuQlXckoxS4hPsJGrhwaEsT-fg2X3RIZiVjZ0Tb9Q0cIP8q-qX0341N~texDBuOqHDgKvYf~jb1RkwEnWNXvVOc4lkNHLsmcCvoSSsTrcXuEo7c6bci3GLyFWe3ox7NZ9ZA-DiqEu6w__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 1000.5π. Distinct panels correspond to: A: |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$|, B: |${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$|, and C: |${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$|. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.
![Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 3000.5π. Distinct panels correspond to: A: ${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$, B: ${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$, and C: ${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig23.jpeg?Expires=1747856791&Signature=QfHMnVtyTCk7iO81f1tTgV9pRo1NlQLxYvzkyUqAOkxhS0nbD2t-0arUiegE~kEVPb7qSed1hNJpOKI3w1NdbEja3SjnhzrXikLXM8z55f5IXiRRdRUAXkqRywoWako1vJs0A7WUkBFqi~~9zSNZvxzm8k6D4jQjGVOc4CL-472TFU0RhTNGEtrMTkdlEO33CFUJscW04w3QZAHK28ryfSiyJ~VgjQEhbtnU30Vd8ox4TKaT2T8Up2xoqm4YFvB40y-SPxDS334cONLxQiLN4rJQIy9Kx-tOcGN3cIB-bxKwcj4l0ogjummmoO8eG1rlPb2fV061eb0DldusK06qdw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Normal and anomalous Green’s functions as functions of x′ in a finite 1D spin-polarized p-wave superconductor at LkF = 3000.5π. Distinct panels correspond to: A: |${\rm Re}[\mathcal {G}(x_{b},x^{\prime },E)]$|, B: |${\rm Im}[\mathcal {F}_{2}(x_{b},x^{\prime },E)]$|, and C: |${\rm Re}[\mathcal {F}_{1}(x_{b},x^{\prime },E)]$|. Parameters: xb = L/1113, E/Δ0 = 0.015Δ0, Δ0 = 0.002μ, and δ = 0.001Δ0.
4.4.1. At resonance
By focusing on the regime without splitting in the LDOS, namely, at resonance, several simplifications are possible. In particular, this regime implies that we consider
a condition that reflects having a zero-energy bound state Eb = 0 and also reveals the formation of a resonant state. In this case, the following relations are also satisfied:
Then, the retarded Green’s function Gr(x, x′, E) is given by a compact expression that reads
Thus, for L → ∞, we can see that Eq. (142) reduces to the expression for a semi-infinite p-wave superconductor presented in Eq. (108). Furthermore, it is important to remark that, at special lengths L satisfying Eq. (140), the Green’s function for a finite p-wave superconductor can be expressed by the semi-infinite Green’s function defined for x > 0 and that for x < L.
4.5. Summary
We have presented a detailed discussion about the properties of the retarded Green’s functions in 1D spin-polarized p-wave superconductors. In particular, we have focused on bulk, semi-infinite, and finite length systems, and demonstrated that their respective retarded Green’s functions reveal interesting functionalities due to the formation of MZMs. In the case of the bulk system, we have shown that the symmetry of the superconducting pair correlations exhibit the expected even-frequency spin-triplet p-wave symmetry that corresponds to the ETO pair symmetry class, which reflects the nature of the pair potential. Interestingly, for a semi-infinite system, we have found that the LDOS and the OTE pairing at the boundary exhibit a large peak at zero energy due to the emergence of an MZM, exhibiting an exponential and oscillatory decay, and can be used as detection probes. We have also highlighted that having a largely dominant OTE pairing means that the superconducting pair symmetries in the presence of MZMs are distinct from those of the bulk, which reveals yet another signature that can be further exploited for detecting the still elusive Majorana quasiparticles in condensed matter systems. In the case of a finite system, we have further shown that having a finite length produces MZMs decaying from both ends, which can then spatially overlap and produce finite oscillations around zero energy in both the OTE and LDOS. Even though very large systems are desirable for achieving MZMs and for their potential applications, having finite systems, although unavoidable, seems to be the key for probing the spatial nonlocality of MZMs and thus for distinguishing them from other topologically trivial zero-energy states.
5. Detecting Majorana zero modes via tunneling spectroscopy
This part is devoted to an analysis of how MZMs and generic ZESABSs manifest when electrons tunnel across a normal–superconductor (NS) junction, where N represents a normal metal lead while S is a superconductor with a given pair potential. In particular, we will focus on conductance, which is perhaps the simplest transport observable and represents the central concept in tunneling spectroscopy [48,69,281–284]. For completeness, we first address tunneling spectroscopy when S represents a conventional superconductor with a spin-singlet s-wave pair potential. Then, we inspect tunneling spectroscopy for superconductors with spin-singlet d-wave and spin-triplet p-wave pair potentials. Before we go further, we note that to discuss tunneling effects, we must take into account elementary scattering processes occurring at the NS interface due to the incidence of quasiparticles. To treat these scattering processes, we will employ the BdG equations given by Eq. (12), where superconducting systems with a spatially nonuniform pair potential, as junctions, can be investigated.
5.1. Tunneling spectroscopy in conventional superconductors
Conventional superconductors exhibit a pair potential that is spin-singlet and isotropic in space (s-wave). Thus, considering the 1D case for simplicity but without loss of generality, the BdG equations given by Eq. (12) can be written as
where u(v) are the electron (hole) components of the wavefunction, μ is the chemical potential, U(x) represents a one-body potential, and Δ(x) is the superconducting spin-singlet s-wave pair potential. Since we are interested in inspecting tunneling across an NS junction, the pair potential is chosen to be
which shows that the NS interface is located at x = 0 and that superconductivity in N is absent. Moreover, the nature of the interface at x = 0 is determined by the profile of the one-body potential, which we here consider to be just a delta barrier, namely, U(x) = Hbδ(x), with Hb the strength of the barrier.
We are now in a position to discuss transport across the NS junction, where multiple scattering processes occur and are determined by incident particles at the interface [285]. For instance, an incident electron from N can be reflected back to N as an electron or hole, processes known as normal and Andreev reflections. The Andreev reflection is the most fundamental transport phenomenon in superconducting junctions [286,287] and the key concept to understand the superconducting proximity effect [288,289]; see also Ref. [290]. Moreover, the incident electron can also be transmitted to S into electron-like and hole-like quasiparticles. These processes are schematically shown in Fig. 24; similar processes happen for incident holes from N or quasielectrons (quasiholes) from S. Thus, for an incident electron from N, the corresponding wavefunction is given by
where the coefficients a(E), b(E), c(E), and d(E) represent the amplitudes of Andreev reflection, normal reflection, electron transmission, and hole transmission, respectively. These coefficients are found from the boundary conditions by integrating the BdG equations around x = 0, which give
where |$k^{\pm }=\sqrt{2m(\mu \pm \Omega )/\hbar ^{2}}$|, |$p^{\pm }=\sqrt{2m(\mu \pm E)/\hbar ^{2}}$|, |$u_{0}=\sqrt{(1 + \Omega /E)/2}$|, |$v_{0}=\sqrt{(1 - \Omega /E)/2}$|, and
with δ > 0 being an infinitesimal positive number. For simplicity and also for pedagogical purposes, we now focus on the Andreev approximation, where the chemical potential is the largest energy scale, namely, μ ≫ |E| and μ ≫ |Ω|. This approximation then implies that p+ ∼ p− ∼ k+ ∼ k− ∼ kF, which then allow us to obtain simple and compact expressions for the coefficients of the wavefunction given by Eq. (145), which are given by [285]
where Γ = v0/u0 and Z = mHb/ℏ2kF. We can now obtain the probabilities of the scattering processes, which are defined by
where A, B, C, and D characterize the probabilities of Andreev reflection, of normal reflection, of transmission as an electron-like quasiparticle, and of transmission as a hole-like quasiparticle (or crossed Andreev reflection). The conservation of current then implies that these probabilities satisfy the following condition
which also reflects the unitarity of the so-called scattering matrix [291]. We can now write the scattering probabilities as
where we have used the transmissivity in the normal state σN = 1/(1 + Z2) with Z characterizing the strength of the barrier at the NS interface. To visualize the energy dependence of the scattering probabilities given by Eqs. (151)–(154), in Fig. 25 we plot them as a function of E at Z = 0 and Z = 2. For Z = 0, we have σN = 1, which represents full transmissivity and the interface is fully transparent. In this full transparent regime, the normal reflection and hole transmission probabilities vanish for all energies, namely, B = D = 0, as seen by the dotted black (b) and dotted red (d) curves in Fig. 25(A). We also see that perfect Andreev reflection occurs for subgap energies, where A = 1 and C = 0 for |E| ≤ Δ0. In the case of low transparency, Z ≠ 0 induces reduced values of σN, which has serious implications on transport across the NS interface. For instance, the Andreev reflection probability A, which is proportional to |$\sigma _{N}^{2}$|, is seriously suppressed for Z = 2 except at the gap edges E = ±Δ0; see the blue curve in Fig. 25(B). In contrast, the normal reflection probability B is dominant at almost all energies except for E = ±Δ0; see the black dotted curve in Fig. 25(B). We finally point out that in both cases of Z = 0 and Z ≠ 0, the electron and hole transmission probabilities vanish, where C = 0 and D = 0 for |E| < Δ0, an effect that occurs because quasiparticles cannot enter the superconductor region as a traveling wave; see the red and green curves in Figs. 25(A), (B).
Scattering processes at the interface of a 1D NS junction along x, where the direction of the arrow denotes the direction of the group velocity and the interface is located at x = 0: an incident electron (solid dark blue arrow) from the normal metal can be reflected back into the normal metal as an electron (solid dark blue arrow) and also as a hole (dotted black arrow), known as normal and Andreev reflections, respectively. Moreover, the incident electron can be transmitted into the superconductor as electron-like (solid light blue arrow) and hole-like (dotted red arrow) quasiparticles, processes known as normal transmission and crossed Andreev reflection, respectively.
Scattering probabilities across an NS interface with S being a conventional spin-singlet s-wave superconductor. (a) Andreev reflection probability, (b) normal reflection probability, (c) transmission probability of an electron-like quasiparticle, (d) transmission probability of a hole-like quasiparticle. (A) Z = 0, (B) Z = 2.
In what follows, we inspect the tunneling conductance, which can be obtained using the probabilities given by Eqs. (151) and (152) using the so-called BTK formula [285]. The essence of this formula is that the conductance of the junction σS is determined by the interface resistance. Thus, if the number of channels through the interface is fixed, the resulting conductance is determined by the effective transparency at the interface. At sufficiently low temperatures, the differential conductance dINS/dV is determined by the tunneling conductance E = eV is given by
where the conductance is obtained as σS(E) = 1 + A − B. Then, by using A and B from Eqs. (151) and (152), we get
where Ω is given by Eq. (147). In Fig. 26, we present the normalized conductance obtained here as a function of eV for different values of Z that describe full, intermediate, and tunnel transparencies of the interface. For σN = 1 (Z = 0), the conductance is quantized for subgap energies such that σS(E) = 2 is satisfied for |E| ≤ Δ0; see the solid black curve in Fig. 26. By taking account of the spin degeneracy, having σS(E) = 2 means that the differential conductance in Eq. (155) reaches 4e2/h. As Z takes finite values, σN reduces, which produces lower conductance values, as depicted by the blue dotted and red dashed curves in Fig. 26. For sufficiently large Z, the normal transmissivity σN is very small inside the superconducting gap, as a result of σS(E) being proportional to |$\sigma _{N}^{2}$|. Thus, the differential conductance at subgap energies can acquire values within the range of 0 up to 4e2/h. Note also that for energies outside the superconducting gap, |E| > Δ0, the conductance σS(E) approaches the normal state value σN. To summarize our discussion for very large Z, and thus low normal transmissivity (σN ≪ 1), it is possible to obtain a simple expression for σS(E), which reads
This expression means that tunneling conductance normalized by its value in the normal state expresses the DOS of a bulk superconductor in the low transparent junction and is insensitive to the phase of the pair potential. In what follows, we will see that this situation is dramatically changed for NS junctions with unconventional superconductors.
Normalized conductance across an NS junction with S being a conventional spin-singlet s-wave superconductor. The different curves correspond to (a) Z = 0, (b) Z = 0.5, and (c) Z = 2.
5.2. Tunneling spectroscopy in unconventional superconductors
In this part we analyze conductance in NS junctions with N being a semi-infinite normal metal and S being a semi-infinite unconventional superconductor. In general, we proceed as in the previous subsection taking into account that unconventional superconductors have anisotropic pair potentials. We consider a 2D NS junctions along x, with a flat interface located at x = 0 and N (S) defined for x < 0 (x > 0); see Fig. 27. Moreover, we assume that the interface is determined by an insulating barrier, which here we model by a delta barrier potential Hbδ(x); the perpendicular direction along y is translationally invariant. With this information, we are now in a position to employ the methodology presented in Section 2 to obtain the scattering states whose coefficients will determine the charge transport in the same spirit as in the previous subsection.
Scattering processes at a 2D normal metal/unconventional superconductor junction along x, where the direction of the arrow denotes the direction of the group velocity and the interface is located at x = 0: an incident electron with incident angle θ from the normal metal can be reflected back into the normal metal as an electron and also as a hole, processes known as normal and Andreev reflections. Moreover, the incident electron can be transmitted into the superconductor as electron-like and hole-like quasiparticles, processes known as normal transmission and crossed Andreev reflection, respectively. The pair potential of the unconventional superconductor is parametrized by the incident angle as θ+ = θ, θ− = π − θ, which represent pair potentials for electron-like and hole-like quasiparticles.
Since the NS junction now is a 2D system, there is an incident angle with respect to the x-axis, denoted by θ, at which quasiparticles hit the NS interface; see Fig. 27. Now, we use the BdG equations in the quasiclassical limit to obtain the scattering states and inspect conductance across the NS interface; see Eqs. (32) and the discussion at the beginning of Section 2.4. In this limit, for an incident electron, the total wavefunction becomes
where u±, v±, γ0, and γ± are given by Eqs. (34), (35), and (36). The coefficients in Eqs. (158) represent the following processes: a, b, c, and d correspond to Andreev reflection, normal reflection, normal transmission, and Andreev transmission (or crossed Andreev reflection), respectively. To fully determine the wavefunction Ψ(θ, x), we still need to find the coefficients a, b, c, and d, which is done by employing the corresponding boundary condition at the NS interface. This condition is given by
Then, by taking kFx = kFcos θ ≫ |γ0|, kFx ≫ |γ±|, the coefficients a(E, θ), b(E, θ), c(E, θ), and d(E, θ) acquire simpler expressions that read
where Zθ = Z/cos θ, with Z = mHb/(ℏ2kF) and
After finding the coefficients of the total wavefunction in Eq. (158), we are now in a position to obtain the tunneling conductance, which now needs to be averaged over the angles θ. Thus, tunneling conductance, normalized by its value in the normal state, is obtained as
where σS and σR are obtained as
and
with E = eV [48]. In this formula, σS(E, θ) expresses the effective transmissivity in the superconducting state. An interesting fact to note in Eq. (164) is that it reveals the formation of the SABSs that we studied in Section 2.4. In fact, the condition for the formation of SABSs is given by the zeros of the conductance denominator in Eq. (164) when the N and S regions are separated σN(θ) → 0 [41,48,49]. This permits us to find the bound state condition equal to 1 = Γ+Γ−, which is exactly the same as the one given by Eq. (40) found in Section 2.4. Thus, the bound state energy found from the zeros in the denominator of the conductance is |$E_{b}=\Delta \, {\rm cos}(\delta _{\phi }/2)\, {\rm sgn}[{\rm sin}(\delta _{\phi }/2)]$|, which is the same as the expression given by Eq. (41) and we have assumed |Δ±| = Δ0, with Δ± = Δ(θ±); here also δϕ = ϕ− − ϕ+ and exp(iϕ±) = Δ±/|Δ±|. See Section 2.4 for a detailed discussion on the SABSs.
It is useful now to express the conductance σS given by Eq. (164) in terms of the energy of the SABS given by Eq. (41). For simplicity, we denote |Δ±| = Δ0 and consider the regime with |E| ≤ Δ0. Hence, we obtain
where |$\gamma =\sigma _{N}(\theta )/(2\sqrt{1-\sigma _{N}(\theta )})$| [248]. Then, for a SABS with zero energy, namely, Eb = 0, the conductance given by Eq. (165) becomes
Therefore, when the energy E ≡ eV is tuned to the energy of the zero-energy SABS, namely, when E = 0, we see that Eq. (166) becomes σS(E = 0, θ) = 2, irrespective of the value of σN(θ). Moreover, the width of resonance at zero energy decreases with the decrease of σN(θ) [41,48,49]. Thus, the conductance expression given by Eq. (166) reveals a remarkable property of tunneling spectroscopy of unconventional superconductors in the presence of ZESABSs. We also anticipate that this expression also captures the essence of tunneling spectroscopy in the presence of MZMs [292,293], as we will discuss later.
5.2.1. For spin-polarized p-wave superconductors
In the case of spin-polarized p-wave superconductors, we now know that they host MZMs at their edges in the 1D case [31]. Thus, in an NS junction formed by these types of superconductors and a normal metal, tunneling occurs into an MZM at the interface. To the impact of the MZM on conductance, we only consider the contribution from θ = 0, which simplifies Eqs. (160) and leads to scattering probabilities given by
In contrast to the obtained probabilities for NS junctions with conventional spin-singlet s-wave pair potential and given by Eqs. (151)–(154), for p-wave superconductors in Eqs. (167) there is an additional − in front of the factor (1 − σN) in the denominator. Moreover, the sign in front of ∣Γ∣2 in the numerator of B is changed from − to +. As we will see below, the impact of the MZM located at the NS interface has its origin in these changes of signs. In Fig. 28, we plot the coefficients A, B, C, and D as a function of energy for distinct values of Z. The first observation is that the probabilities for Z = 0 in Fig. 28 behave exactly the same as in NS junctions with s-wave superconductors shown in Fig. 25. Interestingly, at Z = 2, the Andreev and normal reflection probabilities exhibit a sharp peak and dip around zero energy eV = 0, respectively. This behavior is in stark contrast to what occurs in NS junctions with s-wave superconductors. We highlight that irrespective of the value of Z, there is a perfect Andreev probability with A = 1 at eV = 0, and a zero normal reflection probability B = 0.
Scattering probabilities for an NS junction with a p-wave superconductor (167) for Z = 0 (A) and Z = 2 (B). The different curves represent the following processes: (a) Andreev reflection probability, (b) normal reflection probability, (c) transmission probability of an electron-like quasiparticle, (d) transmission probability of a hole-like quasiparticle for the spin-triplet p-wave case. For Z = 0, B = D = 0 independent of E.
Next, we calculate the conductance σS(eV) = 1 + A − B, from which, by making use of the probabilities A and B given by Eqs. (167), we obtain
This expression represents the conductance across NS junctions formed by a metallic lead N and a semi-infinite spin-polarized p-wave superconductor [25], which was initially predicted in the tunneling conductance of unconventional pairing with ZESABSs in Refs. [41,48]. See also Refs. [75–77]. To see the properties that Eq. (168) reveals, it is instructive to compare Eq. (168) with the conductance expression for NS junctions with s-wave superconductors given by Eq. (156). For σN = 1, we find no difference between the two expressions given by Eq. (156) and Eq. (168). Remarkably, the difference appears with the decrease of σN. In the case of an NS junction with p-wave superconductors, the conductance σS(E) develops a peak around E = eV = 0 for σN < 1, giving rise to a quantized value of σS(E = 0) = 2; see Fig. 29. The fact that there is a factor of 2 is a consequence coming from the spin-polarized nature of the junction. This zero-energy peak (ZEP), sometimes referred to as zero-bias peak (ZBP), is a manifestation of the emergence of an MZM at the interface and, therefore, it represents a strong Majorana signature. It is worth pointing out that this zero-bias conductance peak (ZBCP) stems from having a perfect Andreev reflection probability, as shown by the blue curve in Fig. 28(B). Moreover, the height of the ZBCP of normalized conductance σT(eV) = σS/σN increases with the increase of Z (decrease of σN), as depicted by the red curve in Fig. 29. We can therefore conclude that MZMs at the edge of topological superconductors produce ZBCPs and these peaks can be seen as robust evidence of Majorana physics. A similar ZBCP is expected at the right end of the superconductor, which should emerge as well due to the spatial nonlocality of MZMs. The detection of ZBCPs at both ends is thus necessary and it must comply with the nonlocality effects shown by MZMs. These features are of course complementary to the peaks at zero energy shown in the LDOS and odd-frequency pairing in Section 3. We further point out that ZBCPs can also be taken as evidence of having the dominant OTE pairing [116]; see also Refs. [42,94,113].
Normalized conductance across an NS junction with a p-wave superconductor by its value in the normal state. (a) Z = 0, (b) Z = 0.5, (c) Z = 2.
It is important to mention that, after conceiving the physical realization of spin-polarized p-wave superconductivity in superconductor–semiconductor hybrid systems, there has been an impressive number of theoretical works predicting ZBCPs as signatures of MZMs and topological superconductivity (see e.g. Refs. [146–156]) and nonsinusoidal current–phase curves [157–165]. Along these lines, experimental advances have also reported interesting results [153,166–174], although the unambiguous identification of MZMs still remains a challenging task [154]. The main difficulty in this case is that ZBCPs similar to those expected due to MZMs in Eq. (168) also appear in the presence of topologically trivial ZESABSs [147,154,176,179–214]. These trivial ZEABSs and their associated ZBCPs, however, do not share all the properties of MZMs. Thus, exploiting intrinsic Majorana properties, such as their spatial nonlocality, could open a way to distinguish them from trivial zero-energy states [199,202,204,206,215–221]. Despite the evident open questions, nevertheless, the experimental progress made in recent years, as well as the ongoing efforts following a bottom-up approach to fabricate minimal Kitaev chains [224–226], is very likely to have a positive impact on the realization, detection, and application of MZMs in the near future.
5.2.2. For spin-singlet d-wave superconductors
For NS junctions with a superconductor having a spin-singlet d-wave pair potential, we consider that the d-wave pair potential is expressed by the angle measured from the a-axis of the basal plane since that of the lobe direction coincides with that of the a-axis. We also consider a junction with a misorientation angle α, which denotes the angle between the normal to the interface and the a-axis of a d-wave superconductor. A schematic of the surface of α = 0 and α = π/4 is shown in Fig. 3. In this case, the pair potential Δ± felt by quasiparticles with an injection angle θ is given by
Having defined the pair potential, we plot the normalized tunneling conductance σT(eV) in Eq. (162) as a function of eV in Fig. 30 for distinct misorientation angles α and finite Z (small σN(θ)) [48]. For α = 0, the normalized conductance exhibits a finite value having a V shape and reaching its minimum at eV = 0; see Fig. 30(A). While intermediate values of Z induce finite minimum in the conductance at zero energy, large values of Z produce almost vanishing conductance at zero energy; in this case, the line shape of σT(eV) is almost similar to that by the LDOS of a bulk d-wave superconductor with the increase of Z, as can be seen by comparing the different curves in Fig. 30(A). Thus, no ZBCP appears, which is consistent with the absence of ZESABS at α = 0; see Section 2.4.1. Interestingly, for nonzero α, a ZBCP in σT(eV) always appears and it becomes prominent with the increase of Z. Its peak height and width becomes maximum for α = π/4, as can be clearly seen in Figs. 30(B), (C); see Refs. [41,48,49]. The ZBCP α = π/4 is a consequence of the formation of a ZESABS, which, for this misorientation angle, forms a zero-energy flat band for any θ. It is worth noting that a similar zero-energy flat band appears for the spin-triplet px-wave pairing for both the Sz = 0 and Sz = 1 cases where Δ(θ±) is given by Δ(θ±) = ±Δ0cos θ; see Refs. [70,87,294]. Before closing this part, we would like to briefly mention that the ZBCP discussed in this part has been observed in many experiments on tunneling spectroscopy of high-Tc cuprates, as reported in Refs. [16,50–54,56,57,295–297], and also in other strongly correlated superconductors like CeCoIn5 [60] and PuCoIn5 [61].
Normalized conductance across a junction formed with a d-wave superconductor. The different panels correspond to different misorientation angles such that A: α = 0, B: α = π/8, C: α = π/4. Moreover, the solid blue, dashed red, and solid black lines correspond to distinct values of Z, where (a) Z = 0.5, (b) Z = 2, and (c) bulk DOS.
5.3. Summary
This section has focused on analyzing tunneling conductance across NS junctions based on unconventional superconductors. In particular, we discussed that tunneling into a spin-polarized p-wave superconductor produces a quantized ZBCP with high 2e2/h as a strong signature of the emergence of an MZM at the end of the system in the topological phase. We have shown that this conductance quantization at zero energy is robust against interface imperfections and results from perfect Andreev reflections, an effect that does not occur in spinful junctions with conventional spin-singlet superconductors where conductance can take values between 0 and 4e2/h. Then, we have addressed how tunneling into the surface of a spin-singlet d-wave superconductor can also give rise to a pronounced peak at zero energy due to the formation of a ZESABS, which, however, depends on the misorientation angle. We have also highlighted the intense theoretical and experimental efforts in superconductor–semiconductor hybrids aiming to detect ZBCPs, which, despite the challenges, remain one of the most promising scenarios to realize and detect Majorana physics and topological superconductivity. In this section, we study tunneling conductance for spin-triplet and spin-singlet pairing separately. The existence of ZBCP by ZESABS has been predicted in noncentrosymmetric superconductor junctions with mixed parity pairing states [298–301]. ZBCP has been reported in a recent experiment with AuSn4 where mixed parity surface superconductivity is promising [302].
6. Detecting Majorana zero modes via the anomalous proximity effect
In this section, we consider NS junctions formed by a finite diffusive N region, denoted as DN, and a semi-infinite unconventional superconductor, and then investigate signatures of MZMs and ZESABSs in the anomalous proximity effect. Before going any further, it is important to clarify what we mean by anomalous proximity effect. Generally, the superconducting proximity effect occurs when a superconductor, placed in proximity to another material such as a normal metal, induces superconducting correlations or Cooper pairs, which can have completely different symmetries associated with the nature of the normal metal and interface. Thus, nonsuperconducting materials can acquire superconducting properties, with the coupling between the two systems determining the amount of proximity-induced Cooper pairs. Thus, the proximity effect is present in the NS junctions discussed in Section 5, where proximity-induced Cooper pairs in the ballistic N region are naturally characterized by Andreev reflections [288,289,303,304]. To detect the proximity effect, conductance and LDOS in N are commonly carried out, which are then able to reveal the properties of the induced superconductivity. At this point, we note that, when the superconductor is coupled to diffusive normal metal DN, the resistance of the DN is influenced by the penetration of Cooper pairs due to the proximity effect. Thus, the LDOS and conductance exhibit a behavior that is different from junctions with ballistic N metals. For these reasons, the changes in the penetration of Cooper pairs, LDOS, and conductance in DNS junctions are referred to as the anomalous proximity effect [80–82], which is the focus of this part for unconventional superconductors [305]. For pedagogical purposes, we first address the proximity effect in conventional spin-singlet s-wave superconductors and then inspect the anomalous proximity effect in unconventional superconductors.
6.1. Proximity effect in conventional superconductors
We consider a DNS junction along x formed by a DN metal of length L and a semi-infinite conventional superconductor with uniform spin-singlet s-wave pair potential in S and then explore charge conductance with an N electrode attached to DN. The electrode, the DN metal, and the superconductor are located at x < −L, −L < x < 0, and x > 0, respectively, with an insulating interface barrier at x = 0 given by a δ-function. Here, DN is assumed to be in the dirty limit where the mean free path ℓ is much smaller than the thermal diffusion length |$\sqrt{\hbar D/(2\pi k_{B} T)}$| with Boltzmann constant kB, where D is the diffusion constant and T is the temperature. Within this energy scale, the coherence between electrons and holes is good and the retroreflectivity of the Andreev process is satisfied. Moreover, D introduces an energy scale known as the Thouless energy ETh = ℏD/L2. Moreover, we consider L ≫ ℓ, and its resistance in the normal state is denoted as Rd. Taking all this information into account, the calculation of conductance is carried out using the quasiclassical Green’s function framework [305].
Under these conditions, the conductance is obtained from the charge current [306,307]. In our case, probed at the left side of DN, it is given by [305]
with R being the resistance of the junction at any temperature and ft is given by
where V is the bias voltage. Here, ft(x = −L) is the Fermi distribution function at the left interface where the electrode is attached without any barrier. At low temperatures, the resistance of the junction R is given by
where the first (second) term in the expression for R describes the resistance contributions of the interface (diffusive normal metal DN). Here, Rb and Rd represent the normal state resistances at the interface and in DN, respectively, while ζIm is the imaginary part of ζ(x) that characterizes the proximity effect, namely, ζIm(x) = Im[ζ(x)]. It is noted that cos ζ(x) and sin ζ(x) express the quasiclassical Green’s function. The current IK and ζ(x) are found from
where ζS(N) are solutions ζ(x) in S(N); note that ζ(x) is a function characterizing the proximity effect and is obtained by solving the Usadel equation described above Eq. (172); see also Refs. [305,308] for more details. Then, using the knowledge about resistance in Eqs. (171), conductance in the junction can be found as its inverse. In particular, the conductance, normalized by its value in the normal state, is given by
where RN(S) represents the resistance of the junction in the normal (superconducting) state at sufficiently low temperatures obtained from Eqs. (171).
To visualize the behavior of the conductance given by Eq. (173), in Fig. 31 we plot it as a function of eV for distinct values of Rd/Rb at two representative large and low values of ETh. For completeness, in Fig. 32 we also show the normalized LDOS ρ normalized by its value in the normal state in the middle of DN at x = −L/2 as a function of E = eV for the same cases of Fig. 31. We note that ρ is obtained by ρ = Real[cos ζ] using ζ(x) from the same quasiclassical Green’s functions employed to find the charge current in Eq. (170); see Ref. [305] for more details. At large Thouless energies with ETh = 0.5Δ0, the conductance σT(eV) has an almost U-shaped structure, acquiring lower values when the ratio Rd/Rb is small; see Fig. 31(A). Note that the conductance at Rd/Rb = 0.2 in Fig. 31(A) has a similar behavior to that in ballistic junctions at low transparencies presented in Fig. 26. As the ratio Rd/Rb increases, the subgap conductance values also acquire higher values but still maintain a U shape, with a minimum near eV = 0 that is only developed for very large Rd/Rb; see the blue and green curves in Fig. 31(A). When it comes to the LDOS in this regime with ETh = 0.5Δ0, it has a clear U-shaped profile at small ratio Rd/Rb but, as the ratio increases, the LDOS gets suppressed at low voltages and develops a clear minimum around eV = 0; see Fig. 32(A). We can thus conclude that in this regime, the penetration of Cooper pairs into DN increases with the increase of Rd/Rb that then enhances conductance, which is, however, not visible in the LDOS.
Normalized conductance in DNS junctions where S is a spin-singlet s-wave superconductor and DN is a diffusive normal metal. Left and right panels correspond to distinct values of the Thouless energy: A: ETh = 0.5Δ0. B: ETh = 0.05Δ0. The different curves correspond to distinct values of Rd/Rb: (a) Rd/Rb = 0.2, (b) Rd/Rb = 1, (c) Rd/Rb = 2.
LDOS of the diffusive normal metal DN at x = −L/2 for a DNS junction, where S is a spin-singlet s-wave superconductor. Left and right panels correspond to distinct values of the Thouless energy: A: ETh = 0.5Δ0. B: ETh = 0.05Δ0. The LDOS is normalized by its value in the normal state. The different curves correspond to distinct values of Rd/Rb: (a) Rd/Rb = 0.2, (b) Rd/Rb = 1, (c) Rd/Rb = 2.
In contrast to large Thouless energies ETh = 0.5Δ0 discussed in the previous paragraph, the conductance σT(eV) at low Thouless energies ETh = 0.05Δ0 develops a very different behavior, as can be seen in Fig. 31(B). In fact, σT(eV) exhibits a peak at zero voltage bias for Rd/Rb = 0.2 and 1, as depicted by the red and blue curves in Fig. 31(B). It is interesting to note that this ZBCP becomes prominent at Rd/Rb = 1, although it does not exhibit any robustness as its high strongly depends on the ratio Rd/Rb [306,307]. The ZBCP can be understood as follows: When the energy of an injected electron in DN is near the Fermi surface, the reflected hole by Andreev reflection interferes with the electron due to retroreflectivity, forming a standing wave that helps electrical conduction [309,310]. As a result, the conductance increases within the range in which the energy of the incident electron is smaller than the Thouless energy as measured from the Fermi surface. In contrast, when the energy of the incident electron exceeds the Thouless energy, the conductance rapidly decreases because standing waves are no longer formed. Moreover, the height of the ZBCP has been shown to never exceed unity [306,307] and it is very sensitive to the applied magnetic field [306]. We can thus see that the ZBCP in Fig. 31(B) has a distinct origin to those ZBCPs due to MZMs and ZESABSs of unconventional superconductor junctions discussed in Section 5. Moreover, the corresponding LDOS ρ intriguingly shows a dip-like (gap-like) structure around eV = 0 for all cases of the ratio Rd/Rb, reaching a minimum at eV = 0 independent of the value of Rd/Rb. We can thus conclude that strong Thouless energies promote an enhanced penetration of Cooper pairs around zero energy, which is reflected in the large value of conductance at zero voltage. This enhancement is, however, not seen in the LDOS, which we interpret to happen because the LDOS does not capture the interference effects among quasiparticle paths in DN, in contrast to what conductance does.
6.2. Anomalous proximity effect in unconventional superconductors
In the case of DNS junctions with unconventional superconductors, special care needs to be taken due to the injection angles θ. As before, conductance is obtained from the resistance R, which here reads [81]
where ζIm(x) is a function characterizing the proximity effect, obtained using Eq. (172) with appropriate boundary conditions, while 〈IK〉 represents the angular average over many channels of IK. For spin-singlet superconductors, IK is given by [311]
with
On the other hand, for spin-triplet superconductors it is given by [80,81,311]
with C1 given by
Here, σN = 1/(1 + Z2), Γ± = Δ±/(E + Ω±), ζNr, ζNi denotes the real and imaginary parts of ζN, which characterizes the proximity-induced superconducting correlations in N; see Eq. (172) and Refs. [305,311] for details. We remember that having ζN = 0 indicates an absence of the proximity effect in DN. It can be verified that for a DN region of zero length, L = 0, then Rd = 0, ζN = 0, ζNr = 0, and ζNi = 0 are satisfied. Then, the resistance of the junction is determined only at the interface and we reproduce the results in ballistic junctions in Eqs. (162)–(164). Therefore, in general, Eq. (174), combined with Eq. (176) and Eq. (178), also including Eqs. (172), enables the calculation of conductance for DNS junctions with spin-singlet and spin-triplet unconventional superconductors as given by Eq. (173) but with the resistance here obtained from Eqs. (174).
6.2.1. Anomalous proximity effect in 1D spin-polarized p-wave superconductors
In this subsection, we consider the 1D limit of p-wave superconductor junctions [311]. This corresponds to the Kitaev chain in the continuous limit of the topological phase, where the pair potential is expressed by Δ+ = −Δ− = Δ0. Since Γ+ = −Γ− = Γ with Γ = Δ0/(E + Ω), the expression for IK from Eq. (178) acquires a simpler form given by [311]
with
Then, 〈IK〉 is obtained as 〈IK〉 = IK/σN and is used in Eq. (174) to finally find the normalized conductance σT(eV) following Eq. (173). In Fig. 33 we present σT(eV) at Z = 1.5 for distinct values of the ratio Rd/Rb but fixed Rb. As a reference, we plot σT(eV) with Rd = 0 corresponding to the ballistic junction case discussed in Fig. 29. First, at low Thousless energies such as ETh = 0.05Δ0 in Fig. 33(A), we identify three clear voltage regions with distinct conductance profiles. At eV ∼ 0, the normalized conductance σT(eV) exhibits a peak that increases as Rd/Rb increases, reaching a maximum value at Rd/Rb = 1. For voltages in the interval 0.05Δ0 < |eV| ≲ 0.4Δ0, the conductance develops a maximum value for Rd/Rb = 0 while for Rd/Rb = 1 it is seriously suppressed as compared to the ballistic case with Rd/Rb = 0. For |eV| > 0.4Δ0, the conductance becomes maximum for Rd/Rb = 1 again, as depicted by the red curve in Fig. 33(A). In the case of larger Thouless energies shown in Fig. 33(B) for ETh = 0.5Δ0, we observe an overall similar behavior at lower ETh. For instance, by a direct inspection, we observe that in the ballistic regime Rd/Rb = 0, there is no change between the conductance at ETh = 0.05Δ0 and ETh = 0.5Δ0 shown by the green curves in Figs. 33(A) and (B) respectively. For Rd/Rb ≠ 0, the ZBCP at ETh = 0.5Δ0 undergoes a slight reduction in its sharpness in contrast to what is seen at ETh = 0.05Δ0 but the overall profile is more or less preserved; see curves (a)–(c) in Figs. 33(A), (B).
Normalized conductance σT(eV) in a DNS junction where S is a 1D spin-polarized p-wave superconductor and DN is a diffusive normal metal. Here, σT(eV) is normalized by its value in a normal state for Z = 1.5. Left and right panels correspond to distinct values of the Thouless energy: A: ETh = 0.05Δ0, B: ETh = 0.5Δ0. The different curves in each panel correspond to distinct values of Rd/Rb: (a) Rd/Rb = 0.1, (b) Rd/Rb = 0.3, (c) Rd/Rb = 1, and (d) Rd/Rb = 0 (ballistic).
In order to further understand the proximity effect in the diffusive region DN, we plot the LDOS ρ(E) normalized by its value in the normal state at the middle of DN in x = −L/2 in Fig. 34. For all cases, the resulting ρ(E) has a sharp ZEP [80,81,311] and is completely different from the standard proximity effect in a conventional superconductor junction, which commonly has a zero-energy dip (gap) as seen in Fig. 32. For the ballistic regime with Rd = 0, the normalized LDOS ρ(E) = 1 is satisfied due to ζ(x) = 0, which indicates the absence of the proximity effect. The height of the ZEP becomes larger with the increase of Rd/Rb, as seen in curve (c) of Figs. 34(A) and (B). The functionalities of the ZEP can be understood by solving the Usadel equations (172) for ζ(x), with appropriate boundary conditions, which then gives the normalized LDOS ρ(E) = Re[cos ζ(x)] [80,81]. Thus, E = 0 gives ρ(0) = cosh[2Rd/(RbσN)], which clearly demonstrates that the ZEP in the LDOS depends on Rd, Rb, and σN. It is therefore evident that ρ(0) gets enhanced with the increase of Rd/Rb at fixed values of σN as happens in Fig. 34. Furthermore, it is worth pointing out that the LDOS for ETh = 0.05Δ0 exhibits a dip-like structure at E = ±Ed with 0.05Δ0 < Ed < 0.1Δ0; see Fig. 34(A), with a prominent size for Rd/Rb = 1, but this, however, tends to disappear for ETh = 0.5Δ0. The fact that the LDOS develops a ZEP in the DN region as a result of the proximity effect due to the p-wave superconductor makes it different from the dip-like structure seen in conventional junctions and this is why this effect is called the anomalous proximity effect [80–82,88].
LDOS ρ(E) at x = −L/2 in a DNS junction where S is a 1D spin-polarized p-wave superconductor and DN is a diffusive normal metal. Here, ρ(E) is normalized by its value in the normal state for Z = 1.5. Left and right panels correspond to distinct values of the Thouless energy: A: ETh = 0.05Δ0, B: ETh = 0.5Δ0. The different curves in each panel correspond to distinct values of Rd/Rb: (a) Rd/Rb = 0.1, (b) Rd/Rb = 0.3, (c) Rd/Rb = 1.
The anomalous proximity effect in this case is also characterized by having an induced pair amplitude in the DN region with odd-frequency spin-triplet s-wave symmetry [42,85,86]. Since the edge of the p-wave superconductor hosts an MZM, the anomalous proximity effect can be seen as another Majorana effect that, in turn, can help detect MZMs [42,311]. The odd-frequency pair amplitude s1 is plotted as a function of E in Figs. 35 and 36 for ETh = 0, 05Δ0 and ETh = 0.5Δ0, respectively. In this case, the real part Re(s1) becomes an odd function of E while the imaginary part Im(s1) becomes an even function of E. The magnitudes of both Re(s1) and Im(s1) are enhanced around E = 0 [85]. This means that the ZEP and MZM always accompany the formation of odd-frequency pairing [42,85,86,96,116].
Energy dependence of the OTE pair amplitude at ETh = 0.05Δ0 in a DNS junction where S is a spin-polarized p-wave superconductor and DN a diffusive normal metal. Panels A and B show the real and imaginary parts, respectively. The different curves correspond to different values of Rd/Rb: (a) Rd/Rb = 0.1, (b) Rd/Rb = 0.3, (c) Rd/Rb = 1.
Energy dependence of the OTE pair amplitude at ETh = 0.5Δ0 in a DNS junction where S is a spin-polarized p-wave superconductor and DN a diffusive normal metal. Panels A and B show the real and imaginary parts, respectively. The different curves correspond to different values of Rd/Rb: (a) Rd/Rb = 0.1, (b) Rd/Rb = 0.3, (c) Rd/Rb = 1.
Finally, we address the perfect resonance of charge conductance at zero voltage. For this purpose, we look at σS(eV) = 1/R, which is not a normalized conductance, given by
As shown in Fig. 37, the zero bias conductance σS(0) is quantized and its value is 2e2/h independent of any value of Rd/Rb and ETh. This is a very interesting and robust result that can be shown analytically [80,81,116].
Charge conductance σS(eV) in DN/p-wave superconductor junctions. (a) Rd/Rb = 0.1, (b) Rd/Rb = 0.3, (c) Rd/Rb = 1. A: ETh = 0.05Δ0, B: ETh = 0.5Δ0.
In fact, for a spin-triplet p-wave superconductor, it is possible to show that
where |$\zeta _{Ni}=\frac{2R_{d}}{\sigma _{N}R_{b}}$|. Then, the total resistance of the junction R at eV = 0 becomes
with R0 = RbσN. In the present case, R0 becomes h/e2 since the number of channels is only one. The total resistance of the junction then acquires the value R = R0/2 and we obtain
In Eq. (183), the first term comes from the resistance in DN and it reduces to R0/2 for Rd → ∞. On the other hand, the second term is from the resistance at the interface and it reduces to zero at the large-Rd limit. It is remarkable that the obtained conductance (resistance) at zero voltage is completely independent of Rd and Rb. This is a strong feature of the proximity effect of MZM into DN [42,80,95,116]. The quantization of σS(eV) can also be explained from the viewpoint of the index theorem [312].
Finally, we would like to address the question that the observed quantization is robust against the boundary condition at the left interface between the electrode and DN. In the present analysis, the resistance between the electrode and the DN is set to be zero. It is possible to extend this in the case where the resistance at the interface at x = −L is |$R_{b^{\prime }}$| and the trasmissibity at the interface is σN2, given by
with |$R_{b^{\prime }}=R_{0}/\sigma _{N2}$|.
The total resistance of the junction is then given by
and, after a cumbersome calculation, it is possible to obtain, at E = 0,
Then, by plugging Eqs. (186) into Eq. (185), we obtain R = R0/2 and from Eq. (184) it is not difficult to see that the perfect conductance quantization due to MZMs is verified. We note that the anomalous proximity effect is absent in DNS when S is a d-wave superconductor [78,79]; in this case, the resulting proximity effect becomes the conventional proximity effect realized in junctions between DN and spin-singlet s-wave superconductors. It is worth noting that the anomalous proximity effect has also been studied in noncentrosymmetric superconductor junctions with coexistent spin-singlet s-wave and spin-triplet p-wave pair potentials [311,313], which showed that the anomalous proximity effect remains when the p-wave component is dominant. It is thus possible to establish a direct relationship between the robustness of the anomalous proximity effect and the p-wave nature of superconductivity that gives rise to MZMs.
Before closing this section, we would like to briefly mention the appearance of the anomalous proximity effect in the zero-voltage charge conductance in 1D p-wave superconductors, where an MZM appears at the edge. Although ZBCP appears in the nontopological superconductor, the present ZBCP by the anomalous proximity effect shows remarkable features. The height of ZBCP can become extremely large as compared to the background high voltage value above the energy gap. Here, we denote the conductance as a function of eV and T as σS(eV, T) = σS(eV). First, we focus on the voltage dependence for T = 0. As seen from the perfect resonance independent of the resistance in diffusive normal metal RD and that at the insulating barrier Rb [80,81,116],
is satisfied. Here, eVc is sufficiently large compared to Δ0. With the increase of RB + RD, the ratio in Eq. (187) increases monotonically, which is completely different from conventional junctions and can be taken as a signature of MZMs. Interestingly, if we look at the temperature dependence of σS(eV = 0, T),
for T > Tc, we obtain the same expression as Eq. (187). Here, Tc represents the critical temperature of the superconductor where the resistance of the junction is RD + RB. In this regard, temperature dependence is also a good indicator of MZMs. We emphasize that the ZBCP cannot be taken as the only indicator of MZMs; the temperature and voltage dependences of the normalized conductance discussed here also provide complementary evidence of MZMs. The ZBCP and the enhancement of conductance at zero voltage are results of topologically protected MZMs accompanying odd-frequency spin-triplet s-wave pairing [42,80,85].
6.3. Summary
In this part we have discussed charge transport and LDOS in NS junctions formed by a finite diffusive normal metal N and a superconductor S. We have shown that when the superconductor has a spin-singlet s-wave pair potential, the resulting conductance as a function of energy can exhibit a U-shaped profile but it can also have a peak at zero energy due to impurity scattering, while the LDOS in this case can produce a dip-like feature. The ZBCP, however, is always below perfect quantization and it is not robust under the variation of parameters. Interestingly, for spin-polarized p-wave superconductors, the conductance develops a ZBCP that is robust against impurity scattering. Moreover, this quantized ZBCP is also accompanied by a zero-energy peak in the LDOS and dominant odd-frequency spin-triplet pair correlations, which, together with the ZBCP, define the anomalous proximity effect. Thus, the anomalous proximity effect could be taken as an alternative signature to further explore the detection of MZMs and topological superconductivity [42,96].
7. Detecting Majorana states via the Josephson effect
The Josephson effect is one of the most fundamental phenomena in superconductivity and was predicted by Brian Josephson more than 60 years ago [314]. The Josephson effect reveals the flow of Cooper pairs, known as supercurrent, between two superconductors forming a Josephson junction simply due to a finite phase difference between their pair potentials and without any applied voltage bias [315–321]. In 1991, Akira Furusaki demonstrated that the d.c. Josephson effect can be understood as a result of Andreev reflection processes happening on both interfaces, which then naturally account for the transfer of Cooper pairs between superconductors [322]; see also Refs. [323–325]. This study has paved the way to understand the Josephson effect as a result of phase-dependent ABSs [322–325,325–328], which are formed between the superconductors and exhibit unique profiles for superconductors with unconventional pair potentials [41,67,329–332]. Thus, since Cooper pairs are involved in the transport across the junction [12,321], the Josephson effect has the potential to reveal the type of emergent superconductivity. It is worth mentioning that the Josephson effect is also at the core of interesting quantum applications; see Refs. [40,333–340] for more details. In this part, we are interested in exploring the d.c. Josephson effect in junctions based on unconventional superconductors and how ZESABSs and MZMs can be identified using phase-biased Josephson transport.
7.1. Josephson current in 1D spin-singlet s-wave superconductor junctions
To set the pedagogical groundwork, we first focus on the Josephson effect in a 1D Josephson junction formed by a superconductor–insulator–superconductor (SIS) structure. To describe the Josephson effect, we will first obtain the ABSs and then the supercurrent; see Refs. [163,327]. Thus, to find the ABSs we need to employ the BdG equations discussed in Section 2.1 with a spin-singlet s-wave pair potential. The insulating region is modeled by a delta function barrier potential without superconductivity, while the S regions have a finite and constant spin-singlet s-wave pair potential. Thus, the delta potential and the pair potential are given by
where Hb is the strength of the delta barrier, ΔL(R) = Δ0exp (iφL(R)) the pair potential of the left (right) superconductor with φL(R) being its associated superconducting phase, and Δ0 the isotropic spin-singlet s-wave pair potential amplitude. Without loss of generality, we take φL = φ and φR = 0, such that there is a finite phase difference φ across the junction. We focus on subgap ABS energies, which means that wavefunctions have to exponentially decay as x moves away from either side of the junction. In this case, the wavefunction of the BdG equations is given by
where |$k^{\pm }=\sqrt{\frac{2m}{\hbar ^{2}}(\mu \pm \Omega )}$|, and Γ = E/(E + Ω). The coefficients A, B, C, and D are then found from the boundary conditions for Ψ(x) given by
The ABS energies are then obtained from the nontrivial solutions of Eq. (191), in the same way as in Eq. (41) in Section 2.4. Moreover, we consider the limit of μ ≫ Δ0 where wavevectors as simply approximated as the Fermi wavevector k± ∼ kF. Thus, by taking these conditions into account, we obtain the ABS energies located at x = 0 to be given by [323]
where σN is the normal transparency at the interface given by σN = 1/(1 + Z2), and Z = mHb/(ℏ2kF). These ABS energies are degenerate in spin due to the model considered. We note that there are two ABSs, with positive and negative energies. When Z ≪ 1, which corresponds to σN = 1 and defines the fully transparent regime, the ABSs energies read Eb(φ) = ±Δ0cos(φ/2): this means that the ABSs reach zero energy at φ = π. In contrast, for nonzero barrier strengths Z ≠ 0, the ABSs energies φ = π become |$E_{b}(\varphi )=\pm \Delta _{0}\sqrt{1-\sigma _{N}}$|, leading to a finite energy gap.
The supercurrent across the SIS Josephson junction can be obtained from the thermodynamic relation [10]
where Φ = ℏ/2e is the reduced superconducting flux quantum, while F(φ) is the free energy that in general depends on other parameters of the system. Assuming spin degeneracy, the free energy can be calculated as
with log representing the natural logarithm and T the temperature. Thus, by calculating the derivative of F(φ) with respect to φ, we arrive at [322]
where RN = πℏ/(e2σN) is the resistance of the junction in the normal state. For a sufficiently low transparent junction with large Z, RNI(φ) becomes
where the current–phase relation becomes the simple sinusoidal relation given by Ambegaokar and Baratoff [341]. On the other hand, for a fully transparent limit with Z = 0 and σN = 1, RNI(φ) becomes
where the current–phase relation is proportional to sin (φ/2) for −π < φ < π at sufficiently low temperature [342]. Before going further, it is useful to note that Δ0 depends on temperature T obeying
where Δ0(T = 0) is the zero temperature pair potential and TC the critical temperature. Now, to visualize the behavior of the Josephson current, in Fig. 38(A) we show its phase dependence in the tunnel (Z = 3) and transparent (Z = 0) regimes. We clearly observe that I(φ) has a simple sinusoidal behavior at low transparencies, while it develops a strong skewness in the transparent regime, similar to a sawtooth profile; see the red and blue curves in Fig. 38(A). Moreover, the maximum supercurrents IC, known as critical currents, are plotted in Fig. 39(A) as a function of temperature for the tunnel and transparent regimes. Here we notice that RNIC decreases as the temperature increases, acquiring a similar decreasing behavior near the critical temperature; at low temperatures, however, the transparent regime exhibits larger values.
Phase-dependent supercurrents I(φ) in 1D Josephson junctions for A: s-wave superconductors and B: p-wave superconductors. The curves in each panel correspond to (a) Z = 3, (b) Z = 0. Δ0(0) is the magnitude of the pair potential Δ0 at zero temperature.
Temperature dependence of the critical currents IC in 1D Josephson junctions for A: s-wave pair superconductors and B: p-wave superconductors. The curves in each panel correspond to (a) Z = 3 and (b) Z = 0. Δ0(0) is the magnitude of the pair potential Δ0 at zero temperature.
7.2. Josephson current in 1D spin-polarized p-wave superconductor junctions
In this part, we explore the Josephson effect in Josephson junctions between 1D superconductors having spin-triplet p-wave pair potentials. The structure of the junction is the same as in previous subsection described by Eqs. (189) but now with the pair potentials for a 1D p-wave superconductor [294], namely, |$\Delta _{L(R)}=\hat{k}\Delta _{0} \exp (i\varphi _{L(R)})$| with |$\hat{k}=k/k_{\rm F}$|. We then obtain the ABS energies and supercurrent following the same steps as described in the previous subsection but using the BdG equations within the quasiclassical approximation described in Section 2.3. Then, the corresponding wavefunction can be written as
Hence, by imposing the boundary conditions given by Eqs. (191) and for subgap energies with μ ≫ Δ0, the energy of the ABS is given by [294]
where σN is the normal transparency. Interestingly, the junction in the regime under investigation is already topological; see Section 3, which implies that the ABS energies found in Eq. (200) correspond to topological ABSs. These ABSs are gapless no matter what value the normal transparency σN acquires, thus developing a robust zero-energy crossing at φ = π unlike the ABSs in conventional Josephson junctions of the previous subsection; the only effect of σN here is to reduce the ABSs from the gap edges. The zero-energy protected crossing signals the emergence of a pair of MZMs at the inner sides of the junction. These topological ABSs Eb given by Eq. (200) have also been explored due to the hybridization of ZESABSs in d-wave superconductor junctions [41,249,343]. A consequence of the zero-energy crossing is that these topological ABSs exhibit a 4π periodicity with φ, a remarkable property that was initially derived in the context of d-wave superconductor junctions in 1996 [249]. Then, the supercurrent due to the topological ABSs is given by
where RN = πℏ/(e2σN) is the junction’s resistance in the normal state. For σN = 1, Eq. (201) coincides with Eq. (197) for SIS junctions with s-wave superconductors. The phase-dependent supercurrents are presented in Fig. 38(B) in the tunnel and transparent regimes. In Fig. 39(B) we also show the temperature dependence of their associated critical currents. The first observation is that the supercurrent exhibits a strong skewness near φ = ±π as a result of the robust zero-energy crossing of the ABSs and does not depend on the junction transparency, unlike what occurs in conventional Josephson junctions. Thus, the current–phase relation seriously deviates from a simple sinusoidal one and can be approximated by sin φ/2 within −π ≤ φ ≤ π. Nevertheless, the supercurrent remains 2π-periodic as a function of φ; under fermion parity conservation, however, the zero-energy crossing at φ = π is parity protected and the supercurrent becomes 4π-periodic with φ [31,344]. In relation to the temperature dependence, the supercurrent RNI(φ) at sufficiently low σN is proportional to 1/T; see the red curve in Fig. 39(B). At low temperatures, however, the supercurrent RNI(φ) is proportional to |$1/\sqrt{\sigma _{N}}$| [294,343]. Thus, the magnitude of the Josephson current for a topological Josephson junction made of 1D spin-polarized p-wave superconductors is seriously enhanced as compared to conventional Josephson junctions. This remarkable feature was initially reported during 1996–1997 in the context of d-wave superconductor junctions; see Refs. [249,343,345,346]. The critical current enhancement at low temperatures seen in p-wave Josephson junctions, opposite to what occurs in conventional Josephson junctions, is due to the presence of MZMs and due to ZESABSs in d-wave Josephson junctions.
To close this part, we briefly discuss 1D SNS phase-biased Josephson junctions with the ingredients used to engineer p-wave superconductivity, namely, nanowire SOC, s-wave superconductivity, and the presence of a Zeeman field; see Section 3.2. These junctions are commonly modeled by using the tight-binding description discussed in Section 3.2, where the 1D lattice is separated into a middle N region without pair potential and two S regions on the left and right with a finite pair potential. Moreover, a finite phase difference is considered across the junction. For details of the modeling, we refer the reader to Refs. [161,163]. We first analyze the evolution of the low-energy phase-dependent spectrum as the Zeeman field increases, which is presented in Fig. 40 for short and long S regions, both with a very short N region. At zero Zeeman B = 0, a pair of ABSs within the gap Δ emerges, developing a finite gap φ = π, as expected since the considered parameters place the junction away from the Andreev approximation. As B increases, the ABSs split, and one of them reduces in energy as B approaches the critical field Bc, which now is given as |$B_{c}=\sqrt{\mu _{S}^{2}+\Delta ^{2}}$|, where μS is the chemical potential in S and Δ is the pair potential in S. For B < Bc the ABSs do not depend on LS, as seen by comparing Figs. 40(a) and (d). At B = Bc, the gap closes and the S regions undergo a topological phase transition into a topological phase for B > Bc. We see that the gap closing is more noticeable for longer S regions. Interestingly, as the S regions enter the topological phase, the spectrum shows four ABSs, which are now topological, within the lowest gap that is different from Δ; see Fig. 40(c). There appear two levels around zero energy with a weak dependence for all φ, while the other two levels approach zero energy φ = π. Quite remarkably, for longer S regions, the energy levels closest to zero become flat at zero energy, while the other subgap pair of levels really reach zero energy at φ = 0; see Fig. 40(c). It has been shown that the dispersionless energy levels define two MZMs at φ = 0 located at the outer ends of the S regions, while the four topological ABSs define four MZMs at φ = π, with two additional MZMs located at the inner sides of the S regions; see Refs. [157,161,163]. We stress that the dependence of the energy splitting at around zero energy on LS reflects its intrinsic nature, that it comes from MZMs located at the ends of the S regions. By increasing LS, the wavefunctions of the MZMs overlap less and become more closer to zero energy. It is thus possible to see the length dependence of the zero-energy splitting as a consequence of the Majorana nonlocality.
Phase-dependent low-energy Andreev spectrum in an SNS Josephson junction based on nanowires with SOC for S region length with LS = 2 μm (a)–(c) and LS = 10 μm (d)–(f). The different panels show the evolution with the Zeeman field B. The red curves in (c), (f) indicate the four lowest ABSs giving rise to four MZMs at φ = π for long S regions. Parameters: αR = 20 meV nm, Δ = 0.25 meV, μN(S) = 0.5 meV, LN = 20 nm, LS = 2 μm.
In relation to the supercurrents, in Fig. 41(a) we present the current–phase curves I(φ) in the topological phase for different LS. In Fig. 41(b) we show the critical current IC (maximum supercurrents) as a function of the Zeeman field for different values of the normal transmission, from the tunnel up to the full transparent regime. The first feature to notice is that the current–phase curves in the topological phase B > Bc are 2π-periodic. Moreover, they exhibit a strong dependence on the length of the S regions, developing a sawtooth profile near φ = π for very long S regions. This sawtooth profile has been shown to be robust against variations of the normal transmission, temperature, and also against scalar disorder [161,163], revealing that measuring the length dependence of I(φ) in realistic nanowire-based Josephson junctions is an unambiguous indicator of Majorana physics. When it comes to the critical currents, at full transparencies, the critical current reduces at B increases, remains finite, and forms a kink at B = Bc, and then develops finite oscillations in the topological phase B > Bc; see the top violet curve in Fig. 41(b). The oscillations above Bc were shown to occur due to the zero-energy splitting at φ = π seen in Fig. 40(c), and are therefore intimately related to the formation of MZMs [161]. Thus, IC traces the gap closing and reveals the emergence of MZMs [161,347], thus offering a way to unambiguously detect them even when the junction hosts trivial ZESs [202,218] and even use them for enhancing the efficiency of Josephson diodes [223]. When the normal transparency of the junction is reduced, the critical currents reduce, with a faster rate in the trivial regime B < Bc. In the tunneling regime, the bottom red curve in Fig. 41(b), the critical current in the trivial regime is vanishingly small, while, interestingly, it exhibits a reentrant behavior at B = Bc and still captures the Majorana oscillations for B > Bc. This intriguing effect has been interpreted to be due to the fact that the contribution to the critical current from MZMs is proportional to |$\sim \sqrt{\sigma _{N}}$|, while from ABS the critical current is proportional to ∼σN, where σN is the normal transmission [158,161]. Thus, for σN ≪ 1, the critical current from MZMs becomes a larger quantity. Another feature of Fig. 41(b) is the change in periodicity of the oscillations for B > Bc, where the period of the oscillations is doubled as the transparency of the junction is reduced. In this case, the tunnel regime promotes the presence of two independent pairs of Majorana state that are not coupled via the N region and oscillate with the same period as a function of B. Again, this period-doubling effect is only attributed to a topological Josephson junction hosting four MZMs and its detection could help distinguish MZMs from trivial ZESs [218]. Signatures of MZMs in supercurrents have also been explored in the nonequilibrium regime, with interesting studies due to multiple Andreev reflections [347,348] and intriguing Shapiro steps [349,350], but experiments, despite the effort [351–354], so far seem inconclusive [355,356]. Before ending this part, we highlight that the Zeeman dependence of the critical currents discussed here has also been reported experimentally in Refs. [175,177], and is compatible with the topological nature of the Josephson junction with MZMs. It is therefore encouraging to further explore phase-biased Josephson transport for the detection of MZMs and topological superconductivity.
![(a) Current–phase curves in the topological phase (B > Bc) of an SNS Josephson junction based on nanowires with SOC for S regions with distinct lengths. (b) Zeeman dependence on the critical currents at LS = 2 μm for different values of the normal transmission, ranging from the tunnel regime (bottom red curve) to the full transparent regime (top violet curve). Parameters: αR = 20 meV nm, Δ = 0.25 meV, μN(S) = 0.5 meV, LN = 20 nm, I0 = eΔ/ℏ. Data adapted from Ref. [161].](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/ptep/2024/8/10.1093_ptep_ptae065/1/m_ptae065fig41.jpeg?Expires=1747856791&Signature=iwCev-4y3AZYdIe7E1KVutc37DurFPNUDGcxroKwvux02W1x8Kgt-MnuVM9yQA7P7shqhgXLTnMJSzTX~M4AOdkAnTwr3kgi34KJObZ0pjgVYUrs11ePH3aOkVYLV6oQMr7CRxmWs2z8kVN~YiI47lMfwHaiz92omP93UH0Jm3Dl2c~60qyPzyNrb54BdfpNdXSJ8ZeqBX~OTejrG5LxtDgvjYkynD~cQUZF7~JlTxbJCtWGJTRddnCGg6P~ROBLhhmUOdXWl~bD33LGHcJBwHp4PRVasADPU8-3XSVikX-Ob6b8wSPhooImQce2529Q51SZtlPQX~nv7GMUBOYQlA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
(a) Current–phase curves in the topological phase (B > Bc) of an SNS Josephson junction based on nanowires with SOC for S regions with distinct lengths. (b) Zeeman dependence on the critical currents at LS = 2 μm for different values of the normal transmission, ranging from the tunnel regime (bottom red curve) to the full transparent regime (top violet curve). Parameters: αR = 20 meV nm, Δ = 0.25 meV, μN(S) = 0.5 meV, LN = 20 nm, I0 = eΔ/ℏ. Data adapted from Ref. [161].
7.3. Josephson current in 2D spin-singlet d-wave superconductor junctions
We here turn our attention to the d.c. Josephson effect in Josephson junctions formed between superconductors with spin-singlet d-wave pair potentials, as schematically shown in Fig. 42. As in previous subsections, here we consider an SIS Josephson junction with the insulating region located at x = 0 and determined by U(x) = Hbδ(x), and the S regions having a d-wave pair potential with a finite phase difference, which, within the quasiclassical approximation, can be expressed as
where |$\Delta _{L(R)}(\theta )=\Delta _{0}(\hat{k}^{2}_{x} - \hat{k}^{2}_{y})$| is a pair potential felt by a quasiparticle where the direction of its motion is determined by the angle θ such that |$\exp (i\theta )=\hat{k}_{x} + i\hat{k}_{y}$|, with |$\hat{k}_{x}=k_{Fx}/k_{F}={\rm cos}(\theta )$|, |$\hat{k}_{y}=k_{Fy}/k_{F}={\rm sin}(\theta )$|, and |$k_{F}=|\boldsymbol {k}_{\rm F}|$|. We consider the situation where the lobe direction of the d-wave pair potential (crystal axis) is tilted from the normal to the interface by an angle α(β) as shown in Fig. 42. Moreover, φL(R) is the superconducting phase of the left (right) superconductor. We are interested in the energies of the ABS and in the supercurrent carried by them. It has been shown that, alternatively to the free energy calculation discussed in the previous subsection, the Josephson current can be calculated using Andreev reflection coefficients as
where ae(θ, iωn, φ) [ah(θ, iωn, φ)] represents the Andreev reflection coefficient of an incident electron (hole) from the left superconductor [343,345], |$\tilde{\Omega }_{n,L,\pm }={\rm sgn} \left(\omega _{n}\right) \sqrt{\mid \Delta _{L}\left(\theta _{\pm }\right) \mid ^{2}+ \omega _{n}^{2}}$|, with ωn being Matsubara frequencies, θ+ = θ, θ− = π − θ, RN is the normal state resistance of the 2D junction, and |$\bar{R}_{N}$| is obtained as
where |$\sigma _{N}(\theta )= 1/(1 + Z_{\theta }^{2})$|, Zθ = Z/cos θ, Z = mHb/(ℏ2kF). In what follows, we set σN(θ) ≡ σN and discuss results for d-wave Josephson junctions. Before doing that, however, it is important to mention that, generally speaking, for the considered d-wave Josephson junctions, it is possible to have distinct types of pair potentials. For instance, for scattering from an electron-like quasiparticle injected from the left superconductor, four types of pair potentials appear, which can be summarized as
where Δ0 is the pair potential amplitude and θ is the angle of incidence with respect to the x-axis, while α and β represent the angles between the lobe directions of the d-wave pair potentials in the left and right superconductors and the normal to the interface of the junction; see Fig. 42 for details.
Sketch of a Josephson junction formed by d-wave superconductors. ELQ and HLQ denote electron-like quasiparticle and hole-like one, respectively. Here, α (β) are angles between the crystal axis of a d-wave superconductor with the normal to the interface, while θ is an incidence angle of an electron-like quasiparticle from the left side of the superconductor.
Now we consider a mirror-type Josephson junction, where α = −β, where the energies of the ABSs can be analytically obtained. In particular, for these types of junctions, the Josephson current has been found to be equal to [343,345]
where φ = φL − φR, |$\bar{R}_{N}$| is given by Eq. (204), and
with ΔL ±(θ) = ΔR∓(θ). To find the ABS, we can search for the zeros of the denominator of Eq. (207) in real energies, iωn → E, whose expression then reads
where ΩL ± are obtained as
Then, by solving Eq. (208) for E, we find the energies of the ABSs given by [41,248]
with ∣Eb(φ)∣ ≤ min(|ΔL +(θ)|, |ΔL −(θ)|). The ABSs energies acquire simpler forms at fixed angles. For instance, for α = 0 and α = π/4, Eb can be simply expressed as [41,249]
From these two expressions, we can already draw some important conclusions. For instance, for α = 0, the obtained phase-dependent ABS Eb(φ) is essentially equivalent to that of the spin-singlet s-wave superconductor junction discussed in Eq. (192). On the other hand, for α = ±π/4, the ABS energy Eb(φ) always exhibits a sin (φ/2) dependence that is independent of the interface transparency and forms a 4π periodicity with φ [249,345]. This 4π periodicity appears in the spin-triplet px-wave superconductor junctions discussed in Eq. (200) and has been understood to be the key for detecting MZMs in the a.c. Josephson current [294]. To understand this discussion, in Figs. 43 and 44, we present the ABS energies as a function of phase Eb(φ) for distinct values of θ, α, and Z. For completeness, in Fig. 45 we also show the ABS energies as a function of θ at fixed φ = 0. For α = π/4, the ZESABS with Eb = 0 appears for −π/2 < θ < π/2 except for the nodal direction with θ = 0, as indicated by the red color in Fig. 45(A). For α = π/8, the ZESABS appears for −π/8 < ∣θ∣ < 3π/8, as depicted by the red color in Fig. 45(B). It is also mentioned that the superconducting diode effect can appear for a d-wave superconductor junction on the surface of the topological insulator where ZESABS plays an important role [357].
Energies of the ABSs in Josephson junctions with d-wave superconductors as a function of φ at fixed θ. The different panels correspond to: (A) θ = 0 and α = 0 and (B) θ = π/4 and α = π/4. Red and blue curves correspond to (a) Z = 3 and (b) Z = 0.
Energies of the ABSs in Josephson junctions with d-wave superconductors as a function of φ at fixed θ. The different panels correspond to: (A) θ = 0 and α = π/8 and (B) θ = π/4 and α = π/8. Red and blue curves correspond to (a) Z = 3 and (b) Z = 0.
Energies Eb of the ABSs in a Josephson junction with d-wave superconductors as a function of θ at Z = 3 and distinct values of φ: (A) φ = 0 and α = π/4, (B) φ = 0 and α = π/8. The solid line is Eb and the dashed curves correspond to ±ΔL +(θ) and ±ΔL −(θ).
Having discussed the ABSs, we now calculate the Josephson current using Eq. (206). It is found that for α = 0 and α = ±π/4, it is possible to sum over Matsubara frequencies analytically in Eq. (207). In the case of α = 0, the phase-dependent supercurrent I(φ) is given by
The temperature dependence of the critical current IC(T) is shown in Fig. 46, which essentially becomes equivalent to the spin-singlet s-wave case found by Ambegaokar–Baratoff [341] and shown in Fig. 39(A).
Temperature dependence of the critical current IC(T) of a mirror-type d-wave Josephson junction where α = −β at Z = 2. The different curves correspond to: (a) α = 0, (b) α = 0.1π, (c) α = 0.25π. Here, Δ0(0) is the magnitude of the pair potential Δ0 at zero temperature.
On the other hand, for α = −β = π/4, the phase-dependent supercurrent RNI(φ) becomes
with Δeff(T, θ) = Δ0sin (2θ). In the limit of |$\sqrt{\sigma _{N}} |\Delta _\text{eff}(T,\theta )|\ll 2k_{B}T$|, the critical current IC(T) is approximated to be
which increases with the decrease of the temperature [343,345,346], as can be seen in curve (c) in Fig. 46. This temperature behavior is also seen in the case of the p-wave Josephson junction discussed in the last subsection.
For sufficiently low temperatures, with |$k_{B}T\ll \sqrt{\sigma _{N}} |\Delta _\text{eff}(T,\theta )|$|, the critical current acquires the following form:
RNIC(T) is proportional to |$\sqrt{R_{N}}$| (the inverse of |$\sqrt{\sigma _{N}}$|), similar to the 1D p-wave superconductor junction discussed in the previous subsection. Furthermore, we also note that the critical currents exhibit a nonmonotonic temperature dependence for α within 0 ≤ α ≤ π/4, which is shown by curve (b) in Fig. 46; see Refs. [248,343,345,346].
In order to understand this exotic behavior of the critical currents discussed in the previous paragraph, we decompose RNI(φ) into negative and positive parts, denoted as Gn(φ) and Gp(φ), respectively, for 0 < φ < π [248,343]. Then, these positive and negative parts are given as
where F(θ, iωn, φ) is given by Eq. (207). Here, the term Gp denotes the contribution to RNI(φ) coming from the domain of θ in π/4 + |α| < |θ| < π/2 and −π/4 + |α| < θ < π/4 − |α|. By taking Eq. (207) into account, it is not difficult to see that, when σN < 1/2, the denominator of F(θ, iωn, φ) becomes small ±π/4 − |α| < θ < ±π/4 + |α| because the signs of |$\tilde{\Omega }_{n,L,+} \tilde{\Omega }_{n,L,-}$| and ΔL +(θ)ΔL −(θ) are different. In particular, in the low transparent limit, when |$\sqrt{\sigma _{N}}|\Delta _\text{eff}(T,\theta )|$|≪ |$\mid$|ωn|$\mid$| ≪ |$|\Delta _\text{eff}(T,\theta )|$|, the negative part Gn(φ) is given as [343]
This expression reveals that the magnitude of Gn(φ) is enhanced with the decrease of temperature [248,343]. On the other hand, Gp(φ) becomes almost constant at low temperatures since the denominator of Gp(φ) is not suppressed with the decrease in temperature due to the absence of ZESABS since ΔL(θ+)ΔL(θ−) ≥ 0 is satisfied. At high temperatures, the phase-dependent supercurrent I(φ) becomes positive for 0 < φ < π because the magnitude of Gp(φ) exceeds that of Gn(φ). On the other hand, at low temperatures, |Gn(φ)| > Gp(φ) is satisfied and the resulting phase-dependent supercurrent I(φ) becomes negative for 0 < φ < π. Another feature to note is that, with the decrease of the temperature, a crossover from 0-to-π-junction occurs [343]. To close, we point out that the nonmonotonic temperature dependence of the Josephson current has been experimentally observed in Josephson junctions of high-Tc cuprates [62,65] and it is qualitatively consistent with curve (b) in Fig. 46.
7.4. Summary
We have discussed the Josephson effect in phase-biased Josephson junctions and have shown that they can be very useful in detecting MZMs. We first showed the ABSs and current–phase curves in Josephson junctions with semi-infinite spin-singlet s-wave superconductors. Then, we addressed the ABSs and current–phase curves in the topological phase of Josephson junctions with semi-infinite p-wave superconductors, showing that the topological ABSs develop a robust zero-energy crossing at φ = π that leads to current–phase curves with a skewed profile that is robust to barrier imperfections. Moreover, in this part, we also showed that the ABSs and supercurrents in finite Josephson junctions based on nanowires with SOC provide crucial information to detect MZMs in realistic systems. Specially, these systems host four MZMs whose spatial nonlocality permits their detection in an unambiguous manner via current–phase curves at φ = π and critical currents, even when topologically trivial ZESs are also present. Lastly, we showed that Josephson junctions made of d-wave superconductors host intriguing ABSs, which can even become completely flat at zero energy and promote large critical currents at low temperatures.
8. Summary and outlook
In this review article, we have summarized some theoretical aspects of MZMs from the viewpoint of SABSs in unconventional superconductors. We started by discussing that MZMs are a special type of ZESABS in spin-polarized p-wave superconductors emerging at zero energy located at the edges and are self-conjugate quasiparticles. We then clarified that the regime where MZMs appear corresponds to a topological phase that can also be realized by using semiconductors with SOC, conventional spin-singlet s-wave superconductivity, and a large Zeeman field. We later showed that, by using the Green’s functions of p-wave superconductors, the LDOS can reveal the emergence of MZMs. Moreover, it also helps to understand that MZMs always accompany the formation of dominant odd-frequency spin-triplet s-wave pair correlations, which represent yet another intriguing but less explored property of MZMs. We also discussed that MZMs produce unique signatures in charge conductance, anomalous proximity effect, and skewed phase-biased Josephson currents. While independently these effects can help to identify MZMs, it is necessary to explore all Majorana properties together, such as their zero-energy nature, their spatial nonlocality, and their odd-frequency pairing.
When exploring conductance, it is necessary to perform simultaneous measurements of ZBCPs at both ends of the system, which should exhibit a strong dependence on the length of the system as a signature of spatial nonlocality. Short systems should exhibit an oscillatory conductance around zero energy, which becomes pinned at zero energy only for long wires and lead to ZBCPs. When trying to detect MZMs in supercurrents, a length-dependent profile of the current–phase curves should be inspected to probe Majorana nonlocality. Along these lines, MZMs and the topological phase can also be revealed by exploiting the Zeeman and system’s length dependences of the critical currents. With these ideas, it might be possible to access the zero-energy nature of MZMs and also their spatial nonlocality. Of course, apart from the predicted Majorana signatures, it is also necessary to test them against effects where it is expected that they do not survive. For instance, the Majorana signatures can be tested in the presence of a Zeeman field parallel to the SOC, which is expected to reduce the induced gap where MZMs appear and thus destroy them. To test the odd-frequency pairing of MZMs, it could be possible to further explore the anomalous proximity effect, perhaps in the same way as the long-range proximity effect in superconductor–ferromagnet junctions [90] because in both cases the spin symmetry of the pair amplitude is spin-triplet. Another less explored signature is the topological phase transition, which is expected to occur in the bulk, and should be carried out along with the measurements of the MZMs at the edges.
While the ideas discussed here might not solve all the difficulties, they could certainly help to overcome some of the problems in relation to the unambiguous identification of MZMs when topologically trivial zero-energy states also appear. Despite the challenges, all the theoretical and experimental studies are clearly advancing our understanding of the phenomena present in Majorana devices.
Conflict of interest statement. None declared.
Acknowledgement
We thank R. Aguado, Y. Asano, O. Awoga, M. Benito, S. Bergeret, A. Black-Schaffer, M. Cuoco, P. Burset, D. Chakraborty, Y. Fukaya, A. Furusaki, P. Gentile, A. A. Golubov, S. Hoshino, S. Ikegaya, T. Kokkeler, S. Kashiwaya, Y. Kawaguchi, S. Kobayashi, B. Lu, S. Matsuo, T. Mizushima, N. Nagaosa, R. Nakai, K. Nomura, S. Nakosai, M. Leijnse, M. Sato, B. Sothmann, P. Stano, S. Suzuki, B. Trauzettel, K. Yada, and A. Yamakage for insightful discussions. Y.T. acknowledges financial support from JSPS with Grants-in-Aid for Scientific Research (KAKENHI Grant Nos. 20H00131, 23K17668, and 24K00583.). S.T. is grateful for the support of the Würzburg–Dresden Cluster of Excellence ct.qmat, EXC2147, project-id 390858490, the DFG (SFB 1170), and the Bavarian Ministry of Economic Affairs, Regional Development and Energy within the High-Tech Agenda Project “Bausteine für das Quanten Computing auf Basis topologischer Materialen”. J.C. acknowledges financial support from the Swedish Research Council (Vetenskapsrådet Grant No. 2021-04121) and the Carl Trygger Foundation (Grant No. 22: 2093).
