Subgap dynamics of double quantum dot coupled between superconducting and normal leads

Dynamical processes induced by the external time-dependent fields can provide valuable insight into the characteristic energy scales of a given physical system. We investigate them here in a nanoscopic heterostructure, consisting of the double quantum dot coupled in series to the superconducting and the metallic reservoirs, analyzing its response to (i) abrupt bias voltage applied across the junction, (ii) sudden change of the energy levels, and imposed by (iii) their periodic driving. We explore subgap properties of this setup which are strictly related to the in-gap quasiparticles and discuss their signatures manifested in the time-dependent charge currents. The characteristic multi-mode oscillations, their beating patters and photon-assisted harmonics reveal a rich spectrum of dynamical features that might be important for designing the superconducting qubits.

The double quantum dots embedded on interfaces between various external leads have been proposed for possible spin 1 and spin-orbit quantum bits 2 . Specifically, the superconducting qubits 3 have been considered as promising candidates, making use of the bound states formed inside the pairing gap 4 . Their implementations could protect the parity of Cooper pairs on proximitized superconducting nonoscopic islands 5 . Further perspectives for the proximitized double quantum dots appeared with the topological superconductors 6 , where the zero energy ingap modes are protected by symmetry reasons. These Majorana-type quasiparticles could be used for constructing the charge qubit in a transmission line resonator (transmon) 7 and may be incorporated in the gate tunable superconducting qubits (gatemons) 8 . Readout by means of a switching-event measurement using the attached superconducting quantum interference devices has revealed quantum-state oscillations with sufficiently high fidelity 9 , that seems appealing for realization of quantum computing.
So far the static properties of in-gap bound states have been throughly investigated for the single and multiple quantum dots 10,11 and recently also for nanoscopic length atomic chains, semiconducting nanowires, and magnetic islands proximitized to bulk superconductors 12 . Their particular realizations in the double quantum dots (DQDs) have been experimentally probed by the tunneling spectroscopy, using InAs [13][14][15][16][17][18] , InSb 19 , Ge/Si 20 and carbon nanotubes 21,22 and by the scanning tunneling microscopy applied to various di-molecules deposited on superconducting substrates [23][24][25][26][27] . Rich properties of such in-gap bound states of the DQDs have been analyzed theoretically by a number groups 11,19,[28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43] . Major features of two quantum dots coupled in series to the superconducting lead(s) originate from the ground state configuration which can vary its even-odd parity, depending on: the energy levels, hybridization with the external reservoirs, the inter-dot coupling, and the Coulomb potential 30,38 . Such parity changes are corroborated by crossings of the in-gap bound states and can be empirically detected by discontinuities of the Josephson current in S-DQD-S junctions [14][15][16] or the subgap Andreev current in N-DQD-S junctions 14,18,19 . The resulting zero-bias conductance as a function the quantum dot levels (tunable by the plunger gates) resembles a honeycomb structure [14][15][16] instead of a diamond shape, typical for the single quantum dot junctions. Influence of the coupling to external reservoirs is also meaningful. For instance in a regime of the strong coupling to superconducting lead(s) the spin of quantum dots would be screened 14 .
In general, various arrangements of two quantum dots enable realization of the on-dot and inter-dot electron pairing, affecting the measurable charge transport properties 36 . In particular, for the singly occupied quantum dots (what can be assured by appropriate gating) the superconducting proximity effect could be blocked. Such triplet blockade effect has been recently reported in S-DQD-S 17 and N-DQD-S 18 nanostructures. As regards (1) H =Ĥ S +Ĥ S−QD 1 +Ĥ DQD +Ĥ N−QD 2 +Ĥ N . Figure 1. Schematics. Two quantum dots (QD 1 and QD 2 ) coupled in series between the superconducting (S) and normal (N) metallic reservoirs whose energy levels ε iσ (t) could be varied by the external gate potential. We also consider dynamical phenomena driven by the time-dependent bias voltage imposed between the external leads. www.nature.com/scientificreports/ We restrict our considerations to the wide-band limit, assuming the constant (energy-independent) auxiliary couplings Ŵ N/S = 2π k/q |V Nk/Sq | 2 δ(ε − ǫ Nk/Sqσ ) . We also treat the pairing gap SC as the largest energy scale, focusing on dynamical processes solely inside in the subgap regime. In the limit of infinite | | the selfenergy of the Nambu-matrix Green's function becomes static and the value Ŵ S /2 appearing in the off-diagonal terms can be interpreted as the proximity induced pairing potential. The resulting low-energy physics can be described by 53 In what follows we discuss the time-dependent charge currents j Nσ (t) , j Sσ (t) and occupancies of the quantum dots imposed by the following types of quantum quenches: (i) abrupt bias potential V sd = µ N − µ S applied between N and S electrodes, (ii) sudden change of the energy levels ε iσ due to the gate potential V g , and (iii) periodic driving of the quantum dot levels with a given amplitude and frequency. Expectation values of the physical observables are computed numerically, solving a closed set of the differential equations for appropriate correlation functions (see "Methods" section). The charge current j Nσ (t) flowing between the normal lead and QD 2 can be derived from the time-dependent number of electrons in the normal lead. For ε Nkσ (t) = ε Nkσ this current is formally given by 54 were . . . denotes the quantum statistical averaging and �n iσ (t)� ≡n iσ (t) . The interdot charge flow j 12σ (t) is expressed as whereas the current j Sσ (t) flowing from the superconducting lead to QD 1 can be obtained from the charge conservation law dn 1σ (t) dt = j 12σ (t) + j Sσ (t) . Using equation (4) for the current j Nσ we can define its time-dependent differential conductance G Nσ (V sd , t) = d dV sd j Nσ (t) as a function of the source-drain voltage V sd . Peaks appearing in the dependence of G Nσ (V sd , t) against V sd can be interpreted as the excitation energies between eigenstates, comprising even and odd number of electrons (dubbed the Andreev bound states). Upon approaching the steady limit, t → ∞ , they emerge in the uncorrelated system at energies E = ± 1 2 4V 2 12 + Ŵ 2 S /4 ± Ŵ S 2 (see "Methods" section) and acquire a finite broadening caused by the relaxation processes on continuous spectrum of the normal lead.
In practical realizations of such N-DQD-S heterostructure ( Fig. 1) one should also take into account the Coulomb repulsion between electrons, i=1,2 U i n i↑ n i↓ , competing with the proximity-induced electron pairing and thereby affecting the bound states. Some aspects of the correlations effects have been previously studied under the stationary conditions for this heterostructure by the numerical renormalization group method 30 . Here we shall address the post-quench dynamics, treating the electron-electron interactions within the Hartree-Fock-Bogoliubov decoupling scheme This approximation applied to the static case of the correlated quantum dot hybridized with superconducting lead(s) can qualitatively describe the parity crossings and the energies of in-gap bound states 55 . We use of this decoupling (6) to provide a preliminary insight into the complicated quench-driven dynamics of the interacting setup, which is effectively described by with the renormalized energy levels ε iσ (t) = ε iσ (t) + U i n iσ (t) and the effective on-dot pairings Such mean-field approximation might be reliable at least for the weak interaction case. More subtle analysis, including the Kondo effect of the strongly correlated system ( U i ≫ Ŵ S ), is beyond a scope of this paper. We have done numerical calculations for U 1 = U 2 ≡ U , considering U/ Ŵ S = 0.5 , 1 and 1.5, respectively. Technically we have adapted for this purpose the algorithm outlined in "Methods" section, extending the previous study of the single dot superconducting junctions 46,54 .
We use the convention e = = 1 , expressing the charge currents, time and frequency ω in units of eŴ S / , / Ŵ S and Ŵ S / , respectively. In realistic experimental situations the value of Ŵ S ∼ 200 µ eV would imply the following typical units of time ∼ 3.3 psec, current ∼ 48 nA and frequency ∼ 0.3 THz. We assume the superconducting lead to be grounded, treating its chemical potential as the convenient reference level ( µ S = 0 ). Our calculations are performed for zero temperature. www.nature.com/scientificreports/ Response to a bias voltage. For computational reasons it is convenient to assume that initially, at t = 0 , the quantum dots are disconnected from both external reservoirs (see "Methods" section). Figure 2a presents the transient currents j Nσ (t) and j Sσ (t) right after forming the N-DQD-S heterostructure. In analogy to the previously discussed N-QD-S case 54 such evolution to the stationary limit is achieved through a sequence of the damped quantum oscillations, whose frequencies coincide with the energies of in-gap bound states. In particular, for ε iσ = 0 the period of such oscillations is equal to T = 4π/ Ŵ S and the relaxation processes (originating from the coupling Ŵ N of QD 2 to the metallic lead) impose the damping via exponential envelope function e −tŴ N /2 . In practice, at times t ≥ 50 , the stationary state seems to be fairly well approached. Let us turn to the dynamical response of N-DQD-S setup induced by its biasing, at t = 60 , when the chemical potentials are detuned by by source-drain voltage µ N − µ S = V sd . Figure 2b presents the charge currents j Nσ (t) and j Sσ (t) obtained for V 12 / Ŵ S = 2 , assuming V sd / Ŵ S = 1.5 , 2 and 20, respectively. For the large bias voltage, |V sd | ≫ V 12 , we observe emergence of the quantum beats with the period T B = π/V 12 superimposed with the higher frequency oscillations. Let us recall that charge transport is provided here solely by the anomalous particleto-hole (Andreev) scattering, which is sensitive to the in-gap bound states. For the particular set of model parameters such in-gap bound states appear at energies ± 1 2 4V 2 12 + Ŵ 2 S /4 ± Ŵ S /4 . It has been previously shown 56 that the single quantum dot placed between both normal electrodes responds to a sudden external voltage by the coherent oscillations of the charge current with frequency ω = |V sd − ε dot | . In the present situation we should replace ε dot by the effective in-gap quasiparticle energies, at which the Andreev scattering is amplified. We have four such in-gap bound states, therefore total current can be viewed as a superposition of sinusoidal waves, oscillating with the frequencies � 1/2 = V sd ± ω 1 and � 3/4 = V sd ± ω 2 , where ω 1/2 = V 12 ± Ŵ S /4 . It can be effectively expressed as 4 i=1 a i e − i t sin(� i t) . Individual terms refer here to the damping processes with different parameters i , whereas the coefficients a i control the contributions from these in-gap bound states. For the large bias |V sd | ≫ V 12 and |V sd | ≫ Ŵ S /4 the quantum beats are superimposed with the faster oscillations. It can be shown 56 that such beating patterns depend on a ratio For the case displayed in Fig. 2b this ratio is r = 8 , therefore for V sd / Ŵ S = 20 the repeated sequences of the beats with the periods π 4 , π 2 , π 2 , π 2 , π 2 , π 2 , π 2 , π 2 , π 4 appearing in the current j Nσ (t) should be observed. For non-integer ratio r the resulting beating pattern is more complicated with the different successive periods. Figure 2b displays that for V sd / Ŵ S = 20 the post-quench current j Nσ (t) indeed exhibits the beats mainly with period T B = π/V 12 superimposed with the faster oscillations, whose frequency is equal to V sd . The steady limit current obtained for V sd / Ŵ S = 2 is larger than for V sd / Ŵ S = 1.5 because of the broader transport window involving all the in-gap bound states. We also notice that j Sσ (t) substantially differs from j Nσ (t) , especially for the large bias V sd . We assign this to the fact that DQD sandwiched between the external leads wash out small fluctuations of the current j Sσ (t) , enforcing the final damped oscillations with period 4π/ Ŵ S . Figure 3 shows the beating structure in the time-dependent current j Nσ (t) after abrupt application of the bias voltage. These beats clearly depend on the interdot coupling V 12 via T B = π/V 12 . The beating structure is superimposed with oscillations whose frequency is also sensitive to the bias voltage. By measuring the period of such beating oscillations one could thus practically evaluate the inter-dot coupling V 12 = π/T B . For a realistic value Ŵ S ∼ 200 µeV , and assuming V 12 / Ŵ S = 0.5 , 1 and 2 the beating period would be T B ∼ 21 , 10 and 5 picoseconds, respectively. This time-scale is currently attainable experimentally. We have also performed similar www.nature.com/scientificreports/ calculations including the electron correlations (within the mean-field approximation assuming ε iσ = −U/2 ) and found, to our surprise, that all conclusions concerning the frequencies and the beating patterns remain valid.

