# Multiterminal ballistic Josephson junctions coupled to normal leads

Régis Mélin  Univ. Grenoble-Alpes, CNRS, Grenoble INP††thanks: Institute of
Engineering Univ. Grenoble Alpes, Institut NEEL, 38000 Grenoble, France

###### Abstract

Multiterminal Josephson junctions have aroused considerable theoretical
interest recently and numerous works aim at putting the predictions of
correlations among Coopers (i.e. the so-called quartets) and simulation of
topological matter to the test of experiments. The present paper is motivated
by recent experimental investigation from the Harvard group reporting
$h/4e$-periodic quartet signal in a four-terminal configuration containing a
loop pierced by magnetic flux, together with inversion controlled by the bias
voltage, i.e. the quartet signal can be larger at half-flux quantum than in
zero magnetic field. Here, we theoretically focus on devices consisting of
finite-size quantum dots connected to four superconducting and to a normal
lead. In addition to presenting numerical calculations of the quartet signal
within a simplified modeling, we reduce the device to a non-Hermitian
Hamiltonian in the infinite-gap limit. Then, relaxation has the surprising
effect of producing sharp peaks and log-normal variations in the voltage-
sensitivity of the quartet signal in spite of the expected moderate
fluctuations in the two-terminal DC-Josephson current. The phenomenon is
reminiscent of a resonantly driven harmonic oscillator having amplitude
inverse proportional to the damping rate, and of the thermal noise of a
superconducting weak link, inverse proportional to the Dynes parameter.
Perspectives involve quantum bath engineering multiterminal Josephson
junctions in the circuits of cavity quantum electrodynamics.

##  I Introduction

BCS superconductors BCS  are characterized by macroscopic phase variable
$\varphi$ and energy gap $\Delta$ between the ground state and the first
quasiparticles. Anderson demonstrated Anderson-RPA ; Anderson-mass  that
Coulomb interaction produces collective modes in bulk superconductors,
yielding the Higgs mechanism as an explanation to the Meissner effect.
Josephson Josephson  demonstrated that gauge invariance implies DC-
supercurrent between two phase-biased superconductors, and AC-current
oscillations with voltage-biasing. The so-called Andreev bound states (ABS)
Andreev ; deGennes ; SaintJames ; Kulik  contribute to the DC-Josephson
supercurrent flowing through any type of weak link, see for instance Refs.
Likharev, ; AVERIN, ; Cuevas, ; Scheer, ; Beenakker, . Recently, the ABS were
probed with microwave spectroscopy Bretheau1 ; Bretheau2 ; Olivares ; Janvier
, but the life-time of this “Andreev qu-bit” Zazunov  turns out not to be
infinite. In general, understanding the mechanisms of quasiparticle poisoning
Martinis2009 ; deVisser2011 ; Lenander2011 ; Rajauria2012 ; Wenner2013 ;
Riste2013 ; LevensonFalk2014 ; Nazarov-qp  is a central issue in the studies
of circuits of quantum engineering Kouznetsov ; Clarke1 ; Clarke2 ; Devoret ;
Martinis . In bulk superconductors, the electron-phonon coupling or the
electron-electron repulsive Coulomb interaction produce finite quasiparticle
life-time Kaplan ; Dynes ; Pekola1 ; Pekola2  captured by adding small
imaginary part $\eta$ to their energy. This small Dynes parameter $\eta$
Kaplan ; Dynes ; Pekola1 ; Pekola2  produces exponential time-decay of the
quasiparticle wave-functions. In the following paper, we investigate
multiterminal Josephson junctions coupled to normal leads, and demonstrate
that small relaxation produced by the coupling to those normal conductors has
drastic effect on the current. Based on previous works related to infinite-gap
Hamiltonians Zazunov ; Meng ; topo1_plus_Floquet , we also here-consider the
infinite-gap limit, where non-Hermitian Hamiltonians emerge for quantum dots
coupled to superconducting and normal leads. Those non-Hermitian Hamiltonians
receive interpretation of describing the coherent oscillations and damping in
the ABS dynamics, and they capture the degrees of freedom of the ABS having
spatial extent set by the BCS coherence length, which is vanishingly small in
the infinite-gap limit. The present paper also suggests the interest of future
quantum bath engineering for multiterminal Josephson junctions connected to
cavity-quantum electrodynamics resonators.

Now, we specifically introduce the considered four-terminal Josephson
junctions, starting with the three-terminal Cooper pair beam splitters that
have been the subject of intense investigations since more than twenty years.
Entangled Andreev pairs can split if two independent voltages are applied in
three-terminal $F_{a}SF_{b}$ ferromagnet-superconductor-ferromagnet or
$N_{a}SN_{b}$ normal metal-superconductor-normal metal Cooper pair beam
splitters theory-CPBS1 ; theory-CPBS2 ; theory-CPBS3 ; theory-CPBS4 ; theory-
CPBS4-bis ; theory-CPBS4-ter ; theory-CPBS5 ; theory-CPBS6 ; theory-CPBS7 ;
theory-CPBS8 ; theory-CPBS9 ; theory-CPBS10 ; theory-CPBS11 ; Yeyati2007 ;
Danneau ; Beckmann2007 ; Schindele2012 ; Schindele2014 ; Tan2015 ;
Borzenets2016a ; exp-CPBS1 ; exp-CPBS2 ; exp-CPBS3 ; exp-CPBS4 ; exp-CPBS5 ;
exp-CPBS6 ; exp-CPBS7 ; exp-CPBS8 . “Nonlocal” or “crossed” Andreev reflection
(CAR) appears if the separation between the $N_{a}S$ and $SN_{b}$ interfaces
is comparable to the zero-energy superconducting coherence length. Spin-up
electron from lead $N_{b}$ can be Andreev reflected as spin-down hole into
$N_{a}$, leaving a Cooper pair in the “central” $S$. This yields nonlocal
current response reflecting production of nonlocally split Cooper pairs. In
addition, the positive current-current cross-correlations of Cooper pair
splitting theory-CPBS7 ; theory-noise1 ; theory-noise2 ; theory-noise3 ;
theory-noise4 ; theory-noise5 ; theory-noise6 ; theory-noise7 ; theory-noise8
; theory-noise9 ; theory-noise10 ; theory-noise11 ; theory-noise12  were
experimentally revealed exp-CPBS7 ; exp-CPBS8 .

The double quantum dot experiments exp-CPBS8 ; exp-CPBS6 ; exp-CPBS5  provide
evidence for nonlocal two-particle Andreev resonance in
$N_{a}$-dot-$S_{c}$-dot-$N_{b}$ devices, on the condition of the opposite
energy levels $\epsilon_{a}=-\epsilon_{b}$ on both quantum dots. In the
following, we demonstrate how “nonlocal two-particle Andreev resonance” in
$N_{a}$-dot-$S_{c}$-dot-$N_{b}$ Cooper pair beam splitters can be generalized
to “nonlocal two-Cooper pair resonance” in all-superconducting three-terminal
Josephson junctions.

Those two-Cooper pair resonances build on nonlocal two-Cooper pair states, the
so-called quartets Cuevas-Pothier ; Freyn ; Melin1 ; Jonckheere ; FWS ; Sotto
; engineering ; paperI ; paperII ; long-distance ; Akh2 ; Pillet ; Pillet2 ;
Nazarov-PRR ; Nazarov-PRB-AM  that appear in $(S_{a},\,S_{c},\,S_{b})$ three-
terminal Josephson junctions, where $S_{a}$ and $S_{b}$ are voltage-biased at
$V_{a}$ and $V_{b}$, the superconducting lead $S_{c}$ being grounded at
$V_{c}=0$. Andreev scattering yields quartet phase-sensitive DC-current
response if the condition $V_{a}=-V_{b}\equiv V$ is fulfilled. On this quartet
line $V_{a}=-V_{b}$, the elementary transport process transfers two Cooper
pairs from $S_{a}$ and $S_{b}$ into the grounded $S_{c}$, while exchanging
partners Cuevas-Pothier ; Freyn ; Melin1 ; Jonckheere ; FWS ; Sotto ;
engineering ; paperI ; paperII ; long-distance ; Akh2 ; Pillet ; Pillet2 ;
Nazarov-PRR ; Nazarov-PRB-AM . By the time-energy uncertainty relation, this
phase-sensitive quartet current is DC, since, on the quartet line, the energy
$E_{i}=2e(V_{a}+V_{b})$ of two Cooper pairs from $S_{a}$ and $S_{b}$ in the
initial state is equal to the final state energy $E_{f}=4eV_{c}$ of the two
Cooper pairs transmitted into $S_{c}$.

Several experiments were recently performed on multiterminal Josephson
junctions Lefloch ; Heiblum ; HGE ; multiterminal-exp1 ; multiterminal-exp2 ;
multiterminal-exp3 ; multiterminal-exp4 ; multiterminal-exp5 ; multiterminal-
exp6 ; multiterminal-exp7 ; multiterminal-exp8 ; multiterminal-exp9 . Some of
those explored the possibility of nontrivial topology vanHeck ; Padurariu ;
multiterminal-exp1 ; Nazarov1 ; Nazarov2 ; topo0 ; topo1 ; topo2 ; topo3 ;
topo4 ; topo4-bis ; topo5 ; Berry ; Feinberg1 ; Feinberg2 ; topo1_plus_Floquet
; Levchenko1 ; Levchenko2 ; Akh2 . On the other hand, three groups reported
compatibility with the quartets Cuevas-Pothier ; Freyn ; Melin1 ; Jonckheere ;
FWS ; Sotto ; engineering ; paperI ; paperII ; long-distance ; Akh2 ; Pillet ;
Pillet2 ; Nazarov-PRR ; Nazarov-PRB-AM : the Grenoble group experiment with
metallic structures Lefloch , the Weizmann Institute group experiment with
semiconducting nanowires Heiblum  and the more recent Harvard group experiment
on ballistic graphene-based four-terminal Josephson junctions HGE . This third
experiment HGE  is summarized in figures 1a and 1b. The latter shows typical
experimental variations for the quartet critical current as a function of the
bias voltage $V$ and the magnetic flux $\Phi$ through the loop. The experiment
features the counterintuitive voltage-$V$-tunable inversion, i.e. the
possibility of stronger quartet critical current at half-flux quantum
$\Phi=\pi$ than in zero field $\Phi=0$, see figure 1b. The present paper is
motivated by exploring interpretation of this experiment, in the continuation
of our previous papers I and II paperI ; paperII . We also underline two
recent experimental preprints reporting the quartet line in ballistic devices
multiterminal-exp8 ; multiterminal-exp9 .

Figure 1: Schematic four-terminal experiment: Panel a represents top-view of
the experimental device, with the superconducting leads $S_{a}$ and $S_{b}$
biased at $V_{a,b}=\pm V$. The loop terminated by $S_{c,1}$ and $S_{c,2}$ is
connected to ground by the contact on top of it, and it is pierced by the
magnetic field flux $\Phi$. The four superconducting leads are evaporated on a
graphene layer gated far from the Dirac point, and forming a two-dimensional
(2D) metal. Panel b summarizes the experimental result HGE  that inspires the
present theoretical paper, i.e. the quartet critical current current
$I_{q,c}(\pi)$ at flux $\Phi=\pi$ can be larger or smaller than $I_{q,c}(0)$
at flux $\Phi=0$, in a way that is sensitive to voltage typically being in the
range $eV\approx\Delta/10$ where $\Delta$ is the superconducting gap.

Zero-dimensional (0D) quantum dots were used in Ref. Park,  to model two-
terminal graphene-superconductor hybrids, in connection with experimental
evidence for the Floquet replica of the Andreev spectrum under microwave
radiation. The goal of the present paper is to similarly describe four-
terminal graphene-based Josephson junctions with effective models based on
quantum dots, in the absence of Coulomb interaction. We argue that multilevel
quantum dots coupled to four superconducting and to a normal lead are “minimal
models” of four-terminal Josephson junctions that are intermediate between the
short- and long-junction limits. In addition, we further simplify the
description in terms of double 0D quantum dots. Four or more ABS are within
the gap in the equilibrium limit, instead of two ABS for a single 0D quantum
dot.

Importantly, nonproximitized regions can appear in the two-dimensional (2D)
metal, typically in the “cross-like region” formed in between the contacts
with the superconductors, see figure 1a. Then, the quantum transport process
can have initial or final states in those normal regions of the circuit.
Current lines can be converted by Andreev reflection as supercurrent flowing
in the grounded superconducting loop. Namely, current conservation for
symmetric coupling to the superconducting leads implies that the 2D metal
chemical potential is at the same zero energy as in the superconducting loop.
Thus, the normal carriers transmitted through the 2D metal can be transferred
into the grounded superconducting loop without taking or giving energy.

Now, we explain the connection between the present work and our previous
papers. Compared to our previous Ref. FWS, , we implement a double 0D quantum
dot, instead of the single quantum dot Dynes parameter Kaplan ; Dynes ;
Pekola1 ; Pekola2  model of Ref. FWS, . We also include four superconducting
leads, instead of three in our previous Ref. FWS, . Our previous paper I
paperI  presented an expansion of the quartet critical current in perturbation
in the tunneling amplitudes between a 2D metal and four superconducting leads.
The present paper addresses resonances that appear beyond the weak coupling
limit. Our previous paper II paperII  treated single quantum dot with
relaxation solely originating from the continua of BCS quasiparticles.
Relaxation then produces nonvanishingly small line-width broadening
$\delta_{0}$ for the Floquet resonances, which behaves like
$\log\delta_{0}\sim-\Delta/eV$ at low bias voltage energy $eV$ compared to the
superconducting gap $\Delta$ FWS ; engineering ; Berry ; paperII . In the
four-terminal device of paper II paperII , Landau-Zener quantum tunneling
reduces the quartet current by coherent superpositions in the dynamics of the
two opposite current-carrying ABS, at voltage values that are close to avoided
crossings in the Floquet spectra. This mechanism is analogous to the reduction
of the DC-Josephson current in microwave-irradiated weak links Bergeret .
Conversely, in the present work, the interplay between the time-periodic
dynamics and finite line-width broadening due to the attached normal lead
restores the expected sharp resonance peaks in the voltage-dependence of the
quartet current. For instance, a resonantly driven harmonic oscillator has
amplitude inverse proportional to the damping rate. In addition, we note that
the zero-frequency noise of a single superconducting weak link at finite
temperature is inverse proportional to the rate set by the Dynes parameter
Kaplan ; Dynes ; Pekola1 ; Pekola2  in some parameter window, see Ref. Yeyati-
noise-equilibrium, . The same scaling holds in this reference for the noise at
frequency $2\omega_{S}$, where $\omega_{S}$ is the ABS energy. The quartets in
four-terminal devices HGE ; paperI ; paperII  offer several parameters to
probe this physics, such as the quartet phase variable and the “knobs” of the
bias voltage and magnetic flux.

The paper is organized as follows. The Hamiltonians are presented in section
II. The mechanism for the two-Cooper pair resonance is discussed in section
III. The mechanism for the inversion is presented in section IV. Section V
shows numerical results. The low-voltage limit is discussed in section VI.
Concluding remarks are presented in the final section VII.

##  II Hamiltonians

In this section, we provide the Hamiltonians on which the paper is based.
Subsection II.1 presents the full Hamiltonians of the superconductors, quantum
dots and normal leads. It goes beyond the goal of the present paper to produce
numerical results for the voltage-biased four-terminal Josephson junctions
described by the general Hamiltonians of subsection II.1 in a realistic
geometry. This is why we present in subsection II.2 simplified model
Hamiltonians which will subsequently be used for numerical calculations.

###  II.1 General Hamiltonians

Now, we introduce in subsection II.1.1 the BCS Hamiltonian of the
superconducting leads. Subsections II.1.2 and II.1.3 introduce the
Hamiltonians of the quantum dots and normal leads respectively. As it is
mentioned above, the numerical calculations presented below are based on
simplifications of those Hamiltonians.

####  II.1.1 Superconductors

The Hamiltonian of a BCS superconductor is the following:

$\displaystyle{\cal H}_{BCS}$ $\displaystyle=$ $\displaystyle-W\sum_{\langle
i,j\rangle}\sum_{\sigma_{z}=\uparrow,\downarrow}\left(c_{i,\sigma_{z}}^{+}c_{j,\sigma_{z}}+c_{j,\sigma_{z}}^{+}c_{i,\sigma_{z}}\right)$
$\displaystyle-$
$\displaystyle\Delta\sum_{i}\left(\exp\left(i\varphi_{i}\right)c_{i,\uparrow}^{+}c_{i,\downarrow}^{+}+\exp\left(-i\varphi_{i}\right)c_{i,\downarrow}c_{i,\uparrow}\right),$

where $\sum_{\langle i,j\rangle}$ denotes summation over pairs of neighboring
tight-binding sites labeled by $i$ and $j$, and $\sigma_{z}$ is the component
of the spin along quantization axis. The bulk hopping amplitude is denoted by
$W$ and the superconducting gap $\Delta$ is taken identical in all
superconducting leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$. The variable
$\varphi_{i}$ denotes the superconducting phase variable at the tight-binding
site labeled by $i$. In the following, the current is weak and $\varphi_{i}$
is approximated as being uniform in space, with the values $\varphi_{a}$,
$\varphi_{b}$, $\varphi_{c,1}$, and $\varphi_{c,2}$ in $S_{a}$, $S_{b}$,
$S_{c,1}$ and $S_{c,2}$ respectively.

We assume short distance between the contact points $c_{1}$ and $c_{2}$ (at
the interfaces between $S_{c,1}$ or $S_{c,2}$ and the quantum dot), and we use
the approximation of the gauge

$\displaystyle\varphi_{c,1}$ $\displaystyle=$ $\displaystyle\varphi_{c}$ (2)
$\displaystyle\varphi_{c,2}$ $\displaystyle=$ $\displaystyle\varphi_{c}+\Phi,$
(3)

where $\Phi$ is the flux enclosed in the loop terminated by $S_{c,1}$ and
$S_{c,2}$.

####  II.1.2 Double quantum dots

A legitimate approximation is to discard the Coulomb electron-electron
repulsion at high transparency. Then the finite-size quantum dot is described
by the Hamiltonian

$\displaystyle{\cal H}_{dot}$ $\displaystyle=$ $\displaystyle-W\sum_{\langle
i,j\rangle}\sum_{\sigma_{z}=\uparrow,\downarrow}\left(c_{i,\sigma_{z}}^{+}c_{j,\sigma_{z}}+c_{j,\sigma_{z}}^{+}c_{i,\sigma_{z}}\right),$
(5)
$\displaystyle-\epsilon_{g}\sum_{i}\sum_{\sigma_{z}=\uparrow,\downarrow}c_{i,\sigma_{z}}^{+}c_{i,\sigma_{z}}$

where the hopping amplitude $W$ is identical to Eq. (II.1.1), and
$\epsilon_{g}$ is proportional to the gate voltage. The summation over pairs
of neighboring tight-binding sites is restricted to the region of finite
dimension defining the quantum dot.

