Controlling Quantum Transport via Dissipation Engineering

TL;DR

This paper presents an open-system framework for controlling quantum transport in strongly interacting fermionic systems through dissipation engineering, enabling manipulation of subgap currents, nonreciprocal transport, and other phenomena.

Contribution

It introduces a novel approach to dissipative control of quantum transport in strongly interacting systems, applicable to solid state and cold atom experiments.

Findings

Controlled subgap currents with Multiple Andreev Reflections.

Engineered dissipation can induce nonreciprocal currents.

Different microscopic origins of currents lead to distinct behaviors.

Abstract

Inspired by the microscopic control over dissipative processes in quantum optics and cold atoms, we develop an open-system framework to study dissipative control of transport in strongly interacting fermionic systems, relevant for both solid state and cold atom experiments. We show how subgap currents exhibiting Multiple Andreev Reflections -- the stimulated transport of electrons in the presence of Cooper-pairs -- can be controlled via engineering of superconducting leads or superfluid atomic gases. Our approach incorporates dissipation within the channel, which is naturally occurring and can be engineered in cold gas experiments. This opens opportunities for engineering many phenomena with transport in strongly interacting systems. As examples, we consider particle loss and dephasing, and note different behaviour for currents with different microscopic origin. We also show how to induce nonreciprocal electron and Cooper-pair currents.

Figures

Equations

$$H_{\mathrm{QD}}=\sum _{s=\downarrow ,\uparrow }\omega c_s^{†}{c}_s+U{c}_{\uparrow }^{†}{c}_{\uparrow }{c}_{\downarrow }^{†}{c}_{\downarrow },$$

$$H_{\mathrm{QD}}^{\mathrm{eff}}(t)=H_{\mathrm{QD}}+\sum_{\ell=L,R}\big{(}g_{\ell}e^{2iV_{\ell}t}c_{\downarrow}c_{\uparrow}+\mathrm{h.c.}\big{)},$$

$$\displaystyle H_{\mathrm{QD}}^{\mathrm{eff}}(t)=\sum_{s=\downarrow,\uparrow}\omega c_{s}^{\dagger}c_{s}+Uc_{\uparrow}^{\dagger}c_{\uparrow}c_{\downarrow}^{\dagger}c_{\downarrow}+\sum_{\ell}\big{(}g_{\ell}e^{2iV_{\ell}t}c_{\downarrow}c_{\uparrow}+\mathrm{h.c.}\big{)},$$

$$\displaystyle H_{\mathrm{lead}}=\sum_{\ell=L,R}\sum_{k}\Big{[}\sum_{s=\downarrow,\uparrow}\left(\omega_{k}+V_{\ell}\right)b_{\ell ks}^{\dagger}b_{\ell ks}+\sum_{k^{\prime}}U_{kk^{\prime}}b_{\ell k\uparrow}^{\dagger}b_{\ell k\downarrow}^{\dagger}b_{\ell-k^{\prime}\uparrow}b_{\ell-k^{\prime}\downarrow}\Big{]},$$

$$H_{\mathrm{I}}=\sum _{\ell =L,R}{\kappa }_{\ell }\sum _{ks}\left(b_{\ell ks}^{†}{c}_s+h.c.\right).$$

$$\displaystyle H_{\mathrm{I}}^{(p)}=\sum_{\ell}\big{(}g_{\ell}S^{\dagger}_{\ell}c_{\downarrow}c_{\uparrow}+\mathrm{h.c.}\big{)},$$

$${\dot{\rho }}_{\mathrm{tot}}^I(t)=-i\left[H_{\mathrm{I}}^I(t),{\rho }_{\mathrm{tot}}^I(t)\right]$$

$$H_I^I(t)=\sum _{\ell }{\kappa }_{\ell }\sum _{ks}\left(b_{\ell ks}^{I†}(t)c_s^I(t)+h.c.\right),$$

$$c_s^I(t)=U(t{)}^{†}{c}_{s}U(t),$$

$$b_{\ell ks}^I(t)=e^{i{H}_{\mathrm{B}}t}{b}_{\ell ks}{e}^{-i{H}_{\mathrm{B}}t},$$

$$U(t)=\mathcal{T}{e}^{-i{\int }_0^t{H}_{\mathrm{QD}}^{\mathrm{eff}}(t^{\prime })d{t}^{\prime }},$$

$$H_{\mathrm{B}}=\sum _{\ell }\sum _{ks}\left({\tilde{\omega }}_{\ell k}{\tilde{b}}_{\ell ks}^{†}{\tilde{b}}_{\ell ks}+V_{\ell }{b}_{\ell ks}^{†}{b}_{\ell ks}\right),$$

$$b_{\ell -k\downarrow }=u_{\ell k}{\tilde{b}}_{\ell k\downarrow }+v_{\ell k}{\tilde{b}}_{\ell k\uparrow }^{†},$$

$$b_{\ell k\uparrow }=u_{\ell k}{\tilde{b}}_{\ell k\uparrow }-v_{\ell k}{\tilde{b}}_{\ell k\downarrow }^{†},$$

$$u_{\ell k}=e^{-i{\varphi }_{\ell }}\sqrt{\frac{1}{2}\left(1+\frac{{\omega }_k}{{\tilde{\omega }}_{\ell k}}\right)},$$

$$v_{\ell k}=\sqrt{\frac{1}{2}\left(1-\frac{{\omega }_k}{{\tilde{\omega }}_{\ell k}}\right)}.$$

$$b_{\ell k\downarrow }^I(t)=\left(u_{\ell k}{\tilde{b}}_{\ell -k\downarrow }{e}^{-i{\tilde{\omega }}_{\ell k}t}+v_{\ell k}{\tilde{b}}_{\ell -k\uparrow }^{†}{e}^{i{\tilde{\omega }}_{\ell k}t}\right)e^{-i{V}_{\ell }t},$$

$$b_{\ell k\uparrow }^I(t)=\left(u_{\ell k}{\tilde{b}}_{\ell k\uparrow }{e}^{-i{\tilde{\omega }}_{\ell k}t}-v_{\ell k}{\tilde{b}}_{\ell k\downarrow }^{†}{e}^{i{\tilde{\omega }}_{\ell k}t}\right)e^{-i{V}_{\ell }t}.$$

$${\rho }_{\mathrm{tot}}^I(t)\approx {\rho }^I(t)\otimes {\rho }_L\otimes {\rho }_R$$

$${\rho }_{\ell }=\frac{e^{-{\beta }_{\ell }\left(H_B-V_{\ell }{N}_{\ell }\right)}}{\mathrm{Tr}\left[e^{-{\beta }_{\ell }\left(H_B-V_{\ell }{N}_{\ell }\right)}\right]}$$

$$⟨{\tilde{b}}_{\ell ks}^{†}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}{⟩}_B\equiv {\mathrm{Tr}}_B\left[{\tilde{b}}_{\ell ks}^{†}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}{\rho }_{\ell }\right]={\delta }_{k{k}^{\prime }}{\delta }_{s{s}^{\prime }}{n}_{\ell }({\tilde{\omega }}_{\ell k}),$$

