Anomalous Josephson current through a driven double quantum dot

Josephson junctions based on quantum dots offer a convenient tunability by means of local gates. Here we analyze a Josephson junction based on a serial double quantum dot in which the two dots are individually gated by phase-shifted microwave tones of equal frequency. We calculate the time-averaged current across the junction and determine how the phase shift between the drives modifies the current-phase relation of the junction. Breaking particle-hole symmetry on the dots is found to give rise to a finite average anomalous Josephson current with phase bias between the superconductors fixed to zero. This microwave gated weak link thus realizes a tunable"Floquet $\varphi_{0}$-junction"with maximum critical current achieved for driving frequencies slightly off-resonance with the energy cost of exciting a sub-gap state on each dot. We provide numerical results supported by an analytical analysis for infinite superconducting gap and weak inter-dot coupling. We identify an interaction driven $0-\pi$ transition of anomalous Josephson current as a function of driving phase difference. Finally, we show that this junction can be tuned so as to provide for complete rectification of the time-averaged Josephson current phase relation.


I. INTRODUCTION
The Josephson junction (JJ) has become a ubiquitous device serving in a wide range of applications, including the superconducting qubits which have lead to impressive advances in quantum computing during the past two decades [1][2][3][4][5] .The weak link coupling the two superconductors can either be a plain insulating tunnel barrier, or it may exhibit internal structure like a normal region, a quantum point contact, a magnetic tunnel barrier, or a quantum dot (QD), which all host sub-gap states which may influence strongly the current phase relation (CPR) of the junction [6][7][8][9][10][11][12] .In this way, electrically gateable links like quantum dots or semiconductors offer a certain tunability of the JJ characteristics [13][14][15][16] , a feature which has been employed in the design of a hybrid gatemon [17][18][19] , adding gate control to the superconducting transmon qubit 20,21 , which has already demonstrated its efficiency in solid state quantum computing [22][23][24] .
Whereas normal Josephson junctions carry no current at zero phase bias, ϕ sc = ϕ L − ϕ R , a weak link which breaks both time-reversal and chiral symmetry may carry an anomalous Josephson current between two superconductors maintained at zero phase bias 25 .A number of proposals [25][26][27][28][29][30][31][32][33][34][35][36][37][38] have been made for such ϕ 0 junctions with an anomalous Josephson current, I(ϕ sc ) = I C sin(ϕ sc + ϕ 0 ), at least two of which have already been realized experimentally 39,40 .Of particular relevance to the present work is the proposal by Zazunov et al. 25 to use a multi-orbital QD with inter-orbital (spin-orbit) tunnelling and an external field.With such a link in the JJ, traversing electrons pick up different phases, depending on the tunnelling direction, giving rise to an anomalous Josephson current.This proposal has since been realized in an experiment by Szombati et al. 40 , using an InSb-wire QD contacted by superconducting NbTiN leads.
Here, we propose a nonequilibrium version of the multi-orbital QD considered in Ref. 25, based on the device illustrated in Fig. 1.In this Josephson junction, the two superconductors are coupled by a serial double quantum dot (DQD) where the two dots are driven by individual AC gate voltages with a common amplitude, A, and microwave frequency, Ω.This endows each of the QDs with Floquet sidebands, which play the roles of the additional spin-orbit coupled orbitals in Ref. 25.As we demonstrate below, the phase difference between the two drive voltages, θ d = θ L − θ R , can have a strong influence on the JJ CPR, and with QD levels tuned away from particle-hole symmetry it gives rise to anomalous current, which in the limit of weak tunnel couplings reduces to a simple ϕ 0 junction, with ϕ 0 = θ d .Since the time-averaged critical current is maximized when the microwave frequency is close to the energy for exciting both of the sub-gap states induced in the two proximitized quantum dots, this device comprises a nonadiabatic Cooper pair pump, or more aptly a "Floquet ϕ 0 junction".
Here we circumvent a number of these complications by replacing each of the dots with a non-interacting resonant level.Whereas this would clearly be a poor description of Coulomb blockaded QDs in many other respects, the two models do share the crucial features of the mechanism we wish to illustrate, namely the presence of sub-gap states with a strong gate dependence.As a weak link for a JJ, the resonant level model behaves much like a quantum point contact (QPC) with a CPR which reflects the phase dispersion of the sub-gap Andreev bound states (ABS) 8,44,51 .A JJ based on a Coulomb blockaded QD, however, is known to exhibit a transition from a ϕ 0 = π to ϕ 0 = 0 phase 9,10,12,14,48,49 , and the results presented below are therefore of greater relevance for a realistic Coulomb blockaded QD in its ϕ 0 = 0 phase stabilized for strong tunnel couplings 52 , or for a long normal junction with a finite dwell time 53 .
Dating back to the seminal work on photon assisted tunnelling by Tien and Gordon 54 , the problem of nonadiabatically (microwave) driven Josephson junctions has been expanded to include also junctions with QPC, QD, DQD or magnetic adatom weak links [55][56][57][58][59][60][61][62][63][64] .Experimentally, the ABS in such junctions have been measured and manipulated using microwave spectroscopy [65][66][67][68][69][70] , and these techniques are by now becoming widely available.Recently, Venitucci et al. 62 demonstrated that phase shifted microwave voltages applied to each of the superconductors in a JJ with a single resonant level as the weak link can give rise to photon assisted Cooper pair transfer and a tunable ϕ 0 -junction.Similarly, Soori et al. 71 have studied a finite-size tight-binding model of an SNS junction and found that a phase shifted drive on the two sites comprising their normal region leads to anomalous Josephson current.The model studied here is similar in spirit but not equivalent to these two studies, and our main focus is the anomalous Josephson current and the modified CPR for the JJ driven at resonance with the sub-gap states.
The paper is organized as follows.In Sec.II we present the model.In Sec.III we define the relevant Nambu-Floquet-Keldysh Green functions and provide an expres-sion for the time-averaged current to be calculated.In Sec.IV we study the limit of infinite gap, in which the main effect of the ϕ 0 -junction can be established analytically in the limit of weak inter-dot tunnel coupling.Sec.V contains the numerical results for the current and the CPR for the driven junction.Finally, the results are discussed in Sec.VI.Appendix A provides a supplementary analysis for the infinite-gap limit using Floquet theory, which allows us to also investigate the effects of local Coulomb interactions, and to confirm the rectification of the time-averaged supercurrent.