Conversely, the Hamiltonian of a single 0D quantum dot $D_{x}$ at location
${\bf x}$ is the following:

${\cal
H}_{D_{x},0D}=\epsilon_{x}\sum_{\sigma_{z}}c_{D_{x},\sigma_{z}}^{+}c_{D_{x},\sigma_{z}},$
(6)

where $\epsilon_{x}$ is the on-site energy. Similarly, the double quantum dot
$(D_{x},D_{y})$ is characterized by the on-site energies $\epsilon_{x}$ and
$\epsilon_{y}$ with tunneling amplitude $\Sigma^{(0)}$ between them:

$\displaystyle{\cal H}_{D_{y},0D}$ $\displaystyle=$
$\displaystyle\epsilon_{y}\sum_{\sigma_{z}}c_{D_{y},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}$
(7) $\displaystyle{\cal H}_{T,0}$ $\displaystyle=$
$\displaystyle-\Sigma^{(0)}\sum_{\sigma_{z}}c_{D_{x},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}+h.c.$
(8)

At the exception of one of the numerical calculations presented in section V,
we use $\epsilon_{x}=\epsilon_{y}=0$ in the paper, thus with

$\displaystyle{\cal H}_{D_{x},0D}$ $\displaystyle=$ $\displaystyle 0$ (9)
$\displaystyle{\cal H}_{D_{y},0D}$ $\displaystyle=$ $\displaystyle 0.$ (10)

####  II.1.3 Normal leads

The Hamiltonian of a normal lead is given by

$\displaystyle{\cal H}_{N}$ $\displaystyle=$ $\displaystyle-W\sum_{\langle
i,j\rangle}\sum_{\sigma_{z}=\uparrow,\downarrow}\left(c_{i,\sigma_{z}}^{+}c_{j,\sigma_{z}}+c_{j,\sigma_{z}}^{+}c_{i,\sigma_{z}}\right)$

where the normal lead chemical potential is vanishingly small, and the hopping
amplitude $W$ is identical to Eq. (II.1.1). Eq. (II.1.3) is similar to the
quantum dot Hamiltonian given by Eq. (5) but now, the summation over the pairs
of neighboring tight-binding sites runs over an infinite lattice.

Figure 2: The four-terminal superconducting multilevel quantum dot: The figure
shows a multilevel quantum dot connected to the four superconducting leads
$S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$ biased at $V_{a,b}=\pm V$ and
$V_{c,1}=V_{c,2}=0$, and to the four grounded normal leads $N_{1}$, $N_{2}$,
$N_{3}$ and $N_{4}$. The tunneling between the multilevel quantum dot and the
superconducting leads is provided in Eqs. (II.2.1)-(15). The tunneling to the
normal leads is provided in Eq. (20). This model will be used in the following
sections III.1, III.2.1, III.2.2, III.3 and III.4.  Figure 3: The four-
terminal phenomenological model of superconducting double 0D quantum dot with
normal lead in parallel: The figure shows a double 0D quantum dot connected to
the four superconducting leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$
biased at $V_{a,b}=\pm V$ and $V_{c,1}=V_{c,2}=0$, and to the four grounded
normal leads $N_{1}$, $N_{2}$, $N_{3}$ and $N_{4}$. The corresponding
Hamiltonian is provided in section II.2.2. The coupling $\Sigma^{(0)}$
directly couples the two quantum dots $D_{x}$ and $D_{y}$, see Eq. (8). The
“direct” coupling $\Sigma^{(2)}$ connects $D_{x}$ to $S_{a}$, $D_{x}$ to
$S_{c,1}$, $D_{y}$ to $S_{b}$ and $D_{y}$ to $S_{c,2}$, see Eqs.
(II.2.2)-(25).The “crossed” coupling $\Sigma^{\prime(2)}$ connects $D_{x}$ to
$S_{b}$, $D_{x}$ to $S_{c,2}$, $D_{y}$ to $S_{a}$ and $D_{y}$ to $S_{c,1}$,
see Eqs. (II.2.2)-(30). The coupling $\Sigma^{(3)}$ connects the 0D quantum
dot $D_{x}$ to $N_{x}$ and $D_{y}$ to $N_{y}$, where $N_{x}$ and $N_{y}$ are
the two tight-binding sites on the normal lead $N$-side of the contacts, see
Eqs. (31)-(33). This model will be used in the following sections III.1,
III.2.1, III.2.2, III.2.3, III.3 and III.4.

###  II.2 Simplified model Hamiltonians

Our previous paper I paperI  was based on perturbation theory in the tunnel
amplitudes between a 2D metal and four superconducting leads, in the $V=0^{+}$
adiabatic limit. In this paper I, the poles of the Green’s functions are at
the gap edge singularities and those perturbative calculations do not yield
small voltage scales in the quartet critical current-voltage characteristics.
In our previous paper II, we addressed the connection between the quartet
critical current and the Floquet spectra in a situation where 0D quantum dot
is connected to four superconducting leads. We then obtained reduction of the
quartet critical current at the bias voltage resonance values. Compared to
this paper II, the here-considered four-terminal double quantum dot Josephson
junction is additionally connected to normal lead.

In this section, we provide the Hamiltonians of the three considered
simplified models, the third one being numerically implemented in the
forthcoming section V.

Subsection II.2.1 deals with a multilevel quantum dot connected to four
superconducting and to normal leads. Subsections II.2.2 and II.2.3 deal with
two simple models of double 0D quantum dots, also attached to superconducting
and normal leads. Subsections II.2.2 and II.2.3 correspond to “double 0D
quantum dot connected to normal lead”, in parallel or in series respectively.
The interest of double 0D quantum dots with respect to single, triple or
quadruple 0D quantum dots is explained in subsection II.2.4.

####  II.2.1 Multilevel quantum dot

Now, we consider the Hamiltonian of a multilevel quantum dot connected to
superconducting and normal leads. We also provide a physical discussion. This
defines a first stage in reducing the four-terminal device to simpler
Hamiltonians.

We assume that the multilevel quantum dot on figure 2 is intermediate between
the short- and long-junction limits, i.e. it has dimension $\gtrsim 2\xi_{0}$,
where $\xi_{0}=\hbar v_{F}/\Delta$ is the BCS coherence length, with $v_{F}$
the Fermi velocity. Thus, more than two ABS are formed at equilibrium.

The Hamiltonian: We consider finite-size discrete levels on the multilevel
quantum dot of figure 2, at the energies $\Omega_{\psi}$. The corresponding
Hamiltonian takes the form

${\cal
H}_{multi}=\sum_{\psi}\Omega_{\psi}\sum_{\sigma_{z}}c_{\psi,\sigma_{z}}^{+}c_{\psi,\sigma_{z}},$
(11)

where $c_{\psi,\sigma_{z}}^{+}$ creates a fermion with spin-$\sigma_{z}$ in
the state $|\psi\rangle$. The superconducting leads $S_{a}$, $S_{b}$,
$S_{c,1}$ and $S_{c,2}$ are described by the BCS Hamiltonian given by Eq.
(II.1.1). The normal lead is described by Eq. (II.1.3). Tunneling between the
leads $(S_{a},S_{b},S_{c,1},S_{c,2})$ biased at $(V,-V,0,0)$ and the
multilevel quantum dot is described by multi-channel contacts at the time
$\tau$:

$\displaystyle{\cal H}_{T,a}(\tau)$ $\displaystyle=$
$\displaystyle-\Sigma^{(1)}_{a}\exp\left(-ieV\tau/\hbar\right)\sum_{s}\sum_{\sigma_{z}}c_{a_{s},\sigma_{z}}^{+}c_{\alpha_{s},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,b}(\tau)$ $\displaystyle=$
$\displaystyle-\Sigma^{(1)}_{b}\exp\left(ieV\tau/\hbar\right)\sum_{t}\sum_{\sigma_{z}}c_{b_{t},\sigma_{z}}^{+}c_{\beta_{t},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,c_{1}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(1)}_{c_{1}}\sum_{u}\sum_{\sigma_{z}}c_{{c_{1,u}},\sigma_{z}}^{+}c_{\gamma_{1,u},\sigma_{z}}+h.c.$
(14) $\displaystyle{\cal H}_{T,c_{2}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(1)}_{c_{2}}\sum_{v}\sum_{\sigma_{z}}c_{{c_{2,v}},\sigma_{z}}^{+}c_{\gamma_{2,v},\sigma_{z}}+h.c.,$
(15)

where the integers $s,\,t,\,u$ and $v$ label the collection of hopping
amplitudes in real space, thus realizing a multichannel interface. As an
approximation, the hopping amplitudes in Eqs. (II.2.1)-(15) were chosen
identical within each contact: $\Sigma^{(1)}_{a,s}\equiv\Sigma^{(1)}_{a}$,
$\Sigma^{(1)}_{b,t}\equiv\Sigma^{(1)}_{b}$,
$\Sigma^{(1)}_{c_{1},u}\equiv\Sigma^{(1)}_{c_{1}}$ and
$\Sigma^{(1)}_{c_{2},v}\equiv\Sigma^{(1)}_{c_{2}}$. The contact transparencies
are parameterized by

$\displaystyle\Gamma_{a}$ $\displaystyle=$
$\displaystyle\left(\Sigma^{(1)}_{a}\right)^{2}/W,$ (16)
$\displaystyle\Gamma_{b}$ $\displaystyle=$
$\displaystyle\left(\Sigma^{(1)}_{b}\right)^{2}/W,$ (17)
$\displaystyle\Gamma_{c_{1}}$ $\displaystyle=$
$\displaystyle\left(\Sigma^{(1)}_{c_{1}}\right)^{2}/W,$ (18)
$\displaystyle\Gamma_{c_{2}}$ $\displaystyle=$
$\displaystyle\left(\Sigma^{(1)}_{c_{2}}\right)^{2}/W.$ (19)

The notations $c_{a_{s},\sigma_{z}}^{+}$, $c_{b_{t},\sigma_{z}}^{+}$,
$c_{{c_{1,u}},\sigma_{z}}^{+}$ and $c_{{c_{2,v}},\sigma_{z}}^{+}$ refer to
creating a spin-$\sigma_{z}$ fermion at the tight-binding sites labeled by
$a_{s}$, $b_{t}$, $c_{1,u}$ and $c_{2,v}$ belonging to the superconducting
leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$, see $(a_{1},a_{2},a_{3})$,
$(b_{1},b_{2})$, $(c_{1,1},c_{1,2},c_{1,3})$ and $(c_{2,1},c_{2,2})$ in figure
2. Similarly, $c_{\alpha_{s},\sigma_{z}}$, $c_{\beta_{t},\sigma_{z}}$,
$c_{\gamma_{1,u},\sigma_{z}}$ and $c_{\gamma_{2,v},\sigma_{z}}$ destroy a
spin-$\sigma_{z}$ fermion at the tight-binding sites labeled by $\alpha_{s}$,
$\beta_{t}$, $\gamma_{1,u}$ and $\gamma_{2,v}$, see
$(\alpha_{1},\alpha_{2},\alpha_{3})$, $(\beta_{1},\beta_{2})$,
$(\gamma_{1,1},\gamma_{1,2},\gamma_{1,3})$ and $(\gamma_{2,1},\gamma_{2,2})$
in figure 2.

We assume in addition that $q_{0}$ normal leads labeled by $n_{q}=1,...,q_{0}$
are connected to the multilevel quantum dot by the following tunneling
Hamiltonian:

${\cal
H}_{T,N_{q}}=-\Sigma^{(1)}_{N_{q}}\sum_{w}\sum_{\sigma_{z}}c_{n_{q,w},\sigma_{z}}^{+}c_{\nu_{q,w},\sigma_{z}}+h.c.,$
(20)

where $n_{q,w}$ is the tight-binding site labeled by $w$ in the lead $N_{q}$,
and $\nu_{q,w}$ is its counterpart on the multilevel quantum dot. We took in
Eq. (20) the same tunneling amplitudes
$\Sigma^{(1)}_{N_{q}}\equiv\Sigma^{(1)}_{N_{q},w}$ at each contact. For
instance, for $q_{0}=4$, the tight-binding sites $(n_{1,1},n_{1,2})$,
$n_{2,1}$, $n_{3,1}$ and $(n_{4,1},n_{4,2})$ belonging to the normal leads
$N_{1}$, $N_{2}$, $N_{3}$ and $N_{4}$ respectively are shown in figure 2,
together with their counterparts $(\nu_{1,1},\nu_{1,2})$, $\nu_{2,1}$,
$\nu_{3,1}$ and $(\nu_{4,1},\nu_{4,2})$ in the multilevel quantum dot.

This model will be used in the following sections III.1, III.2.1, III.2.2,
III.3 and III.4.

Physical remarks: Two antagonist effects appear:

(i) In absence of coupling to the superconductors, increasing the coupling to
the normal leads has the tendency to smoothen the energy-dependence of the
local density of states.

(ii) In the presence of coupling to the superconducting leads, the Floquet
resonances have the tendency to produce sharp peaks in the local density of
states as a function of energy.

Physically, the above item (i) captures the metallic limit of “weak sample-to-
sample fluctuations in the quartet critical current”. The item (ii)
corresponds to “strong sample-to-sample fluctuations in the quartet current”,
where samples differ by the number of channels at the contacts, or by
fluctuations in the shape of the multilevel quantum dot. The metallic limit
(i) is ruled out for compatibility with the experiment in Ref. HGE, , because
it does not produce small voltage scales. This is why we focus here on the
regime (ii) of weak coupling to the normal leads.

In most of the work, we make an additional assumption about the spectrum
$\\{\Omega_{\psi}\\}$ of the multilevel quantum dot, see Eq. (11). Namely, we
assume pairs of levels at opposite energies, i.e. there are values $\psi_{1}$,
$\psi_{2}$ of the label $\psi$ such that
$\Omega_{\psi_{1}}=-\Omega_{\psi_{2}}$. Then, two-Cooper pair resonance
emerges at specific values of the bias voltage, such as
$2eV_{*}=\Omega_{\psi_{1}}-\Omega_{\psi_{2}}=2\Omega_{\psi_{1}}$ in the limit
of weak coupling. However, within the considered models, we will find in
section V robustness of the sharp resonances with respect to detuning from the
condition of opposite energy levels.

Next, we simplify the multilevel quantum dot into phenomenological models of
double 0D quantum dots that seem to be the “minimal models” for the two-Cooper
pair resonance.

####  II.2.2 Phenomenological model of double 0D quantum dot with normal lead
in parallel

Now, we present the Hamiltonian of a double 0D quantum dot with normal lead in
parallel, and additionally connected to superconducting leads. We also provide
a physical discussion. This defines a second stage in proposing models that
can practically be numerically implemented.

The Hamiltonian: Specifically, we consider a second model, i.e. a
“phenomenological model of double 0D quantum dot with normal lead in
parallel”, see figure 3. The quantum dots $D_{x}$ and $D_{y}$ have the
vanishingly small Hamiltonians of Eqs. (9) and (10). The coupling between
$D_{x}$ and $D_{y}$ is given by Eq. (8) with the hopping amplitude
$\Sigma^{(0)}$ shown as a light blue line in figure 3. Again, the
superconducting leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$ are described
by the BCS Hamiltonian, see Eq. (II.1.1). Tunneling between the leads $S_{a}$,
$S_{b}$, $S_{c,1}$ and $S_{c,2}$ and the double 0D quantum dot $(D_{x},D_{y})$
is described by single-channel contacts:

$\displaystyle{\cal H}_{T,D_{x},a}(\tau)$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{x},a}\exp\left(-ieV\tau/\hbar\right)\sum_{\sigma_{z}}c_{a,\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,D_{y},b}(\tau)$
$\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{y},b}\exp\left(ieV\tau/\hbar\right)\sum_{\sigma_{z}}c_{b,\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,D_{x},c_{1}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{x},c_{1}}\sum_{\sigma_{z}}c_{{c_{1}},\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}+h.c.$
(23) $\displaystyle{\cal H}_{T,D_{y},c_{2}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{y},c_{2}}\sum_{\sigma_{z}}c_{{c_{2}},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}+h.c.,$
(24)

and we assume

$\Sigma^{(2)}_{D_{x},a}=\Sigma^{(2)}_{D_{y},b}=\Sigma^{(2)}_{D_{x},c_{1}}=\Sigma^{(2)}_{D_{y},c_{2}}\equiv\Sigma^{(2)}.$
(25)

We parameterize the contact transparency with
$\Gamma=\left(\Sigma^{(2)}\right)^{2}/W$.

In addition, we include the following “crossed tunneling terms” between
$D_{x}$ and $S_{b},\,S_{c,2}$ and between $D_{y}$ and $S_{a},\,S_{c,1}$:

$\displaystyle{\cal H}_{T,D_{x},b}(\tau)$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{x},b}\exp\left(ieV\tau/\hbar\right)\sum_{\sigma_{z}}c_{b,\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,D_{y},a}(\tau)$
$\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{y},a}\exp\left(-ieV\tau/\hbar\right)\sum_{\sigma_{z}}c_{a,\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}$
$\displaystyle+h.c.$ $\displaystyle{\cal H}_{T,D_{x},c_{2}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{x},c_{2}}\sum_{\sigma_{z}}c_{{c_{2}},\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}+h.c.$
(28) $\displaystyle{\cal H}_{T,D_{y},c_{1}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(2)}_{D_{y},c_{1}}\sum_{\sigma_{z}}c_{{c_{1}},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}+h.c.,$
(29)

and we assume

$\Sigma^{(2)}_{D_{x},b}=\Sigma^{(2)}_{D_{y},a}=\Sigma^{(2)}_{D_{x},c_{2}}=\Sigma^{(2)}_{D_{y},c_{1}}\equiv\Sigma^{\prime(2)}.$
(30)

We parameterize the crossed contact transparency with
$\Gamma^{\prime}=\left(\Sigma^{\prime(2)}\right)^{2}/W$.

The Hamiltonian of the normal lead is given by Eq. (II.1.3). Tunneling between
the quantum dots $D_{x},\,D_{y}$ and the normal lead $N$ is realized with the
single-channel contacts $(D_{x},N_{x})$ and $(D_{y},N_{y})$:

$\displaystyle{\cal H}_{T,D_{x},N_{x}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(3)}_{D_{x},N_{x}}\sum_{\sigma_{z}}c_{N_{x},\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}+h.c.$
(31) $\displaystyle{\cal H}_{T,D_{y},N_{y}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(3)}_{D_{y},N_{y}}\sum_{\sigma_{z}}c_{N_{y},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}+h.c.,$
(32)

$N_{x}$ and $N_{y}$ being tight-binding sites belonging to the normal lead,
and connected to $D_{x}$ and $D_{y}$ respectively by Eqs. (31)-(32). We
additionally assume that $\Sigma^{(3)}_{D_{x},N_{x}}$ and
$\Sigma^{(3)}_{D_{y},N_{y}}$ are identical:

$\Sigma^{(3)}_{D_{x},N_{x}}=\Sigma^{(3)}_{N_{x},D_{x}}=\Sigma^{(3)}_{D_{y},N_{y}}=\Sigma^{(3)}_{N_{y},D_{y}}\equiv\Sigma^{(3)}.$
(33)