$$⟨{\tilde{b}}_{\ell ks}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}^{†}{⟩}_B\equiv {\mathrm{Tr}}_B\left[{\tilde{b}}_{\ell ks}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}^{†}{\rho }_{\ell }\right]={\delta }_{k{k}^{\prime }}{\delta }_{s{s}^{\prime }}(1-n_{\ell }({\tilde{\omega }}_{\ell k})),$$

$$⟨{\tilde{b}}_{\ell ks}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}{⟩}_B\equiv {\mathrm{Tr}}_B\left[{\tilde{b}}_{\ell ks}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}{\rho }_{\ell }\right]=⟨{\tilde{b}}_{\ell ks}^{†}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}^{†}{⟩}_B\equiv {\mathrm{Tr}}_B\left[{\tilde{b}}_{\ell ks}^{†}{\tilde{b}}_{\ell k^{\prime }{s}^{\prime }}^{†}{\rho }_{\ell }\right]=0,$$

$${\dot{\rho }}^I(t)=-i\mathrm{T}{\mathrm{r}}_{\mathrm{B}}\left(\left[H_I^I(t),{\rho }_{\mathrm{tot}}(0)\right]\right)-{\int }_0^{t}d{t}^{\prime }\mathrm{T}{\mathrm{r}}_{\mathrm{B}}\left(\left[H_I^I(t),\left[H_I^I(t-t^{\prime }),{\rho }_{\mathrm{tot}}(t-t^{\prime })\right]\right]\right).$$

$$\displaystyle\dot{\rho}(t)=-\sum_{\ell}\sum_{s,s^{\prime}}\int_{0}^{t}\bigg{[}\langle B_{\ell s}^{\dagger}(t)B_{\ell s^{\prime}}(t-t^{\prime})\rangle_{B}c_{s}(t)c_{s^{\prime}}^{\dagger}(t-t^{\prime})\rho(t)+\langle B_{\ell s}(t)B^{\dagger}_{\ell s^{\prime}}(t-t^{\prime})\rangle_{B}c_{s}^{\dagger}(t)c_{s^{\prime}}(t-t^{\prime})\rho(t)\bigg{.}$$

$$-⟨B_{\ell s^{\prime }}(t-t^{\prime })B_{\ell s}^{†}(t){⟩}_B{c}_s(t)\rho (t-t^{\prime })c_{s^{\prime }}^{†}(t-t^{\prime })-⟨B_{\ell s^{\prime }}^{†}(t-t^{\prime })B_{\ell s}(t){⟩}_B{c}_s^{†}(t)\rho (t)c_{s^{\prime }}(t-t^{\prime })$$

$$-⟨B_{\ell s}(t)B_{\ell s^{\prime }}^{†}(t-t^{\prime }){⟩}_B{c}_{s^{\prime }}(t-t^{\prime })\rho (t)c_s^{†}(t)-⟨B_{\ell s}^{†}(t)B_{\ell s^{\prime }}(t-t^{\prime }){⟩}_B{c}_{s^{\prime }}^{†}(t-t^{\prime })\rho (t)c_s(t)$$

$$\displaystyle\hskip 85.35826pt\bigg{.}+\langle B_{\ell s^{\prime}}^{\dagger}(t-t^{\prime})B_{\ell s}(t)\rangle_{B}\rho(t)c_{s^{\prime}}(t-t^{\prime})c_{s}^{\dagger}(t)+\langle B_{\ell s^{\prime}}(t-t^{\prime})B^{\dagger}_{\ell s}(t)\rangle_{B}\rho(t)c_{s^{\prime}}^{\dagger}(t-t^{\prime})c_{s}(t)\bigg{]}dt^{\prime}$$

$$\displaystyle\quad\quad\quad+\sum_{\ell}\sum_{s,s^{\prime}}\int_{0}^{t}\bigg{[}\langle B_{\ell s}^{\dagger}(t)B_{\ell s^{\prime}}^{\dagger}(t-t^{\prime})\rangle_{B}c_{s}(t)c_{s^{\prime}}(t-t^{\prime})\rho(t)+\langle B_{\ell s}(t)B_{\ell s^{\prime}}(t-t^{\prime})\rangle_{B}c_{s}^{\dagger}(t)c_{s^{\prime}}^{\dagger}(t-t^{\prime})\rho(t)\bigg{.}$$

$$-⟨B_{\ell s^{\prime }}^{†}(t-t^{\prime })B_{\ell s}^{†}(t){⟩}_B{c}_s(t)\rho (t)c_{s^{\prime }}(t-t^{\prime })-⟨B_{\ell s^{\prime }}(t-t^{\prime })B_{\ell s}(t){⟩}_B{c}_s^{†}(t)\rho (t)c_{s^{\prime }}^{†}(t-t^{\prime })$$

$$-⟨B_{\ell s}^{†}(t)B_{\ell s^{\prime }}^{†}(t-t^{\prime }){⟩}_B{c}_{s^{\prime }}(t-t^{\prime })\rho (t)c_s(t)-⟨B_{\ell s}(t)B_{\ell s^{\prime }}(t-t^{\prime }){⟩}_B{c}_{s^{\prime }}^{†}(t-t^{\prime })\rho (t)c_s^{†}(t)$$

$$\displaystyle\hskip 85.35826pt\bigg{.}+\langle B_{\ell s^{\prime}}^{\dagger}(t-t^{\prime})B_{\ell s}^{\dagger}(t)\rangle_{B}\rho(t)c_{s^{\prime}}(t-t^{\prime})c_{s}(t)+\langle B_{\ell s^{\prime}}(t-t^{\prime})B_{\ell s}(t)\rangle_{B}\rho(t)c_{s^{\prime}}^{\dagger}(t-t^{\prime})c_{s}^{\dagger}(t)\bigg{]}dt^{\prime},$$

$$⟨B_{\ell \downarrow }^{†}(t)B_{\ell \downarrow }(t-t^{\prime }){⟩}_B$$

$$=\frac{{\gamma }_{\ell }}{\pi }{\int }_{-\infty }^{\infty }d\omega e^{i{V}_{\ell }{t}^{\prime }}\left(|u_{\ell }(\omega ){|}^2{e}^{i{\tilde{\omega }}_{\ell }{t}^{\prime }}n({\tilde{\omega }}_{\ell })+|v_{\ell }(\omega ){|}^2{e}^{-i{\tilde{\omega }}_{\ell }{t}^{\prime }}[1-n({\tilde{\omega }}_{\ell })])\right)$$

$$=\frac{{\gamma }_{\ell }}{\pi }{\int }_0^{\infty }d\omega e^{i{V}_{\ell }{t}^{\prime }}\left(e^{i{\tilde{\omega }}_{\ell }{t}^{\prime }}n({\tilde{\omega }}_{\ell })+e^{-i{\tilde{\omega }}_{\ell }{t}^{\prime }}[1-n({\tilde{\omega }}_{\ell })])\right)$$

