Solution of master equations by fermionic-duality:
Time-dependent charge and heat currents through
an interacting quantum dot proximized by a superconductor

Lara C. Ortmanns1,2,3 , Maarten R. Wegewijs1,2,4 , Janine Splettstoesser3

1 Institute for Theory of Statistical Physics, RWTH Aachen, 52056 Aachen,
Germany

2 JARA-FIT, 52056 Aachen, Germany

3 Department of Microtechnology and Nanoscience (MC2), Chalmers University of
Technology, SE-41296 Göteborg

4 Peter Grünberg Institut, Forschungszentrum Jülich, 52425 Jülich, Germany

August 27, 2024

## Abstract

We analyze the time-dependent solution of master equations by exploiting
fermionic duality, a dissipative symmetry applicable to a large class of open
systems describing quantum transport. Whereas previous studies mostly
exploited duality relations after partially solving the evolution equations,
we here systematically exploit the invariance under the fermionic duality
mapping from the very beginning when setting up these equations. Moreover, we
extend the resulting simplifications –so far applied to the local state
evolution– to non-local observables such as transport currents.

We showcase the exploitation of fermionic duality for a quantum dot with
strong interaction –covering both the repulsive and attractive case–
proximized by contact with a large-gap superconductor which is weakly probed
by charge and heat currents into a wide-band normal-metal electrode. We derive
the complete time-dependent analytical solution of this problem involving non-
equilibrium Cooper pair transport, Andreev bound states and strong
interaction. Additionally exploiting detailed balance we show that even for
this relatively complex problem the evolution towards the stationary state can
be understood analytically in terms of the stationary state of the system
itself via its relation to the stationary state of a dual system with inverted
Coulomb interaction, superconducting pairing and applied voltages.

###### Contents

  1. 1 Introduction
  2. 2 Proximized quantum dot
  3. 3 Transport equations
    1. 3.1 Rate and transport equations
    2. 3.2 Liouville-space formulation
  4. 4 Fermionic duality
  5. 5 Solution using fermionic duality
    1. 5.1 Duality-invariant rate variables
    2. 5.2 Duality-adapted basis
    3. 5.3 Diagonalization of the rate superoperator
    4. 5.4 Stationary observables and their duals
    5. 5.5 Evolving state supervector and duality-invariant physical constraints
    6. 5.6 “Universal” stationary duality relation between expectation values
    7. 5.7 Transport current supercovectors
  6. 6 General analytical result
    1. 6.1 State evolution
    2. 6.2 Transport evolution
    3. 6.3 Components of duality invariants
    4. 6.4 Self-duality and transient transport
  7. 7 Limiting cases
    1. 7.1 Interacting dot without pairing ($\alpha=0$)
    2. 7.2 Proximized dot without interaction ($U=0$)
    3. 7.3 Proximized dot at high bias ($\mu\gg|\alpha|,U,|\delta|,T,\Gamma$)
  8. 8 Conclusion
  9. A Stationary duality and Kolmogorov detailed balance
  10. B Condition for self-duality
    1. B.1 Exact self-duality and its breakdown due to interaction
    2. B.2 Asymptotic self-duality

##  1 Introduction

The solution of dynamical equations of open quantum systems exhibiting many-
body effects is inherently a more difficult problem than for closed systems.
One reason for this is the reduction of the number of symmetries, due to the
contact with the macroscopic, dissipative environment. As a result, it may
seem that in general there is no systematic way to simplify the solution
procedure and the physical analysis of the results of interest. This problem
persists even for weakly coupled systems addressed in prior works on
dissipative symmetries [1, 2] and even when symmetries like detailed balance
can be exploited [3, 4, 5, 6, 7].

Recently, it has been found that _fermionic_ open systems, even when strongly
coupled to wide-band metallic contacts, exhibit an extremely useful “symmetry”
that is quite general. This so-called fermionic duality relation is truly
dissipative in nature due to the explicit involvement of the special fermion-
parity decay rate depending only on the interface properties [8, 9]. Unlike
ordinary symmetries this relation involves a transformation of system
_parameters_ in addition to its dynamical variables. Fermionic duality turns
out to include the so-called PT-symmetry [2] specific to Markovian Lindblad
dynamics and generalizes it to strongly coupled non-Markovian systems. It has
been shown to facilitate stationary as well as time-dependent solutions of
master equations (occupations) and quantum master equations (including
coherences), both in the weak coupling limit [10, 8, 11, 12, 13, 14, 15] and
in regimes of strong coupling [16, 17, 9] and can be combined with detailed
balance [13]. Most importantly, it yields unexpected insights into the
physical parameter dependence of both the time scales and the _amplitudes_ of
the state of a decaying quantum system. These affect local observable
quantities such as charge and energy and their non-local transport currents.
The power of the constraints imposed by fermionic duality was illustrated for
relatively simple cases where it completely dictates –rather than just
restricts– the form of the master equation [11]. Such constraints are
extremely valuable for more complicated systems of practical interest for
which an intuitive analysis of analytical results is essentially ruled out
whenever various physical effects stand in competition. Indeed, such an
analysis of the problem solved here has motivated the present work [18].

Figure 1:  Single-level quantum dot with interaction $U$ and pairing $\alpha$
due to proximity of a superconducting contact with a large gap (infinite, not
indicated). The dot is probed by charge- and heat-currents into a weakly
coupled metal M. Of interest are the indicated transient charge and heat
currents when the dot is initialized, for example, by a change of gate-voltage
to a value (a) far from or (b) close to the superconducting resonance. Here we
consider any given initial state that is a mixed state in the energy
eigenbasis [Eq. (19)], independent of the details of its preparation
considered elsewhere [18].

In the above cited prior works, fermionic duality was applied concretely to
the solution of the master-equation eigenvalue problem by expressing half of
its eigenvectors in relation to the other half, which is assumed to be known
by some standard calculation. It was not systematically exploited to find
suitable variables and a suitable basis that would simplify finding the first
half of the solution, and, more importantly, to automatically bring it into a
form that facilitates its analysis. In this paper we push this idea of
fermionic duality as “generalized hermicity” further: going beyond Ref. [8],
we consider from the very beginning quantities that are invariant under the
duality transformation in a weakly coupled open quantum system and refer to
these as duality invariants. Moreover, we apply this analysis not only to the
transition rates (governing state evolution and expectation values of local
observables) but also to the transport rates (governing non-local observables
monitoring the system). Focusing on a concrete example, we show that we obtain
a compact form for the transient state and for charge- and heat-transport
currents. The obtained new “variables” in which the problem can be expressed
are invaluable for any detailed analysis of the full time dependence of its
solution addressing both the decay rates and the amplitudes in dependence of
the initial state and the observables of interest. This detailed analysis
combined with numerical calculations is presented elsewhere [18]. Here we
focus on the required analytical considerations and calculation of the
results. Remarkably, we find in a duality-adapted basis that a _shifted_
version of the transition-rate matrix with duality-invariant elements takes on
a lower-triangular form. Thus, without diagonalizing, the exact time scales of
the dynamics of the problem can simply be read off. Moreover, the eigenvectors
can be expressed compactly in terms of the decay rates and just a single
additional duality-invariant parameter.

Our example is an interacting single-level quantum dot connected to a
superconducting contact –a proximized quantum dot– probed by a weakly coupled
normal-metal contact, see Fig. 1. So far, applications of fermionic duality
mostly focused on models of interacting quantum dots in weak contact with
featureless _normal_ metals [8, 11, 12], although a possible energy dependence
of the coupling to the metal was addressed [13, 14]. When replacing one of the
contacts by a superconductor, the open quantum system is modified by the
superconducting proximity effect and exhibits a much richer dynamics which is
more difficult to analyze and warrants a duality-based approach. The
understanding of proximized interacting quantum dots has seen much theoretical
progress, see for a review [19]. In particular, in previous work the
stationary charge and heat transport were computed [20], as well as the zero-
[21] and finite-frequency noise [22] and counting statistics [23]. The more
complex driven dynamics of hybrid quantum dots was also analysed and exploited
in the context of pumping [24, 25, 26, 27, 28, 29]. Furthermore, the dynamics
of strongly correlated hybrid quantum dot systems was studied [30, 31, 32] but
mostly resorted to numerical tools.

Master equations governing proximized interacting quantum dots were
microscopically derived using the established real-time approach [33, 34] for
the large-gap limit [35, 36] and for a finite gap [20, 37] including
coherences in the energy-basis. The time-dependent state of the proximized dot
was obtained analytically in Ref. [37], but time-dependent charge and heat
currents were not addressed. Yet, we find these two experimentally accessible
quantities particularly interesting, since they provide experimental access to
the complex dynamic behaviour of the proximized quantum dot state, noting also
the continued interest in heat transport in superconductors [38, 39, 40, 41].
Importantly, for such time-dependent decay the proximity effect requires a
careful consideration of the physical procedure of state preparation which
enters the master equation as an initial condition. Already in the normal case
this warrants a systematic analysis of all possible initial and final gate
voltages [8]. For a proximized system the analysis is more complicated since
the induced pairing leads to two inequivalent ways of initializing the quantum
dot using the gate voltage and requires separate analyses of the distinct
ensuing dynamics. This issue has received little attention so far but is
essential for the physical predictions as shown in Ref. [18]. The duality-
based solution required for that analysis is presented here and applies
independent of the way these initial conditions are realized.

In the present work we show that a duality-based calculation of the time-
dependent transport currents provides some important new insights. In doing so
we emphasize that the relevance of the different decay rates / time scales of
the system evolution [20] is determined by their _amplitudes_ and through
these also by the initial state. One thus needs to understand the full
dependence of the amplitudes on physical parameters. For this purpose, the
transition-rate expressions for the master-equation kernel as standardly used
[20, 37, 42] are not suitable since they do not take the extremely restrictive
fermionic duality into account. It would make the analysis of the problem
easier if one could express the _transient_ evolution towards the stationary
state in terms of expectation values of observables in the _stationary_ state
itself and the decay rates. By inspection of transition-rate expressions this
is not obvious at all. However, it turns out that one can express the
amplitudes of the transients in the expectation values of two duality-dictated
observables in the stationary state of the system and, additionally, their
values in the stationary state of the _dual_ system. These stationary
expectation values provide the key to an exhaustive analysis of the solution
[18]. Furthermore, in combination [13] with the detailed-balance property of
this system [3, 7], we find that the stationary observables of the _dual_
system can in turn even be expressed analytically in terms of actual
stationary observables, thus realizing our initial desire. These functions are
“universal”, i.e., independent of system parameters. We, moreover, show that
the duality-based analysis can be extended to the measurable time-dependent
transport quantities when including duality-invariant parts of the transport
rates. In the limit of vanishing superconducting coupling, we reproduce and
extend the results of Ref. [8], which were obtained by exploiting duality for
a normal-conducting system. Also, part of our compact expressions recover
results of Refs. [20, 37] but in a form more suitable for the decay analysis.
Finally, we emphasize that our results apply not only to quantum dots with
repulsive interaction ($U>0$), but also to effectively attractive ones
($U<0$). Nowadays the latter are also accessible in controlled experiments
with a high interest in the interplay with superconductivity [43, 44, 45, 46].
Here fermionic duality is particularly relevant since one of its key features
is that it inverts the sign of the interaction $U$, thereby connecting the
physically very different behaviour of these two classes of systems [14]. We
show that only at the crossover between these two cases at $U=0$ the system
has a higher symmetry characterized by a _self-duality_ of stationary
observables. This ties Coulomb interaction to the breaking of self-duality
even in the presence of superconducting pairing.

The paper is organized as follows. After introducing the proximized quantum
dot system in Sec. 2, we review the established transport equations in Sec. 3.
We present fermionic duality in Sec. 4 and introduce duality invariance which
we systematically exploit in Sec. 5 to derive the solution and analyse its
general features in Sec. 6. Several limiting cases of the solution are
discussed analytically in Sec. 7 and we conclude in Sec. 8. We use units such
that $k_{\text{B}}=\hbar=1$.

##  2 Proximized quantum dot

We consider a quantum dot described by the Hamiltonian
$H_{\text{D}}^{\prime}=\epsilon N+UN_{\uparrow}N_{\downarrow}$, with a single
level $\epsilon$ controlled by a gate voltage and effective interaction $U$ of
either repulsive or attractive sign. Here $N=N_{\uparrow}+N_{\downarrow}$ is
the electron number operator with $N_{\sigma}=d^{\dagger}_{\sigma}d_{\sigma}$.
The dot is tunnel coupled by
$H_{\text{T}}=\sum_{k\sigma}\sqrt{\Gamma/(4\pi)}(d_{\sigma}^{\dagger}c_{k\sigma}+\text{h.c.})$
to a metallic reservoir
$H_{\text{M}}=\sum_{k\sigma}\omega_{k}c_{k\sigma}^{\dagger}c_{k\sigma}$, which
is held at temperature $T$ and electrochemical potential $\mu$. The energy-
independent dot-metal coupling is assumed weak, $\Gamma\ll T$. The observables
of interest are the electron particle current $I_{N}=\partial_{t}\langle
N_{\text{M}}\rangle$, referred to as charge current for simplicity (the actual
charge current being $-|e|I_{N}$) and the heat current
$I_{Q}=\partial_{t}\langle H_{\text{M}}-\mu N_{\text{M}}\rangle$, both
directed into the metal as in Fig. 1. Here
$N_{\text{M}}=\sum_{k\sigma}c_{k\sigma}^{\dagger}c_{k\sigma}$ is the electron
number operator on the metal. The quantum dot is additionally coupled to a
superconductor at zero electro-chemical potential $\mu_{\text{S}}=0$. In the
limit of large superconducting gap exceeding all energy scales, the influence
of the superconductor on the dot is described [47, 48, 34, 49, 50] by a
Hamiltonian term with a pairing amplitude $\alpha$:

$\displaystyle H_{\text{S}}=-\tfrac{1}{2}\alpha
d^{\dagger}_{\uparrow}d_{\downarrow}^{\dagger}+\text{h.c}\ .$ (1)

We choose the phase $\alpha$ to be _real_ by fixing the gauge on the
superconductor. It could even be chosen positive but it is convenient not to
do so for the duality analysis in Sec. 4. In the resulting Hamiltonian
describing the dot-superconductor system, the proximized dot,

$\displaystyle
H_{\text{D}}=H_{\text{D}}^{\prime}+H_{\text{S}}=\sum_{\tau}E_{\tau}\ket{\tau}\bra{\tau}+E_{1}\sum_{\sigma}\ket{\sigma}\bra{\sigma},$
(2)

the pairing generates an Andreev splitting $\delta_{\text{A}}\geq|\alpha|$ of
the discrete 0- and 2-electron levels [51]. For $\alpha\neq 0$ this splitting
differs from the detuning energy $\delta$ as follows:

$\displaystyle E_{\tau}$
$\displaystyle=\tfrac{1}{2}\big{(}\delta+\tau\delta_{\text{A}}\big{)}\quad,\begin{cases}\delta&=2\epsilon+U\\\
\delta_{\text{A}}&=\sqrt{\delta^{2}+\alpha^{2}}\end{cases}.$ (3)

The corresponding states are hybridized into Andreev states with isospin label
$\tau=\pm$:

$\displaystyle\ket{\tau}=\sqrt{\tfrac{1}{2}\Big{[}1-\tau\frac{\delta}{\delta_{\text{A}}}\Big{]}}\ket{0}-\tau\,\frac{\alpha}{|\alpha|}\sqrt{\tfrac{1}{2}\Big{[}1+\tau\frac{\delta}{\delta_{\text{A}}}\Big{]}}\ket{2}.$
(4)

By including $\alpha/|\alpha|=\text{sign}(\alpha)$ into the Andreev states of
Ref. [35] we ensure these are eigenstates of the Hamiltonian (2) for all
_real_ values of $\alpha$ considered later on.111This labelling has the
advantage that $\tau=-$ labels the ground state for all _real_ $\alpha$, but
it implies a swapping $\bm{|}\tau\bm{)}\to\bm{|}-\tau\bm{)}$ under duality
mapping. The latter is in fact desired, see Sec. 4, and is an advantage over
the alternative labelling
$\ket{\tau^{\prime}}\equiv\ket{\tau\alpha/|\alpha|}$, which leaves the states
invariant.The 1-electron spin states $\ket{\sigma}$ of the dot remain
unaffected by the pairing and have energy

$\displaystyle E_{1}=\epsilon=\tfrac{1}{2}\big{(}\delta-U\big{)}.$ (5)

Although $\mu$ can be arbitrary relative to the induced _pairing gap_
$|\alpha|$ on the dot, the considered model only describes non-equilibrium
physics for bias $\mu$ within the _superconductor gap_ which is taken infinite
here. In this description, quasiparticles of the superconductor are
energetically inaccessible. This implies that coherent Cooper pair transfer
($N=0,2$ superpositions) within the dot-superconductor system is not
internally damped and the only source of dissipation is the weak coupling to
the metal probe M. We furthermore focus on the regime considered in Ref. [35]
where only the occupation probabilities $\rho_{\tau}$ of the Andreev states
$\tau=\pm$ together with the odd-parity occupation $\rho_{1}$ need to be
considered, i.e. at all times the proximized dot is in a _mixture_ of energy
states. This is valid provided the induced pairing dominates the dissipation
by the metal

$\displaystyle|\alpha|\gg\Gamma.$ (6)

Of particular interest is the competition of strong Coulomb interaction and
strong pairing in the regime $|U|,|\alpha|\gg T\gg\Gamma$ in the full time-
dependent transport noting that prior results concern stationary charge [35]
and heat-transport [20]. We address this for arbitrary physical parameters and
arbitrary initial energy-mixture described independently by a density operator
denoted $\rho_{0}$. In fact, the initial state $\rho_{0}$ depends on these
parameters but in a way that varies with the choice of physical preparation
procedure. When analysing and comparing the state and transport dynamics
relevant to possible specific experiments this requires separate attention and
is considered elsewhere [18].

##  3 Transport equations

###  3.1 Rate and transport equations

Given the parameters and an initial state of the proximized dot $\rho_{0}$,
diagonal in the energy basis, the transport quantities of interest can be
obtained from the diagonal elements of the time-evolved state denoted
$\rho(t)$. These elements, denoted by $\rho_{\pm}(t)$ and $\rho_{1}(t)$, are
determined by rate equations derived in Ref. [35] from the microscopic model
of Sec. 2:

$\displaystyle\frac{d\rho_{\tau}}{dt}$ $\displaystyle=W_{\tau
1}\rho_{1}-W_{1\tau}\rho_{\tau},$ $\displaystyle\frac{d\rho_{1}}{dt}$
$\displaystyle=\sum_{\tau}W_{1\tau}\rho_{\tau}-\sum_{\tau}W_{\tau 1}\rho_{1}.$
(7)

Here $\rho_{\tau}$ is the occupation probability of the even-parity state
$\ket{\tau}$ with $\tau=\pm$ and $\rho_{1}$ is the occupation – summed over
spin $\sigma=\uparrow,\downarrow$ – of the odd-parity states. The formula for
the charge current into the metal derived in Ref. [35] is

$\displaystyle
I_{N}=\sum_{\eta\tau}\eta\big{[}W_{1,\tau}^{\eta}\rho_{\tau}+W_{\tau,1}^{\eta}\rho_{1}\big{]}$
(8)