This model will be used in the following sections III.1, III.2.1, III.2.2,
III.2.3, III.3 and III.4.

Physical remarks: Once disconnected from the superconducting or normal leads,
this phenomenological model of double 0D quantum dot with normal lead in
parallel naturally yields the pair of opposite-energy levels that is relevant
to the two-Cooper pair resonance.

Now, we consider in section II.2.3 another phenomenological model that also
includes the possibility of quartet current oscillating at the scale of the
Fermi wave-length as a function a dimension of the device.

Figure 4: The four-terminal phenomenological model of superconducting double
0D quantum dot with normal lead in series: The figure shows a double 0D
quantum dot connected to the four superconducting leads $S_{a}$, $S_{b}$,
$S_{c,1}$ and $S_{c,2}$ biased at $V_{a,b}=\pm V$ and $V_{c,1}=V_{c,2}=0$, and
to the four grounded normal leads $N_{1}$, $N_{2}$, $N_{3}$ and $N_{4}$. The
corresponding Hamiltonian is provided in section II.2.3. The coupling
$\Sigma^{(4)}$ [see Eqs. (34)-(36)] couples the two quantum dots $D_{x}$ and
$D_{y}$ to the conductor $N^{\prime}$ that, in the calculations, is
characterized by the local and nonlocal Green’s functions among the tight-
binding sites $N_{x}$ and $N_{y}$. The “direct” coupling $\Sigma^{(2)}$
connects $D_{x}$ to $S_{a}$, $D_{x}$ to $S_{c,1}$, $D_{y}$ to $S_{b}$ and
$D_{y}$ to $S_{c,2}$, see Eqs. (II.2.2)-(25). The “crossed” coupling
$\Sigma^{\prime(2)}$ connects $D_{x}$ to $S_{b}$, $D_{x}$ to $S_{c,2}$,
$D_{y}$ to $S_{a}$ and $D_{y}$ to $S_{c,1}$, see Eqs. (II.2.2)-(30). This
model will be used in the following sections III.2.1, III.2.2, III.2.3, III.3
and V. .

####  II.2.3 Phenomenological model of double 0D quantum dot with normal lead
in series

Now, we consider the Hamiltonian of a double 0D quantum dot with normal lead
in series, and connected to superconducting leads. We also provide a physical
discussion. The numerical calculations presented in the forthcoming section V
are based on this third stage of the simplifications.

The Hamiltonian: Namely, we consider a third model, i.e. a “phenomenological
model of double 0D quantum dot with normal lead in series” connected by
single-channel contacts through a 2D metal, see figure 4. The difference with
the previous section II.2.2 is in the coupling to the normal lead.

As in the above section II.2.2, the quantum dots $D_{x}$ and $D_{y}$ have the
vanishingly small Hamiltonians of Eqs. (9) and (10), the superconducting leads
$S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$ are described by the BCS Hamiltonian
given by Eq. (II.1.1) and the normal lead is described by Eq. (II.1.3).
Tunneling between the leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$ and the
double quantum dot $(D_{x},D_{y})$ is described by the same single-channel
contacts as in the above Eqs. (II.2.2)-(24). In addition, the crossed
tunneling amplitudes given by Eqs. (II.2.2)-(29) are also included for
completeness.

The coupling to the normal lead is phenomenologically accounted for by
replacing the tunneling amplitude $\Sigma^{(0)}$ between the quantum dots by
propagation through the normal conductor $N^{\prime}$ connected to the dots
$D_{x}$ and $D_{y}$ by the single-channel tunneling amplitudes
$\Sigma^{(0)}_{D_{x},N^{\prime}_{x}}$ and
$\Sigma^{(0)}_{D_{y},N^{\prime}_{y}}$ respectively:

$\displaystyle{\cal H}_{T,D_{x},N^{\prime}_{x}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}\sum_{\sigma_{z}}c_{N^{\prime}_{x},\sigma_{z}}^{+}c_{D_{x},\sigma_{z}}+h.c.$
(34) $\displaystyle{\cal H}_{T,D_{y},N^{\prime}_{y}}$ $\displaystyle=$
$\displaystyle-\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}\sum_{\sigma_{z}}c_{N^{\prime}_{y},\sigma_{z}}^{+}c_{D_{y},\sigma_{z}}+h.c.,$
(35)

where we assume

$\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}=\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}=\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}=\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}\equiv\Sigma^{(4)}.$
(36)

The normal conductor $N^{\prime}$ is characterized by the local Green’s
functions at $N^{\prime}_{x}$ or $N^{\prime}_{y}$ , in addition to nonlocal
ones from $N^{\prime}_{x}$ to $N^{\prime}_{y}$ or from $N^{\prime}_{y}$ to
$N^{\prime}_{x}$ across $N^{\prime}$. The nonlocal Green’s function between
$N^{\prime}_{x}$ and $N^{\prime}_{y}$ generally has to be a complex number.

This model will be used in the following sections III.2.1, III.2.2, III.2.3,
III.3 and V.

Physical remarks: Within this phenomenological model of double 0D quantum dot
with normal lead in series, the quartet critical current at resonance is
sensitive to the Green’s function connecting the two quantum dots $D_{x}$ and
$D_{y}$. This Green’s function is itself sensitive to the microscopic details
of the model, i.e. on the value of the separation $R_{0}$ between the tight-
binding sites $N^{\prime}_{x}$ and $N^{\prime}_{y}$ within the interval
$\left[R^{\prime}_{0}-\lambda_{F}/2,R^{\prime}_{0}+\lambda_{F}/2\right]$,
where $\lambda_{F}$ is the Fermi wave-length. Thus, this phenomenological
model of double 0D quantum dot with normal lead in series produces strong
sample-to-sample fluctuations of the quartet current, where different samples
correspond to different values of $R_{0}$.

Once disconnected from the superconducting or normal leads, this
phenomenological model of double 0D quantum dot with normal lead in series
also yields the pair of energy levels that is central to the two-Cooper pair
resonance.

In section V, we present numerical results for this model in the regime of
strong sample-to-sample fluctuations.

Figure 5: The nonlocal resonances in three-terminal devices. Panels a, d, g
and j show the lowest-order diagrams of elastic cotunneling (EC), crossed
Andreev reflection (CAR), double elastic cotunneling (dEC) and double crossed
Andreev reflection (dCAR) respectively. Panels b, e, h and k show the
energy-$E$-dependence of the density of states $\rho_{x}(E)$ on the quantum
dot $D_{x}$ and panels c, f, i and l show the corresponding $\rho_{y}(E)$ on
the quantum dot $D_{y}$. Panels e, f, k and l also show for completeness by
green dashed lines the complementary resonances that do not directly
contribute to the corresponding transport processes.

####  II.2.4 Further physical remarks on the choice of a model

Now, we explain why we focus on the double quantum dot Hamiltonians presented
in the above subsections II.2.2 and II.2.3, instead of a single, three or
four-quantum dot Hamiltonians.

First, we comparatively examine single and double 0D quantum dots.

At equilibrium, the ABS of a single 0D quantum dot have energy scale $\Gamma$
typically set by $\Gamma=\Sigma^{2}/W$, where $\Sigma$ is the tunneling
amplitude between the dot and the superconducting leads. It turns out that, at
small $\Gamma$, the quasiadiabatic regime $V\lesssim V_{*}$ with
$eV_{*}\sim\Gamma$, does not produce enhanced quartet critical current
according to the forthcoming subsection III.3 where the two-Cooper pair
resonance is put in correspondence with nonadiabatic effects.

At equilibrium, double 0D quantum dots are characterized by levels at the
opposite energies $\pm\Omega$, where $\Omega$ is set by the tunneling
amplitude $\Sigma^{(0)}$ between the dots, see Eq. (8). Then, the voltage
energy $eV_{*}=\Omega$ of the first resonance remains finite in the weak-
coupling limit, which implies nonadiabatic behavior and enhancement of the
quartet critical current at low $V$ if $\Omega$ is small compared to the
superconducting gap $\Delta$.

Double 0D quantum dots are thus better candidates than single 0D quantum dots
to produce large quartet critical current at small bias voltage because the
two-Cooper pair resonance is also there in the weak-coupling regime of
small-$\Gamma$, see also the forthcoming section VI.

The interest of double 0D quantum dots is also related to the observation that
the corresponding four energy levels can accommodate the four fermions of a
quartet as a “real state”. In this sense, the double 0D quantum dots have the
“minimal complexity” in comparison with three- or four-quantum dot devices.

##  III Mechanism for the two-Cooper pair resonance

In this section, we present a general mechanism for emergence of the two-
Cooper pair resonance. We start in subsection III.1 with introducing the
three-terminal $N_{a}$-dot-$S_{c}$-dot-$N_{b}$ Cooper pair beam splitters, and
the $S_{a}$-dot-$S_{c}$-dot-$S_{b}$ three-terminal Josephson junctions. Next,
we calculate effective non-Hermitian self-energies and Hamiltonians in
subsection III.2. The transport formula are discussed in subsection III.3 in
connection with emergence of resonances. Current conservation is discussed in
subsection III.4.

###  III.1 Three-terminal two-Cooper pair resonance

In this subsection, we introduce the two-Cooper pair resonance in three-
terminal Josephson junctions, i.e. in $S_{a}$-dot-$S_{c}$-dot-$S_{b}$ devices,
starting with nonlocal resonances in the three-terminal
$N_{a}$-dot-$S_{c}$-dot-$N_{b}$.

Figures 5a, 5b and 5c show the condition for nonlocal resonance of elastic
cotunneling (EC) Yeyati2007 ; theory-CPBS1 ; theory-CPBS2 ; theory-CPBS3 ;
theory-CPBS4 ; theory-CPBS4-bis ; theory-CPBS4-ter ; theory-CPBS5 ; theory-
CPBS6 ; theory-CPBS7 ; theory-CPBS8 ; theory-CPBS9 ; theory-CPBS10 ; theory-
CPBS11  in $N_{a}$-dot-$S_{c}$-dot-$N_{b}$ Cooper pair beam splitters. EC
transfers single-particle states from $N_{a}$ to $N_{b}$ across $S_{c}$, and
contributes to negative nonlocal conductance ${\cal G}_{a,b}=\partial
I_{a}/\partial V_{b}$ on the voltage biasing condition $V_{a}=-V_{b}$, where
$I_{a}$ is the current transmitted into the normal lead $N_{a}$. EC is
resonant if both quantum dots $D_{x}$ and $D_{y}$ have levels at the same
energy $\epsilon_{x}=\epsilon_{y}$, where $\omega=\epsilon_{x}=\epsilon_{y}$
is the energy of the incoming electron in figure 5a.

Figures 5d, 5e, 5f show CAR Yeyati2007 ; theory-CPBS1 ; theory-CPBS2 ; theory-
CPBS3 ; theory-CPBS4 ; theory-CPBS4-bis ; theory-CPBS4-ter ; theory-CPBS5 ;
theory-CPBS6 ; theory-CPBS7 ; theory-CPBS8 ; theory-CPBS9 ; theory-CPBS10 ;
theory-CPBS11  in $N_{a}$-dot-$S_{c}$-dot-$N_{b}$ Cooper pair beam splitters,
which reflects by nonlocal Andreev reflection spin-up electron impinging from
$N_{a}$ as spin-down hole transmitted into $N_{b}$, leaving a Cooper pair in
the central $S_{c}$. CAR contributes for positive value to the nonlocal
conductance ${\cal G}_{a,b}$ at the bias voltages $V_{a}=V_{b}$ and $V_{c}=0$.
The nonlocal CAR resonance is obtained if the quantum dots $D_{x}$ and $D_{y}$
have levels at the opposite energies $\epsilon_{x}=-\epsilon_{y}\equiv\omega$.

Considering $S_{a}$-dot-$S_{c}$-dot-$S_{b}$ three-terminal Josephson
junctions, figures 5g, 5h, 5i feature double elastic cotunneling (dEC) theory-
CPBS8 ; theory-CPBS9 ; theory-CPBS10 , which transfers Cooper pairs from
$S_{a}$ to $S_{b}$ across $S_{c}$ at the bias voltages $V_{a}=V_{b}$. This
process is resonant if $\omega=E_{0}$ (for the spin-up electron crossing
$D_{x}$ or $D_{y}$) and $-\omega+2eV=-E_{1}$ (for the spin-down hole crossing
$D_{x}$ or $D_{y}$). Thus, dEC is resonant if $eV=(E_{0}-E_{1})/2$ and
$\omega=E_{0}$.

Concerning double crossed Andreev reflection (dCAR) in
$S_{a}$-dot-$S_{c}$-dot-$S_{b}$ three-terminal Josephson junctions in figures
5j, 5k and 5l, two Cooper pairs from $S_{a}$ and $S_{b}$ biased at the
voltages $\pm V$ cooperatively enter the grounded $S_{c}$, producing transient
correlations among four fermions, i.e. the so-called quartets. Then,

$\displaystyle\omega$ $\displaystyle=$ $\displaystyle E_{0}$ (37)
$\displaystyle-\omega+2eV$ $\displaystyle=$ $\displaystyle-E_{1}$ (38)

are obtained at resonance. Conversely, the same Eqs. (37)-(38) are obtained
for resonance of the spin-up electron and spin-down holes crossing $D_{y}$.

Now, we assume that $D_{x}$ and $D_{y}$ are gathered into a single multilevel
quantum dot. Assuming the opposite energies $E_{0}=-E_{1}\equiv\Omega$ implies
$\omega=\Omega$ and $-\omega+2eV=\Omega$, which yields

$eV=\omega=\Omega$ (39)

for the nonlocal quartet resonance.

Overall, the argument leading to Eq. (39) confirms that nonlocal quartet
resonance is produced at voltage energy $eV$ that can be much smaller than the
superconducting gap $\Delta$, if the energy scales $\pm\Omega$ are also within
$\pm\Delta$, see also the preceding subsection II.2.4.

###  III.2 Effective non-Hermitian self-energy and Hamiltonian

In this subsection, we present how effective non-Hermitian self-energies or
Hamiltonians can be obtained in the infinite-gap limit from the models
presented in the above section II.

We start in section III.2 with the Dyson equations for the three models
presented in the above subsections II.2.1, II.2.2 and II.2.3. Emergence of
non-Hermitian self-energy in the infinite-gap limit is discussed in subsection
III.2.2. Non-Hermitian Hamiltonians are next presented for double quantum dots
in the infinite-gap limit, see subsection III.2.3.

####  III.2.1 Closing the Dyson equations

In this subsection, we provide the starting-point Dyson equations. We adopt a
general viewpoint on the three models of subsections II.2.1, II.2.2 and
II.2.3, i.e. multilevel quantum dots, double 0D quantum dots connected in
parallel and in series respectively, see also figures 2, 3 and 4.

The “quantum dot $D$” is coupled to $p_{0}$ superconducting leads labeled by
$n_{p}=1,\,...,\,p_{0}$ and to $q_{0}$ normal leads labeled by
$n_{q}=1\,,...,\,q_{0}$, where $p_{0}=q_{0}=4$ are used in figure 2. The Dyson
equations take the form

$\displaystyle\hat{G}_{D,D}$ $\displaystyle=$
$\displaystyle\hat{g}_{D,D}+\sum_{n_{q}=1}^{q_{0}}\sum_{n_{q^{\prime}}=1}^{q_{0}}\hat{g}_{D,D}\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}\hat{\Sigma}^{(5)}_{N_{n_{q^{\prime}}},D}\hat{G}_{D,D}$
(40)
$\displaystyle+\sum_{n_{p}=1}^{p_{0}}\hat{g}_{D,D}\hat{\Sigma}^{(5)}_{D,S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{S_{n_{p}},D}\hat{G}_{D,D}.$

The notation “$D$” in Eq. (40) refers to the collection of the tight-binding
sites in absence of coupling to the superconducting or normal leads, and to
the direct coupling between them. For instance, $D$ represents finite-size
tight-binding multilevel quantum dot, see subsection II.2.1, or
$D_{x}$-$D_{y}$ double quantum dot, see subsections II.2.2 and II.2.3. The
notation $\Sigma^{(5)}$ denotes generic hopping amplitude, with the normal
lead $N_{n_{q}}$ (corresponding to
$\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}=\hat{\Sigma}^{(5)}_{N_{n_{q}},D}$), or with
the superconducting lead $S_{n_{p}}$ (corresponding to
$\hat{\Sigma}^{(5)}_{D,S_{n_{p}}}=\hat{\Sigma}^{(5)}_{S_{n_{p}},D}$).

The notations $\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}$ and
$\hat{g}_{S_{n_{p}},S_{n_{p}}}$ stand for the bare Green’s functions of the
normal or superconducting leads $N_{n_{q}}$, $N_{n_{q^{\prime}}}$ or
$S_{n_{p}}$ respectively, and $\hat{g}_{D,D}$, $\hat{G}_{D,D}$ are the bare
and fully dressed quantum dot Green’s functions. The bare Green’s functions
are given by

$\displaystyle\hat{g}_{D,D,{\bf x},{\bf y}}^{A,1,1}$ $\displaystyle=$
$\displaystyle\sum_{\psi}\langle{\bf
x}|\psi\rangle\frac{1}{\omega-\epsilon_{\psi}-i\eta}\langle\psi|{\bf
y}\rangle$ (41) $\displaystyle\hat{g}_{D,D,{\bf x},{\bf y}}^{A,2,2}$
$\displaystyle=$ $\displaystyle\sum_{\psi}\langle{\bf
x}|\psi\rangle\frac{1}{\omega+\epsilon_{\psi}-i\eta}\langle\psi|{\bf
y}\rangle,$ (42)

where $\omega$ is the real-part of the energy and “$1,1$”,“$2,2$” refer to the
electron-electron and hole-hole channels respectively. In addition,
$\hat{g}_{D,D,{\bf x},{\bf y}}^{A,1,2}=\hat{g}_{D,D,{\bf x},{\bf
y}}^{A,2,1}=0$ because of the absence of superconducting pairing on the
quantum dot. The Dynes parameter $\eta$ in Eqs. (41)-(42) Kaplan ; Dynes ;
Pekola1 ; Pekola2  makes the distinction between the advanced and retarded
Green’s functions, and it can be used to capture relaxation on the quantum dot
FWS ; Bergeret . As in the previous Eq. (11), the notation $|\psi\rangle$ in
Eqs. (41)-(42) stands for the single-particle states, and $\langle{\bf
x}|\psi\rangle$, $\langle{\bf y}|\psi\rangle$ are the corresponding wave-
functions at the tight-binding sites ${\bf x}$ and ${\bf y}$.

The bare and fully dressed Green’s functions $\hat{g}_{D,D}$ and
$\hat{G}_{D,D}$ have entries in the Nambu labels, and in the tight-binding
sites making the contacts between the dot and the normal or superconducting
leads $N_{n_{q}}$ or $S_{n_{p}}$. In addition, those matrices have entries in
the set of the harmonics of the Josephson frequency.