$$=\frac{{\gamma }_{\ell }}{\pi }{\int }_0^{\infty }d\tilde{\omega }{e}^{i{V}_{\ell }{t}^{\prime }}{D}_{\ell }(\tilde{\omega })\left(e^{i\tilde{\omega }{t}^{\prime }}n(\tilde{\omega })+e^{-i\tilde{\omega }{t}^{\prime }}[1-n(\tilde{\omega })])\right)$$

$$=\frac{{\gamma }_{\ell }}{\pi }{\int }_{-\infty }^{\infty }d\tilde{\omega }{D}_{\ell }(\tilde{\omega })e^{i(\tilde{\omega }+V_{\ell })t^{\prime }}n(\tilde{\omega })$$

$$D_{\ell }(E)=\frac{|E|}{\sqrt{E^2-{\Delta }_{\ell }^2}}\Theta (|E|-{\Delta }_{\ell })$$

$$f_{\ell +}(t^{\prime })$$

$$f_{\ell -}(t^{\prime })$$

Full text

††preprint: APS/123-QED

Introduction.

Understanding and controlling the out-of-equilibrium dynamics of strongly interacting many-body systems constitutes one of the key forefronts in quantum physics across a variety of subfields in experiment and theory. In this context, opportunities to achieve their control via dissipation mechanisms have arisen Mul12 ; Dal14 , as is applied for few-body systems in quantum optics Wis09 ; Gar15 . This is especially true in cold-atom platforms, where large separations between frequency scales allows well-controlled theoretical models and implementations of dissipative processes, as realized for laser cooling and trapping Met01 . The longer timescales of cold atom experiments also allow dynamics to be tracked and potentially controlled time-dependently Blo08 ; Blo12 . Out-of-equilibrium transport dynamics remain a ubiquitous paradigm in the solid state Bee91 , and recent developments in cold atom systems have also made it possible to engineer quantised transport of atoms between reservoirs, as well as quantum point contacts and waveguides Kri15 ; Kri17 ; Hus15 ; Leb18 . Here we explore the emerging new opportunity of using dissipation engineering to achieve control of quantum transport properties, that are relevant for both cold-atom and solid-state platforms.

We study transport in a system of strongly interacting fermions coupled to weakly interacting reservoirs, as can be reaslied with cold atoms using optical tweezers connecting larger superfluids, or with solid-state devices using quantum dots (QD) coupled to superconducting leads (S). In a traditional S-QD-S tunelling junction, subgap transport is known to be suppressed for weak electron tunneling as compared to the gap of the attached leads Def10 . Here we demonstrate that subgap transport can be recovered even in the regime of weak tunnelling. This is done via reservoir engineering that allows for independent control of Cooper-pair and single-electron channels. Such channel separation can be accomplished in the solid state by adding two large-gap superconductors to a traditional S-QD-S junction, producing a four terminal structure, or in cold atoms considering driving from a molecular Bose-Einstein condensate Jia11 .

Subgap currents in this context are produced by Multiple Andreev Reflections (MARs) And64 ; Bee97 ; Mar11 ; Def10 , i.e., stimulated transport of electrons via exchange of Cooper-pairs. MARs have been observed in the solid state Oct83 ; Pil10 ; Bui03 ; Dir11 and cold atoms Hus15 , and their signatures can be used to reveal topological phase transitions related to Majorana bound states formation San13 . We show how to engineer well-resolved MAR peaks under weak electron tunneling, and show how these behave in the presence of dissipation in the channel - providing a diagnostic tool for the microscopic nature of the current. We also show that for asymmetric coupling, the reciprocity of the engineered system is broken, yielding electron and Cooper pair currents dependent on the bias direction. This represents a genuinely new way of generating nonreciprocal transport of electrons and Cooper-pairs.

We investigate the transport properties of the junction with an open system approach, while most of the theoretical works rely on Keldish non-equilibrium Green functions or scattering techniques. These approaches are able to treat the tunnelling rate $\gamma$ between the QD and the leads non-perturbatively, but usually treat the Coulomb interaction $U$ between the QD electrons perturbatively or within a mean-field treatment Yey97 ; Zaz06 ; Del08 ; Kan98 . In contrast, open system approaches such as input-output theories Sea02 ; Gar04 ; Gar04b ; Zha11 or master equations Bre06 ; Kos13 ; Pfa13 work well in the opposite regime: for arbitrary interaction $U$ but weak tunnelling rate $\gamma$, implying that MARs have been left beyond their scope. In our framework, the large gap superconducting leads behave effectively as time dependent coherent drives of Cooper-Pairs on the QD (analogous to laser fields in quantum optics). This dynamical model is naturally cast as a dissipative Floquet system, for which we derive a Floquet-Born-Markov master equation Yan16 ; Koh96 ; Gra94 ; Bla15 capturing MARs up to arbitrary order. Our open-system framework provides an opportunity to study the effects of controlled or uncontrolled dissipation acting on the QD. We thus analyze the response of the currents to fermion losses and dephasing, and show in particular robustness of the currents against dephasing. We use in the following natural units in which $\hbar=k_{B}=e=1$, where $-e$ is the electron charge.

Model.

To represent the separate control of Cooper pair driving, we consider a four-terminal QD connected to two pairs of left (L) and right (R) superconducting leads by tunnel junctions, as depicted in Fig. 1. In each pair, we consider one lead in the single-particle mean-field description with a moderate energy gap $\Delta_{\ell}$ ($\ell=L,R$), and one described only by its condensed fraction of Cooper-pairs, assuming that the gap is so large that single-particle excitations are irrelevant. A bias voltage $V=V_{L}-V_{R}$ is generated between the pairs of superconductors, where $V_{L}$ and $V_{R}$ are the voltages of each side. The QD Hamiltonian reads

[TABLE]

and describes electrons of spin $s$, energy $\omega$, and Coulomb interaction $U$. The QD is an effective 4-level system spanned by the non-occupied, single occupied, and double-occupied states $\left\{|0\rangle,|\downarrow\rangle,|\uparrow\rangle,|\downarrow\uparrow\rangle\right\}$. The coupling of the QD to the large-gap superconducting leads (red superconductors in Fig. 1) gives rise to a pairing of the QD electrons, i.e., the proximity effect Mar11 , and results in an effective time-dependent QD Hamiltonian of the form

[TABLE]

where $g_{\ell}$ is the Cooper-pair tunnelling amplitude between the QD and the large-gap superconducting $\ell=L,R$. Hence, the coupling of the large-gap superconductors with the QD takes the form of a driving of the transition between the non-occupied and double-occupied states $|0\rangle$ and $|\downarrow\uparrow\rangle$ of the QD.

