# Green’s functions of quasi-one-dimensional layered systems and their
application to Josephson junctions

Kiryl Piasotski<EMAIL_ADDRESS>Institut für Theorie der
Kondensierten Materie, Karlsruher Institut für Technologie, 76131 Karlsruhe,
Germany  Institut für Quanten Materialien und Technologien, Karlsruher
Institut für Technologie, 76021 Karlsruhe, Germany  Mikhail Pletyukhov
Institut für Theorie der Statistischen Physik, RWTH Aachen, 52074 Aachen,
Germany  Alexander Shnirman  Institut für Theorie der Kondensierten Materie,
Karlsruher Institut für Technologie, 76131 Karlsruhe, Germany  Institut für
Quanten Materialien und Technologien, Karlsruher Institut für Technologie,
76021 Karlsruhe, Germany

###### Abstract

We develop Green’s function formalism to describe continuous multi-layered
quasi-one-dimensional setups described by piece-wise constant single-particle
Hamiltonians. The Hamiltonians of the individual layers are assumed to be
quadratic polynomials in the momentum operator with matrix-valued
(multichannel) coefficients. This, in particular, allows one to study
transport in heterostructures consisting of multichannel conducting,
superconducting, or insulating components with band structures of arbitrary
complexity. We find a general expression for the single-particle Green’s
function of the combined setup in terms of the bulk (translationally
invariant) Green’s functions of its constituents. Furthermore, we provide the
expression for the global density of states of the combined system and
establish the bound state equation in terms of bulk Green’s functions. We
apply our formalism to investigate the spectrum and current-phase relations in
ordinary and topological Josephson junctions, additionally showing how to
account for the effects of static disorder and local Coulomb interaction.

##  I Introduction

Figure 1: Sketch of the model: A heterostructure consists of $M+1$ layers
described by individual matrix-valued Hamiltonians $H_{m}$, which are
consequently coupled with each other across matrix-valued tunneling barriers
$\mathcal{U}_{m}$.

Description of equilibrium and transport properties of layered quantum systems
is a common problem in the domains of quantum electronics and solid state
theoryDatta (1997); Freericks (2016). Albeit typically a single-particle
quantum mechanical problem, its solution is rarely simple due to the intricate
band structures of the materials forming the layers. Most commonly, nowadays,
these systems are analyzed numerically within the tight-binding approximation
(like e.g. in Ref. [Groth  _et al._ , 2014]), allowing one to get a good grip
on the low-lying excitations, as well as to assess the effects of the static
disorder. The main idea behind this approach consistsAndo (1991) in numerical
studies of the ballistic conductance in mesoscopic structures within a lattice
model, being expressed in terms of the lattice Green’s function. It can be
equally well applied to the study of equilibrium properties like Josephson
current (JC) in Josephson junctions (JJ) of arbitrary width, see e.g. Refs.
[Cheng and Lutchyn, 2012,Yokoyama  _et al._ , 2014].

An alternative to the tight-binding numerics is the scattering matrix approach
Landauer (1957, 1970) applied to continuum ballistic models. One of its key
early-day achievements was the calculation of the two-terminal conductance in
terms of the transmission probabilities. The microscopic justification of this
method relies either on taking the continuum limit of the wavefunction
matching (WFM) of Ref. [Ando, 1991] or on studying the continuum limit of the
atomistic Green’s function (AGF) as in Ref. [Caroli  _et al._ , 1971]. The
equivalence of both WMF and AGF approaches has been fully substantiated in
Ref. [Khomyakov  _et al._ , 2005]. The relation between transmission and
reflection coefficients of the scattering matrix on the one hand and Green’s
function, on the other hand, is widely known as the Fisher-Lee relationFisher
and Lee (1981). Its mode-resolved generalization has been recently proposed in
Ref. [Boumrar  _et al._ , 2020].

In application to superconducting systems, the scattering matrix approach has
been extended by Blonder, Tinkham, and KlapwijkBlonder  _et al._ (1982) (BTK)
on the basis of solutions of the Bogoliubov-de-Gennes equations. The BTK
theory has been further generalized by C. Beenakker to the multichannel case
in Ref. [Beenakker, 1991]. In this work the compact equation for the sub-gap
Andreev bound statesAndreev (1964, 1966) (ABS) and the expression for the
continuous excitation spectrum of the JJ have been established in terms of
normal and Andreev scattering matrices. In practice, however, these matrices
are often treated in the so-called Andreev limit $\Delta_{0}\ll\mu$, where
$\Delta_{0}$ is the absolute value of the superconducting order parameter and
$\mu$ is the chemical potential, essentially neglecting the normal reflection
at the superconducting interface.

An alternative theoretical description of the superconducting tunneling and
proximity effects has been developed by G. Arnold in Refs. [Arnold,
1978,Arnold, 1985] using standard nonequilibrium Green’s functions. His
approach is based on the theory of FeuchtwangFeuchtwang (1974a, b, 1975),
which does not make use of the tunneling Hamiltonian. This theory shares many
common features with the AGF method of Caroli et al [Caroli  _et al._ , 1971]
mentioned above, and these techniques will serve as a starting point for our
present consideration.

Discussing the approaches to interface physics, one cannot but mention that
quasiclassical Green’s functions treated in terms of the Eilenberger and the
Usadel equations are also traditionally used for describing superconducting
proximity effects and the JJs, see Refs. [Golubov  _et al._ , 2004] and
[Eschrig, 2018] for reviews.

The goal of the present manuscript is to generalize and further develop the
Green’s function techniques of Arnold, Feuchtwang, Caroli et al, mentioned
above. The structure of the manuscript is twofold.

First of all, we provide multichannel and multi-interface generalizations of
results in Refs. [Caroli  _et al._ , 1971, Arnold, 1978,Arnold, 1985,
Feuchtwang, 1974a, Feuchtwang, 1974b, Feuchtwang, 1975], essentially building
upon the method of Ref. [Caroli  _et al._ , 1971]. In particular, we derive a
compact closed-form expression for the composite Green’s function
$G(x,x^{\prime})$ of a quasi-one-dimensional heterostructure (sketched in Fig.
1) given merely in terms of the bulk Green’s functions
$G^{(0,m)}(x,x^{\prime})$ of its constituting layers (labeled by
$m=0,\ldots,M$). The knowledge of $G(x,x^{\prime})$ is beneficial for several
reasons: It gives direct access to the system’s spectrum, allows for a
calculation of various physical observables (like e.g. the density of states
(DOS) and the current), may serve as a starting point for low-energy
approximations (like e.g. the Andreev limit), as well as provides essential
input for a perturbative diagrammatic treatment of disorder and many-body
interaction effects. The bulk single-layer Green’s functions
$G^{(0,m)}(x,x^{\prime})$ are obtained by means of the standard Fourier
transformation and are computationally inexpensive. For few-channel models
$G^{(0,m)}(x,x^{\prime})$ may often be calculated analytically, as
demonstrated in our examples below.

We shall point out that ideas of that sort, namely the extraction of the
properties of inhomogeneous quantum systems from their bulk counterpart, are
traditional in solid-state theory. This subject has a long history, and it
often bears different names: the quantum theory of surface statesGarcía-
Moliner and Rubio (1969); García-Moliner  _et al._ (1969); García-Moliner and
Velasco (1991, 1994); Rodriguez-Coppola  _et al._ (1990) or the method of
embeddingInglesfield (1971, 1981); Wortmann  _et al._ (2002), for example.

Secondly, we apply the derived expression for $G(x,x^{\prime})$ to the two
paradigmatic examples: (1) Josephson model of a tunneling barrier between two
$s$-wave superconductors, and (2) the model of a JJvan Heck  _et al._ (2017);
Murthy  _et al._ (2020) between two semiconducting wires with strong spin-
orbit interaction and proximity induced superconductivity, submersed into the
parallel magnetic field (each of which is capable of hosting a Majorana zero
modeLutchyn  _et al._ (2010); Oreg  _et al._ (2010)). In particular, we
demonstrate that both the excitation spectrum (i.e. the global DOS) as well as
the JC can be obtained in terms of a single matrix
$d=\begin{pmatrix}G(0,0)&G(0,W)\\\ G(W,0)&G(W,W)\end{pmatrix}^{-1}$ for the
finite-width junction (with the two interfaces at $x=0$ and $x=W$),
degenerating to $d=[G(0,0)]^{-1}$ in the short junction limit (single
interface at $x=0$). In particular, the correction to the global DOS due to
the tunneling between the layers is given by
$\delta\rho=-\frac{1}{\pi}\frac{\partial}{\partial\omega}\text{Im}\,\ln\det
d$, which resembles analogous expressions for $\delta\rho$ in terms of the
scattering matrixAkkermans  _et al._ (1991); Souma and Suzuki (2002). A
general equation for the bound states acquires a particularly simple form
$\det d=0$, which is equally well applicable to all types of heterostructures.
The relation of our equations for the bound states and for $\delta\rho$ to
their scattering matrix analogs becomes transparent on the basis of the
Fisher-Lee relationFisher and Lee (1981): The matrix $d^{-1}$ appears as a
common factor in expressions for all components of the scattering matrix, and
it can be generally interpreted as a core part of the $T$-matrix in the
coordinate representation. The advantage of computing $d$ directly from the
Green’s function $G(x,x^{\prime})$ consists, though, in its additive form,
combining contributions from adjacent layers (and possibly from a local
contact potential at their interface). This property of $d$ is thus analogous
to that of self-energies.

Using our formulas, we study the excitation spectrum and the JC in various
parametric regimes of the two models discussed above, relaxing the Andreev
approximation (which might be essential for one-dimensional wires) and
allowing for arbitrary values of the junction’s width $W$. The knowledge of
the explicit form of $G(x,x^{\prime})$ allows us to account for the effects of
random potential disorder and local Coulomb interaction at the interface of
the JJ. We find that the effect of static disorder arises only beyond the
Andreev limit. Taking into account the local Coulomb interaction, we predict a
crossover between the $0$\- and the $\pi$-junctions.

##  II Formalism of interface Green’s functions

###  II.1 Model and problem formulation

Let us consider a quasi-one-dimensional system that consists of $M+1$ layers
(labeled by $m=0,\ldots,M$) subsequently connected with each other across
tunneling barriers (also labeled by $m=1,\ldots,M$) along the spatial
$x$-axis, see Fig. 1. We further assume that $m$-th layer extends over the
spatial range $(x_{m},x_{m+1})$, with $x_{1}$ and $x_{M}$ being the
coordinates of the leftmost and the rightmost interfaces, respectively. The
leftmost boundary of the whole system is $x_{0}$ and the rightmost one is
$x_{M+1}$. An important class of models with semi-infinite leads is also
included in the present consideration: They are realized by setting
$x_{0}=-\infty$ and $x_{M+1}=+\infty$. For the $m$-th subsystem we assume the
Hamiltonian of the form

$\displaystyle
H_{m}=\frac{1}{2}\mathcal{A}_{m}p^{2}+\mathcal{B}_{m}p+\mathcal{C}_{m},$ (1)

where $\mathcal{A}_{m},\ \mathcal{B}_{m},\ \mathcal{C}_{m}$ are constant
$N_{c}\times N_{c}$ hermitian matrices, parametrically dependent on the
conserved transverse (quasi)momenta $k_{y}$ and $k_{z}$, with $N_{c}$ being
the number of channels (e.g. spin and/or orbital quantum numbers, etc.), and
$p\equiv p_{x}=-i\partial_{x}$ is the momentum operator in the direction along
the system.

The tunneling barriers between the subsystems are modeled by delta-like
potentials at the corresponding interfaces:

$\displaystyle\mathcal{U}(x)=\sum_{m=1}^{M}\mathcal{U}_{m}\delta(x-x_{m}),$
(2)

where $\mathcal{U}_{m}$ are generally considered as $N_{c}\times N_{c}$
hermitian matrices (which, in particular, allows us to treat magnetically
active contacts as well).

The Hamiltonian of the combined system is written as

$\displaystyle
H=\frac{1}{2}p\mathcal{A}(x)p+\frac{1}{2}\left\\{\mathcal{B}(x),\
p\right\\}+\mathcal{C}(x)+\mathcal{U}(x),$ (3)

where

$\displaystyle\begin{Bmatrix}\mathcal{A}(x)\\\ \mathcal{B}(x)\\\
\mathcal{C}(x)\end{Bmatrix}=\sum_{m=0}^{M}\Theta(x_{m}<x<x_{m+1})\begin{Bmatrix}\mathcal{A}_{m}\\\
\mathcal{B}_{m}\\\ \mathcal{C}_{m}\end{Bmatrix}.$ (4)

The eigenvalue problem for the Hamiltonian (3) is complemented by matching
conditions for the wave function and its derivative at each interface. While
the wavefunction is continuous,

$\displaystyle\Psi(x_{m}^{-})=\Psi(x_{m}^{+})\equiv\Psi(x_{m}),\quad
x_{m}^{\pm}=x_{m}\pm 0^{+},$ (5)

its derivative satisfies the following matching condition at the $m$-th
interface:

$\displaystyle\frac{\mathcal{A}_{m}\Psi^{\prime}(x_{m}^{+})-\mathcal{A}_{m-1}\Psi^{\prime}(x_{m}^{-})}{2}=\left[\mathcal{U}_{m}-i\frac{\mathcal{B}_{m}-\mathcal{B}_{m-1}}{2}\right]\Psi(x_{m}).$
(6)

This condition is traditionally derived by integrating the Schrödinger
equation $H\Psi(x)=E\Psi(x)$ over the infinitesimally small region
$(x_{m}^{-},x_{m}^{+})$ enclosing the contact’s coordinate $x_{m}$. Physically
it expresses the conservation of the current density across the interface.

The main goal of our present consideration is to find an expression for the
(retarded) Green’s function of the whole system111Note that the
heterostructural Green’s function $G(x,x^{\prime})$ is a highly non-trivial
object, containing information about all possible multi-barrier scattering
phenomena and interface-localized bound states., which satisfies the equation

$\displaystyle[z-H]G(x,x^{\prime};z)=\delta(x-x^{\prime})$ (7)

for arbitrary complex-valued spectral parameter $z$ (with $\text{Im}\,z>0$, in
particular for $z=\omega+i0^{+}$, in the case of the retarded function). The
differential equation (7) is complemented by the vanishing boundary conditions
$G(x_{0},x^{\prime};z)=G(x_{M+1},x^{\prime};z)=0$ at the system’s ends. By the
virtue of the Lehmann representation, the Green’s function $G(x,x^{\prime};z)$
must also obey the matching conditions (5) and (6) in the variable $x$. We
note that $G(x,x^{\prime};z)$ can be alternatively defined as a solution of a
reciprocal differential equation with respect to the variable $x^{\prime}$
equipped with the boundary conditions $G(x,x_{0};z)=G(x,x_{M+1};z)=0$, while
the corresponding interface matching conditions are obtained by hermitian
conjugating (5) and (6).

As an input we use the translationally invariant (aka bulk) Green’s functions
$G^{(0,m)}(x,x^{\prime};z)$ satisfying the equation

$\displaystyle[z-H_{m}]G^{(0,m)}(x,x^{\prime};z)=\delta(x-x^{\prime})$ (8)

on the whole spatial axis. They are evaluated by means of the Fourier
transformation

$\displaystyle
G^{(0,m)}(x,x^{\prime};z)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}\frac{e^{ik(x-x^{\prime})}}{z-h_{m}(k)},$
(9)

where $h_{m}(k)$ is the Hamiltonian obtained from $H_{m}$ in (1) by the
substitution $p\to k$. Often the integral in (9) may be evaluated analytically
by means of the integration in the complex $k$-plane: To this end, one needs
to establish the complex-valued roots $k(z)$ of the secular equation
$\det[z-h_{m}(k)]=0$ with positive imaginary parts, $\text{Im}\,k(z)>0$. Then
the result of integration in (9) is represented as a sum of residua at the
corresponding poles $k(z)$ (see Appendix A.1 for details). For few-channel
models, the equation $\det[z-h_{m}(k)]=0$ admits analytical solutions (some of
them will be demonstrated in the following examples), while in general the
roots $k(z)$ have to be determined numerically. As it will be demonstrated
later, the root-searching routine is the only numerical part in solving the
problem of establishing $G(x,x^{\prime};z)$, and for this reason, the Green’s
function approach developed below should provide a considerable speed up in
the study of arbitrary heterostructures.

An opening move in establishing $G(x,x^{\prime};z)$ is a determination of the
set of Green’s functions $G^{m}(x,x^{\prime};z)$ describing every isolated
layer $m$ on the corresponding spatial interval $(x_{m},x_{m+1})$. The Green’s
function $G^{m}(x,x^{\prime};z)$ also satisfies the differential equation (8),
but – in contrast to $G^{(0,m)}(x,x^{\prime};z)$ – it vanishes at the
interval’s ends, that is
$G^{m}(x_{m},x^{\prime};z)=G^{m}(x_{m+1},x^{\prime};z)=0$.

Thankfully, the determination of $G^{m}(x,x^{\prime})$ (the argument $z$ is
omitted for brevity) is conveniently solved by the so-called boundary Green’s
function technique, allowing one to write a simple relation between the
propagators of infinite [$G^{(0,m)}(x,x^{\prime})$] and bounded
[$G^{m}(x,x^{\prime})$] systems (see Appendix A.2, and references therein, for
the summary of the key results).

Next, we find out how the isolated layers are coupled with each other. As
such, this coupling is dictated by the matching conditions (5) and (6). The
key question is then how does one implement these requirements in terms of a
local potential? Once this is accomplished, a set of Dyson’s equations may be
set up, relating the full Green’s function $G(x,x^{\prime})$ to those
$G^{m}(x,x^{\prime})$ of the individual layers. Furthermore, it is expected
that the set of these equations admits a closed-form analytical solution, as
by the locality of the coupling potential the integral equations are expected
to reduce to algebraic ones for a finite number of unknowns.

To give an answer to the key question we underscore the two approaches which
we find especially useful in practice.

The first approach, inspired by the earlier ideas of the seminal paper Ref.
[Caroli  _et al._ , 1971], is based on discretizing the equation (7) on a
lattice with the spacing $a$, solving the obtained tight-binding counterpart,
and taking carefully the continuum limit $a\to 0$. As in the tight-binding
description the space consists of discrete points, the tunneling between the
disjoint parts of the system may be defined unambiguously, and hence the
above-described program with a formal solution of the corresponding Dyson
equations can be successfully executed. Further evaluation of the limit $a\to
0$ allows us to recover the complete Green’s function of the corresponding
continuum theory, expressed in terms of boundary Green’s functions of the
individual layers as well as their spatial derivatives up to the second order.
In this paper, we describe the major steps of implementing this approach for
the models of our present interest.

The second approach is based on (the multichannel generalization of) the
Sturm-Liouville theory for second-order differential operators. It
constructively exploits the matching conditions (5) and (6) to recover
$G(x,x^{\prime})$ staying within the paradigm of continuum models. The results
of its application naturally reproduce those described below in the present
paper. Further details of the second approach will be reported elsewhere
[Piasotski  _et al._ , ].

###  II.2 Single barrier

In this subsection, we extrapolate the construction technique of the composite
Green’s function of a double-layer system of Ref. [Caroli  _et al._ , 1971] to
the multichannel case.

Let us start by considering the simplest case of a single barrier separating
left $L\ (m=0,\ x\in(x_{0},0),\ x_{0}<0)$ and right $R\ (m=1,\ x\in(0,x_{2}),\
x_{2}>0)$ subsystems at $x=x_{1}=0$.

We introduce the following lattice counterpart of the single-barrier continuum
model