II. THE MODEL
We consider a non-interacting serial double quantum dot with on-site energies modulated by individual AC gate voltages and coupled to two (Left/Right) superconducting leads (cf.Fig. 1).The Hamiltonian reads with superconducting leads described by BCS Hamiltonians for α = L, R. The two leads are kept at the same chemical potential and are assumed to have the same gap magnitude, ∆ > 0, with different phases, ϕ L,R = ±ϕ sc /2.Both leads are represented by a featureless band-structure near a common chemical potential, i.e. ξ αk = ε αk − µ, corresponding to a common density of states, ν F , near the Fermi level.The time-dependent Hamiltonian of the double quantum dot system reads with individual AC gate voltages given as ε dα (t) = ε d + A cos(Ωt + θ α ), in terms of common (time) average energies, ε d , driving amplitudes, A, frequencies, Ω, and two independent phase constants, θ α .Here, τ i denotes the i'th Pauli matrix, τ 0 the Kronecker delta and t d is the inter-dot tunneling amplitude.The tunneling Hamiltonian reads Written in terms of Nambu spinors, ψ † αk = (c † αk↑ , c α−k↓ ) and φ † α = (d † α↑ , d α↓ ), the full Hamiltonian reads where the phase of the superconducting leads has been gauged into the tunneling matrix, T α = t α σ z e iσzϕα/2 .For simplicity, we assume below that tunneling amplitudes to the leads are real and equal, i.e.
As discussed in the introduction, we neglect the charging energies of both quantum dots altogether and consider this resonant level model as an effective model for a proximitized QD with a doublet sub-gap state.

III. KELDYSH FLOQUET GREEN FUNCTIONS
To calculate the current through the ac-driven device, we employ the non-equilibrium Green function technique [72][73][74] .Dealing with a harmonic drive, it is convenient to use Floquet Keldysh Green functions 75,76 , which offers a representation of the two-time Green functions, which, besides being convenient for numerical calculations, allows for some degree of physical interpretation of the elementary transport process in terms of Floquet side bands.The time-dependent current out of lead α for this driven junction is found as 74 where the trace is taken in Nambu space, and with Nambu/lead-matrix Green functions for the quantum dots defined as with self-energies, which are exact to second order in dotlead tunnelling, where g αk denotes the Nambu Green function in lead α.
From this self-energy, the dot Green functions can be found by solving the steady-state Dyson equations, with matrix products between Green functions implied.With a periodic drive, it is convenient to transform these two-time Green functions into Floquet matrices 76 defined with ω ∈] − Ω/2, Ω/2].This transformation presumes that the Green functions are periodic in both time arguments, with the driving period T = 2π/Ω, and thereby rests on the assumption that the system has reached a nonequilibrium steady state (NESS).In this way, the time-averaged current may be found as This is the zeroth Floquet component of the current.
Here, the Green functions and self-energies are matrices in Nambu, dot, and Floquet space and the trace is performed over all these spaces.The components of the self-energy in dot, and Floquet space are given by Σ where the momentum-summed lead Nambu Green functions are given explicitly as where n F denotes the Fermi function.Henceforth, temperature is assumed to be zero.Finally, using the Dyson equation (10), the retarded double-dot Green function is found by inverting the following infinite dimensional Floquet matrix of 4 × 4 matrices in Nambu-dot space: From the resulting retarded and advanced Green functions, the lesser function is found from Eq. ( 11) by simple matrix multiplication.

