# Quartet currents in a biased three-terminal diffusive Josephson junction

T. Jonckheere1 J. Rech1 C. Padurariu2 L. Raymond1 T. Martin1 D. Feinberg3,4 1
Aix Marseille Univ., Université de Toulon, CNRS, CPT, IphU, AMUtech,
Marseille, France  2 Institut für Komplexe Quantensysteme and IQST,
Universität Ulm, Albert-Einstein-Allee 11, 89069 Ulm, Germany  3 Centre
National de la Recherche Scientifique, Institut NEEL, F-38042 Grenoble Cedex
9, France  4 Université Grenoble-Alpes, Institut NEEL, F-38042 Grenoble Cedex
9, France

###### Abstract

Biasing a three-terminal Josephson junction (TTJ) with symmetrical voltages
$0,V,-V$ leads to new kinds of DC currents, namely quartet Josephson currents
and phase-dependent multiple Andreev reflection (MAR) currents. We study these
currents in a system where a normal diffusive metallic node $N$ is connected
to three terminals $S_{0,1,2}$ by barriers of arbitrary transparency. We use
the quantum circuit theory to calculate the current in each terminal,
including decoherence. In addition to the stationary combination
$\varphi_{Q}=\varphi_{1}+\varphi_{2}-2\varphi_{0}$ of the terminal phases
$\varphi_{i}$, the bias voltage $V$ appears as a new and unusual control
variable for a DC Josephson current. A general feature is the sign changes of
the current-phase characteristics, manifesting in minima of the quartet
“critical current”. Those sign changes can be triggered by the voltage, by the
junction transparency or by decoherence. We study the possible separation of
quartet currents from MAR currents in different regimes of parameters,
including an ”funnel” regime with very asymmetric couplings to $S_{0,1,2}$. In
the regime of low transparency and asymetric couplings, we provide an analytic
perturbative expression for the currents which shows an excellent agreement
with the full numerical results.

###### pacs:

73.23.-b, 73.63.Kv 74.45.+c

##  I Introduction

The understanding of subgap transport in transparent Superconducting/Normal
Metal/Superconducting ($SNS$) Josephson junctions, either point contact-like,
or made of a diffusive metallic region, has reached an advanced state on both
experimental and theoretical levels. Multiple Andreev reflections (MARs) allow
to describe dissipative quasiparticle motion below the superconducting gap
$\Delta$, where transport is assisted by Cooper pairs. MARs  Electrons and
holes cross the structure, being Andreev-reflected at each $NS$ interface.
Each crossing provides the energy $eV$ for electrons and holes, in such a way
that steps in the $I(V)$ characteristics appear around values of the voltage
$V\equiv 2\Delta/ne$ ($n$ positive integer). The most quantitative
understanding of MAR transport has been obtained for point contact junctions,
where transport is fully phase-coherent. QPC  On the other hand, in diffusive
$SNS$ junctions, another important parameter appears due to the dephasing of
electrons and holes between successive Andreev reflections. This leads to the
so-called Thouless energy $E_{Th}$ as another relevant energy scale. Thouless
$E_{Th}$ is proportional to the inverse of the dwell time $\tau_{d}$ for a
quasiparticle diffusing across the structure. In transparent $SNS$ junctions,
shorter than the superconducting coherence length, one has $\Delta<E_{Th}$,
while $E_{Th}<\Delta$ on the contrary in long diffusive junctions. MARs in the
incoherent regime have been extensively studied Bezuglyi2000  and reproducible
junctions with very good quality are routinely fabricated. A minigap forms at
an energy scale set by $E_{Th}$, which also sets the scale of the junction
critical current, e.g. $eI_{c}\simeq G_{N}E_{Th}$ where $G_{N}$ is the
junction normal conductance.

Recently, multiterminal structures, made of three or more superconducting
terminals contacted by a single normal region, have been considered
theoretically, both at equilibrium EPJB2015 ; Rech ; Mi-Akhmerov ; MelinMAR ;
Riwar ; Meyer ; Pillet  and under voltage bias. Cuevas ; Houzet ; Freyn ;
JonckheerePRB2013 ; AkhmerovTTJ ; KeliriDoucot  Experiments have addressed
either three terminals but two being equipotential, Giazotto  or three
independent terminals i.e. biased with different voltages, the latter
corresponding to the present study. Experiments involve junctions made of a
diffusive metal, Pfeffer  clean nanowires Ronen ; Khan  or 2D electron gas
Pankratova ; Graziano1 ; Graziano2  and multichannel graphene. Draelos ;
Harvard ; Finkelstein21 ; Finkelstein22 ; Zhang  Transport in multiterminal
junctions enables DC phase-coherent transport at nonzero but commensurate
voltages. Cuevas ; Houzet ; Freyn  Indeed, setting to zero both voltage and
phase at, say, terminal $S_{0}$, and biasing terminals $S_{1}$ and $S_{2}$ at
voltages $V_{1}$ and $V_{2}$, due to the Josephson relation
$\frac{d\varphi_{i}}{dt}=\frac{2eV_{i}}{\hbar}$, two relevant variables can be
chosen within the sets $(\varphi_{1},\varphi_{2})$ and $(V_{1},V_{2})$. For
instance, fixing $V_{1}+V_{2}=0$ leaves $V_{1}-V_{2}=2V$ and
$\varphi_{1}+\varphi_{2}$ as independent control variables. This holds for any
commensurate combination $nV_{1}+mV_{2}=0$, $n\varphi_{1}+m\varphi_{2}$ being
the relevant stationary phase variable.

As a result, a stationary current component can appear as a DC resonance,
similar to an equilibrium Josephson current, but also controlled by the
voltage V. For $V_{1}+V_{2}=0$ such a current involves two entangled pairs
emitted from $S_{0}$ crossing simultaneously and transmitted, one in $S_{1}$,
the other in $S_{2}$. Those processes have been dubbed as “quartet currents”
since the elementary “exchange currency” between one terminal and the two
others is made of two Cooper pairs instead of one. Of course, they coexist
with the standard AC Josephson current flowing between any pair of terminals.

Another striking signature of the voltage commensurability is the presence of
phase-dependent ”phase-MAR” currents, where quasiparticle transport is
accompanied by quartets instead of single pairs like in the two-terminal case.
JonckheerePRB2013 ; Padurariu2 ; QuartetSQUID  Contrarily to two-terminal MARs
flowing between any pair of the terminals biased at $(0,V,-V)$, phase-MARs
involve all three terminals like quartets.

Such quartet currents, and higher order resonances, are reminiscent of the
synchronization of AC Josephson currents in junction arrays, well documented
some decades ago, Likharev  and recently revisited in a more general context.
AkhmerovMultipair  In these downmixing phenomena, leading to DC currents in
biased arrays, the underlying mechanism can be understood by the classical
nonlinear dynamics of an equivalent RSJ (resistively shunted junction)
circuit, and it is strongly conditioned by the circuit environment, i.e. by
the impedance of the junction plus the external circuit at the AC Josephson
frequency. In contrast, the phenomena addressed in this work are intrinsically
quantum and fully mesoscopic. Just as MARs, they cannot be described by
classical or quantum phase dynamics in a Josephson junction array, using a
given current-phase characteristics for each junction. Likharev ;
AkhmerovMultipair  Quartet and phase-MAR currents are indeed another
consequence of multiple Andreev reflections, allowed only in a multiterminal
structure. They open new phase-coherent channels in a dissipative structure,
that cannot be deduced in a phenomenological way from independent two-terminal
Josephson currents. Discriminating classical downmixing from mesoscopic multi-
pair resonances and separating quartets from phase-MARs remains an
experimental challenge for the future, that could be tackled with a suitable
SQUID geometry. QuartetSQUID

Experiments in TTJ’s Pfeffer ; Ronen ; Harvard  have revealed clear signatures
of Josephson-like currents in a low-bias regime ($eV\ll\Delta$), beyond the
regime of DC Josephson effects between any pair of terminals, i.e. for
currents higher than the critical currents flowing in any pair of terminals.
In this low-bias range, MAR currents are strongly hampered by a necessarily
large number of Andreev reflections. On the contrary, quartets need minimally
four Andreev reflections, which makes the quartet current visible against
MARs. In fact, it has been shown theoretically that the quartet current does
not decrease like MAR currents when $eV$ goes to $0$ but instead enters a
regime where the current sign changes several times with voltage, due to
nonadiabatic resonances. Regis2 ; Regis3  In Refs. Pfeffer, ; Ronen, ;
Harvard, , robust transport anomalies were found as a function of the two
applied voltages $V_{1,2}$, and were interpreted by three “quartet” modes
corresponding to two Cooper pairs simultaneously crossing the structure. In
Ref. Pfeffer,  the junctions are long and the Thouless energy is much smaller
than the superconducting gap. The observation of strong anomalies in a regime
$E_{Th}<V<\Delta$, where the voltage is too large to allow coherent motion of
individual pairs, is especially striking. In Refs. Ronen, ; Harvard,  instead,
the system is in a coherent regime and MARs can be observed at higher
voltages.