$\displaystyle H$ $\displaystyle=H_{L}+H_{R}+|0\rangle W_{0}\langle 0|$ (10)
$\displaystyle-|0\rangle t_{L}\langle-1|-|-1\rangle t_{L}^{\dagger}\langle
0|-|1\rangle t_{R}\langle 0|-|0\rangle t_{R}^{\dagger}\langle 1|,$ (11)

where the left and right disjoint subsystems are defined on the lattice sites
$n_{0}\leq n\leq-1$ and $1\leq n\leq n_{2}$, respectively. They are described
by the Hamiltonians

$\displaystyle H_{L}=$ $\displaystyle-\sum_{n=n_{0}+1}^{-1}(|n\rangle
t_{L}\langle n-1|+|n-1\rangle t_{L}^{\dagger}\langle n|)$
$\displaystyle+\sum_{n=n_{0}}^{-1}|n\rangle W_{L}\langle n|,$ (12)
$\displaystyle H_{R}=$ $\displaystyle-\sum_{n=1}^{n_{2}-1}(|n+1\rangle
t_{R}\langle n|+|n\rangle t_{R}^{\dagger}\langle n+1|)$
$\displaystyle+\sum_{n=1}^{n_{2}}|n\rangle W_{R}\langle n|,$ (13)

with constant, matrix-valued (in the channel space) hopping amplitudes
$t_{L},t_{R}$ and onsite potentials $W_{L},W_{R}$. The central cite $n=0$ is
characterized by the onsite potential $W_{0}$. It is coupled to both
subsystems by the same nearest-neighbor hopping amplitudes $t_{L}$ and $t_{R}$
as occur in the Hamiltonians (12) and (13), respectively.

Treating the coupling term (11) as a perturbation, we set up the following
Dyson equation in the coordinate representation for the Green’s function
$G=\frac{1}{z-H}$ of the full system

$\displaystyle G_{n,n^{\prime}}$
$\displaystyle=G^{L}_{n,n^{\prime}}+G^{C}_{n,n^{\prime}}+G^{R}_{n,n^{\prime}}$
$\displaystyle-(G_{n,-1}^{L}t_{L}^{\dagger}+G_{n,1}^{R}t_{R})G_{0,n^{\prime}}$
$\displaystyle-\delta_{n,0}G_{0,0}^{C}(t_{L}G_{-1,n^{\prime}}+t_{R}^{\dagger}G_{1,n^{\prime}})$
(14)

in terms of the Green’s functions $G^{L}=\frac{1}{z-H_{L}}$,
$G^{C}=|0\rangle\frac{1}{z-W_{0}}\langle 0|$, $G^{R}=\frac{1}{z-H_{R}}$ of the
three disjoint subsystems. The Green’s function $G_{n,n^{\prime}}^{L}$ is
nonzero only for $n,n^{\prime}\leq-1$, while $G_{n,n^{\prime}}^{R}$ is nonzero
only for $n,n^{\prime}\geq 1$. In the corresponding domains, they are
expressed according to (149) in terms of the Green’s functions
$G^{(L,0)}_{n,n^{\prime}}$ and $G^{(R,0)}_{n,n^{\prime}}$ of the two auxiliary
models defined on the whole lattice and using the constant parameters from the
left and the right subsystems, respectively.

Due to the locality of the perturbation, the Dyson equation (14) admits the
explicit solution

$\displaystyle
G_{n,n^{\prime}}=G_{n,n^{\prime}}^{L}+G_{n,n^{\prime}}^{R}+F_{n}D^{-1}\bar{F}_{n^{\prime}},$
(15)

where

$\displaystyle
D=z-W_{0}-t_{L}G_{-1,-1}^{L}t_{L}^{\dagger}-t_{R}^{\dagger}G_{1,1}^{R}t_{R}$
(16)

and

$\displaystyle F_{n}$
$\displaystyle=-\delta_{n,0}+G_{n,-1}^{L}t_{L}^{\dagger}+G_{n,1}^{R}t_{R},$
(17) $\displaystyle\bar{F}_{n^{\prime}}$
$\displaystyle=-\delta_{0,n^{\prime}}+t_{L}G_{-1,n^{\prime}}^{L}+t_{R}^{\dagger}G_{1,n^{\prime}}^{R}.$
(18)

We particularly note that

$\displaystyle F_{0}$ $\displaystyle=\bar{F}_{0}=-1,$ (19) $\displaystyle
G_{0,0}$ $\displaystyle=D^{-1}.$ (20)

On the basis of the solution (15) it is straightforward to compute the global
DOS. Setting $z=\omega+i0^{+}$, we evaluate

$\displaystyle\sum_{n=n_{0}}^{n_{2}}\text{tr}\left[G_{n,n}-G_{n,n}^{L}-G_{n,n}^{R}\right]$
(21)
$\displaystyle=\text{tr}\left[\frac{1}{D}\right]+\sum_{n=n_{0}}^{-1}\text{tr}\left[\frac{1}{D}t_{L}G_{-1,n}^{L}G^{L}_{n,-1}t_{L}^{\dagger}\right]$
$\displaystyle+\sum_{n=1}^{n_{2}}\text{tr}\left[t_{R}^{\dagger}G_{1,n}^{R}G_{n,1}^{R}t_{R}\frac{1}{D}\right],$
(22)

where the trace operation is performed in the channel space. Using the
identities

$\displaystyle\sum_{n=n_{0}}^{-1}G_{-1,n}^{L}G^{L}_{n,-1}$
$\displaystyle=-\frac{\partial G_{-1,-1}^{L}}{\partial\omega},$ (23)
$\displaystyle\sum_{n=1}^{n_{2}}G_{1,n}^{R}G_{n,1}^{R}$
$\displaystyle=-\frac{\partial G_{1,1}^{R}}{\partial\omega},$ (24)

which are most easily proven in the Lehmann representation, and the Jacobi’s
formula $\text{tr}\,[D^{-1}\frac{\partial
D}{\partial\omega}]=\frac{\partial}{\partial\omega}\ln\det D$, we establish
that

$\displaystyle\rho(\omega)=\rho^{L}(\omega)+\rho^{R}(\omega)-\frac{1}{\pi}\text{Im}\,\frac{\partial}{\partial\omega}\ln\det
D(\omega+i0^{+}).$ (25)

The last term in this expression represents a correction to the global DOS due
to the tunneling between the subsystems. It has a form analogous to that of
the familiar expression in terms of the scattering matrixSouma and Suzuki
(2002). The terms $\rho^{L}(\omega)$ and $\rho^{R}(\omega)$ represent the
global DOS of the left and right disjoint subsystems, respectively.

In order to derive the continuum limit of (15) we make certain assumptions
about the scaling of the Hamiltonian parameters with the lattice constant $a$.
In particular, we define for $m=L,R$

$\displaystyle t_{m}+t_{m}^{\dagger}=\frac{\mathcal{A}_{m}}{a^{2}},$ (26)
$\displaystyle i(t_{m}-t_{m}^{\dagger})=\frac{\mathcal{B}_{m}}{a},$ (27)
$\displaystyle W_{m}=\frac{\mathcal{A}_{m}}{a^{2}}+\mathcal{C}_{m},$ (28)

where $\mathcal{A}_{m},\mathcal{B}_{m},\mathcal{C}_{m}=O(a^{0})$ are constant
matrices.

In addition, we choose

$\displaystyle W_{0}$
$\displaystyle=\frac{\mathcal{A}_{L}+\mathcal{A}_{R}}{2a^{2}}+\frac{\mathcal{U}_{1}}{a}.$
(29)

This choice of the leading $O(\frac{1}{a^{2}})$ term in (29) is importantLévy-
Leblond (1995) for ensuring the fulfillment of the matching conditions (5) and
(6) for the Green’s function in the continuum limit [see also in the end of
the subsection]. In turn, the subleading $O(\frac{1}{a})$ term induces the
impurity delta-potential (2), which is generally matrix-valued.

It is important to emphasize that our limiting procedure essentially differs
from the one frequently used in the wide-band limit treatment of the tunneling
regime $t^{\prime}_{L,R}\ll t_{L,R}$, where the hopping amplitudes $t_{L,R}$
in the leads considerably dominate over the hopping amplitudes
$t^{\prime}_{L,R}$ from the leads onto the quantum dot at the site $n=0$, and
the characteristic tunneling rates
$\Gamma_{L,R}=\pi(t^{\prime}_{L,R})^{2}/t_{L,R}$ giving rise to (the imaginary
part of) the dot’s self-energy are then much smaller than the bandwidths of
the leads $\sim t_{L,R}$. Recall that in our treatment we set
$t^{\prime}_{L,R}=t_{L,R}$.

Defining the continuous variable $x=na$ in the limit $a\to 0$, we recover the
continuum analog (1) of the lattice Hamiltonians (12) and (13). Relating the
continuum and the lattice Green’s functions via

$\displaystyle G(x,x^{\prime})=\lim_{a\to 0}\frac{1}{a}G_{n,n^{\prime}},$ (30)

we derive in Appendix B the continuum analog of (15). It reads

$\displaystyle G(x,x^{\prime})=$ $\displaystyle
G^{L}(x,x^{\prime})+G^{R}(x,x^{\prime})$
$\displaystyle+F(x)d^{-1}\bar{F}(x^{\prime}),$ (31)

where $G^{L}(x,x^{\prime})$ and $G^{R}(x,x^{\prime})$ are the boundary Green’s
functions of the corresponding disjoint regions $(x_{0},0)$ and $(0,x_{2})$
which identically vanish outside of them, respectively. The last term in Eq.
(31) contains non-diagonal contributions in the subsystem’s basis (i.e. it is
generically non-zero for all $x,\ x^{\prime}\in(x_{0},\ x_{2})$), and thereby
it mediates the coupling between the subsystems. It is defined through the
following objects

$\displaystyle F(x)=\lim_{a\to 0}F_{n}=$
$\displaystyle-\Theta(-x)\left[G_{2}^{(0,L)}(x,0^{+})+G_{2}^{L}(x,0^{-})-G_{2}^{(0,L)}(x,0^{-})\right]\frac{\mathcal{A}_{L}}{2}$
$\displaystyle+\Theta(x)\left[G_{2}^{(0,R)}(x,0^{-})+G_{2}^{R}(x,0^{+})-G_{2}^{(0,R)}(x,0^{+})\right]\frac{\mathcal{A}_{R}}{2},$
(32) $\displaystyle\bar{F}(x^{\prime})=\lim_{a\to 0}\bar{F}_{n^{\prime}}=$
$\displaystyle-\Theta(-x^{\prime})\frac{\mathcal{A}_{L}}{2}\left[G_{1}^{(0,L)}(0^{+},x^{\prime})+G_{1}^{L}(0^{-},x^{\prime})-G_{1}^{(0,L)}(0^{-},x^{\prime})\right]$
$\displaystyle+\Theta(x^{\prime})\frac{\mathcal{A}_{R}}{2}\left[G_{1}^{(0,R)}(0^{-},x^{\prime})+G_{1}^{R}(0^{+},x^{\prime})-G_{1}^{(0,R)}(0^{+},x^{\prime})\right],$
(33) $\displaystyle d=\lim_{a\to 0}a\,D=$
$\displaystyle-\mathcal{U}_{1}-\frac{1}{8}\mathcal{A}_{L}\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)\mathcal{A}_{L}-\frac{1}{8}\mathcal{A}_{R}\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{R}(x,x)\mathcal{A}_{R},$ (34)

Above we employed a set of compact notations for the coordinate derivatives of
Green’s functions

$\displaystyle G^{m}_{1}(x,x^{\prime})=\frac{\partial
G^{m}(x,x^{\prime})}{\partial x},\quad G^{m}_{2}(x,x^{\prime})=\frac{\partial
G^{m}(x,x^{\prime})}{\partial x^{\prime}},$ (35)

with analogous abbreviations for the translationally invariant propagators
$G^{(0,m)}(x,x^{\prime})$.

The obtained expressions represent a multichannel generalization of Eqs. (22)
and (30) in Ref. [Caroli  _et al._ , 1971]. It is also worth mentioning that
similar expressions for the two-channel case were previously derived in Refs.
[Arnold, 1978,Arnold, 1985].

The expression (31) can be interpreted in terms of the $T$-matrix. Formally
rewriting it as

$\displaystyle G(x,x^{\prime})$ $\displaystyle=G^{L+R}(x,x^{\prime})+\langle
x|G^{L+R}l^{\dagger}d^{-1}lG^{L+R}|x^{\prime}\rangle$ (36)
$\displaystyle=G^{L+R}(x,x^{\prime})$ $\displaystyle+\int dy\int
dy^{\prime}G^{L+R}(x,y)T(y,y^{\prime})G^{L+R}(y^{\prime},x^{\prime}),$ (37)

where $G^{L+R}=G^{L}+G^{R}$ and $l$ is the formally introduced boundary
hermitian operator in terms of

$\displaystyle\langle x|G^{L+R}l^{\dagger}|y\rangle$
$\displaystyle=i\delta(y)F(x),$ (38) $\displaystyle\langle
y^{\prime}|lG^{L+R}|x^{\prime}\rangle$
$\displaystyle=-i\delta(y^{\prime})\bar{F}(x^{\prime}),$ (39)

we identify $T(y,y^{\prime})=\langle y|l^{\dagger}d^{-1}l|y^{\prime}\rangle$
with the $T$-operator in the coordinate representation.

The functions (32) and (33) have the remarkable properties

$\displaystyle F(0^{+})=F(0^{-})\equiv F(0)$ $\displaystyle=-1,$ (40)
$\displaystyle\bar{F}(0^{+})=\bar{F}(0^{-})\equiv\bar{F}(0)$
$\displaystyle=-1,$ (41)

which follow from the standard jump conditions

$\displaystyle[G_{2}^{(0,m)}(0,0^{+})-G_{2}^{(0,m)}(0,0^{-})]\frac{\mathcal{A}_{m}}{2}$
$\displaystyle=1,$ (42)
$\displaystyle\frac{\mathcal{A}_{m}}{2}[G_{1}^{(0,m)}(0^{+},0)-G_{1}^{(0,m)}(0^{-},0)]$
$\displaystyle=1$ (43)

for the Green’s function derivatives. The expressions (40) and (41) thus
appear to be consistent with their lattice analogues (17) and (18),
respectively.

By the virtue of $G^{L}(0,0)=G^{R}(0,0)$ and (40), (41) we recover

$\displaystyle G(0,0)=d^{-1}.$ (44)

This result is also consistent with its lattice analog (20). Replacing $D\to
d/a$ in the lattice version of the global DOS (25), we immediately obtain its
continuum counterpart

$\displaystyle\rho(\omega)=\rho^{L}(\omega)+\rho^{R}(\omega)-\frac{1}{\pi}\text{Im}\,\frac{\partial}{\partial\omega}\ln\det
d(\omega+i0^{+}).$ (45)

The global DOS $\rho^{m}(\omega)$ of the disjoint subsystem $m$ (here $m=L,R$)
is generally given by

$\displaystyle\rho^{m}(\omega)=-\frac{1}{\pi}\text{Im}\int_{x_{m}}^{x_{m+1}}{dx}\,\text{tr}\,G^{m}(x,x;\omega+i0^{+}).$
(46)

Since $d(\omega)$ is hermitian on the real frequency axis, its determinant
$\det d(\omega)$ is real-valued. The real-valued roots of the equation

$\displaystyle\det d(\omega)=0$ (47)

determine energies $E_{B}$ of bound states induced by the tunneling between
the two subsystems. Their contribution to the global DOS (45) naturally
appears in the form $\sum_{B}\delta(\omega-E_{B})$.

In the infinite-space model (that is, when both the left and the right
subsystems are semi-infinite), the bound states found from the equation (47)
are localized near the interface of the two subsystems, and their energies
$E_{B}$ reside in the bandgaps of the whole system. It is also remarkable that
in this case the matrix $d$ can be alternatively written as

$\displaystyle d=[G^{(0,L)}(0,0)]^{-1}-p_{0},$ (48)

and hence the bound state equation (47) acquires the form

$\displaystyle\det\left[1-G^{(0,L)}(0,0;\omega)p_{0}(\omega)\right]=0.$ (49)

Here

$\displaystyle p_{0}$ $\displaystyle=\mathcal{U}_{1}$ (50)
$\displaystyle+\frac{1}{8}\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}\left[\mathcal{A}_{R}G^{R}(x,x)\mathcal{A}_{R}-\mathcal{A}_{L}G_{0}^{L}(x,x)\mathcal{A}_{L}\right]$

is the (generally energy-dependent) term breaking the translational invariance
of the auxiliary infinite-space model using the parameters of the left
subsystem, i.e. described by $G^{(0,L)}(x,x^{\prime})$. In addition, we
employed the auxiliary Green’s function $G_{0}^{L}(x,x^{\prime})$ which
describes the model to the right from the hard-wall potential at $x=0$ but
uses the parameters of the left subsystem (see Appendix A.2 for its explicit
expression as well as for the adopted conventions regarding notations). To
establish (48) we used the identity

$\displaystyle[G^{(0,L)}(0)]^{-1}$ (51)
$\displaystyle=-\frac{1}{8}\mathcal{A}_{L}\left[\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)+\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G_{0}^{L}(x,x)\right]\mathcal{A}_{L},$

which naturally arises in the translationally invariant case
$\mathcal{U}_{1}=0$ and
$(\mathcal{A}_{L},\mathcal{B}_{L},\mathcal{C}_{L})=(\mathcal{A}_{R},\mathcal{B}_{R},\mathcal{C}_{R})$.

Reminding ourselves that the boundary Green’s functions $G^{L}(x,x^{\prime})$
and $G^{R}(x,x^{\prime})$, entering Eqs. (31)-(34), admit a simple
representation (see Appendix A.2) in terms of $G^{(0,L)}(x,x^{\prime})$ and
$G^{(0,R)}(x,x^{\prime})$, respectively, we assert that the above construction
completes our program of establishing $G(x,x^{\prime})$ in the special $M=1$
case.

To verify that $G(x,x^{\prime})$ given by (31) does provide the solution of
the equation (7) with the matching conditions (5) and (6), we first note that
both $G^{L}(x,x^{\prime})+G^{R}(x,x^{\prime})$ and $F(x)$ satisfy the
differential equation (7). Second, on the basis of the relations
$G^{L}(0,x^{\prime})=G^{R}(0,x^{\prime})=0$ and (40) we establish the
continuity of $G(x,x^{\prime})$ at $x=0$, i.e. the condition (5) is fulfilled.
Third, we express the condition (6) for $G(x,x^{\prime})$ in the form

$\displaystyle\frac{\mathcal{A}_{R}}{2}\left[G_{1}^{R}(0^{+},x^{\prime})+F^{\prime}(0^{+})d^{-1}\bar{F}(x^{\prime})\right]$
$\displaystyle-\frac{\mathcal{A}_{L}}{2}\left[G_{1}^{L}(0^{-},x^{\prime})+F^{\prime}(0^{-})d^{-1}\bar{F}(x^{\prime})\right]$
$\displaystyle=-\left[\mathcal{U}_{1}-i\frac{\mathcal{B}_{R}-\mathcal{B}_{L}}{2}\right]d^{-1}\bar{F}(x^{\prime}),$
(52)

and check whether it is fulfilled for $x^{\prime}\neq 0$ (i.e. for
$x^{\prime}>0^{+}$ and for $x^{\prime}<0^{-}$). Under this condition we
identify
$\bar{F}(x^{\prime})=\frac{1}{2}\mathcal{A}_{R}G_{1}^{R}(0^{+},x^{\prime})-\frac{1}{2}\mathcal{A}_{L}G_{1}^{L}(0^{-},x^{\prime})$,
and it remains to prove that

$\displaystyle
d=-\mathcal{U}_{1}-\frac{\mathcal{A}_{R}}{2}F^{\prime}(0^{+})+\frac{\mathcal{A}_{L}}{2}F^{\prime}(0^{-})+i\frac{\mathcal{B}_{R}-\mathcal{B}_{L}}{2}.$
(53)

