# Matrix-pairing states in the alkaline Fe-selenide superconductors: exotic
Josephson junctions

Emilian M. Nica  Department of Physics Box 871504 Arizona State University
Tempe, Arizona 85287-1504<EMAIL_ADDRESS>Qimiao Si  Department of Physics and
Astronomy, Rice University, 6100 Main St, Houston, TX, 77005, USA  Rice Center
for Quantum Materials, Rice University, 6100 Main St, Houston 77005 TX, USA
Onur Erten  Department of Physics Box 871504 Arizona State University Tempe,
Arizona 85287-1504

###### Abstract

True to their unconventional nature, multi-band alkaline Fe-selenides and,
more recently, the heavy-fermion CeCu2Si2 have shown experimental signatures
of fully-gapped but sign-changing superconductivity. A two-orbital pairing
state, called $s\tau_{3}$, with non-trivial matrix structure, was proposed as
a candidate able to reconcile the seemingly contradictory properties of these
superconductors. Motivated by the non-trivial orbital structure of the
proposed $s\tau_{3}$ state, which has opposite signs for the pairing functions
of the two orbitals, we study prototypical Josephson junctions where at least
one of the leads is in a superconducting state of this kind. An analysis of
these junctions in the limit of two degenerate orbitals (bands) and with a
simple form of junction hybridization reveals several remarkable properties.
One is the emergence of gapless, purely electron- and hole-like bound states
for $s\tau_{3}-N-s\tau_{3}$ junctions with arbitrary global phase difference
between the leads, and likewise for $s\tau_{3}-N-I$ junctions. The other is
the absence of static Josephson currents when both leads are superconducting.
In both of these signatures, $s\tau_{3}$ junctions are dramatically different
from more conventional Josephson junctions. We also find that the gapless
bound states are protected by an orbital-exchange symmetry, although the
protection is not topological. Junctions which break this symmetry, such as
$s\tau_{3}-N-s$, have gapped Andreev bound states. In general, the Josephson
effect also re-emerges once the degeneracy of the two orbitals is lifted. We
support these conclusions via analytical and numerical results for the bound
states, together with microscopic calculations of the Josephson current. Our
results indicate that junctions involving $s\tau_{3}$ pairing in alkaline Fe-
selenides will generically have bound states with a small gap together with a
greatly suppressed Josephson current.

##  I Introduction

Unconventional superconductors (SC’s) have led to some of the most important
questions in condensed matter physics. One such puzzle came to light not too
long ago in a branch of the Fe-based family, the alkaline Fe-selenides Lee
(2017). On one hand these SC’s are fully-gapped, as evidenced by ARPES studies
Mou  _et al._ (2011); X.-P. Wang and T. Qian and P. Richard and P. Zhang and
J. Dong and H.-D. Wang and C.-H. Dong and M.-H. Fang and H. Ding (2011); Xu
_et al._ (2012); Wang  _et al._ (2012), while on the other hand, they exhibit
an in-gap spin-resonance in inelastic neutron scattering experiments Park  _et
al._ (2011); Friemel  _et al._ (2012). The first of these features points
toward an $s$-wave pairing state, while the second implies a pairing state
which changes sign under a $\pi/2$ rotations Eschrig (2006); Stockert  _et
al._ (2011); Maier  _et al._ (2011); Dai (2015); Si  _et al._ (2016), such as
$d$-wave pairing. These seemingly mutually-exclusive, traits were shown to be
reconcilable in a pairing state which transforms as a sign-changing $B_{1g}$
representation of the point group, which nonetheless leads to a fully-gapped
SC state. We called this state $s\tau_{3}$ Nica  _et al._ (2017) since it
consists of a $s_{x^{2}y^{2}}$ form factor multiplied by a $\tau_{3}$ Pauli
matrix in a $d_{xz},d_{yz}$ two-orbital space appropriate to the alkaline Fe-
selendides. Underlying the remarkable properties of this pairing candidate is
the non-trivial $\tau_{3}$ matrix structure in orbital space, which ensures
that it transforms as $B_{1g}$. In the tight-binding appropriate to the
alkaline Fe-selenides this pairing state can also be thought of as an
effective $d+d$ intra- and inter-_band_ pairing Nica and Si (2021).
$s\tau_{3}$ pairing was stabilized in a realistic five-orbital model of the
alkaline Fe-selenides Nica  _et al._ (2017). Remarkably, a similar
experimental landscape has also recently emerged in the venerable heavy-
fermion CeCu2Si2, the first-discovered unconventional SC Steglich  _et al._
(1979). Believed to be a typical $d$-wave for almost it’s entire history,
CeCu2Si2 was recently found to exhibit a small gap for temperatures well below
$T_{c}$ from specific heat Kittaka  _et al._ (2014) and London penetration
depth measurements Pang  _et al._ (2018); Yamashita  _et al._ (2017). In-gap
spin resonances in the inelastic neutron spectrum Stockert  _et al._ (2011) as
well as other indicators single out CeCu2Si2 as another example of gapped but
sign-changing superconductivity. A pairing candidate, which incorporates the
orbital and spin-structure of this compound into a non-trivial matrix pairing
was also recently proposed by two of us Nica and Si (2021). More generally,
the $s\tau_{3}$ pairing state represents one of the most dramatic forms of the
overall notion of orbital-selective superconducting pairing introduced in the
context of 111 iron pnictides Yu  _et al._ (2014), which has also been
subsequently discussed in other iron pnictides Yin  _et al._ (2014); Ong  _et
al._ (2016) and the nematic FeSe Sprau  _et al._ (2017); Hu  _et al._ (2018);
Yu  _et al._ (2018).

Figure 1: Summary of our results for Josephson junctions along $x$ where at
least one of the Left (L) and Right (R) leads is in a $s\tau_{3}$ pairing
state. In that case, the pairing functions of the two orbitals have opposite
signs and thus have non-trivial structure on orbital space. We show the most
important results in the limit where the two orbitals are degenerate and
couple identically to a single orbital in the Center (C) metallic part. (a)
$s\tau_{3}-N-s\tau_{3}$ junctions where the L and R leads are identical up to
a global phase difference. For arbitrary global phase differences between L
and R leads, the bound states are either electron- or hole-like and become
gapless for a set of conserved momenta $k_{y}$ . In this limit, the static
Josephson current is zero, since the contributions of the two orbital sectors
effectively cancel. For almost degenerate bands, as expected for systems with
$s\tau_{3}$ pairing, the bound states acquire a small gap via Andreev
scattering, although exceptions can occur for junctions along $z$. A weak
static Josephson current, determined by band-splitting near the FS is also
expected. $s\tau_{3}-N-I$ junctions, where the R lead is in an insulating
state, has similar bound states. (b) $s\tau_{3}-N-s$ junction, where the R
lead consists of a single orbital in a trivial $s$-wave pairing state. In the
same degenerate limit, there are gapped Andreev bound states, but a vanishing
Josephson current. For small band splitting, a weak Josephson current is
expected.

In Ref. Nica and Si (2021), we also pointed out that the $s\tau_{3}$ pairing
state and a similarly-constructed microscopic candidate in CeCu2Si2 are both
equivalent to $d+d$ intra- and inter-band pairing in the band basis. This
equivalence holds when periodic boundary conditions are imposed along all of
the axes. Moreover, we discussed the similarities between $d+d$ and 3He-B,
where a similar $p+p$ pairing can be realized in spin space. In view of the
topologically-protected edge states of 3He-B, a natural question is whether
$d+d$ pairing can also have non-trivial edge states. This consideration opens
up a new direction in the exploration of both the general $d+d$ pairing and
actual microscopic states such as $s\tau_{3}$. For definiteness, in this work
we focus on Josephson junctions as these show features which are specific to
this type of matrix-pairing and therefore can facilitate its identification by
experiments. We reserve a more detailed discussion of the edge state spectrum
for future work.

Although some of the bulk properties of pairing candidates with non-trivial
matrix structure have already been mentioned, much remains to be explored as
far as characteristic experimental signatures are concerned. One promising
route involves Josephson junctions, which in principle can take advantage of
the inherent phase difference between the two orbital sectors in $s\tau_{3}$
pairing. Motivated by this observation, we study prototypical S-N-S and S-N-I
junctions where at least one of the leads is in a $s\tau_{3}$ pairing state.
Although our study is specifically geared toward junctions with $s\tau_{3}$
pairing states and thus target the alkaline Fe-selenides, we expect that some
of the features discussed here can also manifest for similar pairing
candidates with non-trivial matrix structure as proposed for CeCu2Si2 Nica and
Si (2021).

To isolate the most striking features of this pairing state, we describe the
superconducting leads via effective two-orbital models. Furthermore, in order
to capture the effects of the intrinsic phase difference between the pairing
functions of the two orbitals, we consider junctions where the two orbitals
couple to a single orbital in the metallic part. As we show, the bound states
in these rather unusual setups differ sharply from those in more conventional
junctions, which ignore cross-coupling between the various orbital. More
precisely, we examine setups where the Left (L) lead, Center (C) metallic, and
Right (R) leads are arranged in $s\tau_{3}-N-s\tau_{3}$, $s\tau_{3}-N-I$
(insulator), and $s\tau-N-s$ (trivial single-channel s-wave SC) junctions. We
assume that the coupling to the junction does not destabilize $s\tau_{3}$
pairing. Our salient results are as summarized below.

###  I.1 Summary of main results

Our most striking results occur when the two orbitals of the $s\tau_{3}$ SC
leads are degenerate. In this limit, both intra- and inter-orbital
hybridization terms which lead to band splitting are absent, and therefore the
two orbitals correspond to two degenerate bands. The precise definitions of
the intra- and inter-orbital hybridizations are discussed in Sec. II. In
addition, the two orbitals couple identically to a single orbital of the C
part. Due to the effects of the non-trivial orbital structure of $s\tau_{3}$
pairing and to the equal coupling to a single C orbital, the usual Andreev
scattering processes are effectively “frustrated”, leading to gapless bound
states which are either electron- or hole-like for $s\tau_{3}-N-s\tau_{3}$ for
arbitrary global phase differences between L and R leads. This also holds for
$s\tau_{3}-N-I$ junctions. The bound states are non-trivial since they decay
into the SC leads. We also find that the gapless bound states are also
protected by an orbital-exchange symmetry which is present for these junctions
in the limit considered here, although the protection is not topological in
origin. We find that breaking the orbital-exchange symmetry leads to a
restoration of Andreev scattering which mixes electron- and hole-like states,
which, in turn, leads to a gap. An extreme example of this symmetry breaking
is provided by $s\tau_{3}-N-s$ junctions, which exhibit Andreev bound states.
A second important feature is that the bound state spectra for the
$s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-s$ junctions are invariant to changes
in the _global_ phase difference across the junction. This is also in striking
contrast to the typical single-channel Josephson junction, where the Andreev
bound state spectrum changes with phase, leading to a static Josephson effect.
Our unusual findings are confirmed via microscopic calculations of the
Josephson currents in the tunneling limit. We show in this limit that the
contributions from the two sectors of $s\tau_{3}$ pairing, which have opposite
intrinsic phases, cancel.

We also consider more realistic junctions where the SC leads include intra-
and inter-orbital hybridization. We focus on cases where the band splitting in
the vicinity of the Fermi surface (FS) is much smaller than the pairing
amplitudes. We show that intra-orbital hybridization terms break the orbital-
exchange symmetry, leading to bound states more akin to the usual, gapped
Andreev bound states. The same terms also make the bound state spectrum
sensitive to a global phase difference across the junction. These findings are
further confirmed by microscopic calculation of the Josephson current in the
tunneling limit. As a counter-example to the gapped bound state spectrum, we
consider setups where the junction is along the $z$-direction, while the
intra- and inter-orbital hybridization are predominantly in-plane. We show
that gapless states can still be found for in-plane momenta along the
diagonals, where the orbital-exchange symmetry is preserved. In all cases, we
show that the analytical solutions are consistent with numerical results. Fig.
1 provides a summary of our main results.

As discussed in Refs. Nica  _et al._ (2017); Nica and Si (2021), $s\tau_{3}$
pairing is expected to provide a superconducting state whose single-particle
excitations are gapped everywhere on the FS. The simplest case to consider is
well established for the alkaline Fe-selenides. In this case, the splitting of
the relevant bands near the FS, which is essentially determined by intra- and
inter-orbital hybridization terms, is small when compared to the pairing
amplitudes. This is likely the case for a similar proposal in CeCu2Si2 Nica
and Si (2021). In light of our preceding discussion, the band-splitting in
alkaline Fe-selenides and CeCu2Si2 necessarily leads to gapped bound states
and to a finite Josephson current for arbitrary global phase differences
between the two leads. However, we expect that these effects are small.
Therefore, $s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junctions will
typically exhibit bound states with gaps which are much smaller than the bulk
gap, although exceptions with gapless states are also possible for arbitrary
global phase differences. Likewise, $s\tau_{3}-N-s\tau_{3}$ and
$s\tau_{3}-N-s$ junctions will exhibit static Josephson currents which are
determined by band-splitting near the FS, and are therefore small when
compared to similar junctions with orbitally-trivial pairing states. While
these predictions do not provide unequivocal probes for $s\tau_{3}$ pairing,
their combined signatures can provide significant supporting evidence.

The remainder of the article is divided into the following sections. In Sec.
II, we introduce the microscopic models for the junctions. In Sec. III, we
discuss both analytical and numerical solutions for $s\tau_{3}-N-s\tau_{3}$
junctions, and we briefly compare these to junctions with orbitally-trivial
pairing. $s\tau_{3}-N-I$ junctions have similar bound states, and are briefly
discussed at the end of the section. Section IV is devoted to $s\tau_{3}-N-s$
junctions. In the final Sec. V, we summarize our findings and discuss possible
experimental realizations of the junctions. The appendices contain detailed
discussions of some of the points presented in the main text. In Appendix A,
we present the analytical bound state solutions for $s\tau_{3}-N-s\tau_{3}$
junctions in the important limit of degenerate orbitals. Appendix B briefly
touches on the effect of unequal coupling of the two orbitals to the metallic
part, while Appendix C discusses junctions with orbitally-trivial pairings.
Appendices D and E present the analytical solutions for $s\tau_{3}-N-I$ and
$s\tau_{3}-N-s$ junctions, respectively, both in the limit of degenerate
orbitals. The final Appendix F, presents our results for the Josephson current
for $s\tau_{3}-N-s\tau_{3}$ junctions, as well as for junctions with
orbitally-trivial pairing.

##  II Models

The main distinction of the setups considered here are the superconducting
leads which are in a $s\tau_{3}$ pairing state, which was originally proposed
for an effective $d_{xz},d_{yz}$ two-orbital model for the alkaline Fe-
selenides Raghu  _et al._ (2008); Si and Abrahams (2008); Daghofer  _et al._
(2010). The _spin-singlet_ pairing state consists of a
$s_{x^{2}y^{2}}(k_{x},k_{y})$ form factor multiplied by a $\tau_{3}$ Pauli
matrix in orbital space. Due to the non-trivial $\tau_{3}$ orbital structure,
$s\tau_{3}$ pairing transforms as a $B_{1g}$ irreducible representation of the
tetragonal $D_{4h}$ point-group, and thus it changes sign under a $\pi/2$
rotation about $z$. In the simplest case, with the pairing amplitude being
larger than the band-splitting near the FS, which is due to intra- and inter-
orbital hybridization terms in the normal-state Hamiltonian, the Bogoliubov-de
Gennes (BdG) spectrum of this state is _fully gapped_. In the more general
cases, it is always gapped on the FS. In Ref. Nica  _et al._ (2017),
$s\tau_{3}$ pairing was stabilized in a realistic five-orbital model for the
alkaline Fe-selenides.