A recent experiment used a four-terminal junction where two of them are
grounded and enclose a loop pierced by a flux. Harvard  The flux dependence of
the quartet current displayed the expected oscillating behavior, but with two
new features : the presence of the $\frac{hc}{4e}$ periodicity signaling the
quartet mode, and a phase inversion of the oscillation. The latter was
ascribed to resonances within the effective Andreev spectrum, due to the
underlying running phase originating from the voltage bias. Regis1 ; Regis2 ;
Regis3  This sheds light on the microscopic origin of quartets, viewed as
”self-induced Shapiro steps” in Ref. Cuevas, . Theory has, until now,
addressed the case of a few channels, Regis1 ; Regis2 ; Regis3  describing the
junction with a few-level quantum dot. MARs in a TTJ have also been described
with the help of scattering theory, in the full $(V_{1},V_{2})$ plane.
AkhmerovTTJ  Yet a microscopic description is still lacking for a diffusive
junction with a very large number of channels.

Motivated by recent experiments where a sizable normal metal island (rather
than a system of quantum dots) separates the superconducting leads, we
consider here a short diffusive three-terminal junction (TTJ) (Fig. 1)
described by quantum circuit theory Yuli ; Yulibook  to calculate quartet and
MAR DC currents at any voltage. Dephasing can be introduced and may lead to a
small Thouless energy, mimicking one important aspect of a long diffusive
junction. The dependence of currents with voltage and quartet phase is studied
in detail, both for a symmetric and a strongly asymmetric TTJ. This work
follows previous ones where the spectral normal density of states NDOS was
calculated, first at equilibrium, Padurariu1  then at nonzero voltages
Padurariu2  at which BCS-like resonances appear at MAR voltages. Moreover, the
low-voltage behavior of the NDOS of a symmetric TTJ displays a phase-dependent
pseudogap behavior, allowing to make the link with an adiabatic regime when
$V\rightarrow 0$. In the present work, the spectral quartet current is shown
to possess resonances at voltages coinciding with those of the NDOS, but with
signs alternating between successive resonances. Integration on energy yields
the low-voltage phase-coherent DC current, which in general departs strongly
from an adiabatic behavior, with strong oscillations in the voltage
dependence.

This motivates the study of an asymmetric ”funnel” TTJ where the biased
terminals are much more weakly coupled than the unbiased one (Fig. 1). In this
situation, we show that it is possible to reach a quasi-adiabatic regime,
where the low-voltage current does not exhibit sign changes anymore. The
numerical study presented here is complemented by a perturbative expansion of
the nonlinear circuit theory equations.

The other focus of the paper is the balance between quartet and MAR currents,
important in view of obtaining unambiguous signatures of quartet current in
experiments. As expected, even at low voltages, both have the same order of
magnitude for a symmetric TTJ corresponding to a fully resonant cavity making
the junction. On the contrary, in the funnel case, the phase-dependent MAR
component becomes negligible at low voltages. Yet the phase-independent MAR
remains strong if the junction formed by the two biased terminals (at $V,-V$)
is symmetric. Only in a completely asymmetric TTJ, all MAR processes can be
made smaller than quartets at low voltage. Finally, owing to the odd (for
quartets) and even (for MARs) symmetries in phase, the DC current in all
branches, which combines quartet and MAR contributions, yields a phase shift
in the current-phase relationship and therefore makes TTJs a new source of
tunable ”$\varphi$-junction”.

Figure 1: Sketch of a three-terminal junction (TTJ), with superconductors
$S_{0,1,2}$ at voltages $(0,V,-V)$ respectively, coupled to a node by
couplings $\Gamma_{0,1,2}$ (see text). The thickness of the connection between
the leads and the central node represents the strength of the corresponding
coupling contant. (a) symmetric TTJ ($\Gamma_{0}=\Gamma_{1}=\Gamma_{2}$); (b)
funnel TTJ ($\Gamma_{0}\gg\Gamma_{1}=\Gamma_{2}$); (c) fully asymmetric TTJ
($\Gamma_{0}\gg\Gamma_{1}\gg\Gamma_{2}$).

The paper is structured as follows. The model and basic equations are
presented in Section II. Section III discusses how the symmetry with respect
to voltage inversion allows a practical definition of quartet and MAR currents
for a junction with general asymmetry of the contacts. Section IV addresses a
perturbative expansion when the contacts are strongly asymmetric, in the
”funnel” case. Section V provides results on the spectral currents, that can
be compared to the NDOS structure. Section VI describes the phase and voltage
dependence of the quartet currents, for a symmetric and then for an asymmetric
TTJ. Section VII discusses MAR currents and compares them to quartet currents.
Section VIII provides a conclusion and perspectives.

##  II The model

Let us consider a three-terminal junction formed by connecting three
superconductors with a normal diffusive metallic region (Fig. 1). The
superconducting coherence length being larger than both the Fermi wavelength
and the elastic mean-free-path, one can employ the quasiclassical equations of
nonequilibrium superconductivity, which take the form of a diffusive equation
for the quasiclassical Keldysh-Nambu Green’s function: Larkin

$\check{G}=\left(\begin{array}[]{cc}\hat{G}^{R}&\hat{G}^{K}\\\
0&\hat{G}^{A}\end{array}\right)\ ,\quad\check{G}^{2}=\check{1},$ (1)

where $\hat{G}^{R,A,K}$ correspond to the retarded, advanced and Keldysh
components respectively. In addition to Keldysh-Nambu space (to be denoted
with a check hat, $\check{G}$), the quasiclassical Green’s function generally
depends on two times (or energies) and on spatial coordinates,
$\check{G}(\textbf{E},\textbf{x})$. The diffusive equation for the
corresponding Green’s function is also known as the Usadel equation (see also
Ref. Yulibook, ).

$\displaystyle\frac{\partial}{\partial\textbf{x}}\left({\cal
D}(\textbf{x})\check{G}\frac{\partial}{\partial\textbf{x}}\check{G}\right)-i\left[\check{H},\check{G}\right]=0,\quad\check{H}=\left(\begin{array}[]{cc}\hat{H}&0\\\
0&\hat{H}\end{array}\right);$ (4)
$\displaystyle\hat{H}=E\hat{\sigma}_{z}+\frac{1}{2}\Delta(\textbf{x})(i\hat{\sigma}_{y}+\hat{\sigma}_{x})+\frac{1}{2}\Delta^{*}(\textbf{x})(i\hat{\sigma}_{y}-\hat{\sigma}_{x}).$

Here the products should be understood as a matrix product in Nambu-Keldysh
space along with a convolution in the double time (or energy) representation.
The Pauli matrices in Nambu space are here denoted with a hat by
$\vec{\hat{\sigma}}={\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}}$, and
${\cal D}(\textbf{x})$ denotes the diffusion coefficient.

The Usadel equation can be seen as a conservation of the quasiparticle current
density, $\check{j}(\textbf{x})$.

$\displaystyle\frac{\partial}{\partial\textbf{x}}\check{j}(\textbf{x})+\frac{ie^{2}\nu}{\hbar}\left[\check{H},\check{G}\right]=0;\quad\check{j}=-\sigma(\textbf{x})\check{G}\frac{\partial}{\partial\textbf{x}}\check{G}.$
(5)

Here, $\nu$ is the electronic density of states and $\sigma(\textbf{x})$ is
the conductivity. The two quantities are related by $\sigma=e^{2}{\cal D}\nu$.

Here we employ a discretized version of the quasiclassical equations, the so-
called quantum circuit theory. Yuli ; Yulibook  The superconducting terminals
$S_{i}$ are assumed to be homogeneous and described by zero-dimensional
Keldysh-Nambu matrices $\check{G}_{i}$, and they are connected by a zero-
dimensional “node” described by the unknown matrix $\check{G}_{c}$. At
equilibrium, the terminal Green’s functions are functions of one energy
variable and are given by:

$\displaystyle\hat{G}^{R}_{i}=\frac{1}{\xi}\left(\begin{array}[]{cc}\epsilon&\Delta_{i}\\\
-\Delta_{i}^{*}&-\epsilon\end{array}\right)\
;\quad\hat{G}^{A}_{i}=-\frac{1}{\xi^{*}}\left(\begin{array}[]{cc}\epsilon^{*}&\Delta_{i}\\\
-\Delta_{i}^{*}&-\epsilon^{*}\end{array}\right)\ ,$ (10)

where complex energies have been introduced as $\epsilon=E+i0^{+}$ and
$\xi=\sqrt{\epsilon^{2}-|\Delta|^{2}}$. Here
$\Delta_{i}=|\Delta|e^{i\varphi_{i}}$. The positive vanishing complex part of
$\epsilon$ is essential in view of the branch cut of the square root function
in the complex plane.