The fulfillment of this condition is shown in Appendix C. Thus it is finally
justified that the derived $G(x,x^{\prime})$ satisfies the both matching
conditions.

It is remarkable that the formula (53) along with the expressions (172) and
(173) for $F^{\prime}(0^{+})$ and $F^{\prime}(0^{-})$, respectively, represent
the simplest way of determining the matrix $d$ in the single-barrier case:

$\displaystyle d$
$\displaystyle=-\mathcal{U}_{1}+\mathcal{L}_{R}-\mathcal{L}_{L},$ (54)
$\displaystyle\mathcal{L}_{R}$
$\displaystyle=\frac{\mathcal{A}_{R}}{2}G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1}+\frac{i}{2}\mathcal{B}_{R},$
(55) $\displaystyle\mathcal{L}_{L}$
$\displaystyle=\frac{\mathcal{A}_{L}}{2}G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}+\frac{i}{2}\mathcal{B}_{L}.$
(56)

In this setting, the relation of $d$ to the matching condition (6), expressing
the conservation of the current density, becomes especially transparent. The
importance of the objects $\mathcal{L}_{R/L}$ for the properties of the
boundary charge in the half-space models has been recently elucidated in Ref.
[Piasotski  _et al._ , 2022].

###  II.3 Multiple barriers

A direct generalization of the single-barrier result (31) to the case of
multiple barriers is given by the expression (see Appendix D for its
derivation)

$\displaystyle G(x,x^{\prime})$
$\displaystyle=\sum_{m=0}^{M}G^{m}(x,x^{\prime})$
$\displaystyle+\sum_{m,m^{\prime}=1}^{M}F_{m}(x)(d^{-1})_{m,m^{\prime}}\bar{F}_{m^{\prime}}(x^{\prime}),$
(57)

where $d$ is a block tridiagonal matrix, with the additional barrier indices
$m,m^{\prime}$ labeling the blocks. The diagonal blocks

$\displaystyle d_{m,m}$
$\displaystyle=-\mathcal{U}_{m}-\frac{1}{8}\mathcal{A}_{m-1}\lim_{x\to
x_{m}^{-}}\frac{d^{2}}{dx^{2}}G^{m-1}(x,x)\mathcal{A}_{m-1}$
$\displaystyle-\frac{1}{8}\mathcal{A}_{m}\lim_{x\to
x_{m}^{+}}\frac{d^{2}}{dx^{2}}G^{m}(x,x)\mathcal{A}_{m}$ (58)

represent a generalization of (34): They are defined locally at the interface
positions $x_{m}$. In turn, the off-diagonal blocks

$\displaystyle d_{m,m+1}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{m}G_{12}^{m}(x_{m}^{+},x_{m+1}^{-})\mathcal{A}_{m},$
(59) $\displaystyle d_{m+1,m}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{m}G_{12}^{m}(x_{m+1}^{-},x_{m}^{+})\mathcal{A}_{m}$
(60)

describe the propagation between the two adjacent barriers $x_{m}$ and
$x_{m+1}$ across the $m$-th subsystem. These blocks thereby account for the
quantum interference effects.

The functions

$\displaystyle F_{m}(x)$
$\displaystyle=-\Theta(x_{m-1}<x<x_{m})\left[G_{2}^{(0,m-1)}(x,x_{m}^{+})+G_{2}^{m-1}(x,x_{m}^{-})-G_{2}^{(0,m-1)}(x,x_{m}^{-})\right]\frac{\mathcal{A}_{m-1}}{2}$
$\displaystyle+\Theta(x_{m}<x<x_{m+1})\left[G_{2}^{(0,m)}(x,x_{m}^{-})+G_{2}^{m}(x,x_{m}^{+})-G_{2}^{(0,m)}(x,x_{m}^{+})\right]\frac{\mathcal{A}_{m}}{2},$
(61) $\displaystyle\bar{F}_{m^{\prime}}(x^{\prime})$
$\displaystyle=-\Theta(x_{m^{\prime}-1}<x^{\prime}<x_{m^{\prime}})\frac{\mathcal{A}_{m^{\prime}-1}}{2}\left[G_{1}^{(0,m^{\prime}-1)}(x_{m^{\prime}}^{+},x^{\prime})+G_{1}^{m^{\prime}-1}(x_{m^{\prime}}^{-},x^{\prime})-G_{1}^{(0,m^{\prime}-1)}(x_{m^{\prime}}^{-},x^{\prime})\right]$
$\displaystyle+\Theta(x_{m^{\prime}}<x^{\prime}<x_{m^{\prime}+1})\frac{\mathcal{A}_{m^{\prime}}}{2}\left[G_{1}^{(0,m^{\prime})}(x_{m^{\prime}}^{-},x^{\prime})+G_{1}^{m^{\prime}}(x_{m^{\prime}}^{+},x^{\prime})-G_{1}^{(0,m^{\prime})}(x_{m^{\prime}}^{+},x^{\prime})\right]$
(62)

generalize the expressions (32) and (33), respectively. They also possess the
properties

$\displaystyle F_{m}(x_{m^{\prime}})=-\delta_{m,m^{\prime}},$ (63)
$\displaystyle\bar{F}_{m^{\prime}}(x_{m})=-\delta_{m^{\prime},m},$ (64)

analogous to (19) in the single-barrier case. On their basis we establish that

$\displaystyle G(x_{m},x_{m^{\prime}})=(d^{-1})_{m,m^{\prime}},$ (65)

and thus reveal the physical meaning of the matrix $d$: It is the inverse of
the propagator between the contacts $x_{m}$ and $x_{m^{\prime}}$.

The global single-barrier DOS (45) has a straightforward multi-barrier
generalization (see Appendix D)

$\displaystyle\rho(\omega)=$
$\displaystyle\sum_{m=0}^{M}\rho^{m}(\omega)-\frac{1}{\pi}\text{Im}\frac{\partial}{\partial\omega}\text{ln}\det
d(\omega+i0^{+}),$ (66)

where the determinant of $d(\omega+i0^{+})$ is now evaluated for the
$N_{c}M\times N_{c}M$ matrix.

The bound state equation (47) retains its form in the multi-barrier case.

In the double-barrier case $M=2$ we can further simplify the matrix $d$. We
assign the values $x_{1}=0$ and $x_{2}=W$ to the contact coordinates, such
that $W$ is the width of the central ($C$) region. We also assume that both
the left ($L$) and the right ($R$) subsystems are semi-infinite, and in the
following we use the labeling of the regions $m=L,C,R$ instead of $m=0,1,2$.
Owing to the expressions (217), (218), (219), (222) derived in Appendix E we
state

$\displaystyle d$
$\displaystyle\equiv\left(\begin{array}[]{cc}G(0,0)&G(0,W)\\\
G(W,0)&G(W,W)\end{array}\right)^{-1}$ (69)
$\displaystyle=\left(\begin{array}[]{cc}G^{(C,0)}(0,0)&G^{(C,0)}(0,W)\\\
G^{(C,0)}(W,0)&G^{(C,0)}(W,W)\end{array}\right)^{-1}-\left(\begin{array}[]{cc}p_{0}&0\\\
0&p_{W}\end{array}\right).$ (74)

The last term can be interpreted as the self-energy term breaking the
translational invariance of the auxiliary infinite-space model, which uses the
parameters of the central subsystem. It is expressed in terms of the energy-
dependent matrices

$\displaystyle p_{0}=\mathcal{U}_{1}$ (75)
$\displaystyle+\frac{1}{8}\lim_{x\to
0^{-}}\left[\mathcal{A}_{L}\frac{d^{2}}{dx^{2}}G^{L}(x,x)\mathcal{A}_{L}-\mathcal{A}_{C}\frac{d^{2}}{dx^{2}}G^{C}_{0}(x,x)\mathcal{A}_{C}\right],$
$\displaystyle p_{W}=\mathcal{U}_{2}$ (76)
$\displaystyle+\frac{1}{8}\lim_{x\to
W^{+}}\left[\mathcal{A}_{R}\frac{d^{2}}{dx^{2}}G^{R}(x,x)\mathcal{A}_{R}-\mathcal{A}_{C}\frac{d^{2}}{dx^{2}}G^{C}_{W}(x,x)\mathcal{A}_{C}\right].$

Hereby we introduced the auxiliary Green’s functions $G^{C}_{0}(x,x^{\prime})$
and $G^{C}_{W}(x,x^{\prime})$ describing the models to the left from the hard-
wall potential at $x=0$ and to the right from the hard-wall potential at
$x=W$, respectively, and both using the parameters of the central subsystem
(see Appendix A.2 for their explicit expressions).

Applying the formulas (166) and (167), we further simplify (75) and (76) to

$\displaystyle p_{0}$
$\displaystyle=\mathcal{U}_{1}+\mathcal{L}_{L}-\mathcal{L}_{L\to C},$ (77)
$\displaystyle p_{W}$
$\displaystyle=\mathcal{U}_{2}-\mathcal{L}_{R}+\mathcal{L}_{R\to C},$ (78)

where $\mathcal{L}_{R}$ and $\mathcal{L}_{L}$ are defined in (55) and (56),
and

$\displaystyle\mathcal{L}_{L\to C}$
$\displaystyle=\frac{\mathcal{A}_{C}}{2}G_{1}^{(0,C)}(0^{-},0)[G^{(0,C)}(0,0)]^{-1}+\frac{i}{2}\mathcal{B}_{C},$
(79) $\displaystyle\mathcal{L}_{R\to C}$
$\displaystyle=\frac{\mathcal{A}_{C}}{2}G_{1}^{(0,C)}(0^{+},0)[G^{(0,C)}(0,0)]^{-1}+\frac{i}{2}\mathcal{B}_{C}$
(80) $\displaystyle=\mathcal{L}_{L\to C}+[G^{(0,C)}(0,0)]^{-1}.$

These formulas allow one to express the matrix $d$ in the double-barrier case
given by (74) in terms of the objects $\mathcal{L}$ containing only first
derivatives of the translationally invariant Green’s functions $G^{(0,m)}$.

##  III Applications to Josephson systems

To showcase our formalism, we find it instructive to consider a number of
standard-issue problems in the theory of Josephson junctions (JJ).

By a JJ, one commonly understands a weak link between a pair of
superconductors. When two BCS condensates (labeled by 1 and 2) are brought
together and the tunneling of Cooper pairs between them is then switched on,
the combined system finds a new ground state, in which the difference
$\varphi=\varphi_{1}-\varphi_{2}$ of the phases $\varphi_{1},\ \varphi_{2}$ of
the corresponding superconducting order parameters
$\Delta_{1,2}=|\Delta_{1,2}|e^{i\varphi_{1,2}}$ adjusts itself to a particular
value $\varphi=\varphi_{\text{min}}$. In particular, the most common cases are
$\varphi_{\text{min}}=0$ and $\varphi_{\text{min}}=\pi$, defining the so-
called $0$\- and $\pi$-junctions. This result implies that the ground state
energy of the combined system is a function of the aforementioned phase
difference, with a minimum at $\varphi=\varphi_{\text{min}}$. In this regard,
it is important to study the spectral flow of Josephson systems with the
externally varied phase difference222In real experiments, the phase difference
is typically varied by closing the system into a circular geometry far away
from the tunneling region, and varying the magnetic flux threading the
resulting ring., and this defines the first type of problems for showcasing
our formalism. Specifically, in one of the following examples, we shall see
how local Coulomb interaction at the contact between the condensates may lead
to a crossover between the above-mentioned $\varphi_{\text{min}}=0$ and
$\varphi_{\text{min}}=\pi$ ground states.

As Josephson systems feature, by construction, superconducting components,
local charge conservation is violated resulting in the non-zero persistent
current (known as the Josephson current (JC)) between the condensates
comprising the junction. It turns outAbrikosov (2017) that such a current is
also a $\varphi$-dependent quantity that, quite generally, may be shown to be
the $\varphi$-derivative of the aforesaid ground state energy. The study of
the experimentally measurable JC defines the second problem, which we consider
for demonstrating the potential of Green’s function formalism.

###  III.1 Basic definitions

####  III.1.1 Model Hamiltonian

Before proceeding with concrete examples, we first specify notations of the
model Hamiltonians.

In what follows we restrict our consideration to spin-$\frac{1}{2}$ $s$-wave
superconductors, although – as is apparent from Section II – our formalism
allows including arbitrary matrix structure and momentum dependence of the
order parameters up to $O(p^{2})$.

We consider the following second quantized Hamiltonian

$\displaystyle\mathcal{H}=$
$\displaystyle\int_{-\infty}^{\infty}{dx}\hat{\psi}^{\dagger}(x)\left[p\frac{1}{2m(x)}p+\frac{1}{2}\\{A(x),\
p\\}+V(x)\right]\hat{\psi}(x)$
$\displaystyle+\frac{1}{2}\int_{-\infty}^{\infty}{dx}\left[\hat{\psi}^{\dagger}(x)\hat{\Delta}(x)(\hat{\psi}^{\dagger}(x))^{T}+\text{h.c.}\right],$
(81)

where
$\hat{\psi}(x)=(\hat{\psi}_{\uparrow}(x),\hat{\psi}_{\downarrow}(x))^{T}$ is a
two-component spinor, whose spin components $\sigma=\uparrow,\downarrow$ are
the field operators obeying the standard fermionic anticommutation relations

$\displaystyle\\{\hat{\psi}_{\sigma}(x),\
\hat{\psi}_{\sigma^{\prime}}(x^{\prime})\\}=$ $\displaystyle 0,$ (82)
$\displaystyle\\{\hat{\psi}_{\sigma}(x),\
\hat{\psi}_{\sigma^{\prime}}^{\dagger}(x^{\prime})\\}=$
$\displaystyle\delta_{\sigma,\sigma^{\prime}}\delta(x-x^{\prime}).$ (83)

The effective mass $m(x)$ is a piece-wise constant scalar function of $x$;
$A(x)=A^{\dagger}(x),V(x)=V^{\dagger}(x)$ are piece-wise constant $2\times 2$
hermitian matrices, and
$\hat{\Delta}(x)=-\hat{\Delta}^{T}(x)=\Delta(x)i\sigma_{y}$ is the
antisymmetric $2\times 2$ $s$-wave paring matrix expressed in terms of the
spatially dependent scalar order parameter $\Delta(x)$. We assign the
following spatial dependence to these objects:

$\displaystyle\begin{Bmatrix}m(x)\\\ A(x)\\\ V(x)\\\ \Delta(x)\end{Bmatrix}=$
$\displaystyle\Theta(-x)\begin{Bmatrix}m_{L}\\\ A_{L}\\\ V_{L}\\\
\Delta_{L}e^{i\varphi_{L}}\end{Bmatrix}$ (84) $\displaystyle+$
$\displaystyle\Theta(W>x>0)\begin{Bmatrix}m_{C}\\\ A_{C}\\\ V_{C}\\\
0\end{Bmatrix}+\Theta(x-W)\begin{Bmatrix}m_{R}\\\ A_{R}\\\ V_{R}\\\
\Delta_{R}e^{i\varphi_{R}}\end{Bmatrix}.$

Note that the left ($-\infty<x<0$) and the right ($+\infty>x>W$) regions host
superconductors, while the central – tunneling – region $W>x>0$ is normal
($|\Delta_{C}|=0$). Introducing the labelling of the regions in terms of the
index $\lambda=L,C,R$ (note that in the previous section it corresponds to the
index $m$, but here and below we use $\lambda$ to avoid confusion with the
mass notation), we further specify that $m_{\lambda}$ and $\Delta_{L,R}$ take
positive real values, while the hermitian matrices
$A_{\lambda}=A_{\lambda}^{(0)}+\sum_{j=x,y,z}A_{\lambda}^{(j)}\sigma_{j}\equiv
A_{\lambda}^{(0)}+\vec{A}_{\lambda}\cdot\vec{\sigma}$ and
$V_{\lambda}=V_{\lambda}^{(0)}+\sum_{j=x,y,z}V_{\lambda}^{(j)}\sigma_{j}\equiv
V_{\lambda}^{(0)}+\vec{V}_{\lambda}\cdot\vec{\sigma}$ are spanned by the Pauli
matrices $\sigma_{j}$ and the identity matrix with the real-valued
coefficients $A_{\lambda}^{(0,j)}$ and $V_{\lambda}^{(0,j)}$.

Introducing the extended Nambu spinor

$\displaystyle\hat{\Psi}(x)=\begin{pmatrix}\hat{\psi}(x)\\\
i\sigma_{y}(\hat{\psi}^{\dagger}(x))^{T}\end{pmatrix},$ (85)

we rewrite the Hamiltonian (81) in the form

$\displaystyle\mathcal{H}=\frac{1}{2}\int_{-\infty}^{\infty}{dx}\,\hat{\Psi}^{\dagger}(x)H\hat{\Psi}(x),$
(86)

where

$\displaystyle H=$ $\displaystyle\begin{pmatrix}h^{(0)}&\Delta(x)\\\
\Delta^{*}(x)&-\sigma_{y}h^{(0)*}\sigma_{y}\end{pmatrix},$ (87) $\displaystyle
h^{(0)}=$ $\displaystyle p\frac{1}{2m(x)}p+\frac{1}{2}\\{A(x),\ p\\}+V(x).$
(88)

As per common practice, we find it convenient to define a new set of Pauli
matrices $\tau_{x},\ \tau_{y},\ \tau_{z}$, along with an identity $\tau_{0}$,
acting on the space of particles (upper two components of $\Psi(x)$) and holes
(lower two components of $\Psi(x)$). The Hamiltonian (87) now may be written
as

$\displaystyle H=$
$\displaystyle\tau_{z}\left[p\frac{1}{2m(x)}p+\frac{1}{2}\\{\vec{A}(x)\cdot\vec{\sigma},\
p\\}+V^{(0)}(x)\right]$ (89)
$\displaystyle+\tau_{0}\left[\frac{1}{2}\\{A^{(0)}(x),\
p\\}+\vec{V}(x)\cdot\vec{\sigma}\right]$ (90)
$\displaystyle+\Delta(x)\tau_{+}+\Delta^{*}(x)\tau_{-}.$ (91)

In this decomposition of $H$, the terms (89), that is the kinetic energy, the
spin-orbit interaction, and the scalar potential, as well as the pairing
potential (91) are even under the standard time-reversal operation
$\hat{T}=i\sigma_{y}K$, with $K$ denoting the complex conjugation. In turn,
the terms collected in (90), that is the vector potential and the Zeeman
field, are odd under the time-reversal operation.

####  III.1.2 Observables

It is well-known (see e.g. in Ref. [Chamon  _et al._ , 2010]) that introducing
the extended Nambu representation artificially doubles the Hilbert space
assigned to the quantum system. This redundancy of the description has to be
removed in the calculation of observable quantities by enforcing a pseudo-
reality constraint on the Nambu field operators. Eventually this results in
removing the hole-like part of the spectrum residing at negative energies
$\omega<0$. It follows that the excitation spectrum of the system may be
directly inferred from the Eqs. (66), (74), (77)-(80) by restricting the
spectral range to $\omega>0$. For the JJ models, we choose the parameters

$\displaystyle\mathcal{A}_{\lambda}=\frac{\tau_{z}}{m_{\lambda}},\quad\mathcal{U}_{1,2}=0.$
(92)

Note that the contact potential strengths $\mathcal{U}_{1,2}$ are neglected,
since the physical effect of the tunneling barrier in the finite-$W$ setup is
accounted for by an appropriate tuning of $V_{C}$.