When considering Josephson junctions, the gapped nature of the $s\tau_{3}$
state implies that the relative phases of the pairing of the two
$d_{xz},d_{yz}$ orbital sectors are preserved by the junction. Furthermore,
the common $s_{x^{2}y^{2}}$ form factor has nodes along the $k_{x}=\pm\pi/2$
and $k_{y}=\pm\pi/2$ axes. For FS’s away from these axes, we can ignore the
momentum-dependence of the pairing, which does not introduce any qualitatively
new effects.

###  II.1 Junctions along $x$

For all of the three types of junctions along $x$, we adopt the following
model

$\displaystyle H=$ $\displaystyle
H_{\text{L}}+H_{\text{L-C}}+H_{\text{C}}+H_{\text{C-R}}+H_{\text{R}}.$ (1)

$H_{\text{L}}$, $H_{\text{R}}$, and $H_{\text{C}}$ are the bulk Hamiltonians
for the left, right leads and the center metallic parts, respectively.
$H_{\text{L-C}}$ and $H_{\text{C-R}}$ contain the terms at the lead-center
part interfaces, and by convention, are defined on the last and first sites of
the L and R leads, respectively. The L lead which includes the pairing terms
extends along $x[a]\leq(L_{x}-l)/2$, where $L_{x}$ is the number of sites for
the entire L-C-R system , while $l$ is the number of sites of the metallic C
part, along the $x$-direction, respectively. The length is defined in units of
the lattice spacing $a$. The most general form of the bulk L lead Hamiltonian
for $s\tau_{3}$ pairing is

$\displaystyle H_{\text{L}}=$ $\displaystyle H_{\text{L, TB}}+H_{\text{L,
Pair}},~{}\text{for}~{}x[a]<\left(\frac{L_{x}-l}{2}-1\right).$ (2)

$\displaystyle H_{\text{L, TB}}=$
$\displaystyle\sum_{\alpha\textbf{r}\sigma}\bigg{[}-\left(t_{x\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}},\alpha\sigma}+t_{y\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\alpha\sigma}+\text{H.c.}\right)-\mu
c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r},\alpha\sigma}$
$\displaystyle+$
$\displaystyle\sum_{\beta\neq\alpha}t_{4}\left(c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}}+\hat{\mathbf{y}},\beta\sigma}-c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}}-\hat{\mathbf{y}},\beta\sigma}+\text{H.c.}\right)\bigg{]}$
(3)

$\displaystyle H_{\text{L,
Pair}}=\sum_{\alpha\mathbf{r}}\Delta_{\alpha,\text{L}}\left(c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}-c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}\right)+\text{H.c.}$
(4)

The indices $\alpha\in\\{1,2\\}$ stand for the $d_{xz}$ and $d_{yz}$ orbitals,
respectively, while the spin is represented by
$\sigma\in\\{\uparrow,\downarrow\\}$. $t_{x\alpha}$ are the nearest-neighbor
(NN) hopping coefficients for orbital $\alpha$, while $t_{4}$ is a NN orbital
hybridization, with all coefficients taken to be real. These are based on a
simplified two-orbital model introduced in Ref. Raghu  _et al._ (2008). Note
the absence of terms proportional to a $\tau_{2}$ matrix in the tight-binding
part, since these break time-reversal symmetry. The pairing satisfies
$\Delta_{1}=-\Delta_{2}=\Delta$. As mentioned previously, we neglect the
spatial dependence of the pairing which is not of immediate importance to the
effects considered here.

For periodic boundary conditions (BC’s), the BdG Hamiltonian for the L lead
and for single spin sector reads

$\displaystyle
H_{\text{L}}(\mathbf{k})=\begin{pmatrix}\xi_{0}+\xi_{3}&\xi_{1}&\Delta&0\\\
\xi_{1}&\xi_{0}-\xi_{3}&0&-\Delta\\\
\Delta^{*}&0&-(\xi_{0}+\xi_{3})&-\xi_{1}\\\
0&-\Delta^{*}&-\xi_{1}&-(\xi_{0}-\xi_{3})\end{pmatrix}$ (5)

where

$\displaystyle\xi_{0}(\mathbf{k})=$
$\displaystyle(t_{1}+t_{2})\left[\cos(k_{x}a)+\cos(k_{y}a)\right]-\mu$ (6)
$\displaystyle\xi_{3}(\mathbf{k})=$
$\displaystyle(t_{1}-t_{2})\left[\cos(k_{x}a)-\cos(k_{y}a)\right]$ (7)
$\displaystyle\xi_{1}(\mathbf{k})=$
$\displaystyle-4t_{4}\sin(k_{x}a)\sin(k_{y}a).$ (8)

$D_{4h}$ symmetry restricts

$\displaystyle t_{x1}=$ $\displaystyle t_{y2}=t_{1}$ (9) $\displaystyle
t_{x2}=$ $\displaystyle t_{y1}=t_{2}.$ (10)

Three terms $\xi_{0}\tau_{0},\xi_{1}\tau_{1}$, and $\xi_{3}\tau_{3}$ determine
the normal state, corresponding to a common dispersion, inter- and intra-
orbital hybridization terms, respectively. Note that terms
$\propto\cos(k_{x}a)\cos(k_{y}a)$, which preserve the lattice symmetry can
also be added to the orbital-_diagonal_ $\xi_{0}$ terms Raghu  _et al._
(2008). While these can lead to a change in shape of the FS, they preserve the
orbital structure of $H_{\text{L}}$, i.e. do not induce any additional band
splitting, which is essential to the results of this work. Therefore, we
ignore these additional contributions. The important degenerate-orbital limit,
which we mention in the following, occurs when both inter- and intra-orbital
hybridization terms are set to zero with $\xi_{1}=\xi_{3}=0$. For junctions
along $x$ which break the translation symmetry along this direction, this
corresponds to $t_{1}-t_{2}=t_{4}=0.$

The positive eigenvalues are

$\displaystyle
E_{1,2}=\sqrt{\xi^{2}_{0}+\xi^{2}_{3}+\xi^{2}_{1}+|\Delta|^{2}\pm
2\sqrt{\xi^{2}_{0}\left(\xi^{2}_{1}+\xi^{2}_{3}\right)+\xi^{2}_{1}|\Delta|^{2}}}.$
(11)

As mentioned previously, we consider cases where the band splittings near
either of the two FS’s are smaller than the pairing amplitude

$\displaystyle\sqrt{\xi^{2}_{1}(\mathbf{k})+\xi^{2}_{3}(\mathbf{k})}\big{|}_{\mathbf{k}\in\text{FS}}<|\Delta|^{2}.$
(12)

which ensures that the leads are gapped in the bulk.

The Hamiltonian at the L-C interface is

$\displaystyle H_{\text{L-C}}=$ $\displaystyle\sum_{\alpha
y}\bigg{\\{}\sum_{\sigma}\left[-\left(V_{\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}},\sigma}+t_{y\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\alpha\sigma}+\text{H.c.}\right)-\mu
c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r},\alpha\sigma}\right]+\left[\Delta_{\alpha}\left(c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}-c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}\right)+\text{H.c.}\right]\bigg{\\}}$
$\displaystyle,~{}$ $\displaystyle\text{for}~{}x[a]=\frac{L_{x}-l}{2}$ (13)

As mentioned earlier, the $k$-dependence of the $s(k)$ factor in the
$s\tau_{3}$ pairing is unimportant for our purpose, and this is reflected in
the form of the pairing term of Eq. 13. Crucially, the terms proportional to
$V_{\alpha}$ denote the hybridization of the two orbitals in the L lead to _a
single orbital_ in the C part along $x$. Since the point-group symmetry is
necessarily broken at the lead-center interfaces, this symmetry does not _a
priori_ restrict the values of the $V_{\alpha}$’s. In the following, we allow
these couplings to take arbitrary values and we comment on the effects of the
$d_{xz},d_{yz}$ nature of the orbitals in the leads in Sec. V. The remaining
terms are identical to those of $H_{\text{L}}$ while the summation is over the
$y$ coordinate along the junction. We typically choose $V_{1}=V_{2}=t_{1}=1$,
unless stated otherwise. As already mentioned, we neglect the spatial
dependence of the pairing.

The Hamiltonian for the bulk of the central (C) part reads

$\displaystyle H_{\text{C}}=$
$\displaystyle\sum_{\alpha\mathbf{r}\sigma}\left[-t_{1}\left(c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}+\hat{\mathbf{x}},\sigma}+c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\sigma}+\text{H.c.}\right)-\mu
c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r},\sigma}\right],~{}\text{for}~{}\frac{L_{x}-l}{2}+1\leq
x[a]\leq\frac{L_{x}+l+2}{2}-1.$ (14)

As mentioned previously, the C part involves a _single channel_ without any
pairing. We consider for simplicity nearest-neighbor hopping in this sector,
taken to be equal to $t_{1}$.

In the case of $s\tau_{3}-N-s\tau_{3}$ pairing we choose the Hamiltonians at
the C-R junction $H_{\text{C-R}}$ for $x[a]=(L_{x}+l+2)/2$ and for the bulk of
the R lead $H_{\text{R}}$ for $x[a]\geq(L_{x}+l+2)/2+1$ to have identical form
and coefficients to $H_{\text{L-C}}$ and $H_{\text{L}}$, respectively, with
the exception of a possible _global non-zero phase_ $\phi$ for the pairing
terms:

$\displaystyle\Delta_{\text{L}}=$ $\displaystyle\Delta$ (15)
$\displaystyle\Delta_{\text{R}}=$ $\displaystyle\Delta e^{i\phi}.$ (16)

For $s\tau_{3}-N-s$, $H_{\text{C-R}}$ and $H_{\text{R}}$ reflect the presence
of a single channel in the R lead as

$\displaystyle H_{\text{C-R}}=$ $\displaystyle\sum_{\alpha
y}\bigg{\\{}\sum_{\sigma}\bigg{[}-\left(Vc^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}-\hat{\mathbf{x}},\sigma}+t_{y}c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\sigma}+\text{H.c.}\right)$
$\displaystyle-$ $\displaystyle\mu
c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r},\alpha\sigma}\bigg{]}+\left[\Delta
e^{i\phi}\left(c^{{\dagger}}_{\mathbf{r},\uparrow}c^{{\dagger}}_{\mathbf{r},\downarrow}-c^{{\dagger}}_{\mathbf{r},\downarrow}c^{{\dagger}}_{\mathbf{r},\uparrow}\right)+\text{H.c.}\right]\bigg{\\}}$
$\displaystyle,~{}\text{for}~{}x[a]=\frac{L_{J}+l+2}{2}$ (17)

$\displaystyle H_{\text{R, TB}}=$
$\displaystyle-\sum_{\textbf{r}\sigma}\bigg{[}-t\left(c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}-\hat{\mathbf{x}},\sigma}+c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\sigma}+\text{H.c.}\right)$
$\displaystyle-$ $\displaystyle\mu
c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r},\sigma}\bigg{]}$ (18)

$\displaystyle H_{\text{R, Pair}}=\sum_{\mathbf{r}}\Delta
e^{i\phi}\left(c^{{\dagger}}_{\mathbf{r},\uparrow}c^{{\dagger}}_{\mathbf{r},\downarrow}-c^{{\dagger}}_{\mathbf{r},\downarrow}c^{{\dagger}}_{\mathbf{r},\uparrow}\right)+\text{H.c.}.$
(19)

For simplicity, we always consider $V=t=1$.

###  II.2 Junctions along $z$

For $s\tau_{3}-N-s\tau_{3}$ junctions along $z$ we define the complete
Hamiltonian via

$\displaystyle\tilde{H}=$
$\displaystyle\tilde{H}_{\text{L}}+\tilde{H}_{\text{L-C}}+\tilde{H}_{\text{C}}+\tilde{H}_{\text{C-R}}+\tilde{H}_{\text{R}}.$
(20)

In addition to the terms already present in the pure two-dimensional case, we
consider additional NN hopping along $z$ for the L/R leads as

$\displaystyle\tilde{H}_{\text{L, TB}}=$ $\displaystyle H_{\text{L,
TB}}-\sum_{\alpha\textbf{r}\sigma}t_{z}\left(c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{z}},\alpha\sigma}+\text{H.c.}\right),$
(21)

where $H_{\text{L, TB}}$ was defined in Eq. 3, and the sum over lattice sites
is now along all of the three axes. We take the lengths of the L, C, and R
parts along $z$ to be the same as for the junctions along $x$. For simplicity,
we consider only orbital-diagonal, NN hopping along the $z$-direction.
Additional terms are of course possible, but they do not qualitatively change
our conclusions. For clarity, we ignore the momentum dependence of the two-
dimensional pairing, as for junctions along $x$. Similar expressions hold for
the bulk of the R lead, with the exception of a global phase difference
$\phi$. Likewise, the C part amounts to

$\displaystyle\tilde{H}_{\text{C}}=$ $\displaystyle
H_{\text{C}}-\sum_{\alpha\textbf{r}}t_{z}\left(c^{{\dagger}}_{\mathbf{r},\sigma}c_{\mathbf{r}+\hat{\mathbf{z}},\sigma}+\text{H.c.}\right).$
(22)

We use a NN hopping in the C part which is equal in amplitude to that of the
L/R leads $t_{z}=0.2t_{1}$, unless stated otherwise. The reduced value of
$t_{z}$ reflects the tetragonal anisotropy of our target systems.

Finally, the Hamiltonians for the L-C and C-R interfaces are obtained via a
straightforward generalization to hopping along $z$

$\displaystyle H_{\text{L-C}}=$
$\displaystyle-\sum_{\alpha,xy}\bigg{\\{}\sum_{\sigma}\bigg{[}\left(-V_{\alpha
z}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{z}},\sigma}-t_{x\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}},\alpha\sigma}-t_{y\alpha}c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{y}},\alpha\sigma}+\text{H.c.}\right)-\mu
c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r},\alpha\sigma}$
$\displaystyle+$
$\displaystyle\sum_{\beta\neq\alpha}t_{4}\left(c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}}+\hat{\mathbf{y}},\beta\sigma}-c^{{\dagger}}_{\mathbf{r},\alpha\sigma}c_{\mathbf{r}+\hat{\mathbf{x}}-\hat{\mathbf{y}},\beta\sigma}+\text{H.c.}\right)\bigg{]}+\left[\Delta_{\alpha}\left(c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}-c^{{\dagger}}_{\mathbf{r},\alpha\downarrow}c^{{\dagger}}_{\mathbf{r},\alpha\uparrow}\right)+\text{H.c.}\right]\Bigg{\\}}$
(23)

$\tilde{H}_{\text{C-R}}$ can be obtained via a similar generalization of
$H_{\text{C-R}}$.

Unless stated otherwise, we choose units where the nearest-neighbor hopping
and hybridization at the interfaces are $t_{1}=V_{1}=1$ and $\Delta=0.4$.

###  II.3 Orbital-exchange symmetry

We consider the $s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ setups in the
important limit where the two orbitals have identical tight-binding
coefficients, zero intra-orbital hybridization ($\xi_{3}=0$) and inter-orbital
hybridization ($\xi_{1}=0$), and couple identically to the single orbital of
the C part. This limit corresponds to