####  III.2.2 Effective non-Hermitian self-energy

In this subsection, we explain how effective non-Hermitian self-energy is
obtained in the infinite-gap limit. It was already mentioned in the
introductory section that the infinite-gap limit was introduced and considered
over the last few years, see for instance Refs. topo1_plus_Floquet, ; Zazunov,
; Meng,  to cite but a few.

Specifically, Eq. (40) is rewritten as

$\displaystyle\left[\hat{g}_{D,D}^{-1}(\omega)-\sum_{n_{q}=1}^{q_{0}}\sum_{n_{q^{\prime}}=1}^{q_{0}}\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}\hat{\Sigma}^{(5)}_{N_{n_{q^{\prime}}},D}\right.$
(43)
$\displaystyle-\left.\sum_{n_{p}=1}^{p_{0}}\hat{\Sigma}^{(5)}_{D,S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{S_{n_{p}},D}\right]\hat{G}_{D,D}(\omega)=\hat{I}.$

In addition, $\hat{g}_{S,S}$ is independent on $\omega$ in the considered
infinite-gap limit, and the local Green’s function
$\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}$ in the normal leads is also taken as
being independent on energy.

Then, Eq. (43) can be rewritten as

$\displaystyle\left[\omega-\hat{\Sigma}_{Non\,Her.}(\omega)\right]^{-1}\hat{G}_{D,D}(\omega)=\hat{I},$
(44)

where the self-energy

$\displaystyle\hat{\Sigma}_{Non\,Her.}(\omega)$ $\displaystyle=$
$\displaystyle\omega-\hat{g}_{D,D}^{-1}(\omega)$
$\displaystyle+\sum_{n_{q}=1}^{q_{0}}\sum_{n_{q^{\prime}}=1}^{q_{0}}\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}\hat{\Sigma}^{(5)}_{N_{n_{q^{\prime}}},D}$
$\displaystyle+\sum_{n_{p}=1}^{p_{0}}\hat{\Sigma}^{(5)}_{D,S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{S_{n_{p}},D}$

in non-Hermitian, due to $\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}$.

####  III.2.3 Specializing to double quantum dots

In this subsection, we discuss emergence of non-Hermitian Hamiltonian for the
phenomenological model of double 0D quantum dots. The normal leads in
parallel, see subsection II.2.2, are treated in the same framework as normal
leads in series, see subsection II.2.3.

The retarded and advanced Green’s functions $\sim 1/(\omega-{\cal
H}_{eff}^{R})$ and $\sim 1/(\omega-{\cal H}_{eff}^{A})$ are obtained from the
effective non-Hermitian Hamiltonians ${\cal H}_{eff}^{R}\equiv{\cal H}_{eff}$
and ${\cal H}_{eff}^{A}=\left({\cal H}_{eff}\right)^{+}$ respectively. Their
complex eigenvalues receive interpretation of the real and imaginary parts of
the ABS energies, yielding coherent oscillations and damping in the ABS
dynamics, respectively.

We make use of the same compact notations as in the previous subsections
III.2.1 and III.2.2. Then we show that the effective non-Hermitian self-energy
$\hat{\Sigma}_{Non\,Her.}(\omega)$ in Eq. (III.2.2) becomes energy-independent
for those double 0D quantum dots, thus defining the effective non-Hermitian
Hamiltonian $\hat{\cal H}_{eff}^{(\infty)}\equiv\hat{\Sigma}_{Non\,Her.}$ in
the infinite-gap limit.

The matrices in Eq. (43) have entries in the tight-binding sites at the
boundary of the quantum dot, where the tunneling self-energy between the dot
and the superconducting or normal leads is acting. The boundary is identical
to the bulk for the double 0D quantum dots $D_{x}$-$D_{y}$ isolated from the
leads. Thus, $g_{D,D}$ takes the value

$\hat{g}_{D,D}^{R}=\left(\omega+i\eta-\hat{\cal H}_{dot}\right)^{-1},$ (46)

where $\hat{\cal H}_{dot}$ is the double 0D quantum dot Hamiltonian. The
retarded Green’s functions of the double 0D quantum dots with normal leads in
parallel, see figure 3, can then be expressed as

$\hat{G}_{D,D}^{R}=\left(\omega+i\eta-\hat{\cal
H}_{eff}^{(\infty)}\right)^{-1},$ (47)

where the infinite-gap-limit effective Hamiltonian

$\displaystyle\hat{\cal H}_{eff}^{(\infty)}$ $\displaystyle=$
$\displaystyle\hat{\cal
H}_{dot}+\sum_{n_{q}=1}^{q_{0}}\sum_{n_{q^{\prime}}=1}^{q_{0}}\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}\hat{\Sigma}^{(5)}_{N_{n_{q^{\prime}}},D}$
$\displaystyle+\sum_{n_{p}=1}^{p_{0}}\hat{\Sigma}^{(5)}_{D,S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{D_{n_{p}},D}$

is non-Hermitian, due to $\hat{g}_{N_{n_{q}},N_{n_{q^{\prime}}}}$.

We note that infinite-gap Hamiltonians can more generally be deduced from the
self-energy $\hat{\Sigma}_{Non\,Her.}(\omega)$ in Eq. (III.2.2) by extending
the domain of definition of the tunneling amplitude to all tight-binding sites
of the quantum dots, and associating vanishingly small tunneling amplitudes to
those fictitious superconducting leads.

###  III.3 Emergence of resonances in the transport formula

In this subsection, we demonstrate that the current at resonance is inverse-
proportional to the damping rate.

Using the same notations as in the previous subsections III.2.1, III.2.2 and
III.2.3, the current flowing from the multilevel quantum dot into the
superconducting lead $S_{n_{p}}$ is given by the following energy-integral
AVERIN ; Cuevas ; Caroli ; Yeyati-noise-equilibrium :

$\displaystyle I_{D,S_{n_{p}}}$ $\displaystyle=$
$\displaystyle\frac{e}{\hbar}\int d\omega\times$
$\displaystyle\mbox{Tr}\left\\{\hat{\sigma}_{N}^{z}\left[\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{G}^{+,-}_{S_{n_{p}},D}-\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\hat{G}^{+,-}_{D,S_{n_{p}}}\right]\right\\},$

where the Pauli matrix $\hat{\sigma}_{N}^{z}$ is defined as
$\hat{\sigma}_{N}^{z}=\mbox{diag}(1,-1)$ on the set of the Nambu labels, which
is referred to as the subscript “$N$” for “Nambu”. In Eq. (III.3),
$\hat{G}^{+,-}$ is the fully dressed Keldysh Green’s function, see Refs.
Caroli, ; Cuevas, ; Yeyati-noise-equilibrium, . The bare superconducting
Green’s functions $\hat{g}^{+,-}_{S_{n_{p}},S_{n_{p}}}=0$ are vanishingly
small if the infinite-gap limit is taken, and
$\hat{g}^{+,-}_{N_{n_{q}},N_{n_{q}^{\prime}}}\neq 0$ is nonvanishingly small
in the normal leads, due the corresponding finite value of the normal density
of states.

Figure 6: Panel a schematically shows a quantum dot supporting levels at the
opposite energies $E=\pm\Omega$, coupled to the superconducting leads $S_{a}$
and $S_{b}$ biased at $eV_{a,b}=\pm\Omega$ and to the grounded $S_{c,1}$ and
$S_{c,2}$. Panels b-g represent two levels at the opposite energies
$\pm\Omega$, together with the corresponding wave-functions, with absence of
node at $E=-\Omega$, and a node at $E=\Omega$. Panels b-g show the time-
evolution of the quartet process, involving Cooper pair transmitted into
$S_{c,2}$ (panel c), Cooper pair taken from $S_{b}$ (panel d), Cooper
transmitted into $S_{c,1}$ (panel f) and Cooper pair taken from $S_{a}$ (panel
g). Detailed comments about panels b-g are given in the text.

The two terms $\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{G}^{+,-}_{S_{n_{p}},D}$
and $\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\hat{G}^{+,-}_{D,S_{n_{p}}}$ in Eq.
(III.3) can be expanded as follows:

$\displaystyle\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{G}^{+,-}_{S_{n_{p}},D}=$
(50)
$\displaystyle\sum_{q,q^{\prime}}\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\left[\hat{I}+\hat{G}^{R}\hat{\Sigma}^{(5)}\right]_{S_{n_{p}},N_{q}}\hat{g}^{+,-}_{N_{q},N_{q^{\prime}}}\left[\hat{I}+\hat{\Sigma}^{(5)}\hat{G}^{A}\right]_{N_{q^{\prime}},D}$
$\displaystyle\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\hat{G}^{+,-}_{D,S_{n_{p}}}=$
(51)
$\displaystyle\sum_{q,q^{\prime}}\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\left[\hat{I}+\hat{G}^{R}\hat{\Sigma}^{(5)}\right]_{D,N_{q}}\hat{g}^{+,-}_{N_{q},N_{q^{\prime}}}\left[\hat{I}+\hat{\Sigma}^{(5)}\hat{G}^{A}\right]_{N_{q^{\prime}},S_{n_{p}}},$

where $N_{q}$ and $N_{q^{\prime}}$ run over the interfaces between the quantum
dot and the normal lead.

Then, Eqs. (III.3)-(51) become

$\displaystyle I_{D,S_{n_{p}}}$ $\displaystyle=$
$\displaystyle\frac{e}{\hbar}\int d\omega\times$
$\displaystyle\mbox{Tr}\left\\{\hat{\sigma}_{N}^{z}\left[\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{g}^{R}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\hat{G}^{R}_{D,D}\hat{\Gamma}_{D,D}^{+,-,N}\hat{G}^{A}_{D,D}\right.\right.$
$\displaystyle-\left.\left.\hat{\Sigma}_{S_{n_{p}},D}^{(5)}\hat{G}^{R}_{D,D}\hat{\Gamma}_{D,D}^{+,-,N}\hat{G}^{A}_{D,D}\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{g}^{R}_{S_{n_{p}},S_{n_{p}}}\right]\right\\},$

where

$\hat{\Gamma}_{D,D}^{+,-,N}=\sum_{n_{q}=1}^{q_{0}}\sum_{n_{q^{\prime}}=1}^{q_{0}}\hat{\Sigma}^{(5)}_{D,N_{n_{q}}}\hat{g}^{+,-}_{N_{n_{q}},N_{n_{q^{\prime}}}}\hat{\Sigma}^{(5)}_{N_{n_{q^{\prime}}},D}$
(53)

is diagonal in Nambu and in Floquet.

The retarded and advanced Green’s functions are approximated as engineering ;
Berry

$\displaystyle\hat{G}^{R}_{D,D}$ $\displaystyle\simeq$
$\displaystyle\sum_{k^{\prime},l^{\prime}}\frac{\hat{R}_{k^{\prime},l^{\prime}}^{R}}{\omega-k^{\prime}eV-
E_{l^{\prime}}+i\delta_{l^{\prime}}}$ (54) $\displaystyle\hat{G}^{A}_{D,D}$
$\displaystyle\simeq$
$\displaystyle\sum_{k^{\prime\prime},l^{\prime\prime}}\frac{\hat{R}_{k^{\prime\prime},l^{\prime\prime}}^{A}}{\omega-k^{\prime\prime}eV-
E_{l^{\prime\prime}}-i\delta_{l^{\prime\prime}}},$ (55)

where $k^{\prime}$, $l^{\prime}$, $k^{\prime\prime}$ and $l^{\prime\prime}$
are four integers, $E_{l^{\prime}}$, $E_{l^{\prime\prime}}$ and
$\delta_{l^{\prime}}$, $\delta_{l^{\prime\prime}}$ are the Floquet energies
and line-width broadening, and $\hat{R}_{k^{\prime},l^{\prime}}^{R}$,
$\hat{R}_{k^{\prime\prime},l^{\prime\prime}}^{A}$ are the matrix residues.
Then, inserting Eqs. (54)-(55) into Eq. (III.3) and integrating over the
energy $\omega$ yields

$\displaystyle I_{D,S_{n_{p}}}$ $\displaystyle=$ $\displaystyle\frac{\pi
e}{2\hbar}\sum_{(k,l)}\frac{1}{\delta_{l}}\mbox{Tr}\left\\{\hat{\sigma}_{N}^{z}\times\right.$
$\displaystyle\left.\left[\hat{\Sigma}_{D,S_{n_{p}}}^{(5)}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}_{S_{n_{p}},D}^{(5)},\hat{R}_{k,l}^{R}\hat{\Gamma}_{D,D}^{+,-,N}\hat{R}_{k,l}^{A}\right]_{-}\right\\},$

where $\left[...,...\right]_{-}$ is a commutator and $\sum_{(k,l)}$ denotes
summation over the pairs of labels $k$ and $l$ such that

$k^{\prime}eV+E_{l^{\prime}}=k^{\prime\prime}eV+E_{l^{\prime\prime}}\equiv
keV+E_{l}.$ (57)

Thus, we find that the current is inverse proportional to the damping rate set
by the parameter $\delta_{l}$.

###  III.4 Current conservation

Now, we consider that the effective non-Hermitian Hamiltonian with complex
eigenvalues originates from the Hermitian Hamiltonians presented in the above
section II, see also Appendix A. It turns out that the total current is
conserved once the fraction of the current transmitted into the normal lead
has been taken into account. Current conservation can be demonstrated by
assuming that the average number of fermions $\langle\hat{N}_{dot}\rangle$ on
the quantum dot is stationary. The Hamiltonian ${\cal H}$ can be written as a
sum of the Hamiltonians of the quantum dot, of the lead, and tunneling between
them. Then, $d\langle\hat{N}_{dot}\rangle/dt=0$ is equivalent to current
conservation. Thus, non-Hermitian effective Hamiltonian does not contradict
current conservation. We note that self-consistent algorithms were used in
Ref. FWS,  to impose current conservation in three-terminal Josephson
junctions in the presence of a phenomenological Dynes parameter $\eta$.

##  IV Mechanism for the inversion

We discussed in the previous section III how sharp resonance peaks can appear
in the voltage-dependence of the quartet critical current. Now, we provide
simple arguments for the magnetic flux-$\Phi$ sensitivity of the quartet
critical current in the $V=0^{+}$ adiabatic limit, focusing on how the quartet
critical current $I_{q,c}(\Phi=0)$ in zero field $\Phi=0$ compares to
$I_{q,c}(\Phi=\pi)$ at half-flux quantum $\Phi=\pi$. Inversion corresponds to
larger quartet critical current at $\Phi=\pi$ than at $\Phi=0$, i.e.
$I_{q,c}(\Phi=\pi)>I_{q,c}(\Phi=0)$. Absence of inversion corresponds to
$I_{q,c}(\Phi=\pi)<I_{q,c}(\Phi=0)$.

Figure 6 shows a sequence of microscopic processes for the quartets in the
presence of the two energy levels $E=\pm\Omega$ on the quantum dot. Figure 6b
shows spin-up electron at energy $E=\Omega$ on the left-part of the two-level
quantum dot, and how it moves to the right-part. Figure 6c shows spin-up
electron at energy $E=\Omega$ on the right part of the junction, and how it is
converted by Andreev reflection into spin-down hole at energy $E=-\Omega$ at
the interface with the $S_{c,2}$ superconducting lead. Figure 6d shows the
resulting spin-down hole at energy $E=-\Omega$, together with a Cooper pair
taken from $S_{b}$. Figure 6e shows the resulting spin-up electron at energy
$E=-\Omega$ on the right-part of the quantum dot, and how it moves to the left
part. Figure 6f shows spin-up electron at energy $E=-\Omega$ on the left part
of the junction, and how it is converted into spin-down hole at energy
$E=\Omega$ by Andreev reflection at the interface with $S_{c,1}$. Figure 6g
shows the resulting spin-down hole on the left-part of the junction at the
energy $E=\Omega$, together with absorption of a Cooper pair taken from
$S_{a}$, thus coming back to the initial state in figure 6b.

Figure 7: The $eV/\Delta$-dependence of the quartet critical current. Each
panel of this figure shows $I_{q,c}(eV/\Delta,\Phi=0)$ and
$I_{q,c}(\Phi=\pi,eV/\Delta)$ as a function of the normalized bias voltage
$eV/\Delta$, for $A_{I}=0$ and $A_{R}=0.3/W$ (a), $A_{R}=0.2/W$ (b) and
$A_{R}=0.1/W$ (c). Panels d, e and f correspond to $A_{I}=0.05/W$ and
$A_{R}=0.3/W$ (d), $A_{R}=0.2/W$ (e) and $A_{R}=0.1/W$ (f) and panels g, h, i
show $I_{q,c}(\Phi=0,eV/\Delta)$ and $I_{q,c}(\Phi=\pi,eV/\Delta)$ for
$A_{I}=0.1/W$ and $A_{R}=0.3/W$ (g), $A_{R}=0.2/W$ (h) and $A_{R}=0.1/W$ (i).
We use $\Gamma/\Delta=1$, $\Gamma^{\prime}/\Delta=0$, $B_{R}=B_{I}=0$, and
$\epsilon_{x}=\epsilon_{y}=0$. The notation $W$ is used for the band-width.
Figure 8: The $eV/\Delta$-dependence of the quartet critical current. The
figure shows the reduced voltage-$eV/\Delta$-dependence of the logarithm of
the quartet critical currents
$\log_{10}\left(I_{q,c}(\Phi=0,eV/\Delta)\right)$ and
$\log_{10}\left(I_{q,c}(\Phi=\pi,eV/\Delta)\right)$ at $\Phi=0$ and $\Phi=\pi$
respectively. The panels show increasing values of
$\lambda=\Gamma^{\prime}/\Gamma$, where $\Gamma$ is the “direct” coupling
between $D_{x}$ and $S_{a}$, $S_{c,1}$ and between $D_{y}$ and $S_{b}$,
$S_{c,2}$. The crossed coupling $\Gamma^{\prime}$ is in between $D_{x}$ and
$S_{b}$, $S_{c,2}$ and between $D_{y}$ and $S_{a}$, $S_{c,1}$. We used
$\Gamma/\Delta=1$, $A_{R}=0.3/W$, $A_{I}=0.05/W$, and $B_{R}=B_{I}=0$, where
$W$ is the band-width. The on-site energies $\epsilon_{x}=\epsilon_{y}=0$ are
vanishingly small.  Figure 9: The sign of the inversion. The figure shows the
sign of
$I_{q,c}(\lambda,\Phi=\pi,eV/\Delta)-I_{q,c}(\lambda,\Phi=0,eV/\Delta)$, where
$\lambda=\Gamma^{\prime}/\Gamma$. Inversion corresponds to
$I_{q,c}(\lambda,\Phi=\pi,eV/\Delta)-I_{q,c}(\lambda,\Phi=0,eV/\Delta)>0$. The
parameter $\Gamma$ is the “direct” coupling between $D_{x}$ and
$(S_{a},S_{c,1})$ or between $D_{y}$ and $(S_{b},S_{c,2})$. The parameter
$\Gamma^{\prime}$ is in between $D_{x}$ and $(S_{b},S_{c,2})$ or $D_{y}$ and
$(S_{a},S_{c,1})$. We used $\Gamma/\Delta=1$, $A_{R}=0.3/W$, $A_{I}=0.05/W$,
and $B_{R}=B_{I}=0$, where $W$ is the band-width. The on-site energies
$\epsilon_{x}=\epsilon_{y}=0$ are vanishingly small.