Quench of energy levels.
Let us now consider the dynamics induced by a sequence of quantum quenches imposed on the energy levels ε iσ . The first quench ε iσ → ε iσ + V g is performed at t 1 = 60 , safely after N-DQD-S heterostructure achieves its stationary configuration. Later on, at time t 2 = 120 , we rapidly change the energy levels back to their initial values ε iσ + V g → ε iσ . Such step-like change (reminiscent of the pump-and-probe techniques) could be practically driven by the external gate potential applied to DQDs. For understanding the dynamics of our setup it is helpful to inspect the stationary fillings of both quantum dots for various interdot couplings V 12 . Figure 4 shows the occupancy of QD 2 (the neighbor of the normal lead) with respect to the energy level ε 2σ , assuming ε 1σ = ε 2σ so that occupancies of both dots are nearly identical. We recognize three plateau regions, corresponding to n iσ ≈ 1 , 0.5 and 0, respectively. We also notice, that a width of the half-filling region strongly depends on the inter-dot coupling V 12 . The stationary occupancy n 2σ changes from the nearly complete filling to half-filling or from the half-filled case to nearly empty QDs occur in a vicinity of ε iσ ≈ ±V 12 where the in-gap bound states coincide with the chemical potential µ N = µ S (here V sd = 0 ). Our numerical results obtained for various V 12 and V g indicate that the most prominent changes of the time-dependent observables occur for such quenches when the final value of the energy levels ε iσ coincide with the changeovers of n iσ (t = ∞) illustrated in Fig. 4. We have also checked that postquench evolution for different interdot couplings V 12 preserves the same universal properties, provided that the final value ε iσ corresponds to the tilted part of n iσ (t = ∞) curve.   Figure 4. Charge occupancy. The stationary limit (t = ∞) of the QD 2 occupancy as a function of the energy level ε 2σ = ε 1σ determined for several interdot couplings V 12 . The dashed line is calculated within the meanfield approximation for U = 1 . Other parameters: www.nature.com/scientificreports/ Figure 5 shows the time-dependent n 2σ (t) , j Nσ (t) , and j Sσ (t) after lifting the DQD energy levels, at t = 60 , and their return to initial values, at t = 120 , obtained for the strong interdot coupling, V 12 / Ŵ S = 4 . For t ≤ 60 the considered N-DQD-S system is practically in its stationary state with the half-filled QDs and negligible currents j Nσ (t) , j Sσ (t) . More specifically, we have chosen V g / Ŵ S = 3.2 , 3.8, 4, and 5, respectively. Such values of V g correspond to the stationary occupancies equal to ∼ 0.48 , ∼ 0.4 , ∼ 0.25 and ∼ 0.015 , respectively. Let us consider the postquench evolution corresponding to V g / Ŵ S = 3.2 , when the quantum dot level ε iσ coincides with the middle plateau (Fig. 4). The initial occupancy of QD 2 is 0.5 and its stationary value after the first quench (at t = 60 ) changes to ∼ 0.48 , therefore n 2σ (t) exhibits only very small oscillations. Similarly, the charge currents j Nσ and j Sσ are negligible (see the upper curves in Fig. 5 for t < 120 ). After the second quench (at t = 120 ) the occupancy n 2σ ∼ 0.5 , albeit promptly after the quench we observe some transient phenomena with the beating structure (see the upper curve in Fig. 5 for t > 120 ). This beating structure is more evident for the larger gate potentials V g / Ŵ S = 3.8 and 4 (see Fig. 5). We observe oscillations with the period T = π/V 12 , giving rise to the beating structure with another period 2π/ Ŵ S . Upon increasing the gate potential to V g / Ŵ S = 5 the timedepenence of n 2σ after the first quench substantially changes in comparison with the previous cases. Instead of the damped oscillations we now observe an exponential decay, down to nearly zero. Evolution after the second quench is also different in comparison to the previous ones. We now observe the oscillations of n 2σ and both currents with the period T = 2π/ Ŵ S without any beating structure. Concerning the time-dependent occupancies and currents calculated for V 12 / Ŵ S ≥ 1 , they preserve the qualitative properties discussed above. For the smaller interdot couplings V 12 (for instance V 12 / Ŵ S = 0.5 ) the evolution after the first quench preserves all properties characterized for stronger V 12 , but after the second quench we no longer observe the beating patterns, so that only oscillations with the period 4π/ Ŵ S are present.
We have also performed calculations for the interacting system, assuming U/ Ŵ S = 1 . The stationary limit occupancy of QD 2 is shown by the dashed line in Fig. 4. We can notice that the characteristic points, where the totally filled dot changes to the half-filling and another one where the half-filled dot changes to the empty configuration, are shifted in comparison to the noninteracting case. This effect is caused by rescaling of the in-gap states energies. In analogy to our considerations of uncorrelated system we have imposed such variations of the quantum dot levels by the gate potential V g which coincided with these characteristic points of n iσ (t = ∞) . It turned out that postquench evolution revealed the same qualitative features in the time-dependent occupancy n 2σ (t) and the charge currents as for U = 0 . For brevity, we hence skip such results.
Periodically driven energy levels. We now discuss dynamical response of the N-DQD-S heterostructure driven by a periodic driving of the energy levels ε iσ (t) = A sin(ωt) which can be practically achieved by shining an infrared field on the quantum dots. We assume that amplitude A and frequency ω of the oscillations are identical in both QDs. Figure 6 presents the time-dependent current j Sσ (t) obtained for ω/ Ŵ S = 0.1 and several values of the amplitude A. The left (a) panel refers to the uncorrelated case, U = 0 , and the right (b) panel to U/ Ŵ S = 1 , respectively. As a guide to eye we also display the transient current obtained for the static energy levels ε iσ = 0 (top panel in Fig. 5a) with the characteristic damped oscillations whose period is equal to 4π/ Ŵ S . Such current vanishes in the asymptotic limit t → ∞ (here V sd = 0 ) and similar features, but with different profiles of the quantum oscillations, are observable for small amplitudes of the periodic driving as well. They are displayed for V 12 / Ŵ S = 4 in Fig. 6a. We notice that indeed the time-dependent currents asymptotically vanish for A/ Ŵ S ≤ 3.5 . This situation occurs whenever the amplitude A does not exceed the energies of subgap quasiparticles. Such behavior can be contrasted with the larger amplitude driving (for instance A/ Ŵ S = 4 ) when the current j Sσ (t) is forced www.nature.com/scientificreports/ to flow back and forth all over the time. Periodicity is this behavior is a bit subtle and will be analyzed in more detail underneath. Figure 6b shows the current j Sσ (t) of the correlated system (Coulomb potential U 1 = U 2 = U is expressed in units of Ŵ S ) determined for V 12 / Ŵ S = 3 , A/ Ŵ S = 3 , and V sd = 0 . We have chosen such parameters to enforce the nonvanishing current, up to the asymptotic limit t → ∞ . The correlation effects are here quite evident. Upon increasing U the magnitude of oscillating current j Sσ (t) is gradually suppressed. Such effect can be partly assigned to shifting of the subgap quasiparticles to the higher energies and partly to ongoing transfer of the spectral weights (this behavior is also discussed in next subsection). In presence of the finite source-drain voltage V sd the time-dependent phenomena become even more complicated. Its seems, however, that under such highly non-equilibrium conditions the correlation effects become less important.
Finally we briefly investigate the transient currents imposed by different profiles of the periodically driven energy levels ε iσ (t) = ε iσ (t + T) as depicted by the dashed lines in Fig. 7. For all cases we have assumed the same amplitudes and frequencies. As the reference, the upper panel shows the case of the sinusoidally driven energy level. It appears that abrupt (step-like) variations of QDs energy levels are followed by the damped oscillations of transient current j Sσ (t) after each change of ε iσ . Life-time of the resulting damped oscillations is shorter or comparable to the period of driving. For more smooth variation of ε iσ we can notice gradual suppression of the induced oscillations (see the second panel from the top of Fig. 7).

