# Josephson effect in superconductor-normal dot-superconductor junctions
driven out of equilibrium by quasiparticle injection

Sarath Sankar, Julia S. Meyer, and Manuel Houzet  Univ. Grenoble Alpes, CEA,
Grenoble INP, IRIG, Pheliqs, F-38000 Grenoble, France

###### Abstract

We study theoretically the large variations of the supercurrent through a
normal dot that are induced by a small quasiparticle injection current from
normal leads connected to the dot. We find that the supercurrent decomposes
into a subgap contribution, which depends on the voltages applied to the
normal leads, as well as a contribution with opposite sign from energies
outside the gap, which is insensitive to the voltages. As the voltages
gradually suppress the subgap contribution, a critical voltage exists above
which the contribution from energies outside the gap dominates, leading to a
sign reversal of the current-phase relation, namely a transition to a so-
called $\pi$-junction behavior. We determine the critical voltage and analyze
the robustness of the effect with respect to temperature and inelastic
relaxation in the dot.

##  I Introduction

The fascinating properties of superconducting devices are in large part
associated with the presence of an excitation gap for quasiparticles such
that, at low temperatures, the system can be described in terms of a
macroscopic variable, the superconducting phase, while being protected from
dissipation. It turns out, however, that quasiparticles still play an
important role in many situations. Nonequilibrium quasiparticles may be
present due to the breaking of Cooper pairs, e.g., by the absorption of stray
photons or cosmic rays. They have been shown to be very difficult to get rid
off and to be detrimental for the coherence of superconducting qubits (see,
e.g., Ref. Glazman and Catelani (2021) for a recent review). On the other
hand, trapped quasiparticles in the Andreev bound states that form in
Josephson junctions can lead to interesting novel phenomena, such as the
realization of an Andreev spin qubit Padurariu and Nazarov (2010); Hays et al.
(2021).

Recent experiments on gate control of the supercurrent in metallic Josephson
junctions De Simoni et al. (2018); Paolucci et al. (2018); De Simoni et al.
(2019); Paolucci et al. (2019); Rocci et al. (2020) have revived the interest
in a better understanding of the role of quasiparticle injection from normal
parts of the circuit Golokolenov et al. (2021); Ritter et al. (2021a, b);
Basset et al. (2021). Deliberate quasiparticle injection via voltage-biased
normal leads has indeed been studied earlier and shown to have important
effects on the supercurrent Chang and Bagwell (1997); Morpurgo et al. (1998);
Savin et al. (2004). The possible reversal of the sign of the supercurrent
upon increasing the voltage of a normal lead directly coupled to the junction
was first shown, though not emphasized, theoretically in a long ballistic
junction Van Wees et al. (1991). The resulting realization of a so-called
nonequilibrium $\pi$-junction was pointed out in Ref. Volkov (1995), where a
simpler setup consisting of a normal dot connected to two superconductors was
investigated. Experimentally, a nonequilibrium $\pi$-junction was first
realized in a long diffusive junction Baselmans et al. (1999). Also subsequent
work, both theoretical Wilhelm et al. (1998); Yip (1998); Heikkilä et al.
(2002) and experimental Huang et al. (2002); Baselmans et al. (2002),
concentrated on extended junctions. (A short ballistic junction was addressed
in Ref. Heikkilä et al. (2002). However, in that case, the effect is absent.)
More complicated geometries, sometimes called Andreev interferometers, have
been studied as well Samuelsson et al. (2000); Dolgirev et al. (2019a), but
all in the long-junction limit. We note in passing that equilibrium
$\pi$-junctions may be realized in superconductor-ferromagnet-superconductor
junctions Bulaevskii et al. (1977); Buzdin et al. (1982); Ryazanov et al.
(2001); Kontos et al. (2002).

Here we revisit the simple superconductor-normal dot-superconductor setup and
analyze the supercurrent in detail. Our main findings are that the
nonequilibrium $\pi$-transition exists irrespective of the coupling strengths
between the dot and the superconductors, and that it is robust with respect to
temperature and inelastic relaxation due to electron-electron interactions.
Furthermore, we show that the same phenomenon also occurs in multiterminal
Josephson junctions Cohen et al. (2018); Draelos et al. (2019); Pankratova et
al. (2020); Graziano et al. (2020) that are currently under intensive
investigation, as they may be used for Majorana braiding Alicea et al. (2011);
Aasen et al. (2016) and have interesting topological properties Van Heck et
al. (2014); Yokoyama and Nazarov (2015); Riwar et al. (2016)

Using the quasiclassical Usadel equations in the Keldysh formulation Belzig et
al. (1999), we study a chaotic normal dot coupled to various superconducting
and normal leads. We obtain explicit analytical expressions for the currents
flowing from the dot to the different leads. Whereas the links to the normal
leads only carry dissipative currents, the links to the superconducting leads
may carry both dissipative currents and dissipationless supercurrents. We find
that the dissipative and supercurrents can be distinguished by their behavior
under a gobal change of sign of all the voltages: while the dissipative
currents are odd, the supercurrents are even under such a sign reversal.

The distinct behavior under a global change of sign of all voltages allows us
to separately study the supercurrents, which can be expressed as an energy-
integral over the product of a spectral function and a distribution function.
While the spectral function only depends on the spectral properties of the
leads and the couplings between the dot and the leads, the distribution
function also depends on the distribution functions of the leads.
Interestingly the contributions to the supercurrent from subgap energies and
from energies outside the gap show a quite different behavior. The subgap
contributions depend only on the voltage-dependent distribution functions of
the normal leads, which is easy to understand as the gap in the
superconducting leads prevents thermalization at these energies. By contrast,
the contributions from energies outside the gap depend only on the equilibrium
distribution functions of the superconducting leads. Thus the quasiparticle
injection from the normal leads only affects the subgap contributions and, in
particular, suppresses them upon increasing the bias voltages.

Furthermore, we show that, at fixed superconducting phases, the signs of the
subgap contributions and the contributions from energies outside the gap are
opposite, such that the two contributions are in competition. This effect can
be traced back to the negative sign of the spectral current function at
energies above the gap, as emphasized in the two-terminal case in Ref.
Levchenko et al. (2006). In equilibrium, we find that the subgap contributions
always dominate and, therefore, fix the signs of the supercurrents. However,
as the applied voltages to the normal leads are increased, the gradual
suppression of the subgap contributions leads to a supercurrent sign reversal
at a critical value of the voltages, where the contributions from energies
outside the gap become dominant.

In a two-terminal junction, a sign reversal of the supercurrent corresponds to
a transition from a conventional junction to a so-called nonequilibrium
$\pi$-junction. We study this effect in detail in the case of weak and strong
coupling of the dot to the superconducting leads, both at zero temperature and
close to the superconducting critical temperature, $T_{c}$. At zero
temperature, we find that the transition always occurs at voltages
$V^{*}<\Delta_{0}/e$, where $\Delta_{0}$ is the zero-temperature
superconducting gap. Upon further increasing the voltage above $\Delta_{0}/e$,
the critical current saturates to a value that is smaller than or comparable
to the equilibrium critical current. Close to $T_{c}$, the behavior is
qualitatively different. At strong coupling, the transition is pushed up to
voltages $V^{*}>\Delta_{0}/e$. By contrast, at weak coupling, the finite
temperature reduces the value of the critical voltage $V^{*}$. Interestingly,
in that case, the critical current at large voltages in the $\pi$-junction
regime is parametrically larger than the equilibrium supercurrent.

As in any nonequilibrium phenomenon, relaxation plays an important role. We
study the robustness of the effects discussed above with respect to internal
inelastic relaxation inside the dot due to electron-electron interactions.
Modeling this inelastic relaxation by a fictitious fermionic bath Büttiker
(1986), we find that the transition survives even when the inelastic
relaxation rate $\Gamma_{b}$ is fairly large, though it is pushed to large
voltages. Thus the effect is robust.

The paper is organized as follows: In Sec. II, we define the model and provide
the general expressions for the currents in the quasiclassical-Keldysh
formalism. In Sec. III, we analyze the supercurrents in the absence of
inelastic relaxation and obtain the critical voltage at which the
$\pi$-transition happens. In Sec. IV, the effects of inelastic relaxation
within the dot due to electron-electron interactions are considered. Finally
we conclude in Sec. V. Some details and generalizations of the setup can be
found in the Appendices.

##  II The model

We consider a normal dot that is coupled to normal (N) and superconducting (S)
leads through tunnel barriers. In the quasiclassical-Keldysh formalism, the
electric current through the junction connecting the dot to lead $p$ can be
expressed in terms of the quasiclassical Green functions of the dot
($\check{g}$) and the lead ($\check{g}_{p}$), subject to the normalization
condition $\check{g}_{(p)}^{2}=1$, as follows:

$I_{p}=\frac{G_{p}}{16e}\int_{-\infty}^{\infty}dE\,\tilde{I}_{p}(E)$ (1)

with

$\tilde{I}_{p}=\text{Tr}\left[\hat{\tau}_{3}\left(\hat{g}_{p}^{K}\hat{g}^{A}-\hat{g}^{R}\hat{g}_{p}^{K}-\hat{g}^{K}\hat{g}_{p}^{A}+\hat{g}_{p}^{R}\hat{g}^{K}\right)\right],$
(2)

where $G_{p}$ is the tunnel conductance in the normal state, and
$\hat{g}_{(p)}^{R/A/K}$ denotes the retarded, advanced, and Keldysh components
of the quasiclassical Green function, respectively:

$\check{g}_{(p)}=\left(\begin{array}[]{cc}\hat{g}_{(p)}^{R}&\hat{g}_{(p)}^{K}\\\
0&\hat{g}_{(p)}^{A}\end{array}\right).$ (3)

The components $\hat{g}_{(p)}^{R/A/K}$ are $2\times 2$ matrices in Nambu
space, and $\hat{\tau}_{1,2,3}$ are Pauli matrices acting in that space.

If we neglect inelastic relaxation in the dot, the dot Green function is
determined by the equation

$[\check{h},\check{g}]=0\quad{\rm
with}\quad\check{h}=E\check{\tau}_{3}+i\sum_{p}\Gamma_{p}\check{g}_{p}.$ (4)

Here $\check{\tau}_{3}=\hat{\sigma}_{0}\otimes\hat{\tau}_{3}$ with
$\hat{\sigma}_{0}$ being the identity matrix in Keldysh space, $E$ is the
energy measured from the Fermi level in the S leads,
$\Gamma_{p}=G_{p}\delta/(2\pi G_{Q})$, where $\delta$ is the mean level
spacing in the dot and $G_{Q}=e^{2}/\pi$ the conductance quantum (in units
where $\hbar=1$), are the partial level widths due to the tunnel coupling of
the dot to the leads.