We obtain the dissipative dynamics of the QD by coupling the Hamiltonian (2) to the superconductors with moderate gaps $\Delta_{\ell}$ (blue superconductors in Fig. 1) under an open system approach, by deriving a Floquet-Born-Markov master equation Yan16 ; Koh96 ; Gra94 ; Bla15 for the QD. The leads, considered in a mean-field single-particle description, act as baths of Bogoliubov quasiparticles of density of states $D_{\ell}(E)\propto\Theta(|E|-\Delta_{\ell})|E|/(\sqrt{E^{2}-\Delta_{\ell}^{2}})$. The tunnelling of electrons between the leads and the QD is described by a standard tunnelling Hamiltonian of the form $H_{\mathrm{int}}=\sum_{\ell}\kappa_{\ell}\sum_{ks}(b^{\dagger}_{\ell ks}c_{s}+\mathrm{h.c.})$, where $\kappa_{\ell}$ is the electron tunnelling amplitude and $b_{\ell ks}$ the annihilation operator of an electron of spin $s$ and momentum $k$ in the moderate-gap lead $\ell$. The derivation of the master equation, in second-order in $H_{\mathrm{int}}$, results in a single-particle tunnelling rate $\gamma_{\ell}\propto\kappa_{\ell}^{2}$ (the typical line widths of the QD levels) considered as the smallest parameter. Note that while treating perturbatively the single-particle coupling, the master equation describes the QD Coulomb interaction $U$ exactly, as in Kos13 ; Pfa13 . See Supplemental Material for details of the derivation supmat .

Engineering of transport.

From the solutions of the master equation, we calculate the particle current in the leads as a function of the applied bias voltages, taken as opposite from each other for the sake of simplicity ($V_{L}=-V_{R}=V/2$). Figure 2 shows the particle current $I$ in both the moderate and large-gap right leads as a function of $V$. We consider the electron-hole symmetric case $\omega=-U/2$, and vanishing or not Cooper-pair tunnelling $g\equiv g_{\ell}$, taken here real and identical for left and right leads. We also consider identical single-particle tunnelling rates $\gamma\equiv\gamma_{\ell}$ between the QD and the moderate-gap superconducting leads.

When $g=0$, the large-gap superconducting leads are disconnected from the QD, and our system simply consists in a conventional S-QD-S tunnelling junction Kos13 ; Pfa13 . Only one peak of current is observed (see panel A in Fig. 2), whose the shape is related to the superconductor density of states $D_{\ell}(E)\propto\Theta(|E|-\Delta_{\ell})|E|/(\sqrt{E^{2}-\Delta_{\ell}^{2}})$. The peak appears when $E_{i}>\omega>E_{f}$, where $E_{i}$ is the energy of the highest occupied state of the left lead and $E_{f}$ is the energy of the lowest non-occupied state of the right lead (see energy diagram I in Fig. 2). For high bias, the particle current tends to the value $2\gamma$ of normal leads. For low bias ($\gamma\ll V<4\Delta_{\ell}$), i.e. in the subgap region (where no resonance between left-lead occupied and right-lead non-occupied states exist), no current is observed as a result of the weak coupling approximation. Indeed, for $\gamma\ll\Delta$, Andreev reflection at the interface with the moderate-gap superconductor is negligible.

Connecting the large-gap superconducting leads to the QD (i.e., setting $g\neq 0$) allows Andreev reflections to occur. Under such process, an electron (hole) is reflected as a hole (electron) producing the emission (absorption) of a Cooper-pair in the large-gap superconducting leads (see panel B in Fig. 2). After some reflections, electrons of the QD acquire enough energy to tunnel into the moderate-gap superconducting lead. This produces well-resolved single-particle subgap currents more and more pronounced as $g$ increases. These processes are represented in our Floquet-Born-Markov formalism by decay channels corresponding to QD transition (quasi)energies shifted by multiple of Cooper-pair energies (see Supplemental Material). The subgap currents are located at $V=2(|\omega|+\Delta)/(2n+1)$ where $n=1,2,\dots$ denotes the $n^{th}$ MAR (see energy diagrams II and III corresponding respectively to the first and second Andreev reflections). This can be obtained from the condition $E_{i}+nV=\omega=E_{f}-nV$, in which $n$ denotes the number of Cooper-pairs transfer from the left to the right lead. Note however that in general the bias voltage at which a MAR peak appears is a function of both the QD charging energy $U$ and $\omega$ (see Supplemental Material). Hence, while the tunnelling between the QD and the moderate-gap leads is always sequential due to the weak-coupling regime (one electron per time), it can be assisted by transfer of an arbitrary number of Cooper-pairs between the large-gap superconductors, thanks to the stronger tunnelling amplitude $g$. This represents reservoir engineering of subgap single-particle currents. This is our first important result. Interestingly, the Cooper-pair current in the right large-gap superconducting lead is negative outside the subgap region. We attribute this phenomenon to a supercurrent (i.e., Cooper-pair current) reversal, due to the modification of parity of the QD when the voltage exceed the value delimiting the subgap border, as can be seen in the panel C in Fig. 2 Def10 ; Del18 . Note that the sign and amplitude of the supercurrent are dependent of the phases of the superconductors (not shown).

Effects of particle loss and dephasing.

In the previous section, we showed that dissipation induced by reservoir engineering can be used to control subgap transport. Here we examine the robustness of the produced subgap currents against the presence of incoherent processes, that are inherent in real experimental setups. We incorporate these effects into our master equation through an additional dissipator of the form $\mathcal{D}_{I}(\rho)=\gamma_{I}\left(2L\rho L^{\dagger}-\left\{L^{\dagger}L,\rho\right\}\right)$, where $\gamma_{I}$ is the rate of the incoherent process and $L$ the corresponding Lindblad operator (see Supplemental Material). For cold atoms experiments, the dissipation in the channel is often in this Markovian form as can be derived from first principles Dal14 .

We first consider the effects of particle loss (i.e., $\gamma_{I}\equiv\gamma_{\mathrm{loss}}$, $L=c_{s}$) acting on the QD. This occurs naturally in the cold-atom platforms through background gas collisions, and could be engineered using electron beams Ger08 or light scattering quantum gas microscopes with single-site resolution Wei11 ; Bak09 ; She10 (analogous to x-ray scattering in the solid state). In Fig. 3 (panels A-D), we show the particle currents in all the leads as a function of the bias voltage for increasing loss rates $\gamma_{\mathrm{loss}}$. The presence of losses results in competing effects. On the one hand, the additional decay channel tends to empty the QD faster. This results in an increase (decrease) of the currents of electrons entering (reaching) the moderate-gap superconducting leads. On the other hand, pushing the QD towards the non-occupied state $|0\rangle$ increases the effects of the driving (since the driving only affects the QD in the non-occupied or double-occupied states), which favors MARs and thus raises subgap currents. Hence, while source currents (panels A and C) only increase due to electrons losses, drain currents (panels B and D) are subjected to these competing effects, exhibiting amplitude increase or decrease depending on the voltage bias.