Andreev conductance averaged over driving period.
To gain more precise information about the role of amplitude A and frequency ω of the oscillating QDs energy levels we study here the charge currents averaged over a period T = 2π/ω of the driving field. Our main objective is to analyze the spectrum of subgap qua- j Nσ (t)dt induced by the source-drain voltage V sd and, in analogy to the preceding section, assuming the periodically driven energy levels ε iσ (t) = A sin(ωt) . From the differential conductance G Nσ (V sd ) = d dV sd �j Nσ (t)� t 0 one can infer quasienergies of the in-gap bound states 57 . Initially, at t = 0 , the oscillating quantum dot levels ε iσ (t) are imposed simultaneously with the bias voltage µ N − µ S = V sd , assuming both QDs to be empty. We choose the reference time t 0 at which the transient effects become negligible. This choice can be quite arbitrary, because safely after forming our N-DQD-S heterostructure the time-dependent current oscillates with the same period T as enforced on the energy levels (c.f. Figs. 6, 7). Below we discuss the differential conductance G Nσ (V sd ) obtained numerically for a few representative sets of the model parameters. Figure 8 presents the averaged Andreev conductance obtained for two values of the interdot coupling V 12 and several amplitudes A, as indicated. Panels (a-d) display the characteristic features originating from the photonassisted tunneling. We notice that besides the main quasiparticle peaks (for Ŵ N ≪ Ŵ S ) appearing at ± 1 2 4V 2 12 + Ŵ 2 S /4 ± Ŵ S /2 there emerge additional side-peaks originating from the stimulated emission/ absorption of the photon quanta. Their intensity (spectral weight) and avoided-crossing behavior are sensitive to the frequency and amplitude of a microwave field. The main quasiparticle peaks are replicated at multiples of ω and they can be interpreted as higher order harmonics of the initial bound states.
Basic aspects of the photon-assisted tunneling through the quantum dots sandwiched between the normal electrodes have been extensively studied in literature [58][59][60] , predicting the main resonance peaks and their n-th side-bands modulated by the squared Bessel functions of the first kind J 2 n (A/ω) . As regards the specific photonassisted tunneling in the superconducting junctions, it has been observed that the differential conductance G(V sd ) in situations with the single quantum dots can be expressed by G(V sd ) = n J 2 n (kA/ω)G (0) (V sd + nω k ) , where G (0) (V sd ) corresponds to the conductivity without microwave radiation and k denotes the number of electrons transferred in an elementary tunneling process 47,48 . For our N-DQD-S nanostrocture we notice that the main resonant peaks and their side-bands are weighted by the squared Bessel function J 2 0 2A ω . The main resonance peaks and side-bands disappear at such frequencies ω for which the Bessel function vanishes. Figure 8d shows such points for ω/ Ŵ S ∼ 3.3 , 1.45, 0.92, corresponding to the first, second and third zeros of J 0 (2A/ω) . For some given amplitude A the frequency ω at which the main quasiparticle peaks and their higher harmonics disappear is independent of the interdot coupling V 12 (Fig. 8c,d).  www.nature.com/scientificreports/ Let us now consider variation of the averaged Andreev conductance G Nσ with respect to ( V sd , A ) for a few values of the interdot coupling V 12 (Fig. 9). In absence of the microwave field, A = 0 , there exist four peaks in the differential conductance corresponding to two pairs of in-gap bound states. Upon increasing a power of the microwave field (for larger amplitude A) the main quasiparticle peaks loose some part their intensities (spectral weights) at expense of their new higher-order replicas. By varying the amplitude A such replicas appear in the averaged conductance at ±ω , ±2ω , and so on around the main peaks. We can also notice that their spectral weight undergoes substantial redistribution. In particular, at certain values of the amplitude A the spectral weight of individual harmonics vanishes and then reappears.
To check influence of the inter-dot coupling V 12 on the averaged Andreev conductance we present in Fig. 10 the results obtained for ω/ Ŵ S = 1 and two amplitudes A/ Ŵ S = 1 and 2. In the first case the peaks, appearing around ±nω , gradually split into the lower and upper branches with the increasing coupling V 12 . Yet, they never cross each other because of the quantum mechanical interference 61 . For the larger amplitude, A/ Ŵ S = 2 , we clearly notice such avoided-crossing tendency, where each harmonic consists of two nearby located peaks. This is an example of the n-fold fine structure driven in the harmonics, whenever the specific constraint A/ω = n is encountered.
Finally, in Fig. 11we present the averaged conductance G Nσ (V sd ) of the interacting system obtained for V 12 / Ŵ S = 2 , A/ Ŵ S = 2 , ω/ Ŵ S = 2.5 , where panels form top to bottom refer to U/ Ŵ S = 0 , 0.5, 1, and 1.5, respectively. The particle-to-hole scattering mechanism (contributing to the subgap Andreev current) implies the fully symmetric conductance G Nσ (−V sd ) = G Nσ (V sd ) . In the uncorelated system (top panel) the main quasiparticle peaks appear at ± 1 2 4V 2 12 + Ŵ 2 S /4 ± Ŵ S 2 and their higher order replicas are spaced by ±nω . For www.nature.com/scientificreports/ the presently chosen parameters the second-and higher-order harmonics become hardly visible because of their very small spectral weights (see Figs. 8b and 11). Upon increasing the Coulomb potential U the main quasiparticle peaks only slightly change their positions. Major influence of the correlation effects is manifested through noticeable redistribution of the spectral weights, both between the harmonics and between their fine sub-structure. More detailed analysis of the photon-stimulated Andreev transport of the strongly correlated N-DQD-S system would require some sophisticated (nonperturbative) techniques, and such study is beyond the scope of the present work. In addition to the numerical computations of the averaged current directly from the equations of motion, we have also developed the auxiliary procedure based on machine learning algorithm which reliably yields the Andreev conductance for an arbitrary set of the model parameters (see the last subsection of "Methods" section).