The spectral Green functions $\hat{g}_{(p)}^{R/A}$ can be generally expressed
as

$\hat{g}^{R}_{(p)}=\sin\theta_{(p)}\left[\sin\chi_{(p)}\hat{\tau}_{1}+\cos\chi_{(p)}\hat{\tau}_{2}\right]+\cos\theta_{(p)}\hat{\tau}_{3}$
(5)

and
$\hat{g}^{A}_{(p)}=-\hat{\tau}_{3}\hat{g}^{R\,\dagger}_{(p)}\hat{\tau}_{3}$
with complex angles $\theta_{(p)}$ and $\chi_{(p)}$.

For a normal lead, $\theta_{p}=0$ such that
$\hat{g}_{p}^{R/A}=\pm\hat{\tau}_{3}$. Assuming that all the superconducting
leads are grounded and have the same gap amplitude and different
superconducting phases $\phi_{p}$, $\Delta_{p}\equiv\Delta e^{i\phi_{p}}$,
their Green functions are obtained with

$\chi_{p}=\phi_{p}\quad\text{and}\quad\theta_{p}=\theta_{S},\quad\text{where}\quad\tan\theta_{S}=\frac{i\Delta}{E+i0^{+}}.$
(6)

The Green functions $\hat{g}^{R/A}$ of the dot are determined using the
respective blocks of Eq. (4), $[\hat{h}^{R/A},\hat{g}^{R/A}]=0$ with
$\hat{h}^{R/A}=E\hat{\tau}_{3}+i\sum_{p}\Gamma_{p}\hat{g}^{R/A}_{p}$. As
$(\hat{h}^{R/A})^{2}$ is proportional to the identity matrix, it is easy to
see that $\hat{g}^{R(A)}=\hat{h}^{R(A)}/\xi^{R(A)}$ with
$\xi^{R(A)}=\pm\text{sign}(E)\sqrt{(\hat{h}^{R(A)})^{2}}$, where the sign
convention is chosen to match the normal state result
$\hat{g}^{R(A)}=\pm\hat{\tau}_{3}$ in the limit of vanishing coupling to the
leads. Note that $\xi^{A}=-(\xi^{R})^{*}$. Thus the dot parameters $\theta$
and $\chi$ are given as

$\sin\theta=\frac{i\Gamma_{\phi}\sin\theta_{S}}{\xi^{R}},\quad\sin\chi=\frac{\sum_{s}\Gamma_{s}\sin\phi_{s}}{\Gamma_{\phi}}$
(7)

with
$\Gamma_{\phi}=\sqrt{\sum_{s,s^{\prime}}\Gamma_{s}\Gamma_{s^{\prime}}\cos\phi_{ss^{\prime}}}$
and

$\xi^{R}=\text{sign}(E)\sqrt{(E+i\Gamma_{N}+i\Gamma_{S}\cos\theta_{S})^{2}-\Gamma_{\phi}^{2}\sin^{2}\theta_{S}}\,,$
(8)

where the sum over leads $s$ is restricted to the superconducting leads,
$\phi_{ss^{\prime}}=\phi_{s}-\phi_{s^{\prime}}$, and $\Gamma_{N}$ and
$\Gamma_{S}$ are the sums of the partial level widths associated with the N
and S leads, respectively. Note that $\chi$ is real.

The Keldysh part of the equilibrium Green function reads
$\hat{g}^{K}_{p}=(\hat{g}^{R}_{p}-\hat{g}^{A}_{p})f_{L0}(E)$, where
$f_{L0}(E)=\tanh(E/{2T})$ (in units where $k_{B}=1$). For a normal lead $n$
biased at voltage $V_{n}$, one obtains
$\hat{g}^{K}_{n}(E)=2\hat{\tau}_{3}[1-2f(E-eV_{n}\hat{\tau}_{3})]$, where
$f(E)$ is the Fermi-Dirac distribution function. It can also be written in the
form $\hat{g}^{K}_{n}(E)=2\hat{\tau}_{3}f_{Ln}(E)+2\hat{\tau}_{0}f_{Tn}(E)$,
where the longitudinal component of the distribution function,
$f_{Ln}(E)=[f(-E-eV_{n})-f(E-eV_{n})]$, is odd in energy and even in voltage,
whereas the transversal component of the distribution function,
$f_{Tn}(E)=1-f(E-eV_{n})-f(-E-eV_{n})$, is even in energy and odd in voltage.

To obtain the Keldysh part of the dot Green function, following Refs. Snyman
and Nazarov (2009); Catelani et al. (2010), we combine the Keldysh component
of the normalization condition,
$\hat{g}^{R}\hat{g}^{K}+\hat{g}^{K}\hat{g}^{A}=0$, and the Keldysh component
of Eq. (4) to find

$\hat{g}^{K}=\frac{1}{\xi^{R}+\xi^{A}}(\hat{h}^{K}-\hat{g}^{R}\hat{h}^{K}\hat{g}^{A}).$
(9)

Using

$\displaystyle\hat{h}^{K}$ $\displaystyle=$ $\displaystyle
i\sum_{p}\Gamma_{p}\hat{g}_{p}^{K}$ (10) $\displaystyle=$
$\displaystyle(\hat{h}^{R}-\hat{h}^{A})f_{L0}+2i\sum_{n}\Gamma_{n}\left[\hat{\tau}_{3}(f_{Ln}-f_{L0})+\hat{\tau}_{0}f_{Tn}\right],$

where the sum over leads $n$ is restricted to the normal leads, and the
identities derived in Appendix A, it can be written in the form

$\displaystyle\hat{g}^{K}=(\hat{g}^{R}-\hat{g}^{A})f_{L}+(\hat{g}^{R}\hat{\tau}_{3}-\hat{\tau}_{3}\hat{g}^{A})f_{T}$
(11)

with

$\displaystyle f_{L}$ $\displaystyle=$ $\displaystyle f_{L0}+\frac{C}{{\rm
Im}\,\xi^{R}}\sum_{n}\Gamma_{n}(f_{Ln}-f_{L0}),$ (12) $\displaystyle f_{T}$
$\displaystyle=$ $\displaystyle\frac{1}{C\,{\rm
Im}\,\xi^{R}}\sum_{n}\Gamma_{n}f_{Tn},$ (13)

where

$C=\frac{1+\cos\theta\cos\theta^{*}+\sin\theta\sin\theta^{*}}{\cos\theta+\cos\theta^{*}}.$
(14)

Note that for $|E|<\Delta$, $C$ takes a particularly simple form (see Appendix
A), $C={\rm Im}\,\xi^{R}/\Gamma_{N}$, such that the subgap distribution
function only depends on the distribution functions of the normal leads,

$f_{L}(E)=\frac{1}{\Gamma_{N}}\sum_{n}\Gamma_{n}f_{Ln}(E)\quad{\rm
if}\quad|E|<\Delta.$ (15)

This reflects the fact that no thermalization with the superconducting leads
is possible at these energies. At energies $|E|\gg\Delta$, one recovers the
normal state result, where $f_{L}=\sum_{p}\Gamma_{p}f_{Lp}$ is a weighted sum
of the distribution functions of all the leads.

We now have all the elements necessary to evaluate the currents as given by
Eq. (1). The current to a normal lead takes the form

$I_{n}=\frac{G_{n}}{4e}\int dE\,\tilde{I}_{n}(E)\quad{\rm
with}\quad\tilde{I}_{n}=2(f_{T}-f_{Tn}){\rm Re}(\cos\theta).$ (16)

It is a dissipative current that is odd under flipping the signs of all the
voltages (as can be readily deduced from the expression for $f_{Tn})$.

The current to a superconducting lead can be decomposed into a dissipative
current $I_{s}^{dis}$ and a dissipationless supercurrent $I_{s}^{S}$ Likharev
(1979); Barone and Paterno (1982); Volkov (1995), namely,

$\displaystyle I_{s}$ $\displaystyle=$ $\displaystyle
I_{s}^{S}+I_{s}^{dis}=\frac{G_{s}}{2e}\int
dE\,\left(\tilde{I}_{s}^{S}(E)+\tilde{I}_{s}^{dis}(E)\right)$ (17)

with

$\displaystyle\tilde{I}_{s}^{S}(E)$ $\displaystyle=$
$\displaystyle\frac{\sum_{s^{\prime}}\Gamma_{s^{\prime}}\sin\phi_{s^{\prime}s}}{\Gamma_{\phi}}\left[f_{L}{\rm
Im}(\sin\theta){\rm Re}(\sin\theta_{S})+f_{L0}{\rm Re}(\sin\theta){\rm
Im}(\sin\theta_{S})\right],$ (18) $\displaystyle\tilde{I}_{s}^{dis}(E)$
$\displaystyle=$ $\displaystyle f_{T}\left[{\rm Re}(\cos\theta){\rm
Re}(\cos\theta_{S})-\frac{\sum_{s^{\prime}}\Gamma_{s^{\prime}}\cos\phi_{s^{\prime}s}}{\Gamma_{\phi}}{\rm
Re}(\sin\theta){\rm Re}(\sin\theta_{S})\right].$ (19)

Under a global flip in the sign of the voltages in the normal leads we see
that $f_{L}$ is even where as $f_{T}$ is odd which results in the supercurrent
being even and the dissipative current being odd. Thus the dissipative
currents and supercurrents are conserved separately:

$\sum_{n}I_{n}+\sum_{s}I_{s}^{dis}=0\quad{\rm and}\quad\sum_{s}I_{s}^{S}=0.$
(20)

We further note that the two components in our setup satisfy an additional
symmetry property: the supercurrents and disippative currents are odd and
even, respectively, under a global sign flip of the superconducting phases.
(This result may not survive in the case of a finite-length normal region
Dolgirev et al. (2019b).)

Here we are interested in the supercurrents $I_{s}^{S}$. As $\sin\theta_{S}$
is purely real for $|E|<\Delta$ and purely imaginary for $|E|>\Delta$ (see Eq.
(6)), we can distinguish two contributions in Eq. (18): a subgap contribution
that depends on the voltage-dependent dot distribution function $f_{L}$, and a
contribution from energies outside the gap that depends on the equilibrium
distribution function $f_{L0}$. Thus, using Eq. (15) and the observation that
the integrand is even in energy,

$I_{s}^{S}=\frac{\sum_{s^{\prime}}G_{s}\Gamma_{s^{\prime}}\sin\phi_{s^{\prime}s}}{e\Gamma_{\phi}}\left\\{\sum_{n}\frac{\Gamma_{n}}{\Gamma_{N}}\int_{0}^{\Delta}dE\,f_{Ln}{\rm
Im}(\sin\theta)\sin\theta_{S}-i\int_{\Delta}^{\infty}dE\,f_{L0}{\rm
Re}(\sin\theta)\sin\theta_{S}\right\\}.$ (21)