$\displaystyle t_{1}=$ $\displaystyle t_{2}$ (24) $\displaystyle t_{4}=$
$\displaystyle 0$ (25) $\displaystyle V_{1}=$ $\displaystyle V_{2}.$ (26)

The two orbitals in the L (or R if appropriate) lead correspond to degenerate
bands in the normal state. We define a _local_ transformation which acts on
the L (and R) spinors
$\Psi^{T}=(c_{\mathbf{r},1\uparrow},c_{\mathbf{r},2\uparrow},c^{{\dagger}}_{\mathbf{r},1\downarrow},c^{{\dagger}}_{\mathbf{r},2\downarrow})$
as

$\displaystyle\Psi\rightarrow\hat{R}\Psi$ (27)

where

$\displaystyle\hat{R}=\begin{pmatrix}0&i&0&0\\\ i&0&0&0\\\ 0&0&0&-i\\\
0&0&-i&0\end{pmatrix}.$ (28)

This corresponds to an exchange of the two orbitals followed by a gauge
transformation which changes the sign of the pairing $\Delta$. Similarly, in
the C part the transformation acting on the single-orbital spinor
$\psi^{T}=(c_{\mathbf{r},\uparrow},c^{{\dagger}}_{\mathbf{r},\downarrow})$

$\displaystyle\psi\rightarrow\hat{S}\psi$ (29)

where

$\displaystyle\hat{S}=\begin{pmatrix}i&0\\\ 0&-i\\\ \end{pmatrix}.$ (30)

This operation, corresponding to a gauge transformation for the C part, is
required to compensate for the factors of $i$ in $H_{\text{L-C}}$ and
$H_{\text{C-R}}$ if appropriate. Both $\hat{R}$ and $\hat{S}$ are anti-
hermitian operators which obey

$\displaystyle\hat{R}^{{\dagger}}=$ $\displaystyle-\hat{R}$ (31)
$\displaystyle\hat{R}^{2}=$ $\displaystyle-\hat{1}$ (32)

and similarly for $\hat{S}$.

It is straightforward to define an orbital-exchange operator for the entire
lattice models, corresponding to each of the junctions. For
$s\tau_{3}-N-s\tau_{3}$ junctions it reads

$\displaystyle\hat{Q}=\begin{cases}\hat{R},~{}\forall~{}x\in\text{L, R}\\\
\hat{S},~{}\forall~{}x\in\text{C}\end{cases}$ (33)

The $s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junction Hamiltonians in the
degenerate-orbital limit considered here are invariant under $\hat{Q}$. This
orbital exchange symmetry plays an important role in classifying the electron-
like and hole-like solutions encountered in these cases.

##  III $s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junctions

We treat $s\tau_{3}-N-s\tau_{3}$ junctions in some detail in the following and
only briefly discuss $s\tau_{3}-N-I$ junctions at the end of the section and
in Appendix C, since the latter exhibit similar bound states.

###  III.1 Degenerate orbitals

We consider a $s\tau_{3}-N-s\tau_{3}$ junction as defined previously, in the
limit of degenerate orbitals ( Eqs. 24-26 ). We find the bound states in the
continuum limit by linearizing the BdG equations in the vicinity of points
$(\alpha K_{Fx},K_{Fy})$ on the FS Sauls (2018), where $\alpha=\pm 1$. The
detailed solution is presented in Appendix A. Here, we summarize some of the
most important results. In contrast to the typical single-channel junction,
the presence of two orbitals in either leads, which couple to a single orbital
in the C part, imposes additional BC’s. Indeed, in the limit considered here,
only the linear combination $c_{1}+c_{2}$ of the two orbitals couples to the C
part. The corresponding BdG coefficients must vary continuously at the L-C or
C-R interfaces, much like in a single-channel junction. By contrast, the anti-
symmetric combination does not couple to the C part and thus the corresponding
BdG coefficients obey open BC’s at the interface. As shown in Appendix A, the
two channels, corresponding to symmetric and anti-symmetric linear
combinations, do not decouple in the bulk of the leads, due to the non-trivial
matrix structure of $s\tau_{3}$ pairing. Consequently, the continuity and open
BC’s cannot be satisfied simultaneously for solutions which involve both
electron- and hole-like solutions in the C part. The solutions are either
electron- or hole-like in the C part, in contrast to the typical Andreev bound
states which mix the two.

The electron-like bound states with eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon
l}{2v_{Fx}}+\frac{K_{Fx}l}{2}+\frac{(a-m)\pi}{2}\right).$ (34)

have the form

$\displaystyle\Psi^{L/R}_{K_{Fy};\text{Electron}}=$
$\displaystyle\left(u^{L/R}_{1e},u^{L/R}_{2e},v^{L/R}_{1e},v^{L/R}_{2e}\right)^{T}$
(35)

$\displaystyle u^{L/R}_{1e}=$
$\displaystyle\frac{1}{\sqrt{2}}|A||\Delta|e^{i\theta^{0}_{A}}\sin\left[K_{Fx}x\mp\frac{\epsilon
l}{2v_{Fx}}+\frac{(a+m)\pi}{2}\right]$ (36) $\displaystyle u^{L/R}_{2e}=$
$\displaystyle u^{L/R}_{1}$ (37) $\displaystyle v^{L/R}_{1e}=$
$\displaystyle\frac{1}{\sqrt{2}}|A||\Delta|e^{i\theta^{0}_{A}}\sin\left[K_{Fx}\left(x\pm\frac{l}{2}\right)+\gamma_{L/R}\pi\right]$
(38) $\displaystyle v^{L/R}_{2e}=$ $\displaystyle-v^{L/R}_{1}$ (39)
$\displaystyle\gamma_{L}=$ $\displaystyle a$ (40) $\displaystyle\gamma_{R}=$
$\displaystyle m$ (41)

in the L/R leads and

$\displaystyle\Psi^{C}_{K_{Fy};\text{Electron}}=$
$\displaystyle\left(u^{C}_{e},v^{C}_{e}\right)^{T}$ (42)

$\displaystyle u^{C}_{e}=$ $\displaystyle
2|A||\Delta|e^{i\theta^{0}_{A}}\sin\left[K_{Fx}x+\frac{\epsilon
x}{v_{Fx}}+\frac{(a+m)\pi}{2}\right]$ (43) $\displaystyle v^{C}_{e}=$
$\displaystyle 0.$ (44)

in the C part. $|A|$ is a normalization constant, $\theta^{0}_{A}$ is an
arbitrary phase, and $a,m$ are arbitrary integers. Crucially, the C $v$ BdG
coefficient is identically zero.

Similarly, the hole-like solutions with eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon
l}{2v_{Fx}}-\frac{K_{Fx}l}{2}+\frac{(n-b)\pi}{2}\right).$ (45)

are of the form

$\displaystyle\Psi^{L/R}_{K_{Fy};\text{Hole}}=$
$\displaystyle\left(u^{L/R}_{1h},u^{L/R}_{2h},v^{L/R}_{1h},v^{L/R}_{2h}\right)^{T}$
(46)

$\displaystyle u^{L/R}_{1h}=$
$\displaystyle\frac{1}{\sqrt{2}}|B|e^{i\theta^{B}_{0}}\sin\left[K_{Fx}\left(x\pm\frac{l}{2}\right)+\gamma_{L/R}\pi\right]$
(47) $\displaystyle u^{L/R}_{2h}=$ $\displaystyle-u^{L}_{1h}$ (48)
$\displaystyle v^{L/R}_{1h}=$
$\displaystyle\frac{1}{\sqrt{2}}|B|e^{i\theta^{B}_{0}}\sin\left[K_{Fx}x\pm\frac{\epsilon
l}{2v_{Fx}}+\frac{\pi}{2}+\frac{(n+b)\pi}{2}\right]$ (49) $\displaystyle
v^{L/R}_{2h}=$ $\displaystyle v^{L}_{1h}$ (50) $\displaystyle\gamma_{L}=$
$\displaystyle b$ (51) $\displaystyle\gamma_{R}=$ $\displaystyle n$ (52)

for the L/R leads and

$\displaystyle\Psi^{C}_{K_{Fy};\text{Hole}}=$
$\displaystyle\left(u^{C}_{h},v^{C}_{h}\right)^{T}$ (53)

$\displaystyle u^{C}_{h}=$ $\displaystyle 0$ (54) $\displaystyle v^{C}_{h}=$
$\displaystyle\sin\left[K_{Fx}x-\frac{\epsilon
x}{v_{Fx}}+\frac{(n+b)\pi}{2}\right]$ (55)

in the C part. As before, $|B|$ is a normalization constant, $\theta^{0}_{B}$
is an arbitrary phase, while $b,n$ are integers. In contrast to the electron-
like solutions, $u^{C}=0$.

These detailed analytical solutions present a number of important features.
Firstly, the solutions are either electron- or hole-like since the
corresponding BdG coefficient $v$ or $u$ vanishes in the C part. This is in
clear contrast to a single-channel junction, where electron-like and hole-like
solutions are mixed via Andreev scattering. The presence of both open and
continuity BC’s at L-C and C-R interfaces, as well as the non-trivial
structure of the pairing, prevents the existence of solutions which mix
electron and hole states.

Secondly, electron- and hole-like solutions become gapless when

$\displaystyle K_{Fx}l=$ $\displaystyle(m-a)\pi$ (56) $\displaystyle K_{Fx}l=$
$\displaystyle(n-b)\pi,$ (57)

respectively, as indicated by Eqs. 34 and 45. This happens when $K_{Fx}l$ is
either 0 or a multiple of $\pi$. Recall that the solutions are labeled by the
momenta $(K_{Fx},K_{Fy})$ which vary continuously along the FS. Consequently,
the gapless conditions can be realized in multiple instances along the FS if
the latter and/or junction length $l$ are sufficiently large. Linearizing Eqs.
34 and 45 about these points indicates that the electron- and hole-like states
are counter-propagating, as a function of the conserved momentum $K_{Fy}$.

Thirdly, the gaplessness of the electron- and hole-like states is protected by
the orbital-exchange symmetry defined in Sec. II.3. Indeed, it is
straightforward to check that

$\displaystyle\hat{R}\Psi^{L/R}_{K_{Fy};\text{Electron}}=$ $\displaystyle
i\Psi^{L/R}_{K_{Fy};\text{Electron}}$ (58)
$\displaystyle\hat{R}\Psi^{L/R}_{K_{Fy};\text{Hole}}=$ $\displaystyle
i\Psi^{L/R}_{K_{Fy};\text{Hole}},$ (59)

together with

$\displaystyle\Psi^{L/R,{\dagger}}_{K_{Fy};\text{Electron}}\hat{R}=$
$\displaystyle i\Psi^{L/R,{\dagger}}_{K_{Fy};\text{Electron}}.$ (60)

The last relation follows from the anti-hermitian nature of $\hat{R}$. Similar
relations hold for $\hat{S}$ acting on the C spinors. Together, these imply
that any operator $\hat{O}$ added to the Hamiltonian of the junction, which
commutes with $\hat{R}$ and $\hat{S}$, will not mix electron- and hole-like
states via

$\displaystyle\braket{\Psi^{{\dagger}}_{\text{Electron}}}{\hat{O}\hat{R}}{\Psi_{\text{Hole}}}=\braket{\Psi^{{\dagger}}_{\text{Electron}}}{\hat{R}\hat{O}}{\Psi_{\text{Hole}}}$
(61)

which implies that

$\displaystyle(-i)\braket{\Psi^{{\dagger}}_{\text{Electron}}}{\hat{O}}{\Psi_{\text{Hole}}}=$
$\displaystyle(+i)\braket{\Psi^{{\dagger}}_{\text{Electron}}}{\hat{O}}{\Psi_{\text{Hole}}}$
(62)

and subsequently that

$\displaystyle\braket{\Psi^{{\dagger}}_{\text{Electron}}}{\hat{O}}{\Psi_{\text{Hole}}}=0.$
(63)

Note that $\hat{O}$ is not necessarily local at the lattice level. A similar
reasoning can be applied to $\hat{S}$. Therefore, the orbital-exchange
symmetry protects the gapless electron- and hole-like states, in analogy to
the time-reversal operator for spin-polarized edge states in quantum spin Hall
systems Kane and Mele (2005), although we stress that the gapless states are
not due to any topological property of the bulk in our cases. In practice, the
orbital-exchange symmetry is preserved when no intra-orbital hybridization
terms (such as $\xi_{3}$ defined in Eq. 5) are present, although inter-orbital
hybridization terms (such as $\xi_{1}$ in Eq. 5) preserve the orbital-exchange
symmetry.

Fourthly, the bound state spectrum is independent of the global relative phase
$\phi$ between L and R leads, as indicated by Eqs. 34 and 45. This is in clear
contrast to the single-channel junctions, where the change in the Andreev
bound state spectrum with relative phase is proportional to the static
Josephson current Sauls (2018). In the $s\tau_{3}-N-s\tau_{3}$ case, a bound
state spectrum which is insensitive to the relative phase indicates the
absence of a static Josephson current. This is confirmed by microscopic
calculations of the Josephson current in the tunneling approximation,
discussed in detail in Appendix F, in the degenerate orbital limit considered
here. As shown there, the Josephson currents due to each orbital sector cancel
due essentially to the $\pi$ phase difference between the two components of
$s\tau_{3}$ pairing.

Next, we proceed to confirm the analytical solutions via numerical solutions
of the lattice model introduced in Sec. II, in the degenerate-orbital limit.
The results were obtained for a lattice of 100 sites along $x$ with a C part
of 5 sites extending along $49\leq x\leq 53$, in units of the lattice constant
$a$.

In Fig. 2, we show the bound state spectrum as a function of the conserved
momentum $k_{y}$, for varying chemical potential $\mu$ and global relative
phase difference between the two superconducting leads $\phi$. For $\mu=-3.0$
in panel (a), we illustrate that the spectrum is independent of $\phi$, in
accordance with the eigenvalues in Eqs. 34 and 45. This implies the absence of
static Josephson current, as confirmed in Appendix F by a microscopic
calculation of the Josephson current in the tunneling limit. Also note that
the bound states become gapless and cross twice near $k_{y}a\approx 0.3\pi$.
This is consistent with the condition for gapless states as determined by the
analytical solution (Eqs. 56, 57). Indeed, for increasing chemical potential
$\mu=-2.0$, as shown in panel (b), an additional crossing occurs near
$k_{y}=0$. This is due to an $K_{Fx}$ (near $K_{Fy}=0$) which increases by
$\pi$ as the FS expands. The crossings in (a) likewise shift to higher
$k_{y}$. A similar picture presents itself with increasing chemical potential
$\mu=-0.25$ in panel (c). Together, these results are consistent with the
analytical solution.

We also consider the nature of the bound states. In Fig. 3, we show a close-up
of panel (a) of Fig. 2 for $\phi=0$. One pair of hole- and electron-like
states in the vicinity of a crossing are marked by a square and triangle,
respectively. To elucidate the nature of these states, in Fig. 4 we illustrate
the real parts of the BdG coefficients of the hole-like state marked in Fig. 3
by a square. Note that the imaginary pars are trivially zero and are not
shown. The $u$ BdG coefficients shown in panel (a) are consistent with the
analytical solution in Eq. 46. Indeed, $u^{C}$ vanishes in the C part for
$49\leq x[a]\leq 53$, as expected for a hole-like solution. Furthermore,
$u^{L/R}_{1}$ and $u^{L/R}_{2}$ both vanish at the L-C and C-R interfaces,
respectively, and their signs are consistent with Eq. 46 for $n-b$ an odd
multiple of $\pi$. In panel (b), the $v^{C}$ is finite and continuous across
the junction. We conclude that the numerical solution is consistent with the
hole-like solution indicated by the analytical results.