where $\eta=\pm$ counts the number of electrons transferred to the metal (in
units of $-e$). Counting instead the energy transferred, gives a similar
formula for the energy current

$\displaystyle
I_{E}=-\sum_{\eta\tau}\big{[}(E_{1}-E_{\tau})W_{1,\tau}^{\eta}\rho_{\tau}+(E_{\tau}-E_{1})W_{\tau,1}^{\eta}\rho_{1}\big{]}$
(9)

from which the heat current $I_{Q}=I_{E}-\mu I_{N}$ into the metal follows.
The transition rates in the master equation (7) read

$\displaystyle W_{1,\tau}=\sum_{\eta}W_{1,\tau}^{\eta},\qquad
W_{\tau,1}=\sum_{\eta}W_{\tau,1}^{\eta},$ (10)

whereas the transport rates featuring in the current formulas (8)-(9),

$\displaystyle W_{1,\tau}^{\eta}$
$\displaystyle=\phantom{\tfrac{1}{2}}\Gamma_{\eta\tau}f^{-\eta}(E_{\eta,\tau}-\mu),$
$\displaystyle W_{\tau,1}^{\eta}$
$\displaystyle=\tfrac{1}{2}\Gamma_{\bar{\eta}\tau}f^{+\bar{\eta}}(E_{\bar{\eta},\tau}-\mu),$
(11)

explicitly keep track of whether an electron, is transferred to ($\eta=-$) or
from ($\eta=+$) the proximized dot with probability
$f^{-\eta}(\omega)=(e^{-\eta\omega/T}+1)^{-1}$. Importantly, spin degrees of
freedom have already been eliminated being incorporated into the factor $1/2$
distinguishing the rates (11). Thus, $\eta$ indicates the direction relative
to the positive particle current into the metal as in Fig. 1 and we denote the
opposite by $\bar{\eta}\equiv-\eta$. The addition energies are denoted by222
Our notation merely relabels Eq. (2.6) of Ref. [35],
$E_{\eta,\tau}^{\text{\cite[cite]{[\@@bibref{}{Sothmann2010Sep}{}{}]}}}\equiv\eta(E_{\eta\tau}-E_{1})=E_{\eta,(\eta\tau)}$
by setting $\tau\to\eta\tau$.

$\displaystyle
E_{\eta,\tau}=\eta\Big{(}E_{\tau}-E_{1}\Big{)}=\tfrac{1}{2}\Big{(}\eta\tau\delta_{\text{A}}+\eta
U\Big{)}$ (12)

and the effective rates leading to these transitions on the proximized dot are
proportional to the probabilities of $\ket{0}$ and $\ket{2}$ in the
superposition $\ket{\tau}$ [Eq. (4)]:

$\Gamma_{\eta\tau}=\gamma_{p}\tfrac{1}{2}\Big{(}1+\eta\tau\frac{\delta}{\delta_{\text{A}}}\Big{)}.$
(13)

Importantly, $\gamma_{p}=\Gamma$ denotes the lump sum of the tunnel rates that
arises when eliminating spin. The advantage of the formal analysis presented
next lies in avoiding the consideration of these individual processes ($\eta$,
$\tau$) and in working with (anti-)symmetric combinations of transport rates
dictated by fermionic duality [Eq. (49)], instead.

###  3.2 Liouville-space formulation

Although equations (7)-(8) were derived in Ref. [35] the time-dependent
solutions for transient transport quantities (8)-(9) of interest here were not
given. To obtain the simplest form of these solutions by exploiting
symmetries, it is crucial to start from equations in basis-independent form in
contrast to Eqs. (7)-(9) where the basis was fixed. Using Liouville-space
notation, any operator $x$ is written as a supervector $\bm{|}x\bm{)}=x$ and
its adjoint, the left supercovector
$\bm{(}x\bm{|}=\operatorname{tr}\big{(}x^{\dagger}\bullet\big{)}$ denotes a
function on operators such that
$\bm{(}y\bm{|}x\bm{)}=\operatorname{tr}(y^{\dagger}x)$. The rate equations
refer to a basis of pure Andreev states for the parity $+1$ sector ($N=0,2$)
and a uniform spin-mixed state for the parity $-1$ sector ($N=1$),
respectively,

$\displaystyle\bm{|}\tau\bm{)}\equiv\ket{\tau}\bra{\tau},\quad\tau=\pm,\qquad\bm{|}1\bm{)}\equiv\tfrac{1}{2}\sum_{\sigma}\ket{\sigma}\bra{\sigma}.$
(14)

Note that $\bm{(}\tau\bm{|}\tau\bm{)}=1$ but
$\bm{(}1\bm{|}1\bm{)}=\tfrac{1}{2}$ is not normalized due to the elimination
of spin. As a result, the completeness relation for the identity superoperator
reads

$\displaystyle\mathcal{I}$
$\displaystyle=\sum_{\tau}\bm{|}\tau\bm{)}\bm{(}\tau\bm{|}+2\bm{|}1\bm{)}\bm{(}1\bm{|}.$
(15)

Applied to a state, one obtains its expansion in terms of probabilities
$\rho_{\lambda}\equiv\bm{(}\lambda\bm{|}\rho\bm{)}/\bm{(}\lambda\bm{|}\lambda\bm{)}$

$\displaystyle\bm{|}\rho\bm{)}=\sum_{\tau=\pm}\rho_{\tau}\bm{|}\tau\bm{)}+\rho_{1}\bm{|}1\bm{)}.$
(16)

The rate equations (7) can now be written as a basis-independent master
equation

$\frac{d}{dt}\bm{|}\rho(t)\bm{)}=W\bm{|}\rho(t)\bm{)},\qquad\bm{|}\rho(0)\bm{)}=\bm{|}\rho_{0}\bm{)},$
(17)

for the state supervector when using the rate superoperator

$\displaystyle W=$
$\displaystyle-\sum_{\tau}W_{1,\tau}\bm{|}\tau\bm{)}\bm{(}\tau\bm{|}+\sum_{\tau}W_{1,\tau}\bm{|}1\bm{)}\bm{(}\tau\bm{|}+\sum_{\tau}2W_{\tau,1}\bm{|}\tau\bm{)}\bm{(}1\bm{|}-\sum_{\tau}2W_{\tau,1}\bm{|}1\bm{)}\bm{(}1\bm{|},$
(18)

as one verifies by inserting $\mathcal{I}$ into Eq. (17) and using Eq. (14).
To solve the evolution problem we have to diagonalize $W$ and explicitly
express the formal solution

$\displaystyle\bm{|}\rho(t)\bm{)}=e^{Wt}\bm{|}\rho_{0}\bm{)}$ (19)

in the _right_ eigenvectors of $W$. The final step is to compute the
amplitudes as functions of the initial state $\bm{|}\rho_{0}\bm{)}$.
Similarly, the charge current (8) and energy current (9) can be expressed as
scalar products with two basis-independent left supervectors:

$\displaystyle I_{N}(t)$ $\displaystyle=\bm{(}I_{N}\bm{|}\rho(t)\bm{)},$
$\displaystyle\bm{(}I_{N}\bm{|}=\sum_{\eta\tau}\eta\Big{[}W_{1,\tau}^{\eta}\bm{(}\tau\bm{|}+2W_{\tau,1}^{\eta}\bm{(}1\bm{|}\Big{]},$
(20a) $\displaystyle I_{E}(t)$ $\displaystyle=\bm{(}I_{E}\bm{|}\rho(t)\bm{)},$
$\displaystyle\bm{(}I_{E}\bm{|}=-\sum_{\eta\tau}E_{+,\tau}\Big{[}-W_{1,\tau}^{\eta}\bm{(}\tau\bm{|}+2W_{\tau,1}^{\eta}\bm{(}1\bm{|}\Big{]}.$
(20b)

To solve the transport problem, we need to express these in the _left_
eigenvectors of $W$, compute their coefficients and evaluate the scalar
products with the evolving state (19).

##  4 Fermionic duality

To diagonalize the rate superoperator $W$ and analyse the non-trivial
behaviour of the resulting solution, the above form (18) based on _energy
basis_ matrix elements of $W$ ($\propto W_{1\tau}$ or $W_{\tau 1}$) is not a
good starting point (even though it is suitable for the microscopic derivation
of $W$ [35]). It is important to see why this is the case: for a large class
of models including the present system333 Fermionic duality was established in
Ref. [8] for strong _bilinear_ coupling to a flat-band metal [Supplement A.5
loc. cit.] and in Ref. [13] for weak but energy-dependent bilinear coupling.
It holds for any quantum dot Hamiltonian $H_{\text{D}}$ which commutes with
the fermion _parity_ operator, $[H_{\text{D}},p]=0$. It thus allows for adding
pairing terms such as in Eq. (1) to describe large-gap superconductors. We
stress that duality is a property of the _dynamics_ and thus unrelated to the
initial state $\rho_{0}$.  the rate superoperator of the master equation obeys
a fermionic duality relation [8, 15, 9]:

$\displaystyle
W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}=-\big{[}\mathcal{P}\overline{\big{(}W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}\big{)}}\mathcal{P}\big{]}^{\dagger}.$
(21)

This symmetry is truly dissipative since it explicitly features the lump sum
rate $\gamma_{p}$ [Eq. (13)]. As written here, it states that the _shifted_
rate superoperator is _invariant_ under a duality mapping up to a sign. The
duality mapping is a “generalized hermitian adjoint” consisting of three
parts: the ordinary hermitian adjoint ${\dagger}$, the superoperator
$\mathcal{P}=p\,\bullet$ which left-multiplies its argument $\bullet$ with the
fermion-parity operator $p=(-\mathds{1})^{N}$ and the over-bar which denotes
an inversion of parameters. For our model the latter reads

$\displaystyle\overline{X(\epsilon,U,\alpha,\mu)}=X(\bar{\epsilon},\bar{U},\bar{\alpha},\bar{\mu})$
(22)

where $X$ denotes some function and $\bar{x}=-x$ for
$x=\epsilon,U,\alpha,\mu$. Temperature $T$ and coupling $\Gamma$ are left
untouched by this inversion444In our weak coupling limit inversion of $\Gamma$
can be usefully avoided by a convention [8]. Temperature, however, plays a
fundamentally different role in fermionic duality [8, 13, 9], in particular,
the $T\to\infty$ limit.  and were therefore not indicated in Eq. (22).

Unlike earlier discussions of duality, we have to deal with _parameter-
dependent_ bases’ due to the superconducting pairing $\alpha$. For example,
for the energy eigenbasis [Eqs. (4),(14)]:

$\displaystyle\Big{[}\mathcal{P}\overline{\bm{|}\tau\bm{)}}\Big{]}^{\dagger}=\bm{(}-\tau\bm{|},\quad\Big{[}\mathcal{P}\overline{\bm{|}1\bm{)}}\Big{]}^{\dagger}=-\bm{(}1\bm{|},$
(23)

since $\mathcal{P}\bm{|}1\bm{)}=-\bm{|}1\bm{)}$ and
$\mathcal{P}\bm{|}\tau\bm{)}=\bm{|}\tau\bm{)}$ and
$\overline{\bm{|}1\bm{)}}=\bm{|}1\bm{)}$ while the energy eigenstates created
by pairing are swapped $\overline{\bm{|}\pm\bm{)}}=\bm{|}\mp\bm{)}$, when
inverting the signs of $\delta$ and $\alpha$ in Eq. (4). Fermionic duality
makes explicit why working in this basis is neither suitable for finding the
solution nor for performing the analysis of its properties: From Eq. (21) we
see that the rate superoperator $W$ shifted by $\gamma_{p}/2$ is antisymmetric
with respect to the three-fold duality mapping even though $W$ itself is not
(anti-)hermitian $W\neq\pm W^{\dagger}$. This relation immediately implies
that left and right eigenvectors of $W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$
–and thus of $W$– for in general _different_ eigenvalues labelled by their
operators $x$, $y$ are _cross_ -related in pairs

$\displaystyle\Big{[}\mathcal{P}\overline{\bm{|}x\bm{)}}\Big{]}^{\dagger}\propto\bm{(}y\bm{|},\qquad\Big{[}\overline{\bm{(}x\bm{|}}\mathcal{P}\Big{]}^{\dagger}\propto\bm{|}y\bm{)}.$
(24)

Both are equivalent to the operator equation $y=p\bar{x}$. The corresponding
eigenvalues obey

$\displaystyle\lambda_{x}=-\bar{\lambda}_{y}.$ (25)

In order to take advantage of these constraints imposed by the fermionic
duality one should work in an orthogonal basis (i) which has the same property
(24) as the desired eigenvectors, for which (ii) the matrix elements of $W$
transform in the same way as the desired eigenvalues under the scalar duality
mapping, the overbar (22). While the energy basis does have property (i) by
Eq. (23), it fails to have (ii). In particular, the basis-independent duality
(21) in the energy basis implies cross-relations between off-diagonal matrix
elements which involves parameter inversion (22), such as
$W_{1\tau}^{\eta}\propto\bar{W}_{\bar{\tau}1}^{\eta}$ [Eq. (26)]. Instead, in
a duality-adapted basis [Eq. (36)] matrix elements of $W$ are not only cross-
related, but even invariant under the map (22).

##  5 Solution using fermionic duality

Thus guided by duality, we now first identify rate variables which are
_invariant_ under the mapping Eq. (22) up to a sign and then identify cross-
related basis supervectors for which matrix elements of $W$ are proportional
to these variables. Although not smaller in number, these variables
significantly simplify the derivation and analysis of the solution as compared
to a brute-force linear algebra approach. In particular, since they are not
mixed up under a duality mapping they can occur on their own (i.e. without
“counter terms” required by invariance) thus ensuring compact formulas. We
also stress that the brute-force approach cannot achieve such simplifications
since fermionic duality exploits the _functional parameter dependence_ of rate
expressions and not just their mutual linear dependence.

###  5.1 Duality-invariant rate variables

As illustrated in Ref. [8] for a simple, weak-coupling example, fermionic
duality is based on the behaviour of the fermionic reservoir distribution
function $f^{\eta}(\omega)=(e^{\eta\omega}+1)^{-1}=f^{-\eta}(-\omega)$ under
inversion of the energy: $f^{\eta}(\omega)=1-f^{\eta}(-\omega)$. It implies
relations between actual and dual transport rates (11):

$\displaystyle W_{1,\tau}^{\eta}$
$\displaystyle=\Gamma_{\eta\tau}-\bar{W}_{1,\bar{\tau}}^{\eta},$
$\displaystyle 2W_{\tau,1}^{\eta}$
$\displaystyle=\Gamma_{\bar{\eta}\tau}-2\bar{W}_{\bar{\tau},1}^{\eta},$
$\displaystyle W_{1,\tau}^{\eta}$
$\displaystyle=2\bar{W}_{\bar{\tau},1}^{\bar{\eta}}.$ (26)

In the last relation, inverting the direction of the transition from $\tau\to
1$ to $1\to\tau$ inverts both the state ($\bar{\tau}=-\tau$) and the direction
of electron transfer ($\bar{\eta}=-\eta)$. Note how our convention for
inverting indices [Eq. (11) ff.] naturally combines with the notation for
parameter inversion of functions [Eq. (22)]. These relations allow to
eliminate transport / transition rates in one direction in favour of those for
the opposite direction by a sum rule obtained by combining the second and
third relation of (26):

$\displaystyle 2W_{\tau,1}^{\eta}+W_{1,\tau}^{\bar{\eta}}$
$\displaystyle=\Gamma_{\bar{\eta}\tau},$ $\displaystyle
2W_{\tau,1}+W_{1,\tau}$ $\displaystyle=\gamma_{p},$ (27)

where the second expression results from summing over $\eta$. We can thus
focus on parametrizing the rates for one direction, say $\tau\to 1$. Using
relations (26) duality-invariant variables can be obtained as follows by
considering (anti-) symmetry with respect to inversion of the energy-state
($\tau$) or electron-transfer direction ($\eta$) or both. We first decompose
the right hand side of the first sum-rule (27) into even and odd parts, using
Eq. (13) with $\eta,\tau=\pm$,

$\displaystyle\Gamma_{\eta\tau}=\tfrac{1}{2}\big{(}\gamma_{p}+\eta\tau\gamma^{\prime}_{p}\big{)},\qquad\gamma_{p}\equiv\Gamma,\qquad\gamma^{\prime}_{p}\equiv\gamma_{p}\frac{\delta}{\delta_{\text{A}}}.$
(28)

Here $\gamma_{p}$ is the simple lump sum rate of Eq. (13) featuring in the
duality relation (21). The coefficient $\gamma^{\prime}_{p}$ is also
relatively simple, adding only a dependence on the detuning relative to the
pairing,
$\gamma^{\prime}_{p}/\gamma_{p}=\delta/\delta_{\text{A}}=(\delta/|\delta|)/\sqrt{1+(\alpha/\delta)^{2}}$.
Next we decompose the transport rates by making an ansatz for
$W_{1,\tau}^{\eta}$ and then use the first sum rule (27) for
$W_{\tau,1}^{\eta}$

$\displaystyle W_{1,\tau}^{\eta}$
$\displaystyle=\tfrac{1}{2}\Gamma_{\eta\tau}+\tfrac{1}{2}\big{[}\gamma_{C}+\eta\gamma^{\prime}_{c}+\tau\gamma_{s}+\eta\tau\gamma^{\prime}_{S}\big{]},$
(29a) $\displaystyle 2W_{\tau,1}^{\eta}$
$\displaystyle=\tfrac{1}{2}\Gamma_{\bar{\eta}\tau}-\tfrac{1}{2}\big{[}\gamma_{C}-\eta\gamma^{\prime}_{c}+\tau\gamma_{s}-\eta\tau\gamma^{\prime}_{S}\big{]}.$
(29b)

Inserting the first ansatz into the first of Eqs. (26) and summing over
$\eta,\tau=\pm$ as either $\sum_{\eta\tau}$ or $\sum_{\eta\tau}\eta$ or
$\sum_{\eta\tau}\tau$ or $\sum_{\eta\tau}\eta\tau$, we see that indeed the new
variables are “$(\pm)$ invariant” under the scalar duality mapping of physical
parameters:

$\displaystyle\gamma_{p}$ $\displaystyle=+\bar{\gamma}_{p},$
$\displaystyle\qquad\gamma_{C}$ $\displaystyle=-\bar{\gamma}_{C},$
$\displaystyle\qquad\gamma_{s}$ $\displaystyle=+\bar{\gamma}_{s},$ (30a)
$\displaystyle\gamma^{\prime}_{p}$ $\displaystyle=-\bar{\gamma}^{\prime}_{p},$
$\displaystyle\qquad\gamma^{\prime}_{c}$
$\displaystyle=-\bar{\gamma}^{\prime}_{c},$
$\displaystyle\qquad\gamma^{\prime}_{S}$
$\displaystyle=+\bar{\gamma}^{\prime}_{S}.$ (30b)

Therefore these are the appropriate variables in which we should express the
shifted rate superoperator $W+(\gamma_{p}/2)\mathcal{I}$ from the very
beginning to optimally exploit the duality relation (21). The prime notation
is physically motivated: Unprimed ($=$ $\eta$-symmetric) invariants determine
the _state_ evolution (insensitive to electron transfer direction) whereas the
primed ($=$ $\eta$-antisymmetric) invariants enter only into _transport_
quantities (sensitive to electron-transfer direction, cf. Eq. (64b)). The
motivation of subscripts $p$ (parity), $C$ (charge) and $S$ (isospin
polarization)555None of the duality invariants is related to the real spin
which was eliminated already in Eq. (7)-(9).  will become clear below.

Thus, using only duality relations (26), we can already conclude that the most
non-trivial parameter dependence of the solution of the time-dependent problem
must be contained in the 4 invariants $\gamma_{C}$, $\gamma^{\prime}_{c}$,
$\gamma_{s}$, $\gamma^{\prime}_{S}$. They can be explicitly expressed as
linear combinations of transport rates (anti)symmetrized with respect to
$\eta$ and / or $\tau$ [Eq. (45),(49),(74)] but these should be avoided until
_after_ we have solved the problem posed in Sec. 6.2 and turn to analysis of
the parameter dependence.

###  5.2 Duality-adapted basis

Having exploited the duality for the scalar rate coefficients we should also
make use of it for the choice of supervector basis. To this end, we use that
in addition to the mapping of scalars to a dual model (indicated by the
overbar $\overline{\phantom{x}}$ ), the duality additionally involves the
mappings $\bm{|}\bullet\bm{)}^{\dagger}=\bm{(}\bullet\bm{|}$ and
$\mathcal{P}\bm{|}\bullet\bm{)}=\bm{|}p\,\bullet\bm{)}$ involving
super(co)vectors. To decompose the shifted rate superoperator
$W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$ whose support is spanned by the
energy-diagonal supervectors $\\{\bm{|}+\bm{)},\bm{|}-\bm{)},\bm{|}1\bm{)}\\}$
we introduce a new basis of right vectors