The advanced and retarded Green’s functions are related by
$\hat{G}_{i}^{A}=-\hat{\sigma}_{z}\left(\hat{G}_{i}^{R}\right)^{\dagger}\hat{\sigma}_{z}$
and the Keldysh Green’s function $\hat{G}_{i}^{K}$ is obtained from:

$\hat{G}_{i}^{K}=h(E)\left(\hat{G}_{i}^{R}-\hat{G}_{i}^{A}\right),\quad
h(E)=\tanh\left(\frac{\beta E}{2}\right)$ (11)

In the present work, terminals $S_{0,1,2}$ are respectively set to voltages
$V_{0}=0,V_{1}=V,V_{2}=-V$, and the corresponding phases obey the Josephson
relation $\varphi_{i}=\varphi_{0,i}+\frac{2eV_{i}t}{\hbar}$. Here the
reference phase $\varphi_{0}$ is set to zero. While the Green’s function of
terminal $S_{0}$ is given by Eqs. (10)-(11), those of terminals $S_{1,2}$
correspond to an out-of-equilibrium situation and must be described by a two-
time Green’s function $\check{G}_{j}(t_{1},t_{2})$. The latter is connected to
its equilibrium counterpart by the following gauge transformation:

$\check{G}_{i}(t_{1},t_{2})=e^{i\hat{\sigma}_{z}eV_{i}t_{1}}\check{G}_{i}(t_{1}-t_{2})e^{-i\hat{\sigma}_{z}eV_{i}t_{2}}$
(12)

where $\check{G}_{i}(t_{1}-t_{2})$ corresponds to the equilibrium Green’s
function given by Eqs. (1) and (10)-(11) (upon Fourier transforming back from
energy to time variable).

Due to the time-dependence in the Hamiltonian, and periodicity with a single
Josephson frequency $\omega_{0}=\frac{2eV}{\hbar}$, one can use double-time
Fourier transform and represent the Green’s functions in terms of a single
energy variable (within the interval
$[-\frac{\omega_{0}}{2},\frac{\omega_{0}}{2}]$) along with two harmonic
indices,
$\check{G}(E_{1}=E+n\omega_{0},E_{2}=E+m\omega_{0})\to\check{G}(E,n,m)$
Padurariu2 . Notice that this definition is redundant, as
$\check{G}(E,n,m)=\check{G}(E-p\omega_{0},n+p,m+p)$ for any integer $p$, such
that calculations can be performed in the energy interval
$[-\frac{\omega_{0}}{2},\frac{\omega_{0}}{2}]$. Jonckheere2009  Convolution
products in energy are then obtained as matrix products in the space of
harmonics. An alternative representation with only one harmonic index has been
used in Ref. Bezuglyi2, , in this case the matrix products appear as a
recursion relation in the harmonic indices.

Let us now write the circuit theory equations. The node is separated from each
terminal by a connector. Those express the conservation of the Keldysh matrix
current, as a discretized version of Eq. (5). As shown by Nazarov, Yuli  the
(matrix) spectral current flowing between terminal $S_{i}$ and the node takes
the form:

$\displaystyle\check{I}_{ic}=\frac{2e^{2}}{\pi\hbar}\sum_{n}\frac{T^{(i)}_{n}[\check{G}_{i},\check{G}_{c}]}{4+T^{(i)}_{n}(\\{\check{G}_{i},\check{G}_{c}\\}-2)}$
(13)

where the $T^{(i)}_{n}$ are the transmission coefficients of channel $n$
across each connector. Note that since $\check{G}^{2}=\check{1}$, both
matrices $\check{G}_{i}$ and $\check{G}_{c}$ commute with
$\\{\check{G}_{i},\check{G}_{c}\\}$. This explains the notation using the
fraction for matrix inversion in Eq. (13). To simplify, one can formally take
$N$ equivalent channels and define $T^{(i)}_{n}=T\Gamma_{i}/N$, where $T$ is a
scaling transparency factor for all contacts, and the factors $\Gamma_{i}$ are
used to describe contact asymmetries.

The physical DC current $I_{i}$ in terminal $i$ (index $c$ is omitted), as
measured in an experiment, is then given by the Nambu trace of the Keldysh
component of the matrix current $\hat{I}_{ic}^{K}$ as

$I_{i}=\int_{-\omega_{0}/2}^{\omega_{0}/2}dE\sum_{n}\text{Tr}_{N}\left[\hat{\sigma}_{z}\hat{I}_{ic}^{K}(E,n,n)\right].$
(14)

The circuit theory formulation allows to take into account the dephasing which
takes place in a normal metal diffusive conductor. In the single node model
the only source of coherence loss is due to the difference in wave-vectors of
electrons and holes at the same energy. The loss of coherence is taken into
account by connecting a fictitious terminal to the node. The current to this
terminal is given by,

$\displaystyle\check{I}_{fc}=\frac{2e^{2}}{\pi\hbar}\sum_{i}\frac{T\Gamma_{i}}{4}[\check{G}_{f},\check{G}_{c}]\
,\quad\check{G}_{f}=-i\frac{E\tau_{d}}{\hbar}\left(\begin{array}[]{cc}\hat{\sigma}_{z}&0\\\
0&\hat{\sigma}_{z}\end{array}\right).$ (17)

where $\hat{\sigma}_{z}$ is the Pauli matrix in Nambu space, and $\tau_{d}$ is
the dwell time of quasiparticles in the node, including the connectors. Notice
that no physical current flows towards this fictitious terminal.

The basic equation of circuit theory which determines $\check{G}_{c}$ stems
from the law of current conservation:

$\displaystyle\sum_{i}\check{I}_{ic}+\check{I}_{fc}=0\ .$ (18)

The same equation can be re-written as a commutation relation between the
Green’s function of the central node and a matrix $\check{M}$ that depends
both on the Green’s functions of the terminals and that of the central node:

$[\check{G}_{c},\check{M}]=0,$ (19)
$\check{M}=\sum_{i=0,1,2}T\Gamma_{i}\left(\frac{\check{G}_{i}}{1+(T\Gamma_{i}/4)(\\{\check{G}_{i},\check{G}_{c}\\}-2)}+\check{G}_{f}\right).$
(20)

The so-called ”tunnel” approximation to Eq. (13) amounts to neglecting the
term in $T_{n}^{(i)}$ in the denominator, yielding:

$\check{M}=\displaystyle{\sum_{i=\\{0,1,2\\}}}T\Gamma_{i}\left(\check{G}_{i}+\check{G}_{f}\right).$
(21)

This approximation requires that $T\ll 1$. The solution for the node Green’s
function is in this case:

$\check{G}_{c}=\frac{\check{M}}{\sqrt{\check{M}^{2}}}.$ (22)

One must stress that this approximation is much more than a simple
perturbative calculation in the tunneling transmissions. Indeed, the circuit
theory entails multiple scattering within the node. For instance, in the
absence of decoherence, $\tau_{d}=0$, and for a symmetric junction with two
terminals, the tunneling approximation turns out to be exact whatever the
$T_{n}^{(i)}=T_{n}$. VanevicBelzig2006  This is due to the fact that
$\\{G_{c},G_{i}\\}$ becomes independent of $i$ in this case and the
denominator can be factorized in Eq. (19). This property does not hold for
more than two terminals.

In practice, the Green’s functions and matrix currents are expressed in Nambu-
Keldysh-harmonics space and thus have dimension $4(2N_{m}+1)$ where $N_{m}$ is
the maximum number of harmonics of the Josephson frequency retained in the
calculation. Satisfactory convergence is obtained by using
$N_{m}=\frac{2\Delta}{eV}$ as the matrix entries decay quickly at large
harmonics value. The procedure to solve for $\check{G}_{c}$ in the linear
approximation, for each energy, is the following : $\check{M}$ is
diagonalized, its diagonal form is then normalized, so as to contain only $\pm
1$ eigenvalues, then the transfer matrix is applied back to obtain
$\check{G}_{c}$. In the general case, the nonlinear Eq. (19) is solved by
adding an iteration procedure. In practice, the current $I_{i}(\varphi_{Q},V)$
in terminal $i$ is obtained from Eq. (14), where the summation is truncated,
running from $n=-N_{m}$ to $+N_{m}$, where $N_{m}$ is the number of harmonics
retained in the calculation. Alternatively, one can integrate the traced
Keldysh component of the matrix current $\hat{I}_{ic}^{K}(E,0,0)$ over the
interval $[-E_{m},E_{m}]$, where $E_{m}$ is a large energy cutoff
(corresponding here to $E_{m}=(N_{m}+1)\omega_{0}$).

##  III Quartet and MAR currents