The contribution from energies outside the gap vanishes when the dot is
perfectly coupled to the superconductors ($\Gamma_{S}\to\infty$) such that a
short ballistic junction is formed. In that case, the supercurrent is carried
by Andreev states only. A finite value of $\Gamma_{S}$ has a similar effect as
a non-zero length or a finite Thouless energy, yielding a continuum
contribution to the supercurrent Kulik (1970); Levchenko et al. (2006). The
presence of this continuum contribution is crucial for the phenomena described
here.

Equation (21) can be further simplified using Eqs. (7) and (8), namely

$I_{s}^{S}=\frac{1}{e}\sum_{s^{\prime}}G_{s}\Gamma_{s^{\prime}}\sin\phi_{s^{\prime}s}\left\\{\sum_{n}\frac{\Gamma_{n}}{\Gamma_{N}}\int_{0}^{\Delta}dE\,f_{Ln}(E)j(E)+\int_{\Delta}^{\infty}dE\,f_{L0}(E)j(E)\right\\}$

with

$j(E)=\frac{{\rm Re}\,\xi^{R}}{|\xi^{R}|^{2}}\sin^{2}\theta_{S}=\frac{{\rm
Re}\,\xi^{R}}{|\xi^{R}|^{2}}\frac{\Delta^{2}}{\Delta^{2}-E^{2}}.$ (22)

(We will see in the next section that the singularity at $|E|=\Delta$ is
integrable.)

As ${\rm Re}\,\xi^{R}(E>0)$ is positive, we make the important observation
that the supercurrent,

$I_{s}^{S}=\frac{1}{e}\sum_{s^{\prime}}G_{s}\Gamma_{s^{\prime}}\sin\phi_{s^{\prime}s}\left\\{\sum_{n}\frac{\Gamma_{n}}{\Gamma_{N}}K_{n}^{<}+K^{>}\right\\},$
(23)

results from a competition between positive subgap contributions determined by

$K_{n}^{<}=\int_{0}^{\Delta}dE\;f_{Ln}(E)j(E)>0$ (24)

and negative contributions from energies outside the gap determined by

$K^{>}=\int_{\Delta}^{\infty}dE\;f_{L0}(E)j(E)<0.$ (25)

As we will see in the following, it is this competition that leads to a
nonequilibrium $\pi$-junction. In particular, we find that, in equilibrium,
when $f_{Ln}=f_{L0}$, the supercurrent is dominated by the positive subgap
contributions $K_{n}^{<}$. Out of equilibrium, the modified distribution
functions $f_{Ln}$ suppress the subgap contributions such that, eventually,
the supercurrent will be dominated by the negative contribution $K^{>}$ from
energies outside of the gap. The resulting sign change signals the transition
to a $\pi$-junction behavior. In the following section, we determine the
supercurrent as a function of the applied voltages in various regimes.

##  III Voltage-dependent supercurrent in the absence of inelastic relaxation

For the main part of this paper, we will consider a specific setup with two
superconducting leads, phase-biased at a phase difference
$\phi_{2}-\phi_{1}=\phi$, and normal leads, voltage-biased at voltages with
the same absolute value, $|V_{n}|=V$. As $f_{Ln}$ is an even function of the
voltage, the distribution functions of all the normal leads are the same, such
that all the $K_{n}$ in Eq. (26) are the same, and the expression for the
supercurrent simplifies to

$I_{1}^{S}=-I_{2}^{S}=\frac{\Gamma_{S}}{eR}\sin\phi(K^{>}+K^{<}),$ (26)

where $K^{<}$ is given by Eq. (24) and $R=(G_{1}+G_{2})/G_{1}G_{2}$. This
formula allows us to numerically evaluate the supercurrent in all parameter
regimes. In the following, to get a better understanding of the results, we
discuss limiting cases where an analytical expression for the supercurrent can
be obtained. Note that the dissipative currents depend on the signs of the
voltages in the normal leads. The dissipative currents in the superconducting
leads vanish in a symmetric setup in which $\Gamma_{N>}=\Gamma_{N<}$, where
$\Gamma_{N>}(\Gamma_{N<})$ is the sum over the $\Gamma_{n}$ of the leads
biased at $+V(-V)$. This follows from the separate conservation of
supercurrents and dissipative currents. In an asymmetric setup, the
dissipative currents would contribute to the measured critical current
Kutchinsky et al. (1999); Seviour and Volkov (2000); Bezuglyi et al. (2003).
However, in the regime, where the normal leads that drive the system out of
equilibrium are weakly coupled to the dot, $\Gamma_{N}\ll\Gamma_{S}$, this
effect is negligible. In the following, we will concentrate on that regime.

While a non-vanishing coupling is necessary to establish the out-of-
equilibrium distribution function $f_{L}$, we can see from Eqs. (15) that
$f_{L}(E<\Delta)$, which enters the expression for the supercurrent, does not
depend on the absolute magnitude of the couplings. Namely, the value of
$\Gamma_{N}$ affects the supercurrent, Eq. (26), only via the spectral current
$j(E)$. As $j(E)$ is non-vanishing in the absence of a coupling to the normal
leads, we may evaluate it at $\Gamma_{N}=0$ to obtain the result in leading
order. Then the expression for $\xi^{R}$ determining $j(E)$ takes the form

$\xi^{R}=\begin{cases}\sqrt{\left(E+\frac{\Gamma_{S}E}{\sqrt{\Delta^{2}-E^{2}}}\right)^{2}-\frac{\Gamma_{\phi}^{2}\Delta^{2}}{\Delta^{2}-E^{2}}},\quad
E<\Delta,\\\ \,\\\
\sqrt{\left(E+i\frac{\Gamma_{S}E}{\sqrt{E^{2}-\Delta^{2}}}\right)^{2}+\frac{\Gamma_{\phi}^{2}\Delta^{2}}{E^{2}-\Delta^{2}}},\quad
E>\Delta.\end{cases}$ (27)

Note that $\xi^{R}$ vanishes for some $E=E_{g}$, where $E_{g}$ satisfies the
equation,

$E_{g}=\frac{\Gamma_{\phi}\Delta}{\Gamma_{S}+\sqrt{\Delta^{2}-E_{g}^{2}}}.$
(28)

$\xi^{R}$ is purely imaginary for $E<E_{g}$, and consequently $j(E<E_{g})$
vanishes. Furthermore, in the interval, $E_{g}<E<\Delta$, $\xi^{R}$ is real
such that the spectral current simplifies to
$j(E)=\\{\xi_{R}[1-(E/\Delta)^{2}]\\}^{-1}$.

To proceed further we will study two limiting cases: weak coupling
$\Gamma_{S}\ll\Delta_{0}$ in Sec. III.1 and strong coupling
$\Gamma_{S}\gg\Delta_{0}$ in Sec. III.2. We start by considering the zero-
temperature case, where $f_{L0}(E>0)=1$ and $f_{L}(E>0)=\Theta(E-eV)$. We
determine the critical current as a function of voltage and, in particular,
determine the voltage $V^{*}$ at which a switch from a conventional junction
to a $\pi$-junction takes place due to the competition between $K^{>}$ and
$K^{<}$. We then consider the effect of finite temperature in Sec. III.3. Here
analytical results can be obtained in the regime $T\lesssim T_{c}$.

###  III.1 Weak coupling $\Gamma_{S}\ll\Delta_{0}$ at $T=0$

At $T=0$, the gap in the leads is $\Delta=\Delta_{0}$. Let us first consider
the contributions to the supercurrent from energies $E>\Delta_{0}$. In that
regime, we find

$\xi^{R}\approx\begin{cases}E,&E-\Delta_{0}\gg\frac{\Gamma_{S}^{2}}{\Delta_{0}},\\\
\left(\frac{\Gamma_{S}}{\sqrt{\Gamma_{S}^{2}-\Gamma_{\phi}^{2}}}+i\frac{\sqrt{\Gamma_{S}^{2}-\Gamma_{\phi}^{2}}}{\sqrt{2\Delta_{0}(E-\Delta_{0})}}\right)\Delta_{0},&E\to\Delta_{0}^{+}.\end{cases}$
(29)

Thus, as $E\to\Delta_{0}^{+}$, the spectral current saturates at
$j(\Delta_{0}^{+})=-\Gamma_{S}\Delta_{0}(\Gamma_{S}^{2}-\Gamma_{\phi}^{2})^{-3/2}$,
and $K^{>}$ can be approximated as

$K^{>}\approx-\int_{\Delta_{0}+{\Gamma_{S}^{2}}/{\Delta_{0}}}^{\infty}dE\;\frac{1}{E(E^{2}-\Delta_{0}^{2})}\approx-\ln\frac{\Delta_{0}}{\Gamma_{S}}$
(30)

with logarithmic accuracy. (Note that the case
$\Gamma_{\phi}\approx\Gamma_{S}$, which is realized when $\phi\approx 2\pi n$
would require more careful consideration. However, as the critical current is
realized at phases $\phi\approx\pi(n+\frac{1}{2})$, we will not detail it
here.)

Let us now turn to the subgap contributions. At weak coupling, the spectral
gap $E_{g}$ is small. Namely, Eq. (28) yields
$E_{g}\approx\Gamma_{\phi}\ll\Delta_{0}$, varying between
$E_{g}=|\Gamma_{S_{1}}-\Gamma_{S_{2}}|$ at $\phi=\pi n$ and $E_{g}=\Gamma_{S}$
at $\phi=\pi(n+1/2)$. Furthermore,

$\xi^{R}\approx\begin{cases}\sqrt{E^{2}-E_{g}^{2}},&\Delta_{0}-E\gg\frac{\Gamma_{S}^{2}}{\Delta_{0}},\\\
\sqrt{\frac{(\Gamma_{S}^{2}-\Gamma_{\phi}^{2})\Delta_{0}}{2(\Delta_{0}-E)}},&E\to\Delta_{0}^{-}.\end{cases}$
(31)

For $\Delta_{0}-V<\frac{\Gamma_{S}^{2}}{\Delta_{0}}$, we can thus approximate

$\displaystyle K^{<}~{}$ $\displaystyle\approx$
$\displaystyle\int_{\max(E_{g},eV)}^{\Delta_{0}-\Gamma_{S}^{2}/\Delta_{0}}dE\;\frac{\Delta_{0}^{2}}{\sqrt{E^{2}-E_{g}^{2}}(\Delta_{0}^{2}-E^{2})}$
$\displaystyle\approx$
$\displaystyle\begin{cases}\ln\frac{\Delta_{0}^{2}}{\Gamma_{S}E_{g}},&eV<E_{g},\\\
\ln\frac{\Delta_{0}\sqrt{\Delta_{0}^{2}-(eV)^{2}}}{\Gamma_{S}V},&eV\gg
E_{g}.\end{cases}$