In the short junction limit, that is when the junction width $W$ is much
smaller than all other physical length scales in the system, we impose the
scaling $V_{C}=\frac{\mathcal{U}_{1}}{W}$ on the potential of the central
region, while all other terms in $H_{C}$ are supposed to be of $O(W^{0})$.
Then $H_{C}=[\mathcal{U}_{1}+O(W)]\,\delta_{W}(x)$, with the nascent
$\delta$-function $\delta_{W}(x)=\frac{1}{W}\Theta(W>x>0)$. In the limit $W\to
0$, we get the delta-distribution $\delta(x)=\lim_{W\to 0}\delta_{W}(x)$ and
neglect the $O(W)$ terms accompanying $\mathcal{U}_{1}$. The resulting model
has the single barrier at $x=0$ with the contact potential
$\mathcal{U}_{1}\delta(x)$, and the corresponding formulas for the spectral
density (45), (54)-(56) become applicable.

In the expressions (45) and (66) for the spectral density $\rho(\omega)$ of
the composite system there are terms $\rho^{m}(\omega)$ expressing the
spectral density of the isolated subsystems and the term expressed via $d$.
This last term represents the correction due to the tunneling between the
subsystems. It is the only term containing the dependence on the phase
difference $\varphi$, which we indicate explicitly:

$\displaystyle\delta\rho(\omega,\varphi)=-\frac{1}{\pi}\text{Im}\frac{\partial}{\partial\omega}\ln\det\,d(\omega+i0^{+},\varphi).$
(93)

The Josephson current may be expressed as the derivative of the Gibbs free
energy $F$ with respect to the phase difference across the superconducting
leads, $J(\varphi)=\frac{2e}{\hbar}\frac{dF}{d\varphi}$. It is then given in
terms of (93) by

$\displaystyle
J(\varphi)=-\frac{2e}{\hbar\beta}\int_{0}^{\infty}d\omega\ln\left(2\cosh\frac{\beta\omega}{2}\right)\frac{\partial}{\partial\varphi}\delta\rho(\omega,\varphi),$
(94)

where $\beta=\frac{1}{k_{B}T}$ is the inverse temperature. Integrating by
parts we reveal an alternative representation

$\displaystyle
J(\varphi)=-\frac{1}{\Phi_{0}}\text{Im}\int_{0}^{\infty}d\omega\,\tanh\frac{\beta\omega}{2}\frac{\partial}{\partial\varphi}\ln\det
d(\omega+i0^{+},\varphi),$ (95)

where $\Phi_{0}=\frac{h}{2e}$ is the superconducting magnetic flux quantum.

With the help of (95), the identity

$\displaystyle\frac{\partial}{\partial\varphi}\ln\det
d(z,\varphi)=\text{tr}\left\\{\left[d(z,\varphi)\right]^{-1}\frac{\partial
d(z,\varphi)}{\partial\varphi}\right\\}$ (96)

and the granted particle-hole symmetry, one can also convert the integral over
the real frequency axis into the (quickly convergent) Matsubara sum

$\displaystyle J(\varphi)=$
$\displaystyle-\frac{\pi}{\Phi_{0}\beta}\sum_{i\omega_{n}}\text{tr}\left\\{\left[d(i\omega_{n},\varphi)\right]^{-1}\frac{\partial
d(i\omega_{n},\varphi)}{\partial\varphi}\right\\}$ (97)
$\displaystyle\stackrel{{\scriptstyle T=0}}{{\rightarrow}}$
$\displaystyle-\frac{1}{2\Phi_{0}}\int_{-\infty}^{\infty}d\omega\,\text{tr}\left\\{\left[d(i\omega,\varphi)\right]^{-1}\frac{\partial
d(i\omega,\varphi)}{\partial\varphi}\right\\},$ (98)

which in the zero-temperature limit goes over into the integral (98) over the
imaginary frequency axis. These imaginary frequency representations appear
very efficient in numerical calculations of $J(\varphi)$.

###  III.2 Single barrier examples

First, we study the short junction case, in which the tunneling region is
described by the contact potential $\propto\delta(x)$, and the results of
Section II.2 are then employed.

####  III.2.1 Benchmark example: The Josephson model

Let us start by considering the paradigmatic problem of a pair of $s$-wave
superconductors tunnel-coupled with each other across a barrier at $x=0$. In
this case, the extended Nambu representation (85) is not needed, and the
matrix Hamiltonians are expressed in the reduced Nambu basis
$\hat{\Psi}(x)=(\hat{\psi}_{\uparrow}(x),\hat{\psi}_{\downarrow}^{\dagger}(x))^{T}$
as

$\displaystyle
H_{\lambda}=\begin{pmatrix}\frac{p^{2}}{2m}-\mu&\Delta_{\lambda}\\\
\Delta_{\lambda}^{*}&-\frac{p^{2}}{2m}+\mu\end{pmatrix},\quad\lambda=R,L,$
(99)

where the order parameter is assumed to be homogeneous in magnitude throughout
the sample $|\Delta_{\lambda}|=\Delta_{0}>0$. Choosing the symmetric gauge, we
express the phases of $\Delta_{\lambda}$ via the externally induced phase
difference $\varphi=\varphi_{R}-\varphi_{L}$ such that

$\displaystyle\Delta_{\lambda}=\Delta_{0}e^{i\varphi_{\lambda}},\quad\varphi_{\lambda}=\lambda\frac{\varphi}{2}.$
(100)

Hereby we conveniently re-labeled right $\lambda=+$ and left $\lambda=-$
subsystems. As is traditionally done, the nature of the barrier is modeled by
the contact potential

$\displaystyle\mathcal{U}(x)=\delta(x)V_{0}\tau_{z}.$ (101)

Figure 2: The phase-dispersion of the ABS (top panel) and zero-temperature JC
(bottom panel). The results for a clean system are shown in solid lines and
are compared to their disordered counterparts plotted in dashed lines. Various
colors correspond to different values of the chemical potential, as indicated
in the legends. To achieve the Andreev limit universality, we fix the values
(see above the top panel) of the barrier transparency $D$ as well as of the
energy $n_{I}v_{0}$ and momentum $mv_{0}$ scales associated with the impurity
scattering.

The $d$-matrix for this model is calculated in Appendix F.1. For
$\text{Im}\,z>0$ it reads

$\displaystyle
d(z,\varphi)=-V_{0}\tau_{z}+i\frac{k^{(+)}}{m}\tau_{z}+i\frac{k^{(-)}}{m}\frac{z{\color[rgb]{0,0.88,0}-}\Delta_{0}\cos\frac{\varphi}{2}\tau_{x}}{z\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}},$
(102)

where

$\displaystyle
k^{(\pm)}=k_{F}\frac{\sqrt{1+\frac{z}{\mu}\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}\mp\sqrt{1-\frac{z}{\mu}\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}}{2},$
(103)

and $k_{F}=\sqrt{2m\mu}$. The general bound state equation (47) yields

$\displaystyle 0=\frac{m^{2}}{k_{F}^{2}}\det d(\omega+i0^{+},\varphi)$ (104)
$\displaystyle=\left(\frac{k^{(+)}}{ik_{F}}+\frac{mV_{0}}{k_{F}}\right)^{2}+\frac{k^{(-)2}}{k_{F}^{2}}\left(1-\frac{\Delta_{0}^{2}}{\Delta_{0}^{2}-\omega^{2}}\sin^{2}\frac{\varphi}{2}\right),$

where $D=1/[1+(mV_{0}/k_{F})^{2}]$ is the transparency of the tunneling
barrier. It can be satisfied in the gap $|\omega|<\Delta_{0}$, where it holds

$\displaystyle
z\sqrt{1-\left(\frac{\Delta_{0}}{z}\right)^{2}}\bigg{|}_{z=\omega+i0^{+}}=i\sqrt{\Delta_{0}^{2}-\omega^{2}},$
(105)

and the corresponding replacements are to be made in $k^{(\pm)}$. It follows
that equation (104) defines the exact energy-phase relation, valid for all
values of $\mu$ and $D$.

Evaluating the JC on the basis of (95) [additionally multiplying it by the
factor 2 to account the two copies of (102) needed to reproduce the extended
Nambu representation] we note that the sub-gap contribution appears as the
residuum

$\displaystyle J_{\text{ABS}}(\varphi)$
$\displaystyle=\frac{2\pi}{\Phi_{0}}\tanh\frac{\beta\omega_{A}}{2}\frac{\partial_{\varphi}\det
d(\omega_{A},\varphi)}{\partial_{\omega}\det d(\omega_{A},\varphi)}$ (106)
$\displaystyle=-\frac{2\pi}{\Phi_{0}}\tanh\frac{\beta\omega_{A}(\varphi)}{2}\frac{d\omega_{A}(\varphi)}{d\varphi}$
(107)

at the pole $\omega=\omega_{A}$, given by the ABS energy found from the
equation (104). In turn, the continuum contribution to (95) equals

$\displaystyle J_{\text{cont}}(\varphi)=$
$\displaystyle-\frac{\sin\varphi}{\Phi_{0}}\int_{\Delta_{0}}^{\infty}d\omega\,\tanh\frac{\beta\omega}{2}\frac{\Delta_{0}^{2}}{\omega^{2}-\Delta_{0}^{2}}$
$\displaystyle\times\text{Im}\frac{k^{(-)2}}{m^{2}\,\det
d(\omega+i0^{+},\varphi)}.$ (108)

To obtain the known ABS expression

$\displaystyle\omega_{A}(\varphi)=\Delta_{0}\sqrt{1-D\sin^{2}\frac{\varphi}{2}},$
(109)

we invoke the so-called Andreev approximation relying on $\mu\gg\Delta_{0}$.
In this limit we find that $k^{(+)}\to 0$ and $k^{(-)}\to k_{F}$, simplifying

$\displaystyle d(z,\varphi)\approx-
V_{0}\tau_{z}+\frac{ik_{F}}{m}\frac{z-\Delta_{0}\cos\frac{\varphi}{2}\tau_{x}}{z\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}},$
(110)

as well as the equation (104) to the form which admits the solution (109). In
addition, we notice that the JC is completely mediated by the ABS, since
$k^{(-)}$ and $\det d(\omega+i0^{+})$ in (108) become purely real, and
therefore $J_{\text{cont}}(\varphi)$ vanishes in the Andreev limit. Thus, Eq.
(107) yields the following universal relationHaberkorn  _et al._ (1978);
Beenakker (1991)

$\displaystyle J(\varphi)=$
$\displaystyle\frac{\pi\Delta_{0}}{2\Phi_{0}}\frac{D\tanh\left[\frac{\beta\Delta_{0}\sqrt{1-D\sin^{2}\frac{\varphi}{2}}}{2}\right]}{\sqrt{1-D\sin^{2}\frac{\varphi}{2}}}\sin\varphi.$
(111)

In Fig. 2, in solid lines, we show how the energy-phase and the corresponding
current-phase relations approach the universal results of Eqs. (109), (111)
upon an increase in the chemical potential $\mu$.

It is also worth mentioning that the expression (110) multiplied by $-1$ is
analogous by virtue of (44) to the self-energy of a quantum dot coupled to
superconducting leads in the weak tunneling regime which is further
approximated in the wide-band limit, see e.g. Ref. [Žonda  _et al._ , 2015]
for that self-energy expression. One has to replace the tunneling rate
$\Gamma$ of that model with $\frac{k_{F}}{m}$ to gain the formal analogy with
(110). In our treatment, however, an effective $\Gamma$ is itself of the
bandwidth’s order of magnitude, that is $\sim\mu$ [cf. the general discussion
in Section II.2].

Figure 3: Panels a) and e): Dispersion of sub-gap states as a function of
$\varphi$ and $u$ for a system with $\mu=3\Delta_{0},\ D=0.9$. Panels b) and
f): Josephson current as a function of $\varphi$ and $u$ for
$\mu=3\Delta_{0},\ D=0.9$. Panels c) and g): derivative of Josephson current
with respect to phase difference $\varphi$ at $\varphi=0$ for a system with
$\mu=3\Delta_{0}$. The orange line in panel c) indicates the sign change of
$\frac{d}{d\varphi}J(\varphi=0)$. Panels d) and h): The critical current as a
function of $u$ and $D$ for $\mu=3\Delta_{0}$. The orange line in panel d)
indicates the location of the cusp in the critical current [explicitly seen in
panel h)] as a function of Coulomb interaction strength $u$.

One may further ask whether the JC value beyond the Andreev limit in our model
is robust against static disorder, for example. Our approach allows one to get
analytical insights into such questions at a modest expense. In particular, we
may ignore the effects of disorder on the barrier tunneling, assuming that its
imperfection is completely accounted for by the contact potential (101). Under
such an assumption, we put the self-energy insertions into the bulk
propagators $G^{(0,R/L)}$ alone.

Let us consider random non-magnetic impurities, characterized by the strength
$v_{0}$ of the short-range impurity potential and the impurity density
$n_{I}$. Employing the standard $T$-matrix approximation [Altland and Simons,
2010], we obtain the following disorder-induced self-energy

$\displaystyle\Sigma^{\lambda}(z)=$ $\displaystyle
U_{\lambda}\Sigma(z)U_{\lambda}^{\dagger},\quad
U_{\lambda}=e^{\frac{i}{4}\tau_{z}\lambda\varphi},$ (112)
$\displaystyle\Sigma(z)=$ $\displaystyle
n_{I}v_{0}\tau_{z}\frac{1}{1-G^{(0)}(0,0;z)v_{0}\tau_{z}}+O(n_{I}^{2}),$ (113)

dressing the momentum-space propagators of the bulk superconductors:

$\displaystyle G^{(0,\lambda)}_{k}(z)$
$\displaystyle=U_{\lambda}G^{(0)}_{k}(z)U_{\lambda}^{\dagger}\rightarrow\tilde{G}^{(0,\lambda)}_{k}(z)=U_{\lambda}\tilde{G}^{(0)}_{k}(z)U_{\lambda}^{\dagger},$
$\displaystyle\tilde{G}^{(0)}_{k}(z)$
$\displaystyle=\frac{1}{[G^{(0)}_{k}(z)]^{-1}-\Sigma(z)},$ (114)

where $G^{(0)}_{k}(z)$ is given by (225). Representing

$\displaystyle\tilde{G}^{(0)}_{k}(z)=\frac{1}{\tilde{z}(z)-\tau_{z}[\frac{k^{2}}{2m}-\tilde{\mu}(z)]-\tau_{x}\tilde{\Delta}_{0}(z)}$
(115)

and comparing it with (114), we find that the effect of the random disorder
consists in the following energy-dependent renormalization of the model
parameters

$\displaystyle z\rightarrow$
$\displaystyle\tilde{z}(z)=z\left(1-n_{I}v_{0}k_{1,+}k_{2,+}\frac{imv_{0}k^{(-)}}{(mv_{0}k^{(-)})^{2}+(k_{1,+}k_{2,+}+imv_{0}k^{(+)})^{2}}\frac{1}{z\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}\right),$
(116) $\displaystyle\Delta_{0}\rightarrow$
$\displaystyle\tilde{\Delta}_{0}(z)=\Delta_{0}\left(1-n_{I}v_{0}k_{1,+}k_{2,+}\frac{imv_{0}k^{(-)}}{(mv_{0}k^{(-)})^{2}+(k_{1,+}k_{2,+}+imv_{0}k^{(+)})^{2}}\frac{1}{z\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}\right),$
(117) $\displaystyle\mu\rightarrow$ $\displaystyle\tilde{\mu}(z)=\mu-
n_{I}v_{0}k_{1,+}k_{2,+}\frac{k_{1,+}k_{2,+}+imv_{0}k^{(+)}}{(mv_{0}k^{(-)})^{2}+(k_{1,+}k_{2,+}+imv_{0}k^{(+)})^{2}},$
(118)

where $k_{1,+}$ and $k_{2,+}$ are given in (227) and (228).

Since the $d$-function is eventually expressed via the parameters of the bulk
systems, it is sufficient to make the above replacements in (102) in order to
extract spectral information of the composite system subject to random
potential disorder. Using the invariance

$\displaystyle\frac{z}{\Delta_{0}}=\frac{\tilde{z}(z)}{\tilde{\Delta}_{0}(z)},$
(119)

one finds that the dispersion of Andreev levels in Eq. (104) is modified by
the following replacement

$\displaystyle\frac{k^{(\pm)}}{k_{F}}\to\frac{\tilde{k}^{(\pm)}}{k_{F}}=\sqrt{\frac{\tilde{\mu}(z)}{\mu}}$
(120)
$\displaystyle\times\frac{\sqrt{1+\frac{\tilde{z}(z)}{\tilde{\mu}(z)}\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}\mp\sqrt{1-\frac{\tilde{z}(z)}{\tilde{\mu}(z)}\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}}{2}.$

Returning back to the Andreev approximation
$\tilde{\mu}(z)\approx\mu\gg\Delta_{0}$, we find that the results (109) and
(111) remain intact.

In Fig. 2, in dashed lines, we show how the energy-phase and the corresponding
current-phase relations approach the universal results of Eqs. (109), (111)
upon an increase in the chemical potential $\mu$. In addition, Fig. 2
showcases the effect of disorder on the energy- and current-phase relations
away from the universal Andreev limit, showing that static disorder tends to
push the Andreev states into the continuum decreasing the critical current,
similar to the effect of decreasing the barrier’s transparency.

Next, we discuss the effect of the local Coulomb interaction at the point
contact of the two superconductors. We model it as the Coulomb repulsion
between fluctuations of local spin-up and spin-down densities,

$\displaystyle V_{U}=$ $\displaystyle
U\left\\{\hat{\psi}^{\dagger}_{\uparrow}(0)\hat{\psi}_{\uparrow}(0)-\langle
n_{\uparrow}(0)\rangle\right\\}\left\\{\hat{\psi}^{\dagger}_{\downarrow}(0)\hat{\psi}_{\downarrow}(0)-\langle
n_{\downarrow}(0)\rangle\right\\}.$ (121)

With this term our model resembles the Anderson-Josephson model of a quantum
dot with the on-site Coulomb interaction tunnel-coupled to two superconducting
leads (see the review about various studies of this model in Ref. [Meden,
2019]). However, in our case, there is no well-defined quantum dot, since we
treat the contact far beyond the tunneling regime.

We also note that the subtraction of the density averages in Eq. (121)
eventually leads to the energy renormalization in the Cooper channel alone
(see Eq. (122) below, derived under the assumption of no spontaneous spin-
symmetry breaking). Performing this subtraction, we get rid of an
uninteresting renormalization of the diagonal contact-potential component in
the particle-hole basis, which is largely dominated by high-energy normal-
state contributions. One can alternatively envisage this subtraction as a
result of combining the particle-hole diagonal contribution to the
(restricted) Hartree-Fock self-energy with the bare contact potential, which
leads to an effective contact potential and defines the contact’s physical
transparency renormalized by the local Coulomb interaction.

The Hartree-Fock approximation to the local Coulomb interaction (121) is found
on the basis of the Wick’s theorem

$\displaystyle V_{U}$ $\displaystyle\simeq
U\langle\hat{\psi}_{\downarrow}(0)\hat{\psi}_{\uparrow}(0)\rangle\hat{\psi}_{\uparrow}^{\dagger}(0)\hat{\psi}_{\downarrow}^{\dagger}(0)$
$\displaystyle+U\langle\hat{\psi}_{\uparrow}^{\dagger}(0)\hat{\psi}_{\downarrow}^{\dagger}(0)\rangle\hat{\psi}_{\downarrow}(0)\hat{\psi}_{\uparrow}(0)$
$\displaystyle-U\langle\hat{\psi}_{\uparrow}^{\dagger}(0)\hat{\psi}_{\downarrow}^{\dagger}(0)\rangle\langle\hat{\psi}_{\downarrow}(0)\hat{\psi}_{\uparrow}(0)\rangle$
$\displaystyle=\Delta_{loc}\\{\hat{\psi}_{\uparrow}^{\dagger}(0)\hat{\psi}_{\downarrow}^{\dagger}(0)+\hat{\psi}_{\downarrow}(0)\hat{\psi}_{\uparrow}(0)\\}-\frac{\Delta_{loc}^{2}}{U},$
(122)