The three-terminal setup allows two distinct kinds of DC currents to flow
within the structure. First, dissipative currents generalize Cuevas ; Houzet ;
JonckheerePRB2013 ; AkhmerovTTJ  the MAR currents of a two-terminal junction.
As a novel feature, those currents depend on the quartet phase, due to
interfering scattering paths within the node. JonckheerePRB2013 ; Padurariu2 ;
QuartetSQUID  Second, the special condition $V_{1}=-V_{2}$ drives a DC Cooper
pair component from terminal $S_{0}$ to terminals $S_{1,2}$, as a current of
quartets (entangled pairs of Cooper pairs) $I_{Q}$. This current is
dissipationless, like the excess current stabilized at a Shapiro step in a
junction irradiated with microwaves. Cuevas

Thus one can a priori write the Keldysh currents as:

$\check{I}_{ic}=\check{I}_{Q,i}+\check{I}_{MAR,i}$ (23)

corresponding respectively to the quartets and MAR contributions. The MAR
current can be further split as

$\check{I}_{MAR,i}=\check{I}_{avMAR,i}+\check{I}_{phMAR,i}$ (24)

where $\check{I}_{avMAR,i}$ is the ”ordinary” MAR current, obtained by
averaging $\check{I}_{MAR}$ on the quartet phase, while $\check{I}_{phMAR,i}$
corresponds to the phase-dependent part.

For a symmetric TTJ, such as $\Gamma_{1}=\Gamma_{2}$, the MAR current
components flowing in terminals $S_{1}$ and $S_{2}$ are opposite, so that no
MAR current flows through terminal $S_{0}$, whose current is only due to
quartets. For this reason, we sometimes refer to $I_{Q}$ as the quartet
current from $S_{0}$ in the symmetric case. In this case, the components
$I_{Q}$ and $I_{MAR}$ are obtained by integrating over frequency [see Eq.
(14)] the spectral components defined as follows:

$\displaystyle\check{I}_{Q}$
$\displaystyle=\check{I}_{Q,0}=\check{I}_{0}=\frac{e^{2}}{h}T\Gamma_{0}[\check{G}_{0},\check{G}_{c}],$
(25) $\displaystyle\check{I}_{MAR}$
$\displaystyle=\check{I}_{MAR,1}=-\check{I}_{MAR,2}=\frac{e^{2}}{2h}T\Gamma_{1}\left[\check{G}_{1}-\check{G}_{2},\check{G}_{c}\right]$
(26)

where the products entail matrix products of the
$\check{G}_{i,c}^{\alpha,\beta}(E,n,m)$ with $\alpha,\beta$ running on the
Nambu indices $1$ (for electrons) and $2$ (for holes), and $n,m$ running as
integers.

In the general case of asymmetric contacts, one can no longer access a simple
self-contained expression for the various components of the current. Instead,
one may rely on the symmetries of the physical (energy integrated) current.
Indeed, to separate the total current in its quartet and MAR components in the
general case of asymmetric contacts, one may use the symmetries with respect
to $V$ and $\varphi_{Q}$. Time-reversal symmetry actually imposes that the
total current is odd under a change of sign of both voltage and phases. First,
as any dissipative current, the MAR components are odd in voltage thus even in
$\varphi_{Q}$. On the contrary, the quartet current is even in voltage and odd
in $\varphi_{Q}$. Therefore one has for the DC currents:

$\displaystyle I_{Q}^{i}(\varphi_{Q},-V)$ $\displaystyle=$ $\displaystyle
I_{Q}^{i}(\varphi_{Q},V)$ (27) $\displaystyle I_{Q}^{i}(-\varphi_{Q},V)$
$\displaystyle=$ $\displaystyle-I_{Q}^{i}(\varphi_{Q},V)$ (28) $\displaystyle
I_{MAR}^{i}(\varphi_{Q},-V)$ $\displaystyle=$ $\displaystyle-
I_{MAR}^{i}(\varphi_{Q},V)$ (29) $\displaystyle I_{MAR}^{i}(-\varphi_{Q},V)$
$\displaystyle=$ $\displaystyle I_{MAR}^{i}(\varphi_{Q},V)$ (30)

, This suggests to use symmetrization with respect to $V$ as a definition of
MAR and quartet currents in each contact:

$\displaystyle I_{Q}^{i}(\varphi_{Q},V)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}I^{i}(\varphi_{Q},V)+I^{i}(\varphi_{Q},-V)\big{]}$
(31) $\displaystyle I_{MAR}^{i}(\varphi_{Q},V)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}I^{i}(\varphi_{Q},V)-I^{i}(\varphi_{Q},-V)\big{]}$
(32)

A word of caution is needed about those definitions, if used in the biased
terminals. In fact, as found numerically and can be shown by symmetry
considerations, in the symmetric TTJ case $\Gamma_{1}=\Gamma_{2}$, one obtains
$I_{Q}^{1}=I_{Q}^{2}=-\frac{I_{Q}^{0}}{2}$. This corresponds to the intuition
that a quartet current splits equally into $S_{1}$ and $S_{2}$. Conversely, in
the general case $\Gamma_{1}\neq\Gamma_{2}$, one still finds
$I_{Q}^{0}+I_{Q}^{1}+I_{Q}^{2}=0$ but $I_{Q}^{1}\neq I_{Q}^{2}$, except for
low transparency, as shown by a perturbative expansion of the circuit theory
equations (see next Section). This means that high order multiple Andreev
reflections entangle different kind of processes, and the symmetrizing
procedure is merely a convenient definition separating dissipative and
nondissipative current components. Yet this definition is unambiguous in
terminal $S_{0}$: for arbitrary couplings $\Gamma_{i}$, symmetrization in $V$
allows to isolate the quartet current. Let us remark that this procedure is
convenient in calculations, but less in experiments where extrinsic causes can
perturb the symmetry between $V>0$ and $V<0$.

##  IV Perturbative expansion

Circuit theory is nonlinear, which manifests both in the normalization of the
node Green’s function and in its self-consistent definition through Eq. (19).
It is nevertheless possible to find a perturbative solution when some of the
three contacts are more weakly coupled than others. This is especially useful
in the case where the weaker couplings are with the biased contacts, that we
call a ”funnel” TTJ (Fig. 1b).

Let us consider a TTJ with couplings $\Gamma_{i}$ between the node and the
three terminals at voltages $0,V,-V$. The contact transparencies are
$T_{i}=T\Gamma_{i}$ where $T$ denotes the overall transparency of the TTJ.
When $T\ll 1$, which is assumed here, the equations for the node Green’s
function $\check{G}_{c}$ (in Nambu-Keldysh-harmonics space) and the currents
take the form:

$\displaystyle\check{M}$ $\displaystyle=$ $\displaystyle
T\sum_{i=0,1,2}\Gamma_{i}(\check{G}_{i}+\check{G}_{f})$ (33)
$\displaystyle\left[\check{G}_{c},\check{M}\right]$ $\displaystyle=$
$\displaystyle 0$ (34) $\displaystyle(\check{G}_{c})^{2}$ $\displaystyle=$
$\displaystyle 1$ (35) $\displaystyle\check{I}_{i}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{h}T\Gamma_{i}\left[\check{G}_{i},\check{G}_{c}\right]$
(36)

This leaves the couplings $\Gamma_{i}$ and the decoherence parameter
$\tau_{d}$ (see previous paper) as the only parameters of the problem. We
first take $\check{G}_{f}=0$ (no decoherence).

Since we are focusing on the case of a funnel TTJ, let us define
$\Gamma_{0}=\gamma_{0}$ ($=1$ by definition if unspecified),
$\Gamma_{1,2}=\varepsilon\gamma_{1,2}$ ($\gamma_{1,2}\sim\gamma_{0}$) where
$\varepsilon\ll 1$ is an expansion parameter. We know a priori that to obtain
the quartet current we need at least a ”trajectory” (in terms of products of
Green’s functions) passing twice through terminal $0$ and once through both
terminals $1,2$, thus second order in $\varepsilon$ is required. Let us write:

$\check{G}_{c}=\check{G}_{0}+\varepsilon\check{X}+\varepsilon^{2}\check{Y}$
(37)

and solve at order 2 for Eqs. (1-3). Defining
$\check{G}_{12}=\gamma_{1}\check{G}_{1}+\gamma_{2}\check{G}_{2}$, this yields:

$\displaystyle\\{\check{X},\check{G}_{0}\\}$ $\displaystyle=$ $\displaystyle
0$ (38) $\displaystyle\left[\check{X},\check{G}_{0}\right]$ $\displaystyle=$
$\displaystyle\left[\check{G}_{12},\check{G}_{0}\right]$ (39)
$\displaystyle\\{\check{Y},\check{G}_{0}\\}$ $\displaystyle=$
$\displaystyle-\check{X}^{2}$ (40)
$\displaystyle\left[\check{Y},\check{G}_{0}\right]$ $\displaystyle=$
$\displaystyle\left[\check{G}_{12},\check{X}\right]$ (41)

Owing to the constraint $\check{G}_{0}^{2}=1$, a solution for Eq. (38) is:

$\check{X}=\check{Z}-\check{G}_{0}\check{Z}\check{G}_{0}$ (42)

where $\check{Z}$ is an arbitrary matrix. Inserting in the other equations and
solving finally gives:

$\displaystyle\check{Z}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\check{G}_{12}$ (43) $\displaystyle\check{Y}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\big{[}\check{Z}\check{G}_{0}\check{Z}+\\{\check{Z}^{2},\check{G}_{0}\\}-3\check{G}_{0}\check{Z}\check{G}_{0}\check{Z}\check{G}_{0}\big{]}$
(44)

The current in terminal $0$ is obtained as an expansion in $\varepsilon$,
starting at first order:

$\displaystyle\check{I}_{0}$ $\displaystyle=$
$\displaystyle\varepsilon\check{I}_{0}^{(1)}+\varepsilon^{2}\check{I}_{0}^{(2)}$
(45) $\displaystyle\check{I}_{0}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{h}T\left[\check{G}_{0},\check{G}_{12}\right]$ (46)
$\displaystyle\check{I}_{0}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2h}T^{2}\left[\check{G}_{12},\check{G}_{0}\check{G}_{12}\check{G}_{0}\right]$
(47)

The first order term contains components involving terminal $0$ and one of
terminals $(1,2)$, i.e. quasiparticle currents like two-terminal MARs. One
verifies that the quartet current is contained in the second-order terms that
involve all three terminals. This yields a simple expansion for the second
order (quartet) current component:

$\displaystyle\check{I}_{0}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2h}T^{2}\big{(}\Gamma_{1}\Gamma_{2})^{2}(\check{G}_{1}\check{G}_{0}\check{G}_{2}\check{G}_{0}+\check{G}_{2}\check{G}_{0}\check{G}_{1}\check{G}_{0}$
(48) $\displaystyle-$
$\displaystyle\check{G}_{0}\check{G}_{1}\check{G}_{0}\check{G}_{2}-\check{G}_{0}\check{G}_{2}\check{G}_{0}\check{G}_{1}\big{)}$

Let us show that, in a second order expansion, the quartet current components
in $S_{1}$ and $S_{2}$, obtained by symmetrization in $V$, are equal. One
indeed shows easily that they are given by:

$\displaystyle\check{I}_{1}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2h}T^{2}\big{(}\Gamma_{1}\Gamma_{2})^{2}\big{(}\check{G}_{0}\check{G}_{2}\check{G}_{0}\check{G}_{1}-\check{G}_{1}\check{G}_{0}\check{G}_{2}\check{G}_{0}\big{)}$
$\displaystyle\check{I}_{2}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2h}T^{2}\big{(}\Gamma_{1}\Gamma_{2})^{2}(\check{G}_{0}\check{G}_{1}\check{G}_{0}\check{G}_{2}-\check{G}_{2}\check{G}_{0}\check{G}_{1}\check{G}_{0}\big{)}$
(49)

Figure 2: Spectral current (red) and NDOS (grey) for a two-terminal junction
at equilibrium. Scales are arbitrary. The spectral current peaks at the gap
and minigap edges, especially at low transparency; (top, left)
$T=0.1,\tau_{d}=0.025$; (top, right) $T=0.1,\tau_{d}=2$; (bottom, left)
$T=0.8,\tau_{d}=0.025$; (bottom, right) $T=0.8,\tau_{d}=2$.

Replacing $V$ by $-V$ amounts to exchanging $\check{G}_{1}$ and
$\check{G}_{2}$, as well as the phases $\varphi_{1}$, $\varphi_{2}$. Since the
DC current is a function of $\varphi_{1}+\varphi_{2}$ only, one sees that it
also amounts to exchanging $\check{I}_{1}^{(2)}$ and $\check{I}_{2}^{(2)}$,
once integration on energy is performed. Therefore, symmetrization in $V$
yields the same quartet component in $S_{1}$ and $S_{2}$, equal to half of
that in $S_{0}$ given by Eq. (48). Pushing the above expansion to higher order
($3$ and $4$), one can verify that this property does not hold anymore if
$\Gamma_{1}\neq\Gamma_{2}$. This shows that in the general case, if one is
interested in the current in the biased terminals, symmetrization in $V$
allows to separate a ”quartet” from a ”MAR” DC component only for a symmetric
junction. Contrarily, the symmetrizing procedure can always be used in
terminal 0. Notice that antisymmetrization in $\varphi_{Q}$ is equivalent to
symmetrization in $V$, owing to Eqs. (27)-(30).

Taking into account decoherence, if $\check{G}_{f}$ is nonzero, the zeroth
order in $\varepsilon$ for the node Green’s function is obtained by solving:

$\displaystyle\left[\check{G}_{c}^{(0)},\check{G}_{0}+\check{G}_{f}\right]$
$\displaystyle=$ $\displaystyle 0$ (50)
$\displaystyle\big{(}\check{G}_{c}^{(0)}\big{)}^{2}$ $\displaystyle=$
$\displaystyle 1$ (51)

This yields the same formal solution as in the coherent case, replacing
$\check{G}_{0}$ by $\check{G}_{c}^{(0)}$:

$\displaystyle\check{G}_{c}$ $\displaystyle=$
$\displaystyle\check{G}_{c}^{(0)}+\varepsilon\check{X}+\varepsilon^{2}\check{Y}$
(52) $\displaystyle\check{Z}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\check{G}_{12}$ (53) $\displaystyle\check{X}$
$\displaystyle=$
$\displaystyle\check{Z}-\check{G}_{c}^{(0)}\check{Z}\check{G}_{c}^{(0)}$ (54)
$\displaystyle\check{Y}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\big{[}\check{Z}\check{G}_{c}^{(0)}\check{Z}+\\{\check{Z}^{2},\check{G}_{c}^{(0)}\\}-3\check{G}_{c}^{(0)}\check{Z}\check{G}_{c}^{(0)}\check{Z}\check{G}_{c}^{(0)}\big{]}$

and the current in terminal $0$:

$\check{I}_{0}=\frac{e^{2}}{h}T\left[\check{G}_{0},\check{G}_{c}^{(0)}\right]$
(56)

which does not simplify like in Eq. (48).

This perturbative calculation can be generalized to arbitrary transparencies.

##  V Spectral quartet currents

In this Section, we discuss the spectral structure of the current components
in a symmetric TTJ ($\Gamma_{1}=\Gamma_{2}$) biased at ($0,V,-V$), calculated
in the full energy interval, in units of $\frac{e^{2}}{h}$. They are given by
${\cal I}_{ic}(E)=\Re(\check{I}_{ic}^{24}(E,0,0)-\check{I}_{ic}^{13}(E,0,0))$,
and they are compared to the spectral density (NDOS) ${\cal
N}_{c}(E)=\Re(\check{G}_{c}^{11}(E,0,0)-\check{G}_{c}^{22}(E,0,0))$, as
calculated in Ref. Padurariu2, .

As a benchmark, Fig. 2 shows the spectral current and the NDOS for a two-
terminal junction at equilibrium, at fixed phase, varying the decoherence
($\tau_{d}$) and the transparency $T$. In both spectra, one sees the minigap,
which is reduced by decoherence. The spectral current shows a sharp peak at
the minigap edges, plus a marked peak around the gap in the tunnel case ($T\ll
1$, top panels of Fig. 2). The structures are broadened at high transparency.
All take locally the form of a renormalized BCS NDOS.

Figure 3: Funnel TTJ, low transparency: (top) Spectral current (exact, thick;
perturbative, dotted); (bottom) NDOS (exact and perturbative cannot be
distinguished) ($V=0.6$, $\varphi_{Q}=\frac{2\pi}{3}$). Here
$\Gamma_{0}=1,\Gamma_{1}=0.01,\Gamma_{2}=0.03$, $T=0.02$. Curves are
vertically scaled by the squared transparency $T^{2}$. Figure 4: Funnel TTJ,
low transparency: Spectral current with moderate decoherence (exact, thick;
perturbative, dotted) ($V=0.4$, $\varphi_{Q}=\frac{2\pi}{3}$, $\tau_{d}=0.4$,
$\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.002$, $T=0.001$).

Let us now consider a ”funnel” TTJ out of equilibrium (biased at voltages
$V,-V$) with $\varepsilon\ll 1$ and $T\ll 1$ (see Section IV). Fig. 3 shows
the exact result for the spectral current together with that of the
perturbative calculation (see Sec. IV) in the coherent case ($\tau_{d}=0$),
where the agreement between the two is remarkable. For moderate decoherence
(Fig. 4), the quality of the perturbative calculation degrades. Notice that to
reproduce the exact current peaks of the numerical solution, a very small
parameter $\varepsilon$ is required, because of resonances that are difficult
to capture perturbatively. On the other hand, the NDOS is very close to that
of the strongly coupled terminal $S_{0}$, i.e. a simple BCS spectrum. Indeed,
the proximity effect in the node is here dominated by terminal $S_{0}$.