Figure 10: The sign of the inversion. The figure is similar to the previous
figure 9, but now the sign of
$I_{q,c}(\epsilon_{0}/\Delta,\Phi=\pi,eV/\Delta)-I_{q,c}(\epsilon_{0},\Phi=0,eV/\Delta)$
is shown in the $(\epsilon_{0}/\Delta,eV/\Delta)$ plane. We used
$\Gamma/\Delta=1$, $A_{R}=0.3/W$, $A_{I}=0.05/W$, $B_{R}=B_{I}=0$ (where $W$
is the band-width) and $\lambda=\Gamma^{\prime}/\Gamma=1/2$. The on-site
energies are identical for both quantum dots:
$\epsilon_{x}=\epsilon_{y}=\epsilon_{0}$.

Overall the transition between figures 6b and 6c implies negative sign for the
“split quartets” transmitting a Cooper pair into $S_{c,1}$ and another one
into $S_{c,2}$, and a positive sign for the “unsplit quartets”, corresponding
to two Cooper pairs transmitted into $S_{c,1}$ or into $S_{c,2}$. The
resulting opposite signs in the “split” and “unsplit” channels imply the
inversion, see our previous paper I paperI . This physical picture for
“emergence of inversion in the $V=0^{+}$ limit” is confirmed by the
calculations presented in Appendix B.2.

Thus, we demonstrated that “inversion between $\Phi=0$ and $\Phi=\pi$” or
“absence of inversion” is linked to the quantum dot single-particle wave-
functions, i.e. to the number of nodes in their wave-functions.

##  V Numerical results

In this section, we present a selection of the numerical results for the
phenomenological model of double 0D quantum dot with normal lead in series,
see figure 4 and the Hamiltonian in subsection II.2.3. The quantum dots
$D_{x}$ and $D_{y}$ are coupled to each other by the tunneling amplitude
$\Sigma^{(0)}$ defined in Eq. (8). The coupling
$\Gamma=\left(\Sigma^{(2)}\right)^{2}/W$ is defined for hopping between
$D_{x}$ and $S_{a}$, $D_{x}$ and $S_{c,1}$, $D_{y}$ and $S_{b}$, $D_{y}$ and
$S_{c,2}$, see $\Sigma^{(2)}$ in Eq. (25). We add the coupling
$\Gamma^{\prime}=\left(\Sigma^{\prime(2)}\right)^{2}/W$ between $D_{x}$ and
$S_{b}$, $D_{x}$ and $S_{c,2}$, $D_{y}$ and $S_{a}$, $D_{y}$ and $S_{c,1}$,
see $\Sigma^{\prime(2)}$ in Eq. (30). The Green’s functions of $N^{\prime}$
are given by

$\displaystyle\hat{g}_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}=A_{R}+iA_{I}$
(58)
$\displaystyle\hat{g}_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}=-A_{R}+iA_{I}$
(59)
$\displaystyle\hat{g}_{N^{\prime}_{x},N^{\prime}_{x}}^{1,1}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{y}}^{1,1}=B_{R}+iB_{I}$
(60)
$\displaystyle\hat{g}_{N^{\prime}_{x},N^{\prime}_{x}}^{2,2}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{y}}^{2,2}=-B_{R}+iB_{I},$
(61)

where $(A_{R},\,A_{I})$, and $(B_{R},\,B_{I})$ are four real-valued parameters
for the nonlocal and local Green’s functions
$\hat{g}_{N^{\prime}_{x},N^{\prime}_{y}}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{x}}$
and $\hat{g}_{N^{\prime}_{x},N^{\prime}_{x}}$,
$\hat{g}_{N^{\prime}_{y},N^{\prime}_{y}}$ respectively, within the assumption
$\hat{g}_{N^{\prime}_{x},N^{\prime}_{x}}=\hat{g}_{N^{\prime}_{y},N^{\prime}_{y}}$.
All of the numerical calculations presented below are within the assumption
$B_{R}=B_{I}=0$, which is justified as mimicking perfect transmission for
highly transparent extended interfaces between the superconductors and the
ballistic 2D metal. This assumption yields vanishingly small value for the
local Green’s function involving “U-turn” of the quasiparticle coming back on
the interface, a process that is not possible at high transparency in the
absence of disorder. In addition to being physically motivated, the condition
$B_{R}=B_{I}=0$ allows restricting the parameter space to be explored in the
numerical calculations.

The assumption of 0D quantum dots implies that the quartet critical current
depends on the value of the Green’s function crossing the conductor
$N^{\prime}$ between $N^{\prime}_{x}$ and $N^{\prime}_{y}$, see $N^{\prime}$,
$N^{\prime}_{x}$ and $N^{\prime}_{y}$ in figure 4. The nonlocal Green’s
function crossing $N^{\prime}$ oscillates with the combinations
$\cos(k_{F}R_{0})$ and $\sin(k_{F}R_{0})$, where $R_{0}$ is the separation
between $N^{\prime}_{x}$ and $N^{\prime}_{y}$, and $k_{F}$ is the Fermi wave-
vector. Sharp resonance peaks emerge at fixed $k_{F}R_{0}$ in the voltage-$V$
dependence of the quartet current, see the previous section III and the
numerical results in the present section. The interplay between those
$k_{F}R_{0}$-oscillations in space and the sharp Floquet resonances in the
voltage-dependence of the quartet critical current is expected to produce log-
normal distribution of strong sample-to-sample fluctuations, where different
samples are characterized by different values of $R_{0}$. In this section, we
specifically investigate the regime of strong log-normal distribution of
sample-to-sample fluctuations and make the physically-motivated assumption
that fixing $A_{R}$, $A_{I}$ yields quartet critical current-voltage
dependence that is qualitatively representative of the signal in this regime
of strong sample-to-sample fluctuations.

The codes have been developed over the last few years HGE ; Sotto ; FWS ;
engineering ; paperII ; long-distance ; Heiblum ; Berry . They are based on
recursive calculations as a function of the harmonics of the Josephson
frequency Cuevas ; AVERIN  (see also the Appendix of Ref. Sotto, ) and on
sparse matrix algorithms for matrix products. We specifically adapted the code
of our recent Ref. long-distance,  to include the local and nonlocal Green’s
functions given by Eqs. (58)-(61).

In the forthcoming figures 7 and 8, we present numerical results for the
voltage-dependence of the quartet critical current $I_{q,c}$ defined as

$\displaystyle I_{q,c}(\Phi=0,eV/\Delta)=$
$\displaystyle\mbox{max}_{\varphi_{q}}\left[I_{q}(\Phi=0,\varphi_{q},eV/\Delta)-I_{q}(\Phi=0,-\varphi_{q},eV/\Delta)\right]$
$\displaystyle-$
$\displaystyle\mbox{min}_{\varphi_{q}}\left[I_{q}(\Phi=0,\varphi_{q},eV/\Delta)-I_{q}(\Phi=0,-\varphi_{q},eV/\Delta)\right]$
$\displaystyle I_{q,c}(\Phi=\pi,eV/\Delta)=$
$\displaystyle\mbox{max}_{\varphi_{q}}\left[I_{q}(\Phi=\pi,\varphi_{q},eV/\Delta)-I_{q}(\Phi=\pi,-\varphi_{q},eV/\Delta)\right]$
$\displaystyle-$
$\displaystyle\mbox{min}_{\varphi_{q}}\left[I_{q}(\Phi=\pi,\varphi_{q},eV/\Delta)-I_{q}(\Phi=\pi,-\varphi_{q},eV/\Delta)\right].$

In addition, we present colormaps for the sign of $I_{q,c}$ as a function of
the model parameters (see the forthcoming figures 9 and 10).

Figure 7 shows the normalized bias voltage-$eV/\Delta$-dependence of the
quartet critical currents $I_{q,c}(\Phi=0,eV/\Delta)$ and
$I_{q,c}(\Phi=\pi,eV/\Delta)$ at the flux values $\Phi=0$ and $\Phi=\pi$. The
coupling parameters $\Gamma/\Delta=1$ and $\Gamma^{\prime}/\Delta=0$ are used
in figure 7.

Figures 7a, 7b and 7c correspond to absence of relaxation with $A_{I}=0$,
producing small quartet critical current. Figures 7d, 7e, 7f and 7g, 7h, 7i
correspond to the nonvanishingly small relaxation parameters $A_{I}=0.05/W$
and $A_{I}=0.1/W$ respectively, producing large quartet critical current at
small voltage ratio $eV/\Delta$. In addition, sharp resonances emerge in the
variations of the quartet critical current as a function of $eV/\Delta$ for
$A_{I}=0.05/W$ and $A_{I}=0.1/W$ in figures 7d, 7e, 7f and 7g, 7h, 7i
respectively. The numerical results in figure 7 are in agreement with the
general mechanism of section III.

Figure 8 shows the normalized voltage-$eV/\Delta$-dependence of the quartet
critical current $I_{q,c}$, in semi-logarithmic scale, now with nonvanishingly
small $\lambda=\Gamma^{\prime}/\Gamma$, and with $A_{I}=0.05/W$. Those sharp
peaks in the quartet critical current are expected to change as the value of
the normal metal Green’s functions is varied in Eqs. (58)-(61). The quartet
critical current plotted as a function of the bias voltage $V$ is normal in
log-scale. The corresponding peaks are expected to evolve as $k_{F}R_{0}$ is
varied and thus, sample-to-sample fluctuations are expected to have log-normal
distribution.

In addition, multiple cross-overs are obtained in figure 8 between
$I_{q,c}(\Phi=\pi,eV/\Delta)-I_{q,c}(\Phi=0,eV/\Delta)>0$ (i.e. inversion) and
$I_{q,c}(\Phi=\pi,eV/\Delta)-I_{q,c}(\Phi=0,eV/\Delta)<0$ (i.e. noninverted
behavior) in a low-voltage window.

Figure 9 shows the sign of
$I_{q,c}(\lambda,\Phi=\pi,eV/\Delta)-I_{q,c}(\lambda,\Phi=0,eV/\Delta)$ as a
function of the parameters $\lambda=\Gamma^{\prime}/\Gamma$ (on the $x$-axis)
and $eV/\Delta$ (on the $y$-axis). Inversion
$I_{q,c}(\lambda,\Phi=\pi,eV/\Delta)-I_{q,c}(\lambda,\Phi=0,eV/\Delta)>0$ is
obtained at low $eV/\Delta$ for $\lambda=0$, see the previous section IV and
Appendix B for a discussion of the inversion in the $eV/\Delta=0^{+}$ limit.
In figure 9, increasing the coupling $\Gamma^{\prime}$ between $D_{x}$ and
$(S_{b},S_{c,2})$ or between $D_{y}$ and $(S_{a},S_{c,1})$ [in addition to
$\Gamma$ between $D_{x}$ and $(S_{a},S_{c,1})$ or between $D_{y}$ and
$(S_{b},S_{c,2})$], makes the double 0D quantum dot behave closer to a pair of
single quantum dots. This favors the “noninverted behavior” typical of single
quantum dots, as opposed to the “inverted behavior” appearing at
$\lambda=\Gamma^{\prime}/\Gamma=0$ in a double 0D quantum dot, see Appendix B.

The values $\lambda=\Gamma^{\prime}/\Gamma$ such that
$0.2\lesssim\lambda\lesssim 0.6$ in figure 9 generically yield “inversion” at
low bias-voltage $V$, followed by “absence of inversion” at higher $V$-values.
Now, we show that “absence of inversion” can appear at low bias-voltage $V$ if
the on-site energies $\epsilon_{x}=\epsilon_{y}\equiv\epsilon_{0}$ in the
double quantum dot Hamiltonian given by Eqs. (6)-(7) produce detuning from
perfectly opposite levels. Namely, the parameter
$\lambda=\Gamma^{\prime}/\Gamma$ is set to $\lambda=1/2$ in figure 10. Moving
away from $\epsilon_{0}/\Delta=0$ in figure 10 produces typical voltage-
dependence of
$I_{q,c}(\lambda,\Phi=\pi,eV/\Delta)-I_{q,c}(\lambda,\Phi=0,eV/\Delta)$ with
“noninverted behavior” at low $V<V_{*}$ and “inversion” at higher $V>V_{*}$,
see for instance the moderately small values
$-0.5\lesssim\epsilon_{0}/\Delta\lesssim 0.5$ in figure 10. With the
considered parameters, the ratio $eV_{*}/\Delta$ is in the range
$eV_{*}/\Delta\approx 0.1$. Thus, our model is compatible with the
experimental data HGE  shown in figure 1b.

##  VI Low-voltage limit

In this section, we focus on the low-voltage limit and insert in Eq. (III.3)
the equilibrium fully dressed advanced and retarded Green’s functions. We
present an explanation for why the numerical calculations of the preceding
section V suggest that the line-width broadening $\delta_{l}$ in Eq. (III.3)
appears to be much smaller than $\sim A_{I}$ in Eqs. (58)-(61).

As in the previous numerical calculations presented in section V, we consider
normal leads in series, see figure 4 and the Hamiltonian in subsection II.2.3.
The Dyson Eqs. (40) take the form

$\displaystyle\left((\omega-i\eta)-\hat{\Gamma}_{x,x}^{(a,a)}-\hat{\Gamma}_{x,x}^{(c_{1},c_{1})}-\hat{\Gamma}_{x,x}^{(N^{\prime}_{x},N^{\prime}_{x})}\right)\hat{G}_{x,x}$
$\displaystyle-\hat{\Gamma}_{x,y}^{(N^{\prime}_{x},N^{\prime}_{y})}\hat{G}_{y,x}=\hat{g}_{x,x}$
$\displaystyle\left((\omega-i\eta)-\hat{\Gamma}_{y,y}^{(b,b)}-\hat{\Gamma}_{y,y}^{(c_{1},c_{1})}-\hat{\Gamma}_{y,y}^{(N^{\prime}_{y},N^{\prime}_{y})}\right)\hat{G}_{y,x}$
$\displaystyle-\hat{\Gamma}_{y,x}^{(N^{\prime}_{y},N^{\prime}_{x})}\hat{G}_{x,x}=0,$

with

$\displaystyle\hat{\Gamma}_{x,x}^{(a,a)}$ $\displaystyle=$
$\displaystyle-\Gamma_{a}\left(\begin{array}[]{cc}0&\exp(i\varphi_{a})\\\
\exp(-i\varphi_{a})&0\end{array}\right)$ (68)
$\displaystyle\hat{\Gamma}_{x,x}^{(c_{1},c_{1})}$ $\displaystyle=$
$\displaystyle-\Gamma_{c_{1}}\left(\begin{array}[]{cc}0&\exp(i\varphi_{c,1})\\\
\exp(-i\varphi_{c,1})&0\end{array}\right)$ (71)
$\displaystyle\hat{\Gamma}_{y,y}^{(b,b)}$ $\displaystyle=$
$\displaystyle-\Gamma_{b}\left(\begin{array}[]{cc}0&\exp(i\varphi_{b})\\\
\exp(-i\varphi_{b})&0\end{array}\right)$ (74)
$\displaystyle\hat{\Gamma}_{y,y}^{(c_{2},c_{2})}$ $\displaystyle=$
$\displaystyle-\Gamma_{c_{2}}\left(\begin{array}[]{cc}0&\exp(i\varphi_{c,2})\\\
\exp(-i\varphi_{c,2})&0\end{array}\right),$ (77)

where the couplings $\Gamma_{a}$, $\Gamma_{b}$, $\Gamma_{c,1}$ and
$\Gamma_{c,2}$ are between the double quantum dot and the superconducting
leads $S_{a}$, $S_{b}$, $S_{c,1}$ and $S_{c,2}$, see figure 4 and the
analogous Eqs. (16)-(19). In addition,
$\hat{\Gamma}_{x,x}^{(N^{\prime}_{x},N^{\prime}_{x})}$,
$\hat{\Gamma}_{x,y}^{(N^{\prime}_{x},N^{\prime}_{y})}$,
$\hat{\Gamma}_{y,y}^{(N^{\prime}_{y},N^{\prime}_{y})}$ and
$\hat{\Gamma}_{y,x}^{(N^{\prime}_{y},N^{\prime}_{x})}$ are given by

$\displaystyle\hat{\Gamma}_{x,x}^{(N^{\prime}_{x},N^{\prime}_{x})}=\Gamma_{N^{\prime}_{x},N^{\prime}_{x}}\left(\begin{array}[]{cc}Wg_{N^{\prime}_{x},N^{\prime}_{x}}^{1,1}&0\\\
0&Wg_{N^{\prime}_{x},N^{\prime}_{x}}^{2,2}\end{array}\right)$ (80)
$\displaystyle\hat{\Gamma}_{x,y}^{(N^{\prime}_{x},N^{\prime}_{y})}=\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\left(\begin{array}[]{cc}Wg_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}&0\\\
0&Wg_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}\end{array}\right)$ (83)
$\displaystyle\hat{\Gamma}_{y,y}^{(N^{\prime}_{y},N^{\prime}_{y})}=\Gamma_{N^{\prime}_{y},N^{\prime}_{y}}\left(\begin{array}[]{cc}Wg_{N^{\prime}_{y},N^{\prime}_{y}}^{1,1}&0\\\
0&Wg_{N^{\prime}_{y},N^{\prime}_{y}}^{2,2}\end{array}\right)$ (86)
$\displaystyle\hat{\Gamma}_{y,x}^{(N^{\prime}_{y},N^{\prime}_{x})}=\Gamma_{N^{\prime}_{y},N^{\prime}_{x}}\left(\begin{array}[]{cc}Wg_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}&0\\\
0&Wg_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\end{array}\right),$ (89)

with
$\Gamma_{N^{\prime}_{x},N^{\prime}_{x}}=\left(\Sigma_{x,N^{\prime}_{x}}\right)^{2}/W$,
$\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}=\Sigma_{x,N^{\prime}_{x}}\Sigma_{y,N^{\prime}_{y}}/W$,
$\Gamma_{N^{\prime}_{y},N^{\prime}_{x}}=\Sigma_{x,N^{\prime}_{y}}\Sigma_{x,N^{\prime}_{x}}/W$,
$\Gamma_{N^{\prime}_{y},N^{\prime}_{y}}=\left(\Sigma_{x,N^{\prime}_{y}}\right)^{2}/W$,
where $N^{\prime}_{x}$ and $N^{\prime}_{y}$ are the normal-metal tight-binding
sites making the contacts with the 0D quantum dots $D_{x}$ and $D_{y}$
respectively.