We then consider the effects of dephasing (i.e., $\gamma_{I}\equiv\gamma_{\mathrm{deph}}$, $L=c^{\dagger}_{s}c_{s}$) acting on the QD, which occurs naturally through coupling to additional degrees of freedom in the solid state, and can be engineered in cold atoms through light scattering or noise  Dal14 ; Lus17 ; Sar14 ; Pic10 ; Nie17 . We show that dephasing acting on the QD affects identically the source and drain Cooper-pair currents, whereas leaves unchanged the electron currents. Figure 3 (panels E) shows the current of Cooper-pair leaving the QD to reach the right large-gap superconductor for different dephasing rate $\gamma_{\mathrm{deph}}$. Our results show that increasing the dephasing rate reduces the size of the subgap peaks (see panel F). This can be understood as a consequence of the blurring of the QD energy levels caused by the dephasing. Hence, dephasing tend to destroy Cooper-pair subgap currents, but does not affect the single particle currents. This suggests these latter are robust against phonon/photon scattering in condensed matter/cold atomic systems.

Nonreciprocal subgap transport.

Finally, we show how to generate nonreciprocal subgap transport. Bias-direction-dependent properties is generally a desired feature of nanoscale devices, and are known to result from the presence of asymmetry and nonlinearity. Nonreciprocal transport at the quantum level has been investigated in spin Zala ; Landi1 ; Landi ; Emmanuel ; Emmanuel2 ; Emmanuel3 ; Emmanuel4 ; Valente ; Mas18 and QD systems Mal18 ; Ono02 ; Joh05 ; Tan18 ; Sch08 ; Mal18 . This includes the paradigmatic Pauli Blockade effects in a double-QD junction, where a nonreciprocal electron current has been observed for asymmetric QD energy levels Ono02 . For a single QD, the required asymmetry can be provided by different left and right tunnelling rates. In a S-QD-S junction, in the intermediate coupling regime ($\gamma_{\ell}\sim\Delta_{\ell}$), non reciprocal conductance has been observed and explained as originating from asymmetric Kondo resonance at the contact with the leads Eic07 . Here we show that for asymmetric weak single-particle tunnelling rates $\gamma_{L}\neq\gamma_{R}$, the reciprocity of the transport properties can be broken as soon as the Cooper-pair tunnelling amplitudes $g_{\ell}$ is non-zero. In Figure 4, we plot the current-voltage characteristics for the moderate (panel A) and large gap (panel C) superconducting leads for positive and negative bias voltage (see diagrams B and D) for $\gamma_{L}=3\gamma_{R}$. While the total current (the sum of electron and Cooper-pair currents) is still reciprocal (not shown), its electron and Cooper-pair contributions become dependent of the bias direction, as can be clearly seen in Fig. 4. In particular, the current of electrons (Cooper-pair) is larger (smaller) for negative (positive) bias. We interpret this phenomenon as a Cooper-pair-assisted nonreprocical transport, since it occurs only for non-zero $g_{\ell}$. Indeed, for $g_{\ell}=0$, the electron current is reciprocal (not shown). We believe it is a genuinely new way of breaking reciprocity, since while keeping reciprocal the total current, its electron and Cooper-pair contributions – which could be measured independently in our four-terminal scheme – become asymmetric.

Conclusion.

We developed a quantum-optics-inspired framework to study the dynamics of strongly interacting fermions in tunnelling junctions under the influence of dissipation and driving, relevant for both solid-state and cold-atom platforms. For concreteness, we studied the dynamics of a QD coupled to superconducting leads in a four-terminal configuration, where two large-gap superconducting leads are added to a traditional S-QD-S tunnelling junction. We demonstrate the possibility of controlling subgap transport via dissipation engineering. We showed that the added leads generate subgap transport based on MARs despite weak electron tunnelling, and studied the effects of electron loss and dephasing acting on the QD. Finally, we showed that the Cooper-pair driving provided by the added leads is a new way of breaking the reciprocity of the junction, generating nonreciprocal electron and Cooper-pair subgap currents based on MARs.

Our results could be investigated in both solid-state and cold-atom experiments. They could be generalized to multi-QD tunnelling junction, and to include the presence of measurement and feedback loop to control the transport dynamics of fermions in tunnelling junctions Uch18 . More possible outlooks include reservoir engineering of (Floquet)-Majorana fermions Jia11 ; Ben15 ; Li014 , or studies of the interplay between dissipation and driving in thermodynamics problems such as thermoelectric effects RSan13 or quantum heat engines Pek15 involving superconductors.

I SUPPLEMENTAL MATERIAL

In this Supplemental Material, we derive the Floquet-Born-Markov master equation used to investigate the non-equilibrium dynamics of our four-tunnelling junction. We then explain how to solve the master equation and compute the different particle currents, and how to incorporate the effects of additional incoherent processes acting on the QD. We provide some details of our numerical simulations, and finally show the dependence of the positions of the MAR peaks on the QD energy.

II Derivation of the master equation

We model the junction with the Hamiltonian $H=H_{\mathrm{QD}}^{\mathrm{eff}}(t)+H_{\mathrm{lead}}+H_{\mathrm{I}}$ where

[TABLE]

In Eq. (3), $H_{\mathrm{QD}}^{\mathrm{eff}}(t)$ is the effective time-dependent QD Hamiltonian which describes the intrinsic dynamics of the QD and the contribution of the left and right large-gap superconducting leads. The time-dependent pairing terms can be viewed as coming from an interaction Hamiltonian of the form

[TABLE]

where $S_{\ell}$ is the annihilation operator of a Cooper-pair of energy $2V_{\ell}$ in the large-gap superconducting $\ell$, on which one applies a semi-classical description by considering the fields $S_{\ell}$ in coherent states, making the replacement $S_{\ell}\to e^{-2iV_{\ell}t}$, and rescaling $g_{\ell}$. Note that we consider here $g_{\ell}$ arbitrary, since it can be tuned by adjusting the properties of the large-gap superconducting leads and their coupling with the QD. The Hamiltonian $H_{\mathrm{lead}}$ describes the moderate-gap superconducting leads where $b_{\ell ks}$ is the annihilation operator of an electron of spin $s$ and energy $\omega_{k}$, $V_{\ell}$ the bias voltage, and $U_{kk^{\prime}}$ the interaction between electrons. Finally, the Hamiltonian $H_{\mathrm{I}}$ describes electron tunnelling between the QD and the leads with amplitudes $\kappa_{\ell}$. Our four-terminal model allows us by essence to consider independent electron and Cooper-pair tunnelling. We note however that independent tunnelling can also be the result of modelling superconducting leads with two components, as done in Su17 .

We consider the system-bath (S-B) decomposition $H_{\mathrm{S}}(t)\equiv H_{\mathrm{QD}}^{\mathrm{eff}}(t)$ and $H_{\mathrm{B}}\equiv H_{\mathrm{lead}}$, assuming the electron tunnelling amplitude $\kappa_{\ell}$ to be the smallest frequency scale. In order to take into account the time-dependence of the system Hamiltonian in the dissipative processes, we derive a Floquet-Born-Markov master equation, i.e., a Born-Markov master equation based on the Floquet quasienergy spectrum of $H_{\mathrm{S}}(t)$ Yan16 ; Koh96 ; Gra94 . This procedure (i.e. adding the time-dependent driving first and then calculating dissipation) is justified since the coupling with the large-gap superconductors is larger than the ones with the moderate-gap leads. This allows the dissipation to be assisted by Cooper-pair exchanges, yielding the subgap peak structures observed in the main text. The reverse procedure (i.e. calculating dissipation first and then adding the time-dependent driving) is unable to describe the interplay between the electrons and the Cooper-pairs, and is anyway not justified in our regime of parameters.