Let us now consider a fully symmetric TTJ, $\gamma_{1}=\gamma_{2}=\gamma_{0}$.
For low transparency $T\ll 1$, the density of states displays sharp BCS-like
structures centered at energy values corresponding to the six gap edges, e.g.
given by $E=\pm\Delta,\pm(\Delta+V),\pm(\Delta-V)$. Padurariu2  Due to higher
order MAR processes, weaker structures also occur at $E=\pm(\Delta\pm nV)$
with $|n|>1$. Fig. 5 shows the spectral current in terminal $S_{0}$ (which
reduces to the quartet contribution in this case), together with the NDOS,
represented with a different scale. One sees that the spectral current peaks
at the same energies as the NDOS. Remarkably, its sign changes from one
structure to another. For large transparency (bottom panels of Fig. 5), the
peaks broaden but the spectral currents still oscillate from one structure to
the next. For large decoherence and low transparency, one distinguishes three
structures at $E\sim 0,\pm 2V$, of width $\sim\frac{1}{\tau_{d}}$, the latter
being the mirror of the $E\sim 0$ region through Andreev reflection in
terminals $1,2$.

Figure 5: Spectral current (red) and NDOS (grey) for a fully symmetric TTJ,
$V=0.65$, $\varphi_{Q}=\frac{2\pi}{3}$, $\Gamma_{0}=\Gamma_{1}=\Gamma_{2}=1$
(left, $\tau_{d}=0.025$; right, $\tau_{d}=2$). Top : $T=0.1$. Bottom :
$T=0.8$. Current for $\tau_{d}=2,T=0.1$ is scaled up by a factor $5$ to help
visibility. The spectral current exhibits structures together with the NDOS,
with marked sign changes.

##  VI Phase and voltage dependence of the quartet currents.

The quartet phase $\varphi_{Q}$ is stationary for any bias voltage $V$.
Therefore the voltage appears as a novel control parameter for a DC Josephson
current, in contrast with two-terminal junctions. As shown below and in next
Section, the dependence of the quartet and the phase-MAR currents with $V$ is
highly nontrivial and reflects the resonant dynamics of a driven quantum
system. Here the quartet current $I_{Q}$ is calculated in terminal $S_{0}$ for
a TTJ such that $\Gamma_{1}=\Gamma_{2}$, but with the ratio
$\Gamma_{1}/\Gamma_{0}=\varepsilon$ taking any value.

###  VI.1 Fully symmetric TTJ, $\varepsilon\simeq 1$

Let us first consider the situation of a fully symmetric TTJ (see Fig. 1a)
with $\varepsilon\simeq 1$. We fix the quartet phase and plot the quartet DC
current as a function of $V$, as an integral over energy of the spectral
current. In all figures the current is in units of $\frac{e^{2}\Delta}{h}$.
The behavior of a symmetric TTJ is rather complex, with several sign changes
of the current (Fig. 6). Fig. 6b shows the effect of the transparency at zero
decoherence, Fig. 6c corresponds to strong decoherence. In Fig. 6d
transparency is fixed and decoherence is varied. The calculation is limited by
numerics for small $V$ and other sign changes may occur in the very low $V$
regime (where the current experiences strong oscillations with $V$). A similar
behavior was obtained in a separate quantum dot model calculation. Regis1 ;
Regis2 ; Regis3

Figure 6: (a) Symmetric TTJ quartet current as a function of $V$, low
transparency ($\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=1.5$, $T=0.01$,
$\varphi_{Q}=\frac{2\pi}{3}$, $\tau_{d}=1.25$); (inset) corresponding quartet
current-phase characteristics across a sign change. $V=0.495$ (blue); $0.5095$
(red, upscaled by factor $5$); $0.52$ (black); (b) TTJ quartet current as a
function of $V$ for different transparencies, coherent case,
$\varphi_{Q}=\frac{2\pi}{3}$. Curves are scaled (divided) by the transparency
$T$ ($T=0.01$, black; $T=0.3$, green, $T=0.7$, red). $\varphi_{Q}=2\pi/3$,
$\Gamma_{0}=1,\Gamma_{1}=1.5,\Gamma_{2}=0.6$, $\tau_{d}=0$; (c) TTJ quartet
current as a function of $V$ for different transparencies, incoherent case;
same scaling and parameters as (b), $\tau_{d}=1.25$; (d) TTJ quartet current
as a function of $V$ for different decoherence rates, at low transparency
($\tau_{d}=0$, black; $\tau_{d}=0.3$, green, $\tau_{d}=3$, red).
$\varphi_{Q}=2\pi/3$, $\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.8$, $T=0.01$ (same
scaling as (b)).

Figure 6a shows the quartet current-phase characteristics close to a sign
change, around $V=0.5\Delta$. The anharmonicity is strong very close to
$V=0.5\Delta$ at the ”$\pi-0$” transition. A general trend is as follows:
”$\pi$”-junction at low voltage, ”$0$”-junction at intermediate voltage and
”$\pi$”-junction at higher voltage.

The quartet current vanishes at high voltages (Figs. 6b, 6c, and 6d)
regardless of transparency and decoherence ($V>2\Delta$ signals the
quasiparticle transfer onset between the side and central leads).

Changing the transparency does not affect “much” the overall shape of the
quartet current - voltage characteristics: the voltage at which the current
changes sign then depends on the transparency $T$, especially in the presence
of decoherence (Figs. 6b and 6c), and modifies its amplitude (we stress that
the curves have been rescaled by the transparency for convenience). The
introduction of moderate decoherence $\tau_{d}=1$ reduces somehow the
amplitude of the quartet current. Fig. 6d shows the effect of decoherence in
the tunnel regime, where we see that strong decoherence cancels the sign
change at large voltages.

Let us emphasize that the occurrence of a new control parameter $V$ for the
quartet current offers the following possibility : by changing the overall
transparency of the TTJ, or that of one or two contacts only, one can change
the sign of the current, provided the voltage $V$ is suitably chosen. This
extends the range of phenomena causing sign changes in a Josephson current.
This might be useful in devices where commuting a junction from ”$\pi$” to
”$0$” is required.

###  VI.2 ”Funnel” TTJ

Let us now consider a ”funnel” TTJ with $\varepsilon\ll 1$ (Fig. 1b), and
compare the exact current with that found by second order expansion in
$\varepsilon$ (see Sec. IV) The result is displayed in Fig. 7a. In the
coherent case, the agreement with the perturbative result is qualitatively
good if $\varepsilon$ is small enough. Notice that the total current is the
result of an integral over a strongly oscillating function, and the result is
therefore quite sensitive to approximations. Remarkably, the sign of the
quartet current is always negative, apart for high values of $V$. This
corresponds to a ”$\pi$” quartet junction. Fig. 7b shows the quartet current-
phase characteristics for a given $V$, which is close to a harmonic relation
$I_{Q}=I_{Qc}\sin\varphi_{Q}$, with $I_{Qc}<0$. The quartet critical current
peaks around $V=\Delta$.

Conversely, with moderate decoherence, the agreement with the perturbative
treatment in $\varepsilon$ is less good, especially for intermediate voltages
(Fig. 7c): in this region, the quartet current changes sign, from a
”$\pi$”-junction (low voltage) to a ”$0$”-junction (high voltage). The quartet
characteristics $I_{Q}(\varphi_{Q})$ at fixed $V$ becomes strongly anharmonic
in the region of the sign change (Fig. 7d).

Figure 7: (a) Funnel TTJ quartet current as a function of voltage, coherent
case (exact, thick; perturbative, dotted), $\varphi_{Q}=2\pi/3$, $\tau_{d}=0$;
(b) Quartet current-phase characteristics, $V=0.5$ (blue), $1$ (red), $1.5$
(black); $\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.03$, $T=1$,
$\varphi_{Q}=2\pi/3$, $\tau_{d}=0$; (c) Funnel TTJ quartet current as a
function of voltage (exact, thick; perturbative, thin), $\varphi_{Q}=2\pi/3$,
$\tau_{d}=0.4$; (d) Quartet current-phase characteristics, $V=0.8$ (blue),
$0.855$ (red, vertical upscaling by a factor $500$), $0.9$ (black);
$\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.005$, $T=0.1$, $\varphi_{Q}=2\pi/3$,
$\tau_{d}=0.4$.

Let us end this Section by proposing an indirect experimental detection of the
sign changes, when a phase-sensitive experiment is not possible. Fig. 8 shows
the maximum (”critical”) quartet current, taken on all possible quartet phase
values. It displays sharp minima when the phase-sensitive current changes
sign. Similar trends were observed in ferromagnetic Josephson junctions as a
witness of $0/\pi$ transitions, as a function of the junction length. SFS
Here, such minima observed as a function of $V$ are hardly explainable by
conventional models based on extrinsic synchronization by the outer circuit,
as often advocated in experiments against the mesoscopic quartet mechanism.

