Odd-frequency superconducting pairing and subgap density of states at the edge of a two-dimensional topological insulator without magnetism
Jorge Cayao Annica M. Black-Schaffer
cond-mat.mes-hall



Department of Physics and Astronomy, Uppsala University, Box 516, S-751 20 Uppsala, Sweden
 

We investigate the emergence and consequences of odd-frequency spin-triplet s-wave pairing in superconducting hybrid junctions at the edge of a two-dimensional topological insulator without any magnetism.
More specifically, we consider several different normal-superconductor hybrid systems at the topological insulator edge, where spin-singlet s-wave superconducting pairing is proximity-induced from an external conventional superconductor.
We perform fully analytical calculations and show that odd-frequency mixed spin-triplet s-wave pairing arises due to the unique spin-momentum locking in the topological insulator edge state and the naturally non-constant pairing potential profile in hybrid systems. 
Importantly, we establish a one-to-one correspondence between the local density of states (LDOS) at low energies and the odd-frequency spin-triplet pairing in NS, NSN and SNS junctions along the topological insulator edge; at interfaces the enhancement in the LDOS can directly be attributed to the contribution of odd-frequency pairing. Furthermore, in SNS junctions we show that the emergence of the zero-energy LDOS peak at the superconducting phase ϕ=π is associated purely with odd-frequency pairing in the middle of the junction. 



Odd-frequency superconducting pairing and subgap density of states at the edge of a two-dimensional topological insulator without magnetism
    Jorge Cayao and Annica M. Black-Schaffer
    December 30, 2023
===========================================================================================================================================



§ INTRODUCTION


Superconductivity is strongly characterized by the Cooper pair wave function, or the pairing amplitude, obeying the antisymmetry condition imposed by Fermi-Dirac statistics, which restricts its spin and spatial symmetries. The antisymmetry condition also allows for the emergence of odd-frequency superconductivity, where the pairing amplitude is odd in the time, or equivalently frequency, parameter.
The existence of odd-frequency superconductivity was first postulated by Berezinskii<cit.> in 1974, when he pointed out that there is no symmetry restriction for the existence of odd-frequency spin triplet s-wave pairing in ^3He. 
Subsequent works first took this idea to analyze spin-triplet s-wave superconductivity in disordered systems<cit.> and later extended this suggestion to proposed odd-frequency superconductivity in bulk systems with spin-singlet p-wave pairing.<cit.> More recently, bulk odd-frequency superconductivity has also been found in multiband superconductors<cit.> and in superconductors subjected to a time-dependent drive.<cit.>

Following another route for non-bulk materials, it has been shown that odd-frequency superconductivity also emerges as an induced effect in hybrid systems where superconductivity is induced by proximity effect. For ferromagnet-superconductor (FS) junctions it is now well established, both theoretically and experimentally, that it is odd-frequency spin triplet s-wave pairing that explains the long-range proximity effect into the ferromagnet even when using conventional spin-singlet s-wave superconductors.<cit.> In these systems the spin-rotation symmetry is broken by the ferromagnet, which allows for the formation of the spin-triplet state. 
Afterwards, it was also shown that normal metal-superconductor (NS) junctions also exhibit odd-frequency pairing. 
In this case, however, a conventional superconductor can only induce odd-frequency spin-singlet p-wave pairing, 
since only translational symmetry is broken.<cit.> 
However, according to Anderson's theorem<cit.> only s-wave pairing is intrinsically stable against 
ever present non-magnetic disorder, and therefore this latter odd-frequency p-wave is much less stable than the 
odd-frequency s-wave state present in FS junctions. Moreover, normal metals are also fragile to disorder due to 
Anderson localization, namely they allow for finite elastic backscattering processes from non-magnetic impurities, 
causing dissipation of electric current. 
Odd-frequency spin-triplet s-wave pairs can however appear at interfaces of diffusive normal metals if the external superconductor is unconventional with spin-triplet p-wave symmetry.<cit.> 
More recently, the appearance of odd-frequency pairing was also found in systems with 
Rashba spin-orbit coupling.<cit.>

On the other hand, topological insulators,<cit.> characterized by having an insulating bulk but metallic surface states, have strongly suppressed backscattering in their surface states. In fact, for two-dimensional topological insulators (2DTI) backscattering is completely absent due to the perfect helicity of the edge states. Induced s-wave superconductivity in a 2DTI edge therefore represents a promising and very disorder-insensitive platform for the search of robust odd-frequency pairing.
Indeed, it has recently been demonstrated that the combination of surface state helicity and a finite in-surface gradient in a proximity-induced spin-singlet s-wave superconducting state gives rise to odd-frequency spin-triplet s-wave pairing.<cit.> It is the helicity of the TI surface state that allows a spatial symmetry breaking (finite in-surface gradient) to be effectively converted into a change of the spin symmetry, such that spin-triplet s-wave pairing is generated. Notably this mechanism does not need any presence of ferromagnetic regions or even a magnetic field, and thus does not destroy the topological protection of the TI surface states. 

There exits now a growing body of work focusing in detail on behavior of odd-frequency pairing in superconducting hybrid systems at both the surface of 3DTIs<cit.> and edge of 2DTIs.<cit.> However, in these cases finite ferromagnetic regions has been assumed, which destroys the topological protection of the TI surface states. For the situation of odd-frequency pairing without magnetism in 2DTIs there exists no detailed analysis of the pairing amplitudes and especially not of their relationship to the local density of states (LDOS).
In this work we fill this gap and investigate in detail NS, NSN, and SNS hybrid junctions at the edge of a 2DTI without any magnetism present. Our study is based on retarded Green's functions extracted from scattering states and allow us to analytically extract all pairing amplitudes as well as the LDOS in the junctions. 
Very generally, we show that breaking translation symmetry at the NS interface(s) gives rise to four different symmetry classes at 
the NS interfaces: Even-frequency, spin-Singlet, Even in space (ESE), i.e. the conventional order, Odd-frequency, spin-Singlet, Odd in space (OSO), Even-frequency spin-Triplet, Odd in space (ETO), and Odd-frequency, spin-Triplet, Even in space (OTE). In the absence of magnetism in the system we only generate mixed spin-triplet pairing, such that m_z = 0 for the Cooper pairs. Disorder stability is moreover preferential for the local s-wave states present in the ESE and OTE classes and we therefore focus primarily on these. These local pairing terms are also naturally the contributions most directly connected with the LDOS which is also a local quantity. 
A strong relationship between LDOS and  odd-frequency pairing has been previously already established in systems formed out of normal metals with unconventional superconductors<cit.> or when magnetism is present<cit.>. With this work we extend this relationship also to 2DTI superconducting hybrid junctions without any magnetism.


In both NS and NSN junction we find strongly dominating OTE pairing at very low energies in the superconducting interface region, with an exponential decay into the bulk of S. In fact, at zero energy the ESE contribution is completely suppressed and only OTE pairing is present. We also find that the LDOS in the S region experience the same frequency dependence and exponential decay into the bulk of S as the OTE pairing. This shows that the very low energy contribution to the LDOS arises entirely from the OTE pairing. Moreover, we also find that the conductance in NS junctions follows the same behavior as the OTE component at low energies, suggesting that the main contribution to the conductance is also of OTE nature.
Our results here are strongly connected to the fact that the Andreev reflection magnitude at the NS interface reaches its maximum for energies within the gap, a unique characteristic of the helical edges of 2DTIs.<cit.>

In SNS junctions the LDOS reveals the formation of Andreev bound states (ABSs) for energies within the energy gap of 
the superconductor. We find that the ESE and OTE pairing magnitudes exactly capture their emergence, but they behave 
very differently for different values of the superconducting phase ϕ across the junction. Indeed, for very short 
junctions we obtain that at ϕ=0 the ESE pairing dominates over a completely suppressed OTE contribution. 
At ϕ=π the ESE is instead zero while the OTE becomes dominant and even exhibits a zero-energy peak just as 
the LDOS thanks to the topologically protected zero-energy crossing of the ABSs. Moreover, the 
supercurrent across short junctions exhibits its maximum value at ϕ=π, as a result of the 
resonant zero-energy peak which possesses purely OTE pairing.
In longer junctions, the increasing number of ABSs within the gap is also exactly captured in the ESE and OTE amplitudes. In fact, in the middle part of the N region the situation is very similar to the short junction case with only OTE pairing present at ϕ=π. Here the ABSs again have protected crossings generating resonant peaks in the LDOS, which exactly correspond to resonant peaks in the OTE pairing. Thus the LDOS signatures of the ABSs at ϕ = π is entirely a consequence of OTE pairing. 

On the experimental side, HgCd/HgTe<cit.> and InAs/GaSb<cit.> heterostructures represent two of the most promising 2DTI. In both cases, induced superconductivity has already been demonstrated.<cit.>
Moreover, the superconducting junction geometries we consider comprise of a realistic platform for both LDOS and conductance measurements, as has also been demonstrated in experiments. <cit.> Therefore, all experimental prerequisites already exist for the systems we study.
 
This paper is organized as follows. In Sec. <ref> we present the model and the method based on retarded Green's functions calculated from scattering states. In Sec. <ref> we perform a detailed analysis of the pairing amplitudes and investigate their strong relationship with the LDOS for subgap energies for NS, NSN, as well as SNS junctions along the 2DTI edge. We present our conclusions in Sec. <ref>. For completeness, in Appendices <ref>, <ref>, <ref>, and <ref> we provide a detailed account of the method and analytical calculations reported in this work.



§ MODEL AND METHOD


We consider a 2D topological insulator (2DTI), with its 1D metallic edge in proximity to a conventional spin-singlet s-wave superconductor. This system is modeled by the Bogoliubov-de Gennes (BdG) Hamiltonian

    H_BdG   = v_Fp_xτ_zσ_z-μτ_zσ_0+Δ(x)τ_xσ_0 ,

in the basis Ψ(x)=

 (ψ_↑,
 ψ_↓,
 ψ_↓^†,
 -ψ_↑^†
 )^T, where T denotes the transpose operation, ψ^†_σ(x) adds an electron with spin σ=↑,↓ at position x along the edge. 
 The first term represents the 1D metallic edge of a 2DTI,<cit.> where the spin quantization direction is along the z-axis, p_x=-i∂_x, v_F is the velocity of the edge state, σ_i is the ith Pauli spin matrix and τ_i is the ith Pauli matrix in Nambu space. The chemical potential is represented by μ and determines the filling. 
In the last term Δ(x)=Δ(x)  e^iϕ is the induced superconducting pairing potential at the edge. A good interface enables this pairing potential to be induced into the region(s) of the metallic edge in contact with an external superconductor as a result of proximity effect.<cit.> Notice that the pairing potential Δ(x) generally depends on the position at the edge in hybrid junctions. In the normal state (N) the edge states exhibit a linear energy versus momentum dispersion with Δ = 0, while in superconducting (S) regions a finite Δ mixes the electron and hole branches and opens a gap at the Fermi momenta as shown in Fig. <ref>(a,b).

 

We will consider NS, NSN and SNS junctions formed at the edge of the 2DTI, such that the profile of the induced pairing potential is approximated by the step-like functions indicated in green in Fig. <ref>(c). The NS junction captures the physics of the N-S interface, while the NSN and SNS configurations represent finite S regions and Josephson junctions, respectively. 
Here the length, L_ N and L_ S, respectively, of the middle regions become important parameters. Also, in the SNS junction it is possible to attain a finite phase difference between the two S regions, with pairing potentials Δ_ L=Δ and Δ_ R=Δ  e^iϕ, respectively.



 §.§ Retarded Green's function