where
$\Delta_{loc}=U\langle\hat{\psi}_{\downarrow}(0)\hat{\psi}_{\uparrow}(0)\rangle=U\langle\hat{\psi}_{\uparrow}^{\dagger}(0)\hat{\psi}_{\downarrow}^{\dagger}(0)\rangle$
is a local superconducting order parameter.

The local Green’s function $G_{U}(0,0)$ in the Hartree-Fock approximation
results from the Dyson equation

$\displaystyle d_{U}$ $\displaystyle=[G_{U}(0,0)]^{-1}$ (123)
$\displaystyle=[G(0,0)]^{-1}-\Delta_{loc}\tau_{x}=d-\Delta_{loc}\tau_{x}.$
(124)

Using it in the expression (97) [times factor 2 because of the present usage
of the reduced Nambu basis] for the JC and accounting for the correction
$-\frac{\Delta_{loc}^{2}}{U}$ to the free energy occurring in (122), we
establish the mean-field expression for the JC

$\displaystyle J(\varphi)=$
$\displaystyle-\frac{2\pi}{\Phi_{0}\beta}\sum_{i\omega_{n}}\text{tr}\left\\{\left[d_{U}(i\omega_{n},\varphi)\right]^{-1}\frac{\partial
d_{U}(i\omega_{n},\varphi)}{\partial\varphi}\right\\}$
$\displaystyle-\frac{2\pi}{\Phi_{0}}\frac{d}{d\varphi}\frac{\Delta_{loc}^{2}}{U}.$
(125)

With the help of the self-consistency equation (see in Appendix F.2)

$\displaystyle\Delta_{loc}=\frac{U}{2\beta}\sum_{i\omega_{n}}\text{tr}\left\\{[d_{U}(i\omega_{n},\varphi)]^{-1}\tau_{x}\right\\}$
(126)

the expression (127) is modified to the form

$\displaystyle J(\varphi)=$
$\displaystyle-\frac{2\pi}{\Phi_{0}\beta}\sum_{i\omega_{n}}\text{tr}\left\\{\left[d_{U}(i\omega_{n},\varphi)\right]^{-1}\frac{\partial
d(i\omega_{n},\varphi)}{\partial\varphi}\right\\}.$ (127)

A further application of the self-consistent Hartree-Fock approach to the
present model reaches its limitation, emerging in the form of an ultraviolet
divergence which is caused by the ultra-local form of the Coulomb interaction
(see discussion in Appendix F.2), and one has to refine the method.

A similar problem arises in the bare perturbation theory in the Andreev limit
$\mu\gg\Delta_{0}$: The local superconducting order parameter $\Delta_{loc}$
features the logarithmic behavior
$\propto\frac{k_{F}}{m}\ln\frac{\mu}{\Delta_{0}}$, which hints at the
nontrivial competition between the superconductivity and the Kondo
effect.Meden (2019)

To give a conservative quantitative estimate of the interaction effects on the
JC and the Andreev spectrum, we provide below the results of the first-order
perturbation theory for $\mu\gtrsim\Delta_{0}$ and the dimensionless
interaction parameter $u=\frac{Um^{2}\Delta_{0}}{\pi k_{F}^{2}}\ll 1$. In
particular, we approximate the local self-energy term
$\Delta_{loc}(\varphi)=\frac{k_{F}}{m}\bar{\Delta}_{loc}(\varphi)$ by its
leading $O(u)$ contribution given in (251). The resulting $d_{U}$ is used for
the JC evaluation by means of (127) [or (253) at zero temperature] as well as
in the bound state equation

$\displaystyle\det d_{U}(\omega,\varphi)=0.$ (128)

Numerical data for the ABS and JC for a junction with $\mu=3\Delta_{0}$ are
shown in Fig. 3. The panels a) and e) demonstrate the energy-phase relation of
the sub-gap states for various values of the interaction parameter $u\in[0,\
0.5]$ and contact transparency $D=0.9$, while the panels b) and f) show the
phase-dispersion of the JC at these parameters. We observe that upon
increasing $u$ above $\sim 0.2$ the JC derivative [with respect to $\varphi$]
changes its sign from positive to negative, indicating the $0-\pi$ phase
transition. It is very analogous to the phase transition which is well-
established in superconductor-ferromagnet-superconductor junctions [Robinson
_et al._ , 2007, Chtchelkatchev  _et al._ , 2001]. The panels c) and g) show
the JC derivative at $\varphi=0$ as a function of $u$ and $D$. In particular,
we observe that the phase transition occurs at smaller values of $u$ in more
transparent junctions with $D\to 1$. The panels d) and h) demonstrate the
critical current $J_{c}=\max_{\varphi}|J(\varphi)|$ as a function of $u$ and
$D$. We reveal that around the phase transition, the critical current displays
a non-differentiable cusp in the dependence on the interaction parameter $u$.

####  III.2.2 Junction of two Majorana wires: spectral properties

Let us now apply our formalism to the famous Majorana wire problem [Lutchyn
_et al._ , 2010,Oreg  _et al._ , 2010]. In particular, we consider a pair of
two semiconducting wires with the strong spin-orbit interaction $\alpha$ and
induced superconducting correlations, which are additionally submersed into
the external magnetic field $B$ pointing in the wires’ direction. In the
extended Nambu basis, introduced in Section III.1.1, they are described by the
following Bogoliubov-de-Gennes Hamiltonians

$\displaystyle H_{\lambda}=$
$\displaystyle\left(\begin{array}[]{cc}h_{p}^{(0)}+\alpha
p\sigma_{z}-B\sigma_{x}&\Delta_{\lambda}\\\
\Delta_{\lambda}^{*}&-h_{p}^{(0)}-\alpha
p\sigma_{z}-B\sigma_{x}\end{array}\right),$ (131) $\displaystyle
h_{p}^{(0)}=\frac{p^{2}}{2m}-\mu;\quad\lambda=R,L\equiv+,-.$ (132)

Assuming the wires to be semi-infinite, we bring them in contactvan Heck  _et
al._ (2017); Murthy  _et al._ (2020) at $x=0$ enabling the tunneling across
the contact potential

$\displaystyle\mathcal{U}(x)=\delta(x)\,V_{0}\,\tau_{z}\,\sigma_{0},$ (133)

where $\sigma_{0}$ is the identity matrix in the spin space.

As in the previous consideration, we consider again for simplicity the case of
the isospectral junction with $|\Delta_{\lambda}|=\Delta_{0}>0$, with a
symmetrically induced phase difference
$\Delta_{\lambda}=\Delta_{0}e^{\frac{i}{2}\lambda\varphi}$ across the
interface.

The bulk Green’s functions

$\displaystyle G^{(0,\lambda)}(x,x^{\prime})$
$\displaystyle=U_{\lambda}\,g(x-x^{\prime})\,U_{\lambda}^{\dagger},\quad
U_{\lambda}=e^{\frac{i}{4}\tau_{z}\lambda\varphi},$ (134)

associated with the Hamiltonians (131), are evaluated in Appendix G.1. This
information appears sufficient for establishing the composite Green’s function
$G(x,x^{\prime})$ (31) and the $d$-matrix (54) containing spectral properties
relevant for the study of the ABS and JC. For the present model we obtain

$\displaystyle d=$ $\displaystyle\frac{\tau_{z}}{2m}\left\\{-2mV_{0}\right.$
(135)
$\displaystyle\left.\qquad+U_{+}g^{\prime}(0^{+})g^{-1}(0)U_{-}-U_{-}g^{\prime}(0^{-})g^{-1}(0)U_{+}\right\\}.$

Figure 4: The ABS dispersion in the short junction of the two Majorana wires
at zero chemical potential $\mu=0$ and large spin-orbit energy
$E_{\text{SO}}=\frac{m\alpha^{2}}{2}=60\,\Delta_{0}$ dominating over
$\Delta_{0}$ and $B$. Left and right panels refer to the non-topological
($B=0.75\,\Delta_{0}<\Delta_{0}$) and topological
($B=1.25\,\Delta_{0}>\Delta_{0}$) regimes, respectively. As one can see, the
numerical result based on the exact equation $\det d=0$ perfectly agrees with
the approximate low-energy result derived in Appendix G.2 [see Eqs. (282) and
(283)], which coincides with that of Ref. [Murthy  _et al._ , 2020].

The calculation of the JC by means of Eqs. (95)-(98) requires the additional
derivative $\partial_{\varphi}d$. Since the $\varphi$-dependence of $d$ in the
above formula enters only via the gauging matrices $U_{\pm}$, we note the
following useful relation

$\displaystyle\partial_{\varphi}d=$
$\displaystyle\frac{i\tau_{z}}{8m}\left\\{U_{+}[\tau_{z},g^{\prime}(0^{+})g^{-1}(0)]U_{-}\right.$
$\displaystyle\left.\qquad+U_{-}[\tau_{z},g^{\prime}(0^{-})g^{-1}(0)]U_{+}\right\\}.$
(136)

On the basis of the exact expressions for $g(0)$ and $g^{\prime}(0^{\pm})$
(see in Appendix G.1) one can address various limits of the model’s
parameters. For example, we show in Appendix G.2 how the recent results of
Ref. [Murthy  _et al._ , 2020] for the energy spectrum and the JC in the
regime of dominating spin-orbit energy
$E_{\text{SO}}\equiv\frac{m\alpha^{2}}{2}\gg\Delta_{0},B,\mu$ can be
analytically recovered from our general expression for the $d$-function. We
also use it to numerically evaluate the ABS dispersion at large
$E_{\text{SO}}=60\,\Delta_{0}$ [see in Fig. 4] to benchmark the present exact
interface Green’s function approach versus the scattering approach of Ref.
[Murthy  _et al._ , 2020] relying on low-energy approximations in the spin-
orbit dominated regime.

The demonstration presented in Appendix G.2 serves a more general purpose of
explaining how one can derive low-energy approximations for $G(x,x^{\prime})$
in arbitrary heterostructures. Since $G(x,x^{\prime})$ can be always expressed
according to our present findings in terms of bulk Green’s functions, it
suffices to make a low-energy approximation for these functions. This
approximation typically relies on the bulk spectrum linearization near Fermi
points (see also Ref. [Piasotski  _et al._ , 2022] for a similar discussion),
and the approximate bulk Green’s functions are evaluated much easier than
their exact counterparts.

At the same time, it is no longer needed to investigate how the low-energy
approximation affects the matching condition (6) involving derivatives of the
wavefunctions. This is usually a subtle problem since the spectrum
linearization lowers by one the order of the Schrödinger differential
equation, and this requires relaxing the condition on the first derivatives.
On the other hand, the matching condition (6) expresses the current density
conservation, which must be somehow accounted for in the construction of the
eigenfunctions. Our approach circumvents this problem, since it allows us to
make the low-energy approximation directly for $G(x,x^{\prime})$, skipping any
intermediate approximate treatment of Eq. (6).

The generality of our approach allows one to go beyond the low-energy
approximation relying on the dominance of the spin-orbit interaction energy
$E_{\text{SO}}$ in the present model. Thus it enables exploring arbitrary
parameter regimes of the Majorana junction model, as demonstrated in the
following.

Figure 5: The ABS dispersion in the short Majorana junction at zero chemical
potential $\mu=0$ and moderate spin-orbit interaction strength
$E_{\text{SO}}=\frac{m\alpha^{2}}{2}=4\Delta_{0}$. Left and right panels refer
to the non-topological ($B<\Delta_{0}$) and topological ($B>\Delta_{0}$)
regimes, respectively. Blue, orange, and green colors are used to mark three
different strengths of the contact potential $4m\xi V_{0}=0$, $4m\xi V_{0}=1$,
and $4m\xi V_{0}=5$, respectively, which are expressed in terms of the
Ginzburg–Landau coherence length $\xi=1/\sqrt{2m\Delta_{0}}$.

In particular, in Fig. 5 we display the dispersion of the sub-gap states for
the moderate spin-orbit energy $E_{\text{SO}}=4\,\Delta_{0}$, at zero chemical
potential $\mu=0$. We note that in the non-topological regime $B<\Delta_{0}$,
addressed in the left panel of Fig. 5, Andreev levels are being pushed into
the continuum upon an increase of the contact potential strength $V_{0}$, akin
to the behavior in ordinary JJs. In the topological regime, shown in the right
panel of Fig. 5, we find that the increase of $V_{0}$ results in the
unwrapping of the Andreev mode, such that its tails near $\varphi=0,\ 2\pi$
are pushed towards zero energy.

The results for a yet different parameter regime with
$E_{\text{SO}}=4\Delta_{0}$ and – more essentially – nonzero
$\mu=\sqrt{3}\Delta_{0}$ are shown in Fig. 6. As before, the left and right
panels refer to the non-topological and topological regimes, respectively. The
nonzero $\mu$ replaces the boundary value $B$ of the topological phase
transition from $\Delta_{0}$ to the larger value
$\sqrt{\Delta_{0}^{2}+\mu^{2}}$ ($=2\Delta_{0}$ in our specific example).
Apart form that, we observe the apparent increase in the value of $V_{0}$
required to push the Andreev states into and away from the continuum of the
scattering states in the non-topological and topological phases, respectively.

Figure 6: The ABS dispersion in the short Majorana junction at finite chemical
potential $\mu=\sqrt{3}\Delta_{0}$ and moderate spin-orbit interaction
strength $E_{\text{SO}}=\frac{m\alpha^{2}}{2}=4\Delta_{0}$. Left and right
panels refer to the non-topological ($B<\sqrt{\Delta_{0}^{2}+\mu^{2}}$) and
topological ($B>\sqrt{\Delta_{0}^{2}+\mu^{2}}$) regimes, respectively. Blue,
orange, and green colors are used to mark three different strengths of the
contact potential $4m\xi V_{0}=0$, $4m\xi V_{0}=1$, and $4m\xi V_{0}=5$,
respectively.

####  III.2.3 Junction of two Majorana wires: Josephson current

An example of the zero-temperature JC computed for the model parameters
$E_{\text{SO}}=4\Delta_{0}$, $\mu=\sqrt{3}\Delta_{0}$ is shown in Fig. 7, with
the left and right panels referring to the non-topological and topological
phases, respectively. In the non-topological phase, the JC demonstrates a
sinusoidal behavior with a tilt towards the high-symmetry point $\varphi=\pi$.
On the contrary, we see that the JC exhibits a sharp step-like discontinuity
at the high-symmetry point $\varphi=\pi$ in the topological phase (see the
left panel of Fig. 7). This property relates to the discontinuity in the
derivative of the ground state energy (see Fig. 6) with respect to the phase
difference
$\partial_{\varphi}E_{\text{GS}}(\varphi)\propto\text{sgn}(\varphi-\pi),\
\varphi\approx\pi$. Such a result is nonphysical and arises from the non-
commutativity of the zero-energy and zero-temperature limits, as is most
easily seen in the representation (95):

$\displaystyle\lim_{T\rightarrow 0^{+}}\lim_{\omega\rightarrow
0^{+}}\tanh\left[\frac{\omega}{2k_{B}T}\right]=0,$ (137)
$\displaystyle\lim_{\omega\rightarrow 0^{+}}\lim_{T\rightarrow
0^{+}}\tanh\left[\frac{\omega}{2k_{B}T}\right]=1.$ (138)

This implies that any arbitrarily small but finite temperature will smear the
sharp step (see Fig. 8). Additionally, Fig. 7 demonstrates the effect of the
barrier imperfection $V_{0}\neq 0$ on the JC, quite conventionally implying
the reduction of the critical current $J_{c}=\max_{\varphi}|J(\varphi)|$ with
the increase in the back-scattering strength $V_{0}$.

Figure 7: The zero-temperature ($T=0$) limit of the JC calculated at finite
chemical potential $\mu=\sqrt{3}\Delta_{0}$ and relatively small spin-orbit
energy $E_{\text{SO}}=\frac{m\alpha^{2}}{2}=4\Delta_{0}$. Left and right
panels refer to the non-topological
($B=1.25\Delta_{0}<\sqrt{\Delta_{0}^{2}+\mu^{2}}=2\Delta_{0}$) and topological
phases ($B=2.50\Delta_{0}>2\Delta_{0}$), respectively.

The finite-temperature effects are shown in Fig. 8, with the left and right
panels referring to the non-topological and topological phases, as before. In
the non-topological phase, the key effects of thermodynamic fluctuations are
the simultaneous decrease in the critical current and demolition of the
aforementioned $\pi$-tilt of the current-phase relation. As is mentioned in
the previous paragraph, in the topological phase the non-zero temperature
destroys the sharp step in the current profile, further affecting the current-
phase relation in a manner typical for the non-topological phase.

###  III.3 Double barrier example: Long junction of two Majorana wires

We continue by studying a long Josephson junction between the two Majorana
wires. Now we assume that our system consists of the three layers: the right
($R$) and the left ($L$) superconducting leads, and the central ($C$) normal
region which is free of superconducting correlations, as is outlined in
Section III.1.1. The Hamiltonians of the right and left superconducting leads
are given by Eq. (131), while the Hamiltonian of the central region is

$\displaystyle H_{C}=$
$\displaystyle\tau_{z}h_{p}^{(0,C)}+\tau_{z}\alpha_{C}p\sigma_{z}-B_{C}\sigma_{x},$
(139)

with $h^{(0,C)}_{p}$ defined as in Eq. (132) with $\mu\rightarrow\mu_{C}$ and
$m_{C}=m$.

The bound state equation (47), using (74) in the two-barrier case, can be
conveniently written as

$\displaystyle\det\left(\begin{array}[]{cc}1-\hat{p}_{0}\,g^{C}(0)&-\hat{p}_{0}\,g^{C}(-W)e^{\frac{i}{2}\tau_{z}\varphi}\\\
-\hat{p}_{W}\,g^{C}(W)e^{-\frac{i}{2}\tau_{z}\varphi}&1-\hat{p}_{W}\,g^{C}(0)\end{array}\right)=0,$
(142)

where the translationally invariant Green’s function $g^{C}(x-x^{\prime})$ for
the normal central region is evaluated in Appendix G.3. We also shifted the
phase dependence from $p_{0}=U_{-}\,\hat{p}_{0}\,U_{+}$ and
$p_{W}=U_{+}\,\hat{p}_{W}\,U_{-}$ to the off-diagonal elements, responsible
for the quantum coherence in the normal section, in order to elucidate its
importance for supporting the JC. In the present model, the non-decaying terms
$g^{C}(\pm W)$ are achieved due to the gapless spectrum of the Hamiltonian
(139) (in particular, its outer branches do not gap out due to
$\Delta_{0}=0$). The matrices

$\displaystyle\hat{p}_{0,W}$
$\displaystyle=\pm\frac{\tau_{z}}{2m}g^{\prime}(0^{\mp})[g(0)]^{-1}\mp\frac{\tau_{z}}{2m}g^{C\,^{\prime}}(0^{\mp})[g^{C}(0)]^{-1}$
$\displaystyle\pm\frac{i}{2}(\alpha-\alpha_{C})\tau_{z}\sigma_{z}$ (143)

do not carry the $\varphi$-dependence and expressed with the help of the
Green’s function $g(x-x^{\prime})$ of the bulk superconductors evaluated in
Appendix G.1.

In Fig. 9 we show the phase-dependence of the junction-localised bound states
and of the JC at zero and finite temperatures. We focus on the case of the
isospectral junction with $B_{C}=B,\ \mu_{C}=\mu=\sqrt{3}\Delta_{0},\
\alpha_{C}=\alpha=4\Delta_{0}\xi$ varying the junction’s width $W=0.5\,\xi,\
2.5\xi,\ 5\xi,\ 10\xi$ [in the units of the coherence length
$\xi=1/\sqrt{2m\Delta_{0}}$] at the two values of the magnetic field
$B=1.25\Delta_{0},\ 2.5\Delta_{0}$.