IV. INFINITE-GAP LIMIT WITH WEAK INTER-DOT TUNNEL COUPLING
It is instructive to first consider the limit of an infinite superconducting gap.This prohibits quasiparticle tunnelling altogether and transport takes place only via Cooper pairs.In the infinite-gap limit, the retarded QD self-energy becomes with Γ = πν F |t| 2 , corresponding to an effective Hamiltonian describing a proximitized quantum dot with an induced superconducting gap of Γ: with a matrix of tunnelling amplitudes given by In order to illustrate the basic microwave assisted Cooper pair transport mechanism within this infinite-gap model, we calculate here the weak-coupling tunnelling current to second order in the interdot coupling t d , given by the perturbative expression 77 : The driving enters this expression through the timedependent correlation functions, G <,> αη,αη (t, t ), describing the dynamics of the QD proximitized by lead α = L, R.
The perturbative expression for the current requires the Green's functions for t d = 0, and in this case the Hamiltonian (20) describes two independent quantum dots.It is readily diagonalized by the time-dependent Bogoliubov transformation (suppressing the QD index, α = L, R): with Nambu spinors and time-dependent unitary transformation matrix with and real coherence factors given by Notice that we omit the α = L, R subscript for clarity since it only enters in the two different phase shifts, θ α , and can readily be reinstalled.This transformation diagonalizes the Hamiltonian for each of the two different proximitized levels, and endows the quasiparticles with dynamics governed by the equation of motion, where the last term has been obtained as The corresponding transformation of the correlation functions reads The many-body eigenstates of the uncoupled and undriven QD are the empty QD, |0 , the single-electron doublet, |σ = d † σ |0 , and the doubly-occupied QD, , with energies 0, ε d , ε d and 2ε d , respectively.For the proximitized QD, the BCS-like ground state becomes | 0 = u|0 + v|2 , the excited doublet, |σ remains unchanged, and the highest excited state becomes | 2 = u|2 −v|0 , with energies 0, E d , E d , and 2E d , respectively (cf.Fig. 2a).
In contrast to electric charge, the total DQD parity (odd/even number of quasiparticles) is conserved in the infinite-gap limit.Switching between even-parity states, | 0/ 2 , and odd-parity states, |σ , therefore takes place exclusively by inter-dot tunnelling with amplitude t d .If the undriven system is prepared in its even-parity ground state, this will give rise to a finite Josephson current to second order in t d , found from Eq. ( 21) to be In the driven case, a similar formula for the timeaveraged current valid to second order in t d can be obtained when the system is driven with low amplitude, A E d , close to resonance, i.e.Ω ≈ 2E d , as indicated in Fig. 2b.Since the mixing term ( 28) is already proportional to driving amplitude, A, we shall neglect the time dependence of E d (t) in its denominator and in coherence factors u and v, assuming that A max( d , Γ).This allows us to include the mixing term (28) within a rotating wave approximation (RWA), which leads to the following equation of motion where g = AΩΓ/(2E d ) 2 .This equation is readily solved by with a secondary unitary transformation as defined in terms of Here, δ = E d − Ω/2 is the detuning, and the energy Ẽ = δ 2 + g 2 captures the slow time evolution of the co-rotating Nambu spinor with initial condition ζ µ (0) = e iσ z νν θ/2 Ũµν χ ν (0).For concreteness, we assume that driving is turned on at time t = 0, prior to which each proximitized QD is assumed to be thermalized in its ground state, | 0 .Using the relations, the time-evolved states are found as Reinstating the lead index α on θ and inserting this into Eqs.( 29), one finally arrives at the correlation functions with the phase factor in front introduced merely for convenience in formulas below.The time-dependent current in Eq. ( 21) may now be expressed as involving products like + uv e −i(Ωt+θ L +θ R )/2 (g/ Ẽ) 2 sin 2 ( Ẽt) The last two terms contain fast oscillating phase factors, e ±iΩt/2 , and will be strongly suppressed by the subsequent integrations with respect to t.Notice that it is only these fast oscillating terms which contain information about the average phases of the drives, i.e. the information about the initial time, chosen here as t = 0, from which the initial states are being time-evolved.Retaining only the first slowly oscillating term, which only carries information about the phase difference, θ d , these products reduce to with Finally, introducing κ(θ d , ϕ sc ) = f (θ d )e iϕsc/2 = κ + iκ , the current takes the following form × 2 Ẽt sin(2 Ẽt) − sin 2 (2 Ẽt) − 4 sin 4 ( Ẽt) .
Whereas the first term oscillates around zero, the two last terms give rise to a well-defined long-time average lim which results in the time-averaged current valid to leading order in t d , close to resonance, Ω ≈ 2E d , and with the phase shifted by This current vanishes at resonance, Ω = 2E d , and attains its maximum for Ω = 2E d ± AΓ/( √ 3E d ) with maximum current given by which is not strictly valid, since the maximum is attained where δ ∼ g, and counter rotating terms no longer are negligible.
Tuning the levels away from the Fermi level, i.e. for |ε d | Γ, we have E d ≈ ε d and the current becomes which is the Floquet θ d -junction, in which the phase shift of the sinusoidal current phase relation is set directly by the phase shift of the two driving fields together with the sign of the level energies set by ε d .
In the opposite limit, where the two levels are close to the Fermi levels of the two superconducting leads, i.e. ε d Γ, we have E d ≈ Γ and arrive at which corresponds to a 0-junction right above resonance where Ω > 2Γ, and a π-junction right below resonance for Ω < 2Γ.In this limit, the phase shift of the two drives, θ d , serves only to modulate the amplitude, attaining maximum average current when the drives are shifted by θ d = π, and blocking it altogether for θ d = 0.This average current was calculated under the assumption of an even number of electrons occupying each of the two levels, with the specific initial condition that the system is in its lowest energy state at time zero.In a real system, however, quasiparticle poisoning, and relaxation will cause occasional switching of the parity of each of the two levels.With typical parity flip times of the order of 20-200 µs 66,78-80 , a resonant drive frequency, Ω ∼ 2E d , of the order of 10 GHz, say, will take the system through some 10 6 cycles before the parity is flipped, implying that the average current determined here remains meaningful as long as Ẽ 10 −6 Ω, i.e. √ AΓ 10 MHz.The full problem thus entails the stochastic element of random parity switching between even, and odd parity sectors of the Hilbert space.This poses an interesting problem in itself, but shall not be pursued any further in this work.Instead, we shall analyze the steady state Dyson equation (11), in which the parity is relaxed in the infinite gap limit by a weak coupling to a normal metallic reservoir.For a finite gap, the Floquet sidebands of the continuum provide the same effect and the normal metallic reservoir is no longer needed.
Notice that the full lesser component of the Dyson equation has a second contribution 74 , (1+G R Σ R )G < 0 (1+ Σ A G A ), referring to the initial lesser function, and that this term has been omitted altogether in equation (11).This omission rests on the tacit assumption, that Σ < contains relaxation mechanisms, which will wash out the initial conditions, i.e. that Σ < G R,−1 0 )f 0 , where f 0 denotes an initial distribution function.In the present tunnelling problem, Σ < refers to quasiparticle tunnelling to and from the superconducting leads and to the weak tunnelling of electrons directly between the dots and a normal metal reservoir.The former contribution vanishes altogether in the infinite-gap limit, and the steady state Dyson equation (11) as well as the Floquet Keldysh transformation (12), is therefore justified in the infinite gap limit by the normal metal tunneling rate, Γ m , which is large enough to dominate the finite η = (G R,−1 0 − G A,−1 0 )/(2i) used in our numerical implementation of the bare Green functions of the leads, yet small enough not to affect the result.