(For $\phi\to\pi n$, $E_{g}$ has to be replaced by
max$(\Gamma_{\phi},\Gamma_{N})$ in the above formulas.)

With Eqs. (30) and (III.1), we find the equilibrium ($V=0$) supercurrent

$I_{1}^{S,{\rm\,eq}}\approx\frac{\Gamma_{S}}{eR}\ln\frac{\Delta_{0}}{E_{g}}\;\sin\phi.$
(33)

The result describes a conventional junction with current-phase relation
$I^{S}(\phi)=I_{c}\sin\phi$ (neglecting the non-sinusoidal corrections due to
phase-dependence of $E_{g}$ in the logarithm) and critical current $I_{c}^{\rm
eq}=(\Gamma_{S}/eR)\ln(\Delta_{0}/E_{g}(\pi/2))$.

Increasing the voltage in the subgap regime to $eV\gg E_{g}$ (still
$eV<\Delta_{0}$), we obtain

$I_{1}^{S}(V)\approx\frac{\Gamma_{S}}{eR}\ln\frac{\sqrt{\Delta_{0}^{2}-(eV)^{2}}}{V}\;\sin\phi.$
(34)

The prefactor in that expression changes sign at $eV^{*}=\Delta_{0}/\sqrt{2}$,
signaling the transition to a $\pi$-junction. Namely, at $V>V^{*}$, the
current phase relation has the form $I^{S}(\phi)=-I_{c}\sin\phi$ with critical
current $I_{c}=(\Gamma_{S}/eR)\ln(V/\sqrt{\Delta_{0}^{2}/e^{2}-V^{2}})$.

Finally at $eV>\Delta_{0}$, the supercurrent saturates at

$I_{1}^{S,\,>}\approx-\frac{\Gamma_{S}}{eR}\ln\frac{\Delta_{0}}{\Gamma_{S}}\;\sin\phi,$
(35)

describing a $\pi$-junction with critical current
$I_{c}^{>}=(\Gamma_{S}/eR)\ln(\Delta_{0}/\Gamma_{S})$ of the same order of
magnitude as the equilibrium critical current.

Our results are in agreement with the original work by Volkov Volkov (1995).

###  III.2 Strong coupling $\Gamma_{S}\gg\Delta_{0}$ at $T=0$

As in the case of weak coupling, we start by considering the contributions to
the supercurrent from energies $E>\Delta_{0}$. Here

$\xi^{R}\approx E+i\frac{\Gamma_{S}E}{\sqrt{E^{2}-\Delta^{2}}},$ (36)

such that

$K^{>}\approx-\Delta_{0}^{2}\int_{\Delta_{0}}^{\infty}dE\;\frac{1}{E(E^{2}+\Gamma_{S}^{2})}\approx-\frac{\Delta_{0}^{2}}{\Gamma_{S}^{2}}\ln\frac{\Gamma_{S}}{\Delta_{0}}$
(37)

with logarithmic accuracy.

Let us now turn to the subgap contributions. At strong coupling,
$E_{g}\approx\Gamma_{\phi}\Delta_{0}/\Gamma_{S}$, varying from $E_{g}^{\rm
min}=|\Gamma_{S_{1}}-\Gamma_{S_{2}}|\Delta_{0}/\Gamma_{S}$ at
$\phi=2\pi(n+1/2)$ to $E_{g}^{\rm max}=\Delta_{0}$ at $\phi=2\pi n$. Except
for the narrow regime
$\Gamma_{S}-\Gamma_{\phi}\ll\frac{\Delta_{0}^{2}}{\Gamma_{S}}$ corresponding
to phases $|\phi-2\pi n|\ll\frac{\Delta_{0}}{\Gamma_{S}}$ (that is not
relevant for determining the critical current, see below), we find

$\xi^{R}\approx\Gamma_{S}\sqrt{\frac{E^{2}-E_{g}^{2}}{\Delta_{0}^{2}-E^{2}}},$
(38)

and consequently

$K^{<}\approx\frac{\Delta_{0}^{2}}{\Gamma_{S}}\int_{\max(E_{g},eV)}^{\Delta_{0}}dE\;\frac{1}{\sqrt{E^{2}-E_{g}^{2}}\sqrt{\Delta_{0}^{2}-E^{2}}}=\begin{cases}\frac{\Delta_{0}}{\Gamma_{S}}K\left(1-\frac{E_{g}^{2}}{\Delta_{0}^{2}}\right),&eV<E_{g},\\\
\frac{\Delta_{0}}{\Gamma_{S}}F\left(\arcsin\left(\sqrt{\frac{\Delta_{0}^{2}-(eV)^{2}}{\Delta_{0}^{2}-E_{g}^{2}}}\right)|1-\frac{E_{g}^{2}}{\Delta_{0}^{2}}\right),&E_{g}<eV<\Delta_{0},\end{cases}$
(39)

where $K$ and $F$ are the complete and incomplete elliptic integrals of the
first kind, respectively. ( Similarly to the weak coupling case, for $\phi\to
2\pi(n+1/2)$, $E_{g}$ has to be replaced by max$(\Gamma_{\phi},\Gamma_{N})$ in
the above formulas.)

The equilibrium current-phase relation is given as

$I_{1}^{S,{\rm\,eq}}\approx\frac{\Delta_{0}}{eR}\;K\left(\frac{4\Gamma_{S_{1}}\Gamma_{S_{2}}}{\Gamma_{S}^{2}}\sin^{2}\frac{\phi}{2}\right)\sin\phi.$
(40)

The critical current is $I_{c}^{\rm eq}\sim\Delta_{0}/(eR)$; it is realized at
$\phi=\phi_{c}^{\rm eq}\in[0,\pi]$, and it corresponds to a conventional
junction. (For $\Gamma_{S_{1}}=\Gamma_{S_{2}}$, one finds $\phi_{c}^{\rm
eq}\approx 0.59\pi$.)

The fact that $E_{g}$ reaches $\Delta_{0}$ (at $\phi=2\pi n$) and that the
contributions to the equilibrium supercurrent from energies outside the gap
are parametrically smaller than the contributions from subgap energies, as can
be seen by comparing Eqs. (37) and (39), leads to a qualitatively different
scenario for the current reversal compared with the weak-coupling case. The
current-phase relation does not differ from the equilibrium case until $eV$
reaches $E_{g}^{\rm min}$. As the voltage further increases, the phase
$\phi^{*}$ at which $E_{g}(\phi^{*})=eV$ decreases from $\pi$ to $0$. Once it
reaches $\phi_{c}^{\rm eq}$, the critical current starts to decrease (see also
Ref. Bezuglyi et al. (2003)). Analyzing the phase dependence of the
supercurrent around $E_{g}(\phi^{*})=eV$, we conclude that
$I_{c}(V)=I_{1}^{S,{\rm\,eq}}(\phi^{*}(V))$ for $eV\gtrsim E_{g}(\phi_{c}^{\rm
eq})$.

For $\Delta_{0}-eV\ll\Delta_{0}-E_{g}$, we may approximate the incomplete
elliptic integral of the first kind as

$F\left(\arcsin\left(\sqrt{\frac{\Delta_{0}^{2}-(eV)^{2}}{\Delta_{0}^{2}-E_{g}^{2}}}\right)|1-\frac{E_{g}^{2}}{\Delta_{0}^{2}}\right)\approx\sqrt{\frac{\Delta_{0}^{2}-(eV)^{2}}{\Delta_{0}^{2}-E_{g}^{2}}},$

which leads to a current-phase relation of the form

$I_{1}^{S}(V)\approx\frac{\Delta_{0}}{eR}\left(\sqrt{\frac{\Delta_{0}^{2}-(eV)^{2}}{\Delta_{0}^{2}-E_{g}^{2}}}-\frac{\Delta_{0}}{\Gamma_{S}}\ln\frac{\Gamma_{S}}{\Delta_{0}}\right)\sin\phi.$
(41)

Thus the prefactor changes sign when $eV$ reaches
$\sqrt{\Delta_{0}^{2}-(\Delta_{0}^{2}-E_{g}^{2})\frac{\Delta_{0}}{\Gamma_{S}}\ln\frac{\Gamma_{S}}{\Delta_{0}}}$.
In particular at the phase $\phi=\pi/2$, which gives the critical current at
$V>V^{*}$, the sign changes when
$\Delta_{0}-eV\approx\Delta_{0}^{3}/(4\Gamma_{S}^{2})\ln^{2}(\Gamma_{S}/\Delta_{0})$.
A $\pi$-junction is realized once $|I_{1}^{S}(V,\pi/2)|$ exceeds
$I_{1}^{S,{\rm\,eq}}(\phi^{*}(V))\approx(\pi\Delta_{0}/2eR)\phi^{*}(V)$ at
voltages
$\Delta_{0}-eV^{*}\approx\Delta_{0}^{3}/[2(1+\pi)^{2}\Gamma_{S}^{2}]\ln^{2}(\Gamma_{S}/\Delta_{0})$.
As a consequence, the critical current does not vanish at the transition. Such
a behavior is characteristic of junctions with a non-sinusoidal current-phase
relation. Fig. 1 shows current-phase relations at different voltages to
illustrate the above scenario.

Figure 1:  The current-phase relation in the strong coupling case for
different voltages. The equilibrium current-phase relation (dotted line) is
non-sinusoidal. In the symmetric case, $\Gamma_{S_{1}}=\Gamma_{S_{2}}$, shown
here, the equilibrium critical current is realized at $\phi_{c}^{\rm
eq}\approx 0.59\pi$. The critical current starts to decrease once
$eV=E_{g}(\phi_{c}^{\rm eq})$ (dash-dotted line). The transition takes place
when the two extrema at positive and negative values of $I_{1}^{S}$ within the
interval $[0,\pi]$ have the same magnitude (solid line). Thus the critical
current at the transition is finite. For $eV>\Delta_{0}$, the current-phase
relation takes the form $I_{1}^{S,>}=I_{c}^{>}\sin(\phi+\pi)$ (dashed line).

At $eV>\Delta_{0}$, the supercurrent saturates at

$I_{1}^{S,\,>}=-\frac{\Delta_{0}^{2}}{eR\Gamma_{S}}\ln\frac{\Gamma_{S}}{\Delta_{0}}\;\sin\phi.$
(42)

Here the critical current
$I_{c}^{>}=(\Delta_{0}^{2}/eR\Gamma_{S})\ln({\Gamma_{S}}/{\Delta_{0}})$ is
parametrically smaller than the equilibrium supercurrent $I_{c}^{\rm eq}$.

Note that, in the limit $\Gamma_{S}\to\infty$, we recover the result of Ref.
Heikkilä et al. (2002) that the supercurrent gradually decreases with voltage
and vanishes at $V>\Delta_{0}$.