$\displaystyle\bm{|}\mathds{1}\bm{)}$
$\displaystyle=\sum_{\tau}\bm{|}\tau\bm{)}+2\bm{|}1\bm{)},$
$\displaystyle\quad\bm{|}p\bm{)}$
$\displaystyle=\sum_{\tau}\bm{|}\tau\bm{)}-2\bm{|}1\bm{)},$
$\displaystyle\quad\bm{|}A\bm{)}$
$\displaystyle\equiv\sum_{\tau}\tau\bm{|}\tau\bm{)}.$ (31)

and corresponding adjoint left vectors $\bm{(}\mathds{1}\bm{|}$,
$\bm{(}A\bm{|}$, and $\bm{(}p\bm{|}$ orthogonal to these with respect to the
scalar product $\bm{(}y\bm{|}x\bm{)}=\operatorname{tr}(y^{\dagger}x)$. Here
$\mathds{1}$ is the identity operator, $p$ is the parity operator [featuring
in the duality relation (21)] and $A$ is the polarization operator of the
Andreev states (superconducting Bloch-vector or isospin, see Eq. (56) later).

Although other choices are possible, this basis follows right away from
duality noting only that any rate superoperator has the identity as a left
zero eigenvector by probability normalization, $\bm{(}\mathds{1}\bm{|}W=0$.
This implies that its _duality transform_ ,
$[\overline{\bm{(}\mathds{1}\bm{|}}\mathcal{P}]^{\dagger}=\bm{|}p\bm{)}$, is a
right eigenvector of $W$. Thus we should include both $\bm{|}\mathds{1}\bm{)}$
and $\bm{|}p\bm{)}$ in the orthogonal basis, fixing the remaining basis vector
to be $\bm{|}A\bm{)}$ when normalized as follows

$\displaystyle\tfrac{1}{4}\bm{(}\mathds{1}\bm{|}\mathds{1}\bm{)}=\tfrac{1}{4}\bm{(}p\bm{|}p\bm{)}=\tfrac{1}{2}\bm{(}A\bm{|}A\bm{)}=1.$
(32)

The identity superoperator projecting on the support of the energy basis is
thus

$\displaystyle\mathcal{I}=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}\bm{(}\mathds{1}\bm{|}+\tfrac{1}{2}\bm{|}A\bm{)}\bm{(}A\bm{|}+\tfrac{1}{4}\bm{|}p\bm{)}\bm{(}p\bm{|}.$
(33)

Importantly –like the sets of biorthogonal left and right eigenvectors of $W$
that we seek [Eq. (24)]– the left and right orthogonal basis vectors are
_cross-_ related by the duality mapping: using the operator relations $pA=A$
and $p^{2}=\mathds{1}$

$\displaystyle\overline{\mathcal{P}\bm{|}\mathds{1}\bm{)}}^{\dagger}$
$\displaystyle=\bm{(}p\bm{|}$
$\displaystyle,\quad\overline{\mathcal{P}\bm{|}A\bm{)}}^{\dagger}$
$\displaystyle=-\bm{(}A\bm{|}$
$\displaystyle,\quad\overline{\mathcal{P}\bm{|}p\bm{)}}^{\dagger}$
$\displaystyle=\bm{(}\mathds{1}\bm{|}$ (34a)
$\displaystyle\overline{\bm{(}\mathds{1}\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=\bm{|}p\bm{)}$
$\displaystyle,\quad\overline{\bm{(}A\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=-\bm{|}A\bm{)}$
$\displaystyle,\quad\overline{\bm{(}p\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=\bm{|}\mathds{1}\bm{)}.$ (34b)

The _parameter-dependent_ observable $A$ has the scalar duality mapping
$\bar{A}=-A$ [ Eq. (31) and (23)]. For later reference we note the inverse
formulas of Eq. (31):

$\displaystyle\bm{|}\tau\bm{)}=\tfrac{1}{4}\Big{[}\bm{|}\mathds{1}\bm{)}+\bm{|}p\bm{)}\Big{]}+\tau\tfrac{1}{2}\bm{|}A\bm{)},\quad\tau=\pm,\qquad\bm{|}1\bm{)}=\tfrac{1}{4}\Big{[}\bm{|}\mathds{1}\bm{)}-\bm{|}p\bm{)}\Big{]},$
(35)

By merely changing to the duality-adapted variables (29) and basis (35) the
form of the rate superoperator (18) shifted by the identity (33) to
$W+(\gamma_{p}/2)\mathcal{I}$ is greatly simplified:

$\displaystyle W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}=$
$\displaystyle\tfrac{1}{2}\gamma_{p}\tfrac{1}{4}\Big{[}\bm{|}\mathds{1}\bm{)}\bm{(}\mathds{1}\bm{|}-\bm{|}p\bm{)}\bm{(}p\bm{|}\Big{]}-\gamma_{C}\tfrac{1}{2}\Big{[}\bm{|}p\bm{)}\bm{(}\mathds{1}\bm{|}+\bm{|}A\bm{)}\bm{(}A\bm{|}\Big{]}-\gamma_{s}\tfrac{1}{2}\Big{[}\bm{|}p\bm{)}\bm{(}A\bm{|}+\bm{|}A\bm{)}\bm{(}\mathds{1}\bm{|}\Big{]}.$
(36)

This form manifests fermionic duality (21) term-by-term: the first and last
superoperators are $(-)$ invariant [by Eq. (34)], while their coefficients are
$(+)$ invariant [Eq. (46)] and vice versa for the middle term. Thus
$W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$ is a $(-)$ invariant explicitly
confirming Eq. (21). Therefore this is the appropriate basis in which to
expand $W+(\gamma_{p}/2)\mathcal{I}$ in order to optimally exploit the duality
(21) for its diagonalization.

###  5.3 Diagonalization of the rate superoperator

By changing to the duality-adapted basis we have eliminated 3 out of 9 matrix
elements and obtained a lower-triangular form with only 3 independent duality-
invariant matrix elements. This simplifies the calculation of the diagonal
form in several ways:

$\displaystyle W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$
$\displaystyle=\tfrac{1}{2}\gamma_{p}\bm{|}z\bm{)}\bm{(}z^{\prime}\bm{|}-\gamma_{C}\bm{|}c\bm{)}\bm{(}c^{\prime}\bm{|}-\tfrac{1}{2}\gamma_{p}\bm{|}p\bm{)}\bm{(}p^{\prime}\bm{|}.$
(37)

(i) The eigenvalues can be read off from the diagonal of Eq. (36) giving
$\pm\tfrac{1}{2}\gamma_{p}$ and $-\gamma_{C}$ noting the normalization (33).
(ii) The first left basis vector $\bm{(}\mathds{1}\bm{|}$ and the last right
basis vector $\bm{|}p\bm{)}$ are eigenvectors for different eigenvalues. The
two remaining left (right) eigenvectors are then found by standard forward
(backward) recursion for lower-triangular matrices:

$\displaystyle\bm{(}z^{\prime}\bm{|}$ $\displaystyle=\bm{(}\mathds{1}\bm{|},$
(38a) $\displaystyle\bm{(}c^{\prime}\bm{|}$
$\displaystyle=\bm{(}A\bm{|}+\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}+\gamma_{C}}\bm{(}\mathds{1}\bm{|},$
(38b) $\displaystyle\bm{(}p^{\prime}\bm{|}$
$\displaystyle=\tfrac{1}{4}\bm{(}p\bm{|}+\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}-\gamma_{C}}\tfrac{1}{2}\bm{(}A\bm{|}+\Big{(}\frac{\gamma_{C}}{\gamma_{p}}+\frac{\gamma_{s}}{\gamma_{p}}\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}-\gamma_{C}}\Big{)}\tfrac{1}{2}\bm{(}\mathds{1}\bm{|},$
(38c) $\displaystyle\bm{|}z\bm{)}$
$\displaystyle=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}-\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}+\gamma_{C}}\tfrac{1}{2}\bm{|}A\bm{)}-\Big{(}\frac{\gamma_{C}}{\gamma_{p}}-\frac{\gamma_{s}}{\gamma_{p}}\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}+\gamma_{C}}\Big{)}\tfrac{1}{2}\bm{|}p\bm{)},$
(38d) $\displaystyle\bm{|}c\bm{)}$
$\displaystyle=\tfrac{1}{2}\Big{[}\bm{|}A\bm{)}-\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}-\gamma_{C}}\bm{|}p\bm{)}\Big{]},$
(38e) $\displaystyle\bm{|}p\bm{)}$
$\displaystyle=\bm{|}(-\mathds{1})^{N}\bm{)}.$ (38f)

This compact result reveals that the coefficients of the different
eigenvectors in fact have a very similar functional form. Notably, the eigen
_vectors_ are specified by the eigen _values_ , namely the duality invariants
$\pm\tfrac{1}{2}\gamma_{p}$ and $\gamma_{C}$, and _a single_ additional
parameter, the duality invariant $\gamma_{s}$. (iii) Due to their recursive
construction, the normalization of the eigenvectors is automatically fixed by
that of the basis vectors:

$\displaystyle\bm{(}z^{\prime}\bm{|}z\bm{)}=\tfrac{1}{4}\bm{(}\mathds{1}\bm{|}\mathds{1}\bm{)}=1,\qquad\bm{(}c^{\prime}\bm{|}c\bm{)}=\tfrac{1}{2}\bm{(}A\bm{|}A\bm{)}=1,\qquad\bm{(}p^{\prime}\bm{|}p\bm{)}=\tfrac{1}{4}\bm{(}p\bm{|}p\bm{)}=1.$
(39)

(iv) Finally, the cross-relation of eigenvectors (24) dictated by duality is
manifest,

$\displaystyle\overline{\mathcal{P}\bm{|}z\bm{)}}^{\dagger}$
$\displaystyle=\bm{(}p^{\prime}\bm{|},\,$
$\displaystyle\qquad\overline{\mathcal{P}\bm{|}c\bm{)}}^{\dagger}$
$\displaystyle=-\tfrac{1}{2}\bm{(}c^{\prime}\bm{|},\,$
$\displaystyle\qquad\overline{\mathcal{P}\bm{|}p\bm{)}}^{\dagger}$
$\displaystyle=\bm{(}z^{\prime}\bm{|},$ (40a)
$\displaystyle\overline{\bm{(}z^{\prime}\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=\bm{|}p\bm{)},\,$
$\displaystyle\overline{\bm{(}c^{\prime}\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=-2\bm{|}c\bm{)},\,$
$\displaystyle\overline{\bm{(}p^{\prime}\bm{|}\mathcal{P}}^{\dagger}$
$\displaystyle=\bm{|}z\bm{)},$ (40b)

by merely noting the cross-relation of the duality-adapted basis vectors [Eq.
(34)] and the duality invariance of the scalar coefficients [Eq. (46)]. This
mapping preserves the bi-orthogonality of the left and right eigenvectors. It
also shows how the spectral form (37) of
$W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$ ensures duality $(-)$ invariance: The
mapping cross relates the first and last spectral projectors and it relates
the middle spectral projector to itself666The left and right eigenvectors to
middle eigenvalue $-\gamma_{C}$ cannot be normalized to eliminate the factors
$-1/2$ resp. $-2$ without undesirable effects. The $-$ sign _can_ be
eliminated at the expense of normalization factors which are either complex
(making $c$, $c^{\prime}$ non-self adjoint) or discontinuous in $\delta$, both
unnecessary complications. The factor $1/2$ _can_ be eliminated by
normalization $1/\sqrt{2}$ but the our choice ensures that
$\bm{(}c^{\prime}\bm{|}$ occurs without pre-factors in the charge current
operators (which later on cancel in expectation values). . Together with the
$(+)$ invariance of the eigenvalue $\pm\tfrac{1}{2}\gamma_{p}$ and the $(-)$
invariance of the eigenvalue $\gamma_{C}$ this implies that Eq. (37) simply
inverts its sign as it should by Eq. (21). Note carefully that this inversion
is not simply achieved by inverting the signs of all eigenvalues.

###  5.4 Stationary observables and their duals

We now see that we can entirely express the eigenvectors (38) of $W$ in terms
of _expectation values_ of _physical observables_ $\mathds{1}$, $A$ and $p$
with respect to _two_ stationary states, that of the actual system,
$\bm{|}z\bm{)}$, and of the _dual_ system with inverted parameters,
$\bm{|}\bar{z}\bm{)}$ [Eq. (22)]:

$\displaystyle\bm{(}z^{\prime}\bm{|}$ $\displaystyle=\bm{(}\mathds{1}\bm{|},$
$\displaystyle\qquad\bm{|}z\bm{)}$
$\displaystyle=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
A\rangle_{z}\tfrac{1}{2}\bm{|}A\bm{)}+\langle
p\rangle_{z}\tfrac{1}{4}\bm{|}p\bm{)},$ (41a)
$\displaystyle\bm{(}c^{\prime}\bm{|}$ $\displaystyle=\bm{(}A\bm{|}-\langle
A\rangle_{z}\bm{(}\mathds{1}\bm{|},$ $\displaystyle\bm{|}c\bm{)}$
$\displaystyle=\tfrac{1}{2}\Big{[}\bm{|}A\bm{)}-\langle
A\rangle_{\bar{z}}\bm{|}p\bm{)}\Big{]},$ (41b)
$\displaystyle\bm{(}p^{\prime}\bm{|}$
$\displaystyle=\tfrac{1}{4}\bm{(}p\bm{|}+\langle
A\rangle_{\bar{z}}\tfrac{1}{2}\bm{(}A\bm{|}+\langle
p\rangle_{\bar{z}}\tfrac{1}{4}\bm{(}\mathds{1}\bm{|},$
$\displaystyle\bm{|}p\bm{)}$ $\displaystyle=\bm{|}(-\mathds{1})^{N}\bm{)}.$
(41c)

Here $\langle\bullet\rangle_{z}=\bm{(}\bullet\bm{|}z\bm{)}$ takes the
stationary average and
$\langle\bullet\rangle=\bm{(}\bullet\bm{|}\bar{z}\bm{)}$ takes the average for
the stationary state of the dual system,

$\displaystyle\bm{|}\bar{z}\bm{)}=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
A\rangle_{\bar{z}}\tfrac{1}{2}\bm{|}A\bm{)}+\langle
p\rangle_{z}\tfrac{1}{4}\bm{|}p\bm{)}.$ (42)

In the dual of Eq. (41a), $\bar{A}=-A$ occurs twice such that the signs
cancel. One verifies that applying $\bm{(}\mathds{1}\bm{|}$, $\bm{(}A\bm{|}$,
$\bm{(}p\bm{|}$ to Eq. (41a) indeed gives the expectation values of these
observables as written. These averages can be expressed in terms of only three
invariants

$\displaystyle\langle A\rangle_{z}$
$\displaystyle=-\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}+\gamma_{C}},$
$\displaystyle\,\,\,\,\,\langle p\rangle_{z}$
$\displaystyle=-\frac{2\gamma_{C}}{\gamma_{p}}+\frac{2\gamma_{s}^{2}}{\gamma_{p}(\tfrac{1}{2}\gamma_{p}+\gamma_{C})},$
(43a) $\displaystyle\langle A\rangle_{\bar{z}}$
$\displaystyle=\frac{\gamma_{s}}{\tfrac{1}{2}\gamma_{p}-\gamma_{C}},$
$\displaystyle\langle p\rangle_{\bar{z}}$
$\displaystyle=\frac{2\gamma_{C}}{\gamma_{p}}+\frac{2\gamma_{s}^{2}}{\gamma_{p}(\tfrac{1}{2}\gamma_{p}-\gamma_{C})}.$
(43b)

Clearly, in each column of Eq. (43) the dual expressions transform into each
other by $\gamma_{C}\to-\gamma_{C}$ [Eq. (30)], noting that $p\to\bar{p}=p$
but $A\to\bar{A}=-A$ [Eq. (34)]. The mere form of Eqs. (41) already guarantees
the bi-orthogonality of the eigenvectors with one important exception: the
biorthogonality

$\displaystyle\bm{(}p^{\prime}\bm{|}z\bm{)}=\tfrac{1}{4}\big{[}\langle
p\rangle_{\bar{z}}+\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
A\rangle_{z}\langle A\rangle_{\bar{z}}=0,$ (44)

expresses a non-trivial relation between pairs of stationary expectation
values of the physical system and of the dual system. This is verified to hold
by inserting Eqs. (43), but was overlooked in earlier related work [8], see
Sec. 7.1. That such an additional constraint exists could however be expected
since we have four quantities depending on only three invariants $\gamma_{p}$,
$\gamma_{C}$ and $\gamma_{s}$.

Having fully exploited the duality for diagonalizing
$W+\tfrac{1}{2}\gamma_{p}\mathcal{I}$, we can now return to the rate
superoperator $W$ of interest by shifting back, introducing

$\displaystyle\gamma_{c}\equiv\gamma_{C}+\tfrac{1}{2}\gamma_{p},\qquad\gamma^{\prime}_{s}\equiv\gamma^{\prime}_{S}+\tfrac{1}{2}\gamma^{\prime}_{p}.$
(45)

Although these expressions are no longer strictly invariant [Eq. (30)], they
still obey _shifted_ invariance [Eq. (28)]:

$\displaystyle\gamma_{c}=\gamma_{p}-\bar{\gamma}_{c},\qquad\gamma^{\prime}_{s}=\bar{\gamma}^{\prime}_{s}-\bar{\gamma}^{\prime}_{p}.$
(46)

Inserting Eq. (36) and (37) into
$W=(W+\tfrac{1}{2}\gamma_{p}\,\mathcal{I})-\tfrac{1}{2}\gamma_{p}\,\mathcal{I}$
using Eq. (33) we obtain

$\displaystyle W=\,\,$
$\displaystyle\gamma_{p}\tfrac{1}{4}\Big{[}\bm{|}p\bm{)}\bm{(}\mathds{1}\bm{|}-\bm{|}p\bm{)}\bm{(}p\bm{|}\Big{]}-\gamma_{c}\tfrac{1}{2}\Big{[}\bm{|}p\bm{)}\bm{(}\mathds{1}\bm{|}+\bm{|}A\bm{)}\bm{(}A\bm{|}\Big{]}-\gamma_{s}\tfrac{1}{2}\Big{[}\bm{|}p\bm{)}\bm{(}A\bm{|}+\bm{|}A\bm{)}\bm{(}\mathds{1}\bm{|}\Big{]}$
(47) $\displaystyle=\,\,$
$\displaystyle-\gamma_{c}\bm{|}c\bm{)}\bm{(}c^{\prime}\bm{|}-\gamma_{p}\bm{|}p\bm{)}\bm{(}p^{\prime}\bm{|}.$
(48)

Although the first line also has lower-triangular form as in Eq. (36) with
even an extra zero element, the duality relation (21) has become less clear by
un-doing the shift.777 In the cross relation of spectral projectors the
physically important zero eigenvalue projector is only implicitly defined by
its bi-orthogonality to the non-zero eigenvalue projectors. The zero
eigenvalue thus remains implicit until it appears in the cross-relation of the
spectrum, mapping $(0,\gamma_{c},\gamma_{p})$ to
$(\gamma_{p},\gamma_{p}-\gamma_{c},0)$.  Nevertheless, we can now directly
work with $W$ instead of its shifted version and the appropriate variables are
then $\gamma_{p}$, $\gamma^{\prime}_{p}$, $\gamma_{c}$, $\gamma^{\prime}_{c}$,
$\gamma_{s}$, $\gamma^{\prime}_{s}$. These still have the nice feature that
they are linear combinations of transport rates (anti)symmetrized with respect
to $\eta$ and/or $\tau$:

$\displaystyle\gamma_{p}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}\Big{(}W^{\eta}_{1,\tau}+2W^{\bar{\eta}}_{\tau,1}\Big{)},$
$\displaystyle\qquad\gamma_{c}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}W^{\eta}_{1,\tau},$
$\displaystyle\qquad\gamma_{s}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}\tau W^{\eta}_{1,\tau},$ (49a)
$\displaystyle\gamma^{\prime}_{p}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}\eta\tau\Big{(}W^{\eta}_{1,\tau}+2W^{\bar{\eta}}_{\tau,1}\Big{)},$
$\displaystyle\gamma^{\prime}_{c}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}\eta W^{\eta}_{1,\tau},$
$\displaystyle\gamma^{\prime}_{s}$
$\displaystyle\equiv\tfrac{1}{2}\sum_{\eta\tau}\eta\tau W^{\eta}_{1,\tau}.$
(49b)