Discussion
We have studied the double quantum dot coupled between the superconducting and normal leads, addressing its dynamical response to (i) abrupt application of the bias voltage, (ii) sudden change of the energy levels, and (iii) their periodic driving. These effects can be routinely triggered either by dc or ac external potentials. We have analyzed the time-dependent charge flow between the external reservoirs and the quantum dots, revealing an oscillatory behavior (analogous to the Rabi-type mechanism involving pairs of the in-gap quasiparticle states induced by the superconducting proximity effect) with a damping caused by the relaxation processes on a continuum spectrum of the normal lead.
Inspecting the time-dependent profiles of various physical observables we have found the signatures of such frequency components which coincide with the subgap quasiparticle energies. For the quantum quench imposed by the source-drain voltage and by the gate potential the dynamics of proximitized double quantum dot reveals superposition of the fast and slow oscillatory modes, giving rise to the beating patters. These features are well observable over quite long time interval, �t ∼ 10 / Ŵ N , in contrast to much faster transient phenomena realized in the single quantum dot (N-QD-S) heterostructures 54,62 .
In the case of periodically driven energy levels we have found more complex time-dependent behavior. Response of the N-DQD-S heterostructure depends both on the frequency ω and amplitude A of the periodically varying levels. We have illustrated these phenomena in absence (Figs. 6,7) and in presence of the bias voltage (Figs. 8,9,10). We have predicted that amplitude (related to the power of driving force) has crucial effect on activating the higher-order harmonics of in-gap quasiparticle sates, as evidenced for the unbiased (Fig. 6) and biased ( Fig. 9) heterostructures. The frequency, on the other hand, is manifested by replicas of the main quasiparticle peaks. Similar effects have been recently observed experimentally in the Josephson-type junctions, comprising the single quantum dot 47,48 . In our case the proximitized double quantum dot is characterized by a sequence of the photon-assisted enhancements in the differential conductance with an additional fine-structure appearing in the harmonics due to interference effects. Upon varying the frequency (Fig. 8) or the interdot coupling (Fig. 10) the neighboring harmonics never cross each other because of their quantum mechanical interference, which is feasible also in multi-terminal superconducting junctions 61 .
Our considerations could be verified experimentally by means of the subgap tunneling spectroscopy using the carbon nanotubes, semiconducting nanowires or other lithographically constructed quantum dots embedded between the superconducting and metallic electrodes. Another realization would be possible using the scanning microscope technique, where the conducting tip can probe the dimerized molecules deposited on superconducting substrates. The characteristic time-scales determined in this work might be important for designing logical operations with use of the superconducting qubits 8 . In future studies it would be worthwhile to perform more systematic consideration of the correlation effects and address the dynamics of topologically nontrivial superconducting nanostructures.