A similar picture emerges for the electron-like state marked with a triangle
in Fig. 3. In Fig. 5, we show the real parts of the BdG coefficients of this
state across the junction, and compare these with the analytical solution in
Eq. 35. $u^{C}$ is now finite, while $v^{C}$ vanishes in the C region, as
expected. Furthermore, the signs of $u^{L/R}_{1}$ and $u^{L/R}_{2}$ are
consistent with a solution with $m-a$ an odd multiple of $\pi$. Together with
the BdG coefficient for the hole-like state, these numerical solutions confirm
the analytical results.

We briefly discuss the case for a $s\tau_{3}-N-s\tau_{3}$ junction along $z$
for the model introduced in Sec. II.2. The analytical solution is similar to
the junction along $x$ with the exception that both $k_{x},k_{y}$ are
conserved quantities. We therefore linearize about the pair of points
$(K_{Fx},K_{Fy},\alpha K_{Fz})$ and obtain the bound-state spectrum as in Eqs.
34 and 45 which now extend along the $z$-direction. The eigenvalues are
obtained via the replacement $K_{Fx}\rightarrow K_{Fz}.$

We comment on the effects of unequal coupling of the two orbitals in either
lead to the C part. In this case, both symmetric and anti-symmetric linear
combinations of the two orbitals (Appendix A), couple to the C part, in
contrast to the solutions presented here. As shown by numerical results
presented in Appendix B, unequal couplings induce a gap, as well as a spectrum
which depends on the global phase difference $\phi$. The first can be
understood via a breaking of the orbital-exchange symmetry, while the second
can be confirmed via a direct calculation of the Josephson current, as
discussed in Appendix F.

Finally, we note that the surprising behavior for $s\tau_{3}$ junctions
coupled to a single orbital in the C part differs dramatically from that of
junctions consisting of leads in a pairing state which has trivial orbital
structure. We consider junctions with such a state, which we call $s\tau_{0}$,
where the two orbitals have identical pairing functions, including the signs.
For degenerate orbitals, we can apply the same analysis as in the
$s\tau_{3}-N-s\tau_{3}$ case (Appendix A) and we again find that the symmetric
linear combination of the two orbitals in either leads couples across the
junction. However, the anti-symmetric linear combination is entirely decoupled
along the entire junction, and in contrast to $s\tau_{3}$ junctions, imposes
no additional constraints on the bound-state solutions. Thus $s\tau_{0}$
junctions for degenerate orbitals behave essentially like conventional single-
channel junctions. This is illustrated in Appendix C.

Figure 2: Bound state spectrum of a $s\tau_{3}-N-s\tau_{3}$ junction along
$x$, as defined in Sec. II.1, with degenerate orbitals in the L and R leads
which also couple identically to a single orbital in the N part. (a)
Eigenvalues for a fixed chemical potential $\mu=-3.0$ and varying global
relative phase $\phi$ between the R and L superconducting leads, as functions
of the conserved momentum $k_{y}$. Note that the spectrum is insensitive to
$\phi$, in accordance to the analytical solution in Eqs. 34 and 45. Also note
the presence of gapless states near $k_{y}a\approx 0.3\pi$. (b) Same as (a)
for fixed $\mu=-2.0$ and $\phi=0.0$. Although not shown, the spectrum is also
insensitive to $\phi$. There is an additional crossing of gapless states near
$k_{y}a=0.0$. This is consistent with the condition for gapless states
obtained from the analytical solutions in Eqs. 56-57. With decreasing $\mu$
the FS increases, allowing $K_{Fx}l$ to reach values which are multiples of
$\pi$ near $k_{y}=0$. (c) Same as (b) for $\mu=-0.25$ and $\phi=0$. Note that
the number of gapless crossings increases with the size of the FS, as
predicted by the analytical solution. Figure 3: Close-up of panel (a) of Fig.
2 (a). The square and triangle identify the hole- and electron-like states
illustrated in Fig. 4 and 5, respectively. Figure 4: Real parts of the BdG
coefficients for a hole-like solution as a function of $x$ across the
junction, for the momentum indicated by the square in Fig. 3. These are
obtained from the numerical calculations for $\mu=-3.0,\phi=0$. The imaginary
parts vanish identically and are not shown. (a) The $u$ coefficients are
consistent with Eq. 46, as $u^{C}=0$ (blue squares) in the C part extending
from $49\leq x[a]\leq 53$, as expected for a pure hole-like state. Similarly,
$u^{L/R}_{1}$ and $u^{L/R}_{2}$ vanish at the L-C and C-R interfaces,
respectively, in accordance to the analytical calculations. The opposite signs
of these coefficients within the L and R leads and across the junction are
also consistent with the analytical solution with $n-b$ an odd number. (b) The
$v$ coefficients are finite throughout, including the C part, and are
continuous across the junction. These are also consistent with the analytical
solution in Eq. 46. Figure 5: Real part of the BdG coefficients for an
electron-like solution as a function of $x$ across the junction, for the
momentum indicated by the triangle in Fig. 3. As in Fig. 4, these are obtained
from the numerical calculations for $\mu=-3.0,\phi=0$, and the imaginary parts
vanish identically. (a) $u^{C}$ is finite, as expected for a pure electron-
like state, and in agreement with the analytical solution in Eq. 35. The
remaining $u$ coefficients are also consistent with the latter. (b) $v^{C}=0$
again indicates a pure electron-like state. $v^{L/R}_{1}$ and $v^{L/R}_{2}$
vanish at the L-C and C-R interfaces, and their signs are consistent with
$m-a$ and odd, as indicated by Eq. 35.

###  III.2 Non-degenerate orbitals

So far, we have focused on the special case of degenerate orbitals for a
$s\tau_{3}-N-s\tau_{3}$ junction along $x$. We now consider more realistic
cases, where the L and R leads include terms which lift the degeneracy between
the two orbitals. As discussed in Sec. II.1, we consider two types of
symmetry-allowed terms which can lift the orbital degeneracy: (i) intra-
orbital hybridization terms corresponding to $\xi_{3}\tau_{3}$ in the non-
pairing part of the bulk lead Hamiltonians (Eq. 5), and (ii) inter-orbital
hybridization terms corresponding $\xi_{1}\tau_{1}$ in the same model. We
examine the effects of each of these terms on the bound state spectrum
separately.

We first consider intra-orbital hybridization terms only by fixing $t_{1}=1$
and allowing $\delta t=t_{1}-t_{2}$ to vary (Eq. 7). All other parameters are
identical to those in the top panel of Fig. 2. In Fig. 6 (a) we show the bound
state spectrum for the orbital-degenerate case with $\delta t=0$ together with
a case with significant intra-orbital hybridization corresponding to $\delta
t=0.8$. The gapless states for the degenerate orbital case become gapped with
the inclusion of intra-orbital hybridization terms. These break the orbital-
exchange symmetry and thus mix the electron- and hole-like states of the
orbital-degenerate case, in accordance with Sec. III.1. In panel (b), we show
the spectrum for $\delta t=0.8$ as a function of the global phase difference
$\phi$. The eigenvalues change with $\phi$ and we recover gapless states at
$\pi$ phase difference, as is the case for the typical single-channel
junction. The evolution of the bound state spectrum with $\phi$ is also
confirmed via a calculation of the Josephson current in Appendix F.

We now consider the effect of the inter-orbital hybridization terms only,
corresponding to $\xi_{1}$ in Eq. 5. All other parameters are the same in
porevious cases. As shown in Fig. 7 (a) the bound state spectrum remains
gapless under the inclusion of inter-orbital hybridization, in contrast to
intra-orbital hybridization terms. This can also be understood via the
orbital-exchange symmetry, since the inter-orbital hybridization terms
preserve the former. Similarly, the spectrum is invariant under changes in the
global phase difference $\phi$, as indicated in Fig. 7 (b). It is then clear
that, due to intra-orbital hybridization terms, a general
$s\tau_{3}-N-s\tau_{3}$ junction along $x$ exhibits the usual features of a
single-channel junction.

As a final example, we consider the effects of intra- and inter-orbital
hybridization on the junction along $z$. The spectrum in this case can be
labeled by both $k_{x},k_{y}$ conserved momenta. For general values of the
latter, the spectrum is gapped, in analogy to the junction along $x$. However,
along the diagonals of the two-dimensional Brillouin zone $|k_{x}|=|k_{y}|$,
the intra-orbital hybridization terms vanish, and the bound state states again
become gapless. This is illustrated in Fig. 8, and can be understood via the
orbital-exchange symmetry, which is recovered along the diagonals. Thus, in
contrast to the junction along $x$, the junction along $z$ can exhibit gapless
edge states for a rather general model of $s\tau_{3}$ pairing.