Figure 8: The effect of the finite temperature on the JC, calculated at finite
chemical potential $\mu=\sqrt{3}\Delta_{0}$ and relatively small spin-orbit
energy $E_{\text{SO}}=\frac{m\alpha^{2}}{2}=4\Delta_{0}$. Left and right
panels refer to the non-topological
($B=1.25\Delta_{0}<\sqrt{\Delta_{0}^{2}+\mu^{2}}=2\Delta_{0}$) and topological
phases ($B=2.50\Delta_{0}>2\Delta_{0}$), respectively. Figure 9: The phase
dependence of the ABS and JC for long JJs with
$\alpha_{C}=\alpha=4\Delta_{0}\xi,\ \mu_{C}=\mu=\sqrt{3}\Delta_{0}$, and
$B_{C}=B$, for a variety of the junction’s widths $W=0.5\xi,\ 2.5\xi,\ 5\xi$,
and $W=10\xi$ (four different columns). The two upper rows correspond to the
non-topological regime
$B=1.25\Delta_{0}<\sqrt{\Delta^{2}_{0}+\mu^{2}}=2\Delta_{0}$, while the two
lower rows correspond to the topological phase $B=2.5\Delta_{0}>2\Delta_{0}$.
In the second and fourth rows, different colors are used to mark different
temperatures $T=0,\ 0.05\Delta_{0},\ 0.2\Delta_{0}$, and $0.5\Delta_{0}$.

Specifically, the first and third rows demonstrate the evolution of the energy
dispersion with an increase in the width of the central region. We see that in
the non-topological phase
($B=1.25\,\Delta_{0}<\sqrt{\Delta_{0}^{2}+\mu^{2}}=2\Delta_{0}$) the addition
of the normal segment between the two Majorana wires has the effect of
lowering the energy of the sub-gap states, as well as increasing their number.
In turn, in the topological phase ($B=2.5\,\Delta_{0}>2\Delta_{0}$) the only
bound state mode tends to spread over the whole phase interval featuring the
piece-wise linear branches of its phase dispersion.

As for the JC shown in the second and fourth rows of the same figure, we find
that it features qualitatively similar behavior to its short-junction
counterpart, with the major effect of the finite width being the decrease in
the critical current.

Figure 10: The phase dependence of the ABS and JC for long JJs for the long
($W=2.5\xi$) JJ with $\alpha_{C}=\alpha=4\Delta_{0}\xi$. The two upper rows
correspond to the non-topological phase $B=1.25\Delta_{0}$, while the two
lower rows describe the topological regime $B=2.5\Delta_{0}$. The five
different columns correspond to five different values $B_{C}=0,\ 0.5B,\ B,\
1.5B$, and $B_{C}=2B$ of the Zeeman field in the central region. Three
distinct colors mark three different values $\mu_{C}=0,\ \mu,\ 2\mu$ of the
chemical potential in the central region.

Next, we consider an example of the model’s realization, in which some of the
parameters of the central region are distinct from the corresponding ones in
the superconducting leads. For instance, let us study the effect of modifying
the parameters of Abelian ($\mu_{C}$) and non-Abelian ($B_{C}$) parts of the
scalar potential. Experimentally this may be achieved by applying an
appropriate gate voltage (to affect $\mu_{C}$) and by bringing the wire in
proximity with a ferromagnet (to affect $B_{C}$).

In Fig. 10, we present the results for the ABS and the zero temperature JC in
the model with $\alpha_{C}=\alpha=4\xi\Delta_{0},\ W=2.5\xi$, and
$\mu=\sqrt{3}\Delta_{0}$. The five different columns correspond to five
different values of $B_{C}=0,\ 0.5B,\ B,\ 1.5B,$ and $2B$. Like in Fig. 9, the
first and third rows of Fig. 10 show the phase-dispersion of the ABS in the
non-topological ($B=1.25\Delta_{0}$) and topological ($B=2.5\Delta_{0}$)
regimes of the Majorana wires, respectively, while the second and fourth rows
give the corresponding zero-temperature current-phase relations. Three
distinct colors are used to indicate three different values of the chemical
potential in the normal region: $\mu_{C}=0$ (red), $\mu_{C}=\mu$ (black), and
$\mu_{C}=2\mu$ (yellow).

In the non-topological phase ($B=1.25\Delta_{0}$, first and second rows), we
observe that, by varying the parameters of the central region, it is possible
to push the ABS energies downwards as much as needed for creating the crossing
points at zero energy. They are similar to the one observed in the topological
phase in the short junction. However, these crossing points are not
topological in nature, and arise from the crossing of two particle-like and
hole-like Andreev bands, that are symmetric around the $\varphi=\pi$ point,
and hence come in pairs 333Note that one can envisage the situation in which
the two bands touch one another at exactly zero energy. In this case, however,
the states around zero energy are expected to disperse quadratically, as
opposed to the linear dispersion found in the topological regime.. In turn,
the crossing at zero energy in the topological phase arises from the
intersection of two $\varphi=\pi$-asymmetric bands, belonging to two distinct
parity branches, and is topologically protected [Cheng and Lutchyn, 2012]. We
note that an analogous effect was also observed in superconductor-ferromagnet-
superconductor heterostructures in the presence of spin-orbit coupling in Ref.
[Yokoyama  _et al._ , 2014].

As for the zero-temperature JC in the non-topological phase, we reveal that
the appearance of the zero-energy-touching Andreev branches results in the
discontinuities of the current-phase relations, akin to the topological
regime. Like the touching points themselves, the step-like discontinuities
also come in pairs and are notably falling behind the critical current in
size. The latter effect has to do with the presence of additional strongly
dispersing Andreev states, which significantly contribute to the current.

In the topological phase ($B=2.5\Delta_{0}$, third and fourth rows), we see
that independently of the parameter values inside the central region, there is
always only one bound state touching the zero energy at $\varphi=\pi$. In
addition, for the cases of the equal ($\mu_{C}=\mu$, black curve) and enhanced
($\mu_{C}=2\mu_{C}$, yellow curve) chemical potentials of the central segment
we observe that an additional bound state appears for $B_{C}<B$. As in the
previous consideration of the topological regime, we find that the current-
phase relation again features sharp step-like discontinuity at the high-
symmetry point $\varphi=\pi$. We note that for most parameter values the size
of this discontinuity is inappreciably smaller than the critical current.

##  IV Conclusions and Outlook

In this paper, we provided a detailed account of the interface Green’s
function technique regarding a large class of quasi-one-dimensional models of
heterostructures. In our analysis, we assumed the Hamiltonians of the
individual layers to be $N_{c}\times N_{c}$ matrix-valued quadratic
polynomials in momentum operator conjugate to the composition axis of the
structure. Such Hamiltonians may be seen as arising from the expansion of
microscopic ones around the respective Fermi surfaces, allowing one to model a
broad range of realistic physical systems (multi-band systems, systems with
superconducting correlations, etc). For the considered class of models, we
established the representation of the position space Green’s function of the
entire system in terms of bulk position space Green’s functions of its sub-
parts. By demonstrating that the calculation of latter objects is uninvolved,
we opened up a simple pathway to study the multi-barrier scattering phenomena
and the interface-localized bound states in systems of interest. Furthermore,
our work lays one of the first stones into the analysis of many-body and
disorder effects in layered systems, requiring knowledge of Green’s functions
as input. By further extracting the global DOS of the system from its Green’s
function, we reveal that the spectral information may, largely, be drawn from
a single matrix $d$. For a system with $M$ tunneling barriers, we show that
$d$ is an $N_{c}M\times N_{c}M$ block matrix, comprised of $N_{c}\times
N_{c}$-sized blocks admitting for a simple representation in terms of bulk
Green’s functions of the layers. As the spectral density relates to the
logarithmic derivative of $\det\,d$, we also establish that the bound state
energies of the system lie at its zeroes.

We shall point out that our formalism explicitly deals with ballistic contact
problems (beyond the tunneling regime), commonly studied in modern-day
experiments [Matute-Cañadas  _et al._ , 2022]. As nowadays the veil of
theoretical secrecy over the ballistic regime is just being lifted, we hope
our method finds its extensive use in this regard.

The method is further exemplified on models of normal and topological
Josephson junctions. As a prototypical example, we considered a junction
between two conventional superconductors with an induced phase difference
across the interface. We derived the exact energy- and current-phase relations
for the considered model, showing how they reduce to the commonly known
formulas in the Andreev approximation, for which the chemical potential in
superconducting leads is assumed to be the largest energy scale in the model.
Further, we studied whether the universal Josephson relations, arising in the
Andreev limit, are robust against static disorder. By disorder-dressing the
bulk Green’s functions of superconductors in the $T$-matrix approximation and
using those to assemble the Green’s function of the composite system, we
demonstrated the stability of the Josephson relations in the Andreev
approximation, revealing only small corrections to the ABS energy beyond it.
Next, we studied the effect of the local Coulomb interaction at the contact
point between superconductors, showing how it leads to the $0-\pi$ phase
transition, detectable through the change in the current-phase relation of the
system.

As for the topological Josephson systems, we considered models of short and
long junctions between two Majorana wires. We demonstrated how to calculate
energy-phase and current-phase relations for the considered models beyond the
low-energy approximation, reducing the numerical computation to root-searching
problems.

In addition, we indicate that our method may serve as a convenient starting
point to develop low-energy approximations for considered models. When
projecting the Hamiltonians of different subsystems onto the low-energy window
of interest, one is typically challenged to identify and further implement the
correct matching conditions on the distinct low-energy degrees of freedom. Our
method allows us to elegantly circumvent this problem by developing the low-
energy approximations of bulk single-layer Green’s functions, as in Ref.
[Piasotski  _et al._ , 2022]. In this paper, we demonstrated how this strategy
is realized in the Majorana junction model (see Appendix G.2), reproducing the
low-energy results of Ref. [Murthy  _et al._ , 2020].

In our future work, we are going to generalize our developments to the class
of position-dependent single-layer Hamiltonians, which are important to
account for the effects of inhomogeneous external fields.

## Acknowledgments

KP thanks P. Ostrovsky for valuable comments. MP is grateful to P. Khomyakov,
V. Meden, H. Schoeller, and A. Svetogorov for useful discussions. KP and AS
acknowledge funding from DFG Projects SH 81/6-1 and SH 81/7-1.

##  Appendix A Evaluation of the Green’s functions

###  A.1 Green’s functions of translationally invariant systems

Let us consider the infinite-space model described by the Hamiltonian (1) and
evaluate its Green’s function (9) in the position representation. It can be
conveniently rewritten as

$\displaystyle
G^{(0,m)}(x,x^{\prime};z)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}\frac{e^{ik(x-x^{\prime})}}{Q(k;z)}P(k;z),$
(144)

where $P(k;z)=\text{adj}[z-h_{m}(k)]$ is the adjugate matrix, and
$Q(k;z)=\det[z-h_{m}(k)]$ is a $2N_{c}$-order polynomial of $k$, admitting the
representation

$\displaystyle
Q(k;z)=\det(\textstyle{\frac{1}{2}}\mathcal{A}_{m})\prod_{s=1}^{2N_{c}}(k-k_{s}(z))$
(145)

in terms of the roots $k_{s}(z)$.

By the hermiticity of $h_{m}(k)$ the determinant $Q(k;\omega)$ on the real
frequency axis may only have either a pair of complex conjugated roots or
purely real roots. Upon the shift $\omega\to\omega+i0^{+}$ the latter acquire
infinitesimal imaginary parts, whose signs are given by the signs of
$\frac{dk_{s}(\omega)}{d\omega}$, or equivalently by the signs of the
dispersion slope (group velocity) $\frac{d\omega_{k}}{dk}\big{|}_{k=k_{s}}$.
For the Hamiltonian (1), a number of dispersion branches pointing up (down) at
negative $k$ is the same as a number of dispersion branches pointing up (down)
at positive $k$. This implies that for every $\omega$ the number of real roots
with the positive dispersion slope equals the number of real roots with the
negative dispersion slope. Therefore we can generally classify all roots into
the two equally sized sets $\\{k_{s,+}(\omega+i0^{+})\\}_{s=1}^{N_{c}}$ and
$\\{k_{s,-}(\omega+i0^{+})\\}_{s=1}^{N_{c}}$ according to

$\displaystyle\text{Im}\left[k_{s,+}(\omega+i0^{+})\right]>0,\quad\text{Im}\left[k_{s,-}(\omega+i0^{+})\right]<0.$
(146)

Evaluating (144) at $z=\omega+i0^{+}$ we close the integration contour in the
upper and lower halves of the complex momentum plane for $x>x^{\prime}$ and
$x<x^{\prime}$, respectively, to correspondingly pick up the residues at the
isolated singularities $k_{s,+}$ or $k_{s,-}$. This results in

$\displaystyle G^{(0,m)}(x,x^{\prime};\omega+i0^{+})=$
$\displaystyle\sum_{s=1}^{N_{c}}\frac{i\Theta(x-x^{\prime})e^{ik_{s,+}|x-x^{\prime}|}}{\det(\frac{1}{2}\mathcal{A}_{m})\prod_{s^{\prime}\neq
s}(k_{s,+}-k_{s^{\prime},+})\prod_{s^{\prime}}(k_{s,+}-k_{s^{\prime},-})}P(k_{s,+};\omega)$
$\displaystyle-\sum_{s=1}^{N_{c}}\frac{i\Theta(x^{\prime}-x)e^{-ik_{s,-}|x-x^{\prime}|}}{\det(\frac{1}{2}\mathcal{A}_{m})\prod_{s^{\prime}}(k_{s,-}-k_{s^{\prime},+})\prod_{s^{\prime}\neq
s}(k_{s,-}-k_{s^{\prime},-})}P(k_{s,-};\omega).$ (147)

Thus the evaluation of (9) reduces to a problem of the polynomial
factorization (145) with a root post-selection criteria (146). In many
practical applications, these tasks may be executed analytically, otherwise
one can resort to numerical methods, such as the NumPy’s built-in routine
numpy.roots, for example [Harris  _et al._ , 2020].

###  A.2 Boundary Green’s functions

For completeness of the exposition, we shall quote some basic results on the
propagators of bounded systems. For a detailed account of the subject, we
refer an interested reader to Refs. [Wimmer, 2009, Piasotski, 2023], while for
various interesting applications to transport and solid state phenomena see
Refs. [Piasotski  _et al._ , 2022; Zazunov  _et al._ , 2018; Alvarado  _et
al._ , 2020; Kaladzhyan and Bena, 2019; Pinon  _et al._ , 2020; Komnik and
Heinze, 2017; Zazunov  _et al._ , 2017; Peng  _et al._ , 2017; Zazunov  _et
al._ , 2016; Arrachea  _et al._ , 2009; Müller  _et al._ , 2021].

Consider a model described by Green’s function $G^{(0,m)}(x,x^{\prime})$.
Inserting an ultra-local infinite-height potential at $x=X$, which creates a
hard-wall boundary, we separate our initial system into two disjoint
subsystems $x>X$ and $x<X$. The Green’s functions of the both of them are
given by the expression

$\displaystyle G^{m}_{X}(x,x^{\prime})=G^{(0,m)}(x,x^{\prime})$
$\displaystyle-G^{(0,m)}(x,X)\left[G^{(0,m)}(X,X)\right]^{-1}G^{(0,m)}(X,x^{\prime}),$
(148)

satisfying the boundary conditions $G^{m}_{X}(X,x^{\prime})=G^{m}_{X}(x,X)=0$.
It also holds that $G^{m}_{X}(x,x^{\prime})$ identically vanishes when either
$x>X,x^{\prime}<X$ or $x<X,x^{\prime}>X$.

The expression (148) also has the lattice analog (see e.g. in Ref. [Müller
_et al._ , 2021])

$\displaystyle
G_{n,n^{\prime}}^{m}=G^{(0,m)}_{n,n^{\prime}}-G^{(0,m)}_{n,0}[G_{0,0}^{(0,m)}]^{-1}G_{0,n^{\prime}}^{(0,m)}$
(149)

for a nearest-neighbor tight-binding model with the infinite-height potential
at the site $n=0$.

Adding the second hard-wall boundary at $x=Y$, we obtain the Green’s function

$\displaystyle G^{m}_{X,Y}(x,x^{\prime})=G^{m}_{X}(x,x^{\prime})$
$\displaystyle-G^{m}_{X}(x,Y)\left[G^{m}_{X}(Y,Y)\right]^{-1}G^{m}_{X}(Y,x^{\prime})$
(150) $\displaystyle=G^{m}_{Y}(x,x^{\prime})$
$\displaystyle-G^{m}_{Y}(x,X)\left[G^{m}_{Y}(X,X)\right]^{-1}G^{m}_{Y}(X,x^{\prime}),$
(151)

satisfying the boundary conditions
$G^{m}_{X,Y}(X,x^{\prime})=G^{m}_{X,Y}(x,X)=G^{m}_{X,Y}(Y,x^{\prime})=G^{m}_{X,Y}(x,Y)=0$.
It can be also represented in the symmetrized form

$\displaystyle G^{m}_{X,Y}(x,x^{\prime})=G^{(0,m)}(x,x^{\prime})$
$\displaystyle-G^{(0,m)}(x,X)\left[G^{m}_{Y}(X,X)\right]^{-1}G^{m}_{Y}(X,x^{\prime})$
$\displaystyle-G^{(0,m)}(x,Y)\left[G^{m}_{X}(Y,Y)\right]^{-1}G^{m}_{X}(Y,x^{\prime}).$
(152)

Assuming that $X<Y$, we observe that the Green’s function
$G^{m}_{X,Y}(x,x^{\prime})$ identically vanishes when $x$ and $x^{\prime}$
belong to different intervals (smaller than $X$, between $X$ and $Y$, larger
than $Y$). In addition, the following identities hold:

$\displaystyle G^{m}_{X,Y}(x,x^{\prime})=G^{m}_{X}(x,x^{\prime})\quad
x,x^{\prime}<X,$ (153) $\displaystyle
G^{m}_{X,Y}(x,x^{\prime})=G^{m}_{Y}(x,x^{\prime})\quad x,x^{\prime}>Y.$ (154)

With these definitions, we obtain the following expression for the boundary
Green’s functions introduced in the main text:

$\displaystyle
G^{m}(x,x^{\prime})=\Theta(x_{m+1}>x>x_{m})G^{m}_{x_{m},x_{m+1}}(x,x^{\prime}),$
(155)

where we assumed that neither $x_{0}=-\infty$ nor $x_{M+1}=\infty$. In these
special cases, we get simpler expressions on the basis of (148)

$\displaystyle G^{0}(x,x^{\prime})$
$\displaystyle=\Theta(x_{1}-x)G^{0}_{x_{1}}(x,x^{\prime}),\quad
x_{0}=-\infty,$ (156) $\displaystyle G^{M}(x,x^{\prime})$
$\displaystyle=\Theta(x-x_{M})G^{M}_{x_{M}}(x,x^{\prime}),\quad
x_{M+1}=\infty.$ (157)

##  Appendix B Derivation of the continuum limit in the single-barrier model

To derive the continuum limit $a\to 0$ of the lattice Green’s functions (31)
we expand $G^{L}_{-1,-1}$ and $G^{R}_{1,1}$ by two orders higher than the
leading one: This is necessary to cancel the leading $O(\frac{1}{a^{2}})$ term
in (16) and to correctly extract the subleading $O(\frac{1}{a})$ contribution.
This expansion is facilitated by the properties
$G^{L}_{-1,0}=G^{L}_{0,-1}=G^{L}_{0,0}=0$ and
$G^{R}_{1,0}=G^{R}_{0,1}=G^{R}_{0,0}=0$ which are inherent to the boundary
Green’s functions (149). Using these identities as well as
$G^{(0,L)}_{-1,-1}=G^{(0,L)}_{0,0}$ we have