Equations of motion.
Here, we explicitly present the set of differential equations needed for determination of the time-dependent occupancy n iσ (t) = �ĉ † iσ (t)ĉ iσ (t)� and other functions coupled to it (for U=0). Using the exact formula and applying the wide band limit approximation we derive the following set of equations where α = +(−) , β = exp(−i(t − t 1 )V sd ) , t 1 denotes the time at which the bias voltage V sd is applied and . . . stands for the quantum statistical averaging. At this level there appear the new correlation functions �Â iσ (t)B kσ (0)� , where Â ( B ) corresponds to the creation or annihilation operator of electron in the quantum dots (the normal lead). These functions can be determined from the the following equations of motion where �n kσ (0)� = 1 + exp ((ε Nkσ − µ N )/k B T) −1 is the Fermi distribution function for the normal lead electrons. www.nature.com/scientificreports/ We have solved numerically these coupled differential equations (10)(11)(12)(13)(14)(15)(16)(17)(18)(19) subject to the specific initial conditions. For convenience, we have assumed that at t = 0 both external reservoirs were isolated from the quantum dots. In next steps, we have calculated iteratively the time-dependent observables using the Runge Kutta algorithm with sufficiently dense equidistant temporal points t → t + δt → · · · → t + Nδt ≡ t f .