V. NUMERICAL RESULTS
In this section we present numerical results obtained with the Floquet Keldysh Green functions introduced in Sec.III.We shall focus entirely on the time-averaged quantities, which may be found as the zero'th Floquet components, and we shall narrow down the rather large parameter space to illustrate some of the most interesting time-averaged current phase relations realized by this driven junction.
In practice, the inversion of ( 18) is carried out by truncating to the n max lowest Floquet modes, i.e. working with square matrices of dimension 4(1 + 2n max ).For all numerical results presented below, we ensure that n max is large enough that increasing it further does not affect the results.Furthermore, we use a finite broadening in the lead Green functions, replacing 0 + by η = 10 −4 in Eq. ( 16), which, like all energy and frequency ( ≡ 1) parameters used below (except for the infinite-gap limit), is specified in units of ∆.In order to facilitate the numerical integration over the sharp sub-gap states in the infinite gap limit, both levels are weakly coupled to a normal metallic lead with chemical potential aligned with the two superconducting leads, µ m = 0.This gives rise to a finite imaginary part, Γ m , of the d-electron self-energies (15), which is chosen to be smaller than any other scale in the problem, yet resolved by the discretized numerical integrations.In practice, this corresponds to a finite parity relaxation time, which is longer than any other timescale in the problem.As discussed in the previous section, this also constitutes the formal justification of the steady state Dyson equation (11).For a finite gap, the continuum of the superconducting leads provides the necessary broadening for the numerical calculations, and the normal metallic lead is not needed.