Figure 6: The effects of intra-orbital hybridization terms corresponding to
$\xi_{3}$ in the bulk of the L/R leads (Eq. 5) on the bound state spectrum for
a $s\tau_{3}-N-s\tau_{3}$ junction along $x$. We introduce these terms by
fixing $t_{1}=1$ and allowing $\delta t=t_{1}-t_{2}$ to vary (Eq. 7). The
inter-orbital hybridization corresponding to $\xi_{1}$ in Eq. 5 is set to
zero. All other parameters are the same as in Fig. 2 (a) and $\phi=0$. The two
orbitals couple identically to the C part. (a) A gap develops as intra-orbital
hybridization terms are introduced. For clarity, we show the case for
degenerate orbitals $\delta t=0$ (green rhombi), and a case with extreme
intra-orbital hybridization for $\delta t=0.8$ (red squares). Note that a
finite gap opens for any $\delta t\neq 0$. Intra-orbital hybridization terms
break the orbital-exchange symmetry discussed in Sec. II.3, allowing the
electron- and hole-like states in Eqs. 35 and 46 to mix via Andreev
scattering. (b) Effect of a global phase difference $\phi$ across the junction
on the bound state spectrum with finite intra-orbital hybridization $\delta
t=0.8$. In contrast to the degenerate case shown in Fig. 2 (a), the inclusion
of finite intra-orbital hybridization allows the spectrum to change with
$\phi$. Note that the states become gapless at $\phi=\pi$, much like in a
conventional single-channel junction. The sensitivity of the spectrum with
$\phi$ is also confirmed via a calculation of the Josephson current in
Appendix F. Figure 7: The effects of an inter-orbital hybridization term,
corresponding to $\xi_{1}$ (Eq. 8) in the bulk L and R leads, on the bound
state spectrum for a $s\tau_{3}-N-s\tau_{3}$ junction along $x$. The intra-
orbital hybridization terms corresponding to $\xi_{3}$ are set to zero. The
remaining parameters are the same as in Fig. 6. (a) The bound states remain
gapless as $t_{4}$ (Eq. 8) increases, in contrast to the intra-orbital
hybridization case in Fig. 6. The inter-orbital hybridization preserves the
orbital-exchange symmetry defined in Sec. II.3 and thus does not mix the
electron- and hole-like states, which remain adiabatically connected with the
solutions of the orbital-degenerate junction in Eqs. 35 and 46. (b) The
spectrum is insensitive to a global phase difference $\phi$, in contrast to
leads with intra-orbital hybridization. These results are also confirmed via a
calculation of the Josephson current (Appendix F). Figure 8: Bound state
spectrum for a $s\tau_{3}-N-s\tau_{3}$ junction along $z$, calculated from the
model in Sec. II.2, with finite intra- and inter-orbital hybridization. We use
$\delta t=0.4$ and $4t_{4}=0.4$, $\mu=-3.0$ with a NN $t_{z}=0.2$ along $z$.
Note that both $k_{x}$ and $k_{y}$ are conserved in this case. In contrast to
the junction along $x$, gapless states are possible here even when intra-
orbital hybridization is included. The latter vanishes for $|k_{x}|=|k_{y}|$.
For these momenta, we recover the orbital-exchange symmetry defined in Sec.
II.3. (a) bound state eigenvalues as a function of $k_{y}$ for fixed
$k_{x}=0$. The intra-orbital hybridization terms are finite for these momenta,
and the states are gapped, as when the junction is along $x$ (Fig. 6. (b) Same
as (a) for $k_{x}=0.4k_{y}$. The gaps shrink but remain finite. (c) Same as
(a) and (b) along the diagonal of the 2D Brillouin zone for $k_{x}=k_{y}$. The
gap closes as the model recovers the orbital-exchange symmetry, and the
spectrum is similar to that of the junction along $x$ with inter-orbital
hybridization only (Fig. 7).

The bound-state spectra of $s\tau_{3}-N-I$ junctions are very similar to those
of $s\tau_{3}-N-s\tau_{3}$ junctions. In the limit of degenerate orbitals, the
bound states of the former are also electron- or hole like, and they also
become gapless for a set of conserved momenta. The analytical solutions in
this limit are discussed in Appendix C. The evolution with either intra- or
inter-orbital hybridization is also very similar to that in
$s\tau_{3}-N-s\tau_{3}$ junctions and, as such, will not be discussed here in
detail.

##  IV $s\tau_{3}-N-s$ junctions

We consider a $s\tau_{3}-N-s$ junction, where the L lead is in a $s\tau_{3}$
pairing state, while the R lead is in a single-channel $s$-wave state. As in
the case of $s\tau_{3}-N-s\tau_{3}$ junctions, both orbitals on the L side
couple identically to a single orbital in the C metallic part, which in turn
couples to a single orbital in the trivially-paired R lead. The model is
described in Sec. II.1.

We first solve the system in the limit of degenerate orbitals in the L leads,
which couple identically to the single C orbital. The analytical solution is
presented in detail in Appendix E. Here, we summarize some of the most
important results. In contrast to the $s\tau_{3}-N-s\tau_{3}$ case, only the
L-C interface is subject to both open and continuity BC’s. Furthermore, as
explained in Sec. II.3, this model does not preserve the orbital-exchange
symmetry, due to the presence of a single channel in the R lead, even when the
two orbitals are degenerate. Thus, in contrast to the $s\tau_{3}-N-s\tau_{3}$
case, the $s\tau_{3}-N-s$ junction allows for Andreev scattering mixing
electron and hole-like states, and leading to bound states which are generally
gapped. However, in the degenerate orbital limit, the bound states states of
the $s\tau_{3}-N-s$ are invariant under a change in the global relative phase
$\phi$.

The analytical results are confirmed by the numerical solutions. In Fig. 9
(a), we show the bound state spectrum for a $s\tau_{3}-N-s$ junction, whith
parameters similar to those of Fig. 2, for two values of $\mu$. We see the
presence of a gap in both instances, a gap which also occurs for any value of
$\mu$. An analysis of the eigenstates, not shown here for brevity, likewise
indicates that the C spinors involve a mixture of electron- and hole-like
states. Remarkably, the spectrum is invariant under a change in global phase
difference $\phi$ in the degenerate-orbital limit, as shown in Fig. 9 (b). The
spectrum does change with $\phi$ once intra-orbital hybridization terms are
introduced, as in the case of $s\tau_{3}-N-s\tau_{3}$ junctions. The bound-
state variation with $\phi$ is also confirmed via a calculation of the
Josephson current in Appendix F.

Figure 9: Bound state spectrum for a $s\tau_{3}-N-s$ junction along $x$. In
all cases, the two orbitals in the L lead couple identically to the single C
orbital. All other parameters, unless explicitly stated, are the same as in
Fig. 2. (a) Evolution of the spectrum for degenerate orbitals in the L lead
with chemical potential $\mu$. Note that the states are gapped for both values
shown here, as well as for any $\mu$. This stands in contrast to the
$s\tau_{3}-N-s\tau_{3}$ junction, as shown in Fig. 2. The difference can be
attributed to the orbital-exchange symmetry, which is broken for the
$s\tau_{3}-N-s$, thus allowing Andreev scattering which mixes electron- and
hole-like states. (b) The spectrum for the degenerate-orbital case of panel
(a) with $\mu=-3.0$ is invariant under a change in global relative phase
$\phi$. This is due to the absence of intra-orbital hybridization terms, as is
also the case with the $s\tau_{3}-N-s\tau_{3}$ junctions in Fig. 2 (a). (c)
The spectrum does change with $\phi$ once intra-orbital hybridization terms
are introduced for the L lead, in a manner analogous to Fig. 6. This is also
confirmed by the calculation of the Josephson current in Appendix F.

##  V Summary and Discussion

We now summarize our main results and subsequently discuss possible
experimental realizations.

###  V.1 Summary of results

We studied Josephson junctions where the left lead is in an unconventional
$s\tau_{3}$ gapped pairing state which involves two orbitals with opposite-
sign pairing. Similar pairing states were advanced as promising candidates in
alkaline Fe- selenides Nica  _et al._ (2017), where $d_{xz},d_{yz}$ orbitals
provided the two-dimensional manifold for $s\tau_{3}$, as well as in the
heavy-fermion superconductor CeCu2Si2 Nica and Si (2021). We considered
junctions where both orbitals of the left lead couple to a single orbital in
an intermediate metallic part. Such arrangements, although unusual, allow for
a richer phenomenology when compared to the more typical junctions without
cross-coupling. We considered three different types of junctions, referred to
as $s\tau_{3}-N-s\tau_{3}$, $s\tau_{3}-N-I$, and $s\tau_{3}-N-s$, where the
right lead is in a two-orbital $s\tau_{3}$ pairing, insulating, and trivial
single-orbital $s$-wave pairing states, respectively. We discussed both two-
dimensional arrangements with junctions along the $x$-axis, as well as
junctions along $z$.

We studied the three types of junctions in the important limit where the two
orbitals of the $s\tau_{3}$ pairing state, in the left and, when appropriate,
right leads, are degenerate and couple identically to a single orbital of the
central part. Our most striking results can be grouped under two headings. One
is the emergence of purely electron and hole-like bound states, which become
gapless and degenerate for a set of conserved momenta. The other is that the
bound state spectrum is invariant under a change in the global phase
difference between left and right superconducting leads. In both aspects, the
junctions differ sharply from the typical single- or multi-channel junctions
without cross-coupling. The first of these two effects occurs for
$s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junctions along $x$. Here, the
bound states differ from to the typical single-channel junction, where Andreev
scattering mixes purely electron- and hole-like states. The absence of Andreev
scattering in these cases, which leads to gapless electron- and hole-like
states, is due to the combined effects of non-trivial orbital structure of
$s\tau_{3}$ pairing and coupling to a single C orbital. We also find that
these gapless bound states are protected by an orbital-exchange symmetry.
However, we stress that this protection is not due to any topological property
of the junction, unlike the well-known gapless states which occur in a single-
channel $\pi$ junction. The second property, the invariance of the bound state
spectrum with global phase difference, manifests in $s\tau_{3}-N-s\tau_{3}$
and $s\tau_{3}-N-s$ junctions along $x$. The invariance of the bound state
spectrum points toward a vanishing Josephson current Sauls (2018). We
confirmed that this is the case by calculating the Josephson current directly
in the tunneling limit. We stress that $s\tau_{3}-N-s\tau_{3}$ junctions in
the limit of degenerate orbitals differ dramatically from
$s\tau_{0}-N-s\tau_{0}$ junctions in the same limit. The latter involve
orbitally-trivial pairing states which are identical to both orbitals,
including the sign. In clear contrast to the $s\tau_{3}-N-s\tau_{3}$
junctions, the $s\tau_{0}-N-s\tau_{0}$ bound state spectrum is essentially
that of a typical single-channel junction.

We also considered deviations from the degenerate-orbital limit by introducing
intra- and inter-orbital hybridization terms for the normal part of the leads.
These terms are allowed by tetragonal symmetry for $d_{xz},d_{yz}$ orbitals in
the context of alkaline Fe-selenides Raghu  _et al._ (2008). We found that
intra-orbital hybridization terms break the orbital-exchange symmetry and lead
to Andreev scattering which mixes electron and hole-like states, gapping the
bound state spectrum of $s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junctions
along $x$. By contrast, inter-orbital hybridization terms do not break the
orbital-exchange symmetry and preserve the gapless electron- and hole-like
states in the cases mentioned above. The bound states thus behave as in the
intra-orbital hybridization case when both the former and inter-orbital
hybridization are present. Important exceptions can occur as exemplified by
$s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-I$ junctions along $z$, with dominant
NN hopping along $z$. In this case, the bound states become gapless at the
zeroes of the intra-orbital hybridization terms, along the diagonal of the
two-dimensional Brillouin zone in our case, since the junction recovers
orbital-exchange symmetry at these points.

We also considered the effects of lifting the orbital degeneracy in the bulk
of the leads when a finite global phase difference across the junction is
present. Here again we find that the bound state spectrum becomes sensitive to
the phase when intra-orbital hybridization terms are in effect for both
$s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-s$ junctions. This is also confirmed
by calculations of the Josephson current. In contrast, under the addition of
inter-orbital hybridization terms only, the bound state spectrum remains
invariant, also in accordance with a vanishing Josephson current.

The two crucial aspects behind the striking behavior of all of the junctions
considered here involve two-orbital matrix-pairing $s\tau_{3}$ in the leads,
together with identical coupling of both of these orbitals to a single orbital
in the central metallic parts. These setups ensure that the two orbital
sectors of $s\tau_{3}$ pairing are entangled in a non-trivial way near the
interface with the metallic part, even when the two sectors are decoupled far
into the bulk of the leads. This is to be contrasted with cases where the
pairing has trivial matrix structure, as exemplified by $s\tau_{0}$ pairing.
In those cases, the trivial pairing matrix structure means that two orbitals
which are decoupled in the bulk of the leads remain so near the interfaces
with the metallic part. Likewise, $s\tau_{3}$ junctions with central metallic
regions which likewise involve two orbitals, which couple in a ”one-to-one”
manner with those of $s\tau_{3}$ would fail to fully capture the effects of
the non-trivial matrix structure of $s\tau_{3}$. As exemplified by
$s\tau_{3}-N-I$ junctions, the interplay between non-trivial interfaces and
non-trivial matrix-pairing suggests that exotic edge states could also be
engineered for materials which likely involve $s\tau_{3}$ pairing, such as the
alkaline Fe-selenides.

###  V.2 Experimental signatures

Having given a summary of our main results, we now discuss their potential
experimental signatures. We stress that, in the limit of degenerate orbitals,
junctions involving $s\tau_{3}$ pairing are dramatically different from the
typical single-channel junction as well as from junctions involving pairing
states with trivial orbital structure, such as the $s\tau_{0}$ states.
However, any realizations of the $s\tau_{3}$ junctions proposed here will
involve deviations from the ideal case. We argue that such deviations can be
made small, in the sense discussed below, allowing a partial observation of
the striking properties of the ideal case.

We consider three likely deviations from the ideal case in the form of: (i)
unequal coupling of the two orbitals to the single-orbital of the metallic
part, (ii) lifting of the degeneracy of the two orbitals in the leads via
symmetry-allowed intra- and inter-orbital hybridization, and (iii) weak
disorder in the junction. For (i), unequal coupling to the C part gaps the
electron- and hole-like states and induces a nonzero Josephson effect, even
when the orbitals of the leads are degenerate. A single-orbital intermediate
region amounts to an effectively single-band system which remains non-
superconducting for temperatures above the critical temperature of
$s\tau_{3}$, which can be estimated from the transition temperature in
alkaline Fe-selenides Lee (2017). The condition of equal coupling to the
center part is unlikely to occur for junctions along the in-plane axes of the
tetragonal Fe-selenides i.e. along $x$ or $y$, since $d_{xz}$ and $d_{yz}$
orbitals typically cannot couple identically to any other in-plane orbital.
However, junctions along $z$ provide better candidates in this context, as the
the lobes of $d_{xz}$ and $d_{yz}$ are likely to have comparable overlap with
such orbitals as $p_{z}$, $d_{z^{2}}$, and $d_{xy}$ along the $z$ axis.

A second deviation from the ideal case occurs due to the presence of both
intra- and inter-orbital hybridization terms in the leads for junctions
involving alkaline Fe-selenide leads. As we have shown, such terms will lead
to gapping of the electron- and hole-like states and will furthermore ensure
that a static Josephson effect is present. Note however, that, as proposed in
Refs. Nica  _et al._ (2017); Nica and Si (2021), and as reiterated in our
text, $s\tau_{3}$ pairing induces a full gap in the bulk of the leads when the
pairing amplitude exceeds the band splitting near the FS, which is governed by
the intra- and inter-orbital hybridization terms. We expect that this limits
the effects of band splitting in the case of the junctions, such that the
induced gap for the bound states will be small compared to the bulk gap in the
leads. Similarly, the Josephson current, which depends on the strength of the
intra-orbital hybridization terms, will likely be finite but strongly
suppressed, as compared to a similar setup for pairing states with trivial
orbital structure. However, it should be borne in mind that exceptions to this
behavior can occur for junctions along $z$. The intra-orbital hybridization
terms can vanish by symmetry along certain directions in the two-dimensional
Brillouin zone. Here, the junction recovers the orbital-exchange symmetry and
the gapless bound states. We do expect a finite Josephson current in these
cases, due to the contribution of bulk states away from these points.

The third source of deviation from the ideal case is the presence of disorder
in the junction. The gapless bound states in the ideal degenerate cases are
not due to any non-trivial topology, although they are protected by the
orbital-exchange symmetry. Therefore, these states are not robust against
disorder.

Due to these inherent limitations, the junctions discussed here cannot
singlehandedly indicate $s\tau_{3}$ pairing, but they can nonetheless provide
strong supporting evidence.

##  VI Acknowledgements

OE acknowledge support from National Science Foundation Awards No. DMR
1904716. EMN is supported by ASU startup grant. Work at Rice has been
supported by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences, under Award No. DE-SC0018197 and the Robert A. Welch Foundation
under Grant No. C-1411 (Q.S.). One of us (Q.S.) acknowledges the hospitality
of the Aspen Center for Physics, which is supported by NSF grant No.
PHY-1607611.

##  Appendix A Analytical solutions of the $s\tau_{3}-N-s\tau_{3}$ junction
with degenerate orbitals

We consider a $s\tau_{3}-N-s\tau_{3}$ junction with degenerate orbitals/bands
in the normal states corresponding to $t_{1}=t_{2},t_{4}=0$ in Eqs. 9, 8
together with identical coupling to the single channel of the C part
corresponding to $V_{1}=V_{2}$ in Eq. 13. In this limit, only the symmetric
linear combination of the two orbitals at either L-C and C-R interfaces
couples to the single orbital of the C part.

We apply a local unitary transformation to the L lead

$\displaystyle\begin{pmatrix}c_{\mathbf{r},I,\uparrow}\\\
c_{\mathbf{r},II,\uparrow}\\\ c^{{\dagger}}_{\mathbf{r},I,\downarrow}\\\
c^{{\dagger}}_{\mathbf{r},II,\downarrow}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1&0&0\\\
1&-1&0&0\\\ 0&0&1&1\\\
0&0&1&-1\end{pmatrix}\begin{pmatrix}c_{\mathbf{r},1,\uparrow}\\\
c_{\mathbf{r},2,\uparrow}\\\ c^{{\dagger}}_{\mathbf{r},1,\downarrow}\\\
c^{{\dagger}}_{\mathbf{r},2,\downarrow}\end{pmatrix}.$ (64)

Under this transformation, only channel $I$, corresponding to the symmetric
linear combination of the operators for the two orbitals, couples to the
single orbital of the C part. Furthermore, the transformation leaves the
orbital-diagonal tight-binding part of the Hamiltonian invariant
(corresponding to $\xi_{0}$ in Eq. 5) while it transforms the $s\tau_{3}$
pairing as

$\displaystyle H_{\text{L, Pair}}=$
$\displaystyle\sum_{\mathbf{r}}\Delta\left(c^{{\dagger}}_{\mathbf{r},1\uparrow}c^{{\dagger}}_{\mathbf{r},1\downarrow}-c^{{\dagger}}_{\mathbf{r},2\uparrow}c^{{\dagger}}_{\mathbf{r},2\downarrow}\right)+\text{H.c.}$
$\displaystyle\rightarrow\sum_{\mathbf{r}}\Delta\left(c^{{\dagger}}_{\mathbf{r},I\uparrow}c^{{\dagger}}_{\mathbf{r},II\downarrow}-c^{{\dagger}}_{\mathbf{r},II\uparrow}c^{{\dagger}}_{\mathbf{r},I\downarrow}\right)\text{+H.c.},$
(65)

where for simplicity we only consider one spin sector. The model thus reduces
to solving for channel $I$ as in a typical, single-channel Josephson junction,
while channel $II$ must obey open BC’s at the L-C interface.

We proceed to find solutions of the model in the transformed basis which lie
below the bulk gap, corresponding to solutions which decay into the L lead. We
first apply a Fourier transform along $y$ and linearize the BdG equations in
the vicinity of two points $(\alpha K_{Fx}+q_{x},K_{Fy}+q_{y})$ on the FS,
with $q_{x},q_{y}$ small, and $\alpha=\pm 1$. This is in analogy to the
typical single-channel solution Sauls (2018). After identifying
$q_{x}\rightarrow-i\partial_{x}$ in the continuum limit, we separate the BdG
ansatz of the L lead into fast- and slowly-varying parts as

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};L}\\\
\tilde{u}_{II,k_{y};L}\\\ \tilde{v}_{I,k_{y};L}\\\
\tilde{v}_{II,k_{y};L}\end{pmatrix}=$ $\displaystyle
e^{iK_{Fy}y}e^{iK_{Fx}x}e^{\kappa\left(x+\frac{l}{2}\right)}\begin{pmatrix}u_{\alpha=1,I;L}\\\
u_{\alpha=1,II;L}\\\ v_{\alpha=1,I;L}\\\
v_{\alpha=1,II;L}\end{pmatrix}+e^{iK_{Fy}y}e^{-iK_{F}x}e^{\kappa\left(x+\frac{l}{2}\right)}\begin{pmatrix}u_{\alpha=\bar{1},I;L}\\\
u_{\alpha=\bar{1},II;L}\\\ v_{\alpha=\bar{1},I;L}\\\
v_{\alpha=\bar{1},II;L}\end{pmatrix}.$ (66)

Note that the C part is between $-l/2\leq x\leq l/2$. Since we are looking for
solutions which decay into the L lead, we take $\kappa$ to be real and
positive. Also note that we consider an ansatz which is a linear superposition
of solutions with opposite momenta along $x$. This is due to the BdG
coefficients of channel $II$ which must obey open BC at the L-C junction,
while the momentum along $y$ is conserved. We ignore corrections due to
$q_{y}a$.

The slowly-varying parts obey the BdG equation