Figure 8: Nonmonotonous variation with $V$ of the TTJ maximum quartet current
(red), taken on the quartet phase ($\Gamma_{0}=\Gamma_{1}=\Gamma_{2}=1$,
$T=1$, $\tau_{d}=0$).

##  VII MAR currents and the balance with quartet currents

###  VII.1 Phase symmetry of MAR currents

We first consider a junction which is symmetric in terminals $1,2$, i.e.
$\Gamma_{1}=\Gamma_{2}$. The full MAR current is defined as
$I_{MAR}=\frac{e^{2}}{2h}T\Gamma_{1}(I^{1}-I^{2})$, and there is no MAR
current in terminal $S_{0}$. The phase-dependent part of the MAR current
$I_{phMAR}^{i}$ is defined by subtracting from the full MAR current its
average on the quartet phase $\varphi_{Q}$:

$I_{phMAR}^{i}=I_{MAR}^{i}-\langle I_{MAR}^{i}\rangle_{\varphi_{Q}}$ (57)

The phase-MAR current is represented as a function of the phase $\varphi_{Q}$
on Fig. 9. It is an even function of $\varphi_{Q}$, contrarily to the quartet
current which is odd. As a function of voltage, the phase-dependent MAR
current also presents sign changes, at values of $V$ different from those
making the quartet component change sign.

Figure 9: (top) TTJ quartet current (black) and phase MAR current (red) as a
function of $V$ ($\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.8$,
$\varphi_{Q}=2\pi/3$, $T=0.01$, $\tau_{d}=0$); (bottom) Phase MAR current-
phase characteristics across a sign change. $V=0.25$ (full); $0.76$ (dotted).

###  VII.2 A TTJ is a $\varphi$-junction

Owing to the different symmetries of quartet and phase-MAR currents with
respect to the quartet phase, their superposition yields a phase shift, that
makes a TTJ a tunable ”$\varphi$-junction”: the phase shift indeed depends on
the parameters of the TTJ : bias voltage $V$, transparency, decoherence, and
also on the two relevant asymmetries: that comparing contact $0$ to contacts
$(1,2)$ (see the ”funnel” case), and that mutually comparing contacts $1,2$.

This is illustrated by two examples: first, within $(S_{1},S_{2})$ symmetry,
the current in terminal $S_{1}$ superimposes quartet and MAR components thus
displays a phase shift, which in Figure 10(top) varies with the voltage. As
$V$ increases, it passes from a quartet-dominated current, nearly zero at zero
phase, to a MAR-dominated current, with phase-MAR component maximal near zero
phase. Second, with an $(S_{1},S_{2})$ asymmetry, the current in terminal
$S_{0}$ contains MARs in addition to quartets and displays a phase shift, here
at fixed voltage and varying the asymmetry (Figure 10(bottom)). It goes from a
quartet regime to a MAR regime as asymmetry is increased. Provided phase-MARs
are not negligible, such a phase shift could be observed in a phase-sensitive
experiment as proposed in Ref. QuartetSQUID, .

Figure 10: (top) Phase shift in the current in $S_{1}$ (normalized by a factor
$\gamma_{1}$), with $\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=0.2$, $\tau_{d}=0.5$,
$V=0.2$ (black, quartet-dominated current), $0.45$ (blue), $0.55$ (green),
$0.7$ (red, MAR-dominated current)). Dotted lines feature the phase-
independent MAR current; (bottom) Phase shift in the current in $S_{0}$, with
$V=0.4$, $\tau_{d}=1$, $\Gamma_{0}=1$, and $\Gamma_{1}=\Gamma_{2}=1$ (black,
quartet-dominated current), $\Gamma_{1}=1.2,\Gamma_{2}=0.8$ (blue),
$\Gamma_{1}=1.4,\Gamma_{2}=0.7$ (green), $\Gamma_{1}=2,\Gamma_{2}=0.5$ (red,
MAR-dominated current)).

###  VII.3 Conditions to observe quartets against MARs

In a non phase-sensitive transport experiment, and in the case where MARs are
observed in the full $(V_{1},V_{2})$ conductance map, care must be taken in
interpreting the $V-I$ anomaly observed on the ”quartet line” $V_{2}=-V_{1}$.
For a rather symmetric TTJ one finds that $I_{Q}$ and $I_{phMAR}$ have
comparable magnitudes, even at low voltage (Fig. 9).This is due to the fact
that in the present model simulating a cavity with a dense level spectrum, a
symmetric TTJ behaves as a resonant junction and therefore MARs can be strong
even at low voltage. On the contrary, when
$\Gamma_{1},\Gamma_{2}\ll\Gamma_{0}$, the phase-MAR current becomes very small
at low voltages and is dominated by the quartet current (Fig. 11).

Figure 11: TTJ quartet current (red) and phase MAR current (blue) as a
function of $V$ ($\Gamma_{0}=1,\Gamma_{1}=\Gamma_{2}=\varepsilon\Gamma_{0}$,
$\varphi_{Q}=2\pi/3$, $\tau_{d}=0$); (top) $\varepsilon=0.1$; (bottom)
$\varepsilon=0.02$

As shown above, a funnel asymmetry is sufficient to make phase-MAR very small.
This is because those processes involve coherently all three contacts and
taking $\varepsilon\ll 1$ suppresses multiple Andreev reflections between
$S_{0}$ and both ($S_{1},S_{2}$). It also trivially hampers two-terminal MARs
occurring either between $S_{0}$ and $S_{1}$, or between $S_{0}$ and $S_{2}$,
which participate to the phase-independent MAR current. Nevertheless, if
couplings between the node and the leads $S_{1},S_{2}$ are comparable albeit
small, MARs between those contacts (with voltage difference $2V$) become
resonant and can overcome quartets. Fig. 12 shows that if one instead takes
$\Gamma_{0}\gg\Gamma_{1}\gg\Gamma_{2}$ (or
$\Gamma_{0}\gg\Gamma_{2}\gg\Gamma_{1}$) (full asymmetry, see Fig. 1c), all
two-terminal MARs are suppressed at low $V$ and the total MAR current becomes
negligible compared to the quartet current below $V=0.2\Delta$ in the coherent
case.

Figure 12: Quartet and full MAR currents in a fully asymmetric TTJ,
$\Gamma_{0}=1,\Gamma_{1}=0.2,\Gamma_{1}=0.05$. MAR currents in terminals
$S_{0}$ and $S_{1}$ strongly decay at low voltages, contrarily to quartet
currents. Current sign convention is inverted in terminals $1,2$ to help
comparison. Notice that here $I_{q0}\simeq 2I_{q1}\simeq 2I_{q2}$.

##  VIII Discussion and conclusion

In this work we have studied in detail the voltage and phase dependence of the
two DC current components involving quartets in a three-terminal Josephson
junction biased with opposite voltages : quartet current and phase-MAR
current. The salient feature is the sign changes occurring as the voltage is
varied. The quartet current switches from a $\pi$ to a $0$ junction. These
sign changes are absent at low voltage when the junction is strongly
asymmetric in terms of the coupling to the unbiased ”source” terminal compared
to the two others (”funnel junction”). This, and the agreement between a
perturbative and the exact calculation, points towards non-adiabatic effects
as being responsible for such sign changes. This is similar to the behavior
reported in quantum dot junctions. Harvard ; Regis1 ; Regis2 ; Regis3  In the
funnel junction, a wide pseudogap in the normal density of states protects
from nonadiabatic transitions between negative and positive energy Andreev
bound states, and the finite voltage behavior emerges adiabatically from the
zero-voltage one which manifests a $\pi$-junction.

The effect of decoherence has been also investigated, as well as the overall
transparency of the contacts. This leads to shifts of the sign changes, that
sensitively depend on the multiple Andreev reflections occurring at the
contacts. At low voltages, phase-MAR currents become negligible compared to
quartet currents in a funnel asymmetric TTJ. All MAR currents become
negligible if the low-voltage TTJ is doubly asymmetric. When phase-MARs are
not negligible, which is the generic situation, the TTJ becomes a controllable
$\varphi$-junction.

Probing the sign changes of the quartet current, more generally the
$\varphi$-TTJ, requires a full quartet SQUID geometry. QuartetSQUID  A less
demanding probe is the strongly non-monotonic variation of the quartet
”critical” current $I_{qc}$, that can be observed in a direct transport
experiment, as measured by the width of the Josephson-like anomaly found on
the quartet line. This non-monotonic dependence and the strong sensitivity of
the positions of the $I_{Qc}$ minima with voltage are a hallmark of quartet
processes. They are directly related to the mesoscopic nature of the quartet
process and can thus be distinguished from other synchronization mechanisms
due to an external circuit impedance, stronger than that of the junction
itself : such an impedance has no reason to depend sensitively on the details
(transparency, decoherence, asymmetries) of the TTJ.