A. Infinite, and large-gap results
In order to connect to the results of the previous section, we first consider the infinite-gap limit, in which all current is carried by Cooper pairs, at weak tunnel coupling and close to resonance.The resulting current, J = J 0 L − J 0 R (see Eq. ( 6)), is shown in Fig. 3 for different frequencies around the resonance.It is seen to match the perturbative results very well.We plot in Fig. 4 the dependence of the current near the resonance on the two phases, ϕ sc and θ d .We show this together with two cuts illustrating a good match with the result obtained in Eq. (46).
Increasing the amplitude of the drive and fixing the driving phase shift at θ d = π/2, gives rise to highly non-trivial CPRs, of which a few examples are shown in Fig. 5.For a small driving amplitude the CPR is modified by narrow dips of the current, similar to what is observed in superconducting junctions with only a single drive 57,58,[81][82][83] , the main difference being that in this case the dips are not symmetric around ϕ sc = π and do not reach zero, similar to the phase-shifted results in Ref. 62.For higher driving amplitudes (A = 0.8Γ shown) the current is reduced, as for the junctions with only a single drive, but now the CPR is severely modified with no special significance of neither ϕ sc = 0 nor ϕ sc = π, both exhibiting finite supercurrent.
For comparison, in Appendix A we calculate the current using the same parameters as for the blue curve (A = 0.1Γ) in Fig. 5, but now using Floquet states to determine the time evolution of the non-driven ground state.This is done in the infinite-gap limit and with no coupling to a normal metal (Γ m = 0).The long-time average of the resulting current shows good correspondence with the steady-state current in Fig. 5. Furthermore, interactions are straightforwardly included in this approach, and are shown to remove the sharp dips in the current if U becomes of the order of the driving frequency Ω.
The systematic behavior relies on many parameters, but the interaction seems to merely change the resonance condition for the drive.

B. Finite gap results
Turning to the case of a finite BCS gap, ∆, we first fix the superconductor phase difference to zero (ϕ sc = 0) and calculate the frequency dependence of the timeaveraged current.This is shown in the upper left panel of Fig. 6.The lower left panel shows the corresponding time-averaged density of states on the proximitized dots, exhibiting pronounced peaks at two slightly different ABS energies, together with their weaker first, and even weaker second Floquet sidebands.The small peak in the current at the highest frequency corresponds to a resonance between the first side band of the two ABS and the BCS quasiparticle continuum, at around Ω = ∆ + E ABS 1.5∆.At lower frequencies the current attains its largest magnitude slightly off-resonance, and a node right at resonance, Ω 2E ABS 1.2∆.Here a positive ABS energy matches the first sideband of a negative ABS energy, or vice-versa, as illustrated in Fig. 2 for weak t d .For lower frequencies, crossings of ABS sidebands with each other or with the continuum are again reflected in the current.Apart from this additional structure arising from the finite gap or from a substantial driving amplitude, the overall frequency dependence of the current clearly resembles the resonant structure found in the infinite-gap limit in Fig. 3.For the rest of the paper, we fix the drive frequency to be slightly deviated from the main resonance at Ω 2E ABS so as to focus our attention to this anomalous supercurrent signal.
To further investigate the effect of the continuum in the anomalous current we show in the right panel of Fig. 6 how the current varies with ∆ and Ω. Coming from high ∆, the resonant min-zero-max structure observed in Fig. 3 and in the upper left panel of Fig. 6 persists down to a ration of approximately ∆/Γ 2 − 3, with only a slight shift in the resonance frequency.For low ∆ the current is very small and frequency independent, consistent with a weak adiabatic pumping of normal current, similar to the case with normal leads 84 .The line separating the two regions corresponds to the condition for resonance between the continuum and the first Floquet sideband to the negative energy ABS, i.e.Ω = ∆+E ABS , where the energy of the ABS itself also depends on the superconducting gap 44,51,57 .A second line with twice the slope is also observed: it corresponds to a resonance between the second Floquet sideband and the quasiparticle continuum, beyond which the resonances carrying supercurrent are not modified.Interestingly, this second sideband is observed to anticross with the first sideband of the positive energy ABS, which gives rise to a large enhancement of the negative resonant current peak at frequency just below Ω = 2E ABS .Since this anticrossing involves sidebands crossing with the continuum, this enhancement of the current is most likely due to a dissipative quasiparticle current.On the other hand, the nearly vertical features in this figure, including the most pronounced min-zero-max resonance, correspond to a current of Cooper pairs, which are being pumped across the junction by phase shifted resonances between sub-gap states and their Floquet sidebands like the one indicated in Fig. 2.