###  III.3 Finite temperature

As a next step, we consider the effects of finite temperature and, in
particular, the regime $T\lesssim T_{c}$. In that regime, the equilibrium
supercurrent is reduced due to the fact that the distribution function $f^{\rm
eq}_{L}=\tanh(E/{2T})$ suppresses the low-energy contributions to the
supercurrent. This also changes the competition between contributions from
energies below and above the gap, and therefore affects the voltage $V^{*}$ at
which the $\pi$-transition takes place. As, close to $T_{c}$, the gap is given
as $\Delta/T_{c}\approx\sqrt{8\pi^{2}/{7\zeta(3)}}\sqrt{1-T/T_{c}}$, the
relation $\Gamma_{S}\gg\Delta$ holds at arbitrary coupling.

As before, let us start by considering the contributions to the supercurrent
from energies outside the gap. We find

$\displaystyle K^{>}$ $\displaystyle\approx$
$\displaystyle-\Delta^{2}\int_{\Delta}^{\infty}dE\;\tanh\frac{E}{2T}\frac{1}{E(E^{2}+\Gamma_{S}^{2})}$
$\displaystyle\approx$
$\displaystyle\begin{cases}-\frac{\pi\Delta^{2}}{4\Gamma_{S}T_{c}}+\frac{7\zeta(3)\Delta^{2}}{4\pi^{2}T_{c}^{2}},&\Gamma_{S}\ll
T_{c}\sim\Delta_{0},\\\
-\frac{\Delta^{2}}{\Gamma_{S}^{2}}\log\left(\frac{\Gamma_{S}}{T_{c}}\right),&\Gamma_{S}\gg
T_{c}\sim\Delta_{0}.\end{cases}$

In weak coupling, the result is obtained using

$\displaystyle\int_{\Delta/(2T_{c})}^{\infty}dx\;\frac{\tanh
x}{x(x^{2}+(\frac{\Gamma_{S}}{2T_{c}})^{2})}$ $\displaystyle\approx$
$\displaystyle\int_{0}^{\infty}dx\;\frac{1}{x^{2}+(\frac{\Gamma_{S}}{2T_{c}})^{2}}$
$\displaystyle+\int_{0}^{\infty}dx\;\frac{\tanh x-x}{x^{3}},$

where the second term has to be kept as, in equilibrium, the first term is
canceled by the contributions from subgap energies (see below).

We now turn to the subgap contributions. As $E\ll T\sim T_{c}$, the
distribution function may be approximated as

$f_{L}(V)\approx\frac{E}{2T_{c}}\frac{1}{\cosh^{2}({eV}/{2T_{c}})},$ (44)

which yields

$\displaystyle K^{<}$ $\displaystyle\approx$
$\displaystyle\frac{\Delta^{2}/{2\Gamma_{S}T_{c}}}{\cosh^{2}({eV}/{2T_{c}})}\int_{E_{g}}^{\Delta}dE\;\frac{E}{\sqrt{E^{2}-E_{g}^{2}}\sqrt{\Delta^{2}-E^{2}}}$
(45) $\displaystyle=$
$\displaystyle\frac{\pi\Delta^{2}}{4\Gamma_{S}T_{c}}\frac{1}{\cosh^{2}({eV}/{2T_{c}})}.$

As at $T=0$, the contributions from energies below and above the gap are
comparable in the weak coupling case, $\Gamma_{S}\ll\Delta_{0}$. In that
regime, the supercurrent is given as

$I_{1}^{S}(V)\approx\frac{\Gamma_{S}}{eR}\frac{\Delta^{2}}{4\pi^{2}T_{c}^{2}}\left[7\zeta(3)-\pi^{3}\frac{T_{c}}{\Gamma_{S}}\tanh^{2}\frac{eV}{2T_{c}}\right]\sin\phi.$
(46)

As the dominant contributions cancel at $V=0$, i.e., as the critical current
is parametrically smaller than the individual terms, the sign change happens
at a small voltage,

$eV^{*}\approx\frac{2}{\pi}\sqrt{\frac{7\zeta(3)\Gamma_{S}T_{c}}{\pi}}.$ (47)

Interestingly the supercurrent at $V\gg V^{*}$ parametrically exceeds the
equilibrium supercurrent, namely $I_{c}^{\rm
eq}\sim\Gamma_{S}\Delta^{2}/(eRT_{c}^{2})$, whereas
$I_{c}^{>}\sim\Delta^{2}/(eRT_{c})$. This enhancement occurs in the small
temperature range $T_{c}-T\ll\Gamma_{S}^{2}/T_{c}$.

By contrast, in the strong coupling case, $\Gamma_{S}\gg\Delta_{0}$, the
subgap contributions dominate in equilibrium. Therefore a large voltage is
needed to achieve the $\pi$-transition. In particular,

$I_{1}^{S}(V)\approx\frac{\Gamma_{S}}{eR}\frac{\pi\Delta^{2}}{4\Gamma_{S}T_{c}}\left[\frac{1}{\cosh^{2}({eV}/{2T_{c}})}-\frac{4T_{c}}{\pi\Gamma_{S}}\ln\frac{\Gamma_{S}}{T_{c}}\right]\sin\phi,$
(48)

yielding

$eV^{*}\approx T_{c}\ln\left(\frac{\Gamma_{S}}{T_{c}}\right).$ (49)

Here the critical current at large voltages,
$I_{c}^{>}=\Delta^{2}/(eR\Gamma_{S})$, is parametrically smaller than the
equilibrium critical current, $I_{c}^{\rm eq}=\pi\Delta^{2}/(4eRT_{c})$.

###  III.4 Numerical results

(a)
---
(b)
Figure 2:  Critical current as a function of the voltage applied to the normal
leads in the absence of inelastic relaxation in the dot. The minimum of the
critical current at finite voltage signals the transition from a 0-junction to
a nonequilibrium $\pi$-junction. The critical currents are normalized by the
(temperature-dependent) critical current in equilibrium. (a) Results for $T=0$
at different values of $\Gamma_{S}$. The critical current saturates at
$V=\Delta_{0}$. (b) Results at $T=0.9T_{c}$. Interestingly, in the weak
coupling limit, the critical current at high voltages largely exceeds the
equilibrium critical current.

To visualize the results, we evaluate the critical current numerically. Figure
2 shows the voltage dependence of the critical current for various coupling
strengths, at $T=0$ [Fig. 2(a)] and at $T\lesssim T_{c}$ [Fig. 2(b)]. The
minimum of the critical current at a finite voltage signals the transition to
a $\pi$-junction. The transition happens in all parameter regimes with the
characteristic voltage $V^{*}$ increasing with coupling strength $\Gamma_{S}$.
Numerically $V^{*}$ is obtained by determining the voltage at which the
maximal current $I_{\rm max}(\phi)$ in the interval $\phi\in[0,\pi]$ changes
sign. Figure 3(a) shows the temperature dependence of $V^{*}$ for various
coupling strengths. The critical current vanishes at the transition in the
weak-coupling limit; it increases with $\Gamma_{S}$ as the current-phase
relation becomes non-sinusoidal. The enhancement of the nonequilibrium
critical current close to $T_{c}$ in the weak-coupling limit, as discussed in
Sec. III.3, is illustrated in Fig. 3(b). Here the dependence on $\Gamma_{N}$
is taken into account as well: the effect is seen to diminish as $\Gamma_{N}$
increases.

(a)
---
(b)
Figure 3:  (a) Dependence of the voltage $V^{*}$, at which the transition from
a conventional junction to a nonequilibrium $\pi$-junction occurs, on
temperature. The behavior is qualitatively different at weak coupling, where
$V^{*}$ decreases with temperature, and strong coupling, where $V^{*}$
increases with temperature. The dotted lines show the analytical results that
were obtained close to $T_{c}$. (b) Critical current enhancement in the
$\pi$-phase in the weak coupling limit close to $T_{c}$. The dotted lines show
the enhancement factor $I_{c}^{>}/I_{c}^{\rm eq}=4T_{c}/(\pi\Gamma_{S})$ in
the limit $\Gamma_{N}\to 0$. A reduction at finite $\Gamma_{N}$ is observed.

##  IV Effects of inelastic relaxation

So far we neglected inelastic relaxation within the dot. Such relaxation
processes, if strong enough, tend to establish a Fermi-Dirac distribution in
the dot with an effective temperature and chemical potential determined by the
coupling to the reservoirs. As the non-equilibrium $\pi$-junction relies on
deviations from a Fermi-Dirac distribution in the dot, it is expected that
strong enough inelastic processes will destroy the effect. Here we show that
nevertheless the $\pi$-junction remains robust in a large regime of
parameters.

To determine the effect of internal relaxation, we have to compare it with the
relaxation to the reservoirs. In the subgap regime, relaxation can take place
only with the normal reservoirs. As we assume that $\Gamma_{N}$ is small, this
is a very slow process and internal relaxation should start playing a role as
soon as the corresponding rate exceeds $\Gamma_{N}$.

A simple way to model internal relaxation is to couple the system to a
fictitious fermionic bath Büttiker (1986) whose temperature and chemical
potential are chosen such that the energy and charge currents between the dot
and the bath vanish. We denote the temperature and chemical potential of this
fictitious bath $T_{b}$ and $V_{b}$, respectively. The coupling between the
dot and the bath is characterized by the rate $\Gamma_{b}$. The advantage of
this description is that it is readily described using the general formulas in
Sec. II, extending the sum over normal leads to include the fictitious bath.
We will assume $\Gamma_{b}\ll\Gamma_{S}$.

The condition for the vanishing of the charge current between the dot and the
bath can be deduced from Eq. (16), namely

$\int dE\;(f_{T}-f_{Tb}){\rm Re}(\cos\theta)=0.$ (50)

The energy current $J$ can be written as

$J_{b}=\frac{G_{b}}{16e}\int dE\;E\tilde{J}_{b}(E)$ (51)

with

$\tilde{J}_{b}=\text{Tr}\left(\hat{g}_{b}^{K}\hat{g}^{A}-\hat{g}^{R}\hat{g}_{b}^{K}-\hat{g}^{K}\hat{g}_{b}^{A}+\hat{g}_{b}^{R}\hat{g}^{K}\right),$
(52)

which, using similar considerations as the ones leading to Eq. (16), yields
the condition

$\int dE\;E(f_{L}-f_{Lb}){\rm Re}(\cos\theta)=0.$ (53)

Using the expressions for the distribution functions in the dot, Eqs. (12) and
(13), the two conditions can be rearranged such that the left-hand side only
depends on the parameters $T_{b}$ and $V_{b}$ of the fictitious bath, whereas
the right-hand side only depends on the parameters $T$ and $V_{n}$ of the
normal reservoirs. For simplicity, we will consider only the case $T=0$ here.
For the setup considered in the previous section, this yields