II.1 Interaction picture and Bogoliubov approximation

Our starting point is the usual Liouville-Von Neumann equation Bre06

[TABLE]

for the total density matrix $\rho_{\mathrm{tot}}^{I}(t)$ in interaction picture with respect to $H_{0}(t)=H_{\mathrm{S}}(t)+H_{\mathrm{B}}$, where

[TABLE]

with

[TABLE]

where the propagator $U(t)$ is defined as

[TABLE]

with $\mathcal{T}$ the time-ordering operator.

In order to calculate $b_{\ell ks}^{I}(t)$, we first evaluate the time-dependence due to the bias voltages $\sum_{\ell}V_{\ell}N_{\ell}=\sum_{\ell ks}V_{\ell}b_{\ell ks}^{\dagger}b_{\ell ks}$ where $N_{\ell}$ is the total electron number operator in the lead $\ell$, which yields a factor $e^{-iV_{\ell}t}$. This can be done since $\sum_{\ell}V_{\ell}N_{\ell}$ and $H_{B}-\sum_{\ell}V_{\ell}N_{\ell}$ commute. We then obtain the time-dependence provided by the remaining Hamiltonian $H_{B}-\sum_{\ell}V_{\ell}N_{\ell}$ by performing the Bogoliubov mean-field approximation on $H_{B}-\sum_{\ell}V_{\ell}N_{\ell}$. For this purpose, we define the order parameter $\Delta_{\ell}e^{i\phi_{\ell}}=-\sum_{k}^{\prime}U_{kk^{\prime}}\langle b_{\ell-k^{\prime}\downarrow}b_{\ell k^{\prime}\uparrow}\rangle$ where $\Delta_{\ell}$ is the energy gap of the lead $\ell$ and then diagonalize $H_{B}-\sum_{\ell}V_{\ell}N_{\ell}$ which yields

[TABLE]

where $\tilde{b}_{\ell,ks}$ describes the annihilation of a Bogoliubov quasiparticle of spin $s$ and energy $\tilde{\omega}_{\ell k}=\sqrt{\omega_{k}^{2}+\Delta_{\ell}^{2}}$ obtained from the Bogoliubov transformation

[TABLE]

where

[TABLE]

The time-dependence of the lead electron operators thus read

[TABLE]

II.2 Born and Markov approximation

We consider the Born approximation, i.e., the total system-bath density matrix $\rho_{\mathrm{tot}}^{I}(t)$ in the separable form

[TABLE]

where $\rho^{I}$ is the QD density matrix in interaction picture and $\rho_{\ell}$ a stationary state of the lead $\ell$, taken here as the thermal state

[TABLE]

so that the Fermi levels lie in the middle of the gaps of their respective leads, where $\beta_{\ell}=1/k_{B}T_{\ell}$ with $T_{\ell}$ their temperature. In doing so, we obtain the simple expressions

[TABLE]

where $n_{\ell}(\tilde{\omega}_{\ell k})$ is the Fermi occupation number defined as $n_{\ell}(E)=1/(1+e^{\beta_{\ell}E})$.

Expanding Eq. (5) to the second order in $H_{\mathrm{I}}$, tracing over the bath degrees of freedom, we obtain

[TABLE]

Using the Born approximation (14), we first assume that the first term vanishes and then perform the Markov approximation by setting $\rho^{I}(t-\tau)\approx\rho^{I}(t)$ and extending the upper limit of integration to infinity. Expanding the double commutator of Eq. (17), we obtain

[TABLE]

where we removed the superscript I to not burden the notation and where we define $B_{\ell s}(t)=\kappa_{\ell}\sum_{k}b_{\ell ks}(t)$.

II.3 Bath correlation functions

The first summation of the right-hand-side of Eq. (18) contains non-vanishing terms only for $s=s^{\prime}$. They involve bath correlations (16) with coefficient $|u_{\ell k}|^{2}$ or $|v_{\ell k}|^{2}$. Explicitly, the first term for $s=s^{\prime}=\downarrow$ involves the correlation

[TABLE]

In Eq. (19), $\gamma_{\ell}=\pi\kappa_{\ell}^{2}\rho_{N}$ is the tunnelling rate between the QD and the lead $\ell$, where $\rho_{N}=\mathcal{V}mk_{\ell F}/2\pi^{2}\hbar^{2}$ is the density of states of a normal lead assumed to be constant around the Fermi level, where $\mathcal{V}$ is the volume of the lead $\ell$, $m$ the electron mass and $k_{\ell F}$ the Fermi wavenumber, as in Pfa13 ; and

[TABLE]

is the adimensional density of states of a superconducting lead, which naturally appears when making the substitution $\omega\to\tilde{\omega}=\sqrt{\omega^{2}+\Delta_{\ell}^{2}}$ in the integration. A quick look shows that the right-hand-side of Eq. (18) contains only the following two distinct correlations functions

[TABLE]

All the other correlations functions are up to a minus sign these correlations functions or their complex conjugates.

The second summation of the right-hand-side of Eq. (18) contains non-vanishing terms only for $s\neq s^{\prime}$. They also involve bath correlations (16) but with coefficients $uv$ or $u_{\ell k}^{*}v_{\ell k}^{*}$. There are only two distinct correlations functions $f_{\ell 1}(t^{\prime})$ and $f_{\ell 2}(t^{\prime})$ (all the others are up to a minus sign these correlations functions or their conjugates) which read

[TABLE]

Note that while the correlations (21) depend only on the time difference $t^{\prime}$, Eq. (22) explicitly depend on time $t$.

II.4 Case $g_{\ell}=0$ - Standard Born-Markov master equation

Before deriving the general Floquet-Born-Markov equation for the QD, we first investigate the case of $g_{\ell}=0$ for which the QD Hamiltonian is time-independent. This corresponds to the case of a standard two-terminal junction, where the QD is only coupled to two superconducting leads. This yields a standard Born-Markov master equation, similar to the one derived in Kos13 ; Pfa13 .

The QD Hamiltonian $H_{\mathrm{QD}}=\sum_{s=\downarrow,\uparrow}\omega c_{s}^{\dagger}c_{s}+Uc_{\uparrow}^{\dagger}c_{\uparrow}c_{\downarrow}^{\dagger}c_{\downarrow}$ can be diagonalised in the basis $\left\{|0\rangle,|\downarrow\rangle,|\uparrow\rangle,|\downarrow\uparrow\rangle\right\}$ by writing the operators $c_{s}$ in terms of the projection operators $\sigma_{ij}=|i\rangle\langle j|$ ($i,j=0,\downarrow,\uparrow,\downarrow\uparrow$) satisfying $\sigma_{ij}\sigma_{kl}=\delta_{jk}\sigma_{il}$,