C. Modified and rectifying current phase relations
As established in the Appendix A for single-state time evolution in the infinite-gap limit, the symmetries of the Floquet Hamiltonian guarantee the following symmetries of the time-averaged current: As we shall see below, these symmetries are still obeyed in the steady state Green function calculations with a finite BCS gap.The first symmetry relation, Eq.( 51), is apparent in the numerical Green function result for the time-averaged current shown in Fig. 7.For ε d = 0, the vertical black dashed cut illustrates the usual antisymmetry around ϕ sc = π.This symmetry breaks down for ε d = 0 and, as indicated by the vertical green dashed cut, may even lead to a unidirectional supercurrent, corresponding to complete rectification.The horizontal black dashed cut, on the other hand, illustrates the antisymmetry of the current under inversion of ε d for ϕ sc = 0.
In Fig. 8 we illustrate the second symmetry relation (52) of the current under inversion of both phases, ϕ sc and θ d .The vertical red dashed cut shows the anomalous relation between time-averaged current and phase difference on the drives, θ d , at a superconductor phase-difference fixed at ϕ sc = 0, attaining its maximum near, but not right at θ d = π/2.The three horizontal (black, blue and green) dashed cuts illustrate the the strongly modified current phase relations between the time-averaged current and the superconductor phasedifference.Switching from θ d = 0 to θ d = π, the driven Josephson junction is seen to switch current phase relation from a 0-, to a π-junction, as seen in the black, and the blue curves, respectively, up to a slight anharmonicity in both.
Once again, the green cut realizes a rectified timeaveraged current.Since the BCS gap is finite, there is no guarantee that this completely rectified pump current is exclusively a current of Cooper pairs.Nevertheless, as we show in Fig. 11 in Appendix A, a completely rectified current can also be obtained in the infinite-gap limit where all current must be carried by Cooper pairs, indicating that there is no fundamental obstacle to attaining a unidirectional time-averaged supercurrent for all ϕ sc .