$\displaystyle\int_{0}^{\infty}dE\;f_{Tb}\left(1-\frac{\Gamma_{b}}{C\,\text{Im}\,\xi^{R}}\right){\rm
Re}(\cos\theta)$ $\displaystyle=$
$\displaystyle(\Gamma_{N_{<}}-\Gamma_{N_{>}})\int_{0}^{eV}dE\;\frac{1}{C\,\text{Im}\,\xi^{R}}{\rm
Re}(\cos\theta),$ (54)
$\displaystyle\\!\\!\\!\\!\\!\\!\int_{0}^{\infty}dE\;E(1-f_{Lb})\left(1-\frac{C\Gamma_{b}}{\text{Im}\,\xi^{R}}\right){\rm
Re}(\cos\theta)$ $\displaystyle=$
$\displaystyle\Gamma_{N}\int_{0}^{eV}dE\;E\frac{C}{\text{Im}\,\xi^{R}}{\rm
Re}(\cos\theta),$ (55)

where the parameters $\xi^{R}$ and $C$ are specified in Sec. II.

In a symmetric setup, $\Gamma_{N_{>}}=\Gamma_{N_{<}}$, the right-hand side of
Eq. (54) vanishes. This imposes $V_{b}=0$ such that $f_{Tb}=0$. We will
concentrate on this case to illustrate the effect of relaxation. Extensions to
an asymmetric case are discussed in Appendix C. We find that the
characteristic voltage $V^{*}$ of the $\pi$-transition depends on the
asymmetry only very weakly. Thus the following results are qualitatively valid
also in the extreme case of only one normal lead.

In the absence of superconductivity, the bath temperature $T_{b}$ is readily
obtained from Eq. (55) setting $\Delta_{0}=0$. In that case,
$T_{b}=\sqrt{3\Gamma_{N}/\pi^{2}(\Gamma_{N}+\Gamma_{S})}eV$. Superconductivity
suppresses relaxation to the superconducting leads at low energies. This leads
to a faster rise of the temperature in the subgap regime. The temperature
obtained by solving Eq. (55) numerically is shown in Fig. 4(a), where we
consider the case $\Gamma_{N}\ll\Gamma_{b}\ll\Gamma_{S}$ for different
strengths of $\Gamma_{S}$. The results can be understood qualitatively as
follows: At small $V$, the S leads play no role in the heat balance process.
They do, however, modify the density of states in the dot. Initially only
states at energies $E<E_{g}$ are accessible in the dot. (Their density of
states is finite when taking into account finite $\Gamma_{N}$ and
$\Gamma_{b}$.) In that regime, $T_{b}=\sqrt{3}eV/\pi$. As temperature
increases, the integral on the left hand side of Eq. (55) will contain
contributions $\sim e^{-T_{b}/E_{g}}$ from energies $E>E_{g}$. Due to the
increased density of states at $E>E_{g}$, their contribution can be shown to
become relevant at $T_{b}\sim E_{g}/\ln(\Gamma_{S}/\Gamma_{b})\ll E_{g}$ and
to lead to a slow-down of the increase in temperature. Once $eV>E_{g}$, the
enhanced density of states becomes accessible in the integral on the right
hand side of Eq. (55). This leads to an enhanced power injection and results
in a sharp increase in $T_{b}$. At $eV=\Delta_{0}$, we can approximate Eq.
(55) as

$\displaystyle\frac{\Gamma_{N}}{\Gamma_{b}}\int_{0}^{\Delta_{0}}dE\;Ef_{Lb}\approx\int_{\Delta_{0}}^{\infty}dE\;E(1-f_{Lb}),$
(56)

yielding $T_{b}\sim\Delta_{0}/\ln(\Gamma_{b}/\Gamma_{N})$, i.e., a temperature
that is almost independent of the coupling to the superconducting leads.
Finally at $eV>\Delta_{0}$, the slope is determined by the normal state result
$T_{b}=\sqrt{3\Gamma_{N}/(\Gamma_{N}+\Gamma_{S})}eV/\pi$.

(a)
---
(b)
Figure 4:  (a) Dependence of the bath temperature on voltage in a symmetric
setup, obtained by numerically solving Eq. (55). As discussed in the text, the
temperature at $eV=\Delta_{0}$ is found to be almost independent of the
coupling strength to the superconducting leads. Note that the subgap
temperature-dependence varies with phase due to the phase-dependence of
$E_{g}$. (b) Critical current as a function of voltage (log-scale) for the
same parameters as in (a), obtained numerically from Eq. (57). Inelastic
relaxation weakens the nonequilibrium effects and pushes the $\pi$-transition
to large voltages.

Once $T_{b}$ is determined, the supercurrent can be computed from Eq. (26).
For the specific case considered here and
$\Gamma_{N},\Gamma_{b}\ll\Gamma_{S}$, it takes the form

$I^{S}=\frac{\Gamma_{S}}{eR}\sin\phi\left[K^{>}+\tilde{K}^{<}(T_{b})\right].$
(57)

with

$\displaystyle\tilde{K}^{<}(T_{b})$ $\displaystyle=$
$\displaystyle\Theta(\Delta_{0}-eV)\frac{\Gamma_{N}}{\Gamma_{N}+\Gamma_{b}}\int_{eV}^{\Delta_{0}}\\!\\!\\!dE\;j(E)$
(58)
$\displaystyle+\frac{\Gamma_{b}}{\Gamma_{N}+\Gamma_{b}}\int_{0}^{\Delta_{0}}\\!\\!\\!dE\;f_{Lb}(E)j(E).$

If $\Gamma_{b}\ll\Gamma_{N}$, $\tilde{K}^{<}(T_{b})=K^{<}$ up to corrections
of order $\Gamma_{b}/\Gamma_{N}$.

As the corrections to $K^{<}$ are positive, they tend to increase $V^{*}$. As
they are small, one may expect that the corrections to $V^{*}$ are small. This
turns out to be true only as long as $eV^{*}$ remains smaller than
$\Delta_{0}$. Once $eV^{*}$ reaches $\Delta_{0}$, it increases rapidly upon
further increasing $\Gamma_{b}$ due to the energy exchange with the
superconducting leads. As seen in Sec. III.2, the transition happens very
close to $\Delta_{0}$ in the strong coupling regime such that small
corrections are sufficient to push $eV^{*}$ up to $\Delta_{0}$. Using Eq.
(41), we can estimate that this happens at
$\Gamma_{b}\sim(\Gamma_{N}\Delta_{0}/\Gamma_{S})\ln(\Gamma_{S}/\Delta_{0})\ll\Gamma_{N}$.

In the following, we will concentrate on the opposite regime,
$\Gamma_{b}\gg\Gamma_{N}$, where

$\displaystyle\tilde{K}^{<}(T_{b})$ $\displaystyle\approx$
$\displaystyle\int_{0}^{\Delta_{0}}\\!\\!\\!dE\;f_{Lb}(E)j(E)$ (59)
$\displaystyle=$
$\displaystyle\int_{E_{g}}^{\Delta_{0}}dE\;\frac{1}{\xi^{R}}\frac{\Delta_{0}^{2}}{\Delta_{0}^{2}-E^{2}}\tanh\frac{E}{2T_{b}},$

corresponding to the equilibrium result, but at finite temperature $T_{b}$.

As long as $T_{b}\ll E_{g}$ the effect of the finite temperature is
negligible. Analytic results can be obtained for $E_{g}\ll T_{b}\ll\Delta_{0}$
(relevant for $\Gamma_{S}\ll\Delta_{0}$ only), where

$\tilde{K}^{<}(T_{b})\approx\Delta_{0}^{2}\int_{T_{b}}^{\Delta_{0}-\Gamma_{S}^{2}/\Delta_{0}}dE\;\frac{1}{E(\Delta_{0}^{2}-E^{2})}\approx\ln\frac{\Delta_{0}^{2}}{\Gamma_{S}T_{b}},$
(60)

as well as for $T_{b}\gg\Delta_{0}$, where

$\tilde{K}^{<}(T_{b})\approx\frac{\Delta_{0}^{2}}{2T_{b}}\int_{0}^{\Delta_{0}-\Gamma_{S}^{2}/\Delta_{0}}\;\frac{dE}{\Delta_{0}^{2}-E^{2}}\approx\frac{\Delta_{0}}{2T_{b}}\ln\frac{\Delta_{0}}{\Gamma_{S}}$
(61)

if $\Gamma_{S}\ll\Delta_{0}$ and

$\tilde{K}^{<}(T_{b})\approx\frac{\Delta_{0}^{2}}{2\Gamma_{S}T_{b}}\int_{0}^{\Delta_{0}}\;\frac{dE}{\sqrt{\Delta_{0}^{2}-E^{2}}}\approx\frac{\pi\Delta_{0}^{2}}{4\Gamma_{S}T_{b}}$
(62)

if $\Gamma_{S}\gg\Delta_{0}$. To find the $\pi$-transition, we have to compare
these results with $K^{>}$ computed in the previous section, namely Eq. (30)
at weak coupling and Eq. (37) at strong coupling. We see that in the weak
coupling limit, $\Gamma_{S}\ll\Delta_{0}$, the transition happens when
$T_{b}^{*}$ is of the order of $\Delta_{0}$. On the other hand, in the strong
coupling limit, $\Gamma_{S}\gg\Delta_{0}$, the transition happens at
$T_{b}^{*}\sim\Gamma_{S}/\ln(\Gamma_{S}/\Delta_{0})\gg\Delta_{0}$. In both
cases, the corresponding voltage $V^{*}$ is larger than $\Delta_{0}/e$, such
that we may use the relation
$eV^{*}\sim\sqrt{\frac{\Gamma_{S}}{\Gamma_{N}}}T_{b}^{*}$. At weak coupling,
this yields

$eV^{*}\sim\sqrt{\frac{\Gamma_{S}}{\Gamma_{N}}}\Delta_{0},$ (63)

whereas at strong coupling we find

$eV^{*}\sim\frac{\Gamma_{S}^{3/2}}{\sqrt{\Gamma_{N}}\ln\frac{\Gamma_{S}}{\Delta_{0}}}.$
(64)

Thus, in both cases, the transition from a conventional junction to a
nonequilibrium $\pi$-junction still occurs, though it is pushed to voltages
$eV^{*}\gg\Delta_{0}$. Figure 5 shows $V^{*}$ as a function of $\Gamma_{b}$
for different strengths of $\Gamma_{S}$. The rapid increase in $V^{*}$ once it
has reached $\Delta_{0}$ is clearly seen.