Figure 11: The real and imaginary parts of the ABS energies in the infinite-
gap limit: Panels (a1)-(c1) and (a2)-(c2) show the real and imaginary parts of
the ABS energies respectively, for the $4\times 4$ Hamiltonian given by Eq.
(68)-(90), plotted as a function of the variable $\psi^{\prime}$ defined as
$\varphi_{a}=\varphi_{a}^{(0)}+\psi^{\prime}$,
$\varphi_{b}=\varphi_{b}^{(0)}-\psi^{\prime}$
$\varphi_{c,1}=\varphi_{c}^{(0)}$ and $\varphi_{c,2}=\varphi_{c}^{(0)}$, with
$\varphi_{q}=\varphi_{a}^{(0)}+\varphi_{b}^{(0)}-2\varphi_{c}^{(0)}$. The
parameters take the general values indicated on the figure, with in addition
$\Gamma_{a}=\Gamma_{b}=\Gamma_{c,1}=\Gamma_{c,2}=1$ taken as the unit of
energy for the coupling to the superconducting leads, and $\Gamma_{loc}=0.7$,
$\Gamma_{non\,loc}=0.5$ for the “local” and “nonlocal” electron-electron
channels, see Eqs. (80)-(89). Panels a2, b2 and c2 reveal that
$|\mbox{Im}(\lambda_{i})|$ can be small compared to $A_{I}$ or $B_{I}$ (but
not orders of magnitude smaller), which implies small $\delta_{l}$ in Eq.
(III.3), and compatibility with the large quartet current reported in the
numerical results presented in the above section V.

The Dyson Eqs. (VI)-(VI) define the infinite-gap Hamiltonian

${\cal
H}^{(\infty)}=\left(\begin{array}[]{cc}\hat{\Gamma}_{x,x}^{(a,a)}+\hat{\Gamma}_{x,x}^{(c_{1},c_{1})}+\hat{\Gamma}_{x,x}^{(N^{\prime}_{x},N^{\prime}_{x})}&\hat{\Gamma}_{x,y}^{(N^{\prime}_{x},N^{\prime}_{y})}\\\
\hat{\Gamma}_{y,x}^{(N^{\prime}_{y},N^{\prime}_{x})}&\hat{\Gamma}_{y,y}^{(b,b)}+\hat{\Gamma}_{y,y}^{(c_{2},c_{2})}+\hat{\Gamma}_{y,y}^{(N^{\prime}_{y},N^{\prime}_{y})}\end{array}\right),$
(90)

which is now considered in absence of the local couplings
$\hat{\Gamma}_{x,x}^{(N^{\prime}_{x},N^{\prime}_{x})}=\hat{\Gamma}_{y,y}^{(N^{\prime}_{y},N^{\prime}_{y})}\equiv
0$, corresponding to the use of $B_{R}=B_{I}=0$ in the above numerical
calculations, see the discussion following Eqs. (58)-(61) in section V. Then,
the square of the $4\times 4$ infinite-gap Hamiltonian decouples into a pair
of $2\times 2$ blocks, see also Appendix B:

$\left({\cal
H}^{(\infty)}\right)^{2}=\left(\begin{array}[]{cccc}\epsilon_{a,c_{1}}^{1,1}&0&0&\kappa_{1,2}\\\
0&\epsilon_{a,c_{1}}^{2,2}&\kappa_{2,1}&0\\\
0&\kappa^{\prime}_{1,2}&\epsilon_{b,c_{2}}^{1,1}&0\\\
\kappa^{\prime}_{2,1}&0&0&\epsilon_{b,c_{2}}^{2,2}\end{array}\right),$ (91)

with

$\displaystyle\epsilon_{a,c_{1}}^{1,1}$ $\displaystyle=$
$\displaystyle|\gamma_{D_{x},D_{x}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}$
(92) $\displaystyle\epsilon_{a,c_{1}}^{2,2}$ $\displaystyle=$
$\displaystyle|\gamma_{D_{x},D_{x}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}$
(93) $\displaystyle\epsilon_{b,c_{2}}^{1,1}$ $\displaystyle=$
$\displaystyle|\gamma_{D_{y},D_{y}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}$
(94) $\displaystyle\epsilon_{b,c_{2}}^{2,2}$ $\displaystyle=$
$\displaystyle|\gamma_{D_{y},D_{y}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}$
(95) $\displaystyle\kappa_{1,2}$ $\displaystyle=$
$\displaystyle\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}W\left(g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}\gamma_{D_{x},D_{x}}+g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}\gamma_{D_{y},D_{y}}\right)$
(96) $\displaystyle\kappa_{2,1}$ $\displaystyle=$
$\displaystyle\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}W\left(g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}\overline{\gamma_{D_{x},D_{x}}}+g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}\overline{\gamma_{D_{y},D_{y}}}\right)$
(97) $\displaystyle\kappa^{\prime}_{1,2}$ $\displaystyle=$
$\displaystyle\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}W\left(g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}\gamma_{D_{x},D_{x}}+g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\gamma_{D_{y},D_{y}}\right)$
(98) $\displaystyle\kappa^{\prime}_{2,1}$ $\displaystyle=$
$\displaystyle\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}W\left(g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\overline{\gamma_{D_{x},D_{x}}}+g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}\overline{\gamma_{D_{y},D_{y}}}\right),$
(99)

where

$\displaystyle\gamma_{D_{x},D_{x}}$ $\displaystyle=$
$\displaystyle\Gamma_{a}\exp(i\varphi_{a})+\Gamma_{c_{1}}\exp(i\varphi_{c,1})$
(100) $\displaystyle\gamma_{D_{y},D_{y}}$ $\displaystyle=$
$\displaystyle\Gamma_{b}\exp(i\varphi_{b})+\Gamma_{c_{2}}\exp(i\varphi_{c,2}).$
(101)

Within each $2\times 2$ block, the eigenvalues take the form

$\lambda_{\pm}=\frac{1}{2}\left[\epsilon_{a,c_{1}}^{1,1}+\epsilon_{b,c_{2}}^{2,2}\pm\sqrt{\left(\epsilon_{a,c_{1}}^{1,1}-\epsilon_{b,c_{2}}^{2,2}\right)^{2}+4\kappa_{1,2}\kappa^{\prime}_{2,1}}\right],$
(102)

where the following quantities

$\displaystyle\epsilon_{a,c_{1}}^{1,1}+\epsilon_{b,c_{2}}^{2,2}=$ (103)
$\displaystyle|\gamma_{D_{x},D_{x}}|^{2}+|\gamma_{D_{y},D_{y}}|^{2}+2\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}\left(A_{R}^{2}-A_{I}^{2}\right)$
$\displaystyle\left(\epsilon_{a,c_{1}}^{1,1}-\epsilon_{b,c_{2}}^{2,2}\right)^{2}+4\kappa_{1,2}\kappa^{\prime}_{2,1}=$
(104)
$\displaystyle\left(|\gamma_{D_{x},D_{x}}|^{2}-|\gamma_{D_{y},D_{y}}|^{2}\right)^{2}-16\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}A_{R}^{2}A_{I}^{2}$
$\displaystyle+4\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}A_{R}^{2}\left|\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}\right|^{2}$
$\displaystyle-4\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}A_{I}^{2}\left|\gamma_{D_{x},D_{x}}+\gamma_{D_{y},D_{y}}\right|^{2}$

are both real-valued.

It can also be shown that $\lambda_{\pm}$ is real-valued and positive at small
$A_{I}$ because

$\displaystyle\left|\gamma_{D_{x},D_{x}}\overline{\gamma_{D_{y},D_{y}}}-\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\right|^{2}$
(105)

is positive, which implies

$\displaystyle|\gamma_{D_{x},D_{x}}|^{2}|\gamma_{D_{y},D_{y}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{4}W^{4}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}$
(106) $\displaystyle>$
$\displaystyle\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\gamma_{D_{y},D_{y}}\overline{\gamma_{D_{x},D_{x}}}$
$\displaystyle+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}\gamma_{D_{x},D_{x}}\overline{\gamma_{D_{y},D_{y}}},$

where we used the approximation
$\overline{g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}}\simeq
g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}$ and
$\overline{g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}}\simeq
g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}$ at small $A_{I}$. Thus, we obtain

$\displaystyle\left(|\gamma_{D_{x},D_{x}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}\right)$
$\displaystyle\times\left(|\gamma_{D_{y},D_{y}}|^{2}+\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\right)$
$\displaystyle>$
$\displaystyle\left(\Gamma_{N^{\prime}_{x},N^{\prime}_{y}}\right)^{2}W^{2}\left(g_{N^{\prime}_{x},N^{\prime}_{y}}^{2,2}\gamma_{D_{x},D_{x}}+g_{N^{\prime}_{x},N^{\prime}_{y}}^{1,1}\gamma_{D_{y},D_{y}}\right)$
$\displaystyle\times\left(g_{N^{\prime}_{y},N^{\prime}_{x}}^{2,2}\overline{\gamma_{D_{x},D_{x}}}+g_{N^{\prime}_{y},N^{\prime}_{x}}^{1,1}\overline{\gamma_{D_{y},D_{y}}}\right),$

from what we deduce

$\left(\epsilon_{a,c_{1}}^{1,1}+\epsilon_{b,c_{2}}^{2,2}\right)^{2}>\left(\epsilon_{a,c_{1}}^{1,1}-\epsilon_{b,c_{2}}^{2,2}\right)^{2}+4\kappa_{1,2}\kappa^{\prime}_{2,1}$
(108)

at small $A_{I}$. Then, the eigenvalues $\lambda_{\pm}$ of the squared
infinite-gap Hamiltonian are both real-valued at arbitrary $A_{I}$ [see Eqs.
(103)-(104)] and positive at small $A_{I}$ [see Eq. (108)].

In addition, we carried out complementary numerical calculation for the
eigenvalues of the equilibrium $4\times 4$ infinite-gap Hamiltonian defined by
Eqs. (68)-(90), see figure 9 We find the possibility of imaginary part of the
ABS energies smaller (but not orders of magnitude smaller) than $A_{I}$ and
$B_{I}$ with general values of the parameters, i.e. not within the above-
considered $B_{R}=B_{I}=0$ and small $A_{I}$. Thus, small $\delta_{l}$ has to
be used in Eq. (III.3) at low bias voltage, which is compatible with the large
quartet current in figure 8, similarly to the thermal noise of a
superconducting weak link Yeyati-noise-equilibrium .

##  VII Conclusions

Now, we provide final remarks. The paper is summarized in subsection VII.1,
and perspectives are presented in subsection VII.2.

###  VII.1 Summary of the paper

We proposed models of ballistic multiterminal Josephson junctions made with
superconducting leads evaporated on a 2D metal. The models consist of
multilevel quantum dots connected to superconducting and normal leads. We
argued that the qualitative physics of the two-Cooper pair resonance can be
captured with phenomenological double quantum dots connected to
superconducting and normal leads. The coupling to the nonproximitized regions
of the ballistic conductor produces relaxation, and a fraction of the quartet-
phase sensitive current is transmitted into the normal parts of the circuit.

We found resonances in calculations that account for both the time-periodic
dynamics and small relaxation, if the multilevel or double quantum dot support
levels at opposite energies. A related effect was previously found in the
thermal noise of a superconducting weak link Yeyati-noise-equilibrium . In
addition, we addressed detuning from perfectly opposite energy levels.

Noninverted-to-inverted cross-over numerically emerges as the bias voltage
increases, where “inversion” means “larger quartet critical current at half
flux-quantum $\Phi=\pi$ than in zero field at $\Phi=0$”. The corresponding
cross-over voltage $V_{*}$ is small compared to the superconducting gap,
typically $eV_{*}\approx\Delta/10$ with the parameters of our calculations,
which is compatible with the recent Harvard group experiment HGE .

###  VII.2 Perspectives

A challenge for the theory is to model devices that are quite complex, with
e.g. extended interfaces, four or more superconducting leads, nonequilibrium
voltage biasing conditions, loops connecting superconducting terminals,
possibly with radio-frequency radiation. Direct diagonalizations of the
Bogoliubov-de Gennes Hamiltonian would apparently lead to prohibitive
computational expenses. In the field of mesoscopic superconductivity, Nazarov
and co-workers (see Ref. Nazarov-book,  and references therein) proposed and
developed finite element theory for the superconducting-normal metal circuits
that describe the proximity effect, i.e. the interplay between Andreev
reflection and multiple scattering on disorder. This dirty-limit circuit
theory was proposed for multiterminal Josephson junctions Padurariu ;
Padurariu2 . The dirty limit implies short elastic mean free path, a condition
that is not directly met in the ballistic metals that are currently used in
some experiments on superconducting hybrid structures, such as carbon
nanotubes Pillet-CNT , semiconducting nanowires multiterminal-exp7  or
graphene Bretheau3 ; Park . Specifically, tunneling spectroscopy of carbon
nanotube Josephson junctions Pillet-CNT  revealed discrete ABS. Andreev
molecules were realized with semiconducting nanowires multiterminal-exp7 .
Evidence for superconducting phase difference-sensitive continuum of ABS was
obtained in superconductor-graphene-superconductor Josephson junctions
Bretheau1 ; Bretheau2 . As it is mentioned above, microwave experiments on
short superconductor-graphene-superconductor Josephson junctions were recently
carried out, and modeled with single-level quantum dots Park .

Several issues related to averaging could be investigated in the future:

First, it would be interesting to numerically average over a distribution of
the nonlocal Green’s function for double 0D quantum dots. A related issue is
to implement quantum dots having dimension that is large compared to the Fermi
wave-length $\lambda_{F}$ instead of the 0D quantum dots of the present paper.
On the other hand, the “metallic” regime of weak sample-to-sample fluctuations
of the quartet critical current can be addressed with quasiclassics. For
instance, discretized Usadel equations were used to evaluate multiple Andreev
reflections in two-terminal devices in the dirty limit Cuevas-Usadel  and
three-terminal devices were addressed within assumptions about the interface
transparencies, also with Usadel equations Cuevas-Pothier .

Second, the use of Nazarov’s circuit theory implies that the Green’s function
is uniform within a node, which is also satisfied by ballistic chaotic
cavities Vanevic . Thus Nazarov’s circuit theory is also appropriate to
describe a ballistic device if the coupling to the terminals is made through
interfaces with small cross-sections, to ensure that the metal behaves like a
chaotic cavity.

Third, an interesting perspective is to solve Eilenberger equations for the
four-terminal device shown in figure 1.

Fourth, it would also be interesting to average over voltage fluctuations, the
strength of which being controlled by the electromagnetic environment.

To summarize this discussion, perspectives are about the interplay between:
(i) The “Floquet effects” that produce “non-selfaveraging-like” sharp
resonance peaks in the voltage-dependence of the quartet signal for single-
channel contacts, and (ii) The effect of averaging in space over extended
contacts or in energy over the voltage fluctuations induced by the
electromagnetic environment.

Finally, probability current is conserved in the one-dimensional normal metal-
superconductor junction treated by Blonder, Tinkham and Klapwijk BTK  [see
Eqs. (A7)-(A8) in this paper]. We leave as open the question of rigorously
discussing probability conservation in the considered four-terminal device
connected to normal leads, in connection with introducing the phenomenological
Dynes parameter $\eta\neq 0$ and vanishingly small imaginary part of the local
Green’s function $B_{I}=0$ (as it was done here), or with $\eta=0$ and
$B_{I}\neq 0$ (as it could be done in the future).

To conclude, the present paper suggests interest of “quantum bath engineering”
multiterminal Josephson junctions in the circuits of cavity-quantum
electrodynamics.

## Acknowledgements

The author wishes to thank K. Huang, Y. Ronen and P. Kim for stimulating
discussions about their experiment. The author wishes to thank R. Danneau for
useful discussions and comments on the manuscript, and F. Levy-Bertrand and
her colleagues H. Cercellier, K. Hasselbach, M.A. Measson for useful remarks
during an informal seminar on this topic. The author thanks the Infrastructure
de Calcul Intensif et de Données (GRICAD) for use of the resources of the
Mésocentre de Calcul Intensif de l’Université Grenoble-Alpes (CIMENT). The
author acknowledges support from the French National Research Agency (ANR) in
the framework of the Graphmon project (ANR-19-CE47-0007).

##  Appendix A Spectrum in absence of coupling to the superconductors

In this Appendix, we consider the “phenomenological model of double 0D quantum
dots with normal lead in series”, see subsection II.2.3, and show that pairs
of opposite energies emerge in the spectrum.

We start with the simple limit where the four superconducting leads are
disconnected, i.e. $\Gamma=\Gamma^{\prime}=0$, and consider that the two
quantum dots $D_{x}$ and $D_{y}$ are connected to the tight-binding sites
$N^{\prime}_{x}$ and $N^{\prime}_{y}$ by the tunneling amplitudes
$\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}=\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}$ and
$\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}=\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}$. The
tight-binding sites are self-connected by the Green’s functions
$g_{N^{\prime}_{x},N^{\prime}_{x}}$ and $g_{N^{\prime}_{y},N^{\prime}_{y}}$,
and connected to each other by the nonlocal Green’s functions
$g_{N^{\prime}_{x},N^{\prime}_{y}}=g_{N^{\prime}_{y},N^{\prime}_{x}}$. We
denote by $g_{D_{x},D_{x}}^{A}=g_{D_{y},D_{y}}^{A}=1/(\omega-i\eta)$ the
Green’s functions of the “isolated” quantum dots $D_{x}$ and $D_{y}$, see Eqs.
(9)-(10) where $\epsilon_{x}=\epsilon_{y}=0$ in Eqs. (6)-(7). In addition
$\tilde{g}_{D_{x},D_{x}}$, $\tilde{g}_{D_{y},D_{y}}$ are their counterparts
for the connected double 0D quantum dot, and by $\tilde{g}_{D_{x},D_{y}}$ and
$\tilde{g}_{D_{y},D_{x}}$ are the corresponding nonlocal Green’s functions.
The Dyson equations relate the $g$s to the $\tilde{g}$s according to

$\displaystyle\tilde{g}_{D_{x},D_{x}}$ $\displaystyle=$ $\displaystyle
g_{D_{x},D_{x}}+g_{D_{x},D_{x}}\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}g_{N^{\prime}_{x},N^{\prime}_{x}}\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}\tilde{g}_{D_{x},D_{x}}$
$\displaystyle+g_{D_{x},D_{x}}\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}g_{N^{\prime}_{x},N^{\prime}_{y}}\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}\tilde{g}_{D_{y},D_{x}}$
$\displaystyle\tilde{g}_{D_{y},D_{x}}$ $\displaystyle=$ $\displaystyle
g_{D_{y},D_{y}}\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}g_{N^{\prime}_{y},N^{\prime}_{x}}\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}\tilde{g}_{D_{x},D_{x}}$
$\displaystyle+g_{D_{y},D_{y}}\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}g_{N^{\prime}_{y},N^{\prime}_{y}}\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}\tilde{g}_{D_{y},D_{x}},$

which leads to the secular equation

$\left|\begin{array}[]{cc}\omega-\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{x},D_{x}}&-\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{y},D_{y}}\\\
-\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{x},D_{x}}&\omega-\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{y},D_{y}}\end{array}\right|=0,$
(111)