VI. DISCUSSION
As demonstrated above, bridging two superconductors by a double quantum dot with phase shifted microwave tones on their respective gate voltages, as depicted in Fig. 1, comprises an effective Josephson junction with a highly nontrivial CPR.More specifically, the driving induces an alternating tunnelling current, which may exhibit a well-defined non-zero long-time average, and it is this average current, which exhibits an anomalous and often highly anharmonic relation to the superconductor phase difference.In light of the recent interest in Josephson diodes [85][86][87] , it is worth stressing that this driven junction offers complete rectification of the time-averaged supercurrent.
The supercurrent response to the driving relies on nonadiabatic resonant photon assisted tunnelling.This was established in the infinite-gap limit by means of perturbation theory and by time evolution of the non-driven ground state using Floquet theory (cf.Appendix A).For a finite BCS gap, a steady state time-averaged current was calculated by means of Floquet Keldysh Green functions, for which a weak tunnel coupling of each dot to a normal metal was included to eliminate the transient response and allow for parity relaxation.Whereas the general finite-gap current may include some fraction of BCS quasiparticles, the main resonant pump current arising when the drive frequency is slightly off resonance with the energy difference between the two even-parity sub-gap ABS was argued to be carried mainly by Cooper pairs.
For clarity, we have restricted our analysis to a symmetrically coupled device, where only the microwave phase shift breaks the L/R−inversion symmetry.Even in this case, the mean current exhibits a highly non-trivial behavior on the remaining parameters, such as the common mean gate voltage, the inter-dot tunnel coupling together with amplitude and frequency of the microwave tones.The plots chosen to illustrate the salient features for this work therefore by no means exhaust the many possible behaviors of this driven DQD junction.With two different tunnel couplings to the two superconductors, the induced ABS will have different energies and consequently the resonance frequencies on the two quantum dots will be different.The mechanism underlying the resonant current response will, however, remain viable if the two drive frequencies may be adjusted independently.
As demonstrated for the infinite-gap limit in Appendix A, local Coulomb interactions, reflecting the finite charging energies of the quantum dots were shown to alter the resonance conditions and thereby affect the timeaveraged current.Nevertheless, the anomalous Josephson effect (and the rectification) persisted, and was found to exhibit a 0 − π transition in θ d , as the interaction strength increased past a critical value.In real systems, this analysis pertains to the weak coupling regime, U ∆, whereas the opposite regime of U ∆ leads to the formation of YSR states.In this case, quasiparticles from the BCS continua in the leads form singlet bonds with spinful (odd-occupied) quantum dots, and this incomplete proximity effect must be expected to lead to substantial pumping of quasiparticle current.In the limit of U Γ, t d , however, charge fluctuations on the dot will be strongly impeded and the driving only effective at higher frequencies and amplitudes.It should be interesting to explore this regime further.Such work would also allow extending the results of Ref. 52 to a setup for two-tone spectroscopy 88,89 .
As discussed briefly at the end of Section IV, the timeaveraged currents addressed in this work must be expected to depend strongly on the parity flip dynamics arising from quasiparticle poisoning in a given device.Here, we have restricted our attention to the NESS Green function approach or to parity conserving time evolution of a single state, which in spite of their obvious differences all agree on the salient features of the driven DQD junction.For future work, it would be instructive to add stochastic parity dynamics to the Floquet time evolution carried out in Appendix A and make a comparison with the NESS Green function results.
The microwave enabled DQD Josephson junction studied here offers a highly tunable superconducting circuit element, which clearly links the phase-shifted AC input to a traversing supercurrent.Here we have only addressed the relations between time-averaged currents and superconducting, as well as microwave phase differences.To assess the possible value of such a circuit element, future work should address the time-dependent higher harmonics of the induced current as well as its possible implementation in a superconducting circuit.
In the course of finishing this work, we noticed the appearance of new work by A. Soori 90 , who also points to the Josephson diode aspect present in the driven twosite SNS junction explored also in Ref. 71.
Acknowledgments.The Center for Quantum Devices (Project No. DNRF101) and the Center for Nanostructured Graphene (Project No. DNRF103) are funded by the Danish National Research Foundation.We acknowledge fruitful discussions with G. Steffensen, K. Flensberg and M. Geier.
Appendix A: Floquet analysis of the interacting infinite-gap limit The infinite-gap limit offers relatively easy access to the symmetries of the problem, which are also revealed by the steady-state numerical calculations presented in the main text.In this appendix, we employ Floquet theory to provide a brief supplementary analysis of this more tractable limit, in which Local Coulomb interactions on the quantum dots can readily be included.Furthermore, since no quasiparticle excitations are involved in the infinite-gap limit, all currents calculated below are carried exclusively by Cooper pairs.We choose to consider only the even-parity sector, but a similar analysis is straightforwardly made for the odd-parity sector.
In the even-parity sector, the Hilbert space is spanned by the basis, {|00 with tunnelling matrix elements z = t d e iϕsc/2 , and with the local intra-dot Coulomb interaction, U , now included.From this, one may construct the even-parity Floquet Hamiltonian, H F e corresponding to the harmonic driving term, A cos(Ωt), from the matrix elements [91][92][93] where Î denotes the 6 × 6 unit matrix, and V is defined as the 6 × 6 matrix with diagonal elements The red curve in Fig. 5, corresponding to A = 0.7Γ, displays a finite anomalous Josephson current at ϕ sc = 0.In Fig. 10, we use the same parameters to show that this anomalous Josephson current depends strongly on the phase difference of the two drives, θ d , as found also with a finite BCS gap in the red curve of the right panel of Fig. 8. Here, however, one observes also a sign change of the anomalous Josephson current, corresponding to a transition from a 0-to π-junction behavior in θ d , when increasing the interaction strength.For the chosen parameters, this takes place at a critical interaction strength, U c ∼ Ω, but the more detailed parametric dependence of U c is beyond the scope of this paper.
Finally, with Fig. 11, we demonstrate that nearly complete rectification of the time-averaged current is possible also in the infinite-gap limit, where all current is carried by Cooper pairs.Unlike the finite-gap results shown in Figs.7 and 8, parameters have been fine tuned so as to make the current positive for all phase differences, ϕ sc .
In Fig. 11, we demonstrate also an explicit dependence of the current on the initial conditions as a spread in curves obtained for different Floquet gauges 94 , corresponding to different values of θ R .This is indicated by a set of some 63 gray curves, corresponding to evenly spaced values of θ R between 0 and 2π, which are averaged to obtain the blue curve.A similar spread will be obtained for the curves in Fig. 10 (not shown for clarity), whereas in Fig. 9, the driving amplitude is low enough that the results depend only on the phase difference, θ d .This spread increases with driving amplitude and gives a rough indication of sensitivity of the long time average of the Floquet time evolved current on initial conditions, and thereby whether they can be expected to be valid also within a driven steady state.

Symmetries of the current
The time-dependent current and thereby its longtime average obeys a few basic symmetries, which are most easily revealed by reverting to the timedependent infinite-gap Hamiltonian for the even sector obtained by replacing ε

FIG. 1 .
FIG. 1. Sketch of a Josephson junction with a structured weak link (gray region) based on a driven double quantum dot.The superconductors (blue) are maintained at a fixed phase-bias ϕsc = ϕL − ϕR, and the weak link is driven by two microwave gates with same amplitude and frequency, A, Ω, shifted in phase by θ d = θL − θR.The internal, and the two external tunnelling amplitudes are denoted by t d , tL and tR, respectively.

FIG. 2 .
FIG. 2. (a) A schematic of the four many-body eigenstates of a single proximitized undriven QD.The ground state has zero energy, the excited doublet has energy E d , and the twoquasiparticle state has energy 2E d .(b) Diagram illustrating the path of a Cooper pair through the driven DQD junction in progression from panel 1 through 5. Driving the microwave gates with Ω ∼ 2E d induces a near resonant transition from | 0 to | 2 in the left QD (1-2), followed by a two-step excitation transfer to the right QD via t d (2-4), which finally decays via its own microwave gate (4-5).

FIG. 3 .
FIG. 3. (Dashed lines) Plots of the weak coupling anomalous Josephson (ϕsc = 0) current in Eq. (46) in the infinite-gap limit as a function of the driving frequency, showing the maxima on each side of the node right at resonance, Ω = 2E d .Parameters are chosen such that t d = Γ/100, ε d = Γ/10 and θ d = π/2, together with three different driving amplitudes (see inset).(Full lines) Numerical calculation of the current using Eq.(6) with same parameters as above and with ∆ = 2 × 10 4 , η = 2 × 10 −4 and an additional broadening of the QD states corresponding to a normal metal tunnelling rate Γm = Γ/500.

FIG. 4 .
FIG. 4. (Upper panel) Density plot of the weak-coupling current vs. superconductor phase difference and phase shift of the two drives obtained by numerical evaluation of Eq. (6).(Lower panel) Solid lines correspond to cuts along the dashed blue (θ d = π), and green (ϕsc = 0) lines indicated in the upper panel, together with the corresponding analytical infinitegap weak coupling current from Eq. (46) limit (dashed).In both panels, parameters are t d = 0.02Γ, A = ε d = 0.1Γ and Ω = 2Γ.In the numerical evaluation the infinite gap was replaced by ∆ = 10 4 , while η = 10 −4 and Γm = Γ/300.

FIG. 6 .
FIG. 6. (Left) Upper panel: The time-averaged current (in units of I0 = 2e∆) at zero superconductor phase difference (ϕsc = 0) as a function of drive frequency.Lower panels: the corresponding time-averaged density of states on the two resonant levels.Both panels are evaluated with Γ = ∆/3, t d = 2∆/3, A = 4∆/15, ε d = ∆/15 and θ d = π/2.(Right) Variation with Ω and ∆ of the pumped current with t d = 2Γ, A = 0.8Γ, ε d = 0.2Γ, ϕsc = 0, θ d = π/2 and nmax = 7.For high ∆ resonances are observed between different ABS, while for low ∆ the current is featureless, signaling an adiabatic origin of the current.The line separating these two regions corresponds to the resonance between the ABS and the continuum.A second line with twice the slope is also present, corresponding to the first harmonic of that resonance.The dashed line at ∆ = 3Γ is the cut shown in Fig. 6.

7 .
Variation with ε d and ϕsc of the pumped current with parameters 2Γ = t d = 0.7∆, A = 0.8∆, Ω = 0.9∆, θ d = π/2 and nmax = 7 .The current-phase relation is strongly modified by varying ε d and the antisymmetry for inversion of ϕsc and ε d .As in Fig 6, another case is shown in green.here the current-phase relation does not cross zero current.The blue curve shows, for ϕsc = 0, the variation of the pumped current with ε d /∆.

FIG. 8 .
FIG. 8. Variation with ϕsc and θ d of the pumped current with 2Γ = t d = 0.7∆, A = 0.8∆, ε d = 0.8∆, Ω = 0.9∆ and nmax = 11.The current-phase relation is strongly modified and it shows antisymmetry for inversion of both phases.The current vanishes at the black solid lines.We show four cuts (marked by dashed lines) where interesting features are observed.In black, with θ d = 0, the current is that of a π−junction while for θ d = π, shown in blue, the junction becomes a 0−junction again.Another interesting feature is shown in green, with θ d = 1.6π,where for certain parameters one can obtain a current-phase relation that does not vanish for any value of the superconducting phase.

FIG. 12 .
FIG.12.Diagram illustrating the lowest order transport processes leading to anomalous Josephson current.The processes are similar to those illustrated in Fig.2, but are shown here with two Floquet side bands to each quantum dot level and with current carrying tunnelling paths in both directions, which interfere destructively unless ε d = 0, or with interactions included, ε d + U/2 = 0.