Figure 5:  The effect of inelastic relaxation in the dot, characterized by the
rate $\Gamma_{b}$ on the voltage $V^{*}$, at which the $\pi$-transition
occurs, for the case $T=0$. $eV^{*}$ is found to increase slowly until it
reaches $\Delta_{0}$ followed by a much faster variation, if $\Gamma_{b}$ is
further increased. At $\Gamma_{b}\gg\Gamma_{N}$, the system is always in the
latter regime.

##  V Discussion and outlook

We showed that the supercurrent through a superconductor-normal dot-
superconductor junction is strongly modified in the presence of quasiparticle
injection via normal leads. The effect is most pronounced in the absence of
inelastic relaxation in the dot, when the dot distribution function is very
far from a Fermi-Dirac distribution. In the weak coupling case, the critical
current may be suppressed to zero due a tiny quasiparticle injection current
$I_{N}\sim\Gamma_{N}/{\rm min}\\{\Gamma_{S},\Delta_{0}\\}I_{c}^{\rm eq}\ll
I_{c}^{\rm eq}$ at moderate voltage $eV<\Delta_{0}$. Further increasing $V$
leads to a revival of the supercurrent, though with an inverted sign of the
current-phase relation, corresponding to a $\pi$-junction. In the strong
coupling case, a similar $\pi$-transition is observed, but due to the non-
sinusoidal current-phase relation a finite critical current remains at the
transition. The origin of this $\pi$-transition can be easily understood in
the short junction setup considered here: the supercurrent is determined by a
competition between subgap processes and processes involving energies outside
the gap with opposite sign. Interestingly we find that, in the weak-coupling
case, the critical current at high voltages deep in $\pi$-junction regime may
largely exceed the equilibrium critical current close to $T_{c}$. It is
straightforward to generalize the results to multiterminal junctions as
discussed in Appendix B.

Internal relaxation in the dot leads to a more Fermi-Dirac like distribution
function. This slows down the suppression of the critical current with
increasing injection voltage. As long as the internal relaxation rate
$\Gamma_{b}\ll\Gamma_{S}$, the $\pi$-transition is robust, but it occurs at
much larger voltage. We expect the transition to be completely suppressed at
$\Gamma_{b}\gg\Gamma_{S}$. Furthermore, in our study, we did not consider
relaxation by phonons – the only external relaxation processes are due to the
currents to the leads. This requires the rate $\Gamma_{N}$ to be not too
small. If the main external relaxation process is due to phonons, we also
expect the $\pi$-transition to be absent. A suppression of the critical
current due to the quasiparticle injection, the weaker the larger the phonon
relaxation rate $\Gamma_{\rm ph}$, should remain.

The effect of quasiparticle injection on the critical current in a variety of
setups has been intensively studied in recent years. Here we see in detail how
a tiny quasiparticle injection current may completely modify the system
properties in a very simple setup. Our study further highlights the importance
of the shape of the quasiparticle distribution function with much stronger
effects for a non Fermi-Dirac shape.

###### Acknowledgements.

We thank T. Jalabert and C. Chapelier for useful discussions. We acknowledge
support from the ANR through Grant No. ANR-16-CE30-0019.

##  Appendix A Derivation of some identities for solving the kinetic equation

To obtain Eq. (11) for the Keldysh component of the dot Green function, we use
the identity

$\hat{\tau}_{3}-\hat{g}^{R}\hat{\tau}_{3}\hat{g}^{A}=C(\hat{g}^{R}-\hat{g}^{A}),$
(65)

with $C$ given by Eq. (14). To derive this identity, we start from the
following parameterization of $\hat{g}^{R/A}$,

$\displaystyle\hat{g}^{R}$ $\displaystyle=$
$\displaystyle\sin\theta(\sin\chi\hat{\tau}_{1}+\cos\chi\hat{\tau}_{2})+\cos\theta\hat{\tau}_{3},$
(66) $\displaystyle\hat{g}^{A}$ $\displaystyle=$
$\displaystyle\sin\theta^{*}(\sin\chi\hat{\tau}_{1}+\cos\chi\hat{\tau}_{2})-\cos\theta^{*}\hat{\tau}_{3},$
(67)

where we used that $\chi$ is real. Using trigonometric identities, we easily
obtain

$\displaystyle\hat{g}^{R}-\hat{g}^{A}$ $\displaystyle=$ $\displaystyle
2\cos\frac{\theta+\theta^{*}}{2}\left[\sin\frac{\theta-\theta^{*}}{2}(\sin\chi\hat{\tau}_{1}+\cos\chi\hat{\tau}_{2})+\cos\frac{\theta-\theta^{*}}{2}\hat{\tau}_{3}\right],$
(68) $\displaystyle\hat{\tau}_{3}-\hat{g}^{R}\hat{\tau}_{3}\hat{g}^{A}$
$\displaystyle=$ $\displaystyle
2\cos\frac{\theta-\theta^{*}}{2}\left[\sin\frac{\theta-\theta^{*}}{2}(\sin\chi\hat{\tau}_{1}+\cos\chi\hat{\tau}_{2})+\cos\frac{\theta-\theta^{*}}{2}\hat{\tau}_{3}\right].$
(69)

Thus, Eq. (65) holds with

$C=\frac{\cos\frac{\theta-\theta^{*}}{2}}{\cos\frac{\theta+\theta^{*}}{2}}=\frac{1+\cos\theta\cos\theta^{*}+\sin\theta\sin\theta^{*}}{\cos\theta+\cos\theta^{*}}.$
(70)

The further identity
$\hat{g}^{R}\hat{\tau}_{3}-\hat{\tau}_{3}\hat{g}^{A}=C(\hat{\tau}_{0}-\hat{g}^{R}\hat{g}^{A})$
follows trivially from the normalization condition $\hat{g}^{2}=1$.

For the evaluation of the current, it is useful to show that $C$ simplifies
for $|E|<\Delta$, resulting in a simple expression for $f_{L}$ in the dot as
given in Eq. (15). For $|E|<\Delta$, $\cos\theta_{S}$ is purely imaginary and
$\sin\theta_{S}$ is purely real, such that, using Eqs. (7) and (8), we can
write

$\cos\theta=\frac{a+i\Gamma_{N}}{\xi^{R}}\quad\text{and}\quad\sin\theta=\frac{ib}{\xi^{R}},$
(71)

with $\xi^{R}=\sqrt{(a+i\Gamma_{N})^{2}-b^{2}}$, where $a$ and $b$ are real
numbers. As a consequence,

$C=\frac{|\xi^{R}|^{2}+a^{2}+\Gamma_{N}^{2}+b^{2}}{2a{\rm
Re}\,\xi^{R}+2\Gamma_{N}{\rm Im}\,\xi^{R}}.$ (72)

Using ${\rm Re}\,\xi^{R}{\rm Im}\,\xi^{R}=a\Gamma_{N}$ and $({\rm
Im}\,\xi^{R})^{2}=-(a^{2}-\Gamma_{N}^{2}-b^{2}-|\xi^{R}|^{2})/2$, it can
easily be shown that Eq. (72) reduces to $C={\rm Im}\,\xi^{R}/\Gamma_{N}$.

##  Appendix B Critical current hypersurfaces in multiterminal junctions

---
Figure 6:  Critical current hypersurfaces (CCH) for a setup with three
superconducting leads at $T=0$. The overall shape of the CCH depends on the
values of all the couplings $\Gamma_{s}$. (a) Weak coupling. At $V^{*}$, the
CCH shrinks to a point. (b) Strong coupling. At $V^{*}$, the CCH remains
finite. The slight variations of $V^{*}$ for different leads are not visible
on this scale.

The general formula (26) can easily be applied to multiterminal junctions with
more than two superconducting leads. In that case, there is not a single
critical current, but one can define a critical current hypersurface (CCH)
Pankratova et al. (2020) in the space of $m-1$ independent supercurrents,
where $m$ is the number of superconducting leads. (The remaining current is
determined by current conservation.) The CCH encloses the hypervolume, where a
non-dissipative supercurrent can flow. According to Eq. (26), the
supercurrents are given as

$I_{s}^{S}=\frac{1}{e}\sum_{s^{\prime}}G_{s}\Gamma_{s^{\prime}}\sin\phi_{s^{\prime}s}\left(K^{<}+K^{>}\right).$
(73)

In addition to the explicit phase dependence, $K^{<}$ depends on the phase
difference through
$E_{g}\propto\Gamma_{\phi}=\sqrt{\sum_{s,s^{\prime}}\Gamma_{s}\Gamma_{s^{\prime}}\cos\phi_{ss^{\prime}}}$.

The same competition between $K^{<}$ and $K^{>}$ that was responsible for the
$\pi$-transition in the two-terminal setup will lead to a non-monotonous
dependence of the hypervolume enclosed by the CCH as a function of the voltage
applied to the normal leads.

In the weak coupling case, $\Gamma_{S}\ll\Delta_{0}$, we saw that the phase-
dependence of $K^{<}$ does not play an important role. Thus, the CCH will
evolve with increasing voltage without changing its shape and shrink to a
point at $V^{*}$ before increasing again. Here $V^{*}$ has the same value as
for the two-terminal case, $eV^{*}=\Delta_{0}/\sqrt{2}$. On the other hand, in
the strong coupling case, $\Gamma_{S}\gg\Delta_{0}$, the phase-dependence of
$K^{<}$ does play an important role. Thus, the shape of the CCH will depend on
voltage. Furthermore, as the critical currents never vanish, the CCH does not
shrink to a point: it reaches a minimum at $V^{*}$ before increasing again. As
$V^{*}$ depends on the non-sinusoidal shape of the current-phase
characteristic, one obtains the same order of magnitude as for the two-
terminal case, $eV^{*}\lesssim\Delta_{0}$, but the minima for the critical
currents corresponding to different leads may happen at slightly different
voltages.

Figure 6 shows some examples of critical current hypersurfaces in a setup with
three superconducting leads. They were obtained by evaluating the
supercurrents using Eq. (26) on a grid of $m-1$ independent phases taking
values in the interval $[-\pi,\pi]$.

(a)
---
(b)
Figure 7:  Dependence of (a) the bath temperature $T_{b}$ and (b) the bath
chemical potential $V_{b}$ on voltage in an extreme asymmetric setup with
$\gamma_{-}=1$. The same parameters as in Fig. 4 were used. The results for
$T_{b}$ at $\gamma_{-}=0$ are shown as dashed lines for comparison. (a)
---
(b)
Figure 8:  (a) Critical current as a function of voltage in an extreme
asymmetric setup with $\gamma_{-}=1$. The same parameters as in Fig. 4 were
used. The results for $T_{b}$ at $\gamma_{-}=0$ are shown as dashed lines for
comparison. (b) Critical voltage $V^{*}$ as a function of $\Gamma_{b}$ at
$\gamma_{-}=1$. The results for $\gamma_{-}=0$ (Fig. 5) are shown as dashed
lines for comparison. Deviations can be seen in the weak coupling case at
intermediate $\Gamma_{b}$.