where
$\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{x},D_{x}}=\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}g_{N^{\prime}_{x},N^{\prime}_{x}}\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}$,
$\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{y},D_{y}}=\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}g_{N^{\prime}_{y},N^{\prime}_{y}}\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}$,
$\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{y},D_{y}}=\Sigma^{(4)}_{D_{x},N^{\prime}_{x}}g_{N^{\prime}_{x},N^{\prime}_{y}}\Sigma^{(4)}_{N^{\prime}_{y},D_{y}}$,
$\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{x},D_{x}}=\Sigma^{(4)}_{D_{y},N^{\prime}_{y}}g_{N^{\prime}_{y},N^{\prime}_{x}}\Sigma^{(4)}_{N^{\prime}_{x},D_{x}}$.
Assuming symmetric contacts yields
$\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{x},D_{x}}=\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{y},D_{y}}\equiv\Gamma_{loc}$
and
$\Gamma_{D_{x},N^{\prime}_{x},N^{\prime}_{y},D_{y}}=\Gamma_{D_{y},N^{\prime}_{y},N^{\prime}_{x},D_{x}}\equiv\Gamma_{non\,loc}$.
The energy levels are given by

$\displaystyle\omega_{(\pm)}$ $\displaystyle=$
$\displaystyle\Gamma_{loc}\pm\Gamma_{non\,loc}$ (112) $\displaystyle=$
$\displaystyle\left(\Sigma^{(4)}\right)^{2}\left[B_{R}\pm
A_{R}\right]+i\left(\Sigma^{(4)}\right)^{2}\left[B_{I}\pm A_{I}\right],$ (113)

where $A_{R}$, $A_{I}$, $B_{R}$ and $B_{I}$ are given by Eqs. (58)-(61).

Now, we assume $B_{R}=B_{I}=0$, as in section V. Then, Eqs. (112)-(113) yield
energy levels $\omega_{(\pm)}=\pm\Sigma^{(4)}\left(A_{R}+iA_{I}\right)$ having
opposite real and imaginary parts.

##  Appendix B $V=0^{+}$ adiabatic limit

In this Appendix, we examine the $V=0^{+}$ adiabatic limit of four-terminal
Josephson junctions containing a single or two quantum dots (see subsections
B.1 and B.2 below).

###  B.1 Single quantum dot

We start with single 0D quantum dots in the $V=0^{+}$ adiabatic limit,
summarizing a fraction of the Supplemental Material of our previous paper II
paperII .

The Dyson equations take the following form for the considered 0D quantum dot
connected to $p_{0}$ superconducting leads by the tunnel amplitudes
$\hat{\Sigma}^{(5)}_{D_{x},S_{n_{p}}}=\hat{\Sigma}^{(5)}_{S_{n_{p}},D_{x}}$,
with $n_{p}=1,...,p_{0}$:

$\hat{G}_{D_{x},D_{x}}=\hat{g}_{D_{x},D_{x}}+\hat{g}_{D_{x},D_{x}}\sum_{n_{p}=1}^{p_{0}}\hat{\Sigma}^{(5)}_{D_{x},S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{S_{n_{p}},D_{x}}\hat{G}_{D_{x},D_{x}}.$
(114)

In the infinite-gap limit, Eq. (114) can be expressed with the infinite-gap
Hamiltonian $\hat{\cal H}_{eff,\,single\,dot}^{(\infty)}$:

$\hat{G}_{D_{x},D_{x}}^{A}=\left(\omega-i\eta-\hat{\cal
H}_{eff,\,single\,dot}^{(\infty)}\right)^{-1},$ (115)

where

$\hat{\cal
H}_{eff,\,single\,dot}^{(\infty)}=\sum_{n_{p}=1}^{p_{0}}\hat{\Sigma}^{(5)}_{D_{x},S_{n_{p}}}\hat{g}_{S_{n_{p}},S_{n_{p}}}\hat{\Sigma}^{(5)}_{S_{n_{p}},D_{x}}.$
(116)

Specifically, we obtain the following with $p_{0}=4$ superconducting leads:

$\hat{\cal
H}_{eff,\,single\,dot}^{(\infty)}=\left(\begin{array}[]{cc}0&\gamma_{D_{x},D_{x}}\\\
\left(\gamma_{D_{x},D_{x}}\right)^{*}&0\end{array}\right),$ (117)

where

$\displaystyle\gamma_{D_{x},D_{x}}$ $\displaystyle=$
$\displaystyle\Gamma_{a}\exp(i\varphi_{a})+\Gamma_{b}\exp(i\varphi_{b})$
$\displaystyle+$
$\displaystyle\Gamma_{c,1}\exp(i\varphi_{c,1})+\Gamma_{c,2}\exp(i\varphi_{c,2}),$

and $\Gamma_{n_{p}}=\left(\Sigma_{D_{x},S_{n_{p}}}^{(5)}\right)^{2}/W$
parameterizes the line-width broadening of the quantum dot level in the normal
state.

The $(S_{c,1},S_{c,2})$ superconducting leads can be gathered into the single
$S_{c,eff}$ coupled by

$\Gamma_{c,eff}=\Gamma_{c,1}\exp\left(i\varphi_{c,1}\right)+\Gamma_{c,2}\exp\left(i\varphi_{c,2}\right).$
(119)

Using the identical $\Gamma_{c,1}=\Gamma_{c,2}\equiv\Gamma_{c}$ and the gauge
given by Eqs. (2)-(3) yields

$\Gamma_{c,eff}=\Gamma_{c}\left[1+\exp\left(i\Phi\right)\right]\exp\left(i\varphi_{c,1}\right).$
(120)

It was shown in the Supplemental Material of our previous paper II paperII
that nonsymmetric coupling to the superconducting leads can produce inversion
between $\Phi=0$ and $\Phi=\pi$ in the infinite-gap limit. But Eq. (120)
implies $|\Gamma_{c,eff}|=2\Gamma_{c}$ if $\Phi=0$ and $\Gamma_{c,eff}=0$ if
$\Phi=\pi$ for symmetric couplings to the superconducting leads, i.e. the
quartet current at $\Phi=\pi$ is vanishingly small, thus it is automatically
smaller than at $\Phi=0$.

###  B.2 Double 0D quantum dot

In this subsection, we provide simple argument for emergence of inversion in
the double 0D quantum dot with normal lead in parallel in figure 3, in the
limits $eV/\Delta=0^{+}$ and $\Gamma^{\prime}=0$.

In the infinite-gap limit, the $4\times 4$ Hamiltonian of a double 0D quantum
dot with normal leads in parallel is given by

$\hat{\cal
H}_{eff,\,double\,dot}^{(\infty)}=\left(\begin{array}[]{cccc}0&\gamma_{D_{x},D_{x}}&\Sigma^{(0)}&0\\\
\left(\gamma_{D_{x},D_{x}}\right)^{*}&0&0&-\Sigma^{(0)}\\\
\Sigma^{(0)}&0&0&\gamma_{D_{y},D_{y}}\\\
0&-\Sigma^{(0)}&\left(\gamma_{D_{y},D_{y}}\right)^{*}&0\end{array}\right),$
(121)

where we assumed $\Gamma^{\prime}=0$, and $\gamma_{D_{x},D_{x}}$,
$\gamma_{D_{y},D_{y}}$ are given by the above Eqs. (100)-(101). Squaring the
infinite-gap Hamiltonian given by Eq. (121) leads to

$\displaystyle\left(\hat{\cal H}_{eff,\,double\,dot}^{(\infty)}\right)^{2}=$
(122)
$\displaystyle\left(\begin{array}[]{cccc}\left|\gamma_{D_{x},D_{x}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}&0&0&-\Sigma^{(0)}\left[\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}\right]\\\
0&\left|\gamma_{D_{x},D_{x}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}&\Sigma^{(0)}\left[\left(\gamma_{D_{x},D_{x}}\right)^{*}-\left(\gamma_{D_{y},D_{y}}\right)^{*}\right]&0\\\
0&\Sigma^{(0)}\left[\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}\right]&\left|\gamma_{D_{y},D_{y}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}&0\\\
-\Sigma^{(0)}\left[\left(\gamma_{D_{x},D_{x}}\right)^{*}-\left(\gamma_{D_{y},D_{y}}\right)^{*}\right]&0&0&\left|\gamma_{D_{y},D_{y}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}\end{array}\right),$
(127)

which decouples into the following $2\times 2$ blocks:

$\displaystyle\left[\hat{\cal H}^{2}\right]^{(0)}_{2\times
2}=\left(\begin{array}[]{cc}\left|\gamma_{D_{x},D_{x}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}&-\Sigma^{(0)}\left[\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}\right]\\\
-\Sigma^{(0)}\left[\left(\gamma_{D_{x},D_{x}}\right)^{*}-\left(\gamma_{D_{y},D_{y}}\right)^{*}\right]&\left|\gamma_{D_{y},D_{y}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}\end{array}\right)$
(130)

and

$\displaystyle\left[\hat{\cal H}^{2}\right]^{(2)}_{2\times
2}=\left(\begin{array}[]{cc}\left|\gamma_{D_{x},D_{x}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}&\Sigma^{(0)}\left[\left(\gamma_{D_{x},D_{x}}\right)^{*}-\left(\gamma_{D_{y},D_{y}}\right)^{*}\right]\\\
\Sigma^{(0)}\left[\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}\right]&\left|\gamma_{D_{y},D_{y}}\right|^{2}+\left(\Sigma^{(0)}\right)^{2}\end{array}
\right).$ (133)

Thus,

$\displaystyle\gamma_{D_{x},D_{x}}-\gamma_{D_{y},D_{y}}$ $\displaystyle=$
$\displaystyle\Gamma_{a}\exp\left(i\varphi_{a}\right)+\Gamma_{c,1}\exp\left(i\varphi_{c,1}\right)$
$\displaystyle-$
$\displaystyle\Gamma_{b}\exp\left(i\varphi_{b}\right)-\Gamma_{c,2}\exp\left(i\varphi_{c,2}\right)$

and the coupling to the effective $S_{c,eff}$ is now given by the difference

$\Gamma_{c,eff}=\Gamma_{c,1}\exp\left(i\varphi_{c,1}\right)-\Gamma_{c,2}\exp\left(i\varphi_{c,2}\right)$
(135)

instead of the previous Eq. (119) for a single 0D quantum dot. Eq. (135) goes
to

$\Gamma_{c,eff}(\Phi)=\Gamma\left[1-\exp\left(i\Phi\right)\right]\exp\left(i\varphi_{c,1}\right)$
(136)

in the considered limit $\Gamma_{c,1}=\Gamma_{c,2}\equiv\Gamma$ of symmetric
couplings. Thus, the interference $|\Gamma_{c,eff}|(\Phi=0)=0$ and
$|\Gamma_{c,eff}|(\Phi=\pi)=2\Gamma$ yields inversion between $\Phi=0$ and
$\Phi=\pi$ with symmetric coupling to the superconducting leads. This
contrasts with absence of inversion for single 0D quantum dots, see section
B.1 in this Appendix.