The last two columns follow from Eqs. (28)-(29) by inserting the shifted
invariants (45),

$\displaystyle W_{1,\tau}^{\eta}$
$\displaystyle=\tfrac{1}{2}\big{(}\gamma_{c}+\eta\gamma^{\prime}_{c}+\tau\gamma_{s}+\eta\tau\gamma^{\prime}_{s}\big{)},$
(50a) $\displaystyle 2W_{\tau,1}^{\eta}$
$\displaystyle=\tfrac{1}{2}\big{(}(\gamma_{p}-\gamma_{c})+\eta\gamma^{\prime}_{c}-\tau\gamma_{s}+\eta\tau(-\gamma^{\prime}_{p}+\gamma^{\prime}_{s})\big{)}.$
(50b)

and then summing either as $\sum_{\eta\tau}$ or $\sum_{\eta\tau}\eta$ or
$\sum_{\eta\tau}\tau$ or $\sum_{\eta\tau}\eta\tau$. With respect to the
(shifted) invariants, $\gamma_{p},\gamma_{s},\gamma_{c}$, the expectation
values of operators of Eq. (43) read

$\displaystyle\langle A\rangle_{z}$
$\displaystyle=-\frac{\gamma_{s}}{\gamma_{c}},$
$\displaystyle\,\,\,\,\,\langle p\rangle_{z}$
$\displaystyle=1-2\frac{\gamma_{c}^{2}-\gamma_{s}^{2}}{\gamma_{p}\gamma_{c}},$
(51a) $\displaystyle\langle A\rangle_{\bar{z}}$
$\displaystyle=\frac{\gamma_{s}}{\gamma_{p}-\gamma_{c}},$
$\displaystyle\langle p\rangle_{\bar{z}}$
$\displaystyle=1-2\frac{(\gamma_{p}-\gamma_{c})^{2}-\gamma_{s}^{2}}{\gamma_{p}(\gamma_{p}-\gamma_{c})}.$
(51b)

###  5.5 Evolving state supervector and duality-invariant physical constraints

We can now express the evolving state supervector in the eigenbasis of the
evolution

$\displaystyle\bm{|}\rho(t)\bm{)}$
$\displaystyle=e^{Wt}\bm{|}\rho_{0}\bm{)}=\bm{|}z\bm{)}+\bm{|}c\bm{)}e^{-\gamma_{c}t}\bm{(}c^{\prime}\bm{|}\rho_{0}\bm{)}+\bm{|}p\bm{)}e^{-\gamma_{p}t}\bm{(}p^{\prime}\bm{|}\rho_{0}\bm{)}$
(52)

with amplitudes depending on physical expectation values in the initial state
($\rho_{0}$), the stationary state ($z$) _and_ the dual stationary state
($\bar{z}$):

$\displaystyle\bm{(}z^{\prime}\bm{|}\rho_{0}\bm{)}$ $\displaystyle=1,$ (53a)
$\displaystyle\bm{(}c^{\prime}\bm{|}\rho_{0}\bm{)}$ $\displaystyle=\langle
A\rangle_{\rho_{0}}-\langle A\rangle_{z},$ (53b)
$\displaystyle\bm{(}p^{\prime}\bm{|}\rho_{0}\bm{)}$
$\displaystyle=\tfrac{1}{4}\big{[}\langle p\rangle_{\rho_{0}}+\langle
p\rangle_{\bar{z}}\big{]}+\tfrac{1}{2}\langle A\rangle_{\bar{z}}\langle
A\rangle_{\rho_{0}}$ (53c) $\displaystyle=\tfrac{1}{4}\big{[}\langle
p\rangle_{\rho_{0}}-\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
A\rangle_{\bar{z}}\big{[}\langle A\rangle_{\rho_{0}}-\langle
A\rangle_{z}\big{]},$ (53d)

subtracting relation (44) from Eq. (53c) in the last step. At this point it is
relevant to consider the physical constraints on the duality invariants and
the stationary values of duality-adapted observables which will be used later
on.

First, the decay rates of the time-dependent state (52), the (shifted) duality
invariants $\gamma_{c}$ and $\gamma_{p}$, are clearly non-negative since they
are proportional to the sums (49a) of non-negative transition rates of the
master equation. However, their individual non-negativity imposes stronger
physical constraints which also involve the invariant $\gamma_{s}$ which can
be negative like $\gamma^{\prime}_{c}$ and $\gamma^{\prime}_{s}$. These are
not rates but “rate asymmetries” Eq. (49b) which account for various
_differences_ of the transport rates [electron/hole ($\eta$) and state
asymmetry ($\tau$)]. The non-negativity of $W_{1,\tau}$ and $W_{\tau,1}$ [Eq.
(11)] imposes a bound on the negativity of $\gamma_{s}$:

$\displaystyle\gamma_{p}\geq
0,\quad\gamma_{c},\gamma_{p}-\gamma_{c}\geq|\gamma_{s}|.$ (54)

Thus the decay rate $\gamma_{c}$ should not only be non-negative but even
larger than the magnitude of the rate-asymmetry $\gamma_{s}$. Likewise, the
decay rate $\gamma_{p}$ should not only exceed the decay rate $\gamma_{c}$,
but it must do so by more than the magnitude of $\gamma_{s}$. The latter two
conditions are equivalent to one, duality-invariant condition on
$\gamma_{C}=\gamma_{c}-\gamma_{p}/2$:

$\displaystyle|\gamma_{C}|\,=\,|\gamma_{c}-\tfrac{1}{2}\gamma_{p}|\,\,\leq\,\,\big{|}\,|\gamma_{s}|-\tfrac{1}{2}\gamma_{p}\,\big{|},$
(55)

expressing that the sum of transition rates
$2\gamma_{c}=\sum_{\tau}W_{1,\tau}$ is always closer to $\gamma_{p}$ than
their difference $2\gamma_{s}=\sum_{\tau}\tau W_{1,\tau}$. Interestingly,
whether $\gamma_{c}$ dominates $\gamma_{p}/2$ or vice versa is determined by
the sign of $\gamma_{C}$ which can be either $+$ or $-$ [see Eq. (63) below].

Using this we can explicitly see that the stationary supervector (41a) giving
the long-time limit of (52) is a valid physical state for all parameter values
of the model by rewriting it as a sum of two operators with parity $+1$ and
$-1$:

$\displaystyle\bm{|}z\bm{)}=$
$\displaystyle\Big{\\{}\tfrac{1}{2}\big{[}1+\langle
p\rangle_{z}\big{]}\tfrac{1}{4}\big{[}\bm{|}\mathds{1}\bm{)}+\bm{|}p\bm{)}\big{]}+\tfrac{1}{2}\langle
A\rangle_{z}\bm{|}A\bm{)}\Big{\\}}+\tfrac{1}{2}\big{[}1-\langle
p\rangle_{z}\big{]}\tfrac{1}{4}\big{[}\bm{|}\mathds{1}\bm{)}-\bm{|}p\bm{)}\big{]}.$
(56)

Indeed, using the expectation values (51) and the conditions (54) imply

$\displaystyle 0\leq\tfrac{1}{2}[1-\langle p\rangle_{z}]\leq 1,\quad
0\leq\tfrac{1}{2}[1+\langle p\rangle_{z}]\leq 1,\quad|\langle
A\rangle_{z}|\leq\tfrac{1}{2}[1+\langle p\rangle_{z}].$ (57)

The first two conditions ensure that the probabilities
$\tfrac{1}{2}[1\pm\langle p\rangle_{z}]$ of being in the parity $\pm 1$ sector
lie in the range [0,1]. The third condition ensures that conditional on being
in the parity $+1$ sector, the difference of the individual occupations of its
states $\bm{|}\tau\bm{)}$ cannot become too large: For a valid physical state
in the parity $+1$ sector, the length of the (1-dimensional) Bloch vector
$\langle A\rangle_{z}$ should not exceed the radius of the (1-dimensional)
Bloch sphere set by the total probability $\tfrac{1}{2}[1+\langle
p\rangle_{z}]$ of being in that sector. Applying the duality mapping we find
that the dual quantities also obey the corresponding constraints:

$\displaystyle 0\leq\tfrac{1}{2}[1-\langle p\rangle_{\bar{z}}]\leq 1,\quad
0\leq\tfrac{1}{2}[1+\langle p\rangle_{\bar{z}}]\leq 1,\quad|\langle
A\rangle_{\bar{z}}|\leq\tfrac{1}{2}[1+\langle p\rangle_{\bar{z}}].$ (58)

Thus, the dual and actual stationary vector are simultaneously legitimate
states.

Figure 2:  (a) Stationary polarization $\langle A\rangle_{\bar{z}}$ and (b)
parity $\langle p\rangle_{\bar{z}}$ of the dual system as function of the
stationary parity $\langle p\rangle_{z}$ for various polarizations $\langle
A\rangle_{z}$ of the actual system [Eq. (59a)]. (c) Non-linear rescaling
factor $F$ as function of $\langle A\rangle_{z}$ and $\langle p\rangle_{z}$
[Eq. (59b)]. The self-duality condition $F(\langle A\rangle_{z},\langle
p\rangle_{z})=1$ (green 3D curve) holds along the curve $\langle
p\rangle_{z}=\langle A\rangle_{z}^{2}$ (dashed green 2D curve) in base plane
which includes pure states $\bm{|}\pm\bm{)}$ (blue/red points) and the
maximally mixed state (black point) which is the only state for which
$\bm{|}\bar{z}\bm{)}=\bm{|}z\bm{)}$. Corresponding values are shown in (a) and
(b).

###  5.6 “Universal” stationary duality relation between expectation values

The result (51) expresses the stationary observables for the system and the
dual system in the same set of (shifted-)duality invariants. This already
implied the useful relation (44) but one can go one step further and,
remarkably, express the _dual_ polarization and parity in terms of the
_actual_ polarization and parity:

$\displaystyle\langle\bar{A}\rangle_{\bar{z}}=-\langle A\rangle_{\bar{z}}$
$\displaystyle=F\big{(}\langle A\rangle_{z},\langle
p\rangle_{z}\big{)}\cdot\langle A\rangle_{z},\,\,$
$\displaystyle\tfrac{1}{2}\big{[}1+\langle p\rangle_{\bar{z}}\big{]}$
$\displaystyle=F\big{(}\langle A\rangle_{z},\langle
p\rangle_{z}\big{)}\cdot\tfrac{1}{2}\big{[}1+\langle p\rangle_{z}\big{]}.$
(59a) This _stationary duality relation_ thereby explicitly expresses
$\bm{|}\bar{z}\bm{)}$ in terms of $\bm{|}z\bm{)}$ [Eq. 41a] by a rational
function $F$ which is “universal” in the sense that it is _independent_ of the
physical parameters: $\displaystyle F\big{(}\langle A\rangle_{z},\langle
p\rangle_{z}\big{)}\equiv\frac{\tfrac{1}{2}\big{[}1-\langle
p\rangle_{z}\big{]}}{\tfrac{1}{2}\big{[}1+\langle p\rangle_{z}\big{]}-\langle
A\rangle_{z}^{2}}.$ (59b)

The relation (59) is remarkable: given only polarization and parity
_expectation values_ of the actual system with physical parameters
($\epsilon,U,\alpha,\mu,T,\Gamma$) it allows to compute these values for the
system with dual physical parameters ($-\epsilon,-U,-\alpha,-\mu,T,\Gamma$)
without requiring any reference to these parameters: We plot in Fig. 2(a-b)
the (components of the) dual stationary state as a _non-linear_ function of
the (components of the) stationary state of the actual system. In Fig. 2(c) we
plot the scaling factor (59b) which takes non-negative values on the domain
(57) of physically allowed values since there $\tfrac{1}{2}[1+\langle
p\rangle_{z}]\geq|\langle A\rangle_{z}|\geq\langle A\rangle_{z}^{2}$. It can
be arbitrarily large (even diverging at $\langle p\rangle_{z}=-1$ and $\langle
A\rangle_{z}=0$) while always producing legitimate values $\langle
p\rangle_{\bar{z}},\langle A\rangle_{\bar{z}}\in[-1,1]$. This is possible due
to the non-linearity of the duality relation (59), i.e., the dependence of $F$
on $\langle A\rangle_{z}$ and $\langle p\rangle_{z}$.

The stationary duality (59) has been derived previously in Ref. [13] for
master equations obeying both fermionic duality and detailed balance (which is
the case here). This is highlighted in one of the two derivations of Eq. (59)
given in App. A via the energy eigenbasis. However, here, in our duality-
adapted basis this relation acquires a simpler form of a rescaling by a single
factor $F\big{(}\langle A\rangle_{z},\langle p\rangle_{z}\big{)}$ because we
have written relation (59a) in terms of the probability
$\tfrac{1}{2}[1+\langle p\rangle_{z}]$ for the parity $+1$ sector and the
polarization $\langle A\rangle_{z}$ in that sector. This nicely connects to
our earlier result that the stationary state for the actual system
$\bm{|}z\bm{)}$ and the dual system $\bm{|}\bar{z}\bm{)}$ are simultaneously
legitimate states: conditions (57) and (58) are equivalent since the factor
$F$ drops out in $|\langle A\rangle_{z}|/[\tfrac{1}{2}[1+\langle
p\rangle_{z}]]=|\langle A\rangle_{\bar{z}}|/[\tfrac{1}{2}[1+\langle
p\rangle_{\bar{z}}]]\leq 1$: The duality mapping between stationary states
preserves the magnitude of the (1-dimensional) Bloch vector relative to the
(1-dimensional) Bloch sphere in the parity $+1$ sector.

Relation (59) also highlights the special situation of _self-duality_ of
stationary observables drawn in green in Fig. 2. This is defined as:

$\displaystyle\langle\bar{A}\rangle_{\bar{z}}$ $\displaystyle=-\langle
A_{\bar{z}}\rangle=\langle A\rangle_{z},$ $\displaystyle\langle
p\rangle_{\bar{z}}$ $\displaystyle=\langle p\rangle_{z}.$ (60)

In this case the dual stationary state $\bm{|}\bar{z}\bm{)}$ [Eq. (41a)] and
actual one, $\bm{|}z\bm{)}$ [Eq. (42)], _differ_ , but _only_ by inverting the
polarization. This is expected since the duality mapping inverts the energy
spectrum, swapping the Andreev states
$\overline{\bm{|}\tau\bm{)}}=\bm{|}-\tau\bm{)}$ and thus inverting the
polarization of their occupations. In this case $F(\langle
A\rangle_{z},\langle p\rangle_{z})=1$ which by Eq. (59b) is equivalent to
having parity uniquely fixed by polarization in the stationary actual system
and likewise in the dual system:

$\displaystyle\langle p\rangle_{z}$ $\displaystyle=\langle A\rangle_{z}^{2},$
$\displaystyle\langle p\rangle_{\bar{z}}=\langle A\rangle_{\bar{z}}^{2}.$ (61)

One verifies888Inserting Eq. (51) into Eq. (61) gives two solutions:
$\gamma_{c}=\gamma_{p}/2$ and $|\gamma_{s}|=\gamma_{c}$. For the latter case
the physical bound $|\gamma_{s}|\leq\gamma_{p}-\gamma_{c}$ [Eq. (54)] implies
$\gamma_{p}/2\geq\gamma_{c}$. Noting the bounds (63) we see that for $U\geq 0$
the lower-bound $\gamma_{p}/2\leq\gamma_{c}$ implies
$\gamma_{p}/2=\gamma_{c}$, i.e., the second solution is a special case of the
first one, whereas for $U<0$ the strict upper-bound $\gamma_{p}/2>\gamma_{c}$
rules it out.  that this is equivalent to self-duality of the charge decay
rate, $\bar{\gamma}_{c}\equiv\gamma_{p}-\gamma_{c}=\gamma_{c}$, which is
satisfied if and only if

$\displaystyle\gamma_{c}=\tfrac{1}{2}\gamma_{p},$ (62)

or $\gamma_{C}\equiv\gamma_{c}-\gamma_{p}/2=0$. This condition is satisfied at
many parameter points, see Ref. [18], and for $T\to\infty$ it is even
satisfied for all values of the remaining parameters, including $U$. More
interestingly, at any finite $T$ for $U=0$ self-duality always holds as we
verify in App. B.1: By inspection of
$\gamma_{c}=\tfrac{1}{2}\sum_{\eta\tau}W^{\eta}_{1,\tau}$ [Eq. (11),(49a)]
considering all values of $\delta,\mu,\alpha,\Gamma$