$\displaystyle G^{L}_{-1,-1}$
$\displaystyle=G^{(0,L)}_{0,0}-G^{(0,L)}_{-1,0}-G^{(0,L)}_{0,-1}+G^{(0,L)}_{0,0}$
$\displaystyle-[G^{(0,L)}_{-1,0}-G^{(0,L)}_{0,0}][G_{0,0}^{(0,L)}]^{-1}[G_{0,-1}^{(0,L)}-G_{0,0}^{(0,L)}]$
$\displaystyle\approx
a^{2}G_{1}^{(0,L)}(0^{-},0)-\frac{a^{3}}{2}G_{11}^{(0,L)}(0^{-},0)$
$\displaystyle-a^{2}G_{1}^{(0,L)}(0^{+},0)-\frac{a^{3}}{2}G_{22}^{(0,L)}(0,0^{-})$
$\displaystyle-a^{3}G_{1}^{(0,L)}(0^{-},0)][G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{-}).$
(158)

Hereby we identify $G_{n,n^{\prime}}^{(0,L)}$ with $G^{(0,L)}(x,x^{\prime})$
via (30) and use the subindices 1 and 2 to indicate partial derivatives of the
latter function with respect to the corresponding arguments.

Observing that $t_{m}=\frac{\mathcal{A}_{m}-ia\mathcal{B}_{m}}{2a^{2}}$ and

$\displaystyle\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)=-G_{11}^{(0,L)}(0^{-},0)-G_{22}^{(0,L)}(0,0^{-})$
$\displaystyle-2G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{-}),$
(159)

and using the jump conditions (42),(43), we obtain

$\displaystyle
t_{L}G^{L}_{-1,-1}t_{L}^{\dagger}=\frac{\mathcal{A}_{L}-ia\mathcal{B}_{L}}{2a^{2}}G^{L}_{-1,-1}\frac{\mathcal{A}_{L}+ia\mathcal{B}_{L}}{2a^{2}}$
$\displaystyle=-\frac{\mathcal{A}_{L}}{2a^{2}}+\frac{1}{8a}\mathcal{A}_{L}\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)\mathcal{A}_{L}+O(a^{0}).$ (160)

Analogously we find that

$\displaystyle t_{R}^{\dagger}G^{R}_{1,1}t_{R}$
$\displaystyle=-\frac{\mathcal{A}_{R}}{2a^{2}}+\frac{1}{8a}\mathcal{A}_{R}\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{R}(x,x)\mathcal{A}_{R}+O(a^{0}).$ (161)

Collecting together the terms defining $D$ in (16) we observe the cancellation
of the $O(\frac{1}{a^{2}})$ terms thanks to the choice of the $a$-scaling in
(29). The remaining $O(\frac{1}{a})$ terms give rise to (34) defining the
matrix $d$. Note that the linear $z$-dependence is suppressed in this limit
since it appears in the next-to-subleading order of the $a$-expansion.

To approximate the remaining terms $G^{L}_{n,-1}$, $G^{R}_{n,1}$,
$G^{L}_{-1,n^{\prime}}$, and $G^{R}_{1,n^{\prime}}$ in (15), it suffices to
expand them up to the leading order which turns out to be $O(a^{2})$:

$\displaystyle
G^{L}_{n,-1}=G^{L}_{n,-1}-G^{L}_{n,0}\approx-a^{2}G_{2}^{(0,L)}(x,0^{+})$
$\displaystyle+a^{2}G^{(0,L)}(x)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{-})$
$\displaystyle=-a^{2}\left[G_{2}^{(0,L)}(x,0^{+})+G_{2}^{L}(x,0^{-})-G_{2}^{(L,0)}(x,0^{-})\right],$
(162) $\displaystyle G^{R}_{n,1}=G^{R}_{n,1}-G^{R}_{n,0}\approx
a^{2}G_{2}^{(0,R)}(x,0^{-})$
$\displaystyle-a^{2}G^{(0,R)}(x,0)[G^{(0,R)}(0,0)]^{-1}G_{2}^{(0,R)}(0,0^{+})$
$\displaystyle=a^{2}\left[G_{2}^{(0,R)}(x,0^{-})+G_{2}^{R}(x,0^{+})-G_{2}^{(0,R)}(x,0^{+})\right],$
(163) $\displaystyle
G^{L}_{-1,n^{\prime}}=G^{L}_{-1,n^{\prime}}-G^{L}_{0,n^{\prime}}\approx-a^{2}G_{1}^{(0,L)}(0^{+},x^{\prime})$
$\displaystyle+a^{2}G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}G^{(0,L)}(0,x^{\prime})$
$\displaystyle=-a^{2}\left[G_{1}^{(0,L)}(0^{+},x^{\prime})+G_{1}^{L}(0^{-},x^{\prime})-G_{1}^{(0,L)}(0^{-},x^{\prime})\right],$
(164) $\displaystyle
G^{R}_{1,n^{\prime}}=G^{R}_{1,n^{\prime}}-G^{R}_{0,n^{\prime}}\approx
a^{2}G_{1}^{(0,R)}(0^{-},x^{\prime})$
$\displaystyle-a^{2}G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1}G^{(0,R)}(0,x^{\prime})$
$\displaystyle=a^{2}\left[G_{1}^{(0,R)}(0^{-},x^{\prime})+G_{1}^{R}(0^{+},x^{\prime})-G_{1}^{(0,R)}(0^{+},x^{\prime})\right].$
(165)

Multiplying them with $t_{L}^{\dagger}$, $t_{R}$, $t_{L}$, $t_{R}^{\dagger}$,
respectively, cancels the factor $a^{2}$, and we eventually obtain (32) and
(33).

##  Appendix C Proof of the relation (53)

To prove (53) it suffices to show that

$\displaystyle\mathcal{A}_{R}F^{\prime}(0^{+})-iB_{R}=\frac{1}{4}\mathcal{A}_{R}\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{R}(x,x)\mathcal{A}_{R},$ (166)
$\displaystyle\mathcal{A}_{L}F^{\prime}(0^{-})-iB_{L}=-\frac{1}{4}\mathcal{A}_{L}\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)\mathcal{A}_{L}.$ (167)

First we evaluate

$\displaystyle F^{\prime}(0^{+})=[G_{12}^{(0,R)}(0^{+},0^{-})$ $\displaystyle-
G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1}G_{2}^{(0,R)}(0,0^{+})]\frac{\mathcal{A}_{R}}{2},$
(168) $\displaystyle F^{\prime}(0^{-})=-[G_{12}^{(0,L)}(0^{-},0^{+})$
$\displaystyle-
G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{-})]\frac{\mathcal{A}_{L}}{2}.$
(169)

Using the identities

$\displaystyle G_{12}^{(0,R)}(0^{+},0^{-})$
$\displaystyle=G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1}G_{2}^{(0,R)}(0,0^{-}),$
(170) $\displaystyle G_{12}^{(0,L)}(0^{-},0^{+})$
$\displaystyle=G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{+}),$
(171)

previously derived in Ref. [Piasotski  _et al._ , 2022] [see Eqs. (B2) and
(B3) therein], as well as the jump conditions (42) and (43), we simplify

$\displaystyle F^{\prime}(0^{+})$
$\displaystyle=-G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1},$ (172)
$\displaystyle F^{\prime}(0^{-})$
$\displaystyle=-G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}.$ (173)

Next, we evaluate

$\displaystyle-\frac{1}{8}\mathcal{A}_{R}\left[\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{R}(x,x)-\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G_{0}^{R}(x,x)\right]\mathcal{A}_{R}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{R}G_{1}^{(0,R)}(0^{+},0)[G^{(0,R)}(0,0)]^{-1}G_{2}^{(0,R)}(0,0^{+})\mathcal{A}_{R}$
$\displaystyle-\frac{1}{4}\mathcal{A}_{R}G_{1}^{(0,R)}(0^{-},0)[G^{(0,R)}(0,0)]^{-1}G_{2}^{(0,R)}(0,0^{-})\mathcal{A}_{R}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{R}[G_{12}^{(0,R)}(0^{+},0^{-})-G_{12}^{(0,R)}(0^{-},0^{+})]\mathcal{A}_{R}$
$\displaystyle-\mathcal{A}_{R}F^{\prime}(0^{+})-[G^{(0,R)}(0,0)]^{-1}$ (174)

and

$\displaystyle-\frac{1}{8}\mathcal{A}_{L}\left[\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)-\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G_{0}^{L}(x,x)\right]\mathcal{A}_{L}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{L}G_{1}^{(0,L)}(0^{-},0)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{-})\mathcal{A}_{L}$
$\displaystyle-\frac{1}{4}\mathcal{A}_{L}G_{1}^{(0,L)}(0^{+},0)[G^{(0,L)}(0,0)]^{-1}G_{2}^{(0,L)}(0,0^{+})\mathcal{A}_{L}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{L}[G_{12}^{(0,L)}(0^{-},0^{+})-G_{12}^{(0,L)}(0^{+},0^{-})]\mathcal{A}_{L}$
$\displaystyle+\mathcal{A}_{L}F^{\prime}(0^{-})-[G^{(0,L)}(0,0)]^{-1}.$ (175)

Further using (51) we establish that

$\displaystyle\mathcal{A}_{R}F^{\prime}(0^{+})-\frac{1}{4}\mathcal{A}_{R}\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{R}(x,x)\mathcal{A}_{R}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{R}[G_{12}^{(0,R)}(0^{+},0^{-})-G_{12}^{(0,R)}(0^{-},0^{+})]\mathcal{A}_{R}$
(176)

and

$\displaystyle\mathcal{A}_{L}F^{\prime}(0^{-})+\frac{1}{4}\mathcal{A}_{L}\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G^{L}(x,x)\mathcal{A}_{L}$
$\displaystyle=\frac{1}{4}\mathcal{A}_{L}[G_{12}^{(0,L)}(0^{+},0^{-})-G_{12}^{(0,L)}(0^{-},0^{+})]\mathcal{A}_{L}.$
(177)

The relations (166) and (167) follow from these by virtue of the jump
condition in the mixed derivatives

$\displaystyle\frac{1}{4}\mathcal{A}_{m}[G^{(0,m)}_{12}(0^{+},0^{-})-G^{(0,m)}_{12}(0^{-},0^{+})]\mathcal{A}_{m}=i\mathcal{B}_{m},$
(178)

which is a generalization of the equation (18) in Ref. [Piasotski  _et al._ ,
2022].

##  Appendix D Green’s function derivation for the multi-barrier model

The Green’s function computation in the single-barrier model presented in
Section II.2 can be straightforwardly generalized to the case of arbitrary $M$
barriers. Defining their positions at sites $n=n_{1},\ldots,n_{M}$ and
assigning the local onsite potentials $W_{n_{m}}$ to each of them, we also
introduce $M+1$ disjoint subsystems defined on the intervals $n_{m}+1\leq
n\leq n_{m+1}-1$ and describe them by the Green’s functions $G^{m}$.

Switching on the coupling $\sum_{m=1}^{M}v_{m}$ between the disjoint
subsystems via the hoppings through the barrier sites $n_{m}$,

$\displaystyle v_{m}=$ $\displaystyle-|n_{m}\rangle t_{n_{m}-1}\langle
n_{m}-1|-|n_{m}-1\rangle t_{n_{m}-1}^{\dagger}\langle n_{m}|$
$\displaystyle-|n_{m}+1\rangle t_{n_{m}}\langle n_{m}|-|n_{m}\rangle
t_{n_{m}}^{\dagger}\langle n_{m}+1|,$ (179)

we set up the Dyson equation treating $v_{m}$’s as a perturbation. It has the
form

$\displaystyle G_{n,n^{\prime}}=G^{\\{m\\}}_{n,n^{\prime}}$ (180)
$\displaystyle-\sum_{m^{\prime}=1}^{M}(G_{n,n_{m^{\prime}}-1}^{m^{\prime}-1}t_{n_{m^{\prime}}-1}^{\dagger}+G_{n,n_{m^{\prime}}+1}^{m^{\prime}}t_{n_{m^{\prime}}})G_{n_{m^{\prime}},n^{\prime}}$
$\displaystyle-\sum_{m^{\prime}=1}^{M}\delta_{n,n_{m^{\prime}}}g_{n_{m^{\prime}},n_{m^{\prime}}}^{m^{\prime}}(t_{n_{m^{\prime}}-1}G_{n_{m^{\prime}}-1,n^{\prime}}+t_{n_{m^{\prime}}}^{\dagger}G_{n_{m^{\prime}}+1,n^{\prime}})$

generalizing (14). Hereby $G^{\\{m\\}}_{n,n^{\prime}}$ contains not only all
Green’s functions $G^{m}$ of the finite-ranged disjoint subsystems, but also
those
$g_{n,n^{\prime}}^{m}=\frac{\delta_{n,n_{m}}\delta_{n^{\prime},n_{m}}}{z-W_{n_{m}}}$
of the barrier sites $n_{m}$.

Below we find the solution of (180) following the same routine as in Section
II.2.

Introducing the $M$-component vectors

$\displaystyle\bar{\mathcal{G}}_{m^{\prime}}^{(\pm)}$
$\displaystyle=G_{n_{m^{\prime}}\pm 1,n^{\prime}},$ (181)
$\displaystyle\bar{\mathcal{G}}_{m^{\prime}}^{(0)}$
$\displaystyle=G_{n_{m^{\prime}},n^{\prime}},$ (182)

we first derive the equation

$\displaystyle\left(\begin{array}[]{c}\bar{\mathcal{G}}^{(+)}\\\
\bar{\mathcal{G}}^{(-)}\\\
\bar{\mathcal{G}}^{(0)}\end{array}\right)=\left(\begin{array}[]{ccc}1&0&-T^{+,0}\\\
0&1&-T^{-,0}\\\ -T^{0,+}&-T^{0,-}&1\\\
\end{array}\right)^{-1}\left(\begin{array}[]{c}\bar{\mathcal{G}}_{disj}^{(+)}\\\
\bar{\mathcal{G}}_{disj}^{(-)}\\\
\bar{\mathcal{G}}_{disj}^{(0)}\end{array}\right),$ (192)

where

$\displaystyle\bar{\mathcal{G}}_{disj,m^{\prime}}^{(+)}$
$\displaystyle=G_{n_{m^{\prime}}+1,n^{\prime}}^{m^{\prime}},$ (193)
$\displaystyle\bar{\mathcal{G}}_{disj,m^{\prime}}^{(-)}$
$\displaystyle=G_{n_{m^{\prime}}-1,n^{\prime}}^{m^{\prime}-1},$ (194)
$\displaystyle\bar{\mathcal{G}}_{disj,m^{\prime}}^{(0)}$
$\displaystyle=g_{n_{m^{\prime}},n_{m^{\prime}}}^{m^{\prime}}\delta_{n_{m^{\prime}},n^{\prime}},$
(195)

and

$\displaystyle T_{m,m^{\prime}}^{0,+}$
$\displaystyle=-\delta_{m,m^{\prime}}g^{m^{\prime}}_{n_{m^{\prime}},n_{m^{\prime}}}t^{\dagger}_{n_{m^{\prime}}},$
(196) $\displaystyle T_{m,m^{\prime}}^{0,-}$
$\displaystyle=-\delta_{m,m^{\prime}}g^{m^{\prime}}_{n_{m^{\prime}},n_{m^{\prime}}}t_{n_{m^{\prime}}-1},$
(197) $\displaystyle T_{m,m^{\prime}}^{+,0}$
$\displaystyle=-\delta_{m,m^{\prime}-1}G_{n_{m^{\prime}-1}+1,n_{m^{\prime}}-1}^{n_{m^{\prime}}-1}t_{n_{m^{\prime}}-1}^{\dagger}$
$\displaystyle-\delta_{m,m^{\prime}}G_{n_{m^{\prime}}+1,n_{m^{\prime}}+1}^{m^{\prime}}t_{n_{m^{\prime}}},$
(198) $\displaystyle T_{m,m^{\prime}}^{-,0}$
$\displaystyle=-\delta_{m,m^{\prime}}G_{n_{m^{\prime}}-1,n_{m^{\prime}}-1}^{m^{\prime}-1}t_{n_{m^{\prime}}-1}^{\dagger}$
$\displaystyle-\delta_{m,m^{\prime}+1}G_{n_{m^{\prime}+1}-1,n_{m^{\prime}}+1}^{m^{\prime}}t_{n_{m^{\prime}}}.$
(199)

Finding explicitly the solutions for (181), (182)

$\displaystyle
G_{n_{m},n^{\prime}}=-\sum_{m^{\prime}=1}^{M}(D^{-1})_{m,m^{\prime}}\bar{F}_{m^{\prime},n^{\prime}},$
(200) $\displaystyle
g^{m}_{n_{m},n_{m}}\left(t_{n_{m}-1}G_{n_{m}-1,n^{\prime}}+t^{\dagger}_{n_{m}}G_{n_{m}+1,n^{\prime}}\right)$
$\displaystyle=g_{n_{m},n_{m}}^{m}\delta_{n_{m},n^{\prime}}+\sum_{m^{\prime}=1}^{M}(D^{-1})_{m,m^{\prime}}\bar{F}_{m^{\prime},n^{\prime}},$
(201)

we insert them into (180). Thereby we establish the composite Green’s function
of the whole system

$\displaystyle
G_{n,n^{\prime}}=\sum_{m=0}^{M}G_{n,n^{\prime}}^{m}+\sum_{m,m^{\prime}=1}^{M}F_{n,m}(D^{-1})_{m,m^{\prime}}\bar{F}_{m^{\prime},n^{\prime}},$
(202)

where

$\displaystyle D_{m,m^{\prime}}=\delta_{m,m^{\prime}}\,(z-W_{n_{m}}$
$\displaystyle-
t_{n_{m}-1}G^{m-1}_{n_{m}-1,n_{m}-1}t^{\dagger}_{n_{m}-1}-t_{n_{m}}^{\dagger}G^{m}_{n_{m}+1,n_{m}+1}t_{n_{m}})$
$\displaystyle-\delta_{m,m^{\prime}+1}t_{n_{m}-1}G^{m-1}_{n_{m}-1,n_{m-1}+1}t_{n_{m-1}}$
$\displaystyle-\delta_{m,m^{\prime}-1}t^{\dagger}_{n_{m}}G^{m}_{n_{m}+1,n_{m+1}-1}t_{n_{m+1}-1}^{\dagger},$
(203)

and

$\displaystyle F_{n,m}$
$\displaystyle=-\delta_{n,n_{m}}+G_{n,n_{m}-1}^{m-1}t_{n_{m}-1}^{\dagger}+G_{n,n_{m}+1}^{m}t_{n_{m}},$
(204) $\displaystyle\bar{F}_{m^{\prime},n^{\prime}}$
$\displaystyle=-\delta_{n_{m^{\prime}},n^{\prime}}+t_{n_{m^{\prime}}-1}G_{n_{m^{\prime}}-1,n^{\prime}}^{m^{\prime}-1}+t^{\dagger}_{n_{m^{\prime}}}G_{n_{m^{\prime}}+1,n^{\prime}}^{m^{\prime}}.$
(205)

The formula (202) is a multi-barrier generalization of the single-barrier
formula (15). Observing that

$\displaystyle\sum_{n}\bar{F}_{m^{\prime},n}F_{n,m}$
$\displaystyle=\frac{\partial}{\partial\omega}D_{m^{\prime},m},$ (206)

we also express the global DOS in the multi-barrier setup by the analogy with
the single-barrier formula (25) as