[TABLE]

as in Zha11 . Note that we chose the convention $|\downarrow\uparrow\rangle=c_{\downarrow}^{\dagger}c_{\uparrow}^{\dagger}|0\rangle$ instead of $|\downarrow\uparrow\rangle=c_{\uparrow}^{\dagger}c_{\downarrow}^{\dagger}|0\rangle$ regarding the order of action of creation operators to find their matrix representation (otherwise the sign in front of $|s\rangle\langle\downarrow\uparrow|$ would be reversed). The system operators in interaction picture with respect to $H_{0}$ then simply reads

[TABLE]

From Eq. (18), using Eq (24) and the expressions of the bath correlation functions (21) and (22), we obtain the master equation in the standard form (coming back to the Schrödinger picture)

[TABLE]

where

[TABLE]

is the Lamb-shift Hamiltonian coming from standard single particle tunnelling processes, where

[TABLE]

is a time-dependent Lamb-shift Hamiltonian coming from Cooper-pair tunnelling processes, where

[TABLE]

is a time-independent single-particle dissipator coming from standard single-particle tunnelling processes; and where

[TABLE]

is a time-dependent Cooper-pair dissipator coming from Cooper-pair tunnelling processes. Note that the two first terms on the right-hand side do not have the anticommutator since $\sigma_{s\downarrow\uparrow}\sigma_{0s}=0$ and $\sigma_{0s}^{\dagger}\sigma_{s\downarrow\uparrow}^{\dagger}=0$.

The master equation involves the complex rates

[TABLE]

where $j=1,2$, where $E$ corresponds to system transition energies and where $\gamma_{\ell\pm}$ ($\gamma_{\ell j}$) and $\Omega_{\ell\pm}$ ($\Omega_{\ell j}$) are the real and imaginary parts of $\Gamma_{\ell\pm}$ ($\Gamma_{\ell j}e^{-2iV_{\ell}t}$) which explicitly read

[TABLE]

where $\mathrm{P.V.\ }$ denotes the principal value. Note that the integrands appearing in the expressions of the shifts $\Omega_{\ell\alpha}$ ($\alpha=\pm,1,2$) do not converge for $\omega\to\pm\infty$, and one has to introduce a cutoff frequency $\omega_{c}$ in the integration domain to obtain finite values for the shifts. This cutoff represents the fact that the bandwidths of real leads are finite.

Writing down explicitly the master equation for the density matrix elements $\rho_{ij}=\langle i|\rho|j\rangle$, one can see that the populations $\rho_{00},\rho_{\downarrow\downarrow},\rho_{\uparrow\uparrow},\rho_{\downarrow\uparrow\downarrow\uparrow}$ are only coupled to the coherences $\rho_{0\downarrow\uparrow}$ and $\rho_{\downarrow\uparrow 0}$ through the Lamb-shift Hamiltonian $H_{\mathrm{LS}}^{(p)}(t)$ and dissipator $\mathcal{D}^{(p)}\left(\rho,t\right)$. This corresponds to a proximity effect which vanishes for $\Delta_{\ell}=0$. Note that in Pfa13 , these terms are not present due to the use of a number-conserving Bogoliubov transformation. In Kos13 , they are present (for $U+2\omega=0$) but are time-independent, because the authors Kos13 included the bias voltage by making the replacement $\omega\to\omega-V_{\ell}$ in the end of the calculations, instead of including it directly in the lead Hamiltonian (3). We note however that these terms have an important effect only in the very low bias regime $V_{\ell}\ll\gamma_{\ell}$. For $V_{\ell}\gg\gamma_{\ell}$, one can invoke a secular approximation to neglect them. The master equation then becomes similar to the one derived by Kosov et al Kos13 : no subgap current can be observed.

II.4.1 Large gap limit

In the large gap limit $\Delta_{\ell}\gg|\omega\pm V_{\ell}|,|\omega+U\pm V_{\ell}|$, all the $\gamma_{\ell\alpha}(E)$ ($\alpha=\pm,1,2$) vanish since all the transition energies $E$ lie inside the gap. Also, $\Omega_{\ell j}(E)$ becomes independent of $E$, i.e. $\Omega_{\ell j}(E)\approx\tilde{\Omega}_{\ell}$, with

[TABLE]

As a consequence, we have

[TABLE]

which implies that $\mathcal{D}^{(p)}\left(\rho,t\right)=0$; and

[TABLE]

where $g_{\ell}=-2\tilde{\Omega}_{\ell}$. Overall, the master equation reduces to

[TABLE]

Hence, in the large gap limit, the coupling of the QD to the superconductors reduces to a driving of the transition between the non-occupied and double-occupied states of the QD, as it can be derived using other methods (see Mar11 and references therein). This justifies the use of an Hamiltonian of the form (34) in our model (3) to treat the coupling of the large-gap superconducting leads to the QD.

II.5 Case $g_{\ell}\neq 0$ - Floquet-Born-Markov master equation

We now come back to the general case of $g_{\ell}\neq 0$ (the four-terminal junction), where the QD Hamiltonian is time-dependent. In order to perform the time-integration in Eq. (18), we evaluate the time dependence of the system operators $c_{s}(t)$ given by Eq. (7) using the Floquet theory Gri98 . For that purpose, we suppose in the following that $V_{L}=-V_{R}=V/2$, so that the effective QD Hamiltonian is periodic of period $T=2\pi/V$. If it was not the case, one could simply work in the rotating-frame with respect to one of the driving frequency $2V_{\ell}$, let say $2V_{L}$. This would provide an periodic Hamiltonian of period $\delta=2(V_{R}-V_{L})$, and the same theory would apply.

Since the effective QD Hamiltonian is periodic, the QD wavefunction $|\psi(t)\rangle$ satisfying the Schrödinger equation

[TABLE]

can be written as

[TABLE]

where $|\psi_{a}(t)\rangle=e^{-iE_{a}t}|\phi_{a}(t)\rangle$ are the Floquet states with the periodic Floquet modes $|\phi_{a}(t+T)\rangle=|\phi_{a}(t)\rangle$, quasi-energies $E_{a}$, and $d_{a}=\langle\phi_{a}(0)|\psi(0)\rangle$. By definition of the propagator (9), we have

[TABLE]

showing that $e^{-iE_{a}T}$ are the eigenvalues of $U(T)$, which can be numerically computed using $U(T)\approx\prod_{n=0}^{N}e^{-iH_{\mathrm{QD}}^{\mathrm{eff}}(ndt)dt}$ with $N=T/dt-1$. Solving the eigenvalue problem (38), we obtain $E_{a,k}=E_{a}+k\frac{2\pi}{T}$ with $k\in\mathbb{Z}$, and consider the values of $E_{a,k}$ lying in the first Brillouin zone $[-\pi/T,\pi/T]$ to define the quasienergies $E_{a}$. The eigenvectors correspond to the Floquet modes at initial time $|\phi_{a}(0)\rangle$. The Floquet modes at all times $t$ are obtained from these latters using