$\displaystyle\gamma_{c}$ $\displaystyle>\tfrac{1}{2}\gamma_{p}\Leftrightarrow
U>0,$ $\displaystyle\gamma_{c}$
$\displaystyle=\tfrac{1}{2}\gamma_{p}\Leftrightarrow U=0,$
$\displaystyle\gamma_{c}$ $\displaystyle<\tfrac{1}{2}\gamma_{p}\Leftrightarrow
U<0.$ (63)

Thus, _zero interaction_ is equivalent to requiring _exact_ self-duality (60)
for all values of these remaining parameters. One obtains approximate self-
duality $\gamma_{c}\approx\tfrac{1}{2}\gamma_{p}$ asymptotically by rendering
the interaction ineffective by making one energy scale dominate, e.g., large
bias voltage $|\mu|$, strong pairing $|\alpha|$, large detuning $|\delta|$, or
high temperature $T$, see App. B.2. In Secs. 6.4 and 7 we discuss the impact
of self-duality on the transport in such cases. A system that is self-dual in
the above sense ($U=0$) necessarily has non-negative stationary parity
$\langle p\rangle_{z}\geq 0$ [Eq. (61)] for any value of the remaining
parameters as seen in Fig. 2. Strictly negative parity for some parameters
thus requires nonzero magnitude of the interaction, $U\neq 0$. Such a
violation of self-duality corresponds to a nonzero value of $(+)$ invariant
$\gamma_{C}=\gamma_{c}-\gamma_{p}/2$ [Eq. (63), cf. (55)] with the same sign
as the interaction $U$.

As mentioned, self-duality (60) should not be confused with equality of the
actual and dual stationary state, $\bm{|}\bar{z}\bm{)}=\bm{|}z\bm{)}$. The
latter is a stronger condition and occurs when $\langle
A\rangle_{\bar{z}}=\langle A\rangle_{z}$ and $\langle
p\rangle_{\bar{z}}=\langle p\rangle_{z}$. This is satisfied only999$\langle
A\rangle_{\bar{z}}=\langle A\rangle_{z}$ in Eq. (59) implies $\langle
A\rangle_{\bar{z}}=\langle A\rangle_{z}=0$ since $F\geq 0$ and by Eq. (59b)
$F=[1-\langle p\rangle_{z}])/[1+\langle p\rangle_{z}]$. Requiring $\langle
p\rangle_{\bar{z}}=\langle p\rangle_{z}$ by Eq. (59) gives $[1+\langle
p\rangle_{\bar{z}}]/2=[1+\langle p\rangle_{z}]/2=[1-\langle p\rangle_{z}]/2$
so $\langle p\rangle_{z}=0$.  for $\langle A\rangle_{z}=\langle
p\rangle_{z}=0$ implying $F=1$. It is thus a special case of self-duality
indicated by a black dot in Fig. 2 for which the stationary state is maximally
mixed, $\bm{|}z\bm{)}=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}$, i.e., with uniform
occupation of Andreev states ($z_{\tau}=1/4$) and spin states ($z_{1}=1/2$,
including spin-degeneracy). In addition to $\gamma_{c}=\gamma_{p}/2$ required
by self-duality, this requires vanishing transition-rate asymmetry,
$\gamma_{s}=0$ [Eq. (51a)]. This again can be satisfied at many points101010
In particular, note that $\gamma_{s}=0$ _can_ hold at the often discussed
“symmetry point” $\delta=0$ ($\epsilon=-\tfrac{1}{2}U$) even when $U\neq 0$
_but it need not_ hold, see our numerical analysis in Ref. [18]. , and holds
trivially $T\to\infty$ for all parameters. Unlike self-duality, at $U=0$ it
does not hold for all parameters.

###  5.7 Transport current supercovectors

We can now complete the solution of the transport problem by deriving duality-
adapted expressions for the currents. First, we express the left vector for
the charge current (20a) in terms of the left eigenvectors of $W$. Since we
again from the beginning choose a formulation in terms of scalar invariants
and operators (anti)symmetric under the duality mapping, we find the very
compact form