Potential extensions of this work include the computation of the noise
characteristics of the TTJ within the circuit theory framework. In previous
works where the diffusive node between the superconducting leads is replaced
by a system of quantum dots Regis1 ; Jacquet  showed that at low temperature,
the Fano factors relevant to the TTJ vanish for decreasing voltage, as
expected for a Josephson-like signal (the same should be true in the context
of circuit theory).

###### Acknowledgements.

This work received support from the French government under the France 2030
investment plan, as part of the Initiative d’Excellence d’Aix-Marseille
Université - A*MIDEX (authors J. R., T. J., L.R., and T. M.). Centre de Calcul
Intensif d’Aix-Marseille is acknowledged for granting access to its high
performance computing resources.

## References

  * (1) T. M. Klapwijk, G. E. Blonder, and M. Tinkham, Physica B+C 109-110, 1657 (1982); M. Octavio, M. Tinkham, G. E. Blonder, and T. M. Klapwijk, Phys. Rev. B 27, 6739 (1983).
  * (2) E. Scheer, P. Joyez, D. Esteve, C. Urbina, and M. H. Devoret, Phys. Rev. Lett. 78, 3535 (1997).
  * (3) D. J. Thouless, Phys. Rev. Lett. 39, 1167 (1977).
  * (4) E. V. Bezuglyi, E. N. Bratus, V. S. Shumeiko, G. Wendin, Phys. Rev. Lett. 83, 2050 (1999).
  * (5) E. V. Bezuglyi, E. N. Bratus, V. S. Shumeiko, Phys. Rev. B 83(18), 184517 (2011).
  * (6) D. Feinberg, T. Jonckheere, J. Rech, T. Martin, Douçot and R. Mélin, Eur. J. Phys. B 88, 99 (2015).
  * (7) J. Rech, T. Jonckheere, T. Martin, B. Douçot, R. Mélin and D. Feinberg, Phys. Rev. B 90, 075419 (2014).
  * (8) B. van Heck, S. Mi, A. R. Akhmerov, Phys. Rev. B 90, 155450 (2014).
  * (9) R. Mélin, D. Feinberg, and B. Douçot, Eur. Phys. J. B 89, 67 (2016).
  * (10) R. -P Riwar, M. Houzet, J. S. Meyer, and Yu. V. Nazarov, Nat. Comm. 7, 11167 (2016).
  * (11) J. S. Meyer and M. Houzet, Phys. Rev. Lett. 119, 136807 (2017).
  * (12) J.-D. Pillet, V. Benzoni, G. Griesmar, J.-L Smirr, and Ç. Ö. Girit, Nanolett. 7138 (2019).
  * (13) J. C. Cuevas and H. Pothier, Phys. Rev. B75, 174513 (2007).
  * (14) M. Houzet and P. Samuelsson, Phys. Rev. B82, 060517 (2010).
  * (15) A. Freyn, B. Douçot, D. Feinberg, and R. Mélin, Phys. Rev. Lett. 106, 257005 (2011).
  * (16) T. Jonckheere, J. Rech, T. Martin, B. Douçot, D. Feinberg, and R. Mélin, Phys. Rev. B 87, 214501 (2013).
  * (17) M. P. Nowak, M. Wimmer, A. R. Akhmerov, Phys. Rev. B 99, 075516 (2019).
  * (18) A. Keliri and B. Douçot, Phys. Rev. B 107, 094505 (2023)
  * (19) E. Strambini, S. D’Ambrosio, F. Vischi, F. S. Bergeret, Yu. V. Nazarov, F. Giazotto, Nature Nanotechnology 11, 1055 (2016); F. Vischi, M. Carrega, E. Strambini, S. D’Ambrosio, F. S. Bergeret, Yu. V. Nazarov and F. Giazotto, Phys. Rev. B95, 054504 (2017).
  * (20) A. H. Pfeffer, J. E. Duvauchelle, H. Courtois, R. Mélin, D. Feinberg and F. Lefloch, Phys. Rev. B 90, 075401 (2014).
  * (21) Y. Cohen, Y. Ronen, J.-H. Kang, M. Heiblum, D. Feinberg, R. Mélin and H. Shtrikman, Proc. Natl. Acad. Sci. U. S. A. 115, 6991 (2018).
  * (22) S.A. Khan, L. Stampfer, T. Mutas, J.-H. Kang, P. Krogstrup and T.S. Jespersen, Multiterminal Quantized Conductance in InSb Nanocrosses, Advanced Materials 33, 2100078 (2021).
  * (23) N. Pankratova, H. Lee, R. Kuzmin, K.Wickramasinghe,W. Mayer, J. Yuan, M. G. Vavilov, J. Shabani and V. E. Manucharyan, Phys. Rev. X 10, 031051 (2020),
  * (24) G. V. Graziano, J. S. Lee, M. Pendharkar, C. J. Palmström and V. S. Pribiag, Transport studies in a gate-tunable three-terminal Josephson junction. Phys. Rev. B 101, 054510 (2020)
  * (25) G. V. Graziano, M. Gupta, M. Pendharkar, J. T. Dong, C. P. Dempsey, C. Palmstrøm and V. S. Pribiag, Selective control of conductance modes in multi-terminal Josephson junctions, Nat. Comm. 13, 5933 (2022).
  * (26) K.-F. Huang, Y. Ronen, R. Mélin, D. Feinberg, K. Watanabe, T. Taniguchi and P. Kim, Nat. Comm. 13, 3032 (2022).
  * (27) A. W. Draelos, M.-T. Wei, A. Seredinski, H. Li, Y. Mehta, K. Watanabe, T. Taniguchi, I. V. Borzenets, F. Amet, and G. Finkelstein, Nano Lett. 19, 1039 (2019).
  * (28) E.G. Arnault, T. Larson, A. Seredinski, L. Zhao, H. Li, K. Watanabe, T. Taniguchi, I. Borzenets, F. Amet and G. Finkelstein, The multiterminal inverse AC Josephson effect, Nano Lett. 21, 9668 (2021).
  * (29) 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, Nano Lett. 22, 7073 (2022).
  * (30) F. Zhang, A.S. Rashid, M.T. Ahari, W. Zhang, K.M. Ananthanarayanan, R. Xiao, G.J. de Coster, M.J. Gilbert, N. Samarth and M. Kayyalha, Andreev processes in mesoscopic multi-terminal graphene Josephson junctions, arXiv: 2210.04408v2 (2022).
  * (31) R. Mélin and D. Feinberg, Phys. Rev. B 107, L161405 (2023).
  * (32) C. Padurariu, T. Jonckheere, J. Rech, T. Martin and D. Feinberg, Phys. Rev. B 95, 205437 (2017).
  * (33) A.K. Jain, K.K. Likharev, J.E. Lukens, and J.E. Sauvageau, Physics Reports 109, 309 (1984).
  * (34) A. Melo, V. Fatemi and A. R. Akhmerov, SciPost Phys. 12, 017 (2022).
  * (35) R. Mélin, M. Sotto, D. Feinberg, J.-G. Caputo, and B. Douçot, Phys. Rev. B 93, 115436 (2016).
  * (36) R. Mélin, Phys. Rev. B 102, 245435 (2020); R. Mélin and B. Douçot, Phys. Rev. B 102, 245436 (2020).
  * (37) R. Mélin, Phys. Rev. B 105, 155418 (2022).
  * (38) Yu. V. Nazarov, Superlattices Microstruct. 25, 1221-1231 (1999).
  * (39) Yu. V. Nazarov and Ya. M. Blanter, Quantum transport: Introduction to Nanoscience, (Cambridge University Press, 2009).
  * (40) C. Padurariu, T. Jonckheere, J. Rech, R. Mélin, D. Feinberg, T. Martin and Yu. V. Nazarov, Phys. Rev. B 92, 205409 (2015).
  * (41) T. Jonckheere, A. Zazunov, K.V. Bayandin, V. Shumeiko, and T. Martin, Phys. Rev. B 80, 184510 (2009).
  * (42) E. V. Bezuglyi, E. N. Bratus, V. S. Shumeiko, G. Wendin, and H. Takayanagi, Phys. Rev. B 62, 14439 (2000).
  * (43) A. I. Larkin and Yu. V. Ovchinninkov, Sov. Phys. JETP 41, 960 (1975); A. I. Larkin and Yu. V. Ovchinninkov, Sov. Phys. JETP 46, 155 (1977).
  * (44) M. Vanević and W. Belzig, Europhys. Lett. 75, 604 (2006).
  * (45) A. L Buzdin,L. N. Bulaevskil,andS. V. Panyukov, JETP Lett. 35, 178 (1982).
  * (46) R. Jacquet, A. Popoff, K. -I. Imura, J. Rech, T. Jonckheere, L. Raymond, A. Zazunov, and T. Martin, Phys. Rev. B 102, 064510 (2020)