Machine learning approach.
Results presented in the main part of this paper have been obtained by solving the differential equations derived for N-DQD-S heterostructure. The computational procedure has been rather straightforward (see the preceding section), but required quite a lot of time and resources. For instance to produce the conduction maps (Figs. 8,9,10) with 150 × 150 points resolution it takes approximately one week performing multiprocessing calculations on CPU 2x Xeon E5-2660 2.2GHz 16 cores/32 threads. This problem motivated us to construct a machine learning model for our system.
To train our neural network we have used the collected set of data of 76 different conductance maps (with different resolutions), giving us 971760 conductance data points. Subsequently, we have linearly interpolated every single map to doubly increase a number of the data points, finally giving us 3887040 data points. For this purpose we have used the open-source software for machine learning-Tensorflow with application programming interface-Keras.
This neural network has a character of the densely connected type, with 4 input parameters ( V 12 ,ω,V sd ,A) describing noninteracting N-DQD-S setup and 1 single neuron on the output, specifying the averaged Andreev conductance G Nσ . The neural network is composed of 4 hidden layers consisting of 2048, 1024, 512, 256 neurons, respectively. Every hidden layer has a dropout of 1% neurons (which helps to avoid over-fitting our model) and, as an activation function, we have used sigmoid function. One can notice that this neural network is large, because of non-linearity in the system. To train our neural network we have chosen batch = 1024 and epoch = 600 , giving us the fidelity coefficient R 2 = 0.987 . Fig. 12 compares the calculated G Nσ with respect to the value predicted by our neural network. Predictive strength of the machine learning algorithm is illustrated in Fig. 13, which shows the conductance maps obtained from the direct calculation (panel a) and by the neural network (panel b) for such model parameters which were not used during the training process. This neural network model of N-DQD-S heterostructure is available at the following https:// www. dropb ox. com/ sh/ 0hzs9 im3d3 bf0jr/ AADRr 3kltw 2mOdC Ch8te doIWa? dl=0www.dropbox.com/sh/0hzs9im3d3bf0jr/AADRr3kltw2mOdCCh8tedoIWa?dl=0 webpage.  www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.