$\displaystyle\bm{(}I_{N}\bm{|}$
$\displaystyle=\sum_{\eta\tau}\eta\Big{[}W_{1,\tau}^{\eta}\bm{(}\tau\bm{|}+2W_{\tau,1}^{\eta}\bm{(}1\bm{|}\Big{]}$
$\displaystyle=\gamma^{\prime}_{c}\bm{(}\mathds{1}\bm{|}+\gamma^{\prime}_{s}\bm{(}A\bm{|}$
(64a) $\displaystyle=\big{[}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\big{]}\bm{(}z^{\prime}\bm{|}+\gamma^{\prime}_{s}\bm{(}c^{\prime}\bm{|}.$
(64b)

Here we eliminated $W_{\tau,1}^{\eta}$ in favour of $W_{1,\tau}^{\eta}$ [Eq.
(27)], inserted the decomposition (29) of the transport rates, changed basis
using (35) and then rewrote in terms of the left eigenbasis (41). As
anticipated, the primed (= $\eta$-antisymmetric) invariants
$\gamma^{\prime}_{c}$ and $\gamma^{\prime}_{s}$ occur only in the transport
quantities sensitive to electron transfer direction $\eta$.

For the energy and heat current one can proceed analogously111111 Eq. (9) can
be written
$I_{E}(t)=-\sum_{\eta\tau}E_{+,\tau}\big{[}-W_{1,\tau}^{\eta}\rho_{\tau}(t)+W_{\tau,1}^{\eta}\rho_{1}(t)\big{]}=\bm{(}I_{E}\bm{|}\rho(t)\bm{)}$
with left vector
$\bm{(}I_{E}\bm{|}=-\sum_{\eta\tau}E_{+,\tau}\big{[}-W_{1,\tau}^{\eta}\bm{(}\tau\bm{|}+2W_{\tau,1}^{\eta}\bm{(}1\bm{|}\big{]}$
and one then proceeds from there.  by rewriting the left supervector for the
energy current (20b) into the metal. However, it is more convenient to start
from the energy loss of the proximized dot [8],

$\displaystyle I_{E}(t)$ $\displaystyle=-\partial_{t}\langle
H_{\text{D}}\rangle(t)=-\bm{(}H_{\text{D}}\bm{|}W\bm{|}\rho(t)\bm{)}=\bm{(}I_{E}\bm{|}\rho(t)\bm{)},$
(65)

as the Cooper pair condensate does not contribute to the energy of the
proximized dot, whereas it does contribute to its charge. Inserting the
energies (3) and (5) and performing the same steps, reveals the key feature of
the energy, namely, that it couples to interaction $U$ via the parity:

$\displaystyle\bm{(}H_{D}\bm{|}$
$\displaystyle=2E_{1}\bm{(}1\bm{|}+\sum_{\tau}E_{\tau}\bm{(}\tau\bm{|}$
$\displaystyle=\tfrac{1}{2}\delta_{\text{A}}\bm{(}A\bm{|}+\tfrac{1}{4}U\bm{(}p\bm{|}+\big{[}\ldots\big{]}\bm{(}\mathds{1}\bm{|}$
(66a) $\displaystyle=\tfrac{1}{2}\big{(}\delta_{\text{A}}-U\langle
A\rangle_{\bar{z}}\big{)}\bm{(}c^{\prime}\bm{|}+U\bm{(}p^{\prime}\bm{|}+\big{[}\ldots\big{]}\bm{(}z^{\prime}\bm{|}.$
(66b)

The terms with coefficients $[\ldots]$ are diagonal in the energy basis but
don’t contribute to the energy current since $\bm{(}z^{\prime}\bm{|}W=0$.
Thus, $\bm{(}I_{E}\bm{|}=-\bm{(}H_{D}\bm{|}W$ is given by

$\displaystyle\bm{(}I_{E}\bm{|}=\tfrac{1}{2}\big{(}\delta_{\text{A}}-U\langle
A\rangle_{\bar{z}}\big{)}\gamma_{c}\bm{(}c^{\prime}\bm{|}+U\gamma_{p}\bm{(}p^{\prime}\bm{|}.$
(67)

The heat current $I_{Q}(t)=I_{E}(t)-\mu
I_{N}(t)=\bm{(}I_{Q}\bm{|}\rho(t)\bm{)}$ is obtained from

$\displaystyle\bm{(}I_{Q}\bm{|}=$
$\displaystyle-\mu\Big{[}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\Big{]}\bm{(}z^{\prime}\bm{|}+\Big{[}\tfrac{1}{2}\big{(}\delta_{\text{A}}-U\langle
A\rangle_{\bar{z}}\big{)}\gamma_{c}-\mu\gamma^{\prime}_{s}\Big{]}\bm{(}c^{\prime}\bm{|}+U\gamma_{p}\bm{(}p^{\prime}\bm{|}.$
(68)

##  6 General analytical result

The main results of the paper are the final analytical formula for the state
evolution

$\displaystyle\bm{|}\rho(t)\bm{)}=$
$\displaystyle\big{\\{}\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
A\rangle_{z}\tfrac{1}{2}\bm{|}A\bm{)}+\langle
p\rangle_{z}\tfrac{1}{4}\bm{|}p\bm{)}\big{\\}}+\tfrac{1}{2}\Big{[}\bm{|}A\bm{)}-\langle
A\rangle_{\bar{z}}\bm{|}p\bm{)}\Big{]}e^{-\gamma_{c}t}\big{[}\langle
A\rangle_{\rho_{0}}-\langle A\rangle_{z}\big{]}$
$\displaystyle+\bm{|}p\bm{)}e^{-\gamma_{p}t}\Big{\\{}\tfrac{1}{4}\big{[}\langle
p\rangle_{\rho_{0}}-\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
A\rangle_{\bar{z}}\big{[}\langle A\rangle_{\rho_{0}}-\langle
A\rangle_{z}\big{]}\Big{\\}},$ (69)

and the formulas for the charge and heat current,

$\displaystyle I_{N}(t)=$
$\displaystyle\big{(}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\big{)}+\gamma^{\prime}_{s}e^{-\gamma_{c}t}\big{[}\langle
A\rangle_{\rho_{0}}-\langle A\rangle_{z}\big{]},$ (70) $\displaystyle
I_{Q}(t)=$
$\displaystyle-\mu\big{\\{}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\big{\\}}+\Big{\\{}\tfrac{1}{2}\big{(}\delta_{\text{A}}-U\langle
A\rangle_{\bar{z}}\big{)}\gamma_{c}-\mu\gamma^{\prime}_{s}\Big{\\}}e^{-\gamma_{c}t}\big{[}\langle
A\rangle_{\rho_{0}}-\langle A\rangle_{z}\big{]}$
$\displaystyle+U\gamma_{p}e^{-\gamma_{p}t}\Big{\\{}\tfrac{1}{4}\big{[}\langle
p\rangle_{\rho_{0}}-\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
A\rangle_{\bar{z}}\big{[}\langle A\rangle_{\rho_{0}}-\langle
A\rangle_{z}\big{]}\Big{\\}},$ (71)

where the energy current $I_{E}(t)$ is given by setting $\mu=0$ in the
expression for $I_{Q}(t)$. This solution applies to any initial state
$\bm{|}\rho_{0}\bm{)}$ that is diagonal in the energy basis, specified by the
initial expectation values of the polarization $\langle
A_{0}\rangle_{\rho_{0}}$ and parity $\langle p\rangle_{\rho_{0}}$ subject to
the constraint (54), $|\langle p\rangle_{\rho_{0}}|\leq 1$ and $|\langle
A_{0}\rangle_{\rho_{0}}|\leq\tfrac{1}{2}[1+\langle p\rangle_{\rho_{0}}]$,
guaranteeing that the initial state is physical. where by $A_{0}$ we indicate
the dependence of polarization on initial parameters. This is the solution to
the time-evolution and transport problem of Ref. [35]. In the following we
mention its salient features, discussing some specific limiting cases
analytically in Sec. 7, thereby preparing for the numerical analysis of the
general case in Ref. [18].

###  6.1 State evolution

The state supervector (69) is obtained by inserting the amplitude functions
(53) into Eq. (52). The key advantage of this duality-based solution is that
we obtain the state-evolution _entirely_ in terms of _stationary_ expectation
values of the appropriate observables $p$ and $A$, in particular the all-
important amplitudes of the transient decay with rates $\gamma_{c}$ and
$\gamma_{p}$. Roughly speaking, if we know the initial and stationary values
of these observables we understand the entire intermediate state evolution.
Here, the amplitudes are expressed using result (53d) as the difference
between initial and stationary value of the observables $A$ and $p$, their
“initial excess”. This makes explicit that, as expected,
$\bm{|}\rho(t)\bm{)}=\bm{|}\rho_{0}\bm{)}$ for all $t\geq 0$ if
$\bm{|}\rho_{0}\bm{)}=\bm{|}z\bm{)}$ unlike the form (53c) used in Ref. [8],
where the “excess” of observables with respect to their stationary value does
not explicitly enter see Sec. 7.2.

As in prior work we find that the relevant stationary values are _not_ only
those of the actual system but also those of the _dual_ system. This
remarkable insight has significantly advanced the analysis of the parameter
dependence of the transient state evolution [8, 13, 18, 52]. Our result
extends this to interacting systems proximized by a superconductor but
additionally points out two important refinements121212Equation (69) is thus
more than a rewriting of the the expression for $\bm{|}\rho(t)\bm{)}$ given in
Ref. [37], see Eq. (20-25) of that work for the special case of zero
transverse Bloch vector ($I_{x}=I_{y}=0$). . It is well known that a time-
evolving quantum state can be expanded in a set of fixed observables
(generalized Bloch expansion) with coefficients which are their _time-
dependent_ averages in that state. Here we instead express the time-dependent
state in averages of the _stationary_ state ($t\to\infty$), the initial state,
and decay exponentials.

First, the relation
$\bm{(}p^{\prime}\bm{|}\rho_{0}\bm{)}=\bm{(}p\bar{z}\bm{|}\rho_{0}\bm{)}$ [Eq.
(40)] indicates that the amplitude of the $\gamma_{p}$-decay depends on how
well the initial state $\bm{|}\rho_{0}\bm{)}$ “compares” with the _dual_
stationary state $\bm{|}\bar{z}\bm{)}$ as discussed in Ref. [8, 15]. The above
mentioned form (53c) of this amplitude suggests that one needs to analyse the
dual stationary parity $\langle p\rangle_{\bar{z}}$ for this. However, our
result (53d) explicitly shows that the dual stationary state enters _only_
through its polarization $\langle A\rangle_{\bar{z}}$ and _only_ provided that
there is an initial polarization excess, $\langle
A\rangle_{\rho_{0}}\neq\langle A\rangle_{z}$. This simplifies the analysis and
these results carry over to the case $\alpha=0$, see Sec. 7.1.

Second, we can go one step further and substitute the stationary duality
relation (59a) into Eq. (69): we see that the amplitudes of the transient
decay with rates $\gamma_{c}$ and $\gamma_{p}$ ultimately depend only on the
_stationary_ expectation values of $A$ and $p$ of the _system alone_. Thus,
apart from the decay rates the _evolution towards_ the stationary state can be
understood in detail by the stationary state _itself_. If anything, this _non-
linear_ dependence on its stationary state (59b) highlights the non-trivial
nature of the transient evolution. It does not invalidate the usefulness of
above mentioned analysis via dual stationary values which enter linearly: One
either analyses the actual stationary system and the dual parity or one
analyses only the actual stationary system but uses the non-linear –but
parameter independent– relation (59).

###  6.2 Transport evolution

The transport currents (70)-(71) follow by taking the scalar product of the
state (69) with the corresponding left current supervectors [Eqs. (64b) and
(68)] and using bi-orthonormality of the eigen-supervectors. This leads to the
physically appealing structure

$\displaystyle I_{x}(t)$
$\displaystyle=\bm{(}I_{x}\bm{|}\rho(t)\bm{)}=\bm{(}I_{x}\bm{|}z\bm{)}+\bm{(}I_{x}\bm{|}c\bm{)}e^{-\gamma_{c}t}\bm{(}c^{\prime}\bm{|}\rho_{0}\bm{)}+\bm{(}I_{x}\bm{|}p\bm{)}e^{-\gamma_{p}t}\bm{(}p^{\prime}\bm{|}\rho_{0}\bm{)}$
(72)

for $x=N,E,Q$. Denoting $\gamma_{z}\equiv 0$, each term $y=z,c,p$ contributing
to a transport quantity $I_{x}(t)$ on time-scale $\gamma_{y}^{-1}$ depends on
two factors $\bm{(}I_{x}\bm{|}y\bm{)}$ and
$\bm{(}y^{\prime}\bm{|}\rho_{0}\bm{)}$. The first factor quantifies the
ability of quantity $x$ to “probe” the $y$-part of the state evolution, the
second one captures the initial-state dependence of the latter, see also Ref.
[52]. For example, the charge current (70) does not probe the parity decay,
$\bm{(}I_{N}\bm{|}p\bm{)}=0$, and thus –by duality– depends on the physical
parameters only through the stationary state of the actual system, in
particular, through the value of the polarization observable $\langle
A\rangle_{z}$ in $\bm{(}c^{\prime}\bm{|}\rho_{0}\bm{)}$ [Eq. (53b)]. By
contrast, the energy and heat current do probe the parity decay,
$\bm{(}I_{E}\bm{|}p\bm{)},\bm{(}I_{Q}\bm{|}p\bm{)}\neq 0$, and thus both
_linearly_ depend on the dual stationary state through $\langle
A\rangle_{\bar{z}}$. (As discussed at the end of the previous section, this
dependence can ultimately be brought back to $\langle A\rangle_{z}$ and
$\langle p\rangle_{z}$ relation (59) but only at the price of _non-linear_
dependence.)

As for the state evolution, _non-trivial_ time-dependence of the transport
–which is supposed to probe the state of the system– is expressed explicitly
in terms of the stationary state, i.e., the values of observables $p$ and $A$
that characterize it. In the case of energy transport (71) the additional pre-
factors $\delta$, $U$ and $\mu$ that enter do not complicate the analysis. For
the charge transport (70) the charge-sensitive coefficients
$\gamma^{\prime}_{s}$ and $\gamma^{\prime}_{c}$ do depend non-trivially on the
parameters, but their dependence is again restricted by duality: As we show
below in Sec. 7, their parameter dependence is closely related to that of
$\gamma_{c}$ and $\gamma_{s}$, respectively.

Finally, we note that the explicit dependence of the heat current (71) on the
excess polarization and parity can be used to separately measure its different
contributions provided one has full control over the initial state. For
example, for an initial state with no excess polarization, $\langle
A\rangle_{\rho_{0}}=\langle A\rangle_{z}$, only the second term in Eq. (71)
contributes with single-exponential decay with rate $\gamma_{p}$. By contrast,
for an initial state with no excess parity, $\langle
p\rangle_{\rho_{0}}=\langle p\rangle_{z}$, the decay is still double
exponential with rates $\gamma_{c}$ and $\gamma_{p}$.

###  6.3 Components of duality invariants

We now illustrate how the fermionic duality can be used to investigate the
solutions (69)-(71) in more detail by decomposing the (shifted-) duality
invariants into “components”

$\displaystyle\gamma_{C}$
$\displaystyle=\kappa_{C}+\frac{\delta}{\delta_{\text{A}}}\kappa^{\prime}_{s},$
$\displaystyle\gamma_{c}$
$\displaystyle=\kappa_{c}+\frac{\delta}{\delta_{\text{A}}}\kappa^{\prime}_{s},$
$\displaystyle\gamma^{\prime}_{c}$
$\displaystyle=\kappa^{\prime}_{c}+\frac{\delta}{\delta_{\text{A}}}\kappa_{s},$
(73a) $\displaystyle\gamma_{s}$
$\displaystyle=\kappa_{s}+\frac{\delta}{\delta_{\text{A}}}\kappa^{\prime}_{c},$
$\displaystyle\gamma^{\prime}_{s}$
$\displaystyle=\kappa^{\prime}_{s}+\frac{\delta}{\delta_{\text{A}}}\kappa_{c},$
$\displaystyle\gamma^{\prime}_{S}$
$\displaystyle=\kappa^{\prime}_{s}+\frac{\delta}{\delta_{\text{A}}}\kappa_{C},$
(73b)

while still avoiding the more cumbersome full expressions of the rates. The
decomposition (73) separates the contributions from the two summands in
$\Gamma_{\pm}=\gamma_{p}(1\pm\delta/\delta_{\text{A}})/2$ [Eq. (13)] which
enter the invariants via Eq. (11): Formally each component is obtained by
replacing $\Gamma_{\eta\tau}\to\tfrac{1}{2}\gamma_{p}$ in the rate-expression
of the corresponding invariant (49) after inserting Eq. (11). Explicitly,

$\displaystyle\kappa_{c}$
$\displaystyle=\tfrac{1}{4}\gamma_{p}\sum_{\eta\tau}f^{-\eta}(E_{\eta,\tau}-\mu),$
$\displaystyle\kappa^{\prime}_{c}$
$\displaystyle=\tfrac{1}{4}\gamma_{p}\sum_{\eta\tau}\eta
f^{-\eta}(E_{\eta,\tau}-\mu),$ (74a) $\displaystyle\kappa_{s}$
$\displaystyle=\tfrac{1}{4}\gamma_{p}\sum_{\eta\tau}\tau
f^{-\eta}(E_{\eta,\tau}-\mu),$ $\displaystyle\kappa^{\prime}_{s}$
$\displaystyle=\tfrac{1}{4}\gamma_{p}\sum_{\eta\tau}\tau\eta
f^{-\eta}(E_{\eta,\tau}-\mu).$ (74b)

That these are appropriate variables is established by the fact that they
transform in the same simple way as the corresponding invariants [Eqs. (30),
(46)]:

$\displaystyle\kappa_{C}$ $\displaystyle=-\bar{\kappa}_{C},$
$\displaystyle\kappa_{c}$ $\displaystyle=\gamma_{p}-\bar{\kappa}_{c},$
$\displaystyle\kappa_{s}$ $\displaystyle=\bar{\kappa}_{s},$ (75a)
$\displaystyle\kappa^{\prime}_{c}$
$\displaystyle=\phantom{\kappa_{p}}-\bar{\kappa}^{\prime}_{c},$
$\displaystyle\kappa^{\prime}_{s}$ $\displaystyle=\bar{\kappa}^{\prime}_{s}.$
(75b)

Like the four (shifted-)invariants that they decompose, the behaviour of these
four components as function of the physical parameters provides a complete
understanding of the time-dependence of the solution of the problem. The prime
notation for the components makes the same physical distinction as for the
invariants [a component has (no) prime if it determines charge transfer (state
evolution)] which is, however, reversed when occurring with prefactor
$\delta/\delta_{\text{A}}$. For example, primed invariants
$\gamma_{c}^{\prime},\gamma_{s}^{\prime}$ [Eq. (64a)] associated with
transport [Eq. (64a)] count charge transfer $\eta$ with effective rates
$\Gamma_{\eta\tau}=\gamma_{p}(1+\eta\tau\delta/\delta_{\text{A}})/2$: a
transport component is thus either $\eta$-antisymmetric or $\eta$-symmetric
with factor $\delta/\delta_{\text{A}}$.

The components thus have the advantage that they cleanly separate all resonant
behaviour of the superconductor (explicit $\delta$ dependence relative to
$\mu_{\text{S}}=0$ through $\Gamma_{\eta\tau}$) from the superconductor’s
indirect effect via the Andreev levels (implicit $\delta$ dependence relative
to $\mu$ through $f^{-\eta}(E_{\eta,\tau})$ in Eq. (74)). By doing so, the
decomposition (73) shows that the pairs of duality invariants $\gamma_{c}$,
$\gamma^{\prime}_{s}$ and $\gamma^{\prime}_{c}$, $\gamma_{s}$, respectively,
playing distinct physical roles, actually have similar behaviour since they
depend on the same components which can mix up and cancel out when simplifying
final expressions at the price of losing the physical distinction. This is
nicely illustrated for the expression for stationary charge current
$I_{N}(\infty)=\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle A\rangle_{z}$
given by the first term in Eq. (70):

$\displaystyle
I_{N}(\infty)=\frac{\gamma^{\prime}_{c}\gamma_{c}-\gamma^{\prime}_{s}\gamma_{s}}{\gamma_{c}}=\dfrac{\kappa_{c}\kappa^{\prime}_{c}-\kappa_{s}\kappa_{s}^{\prime}}{\kappa_{c}+\frac{\delta}{\delta_{A}}\kappa^{\prime}_{s}}\,\Big{(}1-\frac{\delta^{2}}{\delta_{A}^{2}}\Big{)}.$
(76)

We automatically extract a resonant Lorentzian $\delta$-dependence reflecting
the non-equilibrium Cooper pair transport to the superconductor while its non-
trivial modification by the Andreev states is contained in the pre-factor
governed by the components of the invariants. The latter is a joint effect of
both transport and state-evolution invariants. This stationary current was
studied in Refs. [33, 34, 35, 22] in the weak-coupling limit as well as in
Refs. [53, 54, 55, 56] in the limit of strong coupling where analytical
results were limited to specific regimes. Note that the stationary energy
current equals zero, $I_{E}(\infty)=0$.

###  6.4 Self-duality and transient transport

In the limiting cases discussed in Sec. 7 the state and transport evolutions
simplify substantially when the system obeys self-duality (60), i.e.,
$\gamma_{c}=\tfrac{1}{2}\gamma_{p}$ [Eq. (62)]. As discussed in Sec. 5.6, this
can either hold exactly for all applied voltages ($U=0$) or asymptotically
(e.g., $\mu\to\infty$). We can then generically simplify the stationary state
(41a) by eliminating parity in favour of polarization [Eq. (61)]:

$\displaystyle\bm{|}z\bm{)}=\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
A\rangle_{z}\tfrac{1}{2}\bm{|}A\bm{)}+\langle
A\rangle_{z}^{2}\tfrac{1}{4}\bm{|}p\bm{)}=\sum_{\tau}\bm{|}\tau\bm{)}\tfrac{1}{4}\big{[}1+\tau\langle
A\rangle_{z}\big{]}^{2}+\bm{|}1\bm{)}\tfrac{1}{2}\big{[}1-\langle
A\rangle_{z}^{2}\big{]}.$ (77)

In the evolving state (69) this is not possible without further consideration
since self-duality need not hold initially, $\langle
p\rangle_{\rho_{0}}\neq\langle A_{0}\rangle_{\rho_{0}}^{2}$, because the
energy mixture $\rho_{0}$ is arbitrary. However, as discussed in Ref. [18], a
relevant class of initial states are those prepared from a stationary state of
the system at some different gate voltage $\delta_{0}$ –which by the above
assumptions is still self-dual– and thus $\langle p\rangle_{z_{0}}=\langle
A_{0}\rangle_{z_{0}}^{2}$. This way one obtains an initial state with $\langle
p\rangle_{\rho_{0}}=\langle p\rangle_{z_{0}}$ and $\langle
A\rangle_{\rho_{0}}=\theta\langle A_{0}\rangle_{z_{0}}$ where
$\theta\in[-1,1]$ captures details of the initialization [18]. In this case
all decay amplitudes in the evolving state are fixed by the stationary and
initial polarization:

$\displaystyle\bm{|}\rho(t)\bm{)}=$
$\displaystyle\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
A\rangle_{z}\tfrac{1}{2}\bm{|}A\bm{)}+\langle
A\rangle_{z}^{2}\tfrac{1}{4}\bm{|}p\bm{)}+\tfrac{1}{2}\big{[}\bm{|}A\bm{)}+\langle
A\rangle_{z}\bm{|}p\bm{)}\big{]}e^{-\tfrac{1}{2}\gamma_{p}t}\big{[}\theta\langle
A_{0}\rangle_{z_{0}}-\langle A\rangle_{z}\big{]}$
$\displaystyle+\bm{|}p\bm{)}e^{-\gamma_{p}t}\Big{\\{}\tfrac{1}{4}\langle
A_{0}\rangle_{z_{0}}^{2}+\tfrac{1}{4}\langle
A\rangle_{z}^{2}-\tfrac{1}{2}\theta\langle A\rangle_{z}\langle
A_{0}\rangle_{z_{0}}\Big{\\}}.$ (78)

Moreover, the coefficient of the parity contribution can be verified to be
non-negative which in general [Eq. (53c)-53d] need not be true even when
accounting for constraints (57)-(58). The currents likewise simplify,

$\displaystyle I_{N}(t)=$
$\displaystyle\big{(}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\big{)}+\gamma^{\prime}_{s}e^{-\tfrac{1}{2}\gamma_{p}t}\big{[}\theta\langle
A_{0}\rangle_{z_{0}}-\langle A\rangle_{z}\big{]}$ (79) $\displaystyle
I_{Q}(t)=$
$\displaystyle-\mu\big{\\{}\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}\big{\\}}+\big{\\{}\tfrac{1}{2}\delta_{\text{A}}\gamma_{c}-\mu\gamma^{\prime}_{s}\big{\\}}e^{-\tfrac{1}{2}\gamma_{p}t}\big{[}\theta\langle
A_{0}\rangle_{z_{0}}-\langle A\rangle_{z}\big{]}.$ (80)

In particular $I_{Q}(t)$ has no parity contribution, either exactly (for
$U=0$) or asymptotically (for large $\mu$ where it is dominated by the terms
$\propto\mu$, in which case $\tfrac{1}{2}\delta_{\text{A}}\gamma_{c}$ must
also be dropped). Finally, whenever we have stationary transport, $\alpha\neq
0$, self-duality $\gamma_{c}=\gamma_{p}/2$, is equivalent to having a constant
component, $\kappa_{c}=\gamma_{p}/2$, and a vanishing one,
$\kappa^{\prime}_{s}=0$ [Eq. (49), (73)]. The stationary current (76) is then
modulated by just one component:

$\displaystyle I_{N}(\infty)=\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}=\kappa^{\prime}_{c}\times\Big{(}1-\frac{\delta^{2}}{\delta_{A}^{2}}\Big{)}.$
(81)

##  7 Limiting cases

Finally, to highlight various physical effects captured by the duality
invariants we present explicit formulas for three different limiting cases,
taking either $\alpha\to 0$, $U\to 0$ or $\mu\to\infty$.

###  7.1 Interacting dot without pairing ($\alpha=0$)

For a quantum dot coupled only to a normal metal ($\alpha=0$) the fermionic
duality expresses a symmetry of the problem which is not at all obvious in the
presence of strong interaction ($U\neq 0$) [8, 14]. It is instructive to see
how this case is recovered noting that by the assumption (6) we _cannot_
simply send $\alpha\to 0$ within the range of applicability of the master
equation (7) derived in Ref. [35]. However, for $\Gamma\ll|\alpha|\ll T$
temperature fluctuations eradicate every effect of the pairing except for a
gate voltage regime around the resonance with the superconductor,
$|\delta|\lesssim|\alpha|$, which becomes vanishingly small as $\alpha\to 0$.

Thus, in this limit our general solution should coincide with that obtained
for $\alpha=0$ in Ref. [8] at every $\delta\neq 0$. Indeed, when closing the
pairing-induced gap $\alpha\to 0$ the transport- and transition rates become
discontinuous functions of the gate voltage $\delta$ through both the
transition energies (12) and the effective rates (13). The origin of the
discontinuity lies in Eq. (3) where for $\alpha\neq 0$ we label the
eigenenergies $H_{\text{D}}$ such that they correspond to $\delta$-continuous
branches of the square root. Correspondingly, the labelling of eigenstates (4)
becomes a discontinuous choice of labelling the obvious eigenstates at
$\alpha=0$:

$\displaystyle\ket{\tau}\bra{\tau}\overset{\alpha\to
0}{\quad=\quad}\ket{0}\bra{0}\text{ for $\tau=-\frac{\delta}{|\delta|}$}\quad
and\quad\ket{\tau}\bra{\tau}\overset{\alpha\to
0}{\quad=\quad}\ket{2}\bra{2}\text{ for $\tau=+\frac{\delta}{|\delta|}$}.$
(82)

This gate-voltage dependence of the states is essential to the problem with
superconducting pairing $\alpha\neq 0$131313 This gate-voltage dependence
warrants the careful discussion of the state initialization in Ref. [18].  and
there is no reason to proceed otherwise. However, for $\alpha\to 0$ no
physical result at $\delta\neq 0$ should depend on the choice of labelling
made for $\alpha\neq 0$ but verifying this is messy since _all_ transport
rates (11) become discontinuous for $\alpha\to 0$. The duality invariants now
have the nice property that for $\alpha\to 0$ they are either continuous or
discontinuous only by the _sign factor_ $\delta/|\delta|$. Duality thus
automatically either cancels all discontinuities or collects them into a
single factor. This is clearly useful for numerical calculations but also
simplifies the analysis of the limiting behaviour of proximized quantum dots
here.

Indeed, for $\alpha\to 0$ the invariants $\gamma_{s}$, $\gamma^{\prime}_{s}$
associated with an asymmetry with respect to energy eigenstates [$\tau$
branches, see Eq. (49a)-(49b)] are discontinuous and are related by the
$\delta$-sign to continuous invariants $\gamma_{c}$, $\gamma^{\prime}_{c}$
which are ignorant of this asymmetry:

$\displaystyle\gamma_{s}\overset{\alpha\to
0}{\quad=\quad}\frac{\delta}{|\delta|}\gamma^{\prime}_{c},\qquad\gamma^{\prime}_{s}\overset{\alpha\to
0}{\quad=\quad}\frac{\delta}{|\delta|}\gamma_{c}.$ (83)

Indeed, in Eqs. (49a)-(49b) only the term $\tau=\eta\delta/|\delta|$
contributes to invariants $\gamma_{c}$, $\gamma^{\prime}_{c}$, giving
continuous functions of the transition energies
$E_{\eta,\eta}=\tfrac{1}{2}(\delta+\eta U)=\epsilon,\epsilon+U$ for
$\eta=\mp$: [Eq. (12)]:

$\displaystyle\gamma_{c}$
$\displaystyle=\tfrac{1}{2}\gamma_{p}\sum_{\eta}f^{-\eta}(E_{\eta,\eta}-\mu),$
$\displaystyle\quad\gamma^{\prime}_{c}$
$\displaystyle=\tfrac{1}{2}\gamma_{p}\sum_{\eta}\eta
f^{-\eta}(E_{\eta,\eta}-\mu).$ (84)

Similarly, the parameter-dependent polarization observable [cf. Eq. (34)] is
discontinuous in $\delta$ but in a simple way, see Eq. (82): with the charge
polarization $N-\mathds{1}=-\ket{0}\bra{0}+\ket{2}\bra{2}$

$\displaystyle A=\sum_{\tau}\tau\ket{\tau}\bra{\tau}\overset{\alpha\to
0}{\,=\,}\frac{\delta}{|\delta|}(N-\mathds{1}).$ (85)

Inserting Eq. (83) and Eq. (85) into the observables (43) we recover the
continuous expressions of Ref. [8] but in the more compact form of duality
invariants,

$\displaystyle\langle N-\mathds{1}\rangle_{z}$
$\displaystyle\overset{\alpha\to
0}{\quad=\quad}-\frac{\gamma^{\prime}_{c}}{\gamma_{c}},$
$\displaystyle\quad\langle p\rangle_{z}$ $\displaystyle\overset{\alpha\to
0}{\quad=\quad}1-2\frac{\gamma_{c}^{2}-\gamma_{c}^{\prime
2}}{\gamma_{p}\gamma_{c}},$ (86a) $\displaystyle\langle
N-\mathds{1}\rangle_{\bar{z}}$ $\displaystyle\overset{\alpha\to
0}{\quad=\quad}\frac{\gamma^{\prime}_{c}}{\gamma_{p}-\gamma_{c}},$
$\displaystyle\langle p\rangle_{\bar{z}}$ $\displaystyle\overset{\alpha\to
0}{\quad=\quad}1-2\frac{(\gamma_{p}-\gamma_{c})^{2}-\gamma_{c}^{\prime
2}}{\gamma_{p}(\gamma_{p}-\gamma_{c})}.$ (86b)

Thus, for “$\alpha\to 0$” our appropriate duality-adapted observables
essentially reduce to the parity $p$ and the charge polarization
$N-\mathds{1}$ and our four duality invariants simplify to just two invariants
$\gamma_{c}$ (state) and $\gamma^{\prime}_{c}$ (transport) both continuous in
gate voltage $\delta$. This result shows that also in Ref. [8] it is possible
to eliminate all dual stationary observables in favour of the actual ones.
Indeed, the stationary duality relation (59) also holds with the substitution
$A\to N-\mathds{1}$ (the $\delta$-signs cancel) and shows the consistency of
the actual and dual positivity constraints $|\langle
N-\mathds{1}\rangle_{z}|\leq\tfrac{1}{2}[1+\langle p\rangle_{z}]$ and
$|\langle N-\mathds{1}\rangle_{\bar{z}}|\leq\tfrac{1}{2}[1+\langle
p\rangle_{\bar{z}}]$.

In the state evolution (69) the $\delta$-signs cancel between duality-adapted
observables and their expectation values recovering the continuous Eqs.
(S3-S5) of Ref. [8] for $\alpha\to 0$:

$\displaystyle\bm{|}\rho(t)\bm{)}\overset{\alpha\to 0}{\quad=\quad}$
$\displaystyle\Big{[}\tfrac{1}{4}\bm{|}\mathds{1}\bm{)}+\langle
N-\mathds{1}\rangle_{z}\tfrac{1}{2}\bm{|}N-\mathds{1}\bm{)}+\langle
p\rangle_{z}\tfrac{1}{4}\bm{|}p\bm{)}\Big{]}$
$\displaystyle+\tfrac{1}{2}\Big{[}\bm{|}N-\mathds{1}\bm{)}-\langle
N-\mathds{1}\rangle_{\bar{z}}\bm{|}p\bm{)}\Big{]}e^{-\gamma_{c}t}\big{[}\langle
N\rangle_{\rho_{0}}-\langle N\rangle_{z}\big{]}$
$\displaystyle+\bm{|}p\bm{)}e^{-\gamma_{p}t}\Big{\\{}\tfrac{1}{4}\big{[}\langle
p\rangle_{\rho_{0}}-\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
N-\mathds{1}\rangle_{\bar{z}}\big{[}\langle N\rangle_{\rho_{0}}-\langle
N\rangle_{z}\big{]}\Big{\\}}.$ (87)

Likewise, for the currents (70)-(71) we recover Eq. (S8) of Ref. [8] with
continuous factors:

$\displaystyle I_{N}(t)\overset{\alpha\to 0}{\quad=\quad}$
$\displaystyle\gamma_{c}e^{-\gamma_{c}t}\big{[}\langle
N\rangle_{\rho_{0}}-\langle N\rangle_{z}\big{]}$ (88) $\displaystyle
I_{Q}(t)\overset{\alpha\to 0}{\quad=\quad}$
$\displaystyle\big{[}\tfrac{1}{2}\big{(}\delta-U\langle
N-\mathds{1}\rangle_{\bar{z}}\big{)}-\mu\big{]}\gamma_{c}e^{-\gamma_{c}t}\big{[}\langle
N\rangle_{\rho_{0}}-\langle N\rangle_{z}\big{]}$
$\displaystyle+U\gamma_{p}e^{-\gamma_{p}t}\Big{\\{}\tfrac{1}{4}\big{[}\langle
p\rangle_{\rho_{0}}-\langle p\rangle_{z}\big{]}+\tfrac{1}{2}\langle
N-\mathds{1}\rangle_{\bar{z}}\big{[}\langle N\rangle_{\rho_{0}}-\langle
N\rangle_{z}\big{]}\Big{\\}}.$ (89)

Compared to the result of Ref. [8] the last term in (87) and (89) has been
simplified using Eq. (44) and (53d). This makes explicit that without initial
excess charge ($\langle N\rangle_{\rho_{0}}=\langle N\rangle_{z}$) the heat
current exhibits strictly single-exponential $\gamma_{p}$ decay, whereas
without initial excess parity ($\langle p\rangle_{\rho_{0}}=\langle
p\rangle_{z}$) one still has double exponential heat decay, cf. Sec. 6.2.
These expressions have been successfully applied to the analysis of both
repulsive ($U>0$) and attractive ($U<0$) quantum dots [14].

A final point to note is that the distinct behaviour $\langle N\rangle_{z}$
and $\langle p\rangle_{z}$ is due to the interaction: only when additionally
sending $U\to 0$ do we have $\gamma_{c}\to\tfrac{1}{2}\gamma_{p}$ since
$E_{\eta,\eta}=\epsilon$ independent of $\eta$. By Eq. (62) the system is then
self-dual such that the considerations of Sec. 6.4 apply: the parity is
completely fixed by the charge-polarization

$\displaystyle-\langle N-\mathds{1}\rangle_{\bar{z}}$
$\displaystyle\overset{\alpha,U\to 0}{\quad=\quad}\langle
N-\mathds{1}\rangle_{z},$ $\displaystyle\langle
p\rangle_{z}\overset{U,\alpha\to 0}{\quad=\quad}$ $\displaystyle\langle
N-\mathds{1}\rangle_{z}^{2},$ $\displaystyle\langle
p\rangle_{\bar{z}}\overset{U,\alpha\to 0}{\quad=\quad}$ $\displaystyle\langle
N-\mathds{1}\rangle_{\bar{z}}^{2}.$ (90)

The other nonzero invariant simplifies to an antisymmetric step function,
$\gamma^{\prime}_{c}\to\tfrac{1}{2}\gamma_{p}\big{[}1-2f(\epsilon-\mu)\big{]}$
which determines the charge-polarization [Eq. (86a)]:

$\displaystyle\langle N-\mathds{1}\rangle_{z}\overset{\alpha,U\to
0}{\quad=\quad}2f(\epsilon-\mu)-1.$ (91)

###  7.2 Proximized dot without interaction ($U=0$)

For a non-interacting dot ($U=0$) proximized by a superconductor ($\alpha\neq
0$), see Ref. [57], fermionic duality also remains a dissipative symmetry but
its implications in this case have not been considered. Since in this case
$\gamma_{c}$ has one vanishing component, $\kappa^{\prime}_{s}\to 0$, and the
other is constant, $\kappa_{c}\to\gamma_{p}/2$, the invariants again simplify:

$\displaystyle\gamma_{c}\overset{U\to
0}{\quad=\quad}\tfrac{1}{2}\gamma_{p},\qquad\gamma^{\prime}_{s}\overset{U\to
0}{\quad=\quad}\tfrac{1}{2}\gamma_{p}\frac{\delta}{\delta_{\text{A}}}.$ (92)

and we have _exact_ self-duality by Eq. (62) and Sec. 6.4 again applies. Here
this happens since $E_{\eta,\tau}\to\tfrac{1}{2}\eta\tau\delta_{\text{A}}$
depends only on the product $\eta\tau\equiv\lambda=\pm$ unlike the interacting
case [Eq. (12)]. For the two remaining invariants
$\gamma^{\prime}_{c}=\kappa^{\prime}_{c}+\frac{\delta}{\delta_{\text{A}}}\kappa_{s}$
and
$\gamma_{s}=\kappa_{s}+\frac{\delta}{\delta_{\text{A}}}\kappa^{\prime}_{c}$
the components read

$\displaystyle\kappa_{s}$ $\displaystyle\overset{U\to
0}{\quad=\quad}\tfrac{1}{2}\gamma_{p}\sum_{\lambda\eta}\tfrac{1}{2}\lambda\eta
f^{-\eta}(\tfrac{1}{2}\lambda\delta_{\text{A}}-\mu),$
$\displaystyle\kappa^{\prime}_{c}$ $\displaystyle\overset{U\to
0}{\quad=\quad}\tfrac{1}{2}\gamma_{p}\sum_{\lambda\eta}\tfrac{1}{2}\eta
f^{-\eta}(\tfrac{1}{2}\lambda\delta_{\text{A}}-\mu).$ (93)

Importantly, for $U=0$ the self-duality (92) holds for all values of the
remaining parameters $\delta,\mu,T,\Gamma$, and in particular, the pairing
$\alpha$, such that we have according to Eq. (51)

$\displaystyle\langle A\rangle_{z}$ $\displaystyle\overset{U\to
0}{\quad=\quad}\langle\bar{A}\rangle_{\bar{z}},$ $\displaystyle\langle
p\rangle_{z}$ $\displaystyle\overset{U\to 0}{\quad=\quad}\langle
p\rangle_{\bar{z}},$ $\displaystyle\langle p\rangle_{z}$
$\displaystyle\overset{U\to 0}{\quad=\quad}\langle A\rangle_{z}^{2},$
$\displaystyle\langle p\rangle_{\bar{z}}$ $\displaystyle\overset{U\to
0}{\quad=\quad}\langle A\rangle_{\bar{z}}^{2}.$ (94)

Applying all simplifications (77)-(81) we obtain for the stationary
polarization

$\displaystyle\langle A\rangle_{z}$ $\displaystyle\overset{U\to
0}{\quad=\quad}-2\frac{\gamma_{s}}{\gamma_{p}}=-\sum_{\lambda}\tfrac{1}{2}\Big{(}\lambda+\frac{\delta}{\delta_{\text{A}}}\Big{)}\sum_{\eta}\eta
f^{-\eta}(\tfrac{1}{2}\lambda\delta_{\text{A}}-\mu),$ (95)

and for the non-interacting stationary current (81):

$\displaystyle I_{N}(\infty)=\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}=\tfrac{1}{2}\gamma_{p}\tfrac{1}{2}\sum_{\lambda\eta}\eta
f^{-\eta}(\tfrac{1}{2}\lambda\delta_{\text{A}}-\mu)\times\Big{(}1-\frac{\delta^{2}}{\delta_{A}^{2}}\Big{)}.$
(96)