In this work we primarily investigate the pairing amplitudes and LDOS in different hybrid junctions. 
For this we follow the formalism based on retarded and advanced Green's functions, corresponding to outgoing and incoming waves, respectively. 
In general, the retarded Green's function can be calculated using<cit.> 

    G^r(x,x',ω)=
       α_1Ψ_1(x)Ψ̃_3^T(x')+α_2Ψ_1(x)Ψ̃_4^T(x')    
    
       + α_3Ψ_2(x)Ψ̃_3^T(x')+α_4Ψ_2(x)Ψ̃_4^T(x'),     x>x', 
    β_1Ψ_3(x)Ψ̃_1^T(x')+  β_2Ψ_4(x)Ψ̃_1^T(x')    
    
       +   β_3Ψ_3(x)Ψ̃_2^T(x')+  β_4Ψ_4(x)Ψ̃_2^T(x'),     x<x',

where Ψ_i, i=1,2,3,4 are the wave functions representing the four scattering processes at the NS interface at x = x' indicated by black arrows in Fig. <ref>(b), i.e. right-moving electron, right-moving hole, left-moving quasi-electron, and left-moving quasi-hole. Ψ̃_i corresponds to the conjugated processes obtained using H̃_BdG(k)=H_BdG^*(-k)=H_BdG^T(-k) instead of Eq. (<ref>). The detailed form of these wave functions is provided in the appendix for NS, NSN, and SNS junctions.

The retarded Green's function satisfies the equation of motion

    [ω-H_BdG(x)]G^r(x,x',ω)=δ(x-x') ,

where its integration around x=x',

    [G^r(x>x')]_x=x'-[G^r(x<x')]_x=x'=σ_zτ_z/iv_f ,

provides a system of equations that allows us to find the coefficients α_i and β_i in Eq. (<ref>). Once the retarded Green's functions are found, the advanced Green's function can be calculated using G^a(x,x',ω)=[G^r(x',x,ω)]^†. 

The Green's functions G^r(a) are 4×4 matrices in Nambu space, 

    G^r(x,x',ω)=
    [ G^r_ee(x,x',ω) G^r_eh(x,x',ω); G^r_he(x,x',ω) G^r_hh(x,x',ω) ] ,

where the diagonal elements are the regular electron-electron and hole-hole Green's functions, while the off-diagonal corresponds to the anomalous electron-hole part. Here each element is a 2 × 2 matrix in spin-space.
The electron-electron part allows us to calculate the LDOS

    ρ_(x,ω)=-lim_x→ x'1/π Im Tr[G^r_ee(x,x',ω)],

while the electron-hole part provides the pairing amplitudes. The LDOS is generally one of the most experimentally accessible quantities in a system and can be measured directly using, for example, different scanning tunneling probes.<cit.>
Notice that the LDOS is proportional to the trace of the imaginary local (x=x') Green's function. Therefore, any connection between LDOS and pairing amplitudes is most naturally occurring for local pair amplitudes, which have s-wave symmetry.

Once the retarded and advanced Green's functions are calculated, we can decompose the spin symmetry of the anomalous electron-hole part as 

    G^r(a)_eh(x,x',ω)=f^r(a)_0(x,x',ω)σ_0+∑_j=1^3f^r(a)_j(x,x',ω)σ_j ,

where f_0^r(a) corresponds to spin-singlet, f_1,2^r(a) to equal spin-triplet and f_3^r(a) to mixed spin-triplet retarded (advanced) pairing amplitudes.  

These components all necessarily obey Fermi-Dirac statistics and are therefore antisymmetric under the simultaneous exchange of positions and frequency, 

    f^r_0(x,x',ω)   =f_0^a(x',x,-ω) ,
    
    f^r_j(x,x',ω)   =-f_j^a(x',x,-ω) ,

with j=1,2,3 labeling the spin-triplet pairing amplitudes. Note that when analyzing the symmetry with respect to frequency, the change ω→-ω requires a pass from the retarded Green's function into the advanced function since each are only defined on parts of the time (frequency) axis, see Appendix <ref> for more details.
From this we can isolate the even- and odd-frequency components as

    f_i^r,E(x,x',ω)   =f_i^r(x,x',ω)+f_i^a(x,x',-ω)/2 ,
    
    f_i^r, O(x,x',ω)   =f_i^r(x,x',ω)-f_i^a(x,x',-ω)/2

for all spin components (i = 0,1,2,3).

Taking the full antisymmetry condition of the anomalous Green's functions into account, which includes spin, spatial, and frequency dependence, we are left with a total of four allowed symmetry classes.<cit.> These symmetries are Even-frequency, spin-Singlet, Even in space (ESE), Odd-frequency, spin-Singlet, Odd in space (OSO), Even-frequency, spin-Triplet, Odd in space (ETO), and Odd-frequency, spin-Triplet Even in space (OTE). 
As we will see, NS structures along 2DTI edges generally host all of these four components even though the superconductor itself only has conventional spin-singlet s-wave symmetry (ESE). 
However, considering the stability against disorder and also favoring a direct comparison with the  LDOS which is a fully local property (extracted at x' = x), we will here primarily be concerned with states with local or equivalently uniform s-wave symmetry with f_i(x,x,ω)≡ f_i(x,ω), which are the ESE and OTE components. 



§ PAIRING AMPLITUDES AND LDOS


We now turn to the main objective in this work which is to analyze the pair amplitudes in NS, NSN, and SNS superconducting hybrid junctions at the edge of a 2DTI and its consequences for the LDOS. In normal NS junctions, i.e. not in TIs, it has previously been shown that when a superposition of spin-singlet and singlet-triplet superconductivity is present, the state with larger pairing magnitude dominates the induced superconductivity<cit.> and with remarkable consequences for the LDOS.<cit.> 

Before going into details, we stress that in this work both the retarded and advanced equal spin-triplet pairing amplitudes are zero since we do not have any magnetic order in our hybrid junctions, i.e.  f^r(a)_1,2(x,x',ω)=0. We will therefore refer to the mixed spin-triplet component f_3^r simply as the spin-triplet component. Moreover, we refer to f^r_i as the pairing amplitude, while |f^r_i|=√(f_i^rf_i^r,*) is the pairing magnitude, unless otherwise specified. 






 §.§ NS junction

We first focus on the simplest situation, which consists of a semi-infinite NS junction, where the N region occupies x<0 and S x>0.
Before proceeding further we point out that NS junctions at the metallic edge have previously been found, both theoretically and experimentally, to allow perfect (local) Andreev reflection for energies within the superconducting gap with totally suppressed normal reflection processes (backscattering).<cit.>


Extracting the Green's function from the interface scattering processes we find in the normal region N for x<0 (see Appendix <ref> for details) the even and odd-frequency pairing amplitudes

    f^r, E_0,  N(x,x',ω)   =a_2(ω)/2iv_F C_xx'  e^-i(x+x')/ξ_ω,
    
    f^r, O_0,   N(x,x',ω)   =-a_2(ω)/2v_F S_xx'  e^-i(x+x')/ξ_ω,
    
    f^r, E_3,  N(x,x',ω)   =a_2(ω)/2v_F S_xx'  e^-i(x+x')/ξ_ω,
    
    f^r, O_3, N(x,x',ω)   =-a_2(ω)/2iv_F C_xx'  e^-i(x+x')/ξ_ω,

where C_xx'= cos[k_μ(x-x')],  S_xx'= sin[k_μ(x-x')], and k_μ=μ/v_F, ξ_ω=v_F/ω. Here the superscripts indicate the retarded Green's function (r), the even (E) and odd (O) frequency components, while subscripts stand for spin-component (0 or 3) and normal region (N).
Moreover, a_2 represents the Andreev reflection amplitude for an incident hole from N at the NS interface. At the metallic edge of a 2DTI a_2 it is the same as the Andreev reflection for an incident electron from N  at the NS interface, a_1. For energies within the superconducting gap, a_1,2=  e^-iη(ω), where η(ω)= arccos(ω/Δ).

Despite the non-obvious frequency dependence of  Eqs. (<ref>) it is straightforward to verify their symmetry with respect to frequency. 
For instance, the local ESE amplitude f^r,E_0(x,ω) ≡ f^r,E_0(x,x,ω)=a_2(ω) e^-2ix/ξ_ω/(2iv_F) and thus reversing the sign of ω leads to f^r,E_0(x,-ω)≡ f^a,E_0(x,-ω)=-a^*_2(-ω) e^-2ix/ξ_ω/(2iv_F)=f^r,E_0(x,ω), since a_2^*(-ω) = -a_2(ω). This confirms the even frequency dependence of the ESE pairing amplitude. 

The amplitudes in Eqs. (<ref>) correspond to ESE, OSO, ETO, OTE amplitudes, respectively, and they are all finite for x≠ x' and  finite Andreev reflection amplitude a_2. 
The emergence of all symmetry classes at interfaces of junctions, as the ones discussed here, is not surprising. Indeed, the junction breaks translational symmetry which allows for mixing between different spatial symmetries. Moreover, the 2DTI also allows for mixing between different spin states, without the presence of a magnetic field. Thus all symmetry classes preserving the antisymmetric Fermi-Dirac statistic will appear at the interface.
Moreover, we observe that the even-frequency (odd-freqeuncy) spin-singlet and odd-frequency (even-frequency) spin-triplet pairing amplitudes are equal to each other in magnitude, namely |f^r, E(O)_0(x,x',ω)|=|f^r, O(E)_3(x,x',ω)|.
Considering only s-wave pairing, this means that ESE and OTE contribute equally to the superconducting state in N. This situation with very strong odd-frequency components is due to the unique nature of the 1D metallic edge.
We also note that Eqs. (<ref>) suggest that the pairing amplitudes in N do not decay into N, but rather remain finite as a pure consequence of the NS interface. This has to be understood to happen only very close to the NS interface, where Andreev reflection is manifested.

Performing the same analysis as above for the superconducting region S, x>0, leads us to naturally distinguish between two different terms in the pairing amplitudes, attributed to the bulk (B) and to the interface (I), with the total pairing amplitude f = f_ B + f_ I. For energies within the superconducting gap, |ω|<Δ, we find

    f^r, E_0, B(x,x',ω)   =Z(ω) e^-κ(ω)|x-x'| C_xx' ,
    
    f^r, E_3, B(x,x',ω)   =iZ(ω) e^-κ(ω)|x-x'| S_xx' ,
    
    f^r, O_i, B(x,x',ω)   =0

and 
    f^r, E_0, I(x,x',ω)   =a_3(ω)/2iv_F e^-κ(ω)(x+x')
    B(ω)
    C_xx' ,
    
    f^r, O_0, I(x,x',ω)   =a_3(ω)/2v_F e^-κ(ω)(x+x') S_xx' ,
    
    f^r, E_3, I(x,x',ω)   =a_3(ω)/2v_F e^-κ(ω)(x+x')
    B(ω)
    S_xx' ,
    
    f^r, O_3, I(x,x',ω)   =a_3(ω)/2iv_F e^-κ(ω)(x+x') C_xx' ,

with Z(ω)=(1/iv_F)/[ e^iη(ω)- e^-iη(ω)], B(ω)=-i cot[η(ω)], κ(ω)=√(Δ^2-ω^2), and a_3 = -a_1.
For details on the derivation we refer to Appendix <ref>.
The division into bulk and interface properties is evident for energies within the superconducting gap. Here the interface components exhibit an exponential decay<cit.> proportional to e^-κ(ω)(x+x'), while the bulk components have an exponential decay proportional to e^-κ(ω)|x-x'|. Thus, locally at x=x', the bulk components are independent of spatial coordinates, confirming their bulk nature, while the interface components display an exponential decay into the bulk of the S region. Notice also that the interface amplitudes are all proportional to the Andreev reflection amplitude a_3 = -a_1 and thus purely a consequence of the NS interface.
Analyzing Eqs. (<ref>-<ref>) we find that only the even-frequency ESE and ETO components are finite in the bulk. At the interface, however, all symmetry classes (ESE, OSO, ETO, OTE) are present. 
 Even limiting the interest to the disorder robust s-wave states leaves both ESE and OTE amplitudes in the interface region. 
 Note that an OTE component does not appear in normal metal NS junctions, but there instead a magnetic field is necessary to generate the OTE component.


In Fig. <ref> we display the frequency dependence of the pairing amplitudes at the interface,  x=x'=0. This confirms the odd frequency dependence of the OTE component, while the ESE is fully even in frequency. While the figure only shows the x=x'=0, we have verified the frequency dependence also for x ≠ x'.
Notice that the real part of OTE exhibits an abrupt transition across zero, which arises due to the discontinuity between f^a and f^r at zero frequency.
Thus, the magnitude of the OTE component is finite even for ω=0, as is shown in Fig. <ref>(b). One could argue that at ω=0 the retarded and advanced Green's functions are not well defined and therefore this discussion might lead to a wrong interpretation of the emergence of OTE. We can rule out such a possibility by pointing out that away from zero frequency, but well below Δ, the OTE component is still finite and even clearly dominant over the ESE component. 
In fact, we find that the OTE amplitude dominates over ESE at the interface for |ω|<Δ/√(2).



In order to establish a connection with the LDOS, we compare the LDOS with the local, s-wave, pairing amplitudes at x=x'.
In Fig. <ref> we present the spatial dependence of the LDOS and the local pairing magnitudes in the normal and interface regions with the interface placed at x=0.
In the normal region, x<0, the LDOS acquires a finite value ρ_N(x,ω) =1/(π v_F) independent of both position and energy, as shown in Fig. <ref>(a). This naturally arises because, at the metallic edge of a 2DTI, an incident electron is purely reflected back as a hole without interference.<cit.>
The equal and finite values of ESE f^r, E_0,  N(x,ω) and OTE f^r, O_0,  N(x,ω) in N, first and last equations in Eqs. (<ref>), suggest that ρ_N arises due to their simultaneous contribution. 

In the superconducting region, x>0, the situation is more interesting.  Within the superconducting gap, |ω|<Δ, the LDOS is 

    ρ_ S(x,|ω|<Δ)=1/π v_F e^-2κ(ω)x .

At the interface, x=0, the LDOS is thus finite for energies within the superconducting gap Δ, acquiring its maximum value and then decaying into the bulk of S with a  decay length given by 1/[2κ(ω)], and also seen in Fig. <ref>(a). 
For |ω|>Δ the inset in Fig. <ref>(a) shows the energy dependence and spatial dependence in the S region. Well within the S region the LDOS is depleted below the gap and shows an oscillatory pattern for |ω|>Δ, which arises due to interference between the incident and reflected quasiparticles,<cit.> unlike in the N region. 

For the local pairing amplitudes in S we simplify Eqs. (<ref>) by setting x = x' and obtain

    f^r, E_0, I(x,ω)   =a_3(ω)/2iv_F e^-2κ(ω)x
    B(ω) ,
    
    f^r, O_3, I(x,ω)   =a_3(ω)/2iv_F e^-2κ(ω)x ,

for the ESE and OTE interface amplitudes. Notice that the ESE  interface component is multiplied by the factor B(ω), which leads to important consequences at low energies. 
Remarkably, at the interface and zero energy the ESE interface component is totally suppressed because B(ω=0)=0, while 
the OTE amplitude is still non-zero and even acquires its maximum value, as seen in Fig. <ref>(b). We can therefore directly associate the enhancement of the LDOS close to the interface in S, Fig. <ref>(a), with the OTE pairing and not with the conventional ESE pairing.[A similar conclusion was achieved in junctions between normal metals and unconventional superconductors,<cit.> which are sensitive to backscattering and disorder unlike our 2DTI system.] Interestingly, this behavior is preserved for energies away from (but close to) zero. For energies close to Δ, however, there is an equal contribution of the ESE component.
We have verified that these results do not change by introducing an insulating barrier at the NS interface, since the Andreev reflection coefficient a_1 remains unchanged. Thus the 1D metallic edges of 2DTI offers more robust predictions unlike its 2D surface counterparts, which are sensitive to an interface barrier.<cit.> We have also only considered a step-like pairing potential Δ. Preliminary results (not shown) using a smoother Δ(x) across the NS interfaces, present due to inverse proximity effect, show that there is a reduction of the pairing magnitudes at the NS interfaces, characterized by the sharpness of the Δ(x) profile. However, the classification of the pairing symmetries as well as the main conclusions remain unchanged. This is consistent with the OTE pairing amplitude being generated by an in-surface gradient,<cit.> which naturally is maximized for step-edge profiles.
 
  
Based on our results we can also relate the conductance across the NS junction to the OTE pairing.
The zero-temperature single mode conductance for an incident electron from N is given by G_ NS(ω)=(e^2/h)(1+|a_1|^2-|b_1|^2), where a_1 and b_1 are the Andreev and normal reflection amplitudes, respectively. Since backscattering is forbidden in the 2DTI edge b_1=0. Therefore, the conductance is fully described  by the Andreev reflection coefficient a_1. At the same time, the interface ESE and OTE pairing amplitudes are both proportional to the Andreev reflection coefficient a_3 = -a_1. At low energies the B(ω) factor, however, dramatically suppresses the ESE pairing and left is therefore only the OTE contribution.
Thus, any change in the Andreev reflection process is directly manifested in both the conductance across the NS junction and OTE pairing magnitude at low energies.
This allows us to conclude that the main contribution to G_ NS is given by the OTE component only. Notice that the conductance across NS junctions at the edge of a 2DTI is independent of the chemical potentials in the two regions, unlike the situation with 3DTIs,<cit.> and also becomes constant within the gap while exhibiting a decay outside,<cit.>
as is evident from the behavior of a_1 in Figs. <ref>(b) and also <ref>(a). 

To conclude this part we especially stress that the contribution of the OTE pairing to the enhancement of LDOS and conductance across the NS junction is a direct consequence of the metallic edge of 2DTI, where the Andreev reflection is perfect for energies below the gap<cit.> and also independent of barrier imperfections at the NS interface.<cit.>
Moreover, we point out that Andreev reflection at the edge of a 2DTI has already been experimentally demonstrated<cit.> and significant effort has also been devoted to investigate induced superconductivity at the edge of 2DTIs. <cit.> Our findings can therefore  help to elucidate the nature of the superconducting pairing in both LDOS and conductance measurements in NS junctions.








 §.§ NSN junction

Having investigated the simplest NS junction we now turn to the situation with a finite S region. 
We consider a NSN junction at the edge of a 2DTI, where the S region is placed within 0<x<L_ S, surrounded by N regions on each side.
We arrive at the pairing amplitudes and the LDOS following the same procedure as above, constructing the retarded Green's function from the allowed scattering processes, see Appendix <ref> for details.

In the two normal regions we find the pairing amplitudes proportional to the Andreev reflection amplitude a_2(ω,L_ S)= sin[iκ(ω)L_ S]/ sin[iκ(ω)L_ S-η(ω)] but otherwise exactly following the same form as for the pairing amplitudes in the N region of NS junctions given in Eqs. (<ref>). Thus the only change for the N pairing amplitudes compared to the NS junction is a_2(ω)→ a_2(ω,L_ S). This is fully consistent with the results in N being purely due to the interfaces. 
In particular, notice that the dependence on L_ S  is only present through a_2(ω,L_ S). The magnitude |a_2(ω,L_ S)| as a function of ω for different L_ S is plotted in Fig. <ref>(a). For |ω|<Δ the magnitude increases as the length of the S region increases and saturates to the value found in NS junctions |a_2(ω)|=|a_2(ω,L_ S→∞)| for L_ S> ξ, where ξ=v_ F/Δ is the superconducting coherence length. For energies above Δ, |a_2(ω,L_ S)| develops an oscillatory decay whose amplitude for ξ<L_ S≪∞, is higher than that in NS junctions. 

Turning to the middle superconducting region and concentrating on energies below the gap Δ, we obtain the following even and odd-frequency pairing amplitudes divided into bulk contributions

    f_0, B^r, E(x,x',ω)   =2C_xx'[β_2 e^-κ(ω)|x-x'|+β_3 e^κ(ω)|x-x'|] ,
    
    f_3,  B^r, E(x,x',ω)   =2iS_xx'[β_2 e^-κ(ω)|x-x'|+β_3 e^κ(ω)|x-x'|] ,
    
    f_0,3, B^r, O(x,x',ω)   =0

and interface contributions

    f_0, I^r, E(x,x',ω)   =C_xx'[β_42^- e^-κ(ω)(x+x')+β_13^- e^κ(ω)(x+x')] ,
    
    f_0, I^r, O(x,x',ω)   =-iS_xx'[β_42^+ e^-κ(ω)(x+x')+β_13^+ e^κ(ω)(x+x')] ,
    
    f_3, I^r, E(x,x',ω)   =iS_xx'[β_42^- e^-κ(ω)(x+x')+β_13^- e^κ(ω)(x+x')] ,
    
    f_3, I^r, O(x,x',ω)   =-C_xx'[β_42^+ e^-κ(ω)(x+x')+β_13^+ e^κ(ω)(x+x')] ,

where β_42(13)^±=β_4(1)±β_2(3), β_4,3=- e^-2iη(ω)β_2,1, β_1=-iZ(ω) e^iη(ω)-κ(ω)L_ S/4 sin[iκ(ω)L_ S-η(ω)] and β_2=-β_1 e^2κ(ω)L_ S.
First of all, we stress that these results are in agreement with the expressions for the NS junction. Indeed, we find that in the bulk only even-frequency pairing amplitudes exist, while all pairing symmetries have finite contributions at the interface. 
Observe also that the interface components develop an exponential decay from both interfaces proportional to e^±κ(ω)(x+x'), while the bulk components are proportional to e^±κ(ω)|x-x'| and thus become independent of space for x=x'. 

Next we proceed to analyze how the local pairing amplitudes for x=x' are related to the LDOS.
The regular Green's functions in the N region are exactly the same as in the NS junction case and the LDOS in this N region is subsequently also the same, ρ_  N = 1/(π v_F). On the other hand, in the S region we find a much more complex expression for the LDOS (see Appendix <ref> for details). However, it acquires a simple expression at zero energy
 
    ρ_ S(x,ω=0)=1/π v_F e^-2L_ S/ξ e^2x/ξ+ e^-2x/ξ/1+ e^-2L_ S/ξ .

In Fig. <ref>(a) we plot the spatial and frequency dependence of the LDOS. At the interfaces, x=0, 10ξ, the LDOS is finite and exhibits its maximum value 1/(π v_F). It then exponentially decays from both interfaces into the S region with a decay length for ω=0 given by ξ/2. This behavior is preserved even for energies away from zero as can be observed in Fig. <ref>(a). For ω≠0 the decay length is given by 1/[2κ(ω)].
 

We compare the LDOS with the local pairing amplitudes at x = x', derived from Eqs. (<ref>) and (<ref>), giving the finite components

    f_0, B^r, E(x,ω)   =β_2-β_1 e^-2iη(ω) ,
    
    f_0, I^r, E(x,ω)   =
    B_1(ω)
    [β_1 e^2κ(ω)x -β_2 e^-2κ(ω)x] ,
    
    f_3, I^r, O(x,ω)   =B_2(ω)
    [β_1 e^2κ(ω)x +β_2 e^-2κ(ω)x] ,

which correspond to the ESE class in the bulk and the ESE and OTE classes at the interface, where B_1,2(ω)=±(1± e^-2iη(ω))/2.
Figure <ref>(b) shows the spatial dependence of the interface and normal region pairing magnitudes for different values of ω.
The pairing amplitudes decay from the interfaces at x=0 and x=L_ S into the bulk of S with a decay length given by 1/[2κ(ω)]. Remarkably, this spatial decay is also observed in the LDOS in Fig. <ref>(a), directly supporting their relationship.
Interestingly, at zero energy the coefficients of the interface pairing amplitudes acquire different prefactors, B_1(ω=0)=0 and B_2(ω=0)=2, which gives rise to purely OTE pairing at the interface. A strongly dominating OTE component is preserved even for energies away from ω=0 but well below Δ as clearly seen in Fig. <ref>(b). 


The previous discussion can be clarified by analyzing the frequency dependence of |f_ I^r|, plotted in Fig. <ref>(b-d) for different positions in S. 
At the left interface, x=0, Fig. <ref>(b) shows the OTE pairing being finite and remaining constant as ω increases within the gap. On the other hand, the ESE pairing is zero at ω=0 and increases only following a V-dependence  as ω increases. 
To derive the condition for the parameter space where OTE pairing dominates over the ESE pairing, we solve |f_0, I^r, E(ω,x)|=|f_3, I^r, O(ω,x)| for ω. As for NS junctions, we find a simple expression at x=0, where the OTE component is larger for |ω|<Δ/√(2).
We also point out that the Andreev reflection amplitude a_1(ω,L_ S) = a_2(ω,L_ S) is directly connected to the pairing amplitudes. To visualize this statement we compare the Andreev reflection and pairing magnitudes at the the left interface x=0 in Figs. <ref>(a) and (b). For energies within Δ the Andreev reflection magnitude is finite and constant for finite L_ S≳ξ. In these junctions it is only the OTE amplitude that is finite and constant, while the ESE amplitude is zero at zero energy, increases linearly with ω, and is subdominant all the way until ω = Δ/√(2). 
As we go inside the S region, at x=ξ, the OTE pairing magnitude is reduced and the ESE pairing enhanced, both taking a U-shape form, see Fig. <ref>(c). Finally in the middle of the S region  at x = 5ξ, the OTE magnitude is totally suppressed, while a BCS-like gap is fully developed for the ESE component, see Fig. <ref>(d). This is also what we expect due to the ESE nature of the induced superconducting pairing in the bulk. Notice that in the inset of Fig. <ref>(a) we see how the LDOS at low energies is finite at the interface, x=0, exhibits the same U-shape observed in the pairing magnitudes at x=ξ, while deep in the S bulk at x=5ξ the expected BCS-like induced gap is totally developed. 
Thus the LDOS directly follows the magnitude of the OTE  component |f_3, I^r, O| and we conclude that the enhancement of the LDOS in S close to the interfaces arises mainly due to the OTE contribution. This relationship is exact at ω=0, as the pairing there is purely OTE. For 0<ω<Δ/√(2), when the ESE pairing magnitude is finite, the OTE magnitude is still dominating close to the interface and thus has the largest influence on the LDOS.
This confirms that the emergence of the exotic odd-frequency superconductivity is intrinsically contained in the LDOS at the interface of NSN junctions at the edge of 2DTI, which supports our findings for NS junctions.







 §.§ SNS junction

We finally turn to a SNS junction along the edge of a 2DTI, where the finite normal region with Δ=0 is restricted to 0<x<L_ N. We further set the order parameters in  the left and right S region to Δ_ L=Δ and Δ_ R=Δ  e^iϕ, respectively. We arrive at the expressions for the LDOS and all pairing amplitudes following the same formalism as used for NS and NSN junctions, see Appendix <ref> for a detailed derivation.

As for NS and NSN junctions, all pairing amplitude symmetry classes are also present in SNS junctions. Here we however also obtain a dependence on the superconducting phase across the junction. As we will see, this phase dependence gives rise to a very strong connection between the LDOS and the local s-wave pairing amplitudes.

Let us first focus on the normal N region. We find the even and odd-frequency components of the pairing amplitudes given by

    f_0, N^r, E(x,x',ω)   =W_+(ω) cos[k_μ(x-x')] ,
    
    f_0, N^r, O(x,x',ω)   =W_-(ω)i sin[k_μ(x-x')] ,
    
    f_3, N^r, E(x,x',ω)   =W_+(ω)i sin[k_μ(x-x')] ,
    
    f_3, N^r, O(x,x',ω)   =W_-(ω) cos[k_μ(x-x')]

where W_±(ω)=[m_5 e^i(x+x')/ξ_ω± m_6 e^-i(x+x')/ξ_ω ]/2. 
The coefficients m_5=p_1r̃_4α_1 and 
m_6=q_2s̃_3α_4 contain important information about the scattering processes across 
the junction. 
In fact, p_1 is the amplitude for electron transmission from left S into N, 
r̃_4 gives the Andreev reflection at the left SN interface of a left-moving electron into a hole, 
q_2 gives the Andreev reflection at the right SN interface of a right-moving hole into an electron, and
s̃_3 is the hole transmission from right S into N.
Moreover, we find 
m_6(ω,L_ N,ϕ)=m_5(ω,L_ N,-ϕ) e^i(2L_ N/ξ_ω+ϕ) 
and for energies within the gap, |ω|<Δ,

    m_5(ω,L_ N,ϕ)   =-1/2v_F e^i(ϕ/2- L_ N/ξ_ω)/ sin[η(ω)-L_ N/ξ_ω+ϕ/2] .

The first general observations are that the pairing amplitudes, given by Eqs. (<ref>), become dependent on the superconducting phase difference ϕ through the coefficients m_5,6. Note also that all symmetry classes have finite amplitudes. By simple inspection we notice that they do not exhibit any decay but, for frequencies away from zero, rather develop an oscillatory behavior which extends over the whole N region. 

For the two S regions we obtain similar expressions for the pairing amplitudes as for S in the NS junction case. In fact, they include two components that represent the contributions from the interface, proportional to e^±κ(ω)(x+x'), 
and two components associated with the bulk and thus proportional to e^-κ(ω)|x-x'|. Here the ± signs indicate that the exponential decay occurs from both interfaces and into the bulk of the left and right S regions, respectively. Exact expressions are given in Appendix <ref>, but they are not needed here for the future discussion.

Let us now focus on the local pairing amplitudes ESE and OTE present at x = x' and their 
relation to the LDOS. We have verified that the results found in N and in the left and right S regions
provide the same information and thus for clarity we only discuss the results in the N region. 
The N region components, given by Eqs. (<ref>), for energies within the gap |ω|<Δ at x =x' reduce to

    f_0, N^r, E(x,ω)   =1/2[m_5 e^2i x/ξ_ω+m_6 e^-2ix/ξ_ω] ,
    
    f_3, N^r, O(x,ω)   =1/2[m_5 e^2ix/ξ_ω-m_6 e^-2ix/ξ_ω] .

These expressions for the pairing amplitudes can be shown to be directly connected to the LDOS and especially its phase dependence.
By analyzing m_5,6 we notice that the zeros of their denominators appear at

    2η(ω)-2L_ N/ξ_ω±ϕ=2π n .
 
But this equation is also the exact quantization condition of the Andreev bound states (ABSs) in SNS junctions at the edge of a 2DTI with a normal region of length L_ N,<cit.> see also Appendix <ref>.  These ABSs naturally give very strong peaks in the subgap LDOS in the N region, but here we, quite remarkably, instead obtained them from the pairing amplitudes given by Eqs. (<ref>). This immediately establish a direct relationship between local pairing amplitudes and the LDOS.

A short SNS junction, L_ N≪ξ, only hosts one pair of ABSs at energies ω_±(ϕ)=±Δ cos(ϕ/2), which are 4π-periodic and develop zero-energy crossings at ϕ=π(2n-1), n=1,2,… protected by time-reversal symmetry.<cit.> 
These zero-energy protected crossings give rise to a bound state with Majorana-like properties.<cit.> In a long SNS junction, L_ N≫ξ, on the other hand, many more ABSs fit within the superconducting energy gap. Here the ABSs appear at  ω_±,n(ϕ)=(v_F/2L)[2π(n+1/2)±ϕ], which gives a linear dispersion with respect to the superconducting phase difference ϕ.



  §.§.§ Short junctions

To further analyze the connection between pairing amplitudes and ABSs energies we first investigate the situation of a short junction, L_ N≪ξ, corresponding to the limit L_ N→ 0. 
From Eq. (<ref>) we obtain  m_5(ω,L_ N=0,ϕ)=-(1/2v_ F) e^iϕ/2/ sin[η(ω)+ϕ/2] and m_6(ω,L_ N=0,ϕ)=m_5(ω,L_ N=0,-ϕ) e^iϕ. By simple inspection we notice that this simplified dependence on ϕ has profound consequences for the pairing amplitudes in Eqs. (<ref>). 

At ϕ=2π n, for n=0,1,2,…, the coefficients are the same, i.e., m_5(ω,L_ N=0,2π n)=m_6(ω,L_ N=0,2π n), and therefore

    f_0, N^r, E(x=0,ω)   =-1/2v_F1/ sin[η(ω)] ,
    
    f_3, N^r, O(x=0,ω)   =0 .
 
Thus only ESE pairing survives, while the OTE term is completely vanishing. 
This can be understood from the meaning of coefficients m_5,6. Although they correspond to different (electron and hole) scattering processes, at zero phase difference, they are equal and destructively interfere, which gives rise to zero OTE  but finite ESE pairing.
To visualize these results, we plot in Fig. <ref>(a) the ESE and OTE pairing magnitudes as a function of ω at ϕ=0. Observe that the OTE magnitude is completely suppressed, while the ESE contribution exhibits resonant peaks at the gap edges ±Δ. These ESE peaks are in full agreement with the ABSs energies, ω_±(2π n)=±Δ for  n=0,1,2,…, which at ϕ = 0 merges with the continuum (indicated by blue arrows).


At ϕ=π(2n-1) for n=1,2,…, the coefficients instead acquire opposite values, m_5(ω,L_ N=0,π(2n-1))=-m_6(ω,L_ N=0,π(2n-1)), which again affects the pairing amplitudes in a crucial way

    f_0, N^r, E(x=0,ω)   =0 ,
    
    f_0, N^r, O(x=0,ω)   =-i/2v_F1/ cos[η(ω)] .

Here the ESE contribution is instead completely suppressed, while the OTE term acquires a finite value, as is also shown in Fig. <ref>(b).  
At ϕ=π(2n-1) the ABSs cross zero energy, ω_±(π(2n-1))=0 for n=1,2,…, developing a protected crossing and also introducing a resonant peak at ω=0  in the LDOS.
Other work has found similar phase dependence for a combination of spin-orbit coupling and magnetism.<cit.> Very importantly, in the 2DTI we obtain that no magnetism is necessary.
This LDOS peak (indicated by red arrow) at zero energy is directly reflected in the OTE pairing amplitude also peaking at zero energy.[A zero energy peak in the odd-frequency pairing was also reported in junctions with unconventional superconductors,<cit.> where the pairing function is peaked due to the sign change of the pairing potential. Moreover, a zero energy peak in the OTE amplitude was also predicted in Rashba double wires,<cit.> where Rashba coupling and the additional wire index allows for OTE pairing.]
The mechanism is here the same as mentioned before: two scattering processes, Andreev reflection and normal transmission for electrons and holes, become exactly opposite at ϕ=π, unlike  at ϕ=0 where they are the same. Therefore it exists a clear interference pattern with distinguishable  signatures for each process. We stress that both results, the total suppression of OTE and ESE at ϕ=0 and ϕ=π, respectively, arise due to the metallic nature of the 2DTI edge. In normal metals finite backscattering processes usually induce finite normal reflection amplitudes.

Before concluding we would like to also point out the connection between OTE and supercurrents I(ϕ). In both short and long junctions the supercurrent  can be calculated from the energy levels ω_n  as<cit.> I(ϕ)=-(2e/ħ) ∑_n>0dω_n(ϕ)/dϕ. In short junctions only the pair of ABSs within the gap contribute to I(ϕ). Thus, for short junctions with the ABSs given by ω_±(ϕ)=±Δ cos(ϕ/2), the individual current contribution becomes I_±(ϕ)= ± I_ c sin(ϕ/2), where I_ c=eΔ/ħ is the maximum supercurrent across the junction. Provided fermion parity conservation I(ϕ)  notably exhibits the well-known 4π-periodic fractional Josephson effect<cit.> in ϕ. Notice that I_±(ϕ) is  zero at ϕ=0 while at ϕ=π it acquires  its maximum value. The former arises because there are no ABSs in the junction, while in the latter case the ABSs develop a protected crossing at zero-energy, giving rise to the resonant peak in the LDOS. This resonant peak is, as shown above, purely due to OTE pairing. Thus, the maximum supercurrent, i.e. the critical current I_ c, exhibits contributions from purely OTE pairing, since at ϕ=π only OTE pairing is present.

We can thus conclude that the pairing amplitudes entirely capture the LDOS through the emergence of ABSs, and, remarkably, the zero energy peak in the LDOS is purely a consequence of the OTE pairing as shown in Fig. <ref>(b). These findings are in agreement with the odd-frequency nature previously found for a single Majorana state in time-reversal symmetry breaking systems.<cit.> However, note that in our setup time-reversal symmetry is preserved and the Majorana mode is therefore necessarily two-fold degenerate. 




  §.§.§ Long junctions

Next we proceed by increasing the normal region length L_ N. Then the number of energy levels, bound in the junction within the energy gap Δ, increases and is proportional to L_ N/ξ. In this situation the ESE and OTE pairing amplitudes, given by Eqs. (<ref>), strongly depends on x. In what follows we primarily only need to consider two  situations of this spatial dependence, at the left interface, x=0, and in the middle of N, x=L_ N/2.

At the left interface, x=0, the pairing amplitudes for ϕ=0 read 

    f_0, N^r, E(x=0,ω)   =-1/2v_F cos(ω L_ N)/ sin[η(ω)-L_ Nω] ,
    
    f_3, N^r, O(x=0,ω)   =i/2v_F sin(ω L_ N)/ sin[η(ω)-L_ Nω]

and for ϕ=π

    f_0, N^r, E(x=0,ω)   =-1/2v_F sin(ω L_ N)/ cos[η(ω)-L_ Nω] ,
    
    f_3, N^r, O(x=0,ω)   =1/2iv_F cos(ω L_ N)/ cos[η(ω)-L_ Nω] .




The magnitude of these ESE and OTE pairing terms is plotted in Fig. <ref>(a,b) at ϕ=0 and ϕ=π, respectively. 
At ϕ=0 the ESE and OTE magnitudes reach zero value when ω L_ N=(2n+1)π/2 and ω L_ N=π n for n=0,1,2,…, respectively. In fact, the zeros of ESE correspond to the maxima of OTE and vice-versa.
Specifically, for ω=0 the OTE pairing magnitude is vanishing, while the ESE magnitude remains finite and proportional to 1/(2v_F). 
On the other hand, at ϕ=π the position of the zeros of the ESE and OTE magnitudes is reversed.
 Remarkably, the OTE magnitude now has a resonant peak at zero energy, while the ESE magnitude remains finite at, again, approximately 1/(2v_F). 
Thus, although the OTE pairing magnitude is dominant around zero energy, the finite value of ESE obscures the pure OTE contribution to the LDOS we found in the short junction limit. 

The coexistence of ESE and OTE pairing in long junctions at the interface x=0 can be avoided by probing the pairing amplitudes in the middle of the N region at x=L/2, as presented in Fig. <ref>(c,d). In this case, Eqs. (<ref>) reduce to

    f_0, N^r, E(x=L_ N/2,ω)   =
    m̅_5+m̅_6 ,
    
    f_3, N^r, O(x=L_ N/2,ω)   =
    m̅_5-m̅_6 ,

where m̅_5(ω,L_ N,ϕ)=(-1/2v_F) e^iϕ/2/ sin[η(ω)-L_ N/ξ_ω+ϕ/2] and m̅_6(ω,L_ N,ϕ)=m̅_5(ω,L_ N,-ϕ) e^iϕ. Thus, in the middle of the N region the coefficients m̅_5,6 lose their length dependence in the numerator, unlike m_5,6. Moreover, Eqs. (<ref>) very much resemble the behavior of the pairing amplitudes in short junctions discussed before.
From Eqs. (<ref>) it is simple to conclude that at ϕ=0, m̅_5(ω,L_ N,0)=m̅_6(ω,L_ N,0) and thus

    f_0, N^r, E(x=L_ N/2,ω)   =
    -1/2v_F1/ sin[η(ω)-L_ N/ξ_ω] ,
    
    f_3, N^r, O(x=L_ N/2,ω)   =0 .

Thus the OTE pairing is now completely suppressed for all energies within the energy gap. The ESE component, on the other hand, is still finite with multiple resonant peaks within the energy gap. These resonant peaks in the ESE magnitude directly correspond to the energies of the ABSs ω_±,n(ϕ=0)=(v_F/2L)[2π(n+1/2)] bound within Δ, as indicated by blue arrows in Fig. <ref>(c). 
Remarkably, the situation at ϕ=π implies m̅_5(ω,L_ N,π)=-m̅_5(ω,L_ N,π), leading to instead

    f_0, N^r, E(x=L_ N/2,ω)   =
    0 ,
    
    f_3, N^r, O(x=L_ N/2,ω)   =1/2iv_F1/ cos[η(ω)-L_ N/ξ_ω] .

Thus, at ϕ=π and in the middle of the N region, the OTE component is the only local pairing amplitude, since ESE is completely vanishing. Furthermore, the OTE resonant peaks reveal the emergence of the ABSs with energies ω_±,n(ϕ=π)=(v_F/2L)[2π(n+1/2)±π], as indicated by red arrows in Fig. <ref>(d). 
Thus the peaks appearing in the subgap LDOS from the ABSs in the middle of the N region  can be directly associated with different pairing amplitudes. At ϕ=0 the LDOS peaks are associated with pure ESE pairing, while at ϕ = π they are solely emerging due to the OTE pairing. Also, importantly, the zero energy bias peak at ϕ signaling the presence of Majorana fermions in the SNS junction is entirely associated with OTE pairing. This extends the results for a single pair of ABS in the short junction regime to include all subgap ABSs in longer junctions. 

Above we found that for both short and long junctions it is only the OTE pairing that is intimately connected to the energy resonance at zero energy.
However, the dominance of the OTE pairing is not restricted to zero energy, ϕ=π phase, or the middle of the junction region. This can clearly be appreciated in Fig. <ref>, where we plot the spatial dependence of the ESE and OTE pairing magnitudes over the whole N region for energies away from zero and phase π. 

Although at the interface x=0, L_N, both ESE and OTE pairing magnitudes usually coexist at all phases, it is in the middle of the N region at x=L_ N/2, where OTE dominates over the completely suppressed ESE component for ϕ=π.
The dominating OTE pairing extends over a notably finite region in the middle of the junction. It also survives to quite large subgap energies, ω≳ 0.4Δ, and also for phase values.  Despite ESE often exhibiting a finite value in these regimes, the dominant behavior of OTE over ESE is clearly observed. Also notice how the spatial and energy dependences give rise to an oscillatory behavior of the pairing magnitudes given by Eqs. (<ref>) as can be seen in Fig. <ref>. 
Overall this allows us to conclude that OTE pairing and its influence on the LDOS is not restricted neither to zero energy nor ϕ=π.







§ CONCLUSIONS


In this work we have analytically studied the emergence of odd-frequency superconductivity in NS, NSN, and SNS junctions at the edge of a 2DTI without magnetism or any other time-reversal symmetry breaking perturbations, where the S regions have spin-singlet s-wave superconducting pairing induced by proximity to an external conventional superconductor.
We have shown that odd-frequency mixed spin-triplet s-wave superconductivity does not require the presence of magnetic order but naturally arises at any NS interface as a result of breaking translation invariance in combination with the helicity of the 2DTI edge state. Moreover, we have clearly extended previous studies<cit.> and  established a one-to-one correspondence between the subgap LDOS and odd-frequency pairing in NS, NSN, and SNS junctions at the edge of 2DTI without any magnetism. These geometries are all suitable for LDOS as well as conductance measurements as has recently been demonstrated in experiments. <cit.>

For NS and NSN junctions, we have derived analytical expressions for both the pairing amplitudes and LDOS. Focusing on local, s-wave, pairing we have shown that both even-frequency, spin-singlet, even-parity (ESE) and odd-frequency, spin-triplet, even-parity (OTE) pairing is generally present in these systems. 
In the normal N region(s) the ESE and OTE magnitudes are equal and proportional to the Andreev reflection magnitude. 
In the superconducting S region the pairing amplitudes have two components arising from either the bulk or the interface(s). In the bulk of the S region, only the ESE pairing magnitude survives. However, there are both ESE and OTE interface components that develop an exponential decay into the bulk of the S region with the decay length set by the superconducting coherence length for small energies. 
Very close to the interfaces in the S region the OTE pairing becomes very dominant over an extremely suppressed ESE amplitude for energies well below the superconducting gap.
Moreover, the behavior of the low-energy LDOS in the S region close to the interface has the same dependence on the 
energy and distance from the interface as that of the OTE pairing. This allows us to associate the induced low-energy 
LDOS in the S region purely with OTE pairing.

In SNS junctions we have demonstrated an even stronger functional relationship between superconducting pairing and 
low-energy LDOS inside the N junction region. In this case the pairing amplitudes become dependent on the 
superconducting phase difference ϕ and on the length of the N region. 
In fact, we find that the condition giving the Andreev bound states (ABSs) energies in the junction is exactly the 
same as the condition generating resonant peaks in either the ESE and OTE pairing amplitudes, independent of junction length.
For short junctions, at ϕ=0 the ESE pairing dominates over a completely suppressed OTE pairing, while at 
ϕ=π the ESE is instead zero and the OTE amplitude large. Notably, the pair of ABSs crosses zero at 
ϕ = π, which is directly reflected in the OTE amplitude which has a resonant peak at zero energy. 
We have also shown that the supercurrent across such short junctions acquires its maximum value  
at ϕ=π as a result of the resonant peak at zero energy which is purely due to OTE pairing.
In long junctions more ABSs fit within the N region, but still all ABS energies correspond to resonant peaks in 
the pairing amplitudes. In this case we find the clearest distinction between ESE and OTE pairing in the middle of 
the N region. The ESE is completely suppressed at ϕ = π, while the OTE is finite and has resonant peaks at 
each ABS energy. At ϕ=0 the relation is reversed with ESE having resonances at the ABS energies. Dominant OTE 
pairing exists also at distances away from the middle of the junction and at phase differences away from ϕ=π.

In summary, our findings show that the finite LDOS, as well as local conductance, at the interfaces of NS or NSN junctions at the edge of a 2DTI is entirely a consequence of pure OTE pairing. This results from the unique nature of  the helical edges of 2DTIs, where the Andreev reflection magnitude, strongly connected to LDOS and conductance, at the NS interface reaches its maximum for energies within the gap.<cit.> In SNS junctions the relationship is even stronger with the ABSs directly corresponding to peaks in the pairing amplitudes. In particular, the zero-energy ABS at ϕ = π, protected by topology, is of complete OTE nature. 




§ ACKNOWLEDGEMENTS

We thank P. Burset, D. Kuzmanovski and C. Triola for motivating and helpful discussions.
This work was supported by the Swedish Research Council (Vetenskapsrådet), the Göran Gustafsson Foundation, the Swedish Foundation for Strategic Research (SSF), and the Knut and Alice Wallenberg Foundation through the Wallenberg Academy Fellows program.













§ RETARDED AND ADVANCED GREEN'S FUNCTIONS


The structure of the Green's function G^r is given by Eq. (<ref>) in the main text.
The elements G_ee,hh^r we refer to as the regular parts of G^r, while G_eh,he^r are the anomalous electron-hole parts. 
Electron-hole symmetry imposes for the BdG Hamiltonian PH_BdG^*P^†=-H_BdG, while for the Green's function P[G^r(x,x',ω)]^*P^†=-G^r(x,x',-ω), where P=σ_yτ_y.<cit.> This therefore connects the two diagonal (off-diagonal) elements of G^r. 
The advanced Green's function, G^a, has the same matrix structure as the retarded function given by Eq. (<ref>)
and can be calculated from incoming boundary conditions. Alternatively, we can use the relation between retarded and advanced Green's functions: G^a(x,x',ω)=[G^r(x',x,ω)]^†.

We are interested in the pairing amplitudes, that is, the anomalous (electron-hole) part of the retarded Green's function, 

    G^r_eh(x,x',ω)=
    [  [G^r_eh(x,x',ω)]_↑↓ -[G^r_eh(x,x',ω)]_↑↑;  [G^r_eh(x,x',ω)]_↓↓ -[G^r_eh(x,x',ω)]_↓↑ ] ,

where the minus signs arise due to the specific choice of basis.
In order to decompose the spin symmetry, we write the anomalous Green's function according to Eq. (<ref>) in the main text. Then we arrive at the pairing amplitudes f^r_i:

    f_0^r(x,x',ω)   =[G^r_eh(x,x',ω)]_↑↓-[G^r_eh(x,x',ω)]_↓↑/2 ,
    
    f_3^r(x,x',ω)   =[G^r_eh(x,x',ω)]_↑↓+[G^r_eh(x,x',ω)]_↓↑/2 ,
    
    f_1^r(x,x',ω)   =-[G^r_eh(x,x',ω)]_↑↑+[G^r_eh(x,x',ω)]_↓↓/2 ,
    
    f_2^r(x,x',ω)   =i[G^r_eh(x,x',ω)]_↑↑+[G^r_eh(x,x',ω)]_↓↓/2 .

From these equations we observe that f_0^r acquires a minus sign under the exchange of spins and is thus an odd function under spin exchange, while f^r_1,2,3 acquire a plus sign and are even under spin exchange. 
 f^r_0 is therefore referred to as the spin-singlet component, while f^R_1,2,3 are the spin-triplet components.
 


 §.§ Antisymmetry


The anomalous Green's function represents the wave function of a two-electron system which must obey antisymmetry upon the simultaneous exchange of spins (↑↔↓), spatial coordinates (x↔ x') and time (or energy/frequency) coordinates (t↔ t' or ω→-ω). 

In this work, we use retarded and advanced Green's functions which do  not respect symmetry under frequency. These Green's functions are only partially defined on the time axis and therefore, when the sign of the frequency changes (or time), we should pass from one to the other,<cit.>

    G^r_s_1s_2(x,x',ω)=-G^a_s_2s_1(x',x,-ω) .

By using Eqs. (<ref>) and (<ref>) we get

    [G^r_eh(x,x',ω)]_↑↓   =[G^a_eh(x',x,-ω)]_↓↑ ,
    
    [G^r_eh(x,x',ω)]_↑↑   =-[G^a_eh(x',x,-ω)]_↑↑ ,
    
    [G^r_eh(x,x',ω)]_↓↓   =-[G^a_eh(x',x,-ω)]_↓↓ ,
    
    [G^r_eh(x,x',ω)]_↓↑   =[G^a_eh(x',x,-ω)]_↑↓ .

Similarly, we arrive at the conditions of antisymmetry for the spin-singlet and triplet amplitudes

    f^r_0(x,x',ω)   =f_0^a(x',x,-ω) ,
    
    f^r_1(x,x',ω)   =-f_1^a(x',x,-ω) ,
    
    f^r_2(x,x',ω)   =-f_2^a(x',x,-ω) ,
    
    f^r_3(x,x',ω)   =-f_3^a(x',x,-ω) .

These equations represent the full antisymmetry conditions under the exchange of spin, spatial coordinates, and frequency imposed by Fermi-Dirac statistics on the pairing amplitudes and they are presented in the main text as Eqs. (<ref>).




§ NS JUNCTION


We model a NS junction at the metallic edge of a 2DTI by considering a step-like profile of the induced pairing potential, with the normal region at x<0 and the superconducting region at  x>0, 

    Δ(x)=θ(x)Δ=
       0  ,      x<0, 
    Δ ,     x>0 ,

where we can set the overall superconducting phase to zero. 



 §.§ Scattering processes

In general, there are four different scattering processes at the NS interface which read

    Ψ_1(x)   =
    ϕ_1^N  e^ik_ex+a_1ϕ_3^N  e^ik_hx +b_1ϕ_2^N  e^-ik_ex, x<0   
    
        c_1ϕ_1^S  e^ik_e^Sx +d_1ϕ_4^S  e^-ik_h^Sx, x>0    
    Ψ_2(x)   =
    ϕ_4^N  e^-ik_hx+a_2ϕ_2^N  e^-ik_ex +b_2ϕ_3^N  e^ik_hx, x<0   
    
        c_2ϕ_4^S  e^-ik_h^Sx +    d_2ϕ_1^S  e^ik_e^Sx,  x>0    
    Ψ_3(x)   =
    
        c_3ϕ_2^N  e^-ik_ex +d_3ϕ_3^N  e^ik_hx,  x<0   
    ϕ_2^S  e^-ik_e^Sx+a_3ϕ_4^S  e^-ik_h^Sx+b_3ϕ_1^S  e^ik_e^Sx, x>0   
    Ψ_4(x)   =
    
        c_4ϕ_3^N  e^ik_hx +    d_4ϕ_2^N  e^-ik_ex, x<0     
    ϕ_3^S  e^ik_h^Sx+a_4ϕ_1^S  e^ik_e^Sx+b_4ϕ_4^S  e^-ik_h^Sx, x>0

while the conjugated processes are

    Ψ̃_1(x')   =
    ϕ̃_1^N  e^ik_ex'+ã_1ϕ̃_3^N e^ik_hx' +b̃_1ϕ̃_2^N e^-ik_ex',  x<0   
    c̃_1ϕ̃_1^S e^ik_e^Sx' +  d̃_1ϕ̃_4^S  e^-ik_h^Sx', x>0     
    Ψ̃_2(x')   =
    ϕ̃_4^N e^-ik_hx'+ã_2ϕ̃_2^N e^-ik_ex' +b̃_2ϕ̃_3^N e^ik_hx', x<0   
    c̃_2ϕ̃_4^S e^-ik_h^Sx' + d̃_2ϕ̃_1^S e^ik_e^Sx', x>0    
    Ψ̃_3(x')   =
    c̃_3ϕ̃_2^N e^-ik_ex' +d̃_3ϕ̃_3^N e^ik_hx',  x<0    
    ϕ̃_2^S e^-ik_e^Sx'+ã_3ϕ̃_4^S e^-ik_h^Sx'+b̃_3ϕ̃_1^S e^ik_e^Sx', x>0   
    Ψ̃_4(x')   =
    c̃_4ϕ̃_3^N e^ik_hx' +    d̃_4ϕ̃_2^N e^-ik_ex',  x<0     
    ϕ̃_3^S e^ik_h^Sx'+ã_4ϕ̃_1^S  e^ik_e^Sx'+b̃_4ϕ̃_4^S e^-ik_h^Sx',  x>0 ,

where 

    ϕ_1^N   =[ 1,
    0,
    0,
    0 ]^T,
    ϕ_2^N   =[ 0,
    1,
    0,
    0 ]^T,
    ϕ_3^N   =[ 0,
    0,
    1,
    0 ]^T,
    ϕ_4^N   =[ 0,
    0,
    0,
    1 ]^T,
        ϕ_1^S   =[ u,
    0,
    v,
    0 ]^T,
    ϕ_2^S   =[ 0,
    u,
    0,
    v ]^T,
    ϕ_3^S   =[ v,
    0,
    u,
    0 ]^T,
    ϕ_4^S   =[ 0,
    v,
    0,
    u ]^T,

are wave functions of H_BdG(k), while 

    ϕ̃_1^N   =[ 0,
    1,
    0,
    0 ]^T,
    ϕ̃_2^N   =[ 1,
    0,
    0,
    0 ]^T,
    ϕ̃_3^N   =[ 0,
    0,
    0,
    1 ]^T,
    ϕ̃_4^N   =[ 0,
    0,
    1,
    0 ]^T,
       ϕ̃_1^S   =
    [ 0,
    u,
    0,
    v ]^T,
    ϕ̃_2^S   =
    [ u,
    0,
    v,
    0 ]^T,
    ϕ̃_3^S   =
    [ 0,
    v,
    0,
    u ]^T,
    ϕ̃_4^S   =
    [ v,
    0,
    u,
    0 ]^T,

are the wave functions of the conjugated Hamiltonian H̃_BdG(k).  
Notice that the conjugated scattering processes are constructed after solving the eigenvalue problem for H̃_BdG(k)=H_BdG^*(-k)=H_BdG^T(-k) instead of Eq. (<ref>) in the main text.
In these equations we have used the following relations:

    k_e,h   =(μ±ω)/v_F ,
    
    k_e,h^S   =(μ±√(ω^2-Δ^2))/v_F=k_μ± k(ω) ,
    
     u,v   =√(1/2[1±√(ω^2-Δ^2)/ω]) .

In this work we mainly focus on energies within Δ and then 
 k_e,h^S=k_μ± iκ(ω), with κ(ω)=√(Δ^2-ω^2)/v_F and the coherence factors can be written as u=√(Δ/2ω)  e^iη(ω)/2 and v=√(Δ/2ω)  e^-iη(ω)/2, where η(ω)= arccos(ω/Δ).
 
At the NS interface, Ψ_1 represents the following process: 
an incoming electron (right-moving with spin up) from the N region with wave function ϕ_1^N  e^ik_ex 
experiences reflection and transmission at the NS interface with certain probabilities.
It can be reflected into a left-moving electron with spin-down with wave function 
ϕ_2^N  e^-ik_ex and amplitude b_1 or Andreev reflected into a left-moving hole with spin down with wave function ϕ_3^N  e^ik_hx and amplitude a_1, or transmitted to into the S region in the form of a right-moving quasielectron with wave function 
ϕ_1^S  e^ik_e^Sx and amplitude c_1 or a right-moving quasihole with wave function ϕ_4^S  e^-ik_h^Sx and amplitude d_1. 
Thus, a_1, b_1, c_1 and d_1 are the amplitudes of reflection into a hole (Andreev reflection), reflection into an electron (normal reflection), transmission into an electron and transmission into a hole. 
Likewise, Ψ_2,3,4 correspond to scattering processes for an incoming hole from the N region and incoming electron or incoming hole from the S region, respectively. 
 
The processes Ψ_i and Ψ̃_i are fully determined after finding the coefficients a_i, b_i, c_i and d_i. These in turn are calculated by matching the functions Ψ_i at the NS interface x=0,

    [Ψ_i(x>0)]_x=0   =[Ψ_i(x<0)]_x=0 .

Each Ψ_i is a four column vector and therefore provides four equations. At the end we have a system of 16 equations for the 16 unknown coefficients a_i, b_i, c_i and d_i, which is generally solvable and completely determines the  scattering states Ψ_i.  
The same holds for the conjugated processes Ψ̃_i. 
Notice that we here have written the scattering wave functions, Eqs. (<ref>-<ref>), 
in a general form. However, the processes with amplitudes b_i and d_i are forbidden by helicity 
of the 2DTI edge; an incident electron can be only reflected as a hole by a superconducting barrier or 
transmitted as an electron through it. 
Thus, normal reflection and non-local Andreev transmission are forbidden by helicity conservation<cit.> and b_i=d_i=0.  For the Andreev reflection amplitudes we generally also obtain a_1,2=v/u=-a_3,4.



 §.§ Green's functions and pairing amplitudes





  §.§.§ Normal region


After finding Ψ_i and Ψ̃_i, given by Eqs. (<ref>) and (<ref>), we construct the retarded Green's functions using Eq. (<ref>) in the main text.
In the normal (N) region, we obtain the following expressions for the regular and anomalous parts

    G_ee^r(x,x',ω)
       =
    1/iv_f[ G_ee,↑↑^r         0;         0 G_ee,↓↓^r ],
    
    G_hh^r(x,x',ω)
       =1/iv_f[ G_hh,↓↓^r         0;         0 G_hh,↑↑^r;           ],
    
    G_eh^r(x,x',ω)
       =
    1/iv_f[                       0                       0;                       0 a_1(ω) e^-i(k_ex-k_hx');                         ],
    
    G_he^r(x,x',ω)
       =
    1/iv_f[ a_1(ω) e^i(k_hx-k_ex')                      0;                      0                      0 ],

with G_ee,↑↑^r=θ(x-x') e^ik_e(x-x'), 
G_ee,↓↓^r=θ(x'-x) e^-ik_e(x-x'), 
G_hh,↓↓^r=θ(x'-x) e^ik_h(x-x'),
G_hh,↑↑^r=θ(x-x') e^-ik_h(x-x') and 
a_1=v/u the Andreev reflection amplitude for an incident electron from N. For energies within the gap a_1(ω)= e^-iη(ω).

Now, by using G^r_ee, we calculate the LDOS in the N region following Eq. (<ref>) in the main text and get

    ρ_ N(x,ω)
       =1/π v_F ,

where we have used lim_x→ x'θ(x-x')=1/2. 
Notice how the LDOS is independent of both energy and position, as is expected for the helical edge state in a 2DTI.

By decomposing the spin symmetry of the anomalous part of the retarded and advanced Green's functions, employing Eq. (<ref>), we get the pairing amplitudes

    f_0, N^r(x,x',ω)   =a_1(ω)/2i v_F e^-iμ(x-x')+ω(x+x')/v_F
    ,
    
    f_3, N^r(x,x',ω)   =-f_0^r(x,x',ω),

and f_1,2, N^r=0.
As discussed in Sec. <ref>, the pairing amplitudes must obey Fermi-Dirac statistics. We use Eqs. (<ref>) and check the antisymmetry of previous pairing functions.
For example, for the singlet component we have

    f^a_0, N(x',x,-ω)   =-a_1^*(-ω)/2iv_F e^iμ(x'-x)-ω(x+x')/v_F ,
       =a_1(ω)/2iv_F e^-iμ(x-x')+ω(x+x')/v_F ,
       =f_0^r(x,x',ω) ,

where we have used a_1^*(-ω)= e^iη(-ω)=- e^-iη(ω)=-a_1(ω). Thus, we conclude that f_0^r is antisymmetric. 
Likewise, the triplet component is fully antisymmetric and obeys the relation given by Eq. (<ref>). Notice that in order to check antisymmetry we have to use the advanced pairing functions as discussed in Appendix <ref>. 

We now use Eqs. (<ref>) in order to decompose into the even and odd-frequency components

    f^r, E_0, N(x,x',ω)   =a_2(ω)/2iv_F C_xx'  e^-iω(x+x')/v_F,
    
    f^r, O_0, N(x,x',ω)   =-a_2(ω)/2v_F S_xx'  e^-iω(x+x')/v_F,
    
    f^r, E_3, N(x,x',ω)   =a_2(ω)/2v_F S_xx'  e^-iω(x+x')/v_F,
    
    f^r, O_3, N(x,x',ω)   =-a_2(ω)/2iv_F C_xx'  e^-iω(x+x')/v_F,

where C_xx'= cos[k_μ(x-x')], S_xx'= sin[k_μ(x-x')], and a_2=a_1 the Andreev reflection amplitude. These are Eqs. (<ref>) given in the main text.

Before going further it is worth to point out the following. To check the antisymmetry of the retarded even- and odd-frequency pairing amplitudes is not trivial and it is important to write down their respective advances functions. We have explicitly verified that the given expressions  follow the antisymmetry relations. For a quick check, however, we can focus on the parity in space which follows directly from cos[k_μ(x-x')] and sin[k_μ(x-x')]. Since the spin is already explicit this means the only symmetry left is that of the frequency. For example, for OSO we see that sin[k_μ(x-x')] makes the amplitude odd in space and since it is a spin singlet, the only possibility is for an odd-frequency dependence. 


 


  §.§.§ Superconducting region


In the superconducting region (S) we proceed similarly as in N. We obtain the retarded Green's function, which for energies below the gap Δ, reads


    G^r_ee,↑↑(x,x',ω)   =Z(ω)
      e^ik_μ(x-x')
    [
    ã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|N(x,x',ω)] ,
    
    G^r_ee,↓↓(x,x',ω)   =
    Z(ω)
      e^-ik_μ(x-x')[ã_3(ω)  e^-κ(ω)(x+x')
    + e^-κ(ω)|x-x'|N(x',x,ω)]
     ,
    
    G^r_hh,↓↓(x,x',ω)   =
    Z(ω)
      e^ik_μ(x-x')[ã_3(ω)  e^-κ(ω)(x+x')
    + e^-κ(ω)|x-x'|N(x',x,ω)]
     ,
    
    G^r_hh,↑↑(x,x',ω)   =Z(ω)
      e^-ik_μ(x-x')
    [
    ã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|N(x,x',ω)] ,
    
    G^r_eh,↑↓(x,x',ω)   =
    Z(ω)
      e^ik_μ(x-x')
    [
    u/vã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|] ,
    
    G^r_eh,↓↑(x,x',ω)   =
    Z(ω)
      e^-ik_μ(x-x')
    [
    v/uã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|] ,
    
    G^r_he,↓↑(x,x',ω)   =
    Z(ω)
      e^ik_μ(x-x')
    [
    v/uã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|] ,
    
    G^r_he,↑↓(x,x',ω)   =
    Z(ω)
      e^-ik_μ(x-x')
    [
    u/vã_3(ω)  e^-κ(ω)(x+x')
    +  e^-κ(ω)|x-x'|] ,
    
    G^r_eh(he),↑↑,↓↓(x,x',ω)   =0 ,  
    G^r_ee(hh),↑↓,↓↑(x,x',ω)=0 ,


where N(x,x',ω)=θ(x-x')(u/v)+θ(x'-x)(v/u), Z(ω)=(1/iv_F)/[ e^iη(ω)- e^-iη(ω)], and ã_3=a_3=-a_2.

Notice that, while in the anomalous parts (eh, he) only mixed spin components are finite, in the regular parts (ee, hh) we also obtain finite equal spin components. If we were to consider finite magnetic order in the system the situation would dramatically change with additional off-diagonal (mixed spin) terms in the regular part and equal spin terms in the anomalous components.<cit.>

By using the regular electron-electron part of G^r, given by the first two expressions of Eqs. (<ref>), we obtain the LDOS in the superconducting region,

    ρ_ S(x,ω)
       = Im{i/π v_F[ ρ̅+(1- ρ̅) e^ik(ω)2x]} ,

where ρ̅(ω)=ω/√(ω^2-Δ^2), and k(ω)=iκ(ω).
Within the gap, ω^2<Δ^2, this reduces to

    ρ_ S(x,|ω|<|Δ|)
       =1/π v_F e^-κ(ω)2x ,

which is given as Eq. (<ref>) in the main text.

The pairing amplitudes are calculated by decomposing the spin components according to Eq. (<ref>) and we arrive at


    f_0,S^r(x,x',ω)   =
    Z(ω) e^-κ(ω)|x-x'| C_xx'
    +
    Z(ω)ã_3(ω)/2 e^-κ(ω)(x+x')[
     e^iη(ω) e^ik_μ(x-x')+ e^-iη(ω) e^-ik_μ(x-x')] ,
    
    f_3,S^r(x,x',ω)   =
    Z(ω) e^-κ(ω)|x-x'| iS_xx'
    +
    Z(ω)ã_3(ω)/2 e^-κ(ω)(x+x')[
     e^iη(ω) e^ik_μ(x-x')- e^-iη(ω) e^-ik_μ(x-x')] ,


with f_1,2,S^r=0.

It is important to notice that the pairing amplitudes given by Eqs. (<ref>) contain two terms, which arise from different parts in the S region. The first term is proportional to e^-κ(ω)|x-x'| and is associated with the bulk deep inside the S region. Indeed, locally (x = x'), such a term becomes independent of the space coordinate. The second term is proportional to the Andreev reflection amplitude a_3 and therefore arises due the presence of the NS interface. This term also keeps a spatial dependence even locally that gives the decay into the S region.
Thus, in the S region, the pairing amplitudes can be written with
as bulk (B)  and interface (I) components


    f_0, B^r(x,x',ω)   =
    Z(ω) e^-κ(ω)|x-x'| C_xx' ,
    
    f_3, B^r(x,x',ω)   =
    Z(ω) e^-κ(ω)|x-x'| iS_xx' ,
    
    f_0, I^r(x,x',ω)   =
    Z(ω)ã_3(ω)/2 e^-κ(ω)(x+x')
    K_+ ,
    
    f_3, I^r(x,x',ω)   =
    Z(ω)ã_3(ω)/2 e^-κ(ω)(x+x')
    K_- ,

where we have used K_±(ω,x,x')=
 e^iη(ω) e^ik_μ(x-x')± e^-iη(ω) e^-ik_μ(x-x'). 

We have explicitly checked that all these pairing amplitudes are fully antisymmetric, obeying Eqs. (<ref>). For example, for f_3, I^r we find

    f_3, I^a(x',x,-ω)   =
    Z^∗(-ω)ã^*_3(-ω)/2 e^-κ(-ω)(x+x')
       × K_-(-ω,x',x) ,
       =-Z(ω)ã_3(ω)/2 e^-κ(ω)(x+x')
       × K_-(ω,x,x') ,
       =-f_3, I^r(x,x',ω)

where we have used ã_3^*(-ω)=-a_3(ω), Z^∗(-ω)=Z(ω) and K_-(-ω,x',x)=K_-(ω,x,x'). We have also used the advanced pairing function f^a_3, I which is  calculated from the relation between retarded and advanced Green's functions G^a(x,x',ω)=[G^r(x',x,ω)]^†.

We now use Eqs. (<ref>) in the main text to isolate the even- and odd-frequency components, which gives in the bulk of the S region

    f^r, O_0,3, B(x,x',ω)   =0 ,
    
    f^r, E_0, B(x,x',ω)   =Z(ω) e^-κ(ω)|x-x'| C_xx' ,
    
    f^r, E_3, B(x,x',ω)   =iZ(ω) e^-κ(ω)|x-x'| S_xx' .

and close to the NS interface on the S side

    f^r, E_0, I(x,x',ω)   =a_3(ω)/2iv_F e^-κ(ω)(x+x')
    B(ω)
    C_xx' ,
    
    f^r, O_0, I(x,x',ω)   =a_3(ω)/2v_F e^-κ(ω)(x+x') S_xx' ,
    
    f^r, E_3, I(x,x',ω)   =a_3(ω)/2v_F e^-κ(ω)(x+x')
    B(ω)
    S_xx' ,
    
    f^r, O_3, I(x,x',ω)   =a_3(ω)/2iv_F e^-κ(ω)(x+x') C_xx' ,

where B(ω)=[ e^iη(ω)+ e^-iη(ω)]/[ e^iη(ω)- e^-iη(ω)]. 
Notice that in the bulk we only obtain trivial even-frequency spin-singlet and -triplet components, namely the ESE and ETO pairing amplitudes.
At the interface, however, hosts all symmetry classes (ESE, OSO, ETO, and OTE). Previous two sets of equations, for the bulk and interface, correspond to Eqs. (<ref>-<ref>) in the main text. 





§ NSN


Here we consider a NSN junction, where the S region has a finite length L_ S restricted to 0<x<L_ S, i.e., 

    Δ(x)=
       0  ,      x<0 , 
    Δ ,     0<x<L_ S ,
    
         0  ,      x>0 .




 §.§ Scattering processes

The four scattering processes in a NSN junction read

    Ψ_1(x)   =
    ϕ_1^N e^ik_ex+a_1ϕ_3^N e^ik_hx+b_1ϕ_2^N e^-ik_ex,  x<0   
    ∑_ip_iϕ_i^S e^ik_i^Sx, 0<x<L_ S   
    
       c_1ϕ_1^N e^ik_ex +d_1ϕ_4^N e^-ik_hx, x>L_ S   
    Ψ_2(x)   =
    ϕ_4^N e^-ik_hx+a_2ϕ_2^N e^-ik_ex+b_2ϕ_3^N e^ik_hx,  x<0   
    ∑_iq_iϕ_i^S e^ik_i^Sx, 0<x<L_ S   
    
       c_2ϕ_4^N e^-ik_hx +d_2ϕ_1^N e^ik_ex, x>L_ S   
    Ψ_3(x)   =
    
     c_3ϕ_2^N e^-ik_ex +d_3ϕ_3^N e^ik_hx,  x<0   
    ∑_ir_iϕ_i^S e^ik_i^Sx, 0<x<L_ S   
    ϕ_2^N e^-ik_ex+a_3ϕ_4^N e^-ik_hx+b_3ϕ_1^N e^ik_ex, x>L_ S   
    Ψ_4(x)   =
    
     c_4ϕ_3^N e^ik_hx +d_4ϕ_2^N e^-ik_ex,  x<0   
    ∑_is_iϕ_i^S  e^ik_i^Sx, 0<x<L_ S   
    ϕ_3^N e^ik_hx+a_4ϕ_1^N e^ik_ex+b_4ϕ_4^N e^-ik_hx, x>L_ S ,

while the conjugated processes are

    Ψ̃_1(x)   =
    ϕ̃_1^N e^ik_ex+ã_1ϕ̃_3^N e^ik_hx+b̃_1ϕ̃_2^N e^-ik_ex, x<0   
    ∑_ip̃_iϕ̃_i^S e^ik_i^Sx, 0<x<L_ S   
    c̃_1ϕ̃_1^N e^ik_ex +d̃_1ϕ̃_4^N e^-ik_hx, x>L_ S ,
    
    Ψ̃_2(x)   =
    ϕ̃_4^N e^-ik_hx+ã_2ϕ̃_2^N e^-ik_ex+b̃_2ϕ̃_3^N e^ik_hx,  x<0   
    ∑_iq̃_iϕ̃_i^S e^ik_i^Sx, 0<x<L_ S   
    c̃_2ϕ̃_4^N e^-ik_hx +d̃_2ϕ̃_1^N e^ik_ex, x>L_ S   
    Ψ̃_3(x)   =
    c̃_3ϕ̃_2^N e^-ik_ex +d̃_3ϕ̃_3^N e^ik_hx,  x<0   
    ∑_ir̃_iϕ̃_i^S e^ik_i^Sx, 0<x<L_ S   
    ϕ̃_2^N e^-ik_ex+ã_3ϕ̃_4^N e^-ik_hx+b̃_3ϕ̃_1^N e^ik_ex, x>L_ S   
    Ψ̃_4(x)   =
    c̃_4ϕ̃_3^N e^ik_hx +d̃_4ϕ̃_2^N e^-ik_ex,  x<0    
    ∑_is̃_iϕ̃_i^S e^ik_i^Sx, 0<x<L_ S   
    ϕ̃_3^N e^ik_hx+ã_4ϕ̃_1^N e^ik_ex+b̃_4ϕ̃_4^N e^-ik_hx, x>L_ S

where i=1,2,3,4, k_1,2^S=± k_e^S, k_3,4^S=± k_h^S 
and the spinors ϕ_i^N,S and ϕ̃_i^N,S are given by Eqs. (<ref>) 
and (<ref>), respectively. 

Notice that, as for NS junctions, in this part we have written the general form of the scattering wave functions in the N and S regions. 
The S region is in general formed by a linear combination of four elements given in terms of the amplitudes p_i, q_i, r_i and s_i. Their meaning is as follows: p_1 represent electron transmission from left N to S, p_2 the normal reflection at the right SN interface, p_3 the Andreev reflection at the SN interface, and p_4 the Andreev reflection at the left NS interface. Similar ideas apply to the amplitudes q_i, r_i and s_i.
The amplitudes of all these processes are calculated after matching the wave-functions at the left NS and right SN interfaces.
Due to  helicity conservation in the 2DTI edge states we can again directly obtain 
b_i=d_i=b̃_i=d̃_i=0 and p_2,4=s_2,4=r_1,3=q_1,3=0. Similar relations hold for the conjugated processes.



 §.§ Green's function and pairing amplitudes



  §.§.§ Normal regions


Finding Ψ_i and Ψ̃_i, given by Eqs. (<ref>) and (<ref>), allows us to construct the retarded Green's functions using Eq. (<ref>) in the main text.
In the two normal regions the Green's functions and pairing amplitudes acquire the same form as in the NS junction, with the only difference that the Andreev reflection amplitude becomes dependent on the length L_ S of the S region and reads  a_1(ω,L_ S)= a_2(ω,L_ S)= sin[iκ(ω)L]/ sin[iκ(ω)L-η(ω)].  
Notice that a_2(ω,L_ S)→ a_2(ω)=v/u for L_ S→∞, as expected.



  §.§.§ Superconducting region


Due to the finite length of the S region in NSN junctions we find quite different results here compared to the NS junction. In the S region the retarded Green's function contains the components



    G^r_ee,↑↑(x,x',ω)   = e^i(k_h^Sx-k_e^Sx')γ_2+
     e^i(k_e^Sx-k_h^Sx')γ_3+θ(x-x')[ e^ik_e^S(x-x')γ_1+ e^ik_h^S(x-x')γ_4]
       -θ(x'-x)[ e^ik_e^S(x-x')γ_2+ e^ik_h^S(x-x')γ_3] ,
    
    G^r_ee,↓↓(x,x',ω)   = e^i(k_e^Sx'-k_h^Sx)γ_3+
     e^i(k_h^Sx'-k_e^Sx)γ_2-θ(x-x')[ e^-ik_e^S(x-x')γ_2+ e^-ik_h^S(x-x')γ_3]
       +θ(x'-x)[ e^-ik_e^S(x-x')γ_1+ e^-ik_h^S(x-x')γ_4] ,
    
    G^r_hh,↓↓(x,x',ω)   = e^i(k_h^Sx-k_e^Sx')γ_2+
     e^i(k_e^Sx-k_h^Sx')γ_3-θ(x-x')[ e^ik_h^S(x-x')γ_2+ e^ik_e^S(x-x')γ_3]
       +θ(x'-x)[ e^ik_h^S(x-x')γ_1+ e^ik_e^S(x-x')γ_4] ,
    
    G^r_hh,↑↑(x,x',ω)   = e^i(k_e^Sx'-k_h^Sx)γ_3+
     e^i(k_h^Sx'-k_e^Sx)γ_2+θ(x-x')[ e^-ik_h^S(x-x')γ_1+ e^-ik_e^S(x-x')γ_4]
       -θ(x'-x)[ e^-ik_h^S(x-x')γ_2+ e^-ik_e^S(x-x')γ_3] ,
    
    G_eh,↑↓^r(x,x',ω)   =- e^i(k_e^Sx-k_h^Sx')β_2-
     e^i(k_h^Sx-k_e^Sx')β_3+ e^i k_e^S(x-x')[θ(x-x')β_2+θ(x'-x)β_3]
       + e^i k_h^S(x-x')[θ(x-x')β_3+θ(x'-x)β_2] ,
    
    G_eh,↓↑^r(x,x',ω)   = e^i(k_h^Sx'-k_e^Sx)β_1+
     e^i(k_e^Sx'-k_h^Sx)β_4+ e^-i k_e^S(x-x')[θ(x-x')β_3+θ(x'-x)β_2]
       + e^-i k_h^S(x-x')[θ(x-x')β_2+θ(x'-x)β_3] ,
    
    G_he,↓↑^r(x,x',ω)   = e^i(k_h^Sx-k_e^Sx')β_1+
     e^i(k_e^Sx-k_h^Sx')β_4+ e^i k_e^S(x-x')[θ(x-x')β_2+θ(x'-x)β_3]
       + e^i k_h^S(x-x')[θ(x-x')β_3+θ(x'-x)β_2] ,
    
    G_he,↑↓^r(x,x',ω)   =- e^i(k_e^Sx'-k_h^Sx)β_2-
     e^i(k_h^Sx'-k_e^Sx)β_3+ e^-i k_e^S(x-x')[θ(x-x')β_3+θ(x'-x)β_2]
       + e^-i k_h^S(x-x')[θ(x-x')β_2+θ(x'-x)β_3] ,

where we have defined

    γ_1   =1/iv_fu^4/(u^2-v^2)[u^2- e^i(k_e^S-k_h^S)Lv^2] ,  γ_2=1/iv_fu^2v^2/(u^2-v^2)[-u^2 e^i(k_h^S-k_e^S)L+v^2] ,
    γ_3   =1/iv_fu^2v^2/(u^2-v^2)[-u^2+ e^i(k_e^S-k_h^S)Lv^2] ,  γ_4=1/iv_fv^4/(u^2-v^2)[u^2 e^i(k_h^S-k_e^S)L-v^2] ,
    β_1   =
    1/iv_fu^3v/(u^2-v^2)[-u^2 e^i(k_h^S-k_e^S)L+v^2] ,  β_2=1/iv_fu^3v/(u^2-v^2)[u^2- e^i(k_e^S-k_h^S)Lv^2] ,
    β_3   =1/iv_fu v^3/(u^2-v^2)[u^2 e^i(k_h^S-k_e^S)L-v^2] ,  β_4=1/iv_fu v^3/(u^2-v^2)[-u^2+ e^i(k_e^S-k_h^S)Lv^2] .

The regular electron-electron part of G^r gives the LDOS in the S region which reads 
    ρ_S(x,ω)=
    -1/π v_F Im{γ_1/i[v^2/u^2 e^2ik(ω)L(1+v^2/u^2-2 e^-2ik(ω)x)+(1+v^2/u^2-2v^2/u^2 e^2ik(ω)x)]  } .

For energies within Δ this expression reduces to 

    ρ_ S(x,|ω|<Δ)=ω^2(1- e^-2κ(ω)L)^2+Δ^2(1+ e^-2κ(ω)L)[ e^-2κ(ω)L( e^2κ(ω)x-ω^2/Δ^2) +( e^-2κ(ω)x-ω^2/Δ^2)]/π v_F[Δ^2(1+ e^-2κ(ω)L)^2-4ω^2 e^-2κ(ω)L] .

At ω=0, this expression can be simplified even further and we obtain 

    ρ_ S(x,0)=1/π v_F e^-2L/ξ e^2x/ξ+ e^-2x/ξ/1+ e^-2L/ξ ,

where we have used ξ = v_F/Δ. This is the expression for the LDOS given in Eq. (<ref>) in the main text.

For the pairing amplitudes, limiting ourselves to energies within the gap, we arrive at

    f_0, S^r(x,x',ω)   =
    -β_2 e^-κ(ω)(x+x')[ e^-ik_μ(x-x')-2iη(ω) + e^ik_μ(x-x')]+β_1 e^κ(ω)(x+x')[ e^-ik_μ(x-x') + e^ik_μ(x-x')-2iη(ω)]
       +2 cos[k_μ(x-x')][β_2 e^-κ(ω)|x-x'|-β_1 e^-2iη(ω) e^κ(ω)|x-x'|] ,
    
    f_3, S^r(x,x',ω)   =
    -β_2 e^-κ(ω)(x+x')[ e^ik_μ(x-x')- e^-ik_μ(x-x')-2iη(ω)]+
    β_1 e^κ(ω)(x+x')[ e^ik_μ(x-x')-2iη(ω) - e^-ik_μ(x-x')]
       +2i sin[k_μ(x-x')][β_2 e^-κ(ω)|x-x'|-β_1 e^-2iη(ω) e^κ(ω)|x-x'|] ,

and f^r_1, S = f^r_2, S = 0, where we have used γ_4,2=- e^-2iη(ω)γ_2,1 and β_4,3=- e^-2iη(ω)β_2,1, 
with 

    β_1=-iZ(ω) e^iη(ω)-κ(ω)L_ S/4 sin[iκ(ω)L_ S+η(ω)] ,

β_2=-β_1 e^2κ(ω)L_ S, γ_1= e^iη(ω)β_2,  and γ_2= e^-iη(ω)β_1.

From these and the related advanced expressions we can extract the even- and odd-frequency components which reads

    f_0, S^r, E(x,x',ω)   =C_xx'[ e^-κ(ω)(x+x')β_42^-+ e^κ(ω)(x+x')β_13^-+2β_2 e^-κ(ω)|x-x'|+2β_3 e^κ(ω)|x-x'|] ,
    
    f_0, S^r, O(x,x',ω)   =-iS_xx'[ e^-κ(ω)(x+x')β_42^++ e^κ(ω)(x+x')β_13^+ ,
    
    f_3, S^r, E(x,x',ω)   =iS_xx'[ e^-κ(ω)(x+x')β_42^-+ e^κ(ω)(x+x')β_13^-+2β_2 e^-κ(ω)|x-x'|+2β_3 e^κ(ω)|x-x'|] ,
    
    f_3, S^r, O(x,x',ω)   =-C_xx'[ e^-κ(ω)(x+x')β_42^++ e^κ(ω)(x+x')β_13^+] ,

where β_42(13)^±=β_4(1)±β_2(3).

These pairing amplitudes in S can be composed into bulk (B) and interface (I) contributions. Bulk contributions we designate terms which are independent of the average distance from either interface, i.e. with an overall |x-x'| spatial dependence. The remaining terms we label interface contributions, as they all have a decay length 1/[2κ(ω)] from the interface.
With this division we directly arrive at the results given in Eqs. (<ref>)-(<ref>) in the main text. 






§ SNS JUNCTION


Finally we consider SNS junctions, where the normal region has a finite length L_ N, while the S regions are semi-infinite such that

    Δ(x)=Δ ,      x<0, 
    
        0 ,     x<0<L_ N,
    Δ  e^iϕ ,      x>0,

where ϕ is the superconducting phase difference across the junction.



 §.§ Scattering processes

In this case the processes read

    Ψ_1(x)   =
    ϕ_1^S_L  e^ik^S_ex+a_1 ϕ_3^S_L  e^ik^S_hx+b_1 ϕ_2^S_L  e^-ik_e^Sx,  x<0   
    ∑_ip_i ϕ_i^N  e^ik_ix,  0<x<L_ N   
    
        c_1 ϕ_1^S_R  e^ik_e^Sx +d_1 ϕ_4^S_R  e^-ik_h^Sx, x>L_ N   
    Ψ_2(x)   =
    ϕ_4^S_L e^-ik_h^Sx+a_2ϕ_2^S_L e^-ik^S_ex+b_2ϕ_3^S_L e^ik_h^Sx,  x<0   
    ∑_iq_iϕ_i^N e^ik_ix,  0<x<L_ N   
    
       c_2ϕ_4^S_R e^-ik^S_hx +d_2ϕ_1^S_R e^ik^S_ex, x>L_ N   
    Ψ_3(x)   =
    
     c_3ϕ_2^S_L e^-ik_e^Sx +d_3ϕ_3^S_L e^ik^S_hx,  x<0   
    ∑_ir_iϕ_i^N e^ik_ix,  0<x<L_ N   
    ϕ_2^S_R e^-ik^S_ex+a_3ϕ_4^S_R e^-ik^S_hx+b_3ϕ_1^S_R e^ik^S_ex, x>L_ N   
    Ψ_4(x)   =
    
     c_4ϕ_3^S_L e^ik^S_hx +d_4ϕ_2^S_L e^-ik^S_ex,  x<0   
    ∑_is_iϕ_i^N e^ik_ix,  0<x<L_ N   
    ϕ_3^S_R e^ik^S_hx+a_4ϕ_1^S_R e^ik_e^Sx+b_4ϕ_4^S_R e^-ik^S_hx, x>L_ N



    Ψ̃_1(x)   =
    ϕ̃_1^S_L e^ik^S_ex+ã_1ϕ̃_3^S_L  e^ik^S_hx+b̃_1ϕ̃_2^S_L e^-ik^S_ex,  x<0   
    ∑_ip̃_iϕ̃_i^N e^ik_ix,  0<x<L_ N   
    c̃_1ϕ̃_1^S_R e^ik_e^Sx +d̃_1ϕ̃_4^S_R e^-ik^S_hx, x>L_ N   
    Ψ̃_2(x)   =
    ϕ̃_4^S_L e^-ik^S_hx+ã_2ϕ̃_2^S_L e^-ik_e^Sx+b̃_2ϕ̃_3^S_L e^ik^S_hx,  x<0   
    ∑_iq̃_iϕ̃_i^N e^ik_ix,  0<x<L_ N   
    c̃_2ϕ̃_4^S_R e^-ik^S_hx +d̃_2ϕ̃_1^S_R e^ik^S_ex, x>L_ N   
    Ψ̃_3(x)   =
    c̃_3ϕ̃_2^S_L e^-ik^S_ex +d̃_3ϕ̃_3^S_L e^ik^S_hx,  x<0   
    ∑_ir̃_iϕ̃_i^N e^ik_ix,  0<x<L_ N   
    ϕ̃_2^S_R e^-ik^S_ex+ã_3ϕ̃_4^S_R e^-ik^S_hx+b̃_3ϕ̃_1^S_R e^ik^S_ex, x>L_ N   
    Ψ̃_4(x)   =
    c̃_4ϕ̃_3^S_L e^ik_h^Sx +d̃_4ϕ̃_2^S_L e^-ik_e^Sx,  x<0   
    ∑_is̃_iϕ̃_i^N e^ik_ix,  0<x<L_ N   
    ϕ̃_3^S_R e^ik_h^Sx+ã_4ϕ̃_1^S_R e^ik_e^Sx+b̃_4ϕ̃_4^S_R e^-ik^S_hx
    , x>L_ N

where i=1,2,3,4 and the spinors ϕ_i^N, ϕ̃_i^N, 
ϕ_i^S_L, ϕ̃_i^S_L acquire the same form as in the NS junction, 
given by Eqs. (<ref>) and (<ref>). In the right superconductor, 
the wave functions of H_ BdG read

    ϕ_1^S_R   =[ u  e^iϕ/2,
    0,
    v  e^-iϕ/2,
    0 ]^T
    ϕ_2^S_R   =[ 0,
    u  e^iϕ/2,
    0,
    v  e^-iϕ/2 ]^T
    ϕ_3^S_R   =[ v  e^iϕ/2,
    0,
    u  e^-iϕ/2,
    0 ]^T
    ϕ_4^S_R   =[ 0,
    v  e^iϕ/2,
    0,
    u  e^-iϕ/2 ]^T,

while the corresponding conjugated wave functions are

    ϕ̃_1^S_R   =
    [ 0,
    u  e^-iϕ/2,
    0,
    v  e^iϕ/2 ]^T
    ϕ̃_2^S_R   =
    [ u  e^-iϕ/2,
    0,
    v  e^iϕ/2,
    0 ]^T
    ϕ̃_3^S_R   =
    [ 0,
    v  e^-iϕ/2,
    0,
    u  e^iϕ/2 ]^T
    ϕ̃_4^S_R   =
    [ v  e^-iϕ/2,
    0,
    u  e^iϕ/2,
    0 ]^T .

 As before,  helicity conservation directly imposes b_i=d_i=b̃_i=d̃_i=0 
 and also q_1,3=r_1,3=p_2,4=s_2,4=q̃_1,3=r̃_1,3=p̃_2,4=s̃_2,4=0.



 §.§ Green's functions and pairing amplitudes



  §.§.§ Normal region

We construct the retarded Green's function from the scattering states in the same was as for SN 
and NSN junctions. In the normal region N the regular and anomalous Green's functions  become

    G_ee^r(x,x',ω)   =
    [ e^ik_e(x-x')
    M_1(x,x')                      0;                      0 e^-ik_e(x-x')M_2(x,x') ] ,
    
    G_hh^r(x,x',ω)   =
    [  e^ik_h(x-x')M_1(x',x)                      0;                      0 e^-ik_h(x-x')M_2(x',x) ] ,
    
    G_eh^r(x,x',ω)   =
    [  e^i(k_ex-k_hx') m_5                    0;                    0 e^-i(k_ex-k_hx') m_6 ] ,
    
    G_he^r(x,x',ω)   =
    [  e^i(k_hx-k_ex') m_7                    0;                    0 e^-i(k_hx-k_ex') m_8 ] ,

where M_1(x,x')=θ(x-x')m_1+θ(x'-x)m_3,  
M_2(x,x')=θ(x-x')m_2+θ(x'-x)m_4, 
m_1=u/vm_5, m_2=v/um_6, m_3=v/um_7, m_4=u/vm_8, and

    m_5   ≡ p_1r̃_4α_1=1/iv_fuv/u^2-v^2  e^i(k_e-k_h)L_ N-iϕ ,
    
    m_6   ≡ q_2s̃_3α_4=1/iv_fuv/u^2  e^i(k_h-k_e)L_ N-iϕ-v^2 ,
    
    m_7   ≡ r_1r̃_2α_1=1/iv_fuv/u^2  e^i(k_h-k_e)L_ N+iϕ-v^2 ,
    
    m_8   ≡ q_4s̃_1α_4=1/iv_fuv/u^2-v^2  e^i(k_e-k_h)L_ N+iϕ .

From these equations, we observe that  it is enough to specify the form of m_5. 
Thus, for energies within the superconducting gap previous expressions reduce to
    m_5(ω,L_ N,ϕ)   =-1/2v_F e^i(ϕ/2- L_ N/ξ_ω)/ sin[η(ω)-L_ N/ξ_ω+ϕ/2] 
    
      m_6(ω,L_ N,ϕ)   =m_5(ω,L_ N,-ϕ) e^i(2L_ N/ξ_ω+ϕ) ,
    
      m^*_7,8(-ω)   =m_6,5(ω) ,

which correspond to Eq. (<ref>) presented in the main text.


Instead of deriving the full expression for the LDOS from the regular part of G^r, we here focus on the main peaks as that is much easier and, for this discussion, as enlightening.
In the N region the low-energy peaks in the LDOS are from the discrete ABSs, whose number depend on the ratio L_ N/ξ. There exist different ways to calculate these energy levels. One efficient way is to locate the poles (zeros of the denominator) of the Andreev reflection coefficients, a_1,2, which are

    a_1(ω,L_ N,ϕ)=uv[ e^i(k_eL_ N-ϕ/2)- e^i(k_hL_ N
    +ϕ/2)]/u^2 e^i(k_hL_ N+ϕ/2)-v^2 e^i(k_eL_ N-ϕ/2) .

For energies within Δ, and using u/v= e^iη(ω), this expression can be written as

    a_1(ω,L_ N,ϕ)= sin[L_ N/ξ_ω-ϕ/2]/ sin[η(ω)-L_ N/ξ_ω+ϕ/2]

and also a_2(ω,L_ N,ϕ)=a_1(ω,L_ N,-ϕ). With these expressions, the poles of a_1,2 can be shown to give rise to the following condition for the ABSs, which is also given in the main text:

    2η(ω)-2L_ N/ξ_ω±ϕ=2π n ,   n=0,1,….

For a short junction, L_ N≪ξ, the condition for ABSs reduces to 2η(ω)±ϕ=2π n, with two Andreev levels at energies

    ω_±(ϕ)=±Δ cos(ϕ/2) ,

while for a long junction, L_ N≫ξ we obtain

    ω_±^n(ϕ)=v_F/2L_ N[2π(n+1/2)±ϕ] .



Notice that the Andreev bound states given by Eqs. (<ref>) and  (<ref>), 
and plotted in Fig. <ref>(a-b), develop crossings (indicated by red arrows) around zero energy 
ω=0 at π(2n-1) for n=1,2,3,…,  which are protected by time-reversal symmetry. 

The Andreev reflection magnitude |a_1(ω,L_ N,ϕ)| at ϕ=0,π is also plotted in Fig. <ref>(c-d) for short and long junctions. Observe that it develops resonances at the energies of the protected crossings. In particular, when ϕ=π a single zero-energy peak emerges in short junctions as shown in Fig. <ref>(c). In longer junctions the number of such resonances (that correspond to ABSs) increases, being proportional to L_ N/ξ, and the zero-energy peak coexists with additional peaks as we observe in Fig. <ref>(d). 
Remarkably, the same discussion was performed in the main text but there instead the ABSs energies was derived from the pairing amplitudes.



In terms of the pairing amplitudes in the N region we obtain after some algebra

    f_0, N^r(x,x',ω)   = e^i(k_ex-k_hx') m_5+ e^-i(k_ex-k_hx') m_6/2 ,
    
    f_3, N^r(x,x',ω)   = e^i(k_ex-k_hx') m_5- e^-i(k_ex-k_hx') m_6/2 ,

and f^r_1,N = f^r_2,N = 0, where m_5,6 are given by Eqs. (<ref>).

We can explicitly check the antisymmetry condition of these expressions.
For example, for the singlet component f_0, N^r(x,x',ω) we have

    f_0, N^a(x',x,-ω)   = e^-i(k_ex-k_hx') m_7^*(-ω)+ e^i(k_ex-k_hx') m_8^*(-ω)/2 ,
       = e^i(k_ex-k_hx') m_5+ e^-i(k_ex-k_hx') m_6/2 ,
       =f_0, N^r(x,x',ω) ,

and therefore we conclude that f_0, N^r(x,x',ω) is antisymmetric according to the relations in Eqs. (<ref>).
Likewise, we have verified the antisymmetry condition f^r_3,N(x,x',ω)=-f_3,N^a(x',x,-ω). 

By using Eq. (<ref>) we write down the even- and odd-frequency components

    f_0, N^r, E(x,x',ω)   =W_+(ω) cos[k_μ(x-x')] ,
    
    f_0, N^r, O(x,x',ω)   =W_-(ω)i sin[k_μ(x-x')] ,
    
    f_3, N^r, E(x,x',ω)   =W_+(ω)i sin[k_μ(x-x')] ,
    
    f_3, N^r, O(x,x',ω)   =W_-(ω) cos[k_μ(x-x')]

where W_±(ω)=[m_5 e^i(x+x')/ξ_ω± m_6 e^-i(x+x')/ξ_ω ]/2. 
These equations correspond to Eqs. (<ref>) presented in the main text. We have checked that Eqs. (<ref>) fulfill the antisymmetry conditions given by Eqs. (<ref>). 




  §.§.§ Superconducting regions

In the left superconducting region we obtain

    G^r_ee,↑↑(x,x',ω)   =Z
    [
    a_1(ω,L_ N,ϕ) e^i(k_h^Sx-k_e^Sx')+P(x,x',ω)] ,
    
    G^r_ee,↓↓(x,x',ω)   =Z
    [a_2(ω,L_ N,ϕ) e^i(k_h^Sx'-k_e^Sx)+P(x',x,ω)] ,
    
    G^r_hh,↓↓(x,x',ω)   =Z
    [
    a_1(ω,L_ N,ϕ) e^i(k_h^Sx-k_e^Sx')+Q(x,x',ω)]
    
    G^r_hh,↑↑(x,x',ω)   =Z
    [a_2(ω,L_ N,ϕ) e^i(k_h^Sx'-k_e^Sx)+Q(x',x,ω)] ,
    
    G^r_eh,↑↓(x,x',ω)   =Z
    [
    a_1(ω,L_ N,ϕ)v/u e^i(k_h^Sx-k_e^Sx')+P̅(x,x',ω)]
    
    G^r_eh,↓↑(x,x',ω)   =Z
    [a_2(ω,L_ N,ϕ)u/v e^i(k_h^Sx'-k_e^Sx)+P̅(x',x,ω)] ,
    
    G^r_he,↓↑(x,x',ω)   =Z
    [
    a_1(ω,L_ N,ϕ)u/v e^i(k_h^Sx-k_e^Sx')+Q̅(x,x',ω)] ,
    
    G^r_he,↑↓(x,x',ω)   =Z
    [a_2(ω,L_ N,ϕ)v/u e^i(k_h^Sx'-k_e^Sx)+Q̅(x',x,ω)] ,

where  


P(x,x',ω)   =θ(x-x')u/ve^ik_e^S(x-x')+θ(x'-x)v/ue^ik_h^S(x-x') ,

 Q(x,x',ω)   =θ(x-x')v/ue^ik_e^S(x-x')
+θ(x'-x)u/ve^ik_h^S(x-x') , 

P̅(x,x',ω)   =[θ(x-x')e^ik_e^S(x-x')
+θ(x'-x)e^ik_h^S(x-x')] ,

 Q̅(x,x',ω)   =[θ(x-x')e^ik_e^S(x-x')
+θ(x'-x)e^ik_h^S(x-x')] ,





with Z=1/iv_f1/(u/v)-(v/u)  and Z̅=1/-iv_f1/(u/v)^*-(v/u)^*.
In this case the Andreev reflection amplitudes obey  a_2(ω,L_ N,ϕ)=a_1(ω,L_ N,-ϕ) 
with a_1(ω,L_ N,ϕ) given by Eqs. (<ref>).
This results in the pairing amplitudes

    f^r_0, S_L(x,x',ω)   =
    Z/2 e^-ik(ω)(x+x')[
    a_1(ω,L_ N,ϕ)v/u e^ik_μ(x-x')+a_2(ω,L_ N,ϕ)u/v e^-ik_μ(x-x')] +Z e^ik(ω)|x-x'|C_xx' ,
    
    f^r_3, S_L(x,x',ω)   =
    Z/2 e^-ik(ω)(x+x')[
    a_1(ω,L_ N,ϕ)v/u e^ik_μ(x-x')-a_2(ω,L_ N,ϕ)u/v e^-ik_μ(x-x')] +Z e^ik(ω)|x-x'|iS_xx' ,

and f^r_1, S_L = f^r_2, S_L = 0.
In a similar way we proceed for the right superconducting region to calculate the pairing amplitudes.

We can then extract the odd- and even-frequency components. For energies within the superconducting gap we obtain

    f_0,±^r, E(x,x',ω)   =Z e^iϕ(1∓1)/2{ e^±κ(ω)(x+x')[a_1(4)(ω,L_ N,ϕ)a_1(ω)+a_2(3)(ω,L_ N,ϕ)a_1^*(ω)]+2 e^-κ(ω)|x-x'|}C_xx' ,
    
    f_0,±^r, O(x,x',ω)   =Z e^iϕ(1∓1)/2{ e^±κ(ω)(x+x')[a_1(4)(ω,L_ N,ϕ)a_1(ω)-a_2(3)(ω,L_ N,ϕ)a_1^*(ω)]}iS_xx' ,
    
    f_3,±^r, E(x,x',ω)   =Z e^iϕ(1∓1)/2{ e^±κ(ω)(x+x')[a_1(4)(ω,L_ N,ϕ)a_1(ω)+a_2(3)(ω,L_ N,ϕ)a_1^*(ω)]+2 e^-κ(ω)|x-x'|}iS_xx' ,
    
    f_3,±^r, O(x,x',ω)   =Z e^iϕ(1∓1)/2{ e^±κ(ω)(x+x')[a_1(4)(ω,L_ N,ϕ)a_1(ω)-a_2(3)(ω,L_ N,ϕ)a_1^*(ω)]}C_xx' ,

where ± subscripts correspond to results for the left and right superconducting regions, respectively.
Here we have used the notation a_1(ω)=v/u= e^-iη(ω) and
a_1(ω,L_ N,φ) given by Eq. (<ref>).
Moreover, we find that a_3,4(ω,L_ N,φ)=a_2,1(ω,L_ N,φ) e^2κ(ω).
As with NS junctions, we associate bulk behavior to elements in the pairing amplitudes that do not depend on the Andreev reflections and have an exponential form e^-κ(ω)|x-x'|. We therefore conclude that the second term in the curly brackets in the ESE and ETO pairing amplitudes emerge in the bulk of S. 
However, at the interface, we observe all different allowed pairing amplitudes. 
Notice that all the pairing amplitudes in the right superconductor acquire a phase factor e^iϕ.