##  Appendix C Inelastic relaxation in an asymmetric setup

In the main text, we considered the effect of inelastic relaxation only in a
symmetric setup such that the voltage $V_{b}$ of the fictitious bath remained
0. Here we address the question what happens in an asymmetric setup, where
$\gamma_{-}=(\Gamma_{N_{>}}-\Gamma_{N_{<}})/\Gamma_{N}\neq 0$. In that case,
the subgap contributions to supercurrent depend on the distribution function

$\displaystyle
f_{Lb}(E)=\frac{1}{2}\left[\tanh\frac{E-eV_{b}}{2T_{b}}+\tanh\frac{E+eV_{b}}{2T_{b}}\right],$
(74)

as can be seen from Eq. (58). Both $T_{b}$ and $V_{b}$ are expected to
increase with voltage and gradually suppress the low-energy contributions to
the supercurrent, eventually leading to a $\pi$-transition. To get an idea
about their magnitude, let us first discuss the normal case, such that Eqs.
(54) and (55) simplify to

$\displaystyle eV_{b}=\gamma_{-}\frac{\Gamma_{N}}{\Gamma_{S}}V,$ (75)
$\displaystyle T_{b}^{2}\left[-{\rm Li}_{2}(-e^{eV_{b}/T_{b}})-{\rm
Li}_{2}(-e^{-eV_{b}/T_{b}})\right]=\frac{\Gamma_{N}}{\Gamma_{S}}\frac{V^{2}}{2},\qquad$
(76)

where $\rm{Li}_{2}$ is the dilogarithm function and we assumed
$\Gamma_{N}\ll\Gamma_{S}$ as for the main part of this paper. This assumption
ensures that $T_{b}\gg eV_{b}$, such that we can approximate

$T_{b}\approx\sqrt{\frac{3\Gamma_{N}}{\pi^{2}\Gamma_{S}}}eV\left(1-\frac{\pi^{2}\gamma_{-}^{2}}{6}\frac{\Gamma_{N}}{\Gamma_{S}}\right).$
(77)

We conclude that the finite $\gamma_{-}$ only leads to small modifications of
the distribution function compared to the symmetric case, corresponding to a
shift of the temperature of the order $\delta
T_{b}/T_{b}\sim\Gamma_{N}/\Gamma_{S}$.

At $\Gamma_{b}\gg\Gamma_{N}$, the $\pi$-transition at $\gamma_{-}=0$ happens
in the regime where $T_{b}$ is given by the normal state result. Thus the
above considerations are sufficient to conclude that a finite $\gamma_{-}$ has
negligible effect. This is further illustrated in Figures (7) and (8).

## References

  * Glazman and Catelani (2021) L. Glazman and G. Catelani, SciPost Physics Lecture Notes p. 031 (2021).
  * Padurariu and Nazarov (2010) C. Padurariu and Y. V. Nazarov, Physical Review B 81, 144519 (2010).
  * Hays et al. (2021) M. Hays, V. Fatemi, D. Bouman, J. Cerrillo, S. Diamond, K. Serniak, T. Connolly, P. Krogstrup, J. Nygård, A. L. Yeyati, et al., Science 373, 430 (2021).
  * De Simoni et al. (2018) G. De Simoni, F. Paolucci, P. Solinas, E. Strambini, and F. Giazotto, Nature Nanotechnology 13, 802 (2018).
  * Paolucci et al. (2018) F. Paolucci, G. De Simoni, E. Strambini, P. Solinas, and F. Giazotto, Nano Letters 18, 4195 (2018).
  * De Simoni et al. (2019) G. De Simoni, F. Paolucci, C. Puglia, and F. Giazotto, ACS Nano 13, 7871 (2019).
  * Paolucci et al. (2019) F. Paolucci, G. De Simoni, P. Solinas, E. Strambini, C. Puglia, N. Ligato, and F. Giazotto, AVS Quantum Science 1, 016501 (2019).
  * Rocci et al. (2020) M. Rocci, G. De Simoni, C. Puglia, D. D. Esposti, E. Strambini, V. Zannier, L. Sorba, and F. Giazotto, ACS Nano 14, 12621 (2020).
  * Golokolenov et al. (2021) I. Golokolenov, A. Guthrie, S. Kafanov, Y. A. Pashkin, and V. Tsepelin, Nature Communications 12, 1 (2021).
  * Ritter et al. (2021a) M. Ritter, A. Fuhrer, D. Haxell, S. Hart, P. Gumann, H. Riel, and F. Nichele, Nature Communications 12, 1 (2021a).
  * Ritter et al. (2021b) M. Ritter, N. Crescini, D. Haxell, M. Hinderling, H. Riel, C. Bruder, A. Fuhrer, and F. Nichele, arXiv preprint arXiv:2106.01816 (2021b).
  * Basset et al. (2021) J. Basset, O. Stanisavljević, M. Kuzmanović, J. Gabelli, C. Quay, J. Estève, and M. Aprili, Physical Review Research 3, 043169 (2021).
  * Chang and Bagwell (1997) L.-F. Chang and P. F. Bagwell, Physical Review B 55, 12678 (1997).
  * Morpurgo et al. (1998) A. Morpurgo, T. Klapwijk, and B. Van Wees, Applied Physics Letters 72, 966 (1998).
  * Savin et al. (2004) A. Savin, J. P. Pekola, J. Flyktman, A. Anthore, and F. Giazotto, Applied Physics Letters 84, 4179 (2004).
  * Van Wees et al. (1991) B. Van Wees, K.-M. Lenssen, and C. Harmans, Physical Review B 44, 470 (1991).
  * Volkov (1995) A. Volkov, Physical Review Letters 74, 4730 (1995).
  * Baselmans et al. (1999) J. Baselmans, A. Morpurgo, B. Van Wees, and T. Klapwijk, Nature 397, 43 (1999).
  * Wilhelm et al. (1998) F. K. Wilhelm, G. Schön, and A. D. Zaikin, Physical Review Letters 81, 1682 (1998).
  * Yip (1998) S.-K. Yip, Physical Review B 58, 5803 (1998).
  * Heikkilä et al. (2002) T. T. Heikkilä, J. Särkkä, and F. K. Wilhelm, Physical Review B 66, 184513 (2002).
  * Huang et al. (2002) J. Huang, F. Pierre, T. T. Heikkilä, F. K. Wilhelm, and N. O. Birge, Physical Review B 66, 020507 (2002).
  * Baselmans et al. (2002) J. Baselmans, T. Heikkilä, B. Van Wees, and T. Klapwijk, Physical Review Letters 89, 207002 (2002).
  * Samuelsson et al. (2000) P. Samuelsson, J. Lantz, V. Shumeiko, and G. Wendin, Physical Review B 62, 1319 (2000).
  * Dolgirev et al. (2019a) P. E. Dolgirev, M. S. Kalenkov, and A. D. Zaikin, Scientific Reports 9, 1 (2019a).
  * Bulaevskii et al. (1977) L. Bulaevskii, V. Kuzii, and A. Sobyanin, JETP Lett 25, 290 (1977).
  * Buzdin et al. (1982) A. I. Buzdin, L. Bulaevskii, and S. Panyukov, JETP Lett 35, 178 (1982).
  * Ryazanov et al. (2001) V. Ryazanov, V. Oboznov, A. Y. Rusanov, A. Veretennikov, A. A. Golubov, and J. Aarts, Physical Review Letters 86, 2427 (2001).
  * Kontos et al. (2002) T. Kontos, M. Aprili, J. Lesueur, F. Genêt, B. Stephanidis, and R. Boursier, Physical Review Letters 89, 137007 (2002).
  * Cohen et al. (2018) Y. Cohen, Y. Ronen, J.-H. Kang, M. Heiblum, D. Feinberg, R. Mélin, and H. Shtrikman, Proceedings of the National Academy of Sciences 115, 6991 (2018).
  * Draelos et al. (2019) 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 Letters 19, 1039 (2019).
  * Pankratova et al. (2020) N. Pankratova, H. Lee, R. Kuzmin, K. Wickramasinghe, W. Mayer, J. Yuan, M. G. Vavilov, J. Shabani, and V. E. Manucharyan, Physical Review X 10, 031051 (2020).
  * Graziano et al. (2020) G. V. Graziano, J. S. Lee, M. Pendharkar, C. J. Palmstrøm, and V. S. Pribiag, Physical Review B 101, 054510 (2020).
  * Alicea et al. (2011) J. Alicea, Y. Oreg, G. Refael, F. Von Oppen, and M. P. Fisher, Nature Physics 7, 412 (2011).
  * Aasen et al. (2016) D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon, M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flensberg, et al., Physical Review X 6, 031016 (2016).
  * Van Heck et al. (2014) B. Van Heck, S. Mi, and A. Akhmerov, Physical Review B 90, 155450 (2014).
  * Yokoyama and Nazarov (2015) T. Yokoyama and Y. V. Nazarov, Physical Review B 92, 155437 (2015).
  * Riwar et al. (2016) R.-P. Riwar, M. Houzet, J. S. Meyer, and Y. V. Nazarov, Nature Communications 7, 1 (2016).
  * Belzig et al. (1999) W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, Superlattices and Microstructures 25, 1251 (1999).
  * Levchenko et al. (2006) A. Levchenko, A. Kamenev, and L. Glazman, Physical Review B 74, 212509 (2006).
  * Büttiker (1986) M. Büttiker, Physical Review B 33, 3020 (1986).
  * Snyman and Nazarov (2009) I. Snyman and Y. V. Nazarov, Physical Review B 79, 014510 (2009).
  * Catelani et al. (2010) G. Catelani, L. Glazman, and K. Nagaev, Physical Review B 82, 134502 (2010).
  * Likharev (1979) K. Likharev, Reviews of Modern Physics 51, 101 (1979).
  * Barone and Paterno (1982) A. Barone and G. Paterno, _Physics and applications of the Josephson effect_ , vol. 1 (Wiley Online Library, 1982).
  * Dolgirev et al. (2019b) P. E. Dolgirev, M. S. Kalenkov, A. E. Tarkhov, and A. D. Zaikin, Physical Review B 100, 054511 (2019b).
  * Kulik (1970) I. Kulik, Sov. Phys. JETP 30, 944 (1970).
  * Kutchinsky et al. (1999) J. Kutchinsky, R. Taboryski, C. B. Sørensen, J. B. Hansen, and P. E. Lindelof, Physical Review Letters 83, 4856 (1999).
  * Seviour and Volkov (2000) R. Seviour and A. Volkov, Physical Review B 61, R9273 (2000).
  * Bezuglyi et al. (2003) E. Bezuglyi, V. Shumeiko, and G. Wendin, Physical Review B 68, 134506 (2003).