###  7.3 Proximized dot at high bias ($\mu\gg|\alpha|,U,|\delta|,T,\Gamma$)

Finally, for an interacting and proximized dot at high bias, $\mu\gg
U,|\alpha|,|\delta|,T,\Gamma$ (for simplicity denoted as “$|\mu|\to\infty$”)
we have _approximate_ self-duality as mentioned in Sec. 6.4. In this case
$|\mu|\gg|E_{\eta,\tau}|$ such that all dependence on the Andreev energies
(and thus on $\eta$ and $\tau$) drops out in the components (74). Therefore
the two components $\kappa_{s},\kappa^{\prime}_{s}\to 0$ vanish, one is
constant $\kappa_{c}\to\gamma_{p}/2$ and one simplifies to
$\kappa^{\prime}_{c}\to\tfrac{1}{2}\gamma_{p}\sum_{\eta}\eta
f^{-\eta}(-\mu)=\tfrac{1}{2}\gamma_{p}\tanh[\mu/(2T)]$. Once more the
invariants have a simple property

$\displaystyle\gamma_{s}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\frac{\delta}{\delta_{\text{A}}}\gamma^{\prime}_{c},$
$\displaystyle\gamma^{\prime}_{s}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\frac{\delta}{\delta_{\text{A}}}\gamma_{c},$
(97)

which is very similar to the $\alpha\to 0$ case (83) except for the fact that
the “sign function” $\delta/|\delta|$ in Eq. (92) gets “broadened” by $\alpha$
in $\delta/\delta_{\text{A}}$. Note that this gate-voltage dependence of the
pre-factor, originating from the state-dependence of the effective rate
$\Gamma_{\eta\tau}$ [Eq. (28)], is not wiped out by the high bias which by our
assumptions has to remain below the (infinite) gap of the superconductor, see
Sec. 2. The remaining invariants read

$\displaystyle\gamma_{c}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\tfrac{1}{2}\gamma_{p},$
$\displaystyle\gamma^{\prime}_{c}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\tfrac{1}{2}\gamma_{p}\tanh\Big{(}\frac{\mu}{2T}\Big{)},$
(98)

such that by Eq. (62) we again have self-duality (60):

$\displaystyle\langle A\rangle_{z}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\langle\bar{A}\rangle_{\bar{z}},$
$\displaystyle\quad\langle p\rangle_{z}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\langle p\rangle_{\bar{z}}$
$\displaystyle\langle p\rangle_{z}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\langle A\rangle_{z}^{2},$
$\displaystyle\langle p\rangle_{\bar{z}}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\langle
A\rangle_{\bar{z}}^{2}.$ (99)

and the simplifications (77)-(81) apply. Now the stationary polarization
consists of factors independently controlled by the applied voltages:

$\displaystyle\langle A\rangle_{z}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}-\frac{\delta}{\delta_{\text{A}}}\tanh\Big{(}\frac{\mu}{2T}\Big{)}.$
(100)

Likewise, in the stationary current (76) the pre-factor $\kappa^{\prime}_{c}$
depends only on the bias $\mu$:

$\displaystyle I_{N}(\infty)=\gamma^{\prime}_{c}+\gamma^{\prime}_{s}\langle
A\rangle_{z}$
$\displaystyle\overset{|\mu|\to\infty}{\quad=\quad}\tfrac{1}{2}\gamma_{p}\tanh\Big{(}\frac{\mu}{2T}\Big{)}\times\Big{(}1-\frac{\delta^{2}}{\delta_{A}^{2}}\Big{)}.$
(101)

##  8 Conclusion

In this paper we have demonstrated how to fully exploit fermionic duality for
the solution of open-system dynamics and transport right from the very first
step of setting up the evolution and current equations. Unlike ordinary
symmetry, this dissipative symmetry aids the construction of a _duality-
adapted_ orthogonal Liouville-space basis which is optimal for finding the
distinct left _and_ right eigenvectors of the master-equation rate
superoperator $W$. The corresponding matrix elements have a highly constrained
functional structure governed by _duality invariance_.

We developed these ideas concretely by deriving the full time-dependent
solution (69)-(71) of the dynamics of a quantum dot proximized by a
superconductor and weakly probed by charge- and heat currents into a normal-
metal contact, a problem of high current interest. The obtained compact
expressions capture a variety of effects when interaction and induced pairing
compete, warranting a separate numerical analysis in Ref. [18]. Here we
instead analytically investigated our all-encompassing solution and also
discussed several covered limiting cases ($\alpha=0$, $U=0$, or
$|\mu|\to\infty$) which may prove useful for comparison with other approaches.
For the general case we showed that the transient approach to the stationary
state can be completely understood from duality-invariant decay rates and
expectation values of two duality-adapted observables in the initial and
_stationary_ state of the system. This unexpected result emerged by first
expressing the solution in the stationary expectation values of these
observables for the actual system _and_ for its dual system with inverted
energies, generalizing prior work [8] to include pairing.

Likewise, extending other work [13], we showed that the stationary values of
observables in the _actual_ model in fact determine their values of the _dual_
model by a stationary duality relation [Eq. (59)]. This relation is
“universal” for this class of system, independent of the values of all
physical parameters: temperature $T$, coupling $\Gamma$, voltage $\mu$, but
also the –attractive or repulsive– interaction $U$ and induced superconducting
pairing $\alpha$. This naturally suggested a notion of _self-duality_ of these
observables which we showed occurs essentially whenever the interaction is
irrelevant. This happens when interaction is either explicitly zero or made
ineffective asymptotically and we illustrated these special cases.

When the interaction is relevant self-duality is violated but duality remains
valid and is intimately linked with interactions within the open system. In
particular, we showed that its _sign_ of the interaction $U$
(repulsive/attractive) equals the sign of one of the duality invariants
($\gamma_{C}$), irrespective of the induced superconducting pairing $\alpha$.
Likewise, the duality invariants clarified the singular connection for
$\alpha\to 0$ to the problem of an interacting quantum dot plus normal metal
without superconductor.

Our results underscore that fermionic duality is a powerful tool for advancing
the quantitative understanding not only of the decay rates / time-scales of
electronic of open quantum systems but importantly also of the amplitudes
which decide their (ir)relevance depending on the initial state. It is an
intriguing open question how to systematically extend our work to even more
complicated transport models within the broad class governed by fermionic
duality. The continued interest in (time-)controlled proximized nanoelectronic
systems [58, 59, 60, 61, 62] provides a strong experimental impetus for such
work.

## Acknowledgements

We thank A. Geresdi, F. Hassler, and J. Schulenborg for discussions.

#### Funding information

L. C. O. acknowledges support by the Deutsche Forschungsgemeinschaft (RTG
1995). J.S. acknowledges financial support from the Swedish VR and the Knut
and Alice Wallenberg Foundation.

##  Appendix A Stationary duality and Kolmogorov detailed balance

Here we present two derivations of the stationary duality relation (59).

_Direct derivation:_ Using Eqs. (43) we first express the parity $-1$
probability $\tfrac{1}{2}[1-\langle p\rangle_{z}]$ in terms of $\langle
A\rangle_{z}$ and $\langle A\rangle_{\bar{z}}$ by eliminating
$(\tfrac{1}{2}\gamma_{p}+\gamma_{C})/\gamma_{p}=[1-\langle
A\rangle_{z}/\langle A\rangle_{\bar{z}}]^{-1}$ and then express the remainder
in $\langle A\rangle_{z}$. One then similarly expresses the parity $+1$
probability:

$\displaystyle\tfrac{1}{2}\big{[}1-\langle p\rangle_{z}\big{]}=\frac{1-\langle
A\rangle_{z}^{2}}{1-\langle A\rangle_{z}/\langle
A\rangle_{\bar{z}}},\qquad\tfrac{1}{2}\big{[}1+\langle
p\rangle_{z}\big{]}=\frac{1-\langle A\rangle_{z}\langle
A\rangle_{\bar{z}}}{1-\langle A\rangle_{\bar{z}}/\langle A\rangle_{z}},$ (102)
$\displaystyle\tfrac{1}{2}\big{[}1-\langle
p\rangle_{\bar{z}}\big{]}=\frac{1-\langle A\rangle_{\bar{z}}^{2}}{1-\langle
A\rangle_{\bar{z}}/\langle A\rangle_{z}},\qquad\tfrac{1}{2}\big{[}1+\langle
p\rangle_{\bar{z}}\big{]}=\frac{1-\langle A\rangle_{z}\langle
A\rangle_{\bar{z}}}{1-\langle A\rangle_{z}/\langle A\rangle_{\bar{z}}}.$ (103)

The second row is obtained in a similar way and corresponds to formally
replacing $z\to\bar{z}$ in the first row. From the first (second) relation in
first (second) row one then finds the dual polarization (parity) in terms of
the actual polarization and parity as given by Eq.(59a).

_Derivation using Eq. (74) of Ref.[13]._ Fermionic duality was combined in
Ref. [13] with the assumption that detailed balance
$P_{i}/P_{j}=W_{ij}/W_{ji}$ holds, thus presupposing a unique stationary state
with strictly positive probabilities and the validity of Kolmogorov condition
[3]. These assumptions apply to the system considered here, see Ref. [13] for
details. There it was shown that this implies a “universal” fermionic duality
relation between the probabilities ($P_{i}$) of the energy eigenstates in the
stationary state of the system and of the dual system ($\bar{P}_{i}$),
separately counting possible degenerate states:

$\displaystyle\bar{P}_{i}=\frac{P_{i}^{-1}}{\sum_{k}P_{k}^{-1}}.$ (104)

This simple relation established an interesting connection to “relative
stationary rareness” via Kac’s lemma, see Ref. [13] for further discussion. In
our notation these relations connect the energy-basis probabilities [Eq. (16)]
of $\bm{|}z\bm{)}=\sum_{\tau}z_{\tau}\bm{|}\tau\bm{)}+z_{1}\bm{|}1\bm{)}$ and
$\bm{|}\bar{z}\bm{)}=\sum_{\tau}\bar{z}_{\tau}\bm{|}\bar{\tau}\bm{)}+\bar{z}_{1}\bm{|}1\bm{)}$
of which we eliminate $z_{1}$ and $\bar{z}_{1}$ by normalization. Accounting
for the degeneracy of the $N=1$ spin states Eq. (104) gives for the Andreev
states $\tau=\pm$

$\displaystyle\bar{z}_{\tau}$
$\displaystyle=\frac{(z_{-\tau})^{-1}}{\sum_{\tau}z_{\tau}^{-1}+2(z_{1}/2)^{-1}}$
$\displaystyle=\frac{z_{1}}{z_{1}\sum_{\tau}z_{\tau}+4z_{+}z_{-}}\,z_{{\tau}}=Fz_{\tau}.$
(105)

with rescaling factor (59b)

$\displaystyle F$
$\displaystyle=\frac{z_{1}}{z_{1}\sum_{\tau}z_{\tau}+4z_{+}z_{-}}=\frac{\tfrac{1}{2}\big{[}1-\langle
p\rangle_{z}\big{]}}{\tfrac{1}{2}\big{[}1+\langle p\rangle_{z}\big{]}-\langle
A\rangle_{z}^{2}}.$ (106)

Changing to duality-adapted variables $z_{\tau}=\tfrac{1}{4}\big{[}1+\langle
p\rangle_{z}\big{]}+\tau\tfrac{1}{2}\langle A\rangle_{z}$ and
$z_{1}=\tfrac{1}{2}\big{[}1-\langle p\rangle_{z}\big{]}$, we obtain relation
(59a)