$\displaystyle\begin{pmatrix}\alpha\left[-iv_{Fx}\kappa\right]-\epsilon&0&0&\Delta\\\
0&\alpha\left[-iv_{Fx}\kappa\right]-\epsilon&\Delta&0\\\
0&\Delta^{*}&-\alpha\left[-iv_{Fx}\kappa\right]-\epsilon&0\\\
\Delta^{*}&0&0&-\alpha\left[-iv_{Fx}\kappa\right]-\epsilon\end{pmatrix}\begin{pmatrix}u_{\alpha,I}\\\
u_{\alpha,II}\\\ v_{\alpha,I}\\\ v_{\alpha,II}\end{pmatrix}=$ $\displaystyle
0,$ (67)

where we introduced the Fermi velocities along $x$ via

$\displaystyle v_{Fx}=$ $\displaystyle 2t_{1}\sin(K_{Fx}a).$ (68)

Note that the first and fourth, and second and third rows respectively
decouple in the bulk Hamiltonian and can be solved independently. We find that
the solutions in the bulk of the L lead are

$\begin{pmatrix}u_{\alpha,I}\\\ u_{\alpha,II}\\\ v_{\alpha,I}\\\
v_{\alpha,II}\end{pmatrix}=\begin{pmatrix}A_{\alpha}\left(\epsilon-i\alpha\Lambda\right)\\\
B_{\alpha}\left(\epsilon-i\alpha\Lambda\right)\\\ B_{\alpha}\Delta^{*}\\\
A_{\alpha}\Delta^{*}\end{pmatrix}$ (69)

where $A_{\alpha},B_{\alpha}$ are coefficients to be determined from the BC’s
and

$\displaystyle\kappa=$ $\displaystyle\frac{\Lambda}{v_{Fx}}$ (70)
$\displaystyle\Lambda=$
$\displaystyle\sqrt{\Delta^{2}-\epsilon^{2}},~{}\text{for}~{}\epsilon<|\Delta|.$
(71)

The general solution in the C part, extending from $-l/2\leq x\leq l/2$ can be
similarly determined to be of the form

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y};C}\\\ \tilde{v}_{k_{y};C}\\\
\end{pmatrix}=$ $\displaystyle
e^{iK_{F}x}\begin{pmatrix}E_{1}e^{i\frac{\epsilon x}{v_{Fx}}}\\\
G_{1}e^{-i\frac{\epsilon
x}{v_{Fx}}}\end{pmatrix}+e^{-iK_{F}x}\begin{pmatrix}E_{\bar{1}}e^{-i\frac{\epsilon
x}{v_{Fx}}}\\\ G_{\bar{1}}e^{i\frac{\epsilon x}{v_{Fx}}}\end{pmatrix}.$ (72)

The solutions in the R lead, which decay away from the C-R interface for
$x>l/2$ are

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};R}\\\
\tilde{u}_{II,k_{y};R}\\\ \tilde{v}_{I,k_{y};R}\\\
\tilde{v}_{II,k_{y};R}\end{pmatrix}=$ $\displaystyle
e^{iK_{F}x}e^{-\kappa\left(x-\frac{l}{2}\right)}\begin{pmatrix}M_{1}(\epsilon+i\Lambda)\\\
N_{1}(\epsilon+i\Lambda)\\\ N_{1}|\Delta|e^{-i\phi}\\\
M_{1}|\Delta|e^{-i\phi}\end{pmatrix}+e^{-iK_{F}x}e^{-\kappa\left(x-\frac{l}{2}\right)}\begin{pmatrix}M_{\bar{1}}(\epsilon-i\Lambda)\\\
N_{\bar{1}}(\epsilon-i\Lambda)\\\ N_{\bar{1}}|\Delta|e^{-i\phi}\\\
M_{\bar{1}}|\Delta|e^{-i\phi}\end{pmatrix}.$ (73)

Note that we have introduced a global phase $\phi$ in the pairing of the R
lead.

For our purposes, it proves convenient to parameterize all of the coefficients
in terms of an overall phase and a relative phase as in

$\displaystyle A_{\alpha}=$
$\displaystyle|A|e^{i\theta^{0}_{A}}e^{i\alpha\theta_{A}},$ (74)

and similarly for all $\alpha$-dependent quantities. We also parameterize the
factors which enter in the general solutions for the leads as

$\displaystyle\epsilon\pm i\alpha\Lambda=$ $\displaystyle|\Delta|e^{\pm
i\alpha\theta}$ (75) $\displaystyle\theta=$
$\displaystyle\arg\left(\epsilon+i\Lambda\right).$ (76)

We now consider the BC’s. As discussed previously, channel $II$ in the L and R
leads does not couple to the C part, and thus must obey open BC’s at the L-C
and C-R interfaces, respectively:

$\displaystyle\tilde{u}_{II,k_{y};L}\left(x=\frac{-l}{2}\right)=$
$\displaystyle 0$ (77)
$\displaystyle\tilde{v}_{II,k_{y};L}\left(x=\frac{-l}{2}\right)=$
$\displaystyle 0$ (78)
$\displaystyle\tilde{u}_{II,k_{y};R}\left(x=\frac{l}{2}\right)=$
$\displaystyle 0$ (79)
$\displaystyle\tilde{v}_{II,k_{y};R}\left(x=\frac{l}{2}\right)=$
$\displaystyle 0.$ (80)

By contrast, channel $I$ couples across the junction and therefore satisfies
the continuity conditions

$\displaystyle\tilde{u}_{I,k_{y};L}\left(x=\frac{-l}{2}\right)=$
$\displaystyle\tilde{u}_{k_{y};C}\left(x=\frac{-l}{2}\right)$ (81)
$\displaystyle\tilde{v}_{I,k_{y};L}\left(x=\frac{-l}{2}\right)=$
$\displaystyle\tilde{v}_{k_{y};C}\left(x=\frac{-l}{2}\right)$ (82)
$\displaystyle\tilde{u}_{I,k_{y};R}\left(x=\frac{l}{2}\right)=$
$\displaystyle\tilde{u}_{k_{y};C}\left(x=\frac{l}{2}\right)$ (83)
$\displaystyle\tilde{v}_{I,k_{y};R}\left(x=\frac{l}{2}\right)=$
$\displaystyle\tilde{v}_{k_{y};C}\left(x=\frac{l}{2}\right)$ (84)

Expressed in term of the $\alpha$\- dependent quantities these amount to

$\displaystyle-\frac{K_{Fx}l}{2}+\theta_{B}-\theta=$
$\displaystyle\frac{\pi}{2}+b\pi$ (85)
$\displaystyle-\frac{K_{Fx}l}{2}+\theta_{A}=$
$\displaystyle\frac{\pi}{2}+a\pi$ (86)
$\displaystyle\frac{K_{Fx}l}{2}+\theta_{N}+\theta=$
$\displaystyle\frac{\pi}{2}+n\pi$ (87)
$\displaystyle\frac{K_{Fx}l}{2}+\theta_{M}=$
$\displaystyle\frac{\pi}{2}+m\pi,$ (88)

for the open BC’s. $a,b,m$, and $n$ are arbitrary integers. The continuity
conditions imply that

$\displaystyle\theta_{A}-\theta+\frac{\epsilon l}{v_{Fx}}=$
$\displaystyle\theta_{M}+\theta$ (89) $\displaystyle\theta_{B}-\frac{\epsilon
l}{v_{Fx}}=$ $\displaystyle\theta_{N}$ (90)

$\displaystyle|A|=$ $\displaystyle|M|$ (91) $\displaystyle|B|=$
$\displaystyle|N|$ (92)

$\displaystyle\theta^{0}_{A}=$ $\displaystyle\theta^{0}_{M}$ (93)
$\displaystyle\theta^{0}_{B}=$ $\displaystyle\theta^{0}_{M}-\phi.$ (94)

Note that the _overall phase_ $\theta^{0}_{B}$ can incorporate the global
phase difference of the pairing $\phi$. This is in contrast to the typical
single-channel case and it leads to an insensitivity of the bound state
spectrum w.r.t. to $\phi$.

Note that there are five unknowns consisting of the five relative phase
$\theta_{A},\theta_{B},\theta_{M},\theta_{N}$ together with the eigenenergies
$\epsilon$. However, due to the open BC’s, there are six equations. Therefore,
non-trivial solutions cannot be found for this system of equations. We
consider instead solutions where either $A_{\alpha},M_{\alpha}$ or
$B_{\alpha},N_{\alpha}$ are trivially zero. In either of these cases, the
system involving the relative phase and $\epsilon$ reduces to three equations
with three unknowns. Importantly, these solutions are found to be either hole-
or electron-like, as either the $u$ or $v$ BdG coefficients vanish in the C
part.

The electron-like solutions obtained in this manner are

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};L-e}\\\
\tilde{u}_{II,k_{y};L-e}\\\ \tilde{v}_{I,k_{y};L-e}\\\
\tilde{v}_{II,k_{y};L-e}\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta^{0}_{A}}\begin{pmatrix}\cos\left[K_{Fx}x-\frac{\epsilon
l}{2v_{Fx}}+\frac{(a+m+1)\pi}{2}\right]\\\ 0\\\ 0\\\
\cos\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{\pi}{2}+a\pi\right]\end{pmatrix}.$
(95)

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y};C-e}\\\ \tilde{v}_{k_{y};C-e}\\\
\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta^{0}_{A}}\begin{pmatrix}\cos\left[K_{Fx}x+\frac{\epsilon
x}{v_{Fx}}+\frac{(a+m+1)\pi}{2}\right]\\\ 0\end{pmatrix}.$ (96)

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};R-e}\\\
\tilde{u}_{II,k_{y};R-e}\\\ \tilde{v}_{I,k_{y};R-e}\\\
\tilde{v}_{II,k_{y};R-e}\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta^{0}_{A}}\begin{pmatrix}\cos\left[K_{Fx}x+\frac{\epsilon
l}{2v_{Fx}}+\frac{(a+m+1)\pi}{2}\right]\\\ 0\\\ 0\\\
\cos\left[K_{Fx}\left(x-\frac{l}{2}\right)+\frac{\pi}{2}+m\pi\right]\end{pmatrix}.$
(97)

Using

$\displaystyle\arctan(\theta)=\frac{\Lambda}{\epsilon}$ (98)

we determine the eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon
l}{2v_{Fx}}+\frac{K_{Fx}l}{2}+\frac{(a-m)\pi}{2}\right).$ (99)

Similarly, the hole-like solutions are

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};L-h}\\\
\tilde{u}_{II,k_{y};L-h}\\\ \tilde{v}_{I,k_{y};L-h}\\\
\tilde{v}_{II,k_{y};L-h}\end{pmatrix}=$ $\displaystyle
2|B|e^{i\theta^{B}_{0}}\begin{pmatrix}0\\\
\cos\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{\pi}{2}+b\pi\right]\\\
\cos\left[K_{Fx}x+\frac{\epsilon
l}{2v_{Fx}}+\frac{\pi}{2}+\frac{(n+b+1)\pi}{2}\right]\\\ 0\end{pmatrix}.$
(100)

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y};C-h}\\\ \tilde{v}_{k_{y};C-h}\\\
\end{pmatrix}=$ $\displaystyle 2|B|e^{i\theta^{B}_{0}}\begin{pmatrix}0\\\
\cos\left[K_{Fx}x-\frac{\epsilon
x}{v_{Fx}}+\frac{(n+b+1)\pi}{2}\right]\end{pmatrix}.$ (101)

$\displaystyle\begin{pmatrix}\tilde{u}_{I,k_{y};R-h}\\\
\tilde{u}_{II,k_{y};R-h}\\\ \tilde{v}_{I,k_{y};R-h}\\\
\tilde{v}_{II,k_{y};R-h}\end{pmatrix}=$ $\displaystyle
2|B|e^{i(\theta^{B}_{0}+\phi)}\begin{pmatrix}0\\\
\cos\left[K_{Fx}\left(x-\frac{l}{2}\right)+\frac{\pi}{2}+n\pi\right]\\\
\cos\left[K_{Fx}x-\frac{\epsilon l}{2v_{Fx}}+\frac{(n+b+1)\pi}{2}\right]\\\
0\end{pmatrix}.$ (102)

with eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon
l}{2v_{Fx}}-\frac{K_{Fx}l}{2}+\frac{(n-b)\pi}{2}\right).$ (103)

The solutions in the original basis, obtained via the inverse transformation
of Eq. 64 can be obtained from the solutions shown above.

##  Appendix B $s\tau_{3}-N-s\tau_{3}$ junctions with unequal coupling to the
C part

In Fig. 10 (a), we show the bound state spectrum for a $s\tau_{3}-N-s\tau_{3}$
junction along $x$ as a function of $\delta V=V_{1}-V_{2}$, the difference
between the coupling constants of the two degenerate orbitals in either lead
to the single orbital of the C part. The results show that, in clear contrast
to the case of equal-coupling $\delta V=0$, a finite $\delta V$ breaks the
orbital-exchange symmetry and gaps the purely electron- and hole-like states.
In panel (b) we show the spectrum as a function of $\phi$, for fixed $\delta
V=0.7$. Again in contrast to the $\delta V=0$ cases, the spectrum becomes
dependent on $\phi$, a feature which is also confirmed via a calculation of
the Josephson current in Appendix F.

Figure 10: Bound state spectrum for a $s\tau_{3}-N-s\tau_{3}$ junction along
$x$, where the couplings to the single orbital of the C part, $V_{1}$ and
$V_{2}$ (Eq. 13), are unequal. We consider degenerate orbitals and compare to
the results for $V_{1}=V_{2}$ in Fig. 2 (a). In contrast to the case with
equal $V_{1}=V_{2}$, the bound states acquire a gap with increasing $\delta
V=V_{1}-V_{2}$. This can be understood via a broken orbital-exchange symmetry
which mixes the purely electron- and hole-like states via Andreev scattering.
(b) The spectrum changes as a function of global phase difference across the
junction $\phi$. The results are for $\delta V=0.3$. This is in contrast to
the case with $\delta V=0$, which is shown in Fig. 2.

##  Appendix C $s\tau_{0}-N-s\tau_{0}$ junctions

In this section, we consider $s\tau_{0}-N-s\tau_{0}$ junctions along $x$,
where the leads are in $s\tau_{0}$ pairing states. In these cases, the pairing
functions for both orbitals are identical. We consider degenerate orbitals
which couple to a single orbital in the C part identically. We show that, in
clear contrast to $s\tau_{3}-N-s\tau_{3}$ junctions, $s\tau_{0}-N-s\tau_{0}$
junctions behave essentially like a typical single-channel junction. In Fig.
11 we illustrate the evolution of the bound state spectrum as a function of
global phase difference $\phi$.

Figure 11: Bound state spectrum for a $s\tau_{0}-N-s\tau_{0}$ junction along
$x$, where the two orbitals have identical pairing functions, as a function of
global phase $\phi$. In contrast to $s\tau_{3}$ pairing, $s\tau_{0}$ is an
orbitally-trivial pairing state. The overall setup of the junction is the same
as in the $s\tau_{3}$ cases, with each of the two degenerate orbitals in
either L and R leads coupling identically to a single orbital in the C part.
All of the parameters of the model are the same as in Fig. 2 (a). As indicated
by the results, the bound state spectrum for a $s\tau_{0}-N-s\tau_{0}$
junction differs dramatically from a $s\tau_{3}-N-s\tau_{3}$ junction in two
main aspects: (i) The spectrum for $\phi=0$ is gapped, and (ii) for varying
$\phi$ the spectrum evolves much like a typical single-channel junction,
becoming gapless at $\phi=\pi$. In $s\tau_{0}$ junctions, the symmetric linear
combination of the two orbitals couples across the junction as in the single-
orbital case, while the remaining antisymmetric linear combination decouples
entirely.

##  Appendix D Analytical solutions of the $s\tau_{3}-N-I$ junction with
degenerate orbitals

The solution in the continuum limit is very similar to that of the
$s\tau_{3}-N-s\tau_{3}$ case. Instead of continuity BC’s at the C-R interface,
the C spinors obey open BC’s. Using the ansatze and the conventions of
Appendix A, these conditions amount to

$\displaystyle\frac{K_{Fx}l}{2}+\frac{\epsilon l}{2v_{Fx}}+\theta_{E}=$
$\displaystyle\frac{\pi}{2}+p\pi$ (104)
$\displaystyle\frac{K_{Fx}l}{2}-\frac{\epsilon l}{2v_{Fx}}+\theta_{G}=$
$\displaystyle\frac{\pi}{2}+r\pi,$ (105)

where $p,r$ are integers. As for the $s\tau_{3}-N-s\tau_{3}$ junctions, we
obtain electron- and hole-like solutions of the form

$\displaystyle\begin{pmatrix}u^{L}_{1e}&u^{L}_{2e}&v^{L}_{1e}&v^{L}_{2e}\end{pmatrix}^{T}$
$\displaystyle=$ $\displaystyle
2|A|e^{i\theta^{A}_{0}}|\Delta|\begin{pmatrix}\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)-\frac{\epsilon
l}{v_{Fx}}+p\pi\right]\\\
\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)-\frac{\epsilon
l}{v_{Fx}}+p\pi\right]\\\
\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+a\pi\right]\\\
-\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+a\pi\right]\end{pmatrix}$ (106)

$\displaystyle\begin{pmatrix}u^{C}_{e}v^{C}_{e}\end{pmatrix}^{T}$
$\displaystyle=$ $\displaystyle
2|A|e^{i\theta_{0}}|\Delta|\begin{pmatrix}\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)+\frac{\epsilon}{v_{Fx}}\left(x-\frac{l}{2}\right)+p\pi\right]\\\
0\end{pmatrix}$ (107)

with eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon l}{v_{Fx}}+K_{Fx}l-(p-a)\pi\right)$
(108)

for the electron-like states. The hole-like solutions are

$\displaystyle\begin{pmatrix}u^{L}_{1h}&u^{L}_{2h}&v^{L}_{1h}&v^{L}_{2h}\end{pmatrix}^{T}$
$\displaystyle=$ $\displaystyle
2|B|e^{i\theta^{B}_{0}}|\Delta|\begin{pmatrix}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+b\pi\right]\\\
-\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+b\pi\right]\\\
\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)+\frac{\epsilon
l}{v_{Fx}}+r\pi\right]\\\
\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)+\frac{\epsilon
l}{v_{Fx}}+r\pi\right]\end{pmatrix}$ (109)

$\displaystyle\begin{pmatrix}u^{C}_{h}v^{C}_{h}\end{pmatrix}^{T}$
$\displaystyle=$ $\displaystyle
2|B|e^{i\theta^{0}_{B}}|\Delta|\begin{pmatrix}0\\\
\sin\left[K_{Fx}\left(x-\frac{l}{2}\right)-\frac{\epsilon}{v_{Fx}}\left(x-\frac{l}{2}\right)+r\pi\right]\end{pmatrix}$
(110)

with eigenvalues

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{\epsilon l}{v_{Fx}}-K_{Fx}l+(r-b)\pi\right).$
(111)

As before, $|A|,|B|$ are normalization constants, $\theta^{0}_{A/B}$ are
arbitrary phases, while $a,b,p,r$ are arbitrary integers.

##  Appendix E Analytical solutions of the $s\tau_{3}-N-s$ junction with
degenerate orbitals

Consider a junction of the $s\tau_{3}-N-s$ along $x$, where the R lead is in a
single-channel, $s$-wave pairing state. The model for this junction was
introduced in Sec. II.1. Here, we tackle this analytically for the case where
the two channels of the L leads, which are in a $s\tau_{3}$ pairing state, are
degenerate, and couple identically to the orbital in the C part.

We proceed along the same lines as the $s\tau_{3}-N-s\tau_{3}$ case in the
continuum limit. As in that case, only the symmetric linear combination of the
two channels of the L lead couples to the C part. We therefore consider the
same ansatz for the BdG coefficients in the L lead which also obey the same
open and continuity BC at the L-C interface as before (Eqs. 81-82). The main
distinction is due to the presence of a single channel in the R lead with
general solutions in the bulk given by

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y};R}\\\
\tilde{v}_{k_{y};R}\end{pmatrix}=$ $\displaystyle
e^{iK_{F}x}e^{-\kappa\left(x-\frac{l}{2}\right)}\begin{pmatrix}M_{1}(\epsilon+i\Lambda)\\\
M_{1}|\Delta|e^{-i\phi}\end{pmatrix}+e^{-iK_{F}x}e^{-\kappa\left(x-\frac{l}{2}\right)}\begin{pmatrix}M_{\bar{1}}(\epsilon-i\Lambda)\\\
M_{\bar{1}}|\Delta|e^{-i\phi}\\\ \end{pmatrix}.$ (112)

Furthermore, these BdG coefficients obey a single continuity BC at the C-R
interface. Adopting the same conventions as in the $s\tau_{3}-N-s\tau_{3}$
case, we summarize the boundary conditions as

$\displaystyle-\frac{K_{Fx}l}{2}+\theta_{B}-\theta=$
$\displaystyle\frac{\pi}{2}+b\pi$ (113)
$\displaystyle-\frac{K_{Fx}l}{2}+\theta_{A}=$
$\displaystyle\frac{\pi}{2}+a\pi$ (114)

$\displaystyle\theta_{A}-\theta+\frac{\epsilon l}{v_{Fx}}=$
$\displaystyle\theta_{M}+\theta$ (115) $\displaystyle\theta_{B}-\frac{\epsilon
l}{v_{Fx}}=$ $\displaystyle\theta_{M}$ (116)

for the relative phases and

$\displaystyle|A|=$ $\displaystyle|M|$ (117) $\displaystyle|B|=$
$\displaystyle|N|$ (118)

$\displaystyle\theta^{0}_{A}=$ $\displaystyle\theta^{M}_{0}$ (119)
$\displaystyle\theta^{0}_{B}=$ $\displaystyle\theta^{0}_{M}-\phi.$ (120)

for the amplitudes and global phases of the BdG coefficients. $a,b$ are
arbitrary integers.

In contrast to the case of the $s\tau_{3}-N-s\tau_{3}$ junction, the absence
of additional open BC for the R lead ensures that the equations for the
relative phases and $\epsilon$ have a unique solution. After straightforward
algebra, we determine the eigenvalues from

$\displaystyle\frac{\epsilon}{|\Delta|}=$
$\displaystyle\pm\cos\left(\frac{2\epsilon
l}{3v_{Fx}}+\frac{(a-b)\pi}{3}\right).$ (121)

The corresponding states are

$\displaystyle\begin{pmatrix}\tilde{u}_{1,k_{y},L}\\\
\tilde{u}_{2,k_{y},L,e}\\\ \tilde{v}_{1,k_{y},L,e}\\\
\tilde{v}_{2,k_{y},L,e}\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta^{A}_{0}}\begin{pmatrix}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)-\frac{2\epsilon
l}{3v_{Fx}}+\frac{(2a+b)\pi}{3}\right]+e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+b\pi\right]\\\
\\\ \sin\left[K_{Fx}\left(x+\frac{l}{2}\right)-\frac{2\epsilon
l}{3v_{Fx}}+\frac{(2a+b)\pi}{3}\right]-e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+b\pi\right]\\\
\\\ e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{2\epsilon
l}{3v_{Fx}}+\frac{(a+2b)\pi}{3}\right]+\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+a\pi\right]\\\
\\\ e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{2\epsilon
l}{3v_{Fx}}+\frac{(a+2b)\pi}{3}\right]-\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+a\pi\right]\end{pmatrix}.$
(122)

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y},C}\\\ \tilde{v}_{k_{y},C}\\\
\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta_{0}}\begin{pmatrix}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{\epsilon
x}{v_{Fx}}-\frac{\epsilon l}{6v_{Fx}}+\frac{(2a+b)\pi}{3}\right]\\\ \\\
e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)-\frac{\epsilon
x}{v_{Fx}}+\frac{\epsilon
l}{6v_{Fx}}+\frac{(a+2b)\pi}{3}\right]\end{pmatrix}.$ (123)

$\displaystyle\begin{pmatrix}\tilde{u}_{k_{y},R}\\\ \tilde{v}_{k_{y},R}\\\
\end{pmatrix}=$ $\displaystyle
2|A||\Delta|e^{i\theta_{0}}\begin{pmatrix}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)+\frac{\epsilon
l}{3v_{Fx}}+\frac{(2a+b)\pi}{3}\right]\\\ \\\
e^{-i\phi}\sin\left[K_{Fx}\left(x+\frac{l}{2}\right)-\frac{\epsilon
l}{3v_{Fx}}+\frac{(a+2b)\pi}{3}\right]\end{pmatrix}.$ (124)

Note that the C spinor is a mixture of electron- and hole-like solutions, as
in the case with typical Andreev bound states.

##  Appendix F Josephson current in the tunneling limit

In Sec. III.1 and IV, we showed that in the limit of degenerate orbitals which
couple identically to a single orbital in the C part, $s\tau_{3}-N-s\tau_{3}$
and $s\tau_{3}-N-s$ junctions exhibit bound state spectra which are
insensitive to changes in the relative global phase $\phi$. In view of the
relation between the Josephson current and the derivative of the ground-state
energy with $\phi$, which is typically determined by the Andreev bound state
spectrum, our results suggest the absence of a static Josephson effect.

In order to confirm these findings, we consider simplified models of the
$s\tau_{3}-N-s\tau_{3}$ and $s\tau_{3}-N-s$ junctions and determine the
Josephson current in the tunneling limit. We follow the standard approach in
Ref. Sauls (2018). More precisely, we introduce a Hamiltonian which includes
the leads but which also involves direct tunneling between the latter:

$\displaystyle H=H_{\text{R}}+H_{\text{L}}+H_{\text{T}}.$ (125)

For $s\tau_{3}-N-s\tau_{3}$ junctions, $H_{\text{L/R}}$ are those of Eq. 5,
while for $s\tau_{3}-N-s$ junctions, $H_{\text{R}}$ is a single-channel
$s$-wave bulk Hamiltonian, determined by the dispersion $\xi_{0}$ and $\Delta$
which also enter the expression for the two-orbital $H_{\text{L}}$. The two
leads are connected via a tunneling Hamiltonian $H_{\text{T}}$. For the
$s\tau_{3}-N-s\tau_{3}$ junctions, this takes the form

$\displaystyle H_{\text{T},s\tau_{3}-N-s\tau_{3}}=$
$\displaystyle\sum_{\mathbf{k}\mathbf{p}\sigma}\bigg{(}T^{11}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{R},1\mathbf{k}\sigma}c_{\text{L},1\mathbf{p}\sigma}+T^{21}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{R},2\mathbf{k}\sigma}c_{\text{L},1\mathbf{p}\sigma}$
$\displaystyle+$ $\displaystyle
T^{12}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{R},1\mathbf{k}\sigma}c_{\text{L},2\mathbf{p}\sigma}+T^{22}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{L},2\mathbf{k}\sigma}c_{\text{L},2\mathbf{p}\sigma}+\text{H.c.}\bigg{)},$
(126)

where $T^{ij}_{\mathbf{k}\mathbf{p}}$ are tunneling matrix elements, and where
we introduced R and L indices for the two leads in addition to the $1,2$
indices for the two orbitals. For the $s\tau_{3}-N-s$ junction $H_{\text{T}}$
has the simpler form

$\displaystyle H_{\text{T},s\tau_{3}-N-s}=$
$\displaystyle\sum_{\mathbf{k}\mathbf{p}\sigma}\bigg{(}T^{(1)}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{R},\mathbf{k}\sigma}c_{\text{L},1\mathbf{p}\sigma}+T^{(2)}_{\mathbf{k}\mathbf{p}}c^{{\dagger}}_{\text{R},\mathbf{k}\sigma}c_{\text{L},2\mathbf{p}\sigma}$
$\displaystyle+$ $\displaystyle\text{H.c.}\bigg{)},$ (127)

where the tunneling matrix elements are labeled by a single orbital index. We
shall neglect the spin indices for simplicity.

The total currents out of the L lead read

$\displaystyle\dot{I}_{\text{Tot},s\tau_{3}-N-s\tau_{3}}=$ $\displaystyle-
ie\left(\dot{N}_{\text{L},1}+\dot{N}_{\text{L},2}\right)$ (128)
$\displaystyle=$ $\displaystyle-
ie\sum_{\textbf{kp}}\sum_{ij}\left(T^{ij}_{kp}c^{{\dagger}}_{\text{R},i\mathbf{k}}c_{\text{L},j\mathbf{p}}-\text{H.c.}\right)$
(129)

$\displaystyle\dot{I}_{\text{Tot},s\tau_{3}-N-s}=$ $\displaystyle-
ie\left(\dot{N}_{\text{L},1}+\dot{N}_{\text{L},2}\right)$ (130)
$\displaystyle=$ $\displaystyle-
ie\sum_{\textbf{kp}}\sum_{i}\left(T^{(i)}_{kp}c^{{\dagger}}_{\text{R},\mathbf{k}}c_{\text{L},j\mathbf{p}}-\text{H.c.}\right),$
(131)

where $N_{\text{L},i}$ is the total charge associated with either orbital in
the L lead.

As shown in Ref. Mahan (2000), the Josephson current to leading order in the
tunneling matrix elements is determined from

$\displaystyle I_{\text{J}}(t)=$ $\displaystyle
2e\text{Im}\left[e^{-2ieVt/\hbar}\theta(eU)\right],$ (132)

where $eU$ is a potential drop across the junction.

###  F.1 $s\tau_{4}-N-s\tau_{3}$ and $s\tau_{3}-N-s$ junction

For these cases the quantity $\theta$ is obtained via analytical continuation
from

$\displaystyle\theta_{s\tau_{3}-N-s\tau_{3}}(eU)=$
$\displaystyle\lim_{i\omega\rightarrow
eU+i\eta}\left[2\sum_{ij}\sum_{mn}\sum_{\mathbf{kp}}T^{ij}_{\mathbf{kp}}T^{mn}_{\mathbf{-k,-p}}\frac{1}{\beta}\sum_{i\nu}F^{{\dagger}}_{\text{R},im}(\mathbf{k},i\nu)F_{\text{L},jn}(\mathbf{p},i\nu-i\omega)\right]$
(133)

$\displaystyle\theta_{s\tau_{3}-N-s}(eU)=$
$\displaystyle\lim_{i\omega\rightarrow
eU+i\eta}\left[2\sum_{ij}\sum_{\mathbf{kp}}T^{(i)}_{\mathbf{kp}}T^{(j)}_{\mathbf{-k,-p}}\frac{1}{\beta}\sum_{i\nu}F^{{\dagger}}_{\text{R}}(\mathbf{k},\nu)F_{\text{L},ij}(\mathbf{p},i\nu-i\omega)\right],$
(134)

where the anomalous Green’s functions are determined from the corresponding L
and R lead Hamiltonians as

$\displaystyle F_{\text{L/R},11}(\mathbf{k},i\nu)=$
$\displaystyle\frac{-\Delta_{\text{L/R}}\left[(i\nu)^{2}-\left(\xi_{0}-\xi_{3}\right)^{2}-|\Delta_{\text{L/R}}|^{2}+\xi^{2}_{1}\right]}{\Gamma(\mathbf{k},i\nu)}$
(135)

$\displaystyle F_{\text{L/R},12}(\mathbf{k},i\nu)=$
$\displaystyle\frac{-2\xi_{1}\Delta_{\text{L/R}}\left(i\nu+\xi_{3}\right)}{\Gamma(\mathbf{k},i\nu)}$
(136)

$\displaystyle F_{\text{L/R},21}(\mathbf{k},i\nu)=$
$\displaystyle\frac{2\xi_{1}\Delta_{\text{L/R}}\left(i\nu-\xi_{3}\right)}{\Gamma(\mathbf{k},i\nu)}$
(137)

$\displaystyle F_{\text{L/R},22}(\mathbf{k},i\nu)=$
$\displaystyle\frac{\Delta_{\text{L/R}}\left[(i\nu)^{2}-\left(\xi_{0}+\xi_{3}\right)^{2}-|\Delta_{\text{L/R}}|^{2}+\xi^{2}_{1}\right]}{\Gamma(\mathbf{k},i\nu)}$
(138)

for $s\tau_{3}$ leads, where

$\displaystyle\Gamma(\mathbf{k},i\nu)=$ $\displaystyle\left(i\nu-
E_{1}\right)\left(i\nu+E_{1}\right)\left(i\nu-
E_{2}\right)\left(i\nu+E_{2}\right)$ (139)

and

$\displaystyle
E_{1,2}=\sqrt{\xi^{2}_{0}+\xi^{2}_{3}+\xi^{2}_{1}+|\Delta|^{2}\pm
2\sqrt{\xi^{2}_{0}\left(\xi^{2}_{1}+\xi^{2}_{3}\right)+\xi^{2}_{1}|\Delta|^{2}}}$
(140)

are the BdG bands corresponding to the eigenvalues of the Hamiltonian in Eq.
5. Note that all $\xi$ and $E$ terms are functions of $\mathbf{k}$. Since we
consider $|\Delta_{\text{L}}|=|\Delta_{\text{R}}|$, $\Gamma$’s are independent
of the lead index.

For the single-channel $s$-wave lead, we have the standard expression

$\displaystyle F_{\text{R}}(\mathbf{k},\nu)=$
$\displaystyle\frac{-\Delta_{\text{R}}}{\left[i\nu-\sqrt{\xi^{2}_{0}+|\Delta|^{2}}\right]\left[i\nu+\sqrt{\xi^{2}_{0}+|\Delta|^{2}}\right]}$
(141)

These expressions simplify considerably if we consider the case where all of
the tunneling coefficients are identical for all channels as in
$T^{ij}_{\mathbf{kp}}=T_{\mathbf{kp}}$ and
$T^{(i)}_{\mathbf{kp}}=T_{\mathbf{kp}}$. These correspond to junctions and the
bound state spectra discussed in sections III and IV. The expressions for the
$\theta$ functions are

$\displaystyle\theta_{s\tau_{3}-N-s\tau_{3}}(eU)=$
$\displaystyle\lim_{i\omega\rightarrow
eU+i\eta}2\sum_{\mathbf{kp}}T_{\mathbf{kp}}T_{\mathbf{-k,-p}}\frac{1}{\beta}\sum_{i\nu}\left[\sum_{im}F^{{\dagger}}_{\text{R},im}(\mathbf{k},i\nu)\right]\left[\sum_{jn}F_{\text{L},jn}(\mathbf{p},i\nu-i\omega)\right]$
(142)

$\displaystyle\theta_{s\tau_{3}-N-s}(eU)=$
$\displaystyle\lim_{i\omega\rightarrow
eU+i\eta}2\sum_{\mathbf{kp}}T_{\mathbf{kp}}T_{\mathbf{-k,-p}}\frac{1}{\beta}\sum_{i\nu}F^{{\dagger}}_{\text{R}}(\mathbf{k},\nu)\left[\sum_{ij}F_{\text{L},ij}(\mathbf{p},i\nu-i\omega)\right],$
(143)

The sums in the square brakets amount to

$\displaystyle\left[\sum_{ij}F_{\text{L/R},ij}(\mathbf{k},i\nu)\right]=$
$\displaystyle\frac{-4\Delta_{\text{L/R}}\xi_{3}\left(\xi_{0}+\xi_{1}\right)}{\Gamma(\mathbf{k},i\nu)}$
(144)

From this expression, it is clear that, in the limit of zero intra-orbital
hybridization $\xi_{3}=0$, the _total_ Josephson currents in the tunneling
limit for either $s\tau_{3}-N-s\tau_{3}$ or $s\tau_{3}-N-s$ junctions vanish
identically. This holds even for a finite inter-orbital hybridization
$\xi_{1}$ on the L (and R lead if appropriate). The immediate reason for this
remarkable result is that the L lead continues to preserve the orbital
exchange symmetry defined in Sec. II.3 _in the bulk_ which ensures that

$\displaystyle F_{\text{L},11}(\mathbf{k},i\nu)=$ $\displaystyle-
F_{\text{L},22}(\mathbf{k},i\nu)$ (145) $\displaystyle
F_{\text{L},12}(\mathbf{k},i\nu)=$ $\displaystyle-
F_{\text{L},21}(\mathbf{k},i\nu)$ (146)

as long as the intra-orbital hybridization terms $\xi_{3}$ vanish, even though
the complete Hamiltonian does not preserve the symmetry for the
$s\tau_{3}-N-s$ junction. In the limit considered here, where all of the
orbitals couple identically across the junction, this symmetry in the bulk of
the L lead ensures that the total current out of the L lead vanishes.

Josephson currents which vanish in the limit of zero intra-orbital
hybridization terms are completely consistent with the results of sections III
and IV where the bound state spectra were shown to be insensitive to the
relative phase precisely when the same intra-orbital hybridization terms were
ignored.

###  F.2 $s\tau_{0}-N-s\tau_{0}$ and $s\tau_{0}-N-s$ junction

The expression for the current and the associated $\theta$ functions in Eqs.
142 is the same as for the $s\tau_{3}$ cases. However, the anomalous Green’s
functions for $s\tau_{0}$ are

$\displaystyle F_{\text{L/R},11}(\mathbf{k},i\nu)=$
$\displaystyle\frac{\Delta\left[-(i\omega)^{2}+\left(\xi_{0}-\xi_{3}\right)^{2}_{2}+\xi^{2}_{1}+|\Delta|^{2}\right]}{\chi(\textbf{k},i\nu)}$
(147)

$\displaystyle F_{\text{L/R},12}(\mathbf{k},i\nu)=$
$\displaystyle\frac{2\Delta\xi_{1}\xi_{0}}{\chi(\textbf{k},i\nu)}$ (148)

$\displaystyle F_{\text{L/R},21}(\mathbf{k},i\nu)=$
$\displaystyle\frac{2\Delta\xi_{1}\xi_{0}}{\chi(\textbf{k},i\nu)}$ (149)

$\displaystyle F_{\text{L/R},22}(\mathbf{k},i\nu)=$
$\displaystyle\frac{\Delta\left[-(i\omega)^{2}+\epsilon^{2}_{1}+\xi^{2}_{1}+|\Delta|^{2}\right]}{\chi(\textbf{k},i\nu)}$
(150)

where

$\displaystyle\chi(\textbf{k},i\nu)=\left(i\nu-
E_{1}\right)\left(i\nu+E_{1}\right)\left(i\nu-E_{2}\right)\left(i\nu-
E_{2}\right)$ (151)

and

$\displaystyle E_{1/2}=$
$\displaystyle\sqrt{\left(\xi_{0}\pm\sqrt{\xi^{2}_{1}+\xi^{2}_{3}}\right)^{2}+|\Delta|^{2}}$
(152)

is the bulk BdG spectrum for $s\tau_{0}$ pairing.

We calculate

$\displaystyle\left[\sum_{ij}F_{\text{L/R},ij}(\mathbf{k},i\nu)\right]=$
$\displaystyle\frac{-2\Delta\left\\{(i\omega)^{2}-\left[(\xi_{0}+\xi_{1})^{2}+\xi^{2}_{3}+|\Delta|^{2}\right]\right\\}}{\chi(\textbf{k},i\nu)}$
(153)

This is clearly distinct from the expression for $s\tau_{3}$ pairing in Eq.
144. Indeed, in the $\xi_{1}=\xi_{3}$ limit corresponding to degenerate
orbitals this reduces to

$\displaystyle\left[\sum_{ij}F_{\text{L/R},ij}(\mathbf{k},i\nu)\right]=$
$\displaystyle\frac{-2\Delta}{\left[(i\omega)^{2}-\xi^{2}_{0}-|\Delta|^{2}\right]}$
(154)

which amounts to twice the result for a single-channel junction, due to the
presence of the two, decoupled orbitals. This is consistent with the
dependence of the bound state spectrum on $\phi$ (Appendix C).

## References

  * Lee (2017) D.-H. Lee, Science 357, 32 (2017).
  * Mou  _et al._ (2011) D. Mou, S. Liu, X. Jia, J. He, Y. Peng, L. Zhao, L. Yu, G. Liu, S. He, X. Dong, J. Zhang, H. Wang, C. Dong, M. Fang, X. Wang, Q. Peng, Z. Wang, S. Zhang, F. Yang, Z. Xu, C. Chen, and X. J. Zhou, Phys. Rev. Lett. 106, 107001 (2011).
  * X.-P. Wang and T. Qian and P. Richard and P. Zhang and J. Dong and H.-D. Wang and C.-H. Dong and M.-H. Fang and H. Ding (2011) X.-P. Wang and T. Qian and P. Richard and P. Zhang and J. Dong and H.-D. Wang and C.-H. Dong and M.-H. Fang and H. Ding, Europhy. Lett. 93, 57001 (2011).
  * Xu  _et al._ (2012) M. Xu, Q. Q. Ge, R. Peng, Z. R. Ye, J. Jiang, F. Chen, X. P. Shen, B. P. Xie, Y. Zhang, A. F. Wang, X. F. Wang, X. H. Chen, and D. L. Feng, Phys. Rev. B 85, 220504 (2012).
  * Wang  _et al._ (2012) X.-P. Wang, P. Richard, X. Shi, A. Roekeghem, Y.-B. Huang, E. Razzoli, T. Qian, E. Rienks, S. Thirupathaiah, H.-D. Wang, C.-H. Dong, M.-H. Fang, M. Shi, and H. Ding, Europhys. Lett. 99, 67001 (2012).
  * Park  _et al._ (2011) J. T. Park, G. Friemel, Y. Li, J.-H. Kim, V. Tsurkan, J. Deisenhofer, H.-A. Krug von Nidda, A. Loidl, A. Ivanov, B. Keimer, and D. S. Inosov, Phys. Rev. Lett. 107, 177005 (2011).
  * Friemel  _et al._ (2012) G. Friemel, J. T. Park, T. A. Maier, V. Tsurkan, Y. Li, J. Deisenhofer, H.-A. Krug von Nidda, A. Loidl, A. Ivanov, B. Keimer, and D. S. Inosov, Phys. Rev. B 85, 140511 (2012).
  * Eschrig (2006) M. Eschrig, Adv. Phys. 55, 47 (2006).
  * Stockert  _et al._ (2011) O. Stockert, J. Arndt, E. Faulhaber, C. Geibel, H. S. Jeevan, S. Kirchner, M. Loewenhaupt, K. Schmalzl, W. Schmidt, Q. Si, and F. Steglich, Nat. Phys. 7, 119 (2011).
  * Maier  _et al._ (2011) T. A. Maier, S. Graser, P. J. Hirschfeld, and D. J. Scalapino, Phys. Rev. B 83, 100515 (2011).
  * Dai (2015) P. Dai, Rev. Mod. Phys. 87, 855 (2015).
  * Si  _et al._ (2016) Q. Si, R. Yu, and E. Abrahams, Nature Rev. Mater. 1, 16017 (2016).
  * Nica  _et al._ (2017) E. M. Nica, R. Yu, and Q. Si, npj Quantum Materials 2, 24 (2017), arXiv:1505.04170 .
  * Nica and Si (2021) E. M. Nica and Q. Si, npj Quantum Materials 6, 3 (2021).
  * Steglich  _et al._ (1979) F. Steglich, J. Aarts, C. D. Bredl, W. Lieke, D. Meschede, W. Franz, and H. Schäfer, Phys. Rev. Lett. 43, 1892 (1979).
  * Kittaka  _et al._ (2014) S. Kittaka, Y. Aoki, Y. Shimura, T. Sakakibara, S. Seiro, C. Geibel, F. Steglich, H. Ikeda, and K. Machida, Phys. Rev. Lett. 112, 067002 (2014).
  * Pang  _et al._ (2018) G. Pang, M. Smidman, J. Zhang, L. Jiao, Z. Weng, E. M. Nica, Y. Chen, W. Jiang, Y. Zhang, W. Xie, H. S. Jeevan, H. Lee, P. Gegenwart, F. Steglich, Q. Si, and H. Yuan, Proc. Nat. Acad. Sci. 115, 5343 (2018).
  * Yamashita  _et al._ (2017) T. Yamashita, T. Takenaka, Y. Tokiwa, J. A. Wilcox, Y. Mizukami, D. Terazawa, Y. Kasahara, S. Kittaka, T. Sakakibara, M. Konczykowski, S. Seiro, H. S. Jeevan, C. Geibel, C. Putzke, T. Onishi, H. Ikeda, A. Carrington, T. Shibauchi, and Y. Matsuda, Sci. Adv. 3, e1601667 (2017).
  * Yu  _et al._ (2014) R. Yu, J.-X. Zhu, and Q. Si, Phys. Rev. B 89, 024509 (2014).
  * Yin  _et al._ (2014) Z. P. Yin, K. Haule, and G. Kotliar, Nat. Phys. 10, 845 (2014).
  * Ong  _et al._ (2016) T. Ong, P. Coleman, and J. Schmalian, Proc. Nat. Acad. Sci. 113, 5486 (2016).
  * Sprau  _et al._ (2017) P. O. Sprau, A. Kostin, A. Kreisel, A. E. Böhmer, V. Taufour, P. C. Canfield, S. Mukherjee, P. J. Hirschfeld, B. M. Andersen, and J. C. S. Davis, Science 357, 75 (2017).
  * Hu  _et al._ (2018) H. Hu, R. Yu, E. M. Nica, J.-X. Zhu, and Q. Si, Phys. Rev. B 98, 220503 (2018).
  * Yu  _et al._ (2018) R. Yu, J.-X. Zhu, and Q. Si, Phys. Rev. Lett. 121, 227003 (2018).
  * Raghu  _et al._ (2008) S. Raghu, X.-L. Qi, C.-X. Liu, D. J. Scalapino, and S.-C. Zhang, Phys. Rev. B 77, 220503 (2008).
  * Si and Abrahams (2008) Q. Si and E. Abrahams, Phys. Rev. Lett. 101, 076401 (2008).
  * Daghofer  _et al._ (2010) M. Daghofer, A. Nicholson, A. Moreo, and E. Dagotto, Phys. Rev. B 81, 014511 (2010).
  * Sauls (2018) J. A. Sauls, Philos. Trans. A Math. Phys. Eng. Sci. 376, 20180140 (2018).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * Mahan (2000) G. D. Mahan,  _Many-particle Physics_ (Plenum, New York, 2000).