[TABLE]

II.5.1 Master equation in the Floquet basis

We now decompose the density matrix in the Floquet mode basis $\{|\phi_{a}(0)\rangle\}$, i.e.

[TABLE]

where we restored for clarity the label I denoting the interaction picture for the density matrix. We derive below the equations of motion for the density matrix element $\rho^{I,ab}(t)\equiv\langle\phi_{a}(0)|\rho^{I}(t)|\phi_{b}(0)\rangle$ from Eq. (18).

In this basis, the matrix elements of the system operator $c_{s}(t)$ reads

[TABLE]

Since $|\phi_{a}(t)\rangle$ is periodic of period $T$, we can rewrite $\langle\phi_{a}(t)|c_{s}|\phi_{b}(t)\rangle$ in the Fourier space as

[TABLE]

which yields

[TABLE]

where $\Delta_{abk}=E_{a}-E_{b}+kV$ and

[TABLE]

Using the Floquet basis, the first term of the right-hand side of the master equation (18) reads

[TABLE]

and all other terms can be written in the same way. Hence, we see that the master equation involves the complex rates (30). All together, the master equation (18) written in the Floquet basis gives us the following set of equations for the matrix elements $\rho^{I,ab}(t)\equiv\langle\phi_{a}(0)|\rho^{I}(t)|\phi_{b}(0)\rangle$, i.e. the non-secular Floquet-Born-Markov master equation

[TABLE]

where

[TABLE]

is the Linbladian coming from standard single-particle tunnelling processes, and where

[TABLE]

is the Linbladian coming from Cooper-pair tunnelling processes.

The master equation for the matrix elements in Schrödinger picture $\rho^{ab}(t)$ can be obtained from Eq. (46) by making the replacement $\rho^{I,ab}(t)=e^{i(E_{a}-E_{b})t}\rho^{ab}(t)$. In doing so, one can see that all terms $e^{i(\Delta_{cdk^{\prime}}+\Delta_{ack})t}$ reduces to $e^{i(k+k^{\prime})Vt}$, showing that the master equation exhibits the same periodicity than the effective QD Hamiltonian. This implies that the steady state of the master equation is also periodic with the same period $T=2\pi/V$ Yud16 .

III Solution of the master equation and particle currents

The master equation (46) can be solved using different methods: either via brute force methods without any further treatment or by first writing Eq. (46) in Schrödinger picture and then exploiting its periodicity. One can indeed vectorize the density matrix as $|\rho^{S}(t)\rangle=(\rho^{S,11}(t),\rho^{S,12}(t),\dotsc,\rho^{S,44}(t))$ to obtain the equation of motion

[TABLE]

where $L(t)$ is a periodic time-dependent matrix of period $T$, and then apply again the Floquet theory on Eq. (49) Yud16 . Another solution consists in targeting directly the periodic steady state of Eq. (49) by writing it in Fourier space and solving the linear set of equations for its Fourier components Mal16 . In the main text, we solved the master equation (46) via brute force, but compared the solutions with these other methods and obtained the same results.

From the solution of the master equation (46) for the matrix elements $\rho^{I,ab}(t)$, one can rebuild the whole density matrix (40) and the expectation values of any system operator $O$ through

[TABLE]

We derive now the expressions of the particle currents $I_{\mathrm{QD}}$, $I_{\ell}^{(s)}$ and $I_{\ell}^{(p)}$ respectively in the QD, the moderate-gap and the large-gap superconducting leads $\ell$. The choice of the symbols $(s)$ and $(p)$ expresses the fact that the current in the moderate-gap superconducting leads is mainly due to single-particle ($s$) tunnelling, while the current in the large-gap superconducting leads is only due to pair ($p$) tunnelling. They are defined as

[TABLE]

where $S_{\ell}$ are the Cooper-pair annihilation operators introduced in Eq. (4) and where we used the Langevin equations for $c_{s}$, $b_{\ell ks}$ and $S_{\ell}$, i.e.

[TABLE]

Note that the total number of particles is conserved, i.e.,

[TABLE]

where the factor $2$ is front of $I_{\ell}^{(p)}$ denotes the fact that a Cooper-pair is made of two electrons.

The expectation value of the particles current in the QD is obtained from the solutions of the master equation and Eq. (50), that is

[TABLE]

Due to the conservation of particles (55), it can be related to the particles current in the leads, as shown below. Applying the derivative and using the fact that $dU(t)/dt=-iH_{\mathrm{QD}}^{\mathrm{eff}}(t)U(t)$, we get

[TABLE]

The two first terms on the right-hand-side can be rewritten as

[TABLE]

where [see Eq. (52)]

[TABLE]

Finally, replacing the derivative in Eq. (57) by the right-hand-side of the master equation (46) yields

[TABLE]

with

[TABLE]

IV Adding dissipation acting on the channel

Additional Lindblad dissipation acting on the QD can be accounted for by adding to Eq. (46) a dissipator of the form

[TABLE]

In the Floquet basis and in interaction picture with respect to $H_{\mathrm{QD}}^{\mathrm{eff}}(t)$, this dissipator reads

[TABLE]

By contrast with (46), there is no spectral dependence of the decay rates, as expected from Lindblad dissipation. By solving the master equation (46) after adding on its right-hand-side the dissipator (63), we can then calculate the associated particle current

[TABLE]

V Numerical details

Here we discuss some numerical details about the resolution of the master equation (46).

First, in order to write the master equation, we computed the quasienergies from Eq. (38) using the procedure stated above with $N=10000$, a number of time steps which insures precise enough values for the quasinergies. To avoid the divergence of the superconducting density of states (20), we modify their expression as

[TABLE]

with a small parameter $\epsilon=10^{-1}$. This procedure introduces finite values for the peaks of the DOS, and thus for the real part of the complex rates (31). Their imaginary parts require another treatment to become finite. Here, we transform the denominators appearing in the integrations as

[TABLE]

and introduce a cutoff $\Omega_{f}=100$ to the integration domain.

We solved the master equation using brute force resolution of the differential equations (46), with random or particular initial states, from a initial time $t_{i}=0$ to a final time $t_{f}=10\gamma^{-1}$, where $\gamma\equiv\gamma_{\ell}=10^{-2}$, in order to reach the steady state. We then averaged the particle currents over one period $T$.

VI QD energy dependence of the subgap currents

In this section, we show how the positions of the MAR peaks depend on the QD energy $\omega$, which can be tuned experimentally by applying a gate voltage on the QD. Figure 5 shows the conductance $dI/dV$ as a function of both the bias voltage $V$ and the QD energy $\omega$, for vanishing or not Cooper-pair tunnelling amplitude $g\equiv g_{\ell}$. For $g=0$, the conductance exhibit the standard diamond structure and vanish in the subgap region. For $g\neq 0$, subgap peaks appear and their positions depend on both the bias voltage and QD energy $\omega$. Hence, the subgap peaks do not just occur at rational fractions of the gap energy $\Delta$: they also depend on the properties of the QD.