$\displaystyle\langle\bar{A}\rangle_{\bar{z}}$
$\displaystyle\equiv\sum_{\tau}\tau\bar{z}_{\tau}$
$\displaystyle=F\sum_{\tau}\tau z_{\tau}$ $\displaystyle=F\langle
A\rangle_{z},$ (107) $\displaystyle\tfrac{1}{2}\big{[}1+\langle
p\rangle_{\bar{z}}\big{]}$ $\displaystyle\equiv\sum_{\tau}\bar{z}_{\tau}$
$\displaystyle=F\sum_{\tau}z_{\tau}$
$\displaystyle=F\tfrac{1}{2}\big{[}1+\langle p\rangle_{z}\big{]},$ (108)

Although this derivation highlights the importance of detailed balance for the
“universality” of this relation, it still relies crucially on fermionic
duality through Eq. (104).

##  Appendix B Condition for self-duality

Here we consider the self-duality condition
$\gamma_{c}=\tfrac{1}{2}\gamma_{p}$ for the charge decay rate [Eq.
(11),(49a)]:

$\displaystyle\gamma_{c}=\tfrac{1}{2}\sum_{\eta\tau}W^{\eta}_{1,\tau}=\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}\Big{(}1+\lambda\frac{\delta}{\delta_{\text{A}}}\Big{)}\sum_{\eta}f^{-\eta}\Big{(}\tfrac{1}{2}(\lambda\delta_{\text{A}}+\eta
U)-\mu\Big{)}.$ (109)

###  B.1 Exact self-duality and its breakdown due to interaction

The strict lower / upper bounds (63),
$\gamma_{c}\gtrless\tfrac{1}{2}\gamma_{p}$ for $U\gtrless 0$ at finite
temperature $T$ follow using
$f(x)=[e^{x/T}+1]^{-1}=\tfrac{1}{2}\big{[}1+(1-e^{x/T})/(1+e^{x/T})\big{]}$
and $f^{\eta}(x)=f(\eta x)$ :

$\displaystyle\sum_{\eta}f^{-\eta}\big{(}x+\eta
y)=1+\frac{1-e^{-2y/T}}{(e^{-(x+y)/T}+1)(e^{(x-y)/T}+1)}\gtrless 1.$ (110)

With $x=\lambda\delta_{\text{A}}-\mu$ and $y=U/2$, Eq. (109) gives
$\gamma_{c}\gtrless\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}(1+\lambda\frac{\delta}{\delta_{\text{A}}})=\tfrac{1}{2}\gamma_{p}$
for $U\gtrless 0$. For $T<\infty$ the condition for self-duality (60),
$\gamma_{c}=\gamma_{p}/2$ [Eq. (62)], holds for all values of parameters other
than $U$ if and only if $U=0$ since Eq. (110) is an equality if and only if
$y=0$. For $T\to\infty$ equality is satisfied trivially for all parameter
values.

###  B.2 Asymptotic self-duality

It is possible to achieve self-duality asymptotically, i.e., by making one
energy scale dominate all others (excluding $\Gamma$ which cancels out
$\gamma_{c}=\gamma_{p}/2$). All that is needed to eliminate the dependence on
$\eta$ in the argument of the Fermi distribution in Eq. (109). This is always
possible by taking high temperature $T\gg|U|,|\alpha|,|\delta|$

$\displaystyle\gamma_{c}\approx\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}\Big{(}1+\lambda\frac{\delta}{\delta_{\text{A}}}\Big{)}\sum_{\eta}f^{-\eta}\big{(}0\big{)}=\tfrac{1}{2}\gamma_{p}.$
(111)

This “trivial” limit in fact plays a crucial role in the derivation of
fermionic duality [8, 9], in the underlying renormalized perturbation theory
[17] and methods based on this [63, 16, 64]. If temperature is not dominant,
to achieve self-duality one needs to make ineffective the $\eta$ dependence
which is tied to the interaction $U$ in Eq. (109) in the transition energy
$E_{\eta,\eta\lambda}=\tfrac{1}{2}(\lambda\delta_{\text{A}}+\eta U)-\mu$. For
large bias voltage $|\mu|\gg|\delta|,|\alpha|,U$ such that
$|\mu|\gg|E_{\eta,\eta\lambda}|$ discussed in Sec. 7.3 this works:

$\displaystyle\gamma_{c}\approx\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}\Big{(}1+\lambda\frac{\delta}{\delta_{\text{A}}}\Big{)}\sum_{\eta}f^{-\eta}\big{(}-\mu\big{)}=\tfrac{1}{2}\gamma_{p}.$
(112)

Similarly, for large pairing $|\alpha|\gg|\delta|$ relative to detuning we
recover the $U\to 0$ limit discussed in Sec. 7.2:

$\displaystyle\gamma_{c}\approx\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}\Big{(}1+\lambda\frac{\delta}{\alpha}\Big{)}\sum_{\eta}f^{-\eta}\big{(}\lambda\tfrac{1}{2}|\alpha|-\mu\big{)}=\tfrac{1}{2}\gamma_{p}.$
(113)

For large detuning $|\delta|\gg|\alpha|,U$ one recovers

$\displaystyle\gamma_{c}\approx\tfrac{1}{2}\gamma_{p}\sum_{\lambda}\tfrac{1}{2}\Big{(}1+\lambda\frac{\delta}{|\delta|}\Big{)}\sum_{\eta}f^{-\eta}\big{(}\lambda\tfrac{1}{2}|\delta|-\mu\big{)}=\tfrac{1}{2}\gamma_{p}$
(114)

which we discussed in Sec. 7.1 for $|\delta|\neq 0$ taking $\alpha\to 0$ and
then $U\to 0$.

## References

  * [1] B. Buča and T. Prosen,  _A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains_ ,  New J. Phys. 14(7), 073007 (2012),  10.1088/1367-2630/14/7/073007.
  * [2] T. Prosen,  _Pt-symmetric quantum liouvillean dynamics_ ,  Phys. Rev. Lett. 109(9), 090404 (2012),  10.1103/PhysRevLett.109.090404.
  * [3] A. Kolmogoroff,  _Zur Theorie der Markoffschen Ketten_ ,  Math. Ann. 112(1), 155 (1936),  10.1007/BF01565412.
  * [4] H. J. Carmichael and D. F. Walls,  _Detailed balance in open quantum Markoffian systems_ ,  Z. Phys. B: Condens. Matter 23(3), 299 (1976),  10.1007/BF01318974.
  * [5] G. S. Agarwal,  _Open quantum Markovian systems and the microreversibility_ ,  Z. Phys. A: Hadrons Nucl. 258(5), 409 (1973),  10.1007/BF01391504.
  * [6] R. Alicki,  _On the detailed balance condition for non-hamiltonian systems_ ,  Rep. Math. Phys. 10(2), 249 (1976),  10.1016/0034-4877(76)90046-X.
  * [7] R. Serfozo,  _Basics of Applied Stochastic Processes_ ,  Springer, Berlin, Germany,  ISBN 978-3-540-89331-8,  10.1007/978-3-540-89332-5 (2009).
  * [8] J. Schulenborg, R. B. Saptsov, F. Haupt, J. Splettstoesser and M. R. Wegewijs,  _Fermion-parity duality and energy relaxation in interacting open systems_ ,  Phys. Rev. B 93(8), 081411 (2016),  10.1103/PhysRevB.93.081411.
  * [9] V. Bruch, K. Nestmann, J. Schulenborg and M. Wegewijs,  _Fermionic duality: General symmetry of open systems with strong dissipation and memory_ ,  SciPost Phys. 11(3), 053 (2021),  10.21468/SciPostPhys.11.3.053.
  * [10] L. D. Contreras-Pulido, J. Splettstoesser, M. Governale, J. König and M. Büttiker,  _Time scales in the dynamics of an interacting quantum dot_ ,  Phys. Rev. B 85(7), 075301 (2012),  10.1103/PhysRevB.85.075301.
  * [11] J. Vanherck, J. Schulenborg, R. B. Saptsov, J. Splettstoesser and M. R. Wegewijs,  _Relaxation of quantum dots in a magnetic field at finite bias – Charge, spin, and heat currents_ ,  Phys. Status Solidi B 254(3), 1600614 (2017),  10.1002/pssb.201600614.
  * [12] J. Schulenborg, A. Di Marco, J. Vanherck, M. R. Wegewijs and J. Splettstoesser,  _Thermoelectrics of Interacting Nanosystems—Exploiting Superselection Instead of Time-Reversal Symmetry_ ,  Entropy 19(12), 668 (2017),  10.3390/e19120668.
  * [13] J. Schulenborg, J. Splettstoesser and M. R. Wegewijs,  _Duality for open fermion systems: Energy-dependent weak coupling and quantum master equations_ ,  Phys. Rev. B 98(23), 235405 (2018),  10.1103/PhysRevB.98.235405.
  * [14] J. Schulenborg, M. R. Wegewijs and J. Splettstoesser,  _Thermovoltage in quantum dots with attractive interaction_ ,  Appl. Phys. Lett. 116(24), 243103 (2020),  10.1063/5.0008866.
  * [15] J. Schulenborg,  _Dynamics of open fermionic nano-systems – a fundamental symmetry and its application to electron transport in interacting quantum dots_ ,  Ph.D. thesis, Chalmers University of Technology, Sweden (2018).
  * [16] R. B. Saptsov and M. R. Wegewijs,  _Fermionic superoperators for zero-temperature nonlinear transport: Real-time perturbation theory and renormalization group for Anderson quantum dots_ ,  Phys. Rev. B 86(23), 235432 (2012),  10.1103/PhysRevB.86.235432.
  * [17] R. B. Saptsov and M. R. Wegewijs,  _Time-dependent quantum transport: Causal superfermions, exact fermion-parity protected decay modes, and Pauli exclusion principle for mixed quantum states_ ,  Phys. Rev. B 90(4), 045407 (2014),  10.1103/PhysRevB.90.045407.
  * [18] L. C. Ortmanns, M. R. Wegewijs and J. Splettstoesser,  Manuscript in preparation.
  * [19] V. Meden,  _The Anderson–Josephson quantum dot—a theory perspective_ ,  J. Phys.: Condens. Matter 31(16), 163001 (2019),  10.1088/1361-648x/aafd6a.
  * [20] M. Kamp and B. Sothmann,  _Phase-dependent heat and charge transport through superconductor–quantum dot hybrids_ ,  Phys. Rev. B 99(4), 045428 (2019),  10.1103/PhysRevB.99.045428.
  * [21] A. Braggio, M. Governale, M. G. Pala and J. König,  _Superconducting proximity effect in interacting quantum dots revealed by shot noise_ ,  Solid State Commun. 151(2), 155 (2011),  10.1016/j.ssc.2010.10.043.
  * [22] S. Droste, J. Splettstoesser and M. Governale,  _Finite-frequency noise in a quantum dot with normal and superconducting leads_ ,  Phys. Rev. B 91(12), 125401 (2015),  10.1103/PhysRevB.91.125401.
  * [23] L. Rajabi, C. Pöltl and M. Governale,  _Waiting Time Distributions for the Transport through a Quantum-Dot Tunnel Coupled to One Normal and One Superconducting Lead_ ,  Phys. Rev. Lett. 111(6), 067002 (2013),  10.1103/PhysRevLett.111.067002.
  * [24] B. Wang and J. Wang,  _Optimal quantum pump in the presence of a superconducting lead_ ,  Phys. Rev. B 66(20), 201305 (2002),  10.1103/PhysRevB.66.201305.
  * [25] M. Blaauboer,  _Charge pumping in mesoscopic systems coupled to a superconducting lead_ ,  Phys. Rev. B 65(23), 235318 (2002),  10.1103/PhysRevB.65.235318.
  * [26] Y. Wei and J. Wang,  _Carbon-nanotube-based quantum pump in the presence of a superconducting lead_ ,  Phys. Rev. B 66(19), 195419 (2002),  10.1103/PhysRevB.66.195419.
  * [27] J. Splettstoesser, M. Governale, J. König, F. Taddei and R. Fazio,  _Pumping through a quantum dot in the proximity of a superconductor_ ,  Phys. Rev. B 75(23), 235302 (2007),  10.1103/PhysRevB.75.235302.
  * [28] B. Hiltscher, M. Governale, J. Splettstoesser and J. König,  _Adiabatic pumping in a double-dot Cooper-pair beam splitter_ ,  Phys. Rev. B 84(15), 155403 (2011),  10.1103/PhysRevB.84.155403.
  * [29] L. Arrachea and R. López,  _Anomalous Joule law in the adiabatic dynamics of a quantum dot in contact with normal-metal and superconducting reservoirs_ ,  Phys. Rev. B 98(4), 045404 (2018),  10.1103/PhysRevB.98.045404.
  * [30] R. Žitko, J. S. Lim, R. López and R. Aguado,  _Shiba states and zero-bias anomalies in the hybrid normal-superconductor Anderson model_ ,  Phys. Rev. B 91(4), 045441 (2015),  10.1103/PhysRevB.91.045441.
  * [31] K. Wrześniewski, B. Baran, R. Taranko, T. Domański and I. Weymann,  _Quench dynamics of a correlated quantum dot sandwiched between normal-metal and superconducting leads_ ,  Phys. Rev. B 103(15), 155420 (2021),  10.1103/PhysRevB.103.155420.
  * [32] R. Taranko, K. Wrześniewski, B. Baran, I. Weymann and T. Domański,  _Transient effects in a double quantum dot sandwiched laterally between superconducting and metallic leads_ ,  Phys. Rev. B 103(16), 165430 (2021),  10.1103/PhysRevB.103.165430.
  * [33] M. G. Pala, M. Governale and J. König,  _Nonequilibrium Josephson and Andreev current through interacting quantum dots_ ,  New J. Phys. 9(8), 278 (2007),  10.1088/1367-2630/9/8/278.
  * [34] M. Governale, M. G. Pala and J. König,  _Real-time diagrammatic approach to transport through interacting quantum dots with normal and superconducting leads_ ,  Phys. Rev. B 77(13), 134513 (2008),  10.1103/PhysRevB.77.134513.
  * [35] B. Sothmann, D. Futterer, M. Governale and J. König,  _Probing the exchange field of a quantum-dot spin valve by a superconducting lead_ ,  Phys. Rev. B 82(9), 094514 (2010),  10.1103/PhysRevB.82.094514.
  * [36] J. Eldridge, M. G. Pala, M. Governale and J. König,  _Superconducting proximity effect in interacting double-dot systems_ ,  Phys. Rev. B 82(18), 184507 (2010),  10.1103/PhysRevB.82.184507.
  * [37] M. Kamp and B. Sothmann,  _Higgs-like pair amplitude dynamics in superconductor–quantum-dot hybrids_ ,  Phys. Rev. B 103(4), 045414 (2021),  10.1103/PhysRevB.103.045414.
  * [38] A. Fornieri and F. Giazotto,  _Towards phase-coherent caloritronics in superconducting circuits_ ,  Nat. Nanotechnol. 12(10), 944 (2017),  10.1038/nnano.2017.204.
  * [39] F. Hajiloo, F. Hassler and J. Splettstoesser,  _Mesoscopic effects in the heat conductance of superconducting-normal-superconducting and normal-superconducting junctions_ ,  Phys. Rev. B 99(23), 235422 (2019),  10.1103/PhysRevB.99.235422.
  * [40] L. González Rosado, F. Hassler and G. Catelani,  _Weak localization corrections to the thermal conductivity in $s$-wave superconductors_,  Phys. Rev. B 101(17), 174202 (2020),  10.1103/PhysRevB.101.174202.
  * [41] J. P. Pekola and B. Karimi,  _Colloquium: Quantum heat transport in condensed matter systems_ ,  Rev. Mod. Phys. 93(4), 041001 (2021),  10.1103/RevModPhys.93.041001.
  * [42] M. Heckschen and B. Sothmann,  _Pair-amplitude dynamics in strongly coupled superconductor–quantum dot hybrids_ ,  Phys. Rev. B 105(4), 045420 (2022),  10.1103/PhysRevB.105.045420.
  * [43] G. Cheng, M. Tomczyk, S. Lu, J. P. Veazey, M. Huang, P. Irvin, S. Ryu, H. Lee, C.-B. Eom, C. S. Hellberg and J. Levy,  _Electron pairing without superconductivity_ ,  Nature 521(7551), 196 (2015),  10.1038/nature14398,  25971511.
  * [44] G. Cheng, M. Tomczyk, A. B. Tacla, H. Lee, S. Lu, J. P. Veazey, M. Huang, P. Irvin, S. Ryu, C.-B. Eom, A. Daley, D. Pekker _et al._ ,  _Tunable Electron-Electron Interactions in ${\mathrm{LaAlO}}_{3}/{\mathrm{SrTiO}}_{3}$ Nanostructures_,  Phys. Rev. X 6(4), 041042 (2016),  10.1103/PhysRevX.6.041042.
  * [45] G. E. D. K. Prawiroatmodjo, M. Leijnse, F. Trier, Y. Chen, D. V. Christensen, M. von Soosten, N. Pryds and T. S. Jespersen,  _Transport and excitations in a negative-U quantum dot at the LaAlO3/SrTiO3 interface_ ,  Nat. Commun. 8(395), 1 (2017),  10.1038/s41467-017-00495-7.
  * [46] A. Hamo, A. Benyamini, I. Shapir, I. Khivrich, J. Waissman, K. Kaasbjerg, Y. Oreg, F. von Oppen and S. Ilani,  _Electron attraction mediated by Coulomb repulsion_ ,  Nature 535(7612), 395 (2016),  10.1038/nature18639.
  * [47] I. Affleck, J.-S. Caux and A. M. Zagoskin,  _Andreev scattering and Josephson current in a one-dimensional electron liquid_ ,  Phys. Rev. B 62(2), 1433 (2000),  10.1103/PhysRevB.62.1433.
  * [48] Y. Tanaka, A. Oguri and A. C. Hewson,  _Kondo effect in asymmetric Josephson couplings through a quantum dot_ ,  New J. Phys. 9(5), 115 (2007),  10.1088/1367-2630/9/5/115.
  * [49] T. Meng, S. Florens and P. Simon,  _Self-consistent description of Andreev bound states in Josephson quantum dot devices_ ,  Phys. Rev. B 79(22), 224521 (2009),  10.1103/PhysRevB.79.224521.
  * [50] C. Karrasch, S. Andergassen and V. Meden,  _Supercurrent through a multilevel quantum dot close to singlet-triplet degeneracy_ ,  Phys. Rev. B 84(13), 134512 (2011),  10.1103/PhysRevB.84.134512.
  * [51] A. V. Rozhkov and D. P. Arovas,  _Interacting-impurity Josephson junction: Variational wave functions and slave-boson mean-field theory_ ,  Phys. Rev. B 62(10), 6687 (2000),  10.1103/PhysRevB.62.6687.
  * [52] J. Monsel, J. Schulenborg, T. Baquet and J. Splettstoesser,  _Geometric energy transport and refrigeration with driven quantum dots_ ,  Phys. Rev. B 106(3), 035405 (2022),  10.1103/PhysRevB.106.035405.
  * [53] R. Fazio and R. Raimondi,  _Resonant Andreev Tunneling in Strongly Interacting Quantum Dots_ ,  Phys. Rev. Lett. 80(13), 2913 (1998),  10.1103/PhysRevLett.80.2913.