## References

  * (1) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of Superconductivity, Phys. Rev. 108, 1175 (1957).
  * (2) P. W. Anderson, Random-Phase Approximation in the Theory of Superconductivity, Phys. Rev. 112, 1900 (1958).
  * (3) P. W. Anderson, Plasmons, Gauge Invariance, and Mass, Phys. Rev. 130, 439 (1963).
  * (4) B.D. Josephson, Possible new effects in superconductive tunnelling, Physics Letters 1, 251 (1962).
  * (5) A.F. Andreev, Thermal conductivity of the intermediate state of superconductors, Soviet Physics JETP 20, 1490 (1965)[J. Exptl. Theoret. Phys. 47, 2222 (1964)].
  * (6) P.G. de Gennes and D. Saint-James, Elementary excitations in the vicinity of a normal metal-superconducting metal contact, Phys. Lett. 4, 151 (1963).
  * (7) D. Saint-James, Excitations élémentaires au voisinage de la surface de séparation d’un métal normal et d’un métal supraconducteur J. Phys. (Paris) 25, 899 (1964).
  * (8) I.O. Kulik, Macroscopic quantization and the proximity effect in SNS junctions, Zh. Eksp. Teor. Fiz. , 57, 1745 (1970); [Sov. Phys. JETP 30, 944 (1970)].
  * (9) K.K. Likharev, Superconducting weak links, Rev. Mod. Phys. 51, 101 (1979).
  * (10) D. Averin and H.T. Imam, Supercurrent noise in quantum point contacts, Phys. Rev. Lett. 76, 3814 (1996).
  * (11) J. C. Cuevas, A. Martín-Rodero, and A. Levy Yeyati, Hamiltonian approach to the transport properties of superconducting quantum point contacts, Phys. Rev. B 54, 7366 (1996).
  * (12) E. Scheer, N. Agraït, J. C. Cuevas, A. Levy Yeyati, B. Ludophk, A. Martín-Rodero, G. Rubio Bollinger, J. M. van Ruitenbeekk, C. Urbina, The signature of chemical valence in the electrical conduction through a single-atom contact, Nature 394, 154 (1998).
  * (13) C.W. Beenakker, in Proceedings of the 14th Taniguchi International Symposium on Transport Phenomena in Mesoscopic Systems, edited by H. Fukuyama and T. Ando, Springer, Berlin (1992).
  * (14) L. Bretheau, Ç.Ö. Girit, D. Esteve, H. Pothier and C. Urbina, Tunnelling spectroscopy of Andreev states in graphene, Nature 499, 319 (2013).
  * (15) L. Bretheau, Ç.Ö. Girit, C. Urbina, D. Esteve, H. Pothier, Supercurrent Spectroscopy of Andreev States, Phys. Rev. X 3, 041034 (2013).
  * (16) D.G. Olivares, A. Levy Yeyati, L. Bretheau, Ç.Ö. Girit, H. Pothier, D. Esteve, C. Urbina, Dynamics of quasiparticle trapping in Andreev levels, Phys. Rev. B 89, 104504 (2014).
  * (17) C. Janvier, L. Tosi, L. Bretheau, Ç.Ö. Girit, M. Stern, P. Bertet, P. Joyez, D. Vion, D. Esteve, M.F. Goffman, H. Pothier, and C. Urbina, Coherent manipulation of Andreev states in superconducting atomic contacts, Science 349, 1199 (2015).
  * (18) A. Zazunov, V. S. Shumeiko, E. N. Bratus’, J. Lantz, and G. Wendin, Andreev Level Qubit, Phys. Rev. Lett. 90, 087003 (2003).
  * (19) J.M. Martinis, M. Ansmann, and J. Aumentado, Energy Decay in Superconducting Josephson-Junction Qubits from Nonequilibrium Quasiparticle Excitations, Phys. Rev. Lett. 103, 097002 (2009).
  * (20) P.J. de Visser, J.J.A. Baselmans, P. Diener, S.J.C. Yates, A. Endo, and T.M. Klapwijk, Number Fluctuations of Sparse Quasiparticles in a Superconductor, Phys. Rev. Lett. 106, 167004 (2011).
  * (21) M. Lenander, H. Wang, R.C. Bialczak, E. Lucero, M. Mariantoni, M. Neeley, A.D. O’Connell, D. Sank, M. Weides, J. Wenner, T. Yamamoto, Y. Yin, J. Zhao, A.N. Cleland, and J.M. Martinis, Measurement of energy decay in superconducting qubits from nonequilibrium quasiparticles, Phys. Rev. B 84, 024501 (2011).
  * (22) S. Rajauria, L.M.A. Pascal, Ph. Gandit, F.W.J. Hekking, B. Pannetier, and H. Courtois, Efficiency of quasiparticle evacuation in superconducting devices, Phys. Rev. B 85, 020505(R) (2012).
  * (23) J. Wenner, Y. Yin, E. Lucero, R. Barends, Y. Chen, B. Chiaro, J. Kelly, M. Lenander, M. Mariantoni, A. Megrant, C. Neill, P.J.J. O’Malley, D. Sank, A. Vainsencher, H. Wang, T.C. White, A.N. Cleland, and J.M. Martinis, Excitation of Superconducting Qubits from Hot Nonequilibrium Quasiparticles, Phys. Rev. Lett. 110, 150502 (2013).
  * (24) D. Ristè, C.C. Bultink, M.J. Tiggelman, R.N. Schouten, K.W. Lehnert, L. DiCarlo, Millisecond charge-parity fluctuations and induced decoherence in a superconducting transmon qubit, Nature Comm. 4, 1913 (2013).
  * (25) E.M. Levenson-Falk, F. Kos, R. Vijay, L. Glazman, and I. Siddiqi, Single-Quasiparticle Trapping in Aluminum Nanobridge Josephson Junctions, Phys. Rev. Lett. 112, 047002 (2014).
  * (26) J.S. Meyer, M. Houzet and A.V. Nazarov, Dynamical spin polarization of excess quasi-particles in superconductors, Phys. Rev. Lett. 125, 097006 (2020).
  * (27) D. Kouznetsov, D. Rohrlich, and R. Ortega, Quantum limit of noise of a phase-invariant amplifier, Phys. Rev. A 52, 1665 (1995).
  * (28) The SQUID handbook Vol. I Fundamentals and Technology of SQUIDs and SQUID Systems, J. Clarke, and A.I. Braginski (Eds.), Wiley-Vch (2004).
  * (29) J. Clarke, and F.K. Wilhelm, Superconducting quantum bits, Nature 453, 1031 (2008).
  * (30) M.H. Devoret and R.J. Schoelkopf, Superconducting Circuits for Quantum Information: An Outlook, Science 339, 1169 (2013).
  * (31) F. Arute et al., Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * (32) S.B. Kaplan, C.C. Chi, D.N. Langenberg, J.J. Chang, S. Jafarey, and D.J. Scalapino, Quasiparticle and phonon lifetimes in superconductors, Phys. Rev. B 14, 4854 (1976).
  * (33) R.C. Dynes, V. Narayanamurti, and J.P. Garno, Direct measurement of quasiparticle-lifetime broadening in a strong-coupled superconductor, Phys. Rev. Lett. 41, 1509 (1978).
  * (34) J.P. Pekola, V.F. Maisi, S. Kafanov, N. Chekurov, A. Kemppinen, Yu.A. Pashkin, O.-P. Saira, M. Möttönen, and J.S. Tsai, Environment-assisted tunneling as an origin of the Dynes density of states, Phys. Rev. Lett. 105, 026803 (2010).
  * (35) O.-P. Saira, A. Kemppinen, V.F. Maisi, and J.P. Pekola, Vanishing quasiparticle density in a hybrid Al/Cu/Al single-electron transistor, Phys. Rev. B 85, 012504 (2012).
  * (36) T. Meng, S. Florens and P. Simon, Self-consistent description of Andreev bound states in Josephson quantum dot devices, Phys. Rev. B 79, 224521 (2009).
  * (37) R. L. Klees, G. Rastelli, J. C. Cuevas, and W. Belzig, Microwave Spectroscopy Reveals the Quantum Geometric Tensor of Topological Josephson Matter, Phys. Rev. Lett. 124, 197002 (2020).
  * (38) N. K. Allsopp, V. C. Hui, C. J. Lambert, and S. J. Robinson, Theory of the sign of multi-probe conductances for normal and superconducting materials, J. Phys.: Condens. Matter 6, 10475 (1994).
  * (39) J. M. Byers and M. E. Flatté, Probing Spatial Correlations with Nanoscale Two-Contact Tunneling, Phys. Rev. Lett. 74, 306 (1995).
  * (40) J. Torrès and T. Martin, Positive and negative Hanbury-Brown and Twiss correlations in normal metal-superconducting devices, Eur. Phys. J. B 12, 319 (1999).
  * (41) G. Deutscher and D. Feinberg, Coupling superconducting-ferromagnetic point contacts by Andreev reflections, Appl. Phys. Lett. 76, 487 (2000).
  * (42) M. S. Choi, C. Bruder, and D. Loss, Spin-dependent Josephson current through double quantum dots and measurement of entangled electron states, Phys. Rev. B 62, 13569 (2000).
  * (43) G. Falci, D. Feinberg, and F. W. J. Hekking, Correlated tunneling into a superconductor in a multiprobe hybrid structure, Europhys. Lett. 54, 255 (2001).
  * (44) P. Recher, E. V. Sukhorukov, and D. Loss, Andreev tunneling, Coulomb blockade, and resonant transport of nonlocal spin-entangled electrons, Phys. Rev. B 63, 165314 (2001).
  * (45) G. B. Lesovik, T. Martin, and G. Blatter, Electronic entanglement in the vicinity of a superconductor, Eur. Phys. J. B 24, 287 (2001).
  * (46) R. Mélin and D. Feinberg, Transport theory of multiterminal hybrid structures, Eur. Phys. J. B 26, 101 (2002).
  * (47) N. M. Chtchelkatchev, G. Blatter, G. B. Lesovik, and T. Martin, Bell inequalities and entanglement in solid-state devices, Phys. Rev. B 66, 161320 (2002).
  * (48) R. Mélin and D. Feinberg, Sign of the crossed conductances at a ferromagnet/superconductor/ferromagnet double interface, Phys. Rev. B 70, 174509 (2004).
  * (49) A. V. Lebedev, G. B. Lesovik, and G. Blatter, Generating spin-entangled electron pairs in normal conductors using voltage pulses, Phys. Rev. B 72, 245314 (2005).
  * (50) K. V. Bayandin, G. B. Lesovik, and T. Martin, Energy entanglement in normal metal–superconducting forks Phys. Rev. B 74, 085326 (2006).
  * (51) A. L. Yeyati, F. S. Bergeret, A. Martín-Rodero, and T. M. Klapwijk, Entangled Andreev pairs and collective excitations in nanoscale superconductors, Nature Phys. 3, 455 (2007).
  * (52) D. Beckmann, H. B. Weber, and H. v. Löhneysen, Evidence for crossed Andreev reflection in Superconductor-Ferromagnet hybrid structures, Phys. Rev. Lett. 93, 197003 (2004).
  * (53) S. Russo, M. Kroug, T. M. Klapwijk, and A. F. Morpurgo, Experimental observation of bias-dependent nonlocal Andreev reflection, Phys. Rev. Lett. 95, 027002 (2005).
  * (54) P. Cadden-Zimansky and V. Chandrasekhar, Nonlocal correlations in normal-metal superconducting systems, Phys. Rev. Lett. 97, 237003 (2006).
  * (55) D. Beckmann, H. B.Weber, and H. v. Löhneysen, Negative four-terminal resistance as a probe of crossed Andreev reflection, Appl. Phys. A 89, 603 (2007).
  * (56) P. Cadden-Zimansky, Z. Jiang, and V. Chandrasekhar, Charge imbalance, crossed Andreev reflection and elastic co-tunnelling in ferromagnet/superconductor/normal-metal structures, New J. Phys. 9, 116 (2007).
  * (57) L. Hofstetter, S. Csonka, J. Nygard, and C. Schönenberger, Cooper pair splitter realized in a two-quantum-dot Y-junction, Nature (London) 461, 960 (2009).
  * (58) L. G. Herrmann, F. Portier, P. Roche, A. Levy Yeyati, T. Kontos, and C. Strunk, Carbon nanotubes as Cooper pair beam splitters, Phys. Rev. Lett. 104, 026801 (2010).
  * (59) J. Wei and V. Chandrasekhar, Positive noise cross-correlation in hybrid superconducting and normal-metal three-terminal devices, Nat. Phys. 6, 494 (2010).
  * (60) A. Das, Y. Ronen, M. Heiblum, D. Mahalu, A. V. Kretinin, and H. Shtrikman, High-efficiency Cooper pair splitting demonstrated by two-particle conductance resonance and positive noise cross- correlation, Nat. Commun. 3, 1165 (2012).
  * (61) J. Schindele, A. Baumgartner, and C. Schönenberger, Near-Unity Cooper Pair Splitting Efficiency, Phys. Rev. Lett. 109, 157002 (2012).
  * (62) J. Schindele, A. Baumgartner, R. Maurand, M. Weiss, and C. Schönenberger, Nonlocal spectroscopy of Andreev bound states, Phys. Rev. B 89, 045422 (2014).
  * (63) Z. B. Tan, D. Cox, T. Nieminen, P. Lähteenmäki, D. Golubev, G. B. Lesovik, and P. J. Hakonen, Cooper Pair Splitting by Means of Graphene Quantum Dots, Phys. Rev. Lett. 114, 096602 (2015).
  * (64) I. V. Borzenets, Y. Shimazaki, G. F. Jones, M. F. Craciun, S. Russo, M. Yamamoto, and S. Tarucha, High Efficiency CVD Graphene-lead (Pb) Cooper Pair Splitter, Sci. Rep. 6, 23051 (2016).
  * (65) P. Pandey, R. Danneau, and D. Beckmann, Ballistic Graphene Cooper Pair Splitter, Phys. Rev. Lett. 126, 147701 (2021).
  * (66) M. P. Anantram and S. Datta, Current fluctuations in mesoscopic systems with Andreev scattering, Phys. Rev. B 53, 16390 (1996).
  * (67) P. Samuelsson and M. Büttiker, Chaotic dot-superconductor analog of the Hanbury Brown–Twiss effect, Phys. Rev. Lett. 89,046601 (2002).
  * (68) P. Samuelsson and M. Büttiker, Semiclassical theory of current correlations in chaotic dot-superconductor systems, Phys. Rev. B 66, 201306 (R) (2002).
  * (69) J. Börlin, W. Belzig, and C. Bruder, Full counting statistics of a superconducting beam splitter, Phys. Rev. Lett. 88,197001 (2002).
  * (70) P. Samuelsson, E.V. Sukhorukov,and M. Büttiker, Orbital entanglement and violation of Bell inequalities in mesoscopic conductors, Phys. Rev. Lett. 91, 157002 (2003).
  * (71) L. Faoro, F. Taddei, and R. Fazio, Clauser-Horne inequality for electron-counting statistics in multiterminal mesoscopic conductors, Phys. Rev. B 69, 125326 (2004).
  * (72) G. Bignon, M. Houzet, F. Pistolesi, and F.W.J. Hekking, Current-current correlations in hybrid superconducting and normal metal multiterminal structures, Europhys. Lett. 67, 110 (2004).
  * (73) R. Mélin, C. Benjamin, and T. Martin, Positive cross correlations of noise in superconducting hybrid structures: Roles of interfaces and interactions, Phys. Rev. B 77, 094512 (2008).
  * (74) A. Freyn, M. Flöser and R. Mélin, Positive current cross-correlations in a highly transparent normal-superconducting beam splitter due to synchronized Andreev and inverse Andreev reflections, Phys. Rev. B 82, 014510 (2010).
  * (75) D.S. Golubev and A.D. Zaikin, Shot noise and Coulomb effects on nonlocal electron transport in normal-metal/superconductor/normal-metal heterostructures, Phys. Rev. B 82, 134508 (2010).
  * (76) M. Flöser, D. Feinberg and R. Mélin, Absence of split pairs in cross correlations of a highly transparent normal metal–superconductor–normal metal electron-beam splitter, Phys. Rev. B 88, 094517 (2013).
  * (77) G. Michalek, B. R. Bulka, T. Domański, and K. I. Wysokiński, Statistical correlations of currents flowing through a proximitized quantum dot, Phys. Rev. B 101, 235402 (2020).
  * (78) J. C. Cuevas and H. Pothier, Voltage-induced Shapiro steps in a superconducting multiterminal structure, Phys. Rev. B 75, 174513 (2007).
  * (79) A. Freyn, B. Douçot, D. Feinberg, and R. Mélin, Production of non-local quartets and phase-sensitive entanglement in a superconducting beam splitter, Phys. Rev. Lett. 106, 257005 (2011).
  * (80) T. Jonckheere, J. Rech, T. Martin, B. Douçot, D. Feinberg, and R. Mélin, Multipair DC Josephson resonances in a biased allsuperconducting bijunction, Phys. Rev. B 87, 214501 (2013).
  * (81) R. Mélin, D. Feinberg, and B. Douçot, Partially resummed perturbation theory for multiple Andreev reflections in a short three-terminal Josephson junction, Eur. Phys. J. B 89, 67 (2016).
  * (82) R. Mélin, M. Sotto, D. Feinberg, J.-G. Caputo and B. Douçot, Gate-tunable zero-frequency current cross-correlations of the quartet mode in a voltage-biased three-terminal Josephson junction, Phys. Rev. B 93, 115436 (2016).
  * (83) R. Mélin, J.-G. Caputo, K. Yang and B. Douçot, Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction, Phys. Rev. B 95, 085415 (2017).
  * (84) R. Mélin, R. Danneau, K. Yang, J.-G. Caputo, and B. Douçot, Engineering the Floquet spectrum of superconducting multiterminal quantum dots, Phys. Rev. B 100, 035450 (2019).
  * (85) J.D. Pillet, V. Benzoni, J. Griesmar, J.-L. Smirr, and Ç.Ö. Girit, Nonlocal Josephson Effect in Andreev Molecules Nano Lett. 19, 7138 (2019).
  * (86) V. Kornich, H.S. Barakov, and Yu.V. Nazarov, Fine energy splitting of overlapping Andreev bound states in multiterminal superconducting nanostructures, Phys. Rev. Research 1, 033004 (2019).
  * (87) R. Mélin and B. Douçot, Inversion in a four terminal superconducting device on the quartet line. I. Two-dimensional metal and the quartet beam splitter, Phys. Rev. B 102, 245435 (2020).
  * (88) R. Mélin and B. Douçot, Inversion in a four terminal superconducting device on the quartet line. II. Quantum dot and Floquet theory, Phys. Rev. B 102, 245436 (2020).
  * (89) J.-D. Pillet, V. Benzoni, J. Griesmar, J.-L. Smirr, and Ç Ö Girit, Scattering description of Andreev molecules, SciPost Phys. Core 2, 009 (2020).
  * (90) V. Kornich, H. S. Barakov and Yu. V. Nazarov, Overlapping Andreev states in semiconducting nanowires: competition of 1D and 3D propagation, Phys. Rev. B 101, 195430 (2020).
  * (91) R. Mélin, Ultralong-distance quantum correlations in three-terminal Josephson junctions, Phys. Rev. B 104, 075402 (2021).
  * (92) A. Melo, V. Fatemi and A.R. Akhmerov, Multiplet supercurrent in Josephson tunneling circuits, arXiv:2104.11239 (2021).
  * (93) A.H. Pfeffer, J.E. Duvauchelle, H. Courtois, R. Mélin, D. Feinberg, and F. Lefloch, Subgap structure in the conductance of a three-terminal Josephson junction, Phys. Rev. B 90, 075401 (2014).
  * (94) E. Strambini, S. D’Ambrosio, F. Vischi, F.S. Bergeret, Yu.V. Nazarov, and F. Giazotto, The $\omega$-SQUIPT as a tool to phase-engineer Josephson topological materials, Nat. Nanotechnol. 11, 1055 (2016).
  * (95) Y. Cohen, Y. Ronen, J.H. Kang, M. Heiblum, D. Feinberg, R. Mélin, and H. Strikman, Non-local supercurrent of quartets in a three-terminal Josephson junction, Proc. Natl. Acad. Sci. U. S. A. 115, 6991 (2018).
  * (96) A.W. Draelos, M.-T. Wei, A. Seredinski, H. Li, Y. Mehta, K. Watanabe, T. Taniguchi, I.V. Borzenets, F. Amet, and G. Finkelstein, Supercurrent flow in multiterminal graphene Josephson junctions, Nano Lett. 19, 1039 (2019).
  * (97) K.F. Huang, Y. Ronen, R. Mélin, D. Feinberg, K. Watanabe, T. Taniguchi, and P. Kim, Quartet supercurrent in a multi-terminal Graphene-based Josephson Junction, arXiv:2008.03419 (2020).
  * (98) N. Pankratova, H. Lee, R. Kuzmin, K. Wickramasinghe,1, W. Mayer,J. Yuan,M. Vavilov,J. Shabani and V. Manucharyan, The multi-terminal Josephson effect, Phys. Rev. X 10, 031051 (2020).
  * (99) G.V. Graziano, J.S. Lee, M. Pendharkar, C. Palmstrom and V.S. Pribiag, Transport Studies in a Gate-Tunable Three-Terminal Josephson Junction, Phys. Rev. B 101, 054510 (2020).
  * (100) E.G. Arnault, T. Larson, A. Seredinski, L. Zhao, H. Li, K. Watanabe, T. Tanniguchi, I. Borzenets, F. Amet and G. Finkelstein, The multiterminal inverse AC Josephson effect, arXiv:2012.15253v1 (2020).
  * (101) S.A. Khan, L. Stampfer, T. Mutas, J.-H. Kang, P. Krogstrup and T.S. Jespersen, Multiterminal Quantized Conductance in InSb Nanocrosses, arXiv:2101.02529 (2021).
  * (102) O. Kürtössy, Z. Scherübl, G. Fülöp, I. E. Lukács, T. Kanne, J. Nygard, P. Makk and S. Csonka, Andreev molecule in parallel InAs nanowires, arXiv:2103.14083 (2021).
  * (103) Gino V. Graziano, Mohit Gupta, Mihir Pendharkar, Jason T. Dong, Connor P. Dempsey, Chris Palmstrom, Vlad S. Pribiag, Selective Control of Conductance Modes in Multi-terminal Josephson Junctions, arXiv:2201.01373 (2022).
  * (104) E.G. Arnault, S. Idris, A. McConnell, L. Zhao, T.F.Q. Larson, K. Watanabe, T. Taniguchi, G. Finkelstein, F. Amet, Dynamical Stabilization of Multiplet Supercurrents in Multi-terminal Josephson Junctions, arXiv:2201.11180 (2022).
  * (105) B. van Heck, S. Mi, and A.R. Akhmerov, Single fermion manipulation via superconducting phase differences in multiterminal Josephson junctions, Phys. Rev. B 90, 155450 (2014).
  * (106) C. Padurariu, T. Jonckheere, J. Rech, R. Mélin, D. Feinberg, T. Martin, and Yu.V. Nazarov, Closing the proximity gap in a metallic Josephson junction between three superconductors, Phys. Rev. B 92, 205409 (2015).
  * (107) R.-P. Riwar, M. Houzet, J.S. Meyer, and Y.V. Nazarov, Multi-terminal Josephson junctions as topological materials, Nat. Commun. 7, 11167 (2016).
  * (108) E. Eriksson, R.-P. Riwar, M. Houzet, J. S. Meyer, and Y. V. Nazarov, Topological transconductance quantization in a four-terminal Josephson junction, Phys. Rev. B 95, 075417 (2017).
  * (109) O. Deb, K. Sengupta and D. Sen, Josephson junctions of multiple superconducting wires, Phys. Rev. B 97, 174518 (2018).
  * (110) V. Fatem, A.R. Akhmerov and L. Bretheau, Weyl Josepshon circuits, arXiv:2008.13758v1 (2020).
  * (111) L. Peyruchat, J. Griesmar, J.-D. Pillet and Ç.Ö Girit, Transconductance quantization in a topological Josephson tunnel junction circuit, arXiv:2009.03291v1 (2020).
  * (112) H. Weisbrich, R.L. Klees, G. Rastelli and W. Belzig, Second Chern Number and Non-Abelian Berry Phase in Topological Superconducting Systems, PRX Quantum 2, 010310 (2021).
  * (113) Y. Chen and Y.V. Nazarov, Weyl point immersed in a continuous spectrum: an example from superconducting nanostructures, Phys. Rev. B 104, 104506 (2021).
  * (114) Y. Chen and Y.V. Nazarov, Spin-Weyl quantum unit: theoretical proposal, Phys. Rev. B 103, 045410 (2021).
  * (115) E.V. Repin and Y.V. Nazarov, Weyl points in the multi-terminal Hybrid Superconductor-Semiconductor Nanowire devices, arXiv:2010.11494v1 (2020).
  * (116) B. Douçot, R. Danneau, K. Yang, J.-G. Caputo and R. Mélin, Berry phase in superconducting multiterminal quantum dots, Phys. Rev. B 101, 035411 (2020).
  * (117) B. Venitucci, D. Feinberg, R. Mélin, B. Douçot, Nonadiabatic Josephson current pumping by microwave irradiation, Phys. Rev. B 97, 195423 (2018).
  * (118) L.P. Gavensky, G. Usaj, D. Feinberg and C.A. Balseiro, Berry curvature tomography and realization of topological Haldane model in driven three-terminal Josephson junctions, Phys. Rev. B 97, 220505 (2018).
  * (119) H.-Y. Xie, M.G. Vavilov and A. Levchenko, Topological Andreev bands in three-terminal Josephson junctions, Phys. Rev. B 96, 161406 (2017).
  * (120) H.-Y. Xie, M.G. Vavilov and A. Levchenko, Weyl nodes in Andreev spectra of multiterminal Josephson junctions: Chern numbers, conductances and supercurrents, Phys. Rev. B 97, 035443 (2018).
  * (121) S. Park, W. Lee, S. Jang, Y.-B. Choi, J. Park, W. Jung, K. Watanabe, T. Taniguchi, G. Y. Cho and G.-H. Lee, Steady Floquet-Andreev states in graphene Josephson junctions, Nature 603, 421 (2022).
  * (122) F. S. Bergeret, P. Virtanen, A. Ozaeta, T. T. Heikkilä, and J. C. Cuevas Phys. Rev. B 84, 054504 (2011).
  * (123) A. Martín-Rodero, A. Levy Yeyati, F.J. García-Vidal, Thermal noise in superconducting quantum point contacts, Phys. Rev. B 53, R8891 (1996).
  * (124) C. Caroli, R. Combescot, P. Nozières and D. Saint-James, Direct calculation of the tunneling current, Jour. Phys. C: Solid State Phys. 4, 916 (1971).
  * (125) Yu.V. Nazarov and Y.M. Blanter, Quantum transport, Cambridge University Press (2012).