$\displaystyle\rho(\omega)=\sum_{m=0}^{M}\rho^{m}(\omega)-\frac{1}{\pi}\text{Im}\,\frac{\partial}{\partial\omega}\ln\det
D(\omega+i0^{+}).$ (207)

To take the continuum limit of $D_{m,m^{\prime}}$ we approximate in (203) [cf.
Eq. (29)]

$\displaystyle
W_{n_{m}}\approx\frac{\mathcal{A}_{m-1}+\mathcal{A}_{m}}{2a^{2}}+\frac{\mathcal{U}_{m}}{a}$
(208)

and

$\displaystyle t_{n_{m}-1}G^{m-1}_{n_{m}-1,n_{m}-1}t^{\dagger}_{n_{m}-1}$
$\displaystyle\approx-\frac{\mathcal{A}_{m-1}}{2a^{2}}+\frac{1}{8a}\mathcal{A}_{m-1}\lim_{x\to
x_{m}^{-}}\frac{d^{2}}{dx^{2}}G^{m-1}(x,x)\mathcal{A}_{m-1},$ (209)
$\displaystyle t_{n_{m}}^{\dagger}G^{m}_{n_{m}+1,n_{m}+1}t_{n_{m}}$
$\displaystyle\approx-\frac{\mathcal{A}_{m}}{2a^{2}}+\frac{1}{8a}\mathcal{A}_{m}\lim_{x\to
x_{m}^{+}}\frac{d^{2}}{dx^{2}}G^{m}(x,x)\mathcal{A}_{m},$ (210)

analogously to the single-barrier case discussed in Appendix B. In addition,
there are two new objects in (203), which are approximated by

$\displaystyle t_{n_{m}-1}G^{m-1}_{n_{m}-1,n_{m-1}+1}t_{n_{m-1}}$
$\displaystyle\approx-\frac{1}{4a}\mathcal{A}_{m-1}G_{12}^{m-1}(x_{m}^{-},x_{m-1}^{+})\mathcal{A}_{m-1},$
(211) $\displaystyle
t^{\dagger}_{n_{m}}G^{m}_{n_{m}+1,n_{m+1}-1}t_{n_{m+1}-1}^{\dagger}$
$\displaystyle\approx-\frac{1}{4a}\mathcal{A}_{m}G^{m}_{12}(x_{m}^{+},x_{m+1}^{-})\mathcal{A}_{m}.$
(212)

Altogether these terms provide the continuum limit of $D$ (203) giving rise to
the expressions (58)-(60) defining the matrix $d=\lim_{a\to 0}a\,D$ in the
multi-barrier case.

A derivation of the continuum limit of (204) and (205) employs the same type
of the $a$-expansion as in (162)-(165). It leads us to the final result (61),
(62) analogous to the result (32), (33) of the single-barrier case.

##  Appendix E Simplification of the Green’s function in the double-barrier
case

On the basis of expressions (148) and (152) we obtain

$\displaystyle G_{0,W}^{C}(x,x^{\prime})$ $\displaystyle=g(x-x^{\prime})$
(213) $\displaystyle-g(x)\left\\{g(0)-g(-W)[g(0)]^{-1}g(W)\right\\}^{-1}$
$\displaystyle\times\left\\{g(-x^{\prime})-g(-W)[g(0)]^{-1}g(W-x^{\prime})\right\\}$
$\displaystyle-g(x-W)\left\\{g(0)-g(W)[g(0)]^{-1}g(-W)\right\\}^{-1}$
$\displaystyle\times\left\\{g(W-x^{\prime})-g(W)[g(0)]^{-1}g(-x^{\prime})\right\\},$

in terms of the translationally invariant counterpart $G^{(C,0)}(x,0)\equiv
g(x)$. According to (155), the boundary Green’s function of the central region
is then given by $G^{C}(x,x^{\prime})=\Theta(W>x>0)G_{0,W}^{C}(x,x^{\prime})$.

Taking in account that $G_{0,W}^{C}(x,x^{\prime})$ identically vanishes when
$x$ and $x^{\prime}$ belong to different spatial regions [either $(-\infty,0)$
or $(0,W)$ or $(W,\infty)$], we get the identities

$\displaystyle\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{-},W^{+})=\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{-},W^{-})$ (214)
$\displaystyle=\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{+},W^{+})=0.$ (215)

They allow us to replace

$\displaystyle\frac{1}{4}\mathcal{A}_{C}G_{12}^{C}(0^{+},W^{-})\mathcal{A}_{C}$
(216) $\displaystyle=$
$\displaystyle\frac{1}{4}\mathcal{A}_{C}\left[\frac{\partial^{2}}{\partial
x\partial x^{\prime}}G_{0,W}^{C}(0^{+},W^{-})-\frac{\partial^{2}}{\partial
x\partial x^{\prime}}G_{0,W}^{C}(0^{-},W^{-})\right.$
$\displaystyle\quad\left.+\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{-},W^{+})-\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{+},W^{+})\right]\mathcal{A}_{C}.$

After the lengthy calculation using (213) and (42), (43) we obtain

$\displaystyle\frac{1}{4}\mathcal{A}_{C}G_{12}^{C}(0^{+},W^{-})\mathcal{A}_{C}$
$\displaystyle=\left\\{g(W)-g(0)[g(-W)]^{-1}g(0)\right\\}^{-1}.$ (217)

Analogously we deduce

$\displaystyle\frac{1}{4}\mathcal{A}_{C}G_{12}^{C}(W^{-},0^{+})\mathcal{A}_{C}$
$\displaystyle=\left\\{g(-W)-g(0)[g(W)]^{-1}g(0)\right\\}^{-1}.$ (218)

To derive

$\displaystyle-\frac{1}{8}\mathcal{A}_{C}\left[\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G_{0}^{C}(x,x)+\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{C}(x,x)\right]\mathcal{A}_{C}$
$\displaystyle=\left\\{g(0)-g(-W)[g(0)]^{-1}g(W)\right\\}^{-1},$ (219)

we use the identities

$\displaystyle\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{+},0^{-})=\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{-},0^{+})=0$ (220)

to replace

$\displaystyle-\frac{1}{8}\mathcal{A}_{C}\left[\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G_{0}^{C}(x,x)+\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G^{C}(x,x)\right]\mathcal{A}_{C}$ (221)
$\displaystyle=$ $\displaystyle-\frac{1}{8}\mathcal{A}_{C}\left[\lim_{x\to
0^{-}}\frac{d^{2}}{dx^{2}}G_{0,W}^{C}(x,x)+\lim_{x\to
0^{+}}\frac{d^{2}}{dx^{2}}G_{0,W}^{C}(x,x)\right.$
$\displaystyle\left.-2\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{+},0^{-})-2\frac{\partial^{2}}{\partial x\partial
x^{\prime}}G_{0,W}^{C}(0^{+},0^{-})\right]\mathcal{A}_{C}.$

After the lengthy calculation using (213) and (42), (43) we obtain (219).

Analogously we deduce

$\displaystyle-\frac{1}{8}\mathcal{A}_{C}\left[\lim_{x\to
W^{+}}\frac{d^{2}}{dx^{2}}G_{W}^{C}(x,x)+\lim_{x\to
W^{-}}\frac{d^{2}}{dx^{2}}G^{C}(x,x)\right]\mathcal{A}_{C}$
$\displaystyle=\left\\{g(0)-g(W)[g(0)]^{-1}g(-W)\right\\}^{-1}.$ (222)

##  Appendix F Josephson system

###  F.1 Position-space Green’s functions of bulk $s$-wave superconductors

The momentum-space Green’s functions of the bulk $s$-wave superconductors are
equal at $\varphi=0$ to

$\displaystyle G_{k}^{(0)}(z)=$
$\displaystyle\frac{1}{z^{2}-(h_{k}^{(0)})^{2}-\Delta^{2}}\left(\begin{array}[]{cc}z+h_{k}^{(0)}&\Delta_{0}\\\
\Delta_{0}&z-h_{k}^{(0)}\end{array}\right),$ (225)

where $h_{k}^{(0)}=\frac{k^{2}}{2m}-\mu$.

For the evaluation of

$\displaystyle
G^{(0)}(x,x^{\prime})=\int_{-\infty}^{\infty}\frac{dk}{2\pi}e^{ik(x-x^{\prime})}G_{k}^{(0)}(z)$
(226)

we note the following poles of $G_{k}^{(0)}(z)$

$\displaystyle
k_{1,\pm}=\pm\sqrt{2m\mu+2mz\sqrt{1-\left(\frac{\Delta_{0}}{z}\right)^{2}}},$
(227) $\displaystyle
k_{2,\pm}=\mp\sqrt{2m\mu-2mz\sqrt{1-\left(\frac{\Delta_{0}}{z}\right)^{2}}}.$
(228)

Considering $\text{Im}(z)>0$ and deforming appropriately the integration
contour in the complex $k$-plane, we obtain by means of the residue theorem

$\displaystyle G^{(0)}(x,x^{\prime})=$
$\displaystyle-\frac{ime^{ik_{1,+}|x-x^{\prime}|}}{2k_{1,+}}\left[e^{\tau_{x}\chi}+\tau_{z}\right]$
(229)
$\displaystyle+\frac{ime^{ik_{2,+}|x-x^{\prime}|}}{2k_{2,+}}\left[e^{\tau_{x}\chi}-\tau_{z}\right],$
(230)

where $e^{\tau_{x}\chi}=\frac{1+\tau_{x}\tanh\chi}{\sqrt{1-\tanh^{2}\chi}}$
and $\tanh\chi=\frac{\Delta_{0}}{z}$.

For $\text{Im}(z)<0$ it is necessary to replace $k_{1,+}\to k_{1,-}$ and
$k_{2,+}\to k_{2,-}$ in the above equations.

Next, we evaluate

$\displaystyle G^{(0)}(0,0)$
$\displaystyle=-\frac{im}{k_{1,+}k_{2,+}}\left[k^{(+)}\tau_{z}-k^{(-)}e^{\tau_{x}\chi}\right],$
(231) $\displaystyle\,[G^{(0)}(0,0)]^{-1}$
$\displaystyle=-\frac{k^{(+)}\tau_{z}+k^{(-)}e^{-\tau_{x}\chi}}{im},$ (232)

where

$\displaystyle k^{(\pm)}=\frac{k_{1,+}\pm
k_{2,+}}{2}\,\text{sgn}\,\text{Im}(z),$ (233)

as well as

$\displaystyle G_{1}^{(0)}(0^{\pm},0)=G_{2}^{(0)}(0,0^{\pm})=\pm m\tau_{z}.$
(234)

With this in hands and noting

$\displaystyle
G^{(0,\lambda)}(x,x^{\prime})=U_{\lambda}G^{(0)}(x,x^{\prime})U_{\lambda}^{\dagger},\quad
U_{\lambda}=e^{\frac{i}{4}\tau_{z}\lambda\varphi},$ (235)

one can show on the basis of Eqs. (54), (55), and (56) that the matrix $d$ of
the Josephson system is given by Eq. (102).

###  F.2 Local Coulomb interaction in the Hartree-Fock approximation

To compute the expectation values
$\langle\psi_{\downarrow}(0)\psi_{\uparrow}(0)\rangle$ and
$\langle\psi_{\uparrow}^{\dagger}(0)\psi_{\downarrow}^{\dagger}(0)\rangle$
with respect to the Hartree-Fock Hamiltonian, we define the following
imaginary-time Matsubara Green’s functions

$\displaystyle G_{\gamma\delta;U}^{\text{Mat}}(x,x^{\prime};\tau)=$
$\displaystyle-\langle
T_{\tau}\hat{\Psi}_{\gamma}(x,\tau)\hat{\Psi}_{\delta}^{\dagger}(x^{\prime},0)\rangle$
(236) $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{i\omega_{n}}e^{-i\omega_{n}\tau}G_{\gamma\delta;U}(x,x^{\prime};i\omega_{n}),$
(237)

where $\hat{\Psi}_{1}(x)=\hat{\psi}_{\uparrow}(x)$ and
$\hat{\Psi}_{2}(x)=\hat{\psi}_{\downarrow}^{\dagger}(x)$. In this terms

$\displaystyle\langle\psi_{\downarrow}(0)\psi_{\uparrow}(0)\rangle$
$\displaystyle=G^{\text{Mat}}_{12,U}(0,0;\tau=0^{-})$
$\displaystyle=\frac{1}{2\beta}\sum_{i\omega_{n}}e^{i\omega_{n}0^{+}}\text{tr}\,\left\\{[d_{U}(i\omega_{n})]^{-1}\tau_{x}\right\\},$
(238)
$\displaystyle\langle\psi_{\uparrow}^{\dagger}(0)\psi_{\downarrow}^{\dagger}(0)\rangle$
$\displaystyle=G^{\text{Mat}}_{21,U}(0,0;\tau=0^{-})$
$\displaystyle=\frac{1}{2\beta}\sum_{i\omega_{n}}e^{i\omega_{n}0^{+}}\text{tr}\,\left\\{[d_{U}(i\omega_{n})]^{-1}\tau_{x}\right\\}.$
(239)

The function $[d_{U}(z)]^{-1}=[d(z)-\Delta_{loc}\tau_{x}]^{-1}$ has the
following off-diagonal components

$\displaystyle\frac{1}{2}\text{tr}\,\left\\{[d_{U}(z)]^{-1}\tau_{x}\right\\}=-\frac{m}{k_{F}}\left(\frac{\frac{k^{(-)}}{k_{F}}\frac{i\Delta_{0}}{z}\cos\frac{\varphi}{2}}{\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}+\bar{\Delta}_{loc}\right)$
(240)
$\displaystyle\times\frac{1}{\frac{\frac{k^{(-)\,2}}{k_{F}^{2}}}{1-(\frac{\Delta_{0}}{z})^{2}}+\left(\frac{\frac{k^{(-)}}{k_{F}}\frac{i\Delta_{0}}{z}\cos\frac{\varphi}{2}}{\sqrt{1-(\frac{\Delta_{0}}{z})^{2}}}+\bar{\Delta}_{loc}\right)^{2}+(\bar{V}_{0}+\frac{k^{(+)}}{ik_{F}})^{2}},$

where $\bar{V}_{0}=\frac{mV_{0}}{k_{F}}$ and
$\bar{\Delta}_{loc}=\frac{m\Delta_{loc}}{k_{F}}$.

At zero temperature the Matsubara sum in (239) turns into the imaginary-
frequency integral. The integrand $\propto\Delta_{0}\cos\frac{\varphi}{2}$ in
(240) decays sufficiently fast at high frequencies, so that the convergence
factor $e^{i\omega 0^{+}}$ for it may be omitted. The other integrand
$\propto\bar{\Delta}_{loc}$ does require the convergence factor, which is
needed to handle the terms $\sim\frac{1}{z}$ at high frequencies. Collecting
all contributions we obtain the following self-consistency equation

$\displaystyle\bar{\Delta}_{loc}=-u\int_{0}^{\infty}d\bar{\omega}\frac{\bar{k}^{(-)}\cos\frac{\varphi}{2}}{\sqrt{\bar{\omega}^{2}+1}}R(\bar{\omega}),$
(241)
$\displaystyle\qquad-u\bar{\Delta}_{loc}\int_{0}^{\infty}d\bar{\omega}[R(\bar{\omega})-R_{0}(\bar{\omega})]$
(242)
$\displaystyle\qquad-u\bar{\Delta}_{loc}\int_{0}^{\infty}d\bar{\omega}R_{0}(\bar{\omega})\cos(\bar{\omega}0^{+}),$
(243)

where

$\displaystyle u=\frac{Um^{2}\Delta_{0}}{\pi k_{F}^{2}}$ (244)

is a dimensionless interaction parameter, and

$\displaystyle R(\bar{\omega})$
$\displaystyle=\frac{1}{\frac{\bar{k}^{(-)\,2}\bar{\omega}^{2}}{\bar{\omega}^{2}+1}+\left(\frac{\bar{k}^{(-)}\cos\frac{\varphi}{2}}{\sqrt{\bar{\omega}^{2}+1}}+\bar{\Delta}_{loc}\right)^{2}+\left(\bar{V}_{0}-i\bar{k}^{(+)}\right)^{2}},$
(245) $\displaystyle\bar{k}^{(\pm)}$
$\displaystyle=\frac{\sqrt{1+i\frac{\Delta_{0}}{\mu}\sqrt{\bar{\omega}^{2}+1}}\mp\sqrt{1-i\frac{\Delta_{0}}{\mu}\sqrt{\bar{\omega}^{2}+1}}}{2},$
(246) $\displaystyle R_{0}(\bar{\omega})$
$\displaystyle=\frac{1}{\bar{k}_{0}^{(-)\,2}+\left(\bar{V}_{0}-i\bar{k}_{0}^{(+)}\right)^{2}},$
(247) $\displaystyle\bar{k}_{0}^{(\pm)}$
$\displaystyle=\frac{\sqrt{1+i\frac{\Delta_{0}}{\mu}\bar{\omega}}\mp\sqrt{1-i\frac{\Delta_{0}}{\mu}\bar{\omega}}}{2}.$
(248)

The perform the integral in (243) we rotate the integration contour to the
real axis

$\displaystyle
u\bar{\Delta}_{loc}\,i\int_{i0}^{i\infty}d\bar{z}R_{0}(-i\bar{z})\frac{e^{\bar{z}0^{+}}+e^{-\bar{z}0^{+}}}{2}$
$\displaystyle=$
$\displaystyle-\frac{i}{2}u\bar{\Delta}_{loc}\int_{0}^{\infty}d\bar{\omega}\,[R_{0}(i\bar{\omega}+0^{+})-R_{0}(-i\bar{\omega}+0^{+})]$
(249) $\displaystyle=$
$\displaystyle-u\bar{\Delta}_{loc}\int_{\mu/\Delta_{0}}^{\infty}\frac{d\bar{\omega}\,e^{-\bar{\omega}0^{+}}\,\sqrt{\frac{\Delta_{0}}{\mu}\bar{\omega}+1}}{\left(\sqrt{\frac{\Delta_{0}}{\mu}\bar{\omega}-1}+\bar{V}_{0}\right)\left(\frac{\Delta_{0}}{\mu}\bar{\omega}+1+\bar{V}_{0}^{2}\right)}.$
(250)

The resulting integral has the ultraviolet logarithmic divergence which is not
cut off at the scale $\mu$. The origin of this divergence is rooted in the
ultra-local form of the Coulomb interaction in our model, and for the
integral’s regularization it is necessary to introduce the cut-off scale
$\omega_{c}\gg\mu$, whose inverse $v_{F}/\omega_{c}$ captures the microscopic
length scale of the interaction. Thus (243) is estimated by
$-u\bar{\Delta}_{loc}\frac{\mu}{\Delta_{0}}\ln\frac{\omega_{c}}{\Delta_{0}}$.
This behaviour signals the necessity to refine the model’s consideration, e.g.
by resorting to renormalization group methods.Meden (2019)

In turn, in the bare perturbation theory the local superconducting order
parameter is given by the convergent integral

$\displaystyle\bar{\Delta}_{loc}\approx-u\int_{0}^{\infty}d\bar{\omega}\frac{\bar{k}^{(-)}\cos\frac{\varphi}{2}}{\sqrt{\bar{\omega}^{2}+1}}R(\bar{\omega})\bigg{|}_{\bar{\Delta}_{loc}=0}.$
(251)

It is remarkable that in the Andreev limit $\mu\gg\Delta_{0}$ it also features
the logarithmic dependence of the high energy scale $\mu$:

$\displaystyle\bar{\Delta}_{loc}\approx-
Du\cos\frac{\varphi}{2}\ln\frac{\mu}{\Delta_{0}},\quad
D=\frac{1}{1+\bar{V}_{0}^{2}}.$